鱼眼拼接之Harris检测特征点

理论部分是根据《鱼眼图像全景拼接系统》

%我用的是鱼眼图像  所以先提取区域校正

A=imread('F:\orl_zhifangtu\s1.jpg');

T=40;

[A,R]=kuaisusaomiao(A,T);

[m,n,k]=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);

k=0.045;

threshold=1000000;

cnt=0;（特征点个数  可以在MATLAB里面看到）

Rmax=0;

result=zeros(m,n);

a=[-1 0 1];

b=a';

Ix=filter2(a,C);

Iy=filter2(b,C);

w=fspecial('gaussian',[5 5],0.8);

E=filter2(w,Ix.^2);

F=filter2(w,Iy.^2);

G=filter2(w,Ix.*Iy);

for x=3:m-2

      for y=3:n-2

         for i=-2:2

             forj=-2:2

                RR(x,y)=E(x,y)*F(x,y)-(G(x,y))^2-k*(E(x,y)+F(x,y))^2;

                RR(x+i,y+j)=E(x+i,y+j)*F(x+i,y+j)-(G(x+i,y+j))^2-k*(E(x+i,y+j)+F(x+i,y+j))^2;

                if(RR(x,y)<RR(x+i,y+j))

                   break;

                else

                   RRmax=RR(x,y);

                end

               end

           end

          if(RRmax>threshold)

              cnt=cnt+1;

              result(x,y)=1;

            end

       end

end

>> [X,Y]=find(result==1);

imshow(C);

hold on

plot(Y,X,'rx');

显示出来效果如第二幅图所示：

当阈值取threshold=500000;时效果如下所示：

效果不好 特征点太多   而且有的特征点根本没有显示出来

当阈值改为threshold=1000000;时 特征点比上幅图减少了  但是也没有SIFT效果好  可能是编的不够好

这个程序要自己修改阈值 而且不知道该修改成多少才合适  所以不好 

后来在百度文库上搜到一个程序  效果挺好的 而且简单易懂  网址就不贴了  忘记了当时

前面的鱼眼处理部分不变

只是改了后面的检测特征点那里

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);

%随便拿一幅图像 都可以从以下开始进行harris特征点检测

C=rgb2gray(C);

>> ori_im=C;

>> fx=[-2 -1 0 1 2];

>> Ix=filter2(fx,ori_im);

>> fy=[-2;-1;0;1;2];

>> Iy=filter2(fy,ori_im);

>> Ix2=Ix.^2;

>> Iy2=Iy.^2;

>> Ixy=Ix.*Iy;

>> h=fspecial('gaussian',[77],0.8);

Ix2=filter2(h,Ix2);

>> Iy2=filter2(h,Iy2);

>> Ixy=filter2(h,Ixy);

>> result=zeros(m,n);

>> R=zeros(m,n);

>> Rmax=0;

>> for i=1:m

     for j=1:n

        M=[Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];

        R(i,j)=det(M)-0.04*(trace(M))^2;

        if(R(i,j)>Rmax)

            Rmax=R(i,j);

        end

    end

end

>> cnt=0;

>> for i=2:m-1

      for j=2:n-1

        if(R(i,j) > 0.0008*Rmax &&R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) >R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1)&& R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) &&R(i,j) > R(i+1,j+1))

           result(i,j)=1;

           cnt=cnt+1;

         end

     end

end

>>[posc,posr]=find(result==1);

>> imshow(ori_im);

>> hold on;

>> plot(posr,posc,'rx');

他这个程序是先遍历所有找到最大的那个值  然后把最大值的0.0008倍设置为阈值  当然这个倍数也可以自己调节

效果还不错 如下：

效果还是挺棒棒哒的。。。最起码比我之前的那个效果好多了