finish SIFT part

DoGdetector.m

@@ -0,0 +1,32 @@
+function [locsDoG, GaussianPyramid] = DoGdetector(im, sigma0, k, ...
+                                                levels, th_contrast, th_r)
+%%DoGdetector
+%  Putting it all together
+%
+%   Inputs          Description
+%--------------------------------------------------------------------------
+%   im              Grayscale image with range [0,1].
+%
+%   sigma0          Scale of the 0th image pyramid.
+%
+%   k               Pyramid Factor.  Suggest sqrt(2).
+%
+%   levels          Levels of pyramid to construct. Suggest -1:4.
+%
+%   th_contrast     DoG contrast threshold.  Suggest 0.03.
+%
+%   th_r            Principal Ratio threshold.  Suggest 13
+%
+%   Outputs         Description
+%--------------------------------------------------------------------------
+%
+%   locsDoG         N x 3 matrix where the DoG pyramid achieves a local extrema
+%                   in both scale and space, and satisfies the two thresholds.
+%
+%   GaussianPyramid A matrix of grayscale images of size (size(im),numel(levels))
+GaussianPyramid = createGaussianPyramid(im, sigma0, k, levels);
+[DoGPyramid, DoGLevels] = createDoGPyramid(GaussianPyramid, levels);
+PrincipalCurvature = computePrincipalCurvature(DoGPyramid);
+locsDoG = getLocalExtrema(DoGPyramid, DoGLevels, PrincipalCurvature, th_contrast, th_r);
+end
+

computeGradient.m

@@ -0,0 +1,28 @@
+function g_histogram = computeGradient(patch)
+[m,n] = size(patch);
+[FX,FY] = gradient(patch);
+g_histogram = zeros(1,8);
+for i=1:m
+    for j=1:n
+        if FX(i,j)>0 && -tan(pi/8)<(FY(i,j)/FX(i,j)) && (FY(i,j)/FX(i,j))<=tan(pi/8)
+            g_histogram(1) = g_histogram(1)+norm([FY(i,j) FX(i,j)]);
+        elseif FX(i,j)>0 && FY(i,j)>0 && tan(pi/8)<(FY(i,j)/FX(i,j)) && (FY(i,j)/FX(i,j))<=tan(3*pi/8)
+            g_histogram(2) = g_histogram(2)+norm([FY(i,j) FX(i,j)]);
+        elseif FY(i,j)>0 && tan(3*pi/8)<=abs((FY(i,j)/FX(i,j)))
+            g_histogram(3) = g_histogram(3)+norm([FY(i,j) FX(i,j)]);
+        elseif FX(i,j)<0 && FY(i,j)>0 && tan(5*pi/8)<=(FY(i,j)/FX(i,j)) && (FY(i,j)/FX(i,j))<=tan(7*pi/8)
+            g_histogram(4) = g_histogram(4)+norm([FY(i,j) FX(i,j)]);
+        elseif FX(i,j)<0 && abs((FY(i,j)/FX(i,j)))<=tan(pi/8)
+            g_histogram(5) = g_histogram(5)+norm([FY(i,j) FX(i,j)]);
+        elseif FX(i,j)<0 && FY(i,j)<0 && tan(9*pi/8)<=(FY(i,j)/FX(i,j)) && (FY(i,j)/FX(i,j))<=tan(11*pi/8)
+            g_histogram(6) = g_histogram(6)+norm([FY(i,j) FX(i,j)]);
+        elseif FY(i,j)<0 && (FY(i,j)/FX(i,j))<=tan(11*pi/8)
+            g_histogram(7) = g_histogram(7)+norm([FY(i,j) FX(i,j)]);
+        else
+            g_histogram(8) = g_histogram(8)+norm([FY(i,j) FX(i,j)]);
+        end
+    end 
+end
+
+g_histogram=g_histogram./sum(g_histogram);
+end

computePrincipalCurvature.m

@@ -0,0 +1,28 @@
+function PrincipalCurvature = computePrincipalCurvature(DoGPyramid)
+%%Edge Suppression
+% Takes in DoGPyramid generated in createDoGPyramid and returns
+% PrincipalCurvature,a matrix of the same size where each point contains the
+% curvature ratio R for the corre-sponding point in the DoG pyramid
+%
+% INPUTS
+% DoG Pyramid - size (size(im), numel(levels) - 1) matrix of the DoG pyramid
+%
+% OUTPUTS
+% PrincipalCurvature - size (size(im), numel(levels) - 1) matrix where each 
+%                      point contains the curvature ratio R for the 
+%                      corresponding point in the DoG pyramid
+[R,C,L] = size(DoGPyramid);
+PrincipalCurvature = zeros(R,C,L);
+for k=1:L
+    [Dx,Dy] = gradient(DoGPyramid(:,:,k));
+    [Dxx,Dxy] = gradient(Dx);
+    [Dyx,Dyy] = gradient(Dy); 
+    for i=1:R
+       for j=1:C
+          H = [Dxx(i,j) Dxy(i,j);Dyx(i,j) Dyy(i,j)];
+          PrincipalCurvature(i,j,k) = trace(H)^2/det(H);
+       end
+    end
+end
+
+end

computeSIFT.m

@@ -0,0 +1,60 @@
+function [locs,desc] = computeSIFT(im, GaussianPyramid, locsDoG, k, levels)
+[R,C,L] = size(GaussianPyramid);
+patchWidth = 9;
+locs = zeros(0,3);
+desc = zeros(0,128);
+norm_patch_size = 64;
+
+for i=1:size(locsDoG,1)
+    level = find(levels==locsDoG(i,3));
+    patch_half_size = floor(8*k^(1+levels(level)-levels(1)));
+    bound = floor(sqrt(2)*patch_half_size);
+    if locsDoG(i,1)>=bound && locsDoG(i,2)>=bound && locsDoG(i,1)<=C-bound && locsDoG(i,2)<=R-bound
+        locs = [locs;locsDoG(i,:)];
+        % compute 128-D discriptor
+
+        patch = GaussianPyramid(locsDoG(i,2)+1-bound:locsDoG(i,2)+bound,locsDoG(i,1)+1-bound:locsDoG(i,1)+bound,level);       
+        % compute major gradient
+        tmp_patch = GaussianPyramid(locsDoG(i,2)+1-patch_half_size:locsDoG(i,2)+patch_half_size,locsDoG(i,1)+1-patch_half_size:locsDoG(i,1)+patch_half_size,level);
+        tmp_patch = imresize(tmp_patch, [norm_patch_size norm_patch_size]);
+        g_histogram = computeGradient(tmp_patch);
+        [~,major_g_index] = max(g_histogram);
+        % rotate patch
+        patch = imrotate(patch,-45*(major_g_index-1));
+        [h,w] = size(patch);
+        h = floor(h/2);
+        w = floor(w/2);
+        patch = patch(h+1-patch_half_size:h+patch_half_size,w+1-patch_half_size:w+patch_half_size);
+        patch = imresize(patch, [norm_patch_size norm_patch_size]);
+
+        desc = [desc;zeros(1,128)];
+        for p=1:4
+            for q=1:4
+                g_histogram = computeGradient(patch(floor((p-1)*norm_patch_size/4+1):floor(p*norm_patch_size/4), floor((q-1)*norm_patch_size/4+1):floor(q*norm_patch_size/4)));
+                desc(end,((p-1)*4+q-1)*8+1:((p-1)*4+q)*8) = g_histogram;
+            end
+        end
+    end
+end
+
+% for i=1:size(locsDoG,1)
+%     level = find(levels==locsDoG(i,3));
+%     patch_half_size = 16;
+%     if locsDoG(i,1)>=patch_half_size && locsDoG(i,2)>=patch_half_size && locsDoG(i,1)<=C-patch_half_size && locsDoG(i,2)<=R-patch_half_size
+%         locs = [locs;locsDoG(i,:)];
+%         % compute 128-D discriptor
+%         
+%         patch = GaussianPyramid(locsDoG(i,2)+1-patch_half_size:locsDoG(i,2)+patch_half_size,locsDoG(i,1)+1-patch_half_size:locsDoG(i,1)+patch_half_size,level);       
+%         patch = imresize(patch, [norm_patch_size norm_patch_size]);
+%         
+%         desc = [desc;zeros(1,128)];
+%         for p=1:4
+%             for q=1:4
+%                 g_histogram = computeGradient(patch(floor((p-1)*norm_patch_size/4+1):floor(p*norm_patch_size/4), floor((q-1)*norm_patch_size/4+1):floor(q*norm_patch_size/4)));
+%                 desc(end,((p-1)*4+q-1)*8+1:((p-1)*4+q)*8) = g_histogram;
+%             end
+%         end
+%     end
+% end
+
+end

createDoGPyramid.m

@@ -0,0 +1,17 @@
+function [DoGPyramid, DoGLevels] = createDoGPyramid(GaussianPyramid, levels)
+%%Produces DoG Pyramid
+% inputs
+% Gaussian Pyramid - A matrix of grayscale images of size
+%                    (size(im), numel(levels))
+% levels      - the levels of the pyramid where the blur at each level is
+%               outputs
+% DoG Pyramid - size (size(im), numel(levels) - 1) matrix of the DoG pyramid
+%               created by differencing the Gaussian Pyramid input
+[R,C,L] = size(GaussianPyramid);
+DoGLevels = levels(2:end);
+DoGPyramid = zeros(R,C,L-1);
+for i=1:L-1
+    DoGPyramid(:,:,i) = GaussianPyramid(:,:,i+1)-GaussianPyramid(:,:,i);
+end
+
+end

createGaussianPyramid.m

@@ -0,0 +1,23 @@
+function [GaussianPyramid] = createGaussianPyramid(im, sigma0, k, levels)
+% function [GaussianPyramid] = createGaussianPyramid(im, sigma0, k, levels)
+%%Produces Gaussian Pyramid
+% inputs
+% im - a grayscale image with range 0 to 1
+% sigma0 - the standard deviation of the blur at level 0
+% k - the multiplicative factor of sigma at each level, where sigma=sigma_0 k^l
+% levels - the levels of the pyramid where the blur at each level is
+% sigma=sigma0 k^l
+% outputs
+% A matrix of grayscale images of size (size(im),numel(levels))
+
+im = im2double(im);
+if size(im,3)==3
+    im= rgb2gray(im);
+end
+
+GaussianPyramid = zeros([size(im),length(levels)]);
+for i = 1:length(levels)
+    sigma_ = sigma0*k^levels(i);
+    h = fspecial('gaussian',floor(3*sigma_*2)+1,sigma_);
+    GaussianPyramid(:,:,i) = imfilter(im,h);
+end

getLocalExtrema.m

@@ -0,0 +1,62 @@
+function locsDoG = getLocalExtrema(DoGPyramid, DoGLevels, ...
+                        PrincipalCurvature, th_contrast, th_r)
+%%Detecting Extrema
+% INPUTS
+% DoG Pyramid - size (size(im), numel(levels) - 1) matrix of the DoG pyramid
+% DoG Levels  - The levels of the pyramid where the blur at each level is
+%               outputs
+% PrincipalCurvature - size (size(im), numel(levels) - 1) matrix contains the
+%                      curvature ratio R
+% th_contrast - remove any point that is a local extremum but does not have a
+%               DoG response magnitude above this threshold
+% th_r        - remove any edge-like points that have too large a principal
+%               curvature ratio
+%
+% OUTPUTS
+% locsDoG - N x 3 matrix where the DoG pyramid achieves a local extrema in both
+%           scale and space, and also satisfies the two thresholds.
+[R,C,L] = size(DoGPyramid);
+locsDoG = zeros(0,3);
+for i=1:R
+    for j=1:C
+        for k=1:L
+            isExtrema = 1;            
+            if PrincipalCurvature(i,j,k)<th_r && PrincipalCurvature(i,j,k)>0 && abs(DoGPyramid(i,j,k))>th_contrast            
+                localmax = max(max(DoGPyramid(max(i-1,1):min(i+1,R),max(j-1,1):min(j+1,C),k)));
+                localmin = min(min(DoGPyramid(max(i-1,1):min(i+1,R),max(j-1,1):min(j+1,C),k)));
+                if DoGPyramid(i,j,k)~=localmax && DoGPyramid(i,j,k)~=localmin
+                    continue
+                end
+                if k>1
+                    if DoGPyramid(i,j,k)==localmin
+                        if DoGPyramid(i,j,k)>DoGPyramid(i,j,k-1)
+                            isExtrema = 0;
+                        end
+
+                    else
+                        if DoGPyramid(i,j,k)<DoGPyramid(i,j,k-1)
+                            isExtrema = 0;
+                        end
+                    end
+                end
+                if k<L
+                    if DoGPyramid(i,j,k)==localmin
+                        if DoGPyramid(i,j,k)>DoGPyramid(i,j,k+1)
+                            isExtrema = 0;
+                        end
+                    else
+                        if DoGPyramid(i,j,k)<DoGPyramid(i,j,k+1)
+                            isExtrema = 0;
+                        end
+                    end
+                end
+
+                if isExtrema
+                    locsDoG = [locsDoG;j i DoGLevels(k)];
+                end
+            end            
+        end
+    end    
+end
+
+end

plotMatches.m

@@ -0,0 +1,31 @@
+function [] = plotMatches(im1, im2, matches, locs1, locs2)
+% INPUT
+% im1, im2 - grayscale images, scaled into [0 1]
+% matches -  the list of matches returned by briefMatch()
+% locs1, locs 2 - the locations of keypoints returne by briefLite() . 
+%                 m1 x 3 and m2 x3 respectively. each row corresponds to the 
+%                 space-scale position: [x_coordinate,y_coordiante,scale]
+
+
+    im1= im2double(im1);
+    im2= im2double(im2);
+
+    width = size(im1,2) + size(im2,2);
+    height = max(size(im1,1), size(im2,1));
+    nchannel = size(im1,3); 
+    img = zeros(height, width,nchannel);
+%     size(img)
+    img(1:size(im1,1),1:size(im1,2),:) = im1;
+    img(1:size(im2,1),size(im1,2)+1:size(im1,2) + size(im2,2),:) = im2; 
+    imshow(img);
+
+    axis equal;
+    hold on;
+    for i = 1:size(matches,1)
+        p1 = locs1(matches(i,1),:);
+        p2 = locs2(matches(i,2),:); 
+        line([p1(1) size(im1,2)+p2(1)],[p1(2),p2(2)], 'Color','r','LineWidth',1); 
+    end
+    hold off; 
+
+end

showInterestPoint.m

@@ -0,0 +1,14 @@
+im = imread('..\data\chickenbroth_01.jpg');
+im = im2double(im);
+if size(im,3)==3
+    im= rgb2gray(im);
+end
+levels = [-1; 0; 1; 2; 3; 4];
+sigma0 = 1;
+k = sqrt(2);
+th_contrast = 0.03;
+th_r = 12;
+[locsDoG, GaussianPyramid] = DoGdetector(im, sigma0, k, levels, th_contrast, th_r);
+imshow(im);
+hold on;
+plot(locsDoG(:,1),locsDoG(:,2),'.','MarkerEdgeColor','g','MarkerSize',10);

siftLite.m

@@ -0,0 +1,18 @@
+function [locs, desc] = siftLite(im)
+% INPUTS
+% im - gray image with values between 0 and 1
+%
+% OUTPUTS
+% locs - an m x 3 vector, where the first two columns are the image coordinates 
+% 		 of keypoints and the third column is the pyramid level of the keypoints
+% desc - an m x n bits matrix of stacked BRIEF descriptors. 
+%		 m is the number of valid descriptors in the image and will vary
+% 		 n is the number of bits for the BRIEF descriptor
+levels = [-1; 0; 1; 2; 3; 4];
+sigma0 = 1;
+k = sqrt(2);
+th_contrast = 0.03;
+th_r = 12;
+[locsDoG, GaussianPyramid] = DoGdetector(im, sigma0, k, levels, th_contrast, th_r);
+[locs,desc] = computeSIFT(im, GaussianPyramid, locsDoG, k, levels);
+end

siftMatch.m

@@ -0,0 +1,26 @@
+function [matches] = siftMatch(desc1, desc2, ratio)
+%function [matches] = briefMatch(desc1, desc2, ratio)
+% Performs the descriptor matching
+% inputs  : desc1 , desc2 - m1 x n and m2 x n matrix. m1 and m2 are the number of keypoints in image 1 and 2.
+%						    n is the number of bits in the brief
+% outputs : matches - p x 2 matrix. where the first column are indices
+%									into desc1 and the second column are indices into desc2
+
+if nargin<3
+    ratio = .8;
+end
+
+% compute the pairwise Hamming distance
+[D,I] = pdist2(desc2,desc1,'cityblock','Smallest',2);
+ix = 1:size(desc1,1);
+
+% suprress match between descriptors that are not distriminative.
+r = D(1,:)./D(2,:); 
+ix = ix(r < ratio | isnan(r)); 
+I2 = I(1,r < ratio | isnan(r));
+
+%output
+matches = [ix' I2'];
+
+
+end

testMatch.m

@@ -0,0 +1,21 @@
+im1 = imread('../data/model_chickenbroth.jpg');
+% im1 = imread('../data/incline_L.png');
+% im1 = imread('../data/pf_scan_scaled.jpg');
+im1 = im2double(im1);
+if size(im1,3)==3
+    im1= rgb2gray(im1);
+end
+[locs1, desc1] = siftLite(im1);
+
+im2 = imread('../data/chickenbroth_01.jpg');
+% im2 = imrotate(im2,180);
+% im2 = imread('../data/incline_R.png');
+% im2 = imread('../data/pf_stand.jpg');
+im2 = im2double(im2);
+if size(im2,3)==3
+    im2= rgb2gray(im2);
+end
+[locs2, desc2] = siftLite(im2);
+
+[matches] = siftMatch(desc1, desc2);
+plotMatches(im1, im2, matches, locs1, locs2)