风机桨叶故障诊断（二） 获取图像几何主方向

    昨天，我将视频资源按帧抽取并筛选得到了可以用来提取样本的图像库。今天还是进行项目的准备工作。当我们拿到一张图片，我们的软件要做的大致可以分为三个步骤：从原图中识别桨叶——从桨叶中寻找故障——诊断故障类型。

    我们第一步要识别的核心物体是桨叶，首先桨叶的运动是有其特点的，它绕着一个轴在旋转，所以我们得到的图像中桨叶的方向是360°任意的。就如下图所示：

            

    我考虑到由于样本数量并不庞大且多样，设计一个算法识别出一个物体是否为桨叶（且提供的桨叶图像方向任意）的方案并不理想。我感觉如果我们对于任意的桨叶图像，能够识别出桨叶所指的方向，然后将其统一旋转到水平方向来归一化样本，应该是个不错的选择！

    现在的问题就是，如何将图像中物体的几何主方向提取出来。

    首先我们想提取桨叶的几何轮廓信息，而要排除背景等因素的干扰。考虑使用边缘提取算法。下面我们均以下图所示桨叶为例，展示算法的结果

        

    下图是采用不同算子进行边缘提取后的结果

    

    注：上图所示结果在边缘提取前进行了高斯滤波，我在试验了批量样本后发现在边缘提取前进行滤波是必要的，可以排除背景中包括云层等因素的影响。下面是某个桨叶图像直接提取的结果（未滤波）

    

    进行了大量的探索之后我们发现：

    1.在边缘提取前有必要进行滤波，我使用的是高斯滤波

    2.边缘提取时sobel算子和roberts算子表现比较好，所以后面的边缘提取我将采用sobel算子，我再赘述。

    到这里，我们得到了桨叶的几何轮廓信息，就像下图所示，怎么从中得到桨叶的方向呢？

     

     我们可以看到桨叶的几何轮廓非常具有特点，大致由两条角度近似的直线相交在一起。可以看出，这样的轮廓，其具有明显的像素值的梯度方向，也就是与两条角度近似直线方向的垂直方向。换个角度想，某个角度的梯度大小与频率域的某个方向的总能量大小时对应的。所以我们可以用二维傅里叶变换将图像变换到频率域，用结果求出对应的能量谱，将频率域沿原点上半平面分180份（因为下半平面与上半平面的角度是一一对应的，所以只考虑上半平面就可以了），分别计算每一份小扇形所包含的区域内的能量总和，也就代表着这个角度上所包含的能量大小。

     上图边缘提取结果再进行二维傅里叶变换后的能量谱如下图所示（图像中心对应于0频率），可以明显看出能量主要集中在水平方向。

     

     我们取180个角度中能量最大的那个作为我们的主方向输出，即可完成算法。算法输出的主方向角度和对应的桨叶图像示例如下：

               41°

                     1°

                         135°

    由上面的结果可以看出，对于识别桨叶所指角度的问题，这个方法十分有效且可靠。这样，我们拿到一个桨叶图像后，即可调用算法识别主方向，然后将其旋转到水平方向，这样我们的样本就得到了归一化，相信之后的预测算法能够表现的更好。

    不过我们还可以注意到一个问题，算法对于负角度的桨叶输出结果是与之对应的正角度，在后面的工作中，应该要再用一个方法区分一下正反。

     识别图像中物体主方向算法的matlab实现如下：

function [mdex] = MainDirec(a)
%根据图像边缘信息估计纹理主方向%若是彩色图像，将其变换为灰度图像
if(size(a,3)~=1)a=rgb2gray(a);
end%使用sobel算子进行边缘提取
[w,h]=size(a);
edge_a=edge(a,'sobel');%注：fft2是2维离散傅立叶变换
% 傅里叶变换后，如果使用fftshift命令，0频率分量将会移到坐标中心
fedge_a=fftshift(fft2(edge_a));%求傅里叶变换的实部、虚部
realf=real(fedge_a);
imagf=imag(fedge_a);
%求傅里叶变换的能量
pfa=sqrt(realf.^2+imagf.^2);
%绘制出能量谱看下
imshow(pfa);%dhist用来记录每个角度上的能量
dhist=zeros(1,180);%0频率分量出的坐标（坐标中心）
x0=w/2; y0=h/2;for i=1:wfor j=1:h% 函数 theat=atan2(imagf,realf); 将 theat 归一化到 [0,pi]temp=atan2(j-y0,i-x0);if(temp<0)temp=temp+pi;end%round 向最近的整数取整%mod(a,b)就是求的是a除以b的余数%下面的一行代码将[0,pi]的角度换回[0,179],这样temp+1的范围就是[1,180]temp=mod(round(temp*180/pi),180);%dhist用来记录每个角度上的能量%将当前位置的能量累加到记录当前位置所属角度的能量的变量dhist中dhist(temp+1)=dhist(temp+1)+pfa(i,j);end
end% %将精确的能量谱简化成分段的形式
% dsp=[0:10:170;10:10:180];
% hist=zeros(1,18);
% for i=1:180
%     for k=1:18
%         if(i>dsp(1,k) && i<=dsp(2,k))
%             hist(k)=hist(k)+dhist(i);
%         end
%     end
% end
% [mhist,mdex]=max(hist);[mhist,mdex]=max(dhist);end