function [Match1 Match2 matchedFP pocetFP] = SIFT( varargin )

if ischar( varargin{1} )
    img1 = imread( varargin{1} );
    img2 = imread( varargin{2} );
else
    img1 = varargin{1};
    img2 = varargin{2};
end

clear varargin

warning off all
if isrgb(img1)
    img1gr = rgb2gray( img1 );
else
    img1gr = img1;
end
if isrgb(img2)
    img2gr = rgb2gray( img2 );
else 
    img2gr = img2;
end

%% Premenné 
%-------------------- Parameters -------------------%
SwitchWaitbars = 'on'; % 'off'/'on'

%--------------- Descriptors 4 FP ------------------%
TypeOfFP = 'Aff'; % Parameter nema ziadny vplyv na proces registracie - Typ výstupu pre vysledne zobrazovanie 'Aff'/'HL' - 'HL' vykresli stvorce
TypeOfDescriptor = 'SIFT'; % !!! NEMENIT !!! - definicia typu deskriptora (implementovany len SIFT)
NOfWindows = 4; % Pocet okien ktore sa pouziju v okoli kazdeho bodu pre deskriptor - NOfWindows^2 s velkostou okna 4x4
AngleBinDescript = 45; % degrees, number of bins
ThreshDescript = 0.2; 
FactorDescript = 1.5; %
StepSampleFunction = @(x) sqrt(1+x^2); % !!! NEMENIT !!! Funkcia pre definovanie podvzorkovania pixlov v okoli kazdeho FP


%--------------------- Harris ----------------------%
TypeOfCornerDetector = 'HarmonicMean'; % Type detektora - 'Harris'/'HarmonicMean'
TypeOfNBHOOD = 'const'; % Typ okolia bodu - 'const'/'dif' | 'dif' - je závislé na mierke
NBHOOD = ones(3); % Scalarna alebo binarna hodnota matice, ktora bude reprezentovat okolie bodu
BorderDistance = 2*NOfWindows;
ThreshType = 'const'; % Typ prahu - 'const'/'percent'
HarrisThresh = 10; % Hodnota prahu | V pripade 'HarmonicMean', prah ~ 10; 'Harris' ~ 1000
k = 0.055; % konstanta pre vypocet Harrisovej funkcie
Dilate = 'no'; %Pouzitie kriteria minimalnej vzdialenosti bodu 'yes'/'no'
radius = 2; % Radius dilatacie bodov. Teda budu odstranene vsetky FP z okolia ineho FP ktore su blizsie ako nastavena hodnota
sigma_nmbr = 9; % Pocet zmensených obrazov - 2 - 17


%%
dispPlaces = 1; % 1/0 - shows only places of FP's, without scale, orientation.
dispMOCS = 0; % 1/0 - displays rectangles with edges proport. to charac. scale, and rotated as main orientation - takes a lot of time 
dispMatches = 1; % 0/1
% need to be or dispMOCS = 1 or dispPlaces = 1, or both = 0; both = 1, isn't supported ( i.e. have to be dispPlaces*dispMOCS = 0 ). 
%%
%--------------- Main Orientation -----------------%
ThreshMainOrient = 0.8;
FactorMainOrient = 1.5;
AngleBinMainOrient = 10; % can be interpretated, as main orientation will be calculated with accuracy = AngleBinMainOrient degrees
TypeOfMainOrient = 'some'; % 'one'/'some'. 'one' - for each FP take only one Main Orientation 
% ( take it, if this orient. is the only one bigger then ThreshMainOrient % of it's value )
% 'some' - take all orientations, that bigger then thresh, described above ( > max_value*ThreshMainOrient ).

%----------------- Matching FPs -------------------%
TypeOfSearch = 'determ'; %'determ' 1/2, 'kmeans','kNN' for not big data preferable 'determ' 2,
% 'kmeans' - just for experiment.
HelpScalarOrVector = 2;
TypeOfMatchThresh = 'first/second'; % 'first'  / 'first/second'  / 'first/second&first'
ThreshFS = 0.78; % thresh for first/second ( 0.7 - ~0.85 )
ThreshF = 0.48; % thresh for first ( 0.2 - 0.5 )

%--------------------------------------------------------------------%
%- Finding FP and they char. scale (Harris-Laplace)  -% 
% if dispPlaces && dispMatches
%     subplot(2,2,1)
% elseif dispPlaces
%     subplot(1,2,1)
% end
[ HrLPoints1 ] = harrislpls( img1gr, TypeOfNBHOOD, NBHOOD, TypeOfCornerDetector, BorderDistance, ThreshType, HarrisThresh, k, Dilate, radius, sigma_nmbr, dispPlaces, SwitchWaitbars );
% if dispPlaces && dispMatches
%     subplot(2,2,2)
% elseif dispPlaces
%     subplot(1,2,2)
% end
[ HrLPoints2 ] = harrislpls( img2gr, TypeOfNBHOOD, NBHOOD, TypeOfCornerDetector, BorderDistance, ThreshType, HarrisThresh, k, Dilate, radius, sigma_nmbr, dispPlaces, SwitchWaitbars );

%----------- Finding Main Orientation -------------%
[ HrLOrntPoints1 ] = mainOrient( HrLPoints1, img1gr, ThreshMainOrient, FactorMainOrient, AngleBinMainOrient, TypeOfMainOrient, SwitchWaitbars );
[ HrLOrntPoints2 ] = mainOrient( HrLPoints2, img2gr, ThreshMainOrient, FactorMainOrient, AngleBinMainOrient, TypeOfMainOrient, SwitchWaitbars );

%------------------ Descripting -------------------%
[ DescriptHrLPoints1 ] = descriptFPoints( HrLOrntPoints1, img1gr, TypeOfDescriptor, NOfWindows, StepSampleFunction, AngleBinDescript, ThreshDescript, FactorDescript, SwitchWaitbars );
[ DescriptHrLPoints2 ] = descriptFPoints( HrLOrntPoints2, img2gr, TypeOfDescriptor, NOfWindows, StepSampleFunction, AngleBinDescript, ThreshDescript, FactorDescript, SwitchWaitbars );

%-------------------- Matching --------------------%
switch TypeOfMatchThresh
    case 'first'
Thresh = ThreshF;
[ MatchedPairs ] = matchFPoints( DescriptHrLPoints1, DescriptHrLPoints2, TypeOfSearch, HelpScalarOrVector, TypeOfMatchThresh, Thresh, SwitchWaitbars );
    case 'first/second'
Thresh = ThreshFS;
[ MatchedPairs ] = matchFPoints( DescriptHrLPoints1, DescriptHrLPoints2, TypeOfSearch, HelpScalarOrVector, TypeOfMatchThresh, Thresh, SwitchWaitbars );
    case 'first/second&first'
Thresh = ThreshFS;
[ MatchedPairs ] = matchFPoints( DescriptHrLPoints1, DescriptHrLPoints2, TypeOfSearch, HelpScalarOrVector, 'first/second', Thresh, SwitchWaitbars );
ind1 = find( MatchedPairs );
ind2 = MatchedPairs( ind1 );
% take only those pairs, that are at distance each from other less then ThreshF
ind = find( sqrt(sum( ((DescriptHrLPoints1( ind1, : ) - DescriptHrLPoints2( ind2, : )).^2), 2 )) <= ThreshF );
ind1 = ind1( ind );
ind2 = ind2( ind );
end

%--------------------- 4 Ploting ------------------%
switch TypeOfMatchThresh
    case 'first/second&first'
Match1 = HrLOrntPoints1( ind1, [ 1 2 3 end ] );
Match2 = HrLOrntPoints2( ind2, [ 1 2 3 end ] );
    otherwise
ind1 = find( MatchedPairs );
ind2 = MatchedPairs( ind1 );
Match1 = HrLOrntPoints1( ind1, [ 1 2 3 end ] );
Match2 = HrLOrntPoints2( ind2, [ 1 2 3 end ] );
end

%------------- Removing the same pairs ------------%
[ b, ind, ind2 ] = unique( [ Match1(:, [1 2]), Match2(:, [1 2]) ], 'rows' ); clear b ind2
Match1 = Match1( ind, : );
Match2 = Match2( ind, : );

matchedFP(:,1:2) = Match1(:,1:2);
matchedFP(:,3:4) = Match2(:,1:2);
pocetFP = size(matchedFP,1);

% for i=1:(size(matchedFP,1))
%     distance = sqrt( ((matchedFP(i,1) - matchedFP(i,3))^2) + ((matchedFP(i,2) - matchedFP(i,4))^2) );
%     
%     if distance > 15
%        matchedFP(i,1:4) = 0;
%         
%     end
%     
%     
% end

% if dispMOCS && dispMatches
%     subplot(2,2,1)
%     PlotFP( img1, HrLOrntPoints1 );
% elseif dispMOCS
%     subplot(1,2,1)
%     PlotFP( img1, HrLOrntPoints1 );
% end
    
% if dispMOCS && dispMatches
%     subplot(2,2,2)
%     PlotFP( img2, HrLOrntPoints1 );
% elseif dispMOCS
%     subplot(1,2,2)
%     PlotFP( img2, HrLOrntPoints1 );
% end

% if dispMatches && ( dispMOCS || dispPlaces )
%     subplot(2,2,[ 3 4 ])
%     PlotMatches( img1, Match1, img2, Match2, 'sigma', TypeOfFP );
% elseif dispMatches
%     figure, PlotMatches( img1, Match1, img2, Match2, 'sigma', TypeOfFP );


end