Project 2: Project 2: Local Feature Matching

Result of a hybrid image.

Objective

The goal of this assignment is to understand a local feature matching algorithm using techniques.

The assignment contains several parts:

  1. Interest point detection.
  2. Local feature description.
  3. Feature matching.

Methods and algorithm

1. Interest point detection

In this assignment, we are taught to use harris corner detector to choose interested points. The picture could be divided into 3 regions: flat, edge and corner, and it is corners that are most interested. Here we use second moment matrix to find corners.

2. Local feature description

In this assignment, the feature width is set to 16X16 , and divded to 16 4X4 cell which detects 8 orientations. So that each discriptor contains 4X4X8 dimentsion histogram that discribe the gradient of pixels around the interested point. Since the connected components will generate severl interested point, we need to find the local maxima in the connected components. In this assignment, function bwconnecomp is used to find connected components and generate one local maxima interested point.

3. Feature matching

In this assignment, distance between two features of different interested point pair is calculated to eliminate wrong pairing caused by having similar pattern. Therefore pair with least distance will be chosen.

Explanation of codes

Interest point detection


    [m,n,L] = size(image);
    image = double(image);
    g1 = fspecial('gaussian',[7,7], 1);
    g2 = fspecial('gaussian',[8,8],0.1);
    [dx,dy]=gradient(g1);
%%Build second moment matrix
    for i = 1:L
        Ix(:,:,i) = imfilter(image(:,:,i),dx); 
        Iy(:,:,i) = imfilter(image(:,:,i),dy);    
        Ix2(:,:,i) = imfilter(Ix(:,:,i).^2, g2);
        Iy2(:,:,i) = imfilter(Iy(:,:,i).^2, g2);
        Ixy(:,:,i) = imfilter(Ix(:,:,i).*Iy(:,:,i), g2);
    end
%%Measure corner response    
    alpha = 0.05;
    cim = (Ix2.*Iy2 - Ixy.^2) - alpha*(Ix2 + Iy2).^2;     
    thresh = 0.000001;
%%Transform to binary image
   	cimMark = cim>thresh;
    list=[];

%%Use bwconnecomp function to generate local maxima
    for k = 1:L
        connect(k)=bwconncomp(cimMark(:,:,k),8);
        cim2 = cim(:,:,k);
        [height,length]=size(connect(k).PixelIdxList);
        for i = 1:length
            cell = cell2mat(connect(k).PixelIdxList(i));
            a = -Inf;
            b = 1;
            for j = 1:size(cell)
                if cim2(cell(j))>a
                    a = cim2(cell(j));
                    b = cell(j);
                end
            end
            if ~any(list==b)
                list = [list b];
            end    
        end
    end
    y = ceil(list./m);
    x = mod(list,m);
    delMask = (ceil(x - (feature_width / 2)) < 1 | ceil(x + (feature_width / 2 - 1)) > m | ceil(y - feature_width / 2) < 1 | ceil(y + feature_width / 2 - 1) > n);
    x(delMask) = [];
    y(delMask) = [];  

Local feature description


x = round(x);
y = round(y);
[M,N,L]=size(image);
g1 = fspecial('gaussian',[feature_width,feature_width],1);
[dx, dy] = gradient(g1);
grad_imX = imfilter(image,dx);
grad_imY = imfilter(image,dy);
n = length(x);
features = zeros(n,128);
for k = 1:n
    his_int{k} = cell(4);
    for i = 1:16
        his_int{k}{i}=zeros(1,8);
    end
%%Generate gradients
    for l = 1:L
        grad_interestX = grad_imX(round(x(k)-feature_width/2):round(x(k)+feature_width/2-1),(round(y(k)-feature_width/2):round(y(k)+feature_width/2-1)));
        grad_interestY = grad_imY(round(x(k)-feature_width/2):round(x(k)+feature_width/2-1),(round(y(k)-feature_width/2):round(y(k)+feature_width/2-1)));
        Mag = sqrt(grad_interestX.^2+grad_interestY.^2);
        [m_grad,n_grad,L_grad]=size(grad_interestX);

%% Assign magnitude to different orientation according to the angle. For convinience, magnitude will be add to both orientation if the angle is located between them.
        for i = 1:m_grad
            for j = 1:n_grad
                if grad_interestY(i,j)==0&&grad_interestX==0
                    bins = [1 5];
                end
                if grad_interestX(i,j)==0&&grad_interestY(i,j)>0
                    bins = [3];
                end
                if grad_interestX(i,j)==0&&grad_interestY(i,j)<0
                    bins = [7];
                end
                if grad_interestX(i,j)>0
                    if grad_interestY(i,j)/grad_interestX(i,j) ==0
                        bins = [1];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)>0 && grad_interestY(i,j)/grad_interestX(i,j) <1
                        bins = [1 2];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)==1
                        bins = [2];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)>1
                        bins = [2 3];              
                    elseif grad_interestY(i,j)/grad_interestX(i,j)<-1
                        bins = [7 8];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)==-1
                        bins = [8];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)>-1 && grad_interestY(i,j)/grad_interestX(i,j)<0
                        bins = [1 8];
                    end
                end
                if grad_interestX(i,j)<0
                    if grad_interestY(i,j)/grad_interestX(i,j) ==0
                        bins = [5];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)>0 && grad_interestY(i,j)/grad_interestX(i,j) <1
                        bins = [5 6];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)==1
                        bins = [6];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)>1
                        bins = [6 7];              
                    elseif grad_interestY(i,j)/grad_interestX(i,j)<-1
                        bins = [3 4];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)==-1
                        bins = [4];
                    elseif grad_interestY(i,j)/grad_interestX(i,j)>-1 && grad_interestY(i,j)/grad_interestX(i,j)<0
                        bins = [4 5];
                    end
                end
%%Assign magnitude 
                his_int{k}{ceil(i*4/feature_width),ceil(j*4/feature_width)}(bins) = his_int{k}{ceil(i*4/feature_width),ceil(j*4/feature_width)}(bins)+Mag(i,j);                   
            end
        end
    end
    for und = 1 : 16
        features(k, (und - 1) * 8 + 1 : und * 8) = his_int{k}{und};
    end
    features(k, :) = features(k, :) .* (1 / norm(features(k, :)));
end

Feature matching


[num,area]=size(features1);
[num2,area2]=size(features2);
matches = zeros(num,2);
confidences = zeros(num,1);
for i = 1:num
    dis = Inf;
    dis_2 = Inf;
    point1 = 1;
    for j = 1:num2
        L = sqrt(sum((features1(i,:)-features2(j,:)).^2));
        if L < dis
            point1 = j;
            dis_2 = dis;
            dis = L;
        end
    end
    matches(i,1) = i;
    matches(i,2) = point1;
    confidences(i)=(dis_2/dis);    
end

Notre Dame Accurancy=80% with 200 points


Mount Rushmore Accurancy=66%


Episcopal Gaudi Accurancy=4%