边缘提取（二）：Steerable filter简单理解

Steerable filter简单理解

（如果要看代码，可以直接跳转到下面，中间的理论推导可能有废话，如果要看steerable的原文，The Design and Use of Steerable Filters  William T，91年的老文章了）

     若多元函数  f在点  存在对所有自变量的偏导数，则称向量gradf 为函数 f 在点  的梯度，记向量gradf 的长度（或模）为[2]

                                                                       [1]

                                                       [2]

（上述这段摘自百度百科）

对于函数来说，变化最剧烈的方向即函数在该点的梯度方向，对于一幅灰度图im(x,y)来说，我们可以将其看成是一个二元函数（p为像素值，im为原图）

那么图中每个点像素值变化最剧烈的方向即是该点的梯度方向（如图1中的橘色箭头，指的就是箭头端点处的梯度方向）。理解上述观点只需要基本的高等数学的知识，如果想要增强理论基础，可以参见                https://blog.csdn.net/myarrow/article/details/51332421）。

                                                       （图 1 橘色箭头指向的是箭头端点处的边缘方向）

边缘是像素值变化相对剧烈的点，且该点的梯度方向是像素值变化最剧烈的方向，我们将与该点梯度方向垂直的方向称为该点的边缘方向。（如图一中的蓝色箭头所示）

      求取图像边缘的过程其实就可以等价于求导，然后筛选导数绝对值较大的位置。图像作为一个二元函数，求导准确来说是求方向导数，沿着不同方向求导，得到的值也不一样。根据多元函数的知识我们可以知道，方向导数其实就是该点梯度与该方向方向余弦的点积（投影）。

                                                                                                                 [3]

所以当求导方向与该点梯度平行时，最能突出该点的变化，如果垂直，则完全忽视了该点的变化。这样的性质在边缘提取中可以用来特意的突出某一方向的边缘，是设计steerable filter 的基础，理解了上面的话之后，设计一个提取特定方向的steerable filter就很简单了。

设图像 ，则该点的梯度gradf为, ，设定要凸出边缘朝向为的边缘，则方向导数的方向应该为，方向余弦为

                                                                         [4]

设gx，gy都是沿x方向和y方向的卷积核，公式4给出的是利用卷积求图像偏导的公式，gx，gy可以是一阶的高斯导数离散化之后的模板，也可以是诸如sobel沿x和y方向的两个卷积核。根据公式3和上面给定的求导方向，我们可以设计出一个新的卷积核（5），专门用于凸出某一个方向的边缘

                                                                                                              [5] 

 

           [6]

在具体实现过程中，如果不需要设计新的卷积模板，那么可以直接对横向边缘图和垂直边缘图按照方向余弦加权（公式6），这是等价的（这是有卷积的运算性质决定的）

 matlabcode:

clearvars;
close all;
im=imread('1.png');
im=rgb2gray(im);
theta=0;%要提取的边缘方向
theta1=theta;
sigma=1;
Wsize=1;
theta = theta/180*pi;%转换为弧度
% 计算二维高斯核在x,y方向的偏导gx,gy
k =-Wsize:Wsize;
%%
%卷积核的计算,利用一阶差分的高斯函数做的，很简单，只需要计算一下一阶导数就可以了
mm=fspecial('gaussian',[3 3],1);
nn=zeros(3,3);
nn(:,1)=1/(sigma^2);
nn(:,3)=-1/(sigma^2);
gx=nn.*mm;
gy=nn'.*mm;%得到离散化的高斯一阶差分算子
%%
% 计算图像Ｉ在x,y方向的滤波结果
Ix = conv2(im,gx,'same');
Iy = conv2(im,gy,'same');
g_new=cos(theta)*gx+sin(theta)*gy;
% 计算图像Ｉ在theta方向的滤波结果
J = cos(theta)*Ix+sin(theta)*Iy;
J2=conv2(im,g_new);
% J=J;%给定一个固定的阈值
% J2=J2;
%J2=cos(pi-theta)*Ix+sin(pi-theta)*Iy;%J2=J2>10;
%J=J+J2;
figure,
subplot(1,3,1),axis image; colormap(gray);imshow(im),title('原图像');
subplot(1,3,2),axis image; colormap(gray);imshow(cos(theta)*gx+sin(theta)*gy),title('滤波模板');
subplot(1,3,3),axis image; colormap(gray);imshow(J),title('滤波结果'); 
figure,subplot(2,1,1);imshow(J,[]);title(num2str(theta1));subplot(2,1,2);imshow(J2,[]);title(num2str(theta1));

效果图：

                                                         

 

 

                                       （可以很清楚的看到不同方向的边缘提取算子突出的边缘是不一样的）