Shi-Tomasi算子的运用 ，用于检测角点

角点检测

当一个窗口在图像上移动，在平滑区域如图(a)，窗口在各个方向上没有变化。在边缘上如图(b)，窗口在边缘的方向上没有变化。在角点处如图(c)，窗口在各个方向上具有变化。Harris角点检测正是利用了这个直观的物理现象，通过窗口在各个方向上的变化程度，决定是否为角点。

将图像窗口平移[u,v]产生灰度变化E(u,v)

由：， 得到：

对于局部微小的移动量 [u,v]，近似表达为：

其中M是 2*2 矩阵，可由图像的导数求得：

E(u,v)的椭圆形式如下图：

 

定义角点响应函数 R 为：

Harris角点检测算法就是对角点响应函数R进行阈值处理：R > threshold，即提取R的局部极大值。

 

Shi--Tomasi角点检测法，如果像素点的最小特征值大于周围像素的特征值，则该点是角点。

代码：

<span style="font-size:18px;"><strong>im=imread('lena.jpg');
tau=100;
im=double(im);
keyXs=[];
keyYs=[];
win=3;
[height,width] = size(im);
result = zeros(height,width);
%Then I will get the gradients of the image along the x and y axises.
sobel_x=1/4*[-1 0 1;-2 0 2;-1 0 1];
sobel_y=1/4*[-1 0 1;-2 0 2;-1 0 1]';
diffx=imfilter(im,sobel_x);         %对图像x方向进行梯度
diffy=imfilter(im,sobel_y);       %对图像y方向的梯度进行计算
%For smoothing the differentiation of the image along the x and y
%direction, the gauss filter of the diffx and diffy is must.
gauss_win=win;
sigma=1;
[x,y]=meshgrid(-gauss_win:gauss_win,-gauss_win:gauss_win);
gauss2D=exp(-(x.^2+y.^2)/(2*sigma.^2));  %产生高斯算子
gauss2D=gauss2D/(sum(sum(gauss2D)));  %对高斯算子进行归一化
%Then calculate the M matrix.
A=imfilter(diffx.*diffx,gauss2D);      %二阶x方向梯度进行高斯滤波
B=imfilter(diffy.*diffy,gauss2D);      %二阶y方向梯度进行高斯滤波
C=imfilter(diffx.*diffy,gauss2D);      %对图像x y方向的梯度进行高斯滤波
supress_win=2;
    threshold=100;
    points_count=0;
    bigger=zeros(height,width);
    smaller=zeros(height,width);
    for x=1:width
        for y=1:height
            M=[A(y,x) C(y,x);C(y,x) B(y,x)];
            %It is too time-consuming.
            %eigenvalue=eig(M);
            %bigger(y,x)=max(eigenvalue);
            %smaller(y,x)=min(eigenvalue);
            temp1=M(1,1)+M(2,2);
            temp2=sqrt((M(1,1)-M(2,2))^2+4*M(1,2)^2);
            bigger(y,x)=(temp1+temp2)/2;
            smaller(y,x)=(temp1-temp2)/2;
        end
    end
    for x=supress_win+1:width-supress_win
        for y=supress_win+1:height-supress_win     
            temp=smaller(y,x);
            if(temp>threshold)
                %Then I will make the non-maximumu suppression to the
                %samller matrix after the threholding.
                flag=0;
                for i=-supress_win:supress_win
                    for j=-supress_win:supress_win
                        if(temp>=smaller(y+j,x+i))
                            flag=flag+1;
                        end
                    end
                end
                if(flag==((2*supress_win+1)*(2*supress_win+1)))
                    result(y,x)=1;
                    points_count=points_count+1;
                    keyXs(points_count)=x;
                    keyYs(points_count)=y;
                end
            end
        end
    end
end
</strong></span>