الگوریتم سیفت: ( scale invariant feature transform )

در مبحث به دست آوردن نقاط ویژه که کاربرد فراوانی در انطباق تصویر دارد ، به دست آوردن نقاط ویژه هریس (harris point  )

از پایه ای ترین روش ها می باشد با این وجود عدم کارایی این روش و روش های مشابه در مقیاس های مختلف ( scale ) از تصویر استفاده از این روش ها رو به چالش می کشد.

با هدف رفع این مشکل ویژگی های SIFT ( Scale Invariant Feature Transfirm)  معرفی شدند. در این روش تصویر به جای محاسبه و استفاده از لاپلاسین گوسین LoG  از تفاضل گوسی DoG استفاده می شود و به جای استفاده از سیگماهای مختلف در فیلتر گوسی ، تفاضل گوسی روی اکتاو octave  مختلف از هرم تصویر اعمال می شود تا نقاط ویژه به دست آید .

در این بین برای حذف نقاط نا مناسب از یک سطح آستانه استفاده می شود که اگر حاصل تفاضل از این مقدار کمتر بود یعنی نقطه قدرتمندی نیست و در نظر گرفته نمی شود. از سوی دیگر لبه ها همچنین برای تشخیص نقاط مناسب مشکل ایجاد می کنند برای منظور در ماسک 2*2 تعریف شده مقادیر ویژه محاسبه می شود اگر دو تفاوت دو مقدار ویژه بیش از یک آستانه مشخص بود این نقطه لبه بوده و در نظر گرفته نمی شود.

همچنین برای هر نقطه زاویه یا جهت تعیین می شود که معمولا 36 جهت در نظر گرفته می شود که بر اساس گرادیان محاسبه می شود.در نهایت برای هر نقطه همسایگی 16*16 به صورت 16 بلوک 4*4 بررسی می شود که در هر بلوک هیستوگرام جهت ها به صورت 8 ستونی محاسبه می شوند پس در انتها به ازائ هر نقطه 128 مقدار بردار توصیف گر نقطه را تشکیل می دهند. در تصویر زیر به ترتیت از چپ به راست تصویر ورودی ، تصویر وسط تصویرنقاط استخراج شده و تصویر سمت راست جهت غالب در هر نقطه را نشان می دهد

نمونه ای از اجرای  الگوریتم sift   به صورت زیر است : 

%% Initialization

% Here, the x-axis correspond to coloum, while the y-axis correspond to row

clear all; close all;

img = imread('cameraman.tif');

[~,~,ColorChannel] = size(img);

if ColorChannel > 1

    img = rgb2gray(img);

end

img = double(img)/255;

ScaleSpaceNum = 3; % number of scale space intervals

SigmaOrigin = 2^0.5; % default sigma

ScaleFactor = 2^(1/ScaleSpaceNum);

StackNum = ScaleSpaceNum + 3; % number of stacks = number of scale space intervals + 3

OctaveNum = 3;

GaussianFilterSize = 21;

OctaveImage = {OctaveNum,StackNum}; % save the Gaussian-filtered results of image

OctaveImageDiff = {OctaveNum,StackNum-1}; % save the Difference of Gaussian-filtered results of image

%% Gaussian Convolution of Images in Each Octave

ImgOctave = cell(OctaveNum,1);

for Octave = 1:OctaveNum

    Sigma = SigmaOrigin * 2^(Octave-1); % when up to a new octave, double the sigma

    ImgOctave{Octave} = imresize(img, 2^(-(Octave-1)));

    for s = 1:StackNum

        SigmaScale = Sigma * ScaleFactor^(s-2);

        % calculate Guassian kernel

        GaussianFilter = fspecial('gaussian',[GaussianFilterSize,GaussianFilterSize],SigmaScale);

        % do convolution with Gassian kernel

        OctaveImage{Octave,s} = imfilter(ImgOctave{Octave}, GaussianFilter,'symmetric');

    end  

end

% calculate difference of Gaussian (in original paper eq.1)

for Octave = 1:OctaveNum

    for s = 1:StackNum-1

        OctaveImageDiff{Octave,s} = OctaveImage{Octave,s+1} - OctaveImage{Octave,s};

    end

end

%% Find the local minima and maxima between stacks

DiffMinMaxMap = cell(OctaveNum,StackNum-3);

for Octave = 1:OctaveNum

    for s = 2:size(OctaveImageDiff,2)-1

        CompareDiffImg = zeros(size(OctaveImage{Octave,1},1)-2,size(OctaveImage{Octave,1},2)-2,27);

        indx = 0; % 3rd dimension indx for CompareDiffImg

        for s2 = s-1:s+1

            for k = 1:9

                [i,j] = ind2sub([3,3],k);

                CompareDiffImg(:,:,indx+k) = OctaveImageDiff{Octave,s2}(i:end-3+i,j:end-3+j);

            end

            indx = indx + 9;

        end

        [~,MinMap] = min(CompareDiffImg,[],3);

        MinMap = (MinMap == 14);

        [~,MaxMap] = max(CompareDiffImg,[],3);

        MaxMap = (MaxMap == 14);

        DiffMinMaxMap{Octave,s-1} = MinMap + MaxMap; % the center indx is 9 + 5 = 14

    end

end

%% Accurate keypoint localization

UnstableExtremaThreshold = 0.03; % set the threshold as 0.03 (the same as the original paper)

DiffMinMaxMapAccurate = DiffMinMaxMap;

for Octave = 1:OctaveNum

    Sigma = SigmaOrigin * 2^(Octave-1);

    for DiffMinMaxMapIndx = 1:size(DiffMinMaxMap,2)

        Map = DiffMinMaxMap{Octave,DiffMinMaxMapIndx};

        SSCindx = find(Map); % Scale-Space-Corner Index

        if isempty(SSCindx)

            continue;

        end

        for ssc = 1:length(SSCindx)

            [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

            if (Row <= 1) || (Row >= size(Map,1)) || (Col <= 1) || (Col >= size(Map,2))

                DiffMinMaxMapAccurate{Octave,DiffMinMaxMapIndx}(Row,Col) = 0; % discard out of matrix boundary

                continue;

            end

            RowDiff = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row+1,Col) - OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row-1,Col);

            ColDiff = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col+1) - OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col-1);

            ScaleDiff = OctaveImageDiff{Octave,DiffMinMaxMapIndx+2}(Row,Col) - OctaveImageDiff{Octave,DiffMinMaxMapIndx}(Row,Col);

            offset = [2; 2; Sigma * ScaleFactor^(DiffMinMaxMapIndx) - Sigma * ScaleFactor^(DiffMinMaxMapIndx-2)];

            DxHat = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col) + 0.5 * ([RowDiff,ColDiff,ScaleDiff] * offset);

            if abs(DxHat) < UnstableExtremaThreshold

                DiffMinMaxMapAccurate{Octave,DiffMinMaxMapIndx}(Row,Col) = 0; % discard unstable extrema

            end

        end

    end

end

%% Eliminating edge responses

gamma = 10; % set the threshold gamma as 10 (the same as the original paper)

DetermineThreshold = (gamma + 1)^2 / gamma;

DiffMinMaxMapNoEdge = DiffMinMaxMapAccurate;

for Octave = 1:OctaveNum

    for DiffMinMaxMapIndx = 1:size(DiffMinMaxMap,2) 

        Map = DiffMinMaxMapAccurate{Octave,DiffMinMaxMapIndx};

        SSCindx = find(Map); % Scale-Space-Corner Index

        if isempty(SSCindx)

            continue;

        end

        for ssc = 1:length(SSCindx)

            [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

            if (Row <= 1) || (Row >= size(Map,1)) || (Col <= 1) || (Col >= size(Map,2))

                DiffMinMaxMapNoEdge{Octave,DiffMinMaxMapIndx}(Row,Col) = 0; % discard out of matrix boundary

                continue;

            end

            Dyy = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row+1,Col) - 2*OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col) + OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row-1,Col);

            Dxx = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col+1) - 2*OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col) + OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row,Col-1);

            Dxy = OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row-1,Col+1) - OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row-1,Col-1) - OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row+1,Col+1) + OctaveImageDiff{Octave,DiffMinMaxMapIndx+1}(Row+1,Col-1);

            TrH = Dxx + Dyy;

            DetH = Dxx*Dyy - Dxy^2;

            if ((TrH^2 / DetH) >= DetermineThreshold)

                DiffMinMaxMapNoEdge{Octave,DiffMinMaxMapIndx}(Row,Col) = 0; % discard unstable extrema

            end

        end

    end

end

%% SIFT feature descriptors generation

% the patch size for dominant orientation calculation is 9;

% the patch size for feature transformation is 16;

DominantOrientation = cell(size(DiffMinMaxMap));

SIFT = cell(size(DiffMinMaxMap));

for Octave = 1:OctaveNum

    Sigma = SigmaOrigin * 2^(Octave-1); % when up to a new octave, double the sigma

    for DiffMinMaxMapIndx = 1:size(DiffMinMaxMap,2)

        stack = DiffMinMaxMapIndx+1;

        SigmaScale = Sigma * ScaleFactor^(stack-2);

        GaussianSmoothedImage = OctaveImage{Octave,stack};

        Map = DiffMinMaxMapNoEdge{Octave,DiffMinMaxMapIndx};

        SSCindx = find(Map); % Scale-Space-Corner Index

        if isempty(SSCindx)

            continue;

        end

        DomOri = zeros(length(SSCindx),2,36); % first column is for sita, second column is for magnitude

        sift = zeros(length(SSCindx),128,36);

        for ssc = 1:length(SSCindx)

            [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

            Row = Row+1; Col = Col+1; % offset = [1,1];

            if (Row <= 10) || (Row >= size(GaussianSmoothedImage,1)-8) || (Col <= 10) || (Col >= size(GaussianSmoothedImage,2)-8)

                % skip if out of boundary

                continue;

            end

            % magnitude and sita of sample points in the patch

            mag = ((GaussianSmoothedImage(Row-8:Row+7,Col-7:Col+8) - GaussianSmoothedImage(Row-8:Row+7,Col-9:Col+6)).^2 + (GaussianSmoothedImage(Row-7:Row+8,Col-8:Col+7) - GaussianSmoothedImage(Row-9:Row+6,Col-8:Col+7)).^2).^0.5;

            sita = atan2((GaussianSmoothedImage(Row-7:Row+8,Col-8:Col+7) - GaussianSmoothedImage(Row-9:Row+6,Col-8:Col+7)),(GaussianSmoothedImage(Row-8:Row+7,Col-7:Col+8) - GaussianSmoothedImage(Row-8:Row+7,Col-9:Col+6)));

            sita = mod(sita + 2*pi, 2*pi); % the range of atan2 function is -pi~pi, map it to 0~2*pi

            % Dominant orientation calculation

            GaussianFilter = fspecial('gaussian',[9,9],1.5*SigmaScale);

            Dmag = mag(5:13,5:13).* GaussianFilter; % magnitude is weighted by gaussian filter

            sitaquantize = mod(sita(5:13,5:13) + pi/36,2*pi);

            sitaquantize = floor(sitaquantize / (2*pi/36));

            sitabin = zeros(36,1);

            for patchindx = 1:9^2

                sitabin(sitaquantize(patchindx)+1) = sitabin(sitaquantize(patchindx)+1) + Dmag(patchindx);

            end

            maxsitabin = max(sitabin);

            DominantOriBin = find((sitabin / maxsitabin) >= 0.8); % duplicate the feature when non-maximum magnitude of orientation is also big

            DominantOriSita = (DominantOriBin-1)*(2*pi/36);

            DomOri(ssc,1,1:length(DominantOriSita)) = DominantOriSita;

            DomOri(ssc,2,1:length(DominantOriSita)) = sitabin((sitabin / maxsitabin) >= 0.8);

            % SIFT feature calculation

            for DomOriIndx = 1:length(DominantOriSita)

                GaussianFilter = fspecial('gaussian',[16,16],8);

                Smag = mag.* GaussianFilter;

                sitaquantize = mod(sita - DominantOriSita(DomOriIndx) + pi/8 + 2*pi,2*pi);

                sitaquantize = floor(sitaquantize / (2*pi/8));

                sitabin = zeros(8,4,4);

                for patchindx = 1:16^2

                    [row,col] = ind2sub([16,16], patchindx);

                    row = floor((row-1)/4)+1;

                    col = floor((col-1)/4)+1;

                    sitabin(sitaquantize(patchindx)+1,row,col) = sitabin(sitaquantize(patchindx)+1,row,col) + Smag(patchindx);

                end

                sitabin = bsxfun(@times,sitabin,sum(sitabin.^2,1).^(-0.5)); % normalize the vector

                sitabin(sitabin > 0.2) = 0.2; % threshold the maximum value as 0.2

                sitabin = bsxfun(@times,sitabin,sum(sitabin.^2,1).^(-0.5)); % renormalize

                sift(ssc,:,DomOriIndx) = sitabin(:);

            end

        end

        DominantOrientation{Octave,DiffMinMaxMapIndx} = DomOri;

        SIFT{Octave,DiffMinMaxMapIndx} = sift;

    end 

end

%% Show results (this part of code is just for verification, and is a mess...)

% close all;

% imgtemp = img*0.5;

% imgtemp(2:end-1,2:end-1) = (imgtemp(2:end-1,2:end-1) + (DiffMinMaxMap{1,1}+DiffMinMaxMap{1,2}+DiffMinMaxMap{1,3}))*255;

% figure,imshow(uint8(imgtemp));

% imgtemp2 = img*0.5;

% imgtemp2(2:end-1,2:end-1) = (imgtemp2(2:end-1,2:end-1) + (DiffMinMaxMapAccurate{1,1}+DiffMinMaxMapAccurate{1,2}+DiffMinMaxMapAccurate{1,3}))*255;

% figure,imshow(uint8(imgtemp2));

% imgtemp3 = img*0.5;

% imgtemp3(2:end-1,2:end-1) = (imgtemp3(2:end-1,2:end-1) + (DiffMinMaxMapNoEdge{1,1}+DiffMinMaxMapNoEdge{1,2}+DiffMinMaxMapNoEdge{1,3}))*255;

% figure,imshow(uint8(imgtemp3));

imgtemp4 = img*0.5;

figure,imshow(uint8(imgtemp4*255));

for i = 1:3

for j = 1:3

for k = 1:2

if isempty(DominantOrientation{i,j})

continue

end

DomOri = DominantOrientation{i,j}(:,:,k);

DOx = cos(DomOri(:,1)).*DomOri(:,2)*200;

DOy = sin(DomOri(:,1)).*DomOri(:,2)*200;

Map = DiffMinMaxMapNoEdge{i,j};

SSCindx = find(Map); % Scale-Space-Corner Index

for ssc = 1:length(SSCindx)

    hold on;

    [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

    Row = Row*2^(i-1)+1; Col = Col*2^(i-1)+1; % offset = [1,1];

    if (Row > size(img,1)) || (Col > size(img,2))

        continue

    end

    quiver(Col,Row,DOx(ssc),DOy(ssc),'r')

end

end

end

end

hold off;

%

imgtemp5 = img*0.5;

for i = 1:3

for j = 1:3

for k = 1:2

if isempty(DominantOrientation{i,j})

continue

end

Map = DiffMinMaxMapNoEdge{i,j};

SSCindx = find(Map); % Scale-Space-Corner Index

for ssc = 1:length(SSCindx)

    [Row,Col] = ind2sub([size(Map,1),size(Map,2)], SSCindx(ssc));

    Row = Row*2^(i-1)+1; Col = Col*2^(i-1)+1; % offset = [1,1];

    if (Row > size(img,1)) || (Col > size(img,2))

        continue

    end

    imgtemp5(Row,Col) = 1;

end

end

end

end

figure,imshow(uint8(imgtemp5*255));