鱼眼拼接之SIFT提取特征点

之前写的那两篇 是我没有完全理解SIFT理论部分 后来看了这位女侠的http://blog.csdn.net/abcjennifer/article/details/7639681终于弄明白了  


%这部分是之前写过的鱼眼提取有效区域及经度坐标校正
A=imread('F:\orl_zhifangtu\s1.jpg');
T=40;
[A,R]=kuaisusaomiao(A,T);
[m,n,H]=size(A);
C=[];
x=n/2;
y=m/2;
for u=1:m
    for v=1:n
        i=u;
        j=round(sqrt(R^2-(y-u)^2)*(v-x)/R+x);
        if(R^2-(y-u)^2<0)
            continue;
        end
        C(u,v,1)=A(i,j,1);
        C(u,v,2)=A(i,j,2);
        C(u,v,3)=A(i,j,3);
    end
end

C=uint8(C);

这是原图

这是校正然后灰度化以后的图

接下来是SIFT高斯卷积 这里之前我用的http://blog.csdn.net/v_JULY_v/article/details/6245939这个理论部分 按照他的高斯函数sigma第一层第一幅图取值是1.6 第二幅图是k*1.6 第三幅是K^2*1.6 第四幅是K^3*1.6 第五幅是K^4*1.6（第一层五幅图的高斯函数已经完毕 接下来第二层）;第二层第一幅图是2*1.6 第二幅图是2*K^2*1.6（是上一层第二幅的两倍） 第三幅是2*K^3*1.6 以此类推。。。第三层第一幅图是2*2*1.6 第二幅图是2*2*K*1.6（即是上一层的两倍）以此类推。。。。但完完全全按照这个大神编  后面会有问题  就是比较极值显示时 会一个点都没有 （我忘记了当时是一个点都没有还是只有肉眼可见的零星的几个点 反正不是我们要的特征点） 所以后来看这个女侠解说的理论 尺度连续性那里总算看懂了  但是编的时候如果完全按照那个参数sigma编  还是会有一样的问题   
7，0.7
所以我总共编了三层金字塔，每一层有六幅图像，每一次的sigma不是按照k变化的  而是随意取的0.3，0.4，0.5，0.6，0.7，0.8这六个参数  可得到六个高斯函数 然后对原图进行卷积 得到 第一层高斯卷积后的图


%先灰度化
C=rgb2gray(C);
C=double(C);
[m1,n1,h1]=size(C);
k=2^(1/3);
threshold=3;
h11=fspecial('gaussian',[5 5],0.3);
C11=imfilter(C,h11,'conv');
h12=fspecial('gaussian',[5 5],0.4);
C12=imfilter(C,h12,'conv');
h13=fspecial('gaussian',[5 5],0.5);
C13=imfilter(C,h13,'conv');
h14=fspecial('gaussian',[5 5],0.6);
C14=imfilter(C,h14,'conv');
h15=fspecial('gaussian',[5 5],0.7);
C15=imfilter(C,h15,'conv');
h16=fspecial('gaussian',[5 5],0.8);
C16=imfilter(C,h16,'conv');
%以上就得到了第一层的六幅高斯卷积后的图C11，C12，C13,C14，C15，C16
然后两两相减 得到第一层的DOG差分空间后的图  只有五幅 D11，D12，D13，D14，D15
D11=C11-C12;
D12=C12-C13;
D13=C13-C14;
D14=C14-C15;
D15=C15-C16;
%接下来就是对第一层进行极值比较 初步得到候选的特征点 只有D12 D13 D14前后都有图 所以才可进行比较 首尾即D11，D15没办法得到26个特征点 就不能进行比较 所以我把D12 D13 D14这三幅图得到的特征点 分别装在E11 E12 E13里 
E11=zeros(m1,n1);
for i=3:m1-2
    for j=3:n1-2
        if(D12(i,j)>D12(i-1,j-1) && D12(i,j)>D12(i,j-1) && D12(i,j)>D12(i+1,j-1) && D12(i,j)>D12(i-1,j) && D12(i,j)>D12(i+1,j) && D12(i,j)>D12(i-1,j+1) && D12(i,j)>D12(i,j+1) && D12(i,j)>D12(i+1,j+1))
            if(D12(i,j)>D11(i,j) && D12(i,j)>D11(i-1,j-1) && D12(i,j)>D11(i,j-1) && D12(i,j)>D11(i+1,j-1) && D12(i,j)>D11(i-1,j) && D12(i,j)>D11(i+1,j) && D12(i,j)>D11(i-1,j+1) && D12(i,j)>D11(i,j+1) && D12(i,j)>D11(i+1,j+1))
                if(D12(i,j)>D13(i,j) && D12(i,j)>D13(i-1,j-1) && D12(i,j)>D13(i,j-1) && D12(i,j)>D13(i+1,j-1) && D12(i,j)>D13(i-1,j) && D12(i,j)>D13(i+1,j) && D12(i,j)>D13(i-1,j+1) && D12(i,j)>D13(i,j+1) && D12(i,j)>D13(i+1,j+1))
                    if(D12(i,j)>threshold)
                         E11(i,j)=1;
                    end
                end
            end
        elseif(D12(i,j)<D12(i-1,j-1) && D12(i,j)<D12(i,j-1) && D12(i,j)<D12(i+1,j-1) && D12(i,j)<D12(i-1,j) && D12(i,j)<D12(i+1,j) && D12(i,j)<D12(i-1,j+1) && D12(i,j)<D12(i,j+1) && D12(i,j)<D12(i+1,j+1))
            if(D12(i,j)<D11(i,j) && D12(i,j)<D11(i-1,j-1) && D12(i,j)<D11(i,j-1) && D12(i,j)<D11(i+1,j-1) && D12(i,j)<D11(i-1,j) && D12(i,j)<D11(i+1,j) && D12(i,j)<D11(i-1,j+1) && D12(i,j)<D11(i,j+1) && D12(i,j)<D11(i+1,j+1))
                if(D12(i,j)<D13(i,j) && D12(i,j)<D13(i-1,j-1) && D12(i,j)<D13(i,j-1) && D12(i,j)<D13(i+1,j-1) && D12(i,j)<D13(i-1,j) && D12(i,j)<D13(i+1,j) && D12(i,j)<D13(i-1,j+1) && D12(i,j)<D13(i,j+1) && D12(i,j)<D13(i+1,j+1))
                    if(D12(i,j)<-threshold)
                        E11(i,j)=-1;
                    end
                end
            end
        end
    end
end
E12=zeros(m1,n1);
for i=2:m1-1
    for j=2:n1-1
        if(D13(i,j)>D13(i-1,j-1) && D13(i,j)>D13(i,j-1) && D13(i,j)>D13(i+1,j-1) && D13(i,j)>D13(i-1,j) && D13(i,j)>D13(i+1,j) && D13(i,j)>D13(i-1,j+1) && D13(i,j)>D12(i,j+1) && D13(i,j)>D13(i+1,j+1))
            if(D13(i,j)>D12(i,j) && D13(i,j)>D12(i-1,j-1) && D13(i,j)>D12(i,j-1) && D13(i,j)>D12(i+1,j-1) && D13(i,j)>D12(i-1,j) && D13(i,j)>D12(i+1,j) && D13(i,j)>D12(i-1,j+1) && D13(i,j)>D12(i,j+1) && D13(i,j)>D12(i+1,j+1))
                if(D13(i,j)>D14(i,j) && D13(i,j)>D14(i-1,j-1) && D13(i,j)>D14(i,j-1) && D13(i,j)>D14(i+1,j-1) && D13(i,j)>D14(i-1,j) && D13(i,j)>D14(i+1,j) && D13(i,j)>D14(i-1,j+1) && D13(i,j)>D14(i,j+1) && D13(i,j)>D14(i+1,j+1))
                    if(D13(i,j)>threshold)
                      E12(i,j)=1;
                    end
                end
            end
        elseif(D13(i,j)<D13(i-1,j-1) && D13(i,j)<D13(i,j-1) && D13(i,j)<D13(i+1,j-1) && D13(i,j)<D13(i-1,j) && D13(i,j)<D13(i+1,j) && D13(i,j)<D13(i-1,j+1) && D13(i,j)<D12(i,j+1) && D13(i,j)<D13(i+1,j+1))
            if(D13(i,j)<D12(i,j) && D13(i,j)<D12(i-1,j-1) && D13(i,j)<D12(i,j-1) && D13(i,j)<D12(i+1,j-1) && D13(i,j)<D12(i-1,j) && D13(i,j)<D12(i+1,j) && D13(i,j)<D12(i-1,j+1) && D13(i,j)<D12(i,j+1) && D13(i,j)<D12(i+1,j+1))
                if(D13(i,j)<D14(i,j) && D13(i,j)<D14(i-1,j-1) && D13(i,j)<D14(i,j-1) && D13(i,j)<D14(i+1,j-1) && D13(i,j)<D14(i-1,j) && D13(i,j)<D14(i+1,j) && D13(i,j)<D14(i-1,j+1) && D13(i,j)<D14(i,j+1) && D13(i,j)<D14(i+1,j+1))
                    if(D13(i,j)<-threshold)
                      E12(i,j)=-1;
                    end
                end
            end
        end
    end
end
E13=zeros(m1,n1);
for i=2:m1-1
    for j=2:n1-1
        if(D14(i,j)>D14(i-1,j-1) && D14(i,j)>D14(i,j-1) && D14(i,j)>D14(i+1,j-1) && D14(i,j)>D14(i-1,j) && D14(i,j)>D14(i+1,j) && D14(i,j)>D14(i-1,j+1) && D14(i,j)>D14(i,j+1) && D14(i,j)>D14(i+1,j+1))
            if(D14(i,j)>D13(i,j) && D14(i,j)>D13(i-1,j-1) && D14(i,j)>D13(i,j-1) && D14(i,j)>D13(i+1,j-1) && D14(i,j)>D13(i-1,j) && D14(i,j)>D13(i+1,j) && D14(i,j)>D13(i-1,j+1) && D14(i,j)>D13(i,j+1) && D14(i,j)>D13(i+1,j+1))
                if(D14(i,j)>D15(i,j) && D14(i,j)>D15(i-1,j-1) && D14(i,j)>D15(i,j-1) && D14(i,j)>D15(i+1,j-1) && D14(i,j)>D15(i-1,j) && D14(i,j)>D15(i+1,j) && D14(i,j)>D15(i-1,j+1) && D14(i,j)>D15(i,j+1) && D14(i,j)>D15(i+1,j+1))
                    if(D14(i,j)>threshold)
                        E13(i,j)=1;
                    end
                end
            end
        elseif(D14(i,j)<D14(i-1,j-1) && D14(i,j)<D14(i,j-1) && D14(i,j)<D14(i+1,j-1) && D14(i,j)<D14(i-1,j) && D14(i,j)<D14(i+1,j) && D14(i,j)<D14(i-1,j+1) && D14(i,j)<D14(i,j+1) && D14(i,j)<D14(i+1,j+1))
            if(D14(i,j)<D13(i,j) && D14(i,j)<D13(i-1,j-1) && D14(i,j)<D13(i,j-1) && D14(i,j)<D13(i+1,j-1) && D14(i,j)<D13(i-1,j) && D14(i,j)<D13(i+1,j) && D14(i,j)<D13(i-1,j+1) && D14(i,j)<D13(i,j+1) && D14(i,j)<D13(i+1,j+1))
                if(D14(i,j)<D15(i,j) && D14(i,j)<D15(i-1,j-1) && D14(i,j)<D15(i,j-1) && D14(i,j)<D15(i+1,j-1) && D14(i,j)<D15(i-1,j) && D14(i,j)<D15(i+1,j) && D14(i,j)<D15(i-1,j+1) && D14(i,j)<D15(i,j+1) && D14(i,j)<D15(i+1,j+1))
                    if(D14(i,j)>threshold)
                       E13(i,j)=-1;
                    end
                end
            end
        end
    end
end

这样就得到了初步得到了第一层的候选特征点

figure
>> imshow(D11,[0 1])
>> figure
>> imshow(D12,[0 1])
>> figure
>> imshow(D13,[0 1])
>> figure
>> imshow(D14,[0 1])
>> figure
>> imshow(D15,[0 1])

这是D11

这是D12

这是D13

这是D14

这是D15

%下面显示E11 即极大极小值 特征点

figure
>> imshow(E11,[-1 1])
>> figure
>> imshow(E12,[-1 1])
>> figure
>> imshow(E13,[-1 1])

白色表示极大值 黑色表示极小值

但是E12 E13这两幅 出来就看不到特征点了 我在想 是不是 一层都只有一幅有特征点的图出来 后面的两幅图都看不到特征点 

下面是E12 和 E13

可以看到 两幅都没有特征点出来 好奇怪明明都是26个像素比较 容我想想  不过既然只有E11有特征点出来  我们就暂时第一层就用E11吧  至少出来了一个

第二层 第三层跟上面步骤道理一样：

m2=round(m1/2);
n2=round(n1/2);
C20=imresize(C11,[m2 n2]);
h21=fspecial('gaussian',[5 5],k*0.3);
C21=imfilter(C20,h21,'conv');
h22=fspecial('gaussian',[5 5],k*0.4);
C22=imfilter(C20,h22,'conv');
h23=fspecial('gaussian',[5 5],k*0.5);
C23=imfilter(C20,h23,'conv');
h24=fspecial('gaussian',[5 5],k*0.6);
C24=imfilter(C20,h24,'conv');
h25=fspecial('gaussian',[5 5],k*0.7);
C25=imfilter(C20,h25,'conv');
h26=fspecial('gaussian',[5 5],k*0.8);
C26=imfilter(C20,h26,'conv');
m3=round(m2/2);
n3=round(n2/2);
C30=imresize(C21,[m3 n3]);
h31=fspecial('gaussian',[5 5],k^2*0.3);
C31=imfilter(C30,h31,'conv');
h32=fspecial('gaussian',[5 5],k^2*0.4);
C32=imfilter(C30,h32,'conv');
h33=fspecial('gaussian',[5 5],k^2*0.5);
C33=imfilter(C30,h33,'conv');
h34=fspecial('gaussian',[5 5],k^2*0.6);
C34=imfilter(C30,h34,'conv');
h35=fspecial('gaussian',[5 5],k^2*0.7);
C35=imfilter(C30,h35,'conv');
h36=fspecial('gaussian',[5 5],k^2*0.8);
C36=imfilter(C30,h36,'conv');
>> D21=C21-C22;
D22=C22-C23;
D23=C23-C24;
D24=C24-C25;
D25=C25-C26;
D31=C31-C32;
D32=C32-C33;
D33=C33-C34;
D34=C34-C35;
D35=C35-C36;
>> figure
>> imshow(D21,[0 1])

这是第二层的第一幅DOG高斯尺度空间的图 依次可以显示第二层第三层的所有图。。。

比较第二层的极值点

E21=zeros(m2,n2);
for i=2:m2-1
    for j=2:n2-1
        if(D22(i,j)>D22(i-1,j-1) && D22(i,j)>D22(i,j-1) && D22(i,j)>D22(i+1,j-1) && D22(i,j)>D22(i-1,j) && D22(i,j)>D22(i+1,j) && D22(i,j)>D22(i-1,j+1) && D22(i,j)>D22(i,j+1) && D22(i,j)>D22(i+1,j+1))
            if(D22(i,j)>D21(i,j) && D22(i,j)>D21(i-1,j-1) && D22(i,j)>D21(i,j-1) && D22(i,j)>D21(i+1,j-1) && D22(i,j)>D21(i-1,j) && D22(i,j)>D21(i+1,j) && D22(i,j)>D21(i-1,j+1) && D22(i,j)>D21(i,j+1) && D22(i,j)>D21(i+1,j+1))
                if(D22(i,j)>D23(i,j) && D22(i,j)>D23(i-1,j-1) && D22(i,j)>D23(i,j-1) && D22(i,j)>D23(i+1,j-1) && D22(i,j)>D23(i-1,j) && D22(i,j)>D23(i+1,j) && D22(i,j)>D23(i-1,j+1) && D22(i,j)>D23(i,j+1) && D22(i,j)>D23(i+1,j+1))
                    if(D22(i,j)>threshold)
                       E21(i,j)=1;
                    end
                end
            end
        elseif(D22(i,j)<D22(i-1,j-1) && D22(i,j)<D22(i,j-1) && D22(i,j)<D22(i+1,j-1) && D22(i,j)<D22(i-1,j) && D22(i,j)<D22(i+1,j) && D22(i,j)<D22(i-1,j+1) && D22(i,j)<D22(i,j+1) && D22(i,j)<D22(i+1,j+1))
            if(D22(i,j)<D21(i,j) && D22(i,j)<D21(i-1,j-1) && D22(i,j)<D21(i,j-1) && D22(i,j)<D21(i+1,j-1) && D22(i,j)<D21(i-1,j) && D22(i,j)<D21(i+1,j) && D22(i,j)<D21(i-1,j+1) && D22(i,j)<D21(i,j+1) && D22(i,j)<D21(i+1,j+1))
                if(D22(i,j)<D23(i,j) && D22(i,j)<D23(i-1,j-1) && D22(i,j)<D23(i,j-1) && D22(i,j)<D23(i+1,j-1) && D22(i,j)<D23(i-1,j) && D22(i,j)<D23(i+1,j) && D22(i,j)<D23(i-1,j+1) && D22(i,j)<D23(i,j+1) && D22(i,j)<D23(i+1,j+1))
                    if(D22(i,j)<-threshold)
                       E21(i,j)=-1;
                    end
                end
            end
        end
    end
end
E22=zeros(m2,n2);
for i=2:m2-1
    for j=2:n2-1
        if(D23(i,j)>D23(i-1,j-1) && D23(i,j)>D23(i,j-1) && D23(i,j)>D23(i+1,j-1) && D23(i,j)>D23(i-1,j) && D23(i,j)>D23(i+1,j) && D23(i,j)>D23(i-1,j+1) && D23(i,j)>D23(i,j+1) && D23(i,j)>D23(i+1,j+1))
            if(D23(i,j)>D22(i,j) && D23(i,j)>D22(i-1,j-1) && D23(i,j)>D22(i,j-1) && D23(i,j)>D22(i+1,j-1) && D23(i,j)>D22(i-1,j) && D23(i,j)>D22(i+1,j) && D23(i,j)>D22(i-1,j+1) && D23(i,j)>D22(i,j+1) && D23(i,j)>D22(i+1,j+1))
                if(D23(i,j)>D24(i,j) && D23(i,j)>D24(i-1,j-1) && D23(i,j)>D24(i,j-1) && D23(i,j)>D24(i+1,j-1) && D23(i,j)>D24(i-1,j) && D23(i,j)>D24(i+1,j) && D23(i,j)>D24(i-1,j+1) && D23(i,j)>D24(i,j+1) && D23(i,j)>D24(i+1,j+1))
                    if(D23(i,j)>threshold)
                      E22(i,j)=1;
                    end
                end
            end
        elseif(D23(i,j)<D23(i-1,j-1) && D23(i,j)<D23(i,j-1) && D23(i,j)<D23(i+1,j-1) && D23(i,j)<D23(i-1,j) && D23(i,j)<D23(i+1,j) && D23(i,j)<D23(i-1,j+1) && D23(i,j)<D23(i,j+1) && D23(i,j)<D23(i+1,j+1))
            if(D23(i,j)<D22(i,j) && D23(i,j)<D22(i-1,j-1) && D23(i,j)<D22(i,j-1) && D23(i,j)<D22(i+1,j-1) && D23(i,j)<D22(i-1,j) && D23(i,j)<D22(i+1,j) && D23(i,j)<D22(i-1,j+1) && D23(i,j)<D22(i,j+1) && D23(i,j)<D22(i+1,j+1))
                if(D23(i,j)<D24(i,j) && D23(i,j)<D24(i-1,j-1) && D23(i,j)<D24(i,j-1) && D23(i,j)<D24(i+1,j-1) && D23(i,j)<D24(i-1,j) && D23(i,j)<D24(i+1,j) && D23(i,j)<D24(i-1,j+1) && D23(i,j)<D24(i,j+1) && D23(i,j)<D24(i+1,j+1))
                    if(D23(i,j)<-threshold)
                       E22(i,j)=-1;
                    end
                end
            end
        end
    end
end
E23=zeros(m2,n2);
for i=2:m2-1
    for j=2:n2-1
        if(D24(i,j)>D24(i-1,j-1) && D24(i,j)>D24(i,j-1) && D24(i,j)>D24(i+1,j-1) && D24(i,j)>D24(i-1,j) && D24(i,j)>D24(i+1,j) && D24(i,j)>D24(i-1,j+1) && D24(i,j)>D24(i,j+1) && D24(i,j)>D24(i+1,j+1))
            if(D24(i,j)>D23(i,j) && D24(i,j)>D23(i-1,j-1) && D24(i,j)>D23(i,j-1) && D24(i,j)>D23(i+1,j-1) && D24(i,j)>D23(i-1,j) && D24(i,j)>D23(i+1,j) && D24(i,j)>D23(i-1,j+1) && D24(i,j)>D23(i,j+1) && D24(i,j)>D23(i+1,j+1))
                if(D24(i,j)>D25(i,j) && D24(i,j)>D25(i-1,j-1) && D24(i,j)>D25(i,j-1) && D24(i,j)>D25(i+1,j-1) && D24(i,j)>D25(i-1,j) && D24(i,j)>D25(i+1,j) && D24(i,j)>D25(i-1,j+1) && D24(i,j)>D25(i,j+1) && D24(i,j)>D25(i+1,j+1))
                    if(D24(i,j)>threshold)
                    E23(i,j)=1;
                    end
                end
            end
        elseif(D24(i,j)<D24(i-1,j-1) && D24(i,j)<D24(i,j-1) && D24(i,j)<D24(i+1,j-1) && D24(i,j)<D24(i-1,j) && D24(i,j)<D24(i+1,j) && D24(i,j)<D24(i-1,j+1) && D24(i,j)<D24(i,j+1) && D24(i,j)<D24(i+1,j+1))
            if(D24(i,j)<D23(i,j) && D24(i,j)<D23(i-1,j-1) && D24(i,j)<D23(i,j-1) && D24(i,j)<D23(i+1,j-1) && D24(i,j)<D23(i-1,j) && D24(i,j)<D23(i+1,j) && D24(i,j)<D23(i-1,j+1) && D24(i,j)<D23(i,j+1) && D24(i,j)<D23(i+1,j+1))
                if(D24(i,j)<D25(i,j) && D24(i,j)<D25(i-1,j-1) && D24(i,j)<D25(i,j-1) && D24(i,j)<D25(i+1,j-1) && D24(i,j)<D25(i-1,j) && D24(i,j)<D25(i+1,j) && D24(i,j)<D25(i-1,j+1) && D24(i,j)<D25(i,j+1) && D24(i,j)<D25(i+1,j+1))
                    if(D24(i,j)<-threshold)
                    E23(i,j)=-1;
                    end
                end
            end
        end
    end
end

下面显示E21

figure
>> imshow(E21,[-1 1])

不知你们看到这幅图的四五个黑白点没 就是极值点 我自己也是醉了 

后面第三层的就不显示了  我估计就是一片灰色 连这四五个点都不会给我显示出来了。。。

经过上面的分析以及显示的图像  我们知道其实只需要第一层就可以得到候选特征点 所以我只要E11  完整的程序如下所示：

 A=imread('F:\orl_zhifangtu\s1.jpg');
T=40;
[A,R]=kuaisusaomiao(A,T);
[m,n,H]=size(A);
C=[];
x=n/2;
y=m/2;
for u=1:m
    for v=1:n
        i=u;
        j=round(sqrt(R^2-(y-u)^2)*(v-x)/R+x);
        if(R^2-(y-u)^2<0)
            continue;
        end
        C(u,v,1)=A(i,j,1);
        C(u,v,2)=A(i,j,2);
        C(u,v,3)=A(i,j,3);
    end
end
C=uint8(C);
C=rgb2gray(C);
H=C;
C=double(C);
[m1,n1,h1]=size(C);
k=2^(1/3);
threshold=3;
h11=fspecial('gaussian',[5 5],0.3);
C11=imfilter(C,h11,'conv');
h12=fspecial('gaussian',[5 5],0.4);
C12=imfilter(C,h12,'conv');
h13=fspecial('gaussian',[5 5],0.5);
C13=imfilter(C,h13,'conv');
h14=fspecial('gaussian',[5 5],0.6);
C14=imfilter(C,h14,'conv');
h15=fspecial('gaussian',[5 5],0.7);
C15=imfilter(C,h15,'conv');
h16=fspecial('gaussian',[5 5],0.8);
C16=imfilter(C,h16,'conv');
D11=C11-C12;
D12=C12-C13;
D13=C13-C14;
D14=C14-C15;
D15=C15-C16;
E11=zeros(m1,n1);
for i=3:m1-2
    for j=3:n1-2
        if(D12(i,j)>D12(i-1,j-1) && D12(i,j)>D12(i,j-1) && D12(i,j)>D12(i+1,j-1) && D12(i,j)>D12(i-1,j) && D12(i,j)>D12(i+1,j) && D12(i,j)>D12(i-1,j+1) && D12(i,j)>D12(i,j+1) && D12(i,j)>D12(i+1,j+1))
            if(D12(i,j)>D11(i,j) && D12(i,j)>D11(i-1,j-1) && D12(i,j)>D11(i,j-1) && D12(i,j)>D11(i+1,j-1) && D12(i,j)>D11(i-1,j) && D12(i,j)>D11(i+1,j) && D12(i,j)>D11(i-1,j+1) && D12(i,j)>D11(i,j+1) && D12(i,j)>D11(i+1,j+1))
                if(D12(i,j)>D13(i,j) && D12(i,j)>D13(i-1,j-1) && D12(i,j)>D13(i,j-1) && D12(i,j)>D13(i+1,j-1) && D12(i,j)>D13(i-1,j) && D12(i,j)>D13(i+1,j) && D12(i,j)>D13(i-1,j+1) && D12(i,j)>D13(i,j+1) && D12(i,j)>D13(i+1,j+1))
                    if(D12(i,j)>threshold)
                         E11(i,j)=1;
                    end
                end
            end
        elseif(D12(i,j)<D12(i-1,j-1) && D12(i,j)<D12(i,j-1) && D12(i,j)<D12(i+1,j-1) && D12(i,j)<D12(i-1,j) && D12(i,j)<D12(i+1,j) && D12(i,j)<D12(i-1,j+1) && D12(i,j)<D12(i,j+1) && D12(i,j)<D12(i+1,j+1))
            if(D12(i,j)<D11(i,j) && D12(i,j)<D11(i-1,j-1) && D12(i,j)<D11(i,j-1) && D12(i,j)<D11(i+1,j-1) && D12(i,j)<D11(i-1,j) && D12(i,j)<D11(i+1,j) && D12(i,j)<D11(i-1,j+1) && D12(i,j)<D11(i,j+1) && D12(i,j)<D11(i+1,j+1))
                if(D12(i,j)<D13(i,j) && D12(i,j)<D13(i-1,j-1) && D12(i,j)<D13(i,j-1) && D12(i,j)<D13(i+1,j-1) && D12(i,j)<D13(i-1,j) && D12(i,j)<D13(i+1,j) && D12(i,j)<D13(i-1,j+1) && D12(i,j)<D13(i,j+1) && D12(i,j)<D13(i+1,j+1))
                    if(D12(i,j)<-threshold)
                        E11(i,j)=-1;
                    end
                end
            end
        end
    end
end

%下面这部分是把特征点显示在一个原图中  而不是像E11一样显示在一个空矩阵中
[X11,Y11]=find(E11==1 | E11==-1);
imshow(H)
hold on;
plot(Y11,X11,'rx');（这里敲键盘的时候千万别敲错啊 我之前敲错了 换了X11 Y11的位置 结果出来是个亮瞎我铝合金狗眼的图  我自己想了两天 还以为是我哪里编错了）

这就是以上SIFT程序得到的特征点结果  勉强可以吧

后来搜了http://blog.csdn.net/abcjennifer/article/details/7639681她写的MATLAB程序  我真的有一巴掌拍死自己的冲动  人家写得多么简洁  而我的。。。让我自己一刀了断。。。所以你们还是直接下她的程序比较好  短小精悍 又容易懂  

我写这个就当是自己的学习笔记好了。。。用来以后 自己边吐血边回忆。。。