纹理分割（一）Gabor滤波器学习

第一个项目根据公司那边提供的学习资料，需要用到Gabor滤波器对图像进行处理，公司那边关于项目的说法比较商业化，叫X-Ray Image Auto Judging System，之前找了很久论文都没有思路，用这个英文查找论文，也是不对路，这让我在前期浪费不少时间，后来查阅大量论文之后，确定关于目前的项目的学术说法应该是“轮胎X射线图缺陷检测”，英文是“X-Ray
Defect Detection”，目前放暑假，不在学校检索论文比较麻烦，只能根据各路大神网友的博客在做这个项目。

这个项目涉及的主要学术方向是“复杂多纹理图像提取”，能查到的最普遍的做法是采用Gabor滤波器对图像进行处理，可能阅读论文量不到，也许有更好的图像处理方法和思路，这里记录下我目前的做法吧，有效果，但是不是特别好。

看到的大神写到的博客，非常详细，大部分信息整理自以下博客：

1. Gabor滤波器详细介绍：http://mplab.ucsd.edu/tutorials/gabor.pdf

2. Gabor滤波器学习：http://blog.csdn.net/jinshengtao/article/details/17797641

3. 国外大学对Gabor滤波器的参数的讲解：http://matlabserver.cs.rug.nl/edgedetectionweb/web/edgedetection_params.html

4. 备份论文，在校外没法下载：http://www.researchgate.net/publication/252703921_Designing_multiple_Gabor_filters_for_multitexture_image_segmentation

5. 德国某公司的宣传网页，看样子做的非常棒：http://www.cyxplus.fr/produit/our-industy-activities/cyxpert-automatic-defect-detection

一、Gabor Filter 

一维Gabor滤波器

Gabor滤波器是处理一维信号（比如音频）最佳的带通滤波器，一个复杂的Gabor滤波器是一个高斯核函数乘以一个复杂的sin函数（A complex Gabor filter is defined as the product of a Gaussian kernel times a complex
sinusoid），比如：



k  theta  fo是滤波器的参数。我们可以把这个复杂的Gabor滤波器想象成两个相位滤波器分成了一个复杂函数的实部和虚部，



二维Gabor滤波器
但是我要用到的是二维Gabor滤波器，也叫空间Gabor滤波器（The Spatial (2-D) Gabor Filter），但是第一条参考中给出的二维Gabor滤波器讲解太过于。。怎么说。

按照普通的解释来吧，二维Gabor函数的数学表达式：



如何得到的这个公式呢？大神博客（http://blog.csdn.net/yanmy2012/article/details/8090400）中给出的详细求解过程，截图如下：







根据第三条参考文献的解释，关于实部和虚部的说明，只有一句话，实部可以对图像进行平滑滤波，虚部可以用来边缘检测，具体的用法可以可以参考相关论文，不在学校下论文不方便，暂且写这点吧。

下面具体看一下Gabor滤波器的参数说明：







二、Gabor滤波器的应用和适应性

根据第二条参考文献给出的说明，可以总结以下几点：

1. Gabor滤波器可以很好的近似单细胞的感受野细胞（光强刺激下的传递函数），在提取目标的局部空间和频率域信息方面具有良好的特性。

2. 虽然Gabor小波本身不能构成正交基，但在特定参数下可构成紧框架。Gabor小波对于图像的边缘敏感，能够提供良好的方向选择和尺度选择特性，而且对于光照变化不敏感,能够提供对光照变化良好的适应性。-------Gabor小波被广泛应用于视觉信息理解

3. 二维Gabor小波变换是在时频域进行信号分析处理的重要工具，其变换系数有着良好的视觉特性和生物学背景。------因此被广泛应用于图像处理、模式识别等领域。

4. 与传统的傅立叶变换相比，Gabor小波变换具有良好的时频局部化特性。即非常容易地调整Gabor滤波器的方向、基频带宽及中心频率从而能够最好的兼顾信号在时空域和频域中的分辨能力.

5. Gabor小波变换具有多分辨率特性即变焦能力。即采用多通道滤波技术，将一组具有不同时频域特性的Gabor小波应用于图像变换，每个通道都能够得到输入图像的某种局部特性，这样可以根据需要在不同粗细粒度上分析图像。

6. 在特征提取方面，Gabor小波变换与其它方法相比：一方面其处理的数据量较少，能满足系统的实时性要求；另一方面，小波变换对光照变化不敏感，且能容忍一定程度的图像旋转和变形，当采用基于欧氏距离进行识别时，特征模式与待测特征不需要严格的对应，故能提高系统的鲁棒性。

为什么能够有这些应用呢？我们需要深入的理解下Gabor的特性。根据大神博客（http://blog.csdn.net/yanmy2012/article/details/8090400）中的详细解释。在基本视觉皮层里的简单细胞的感受野局限在很小的空域范围内，并且高度结构化。

1. Gabor变换所采用的核（Kernels）与哺乳动物视觉皮层简单细胞2D感受野剖面（Profile）非常相似，具有优良的空间局部性和方向选择性，能够抓住图像局部区域内多个方向的空间频率（尺度）和局部性结构特征。这样，Gabor分解可以看作一个对方向和尺度敏感的有方向性的显微镜。

2. 二维Gabor函数也类似于增强边缘以及峰、谷、脊轮廓等底层图像特征，这相当于增强了被认为是面部关键部件的眼睛、鼻子、嘴巴等信息，同时也增强了诸于黑痣、酒窝、伤疤等局部特征，从而使得在保留总体人脸信息的同时增强局部特性成为可能。它的小波特性说明了Gabor滤波结果是描述图像局部灰度分布的有力工具,因此,可以使用Gabor滤波来抽取图像的纹理信息。

3.  由于Gabor特征具有良好的空间局部性和方向选择性，而且对光照、姿态具有一定的鲁棒性，因此在人脸识别中获得了成功的应用。然而，大部分基于Gabor特征的人脸识别算法中，只应用了Gabor幅值信息，而没有应用相位信息，主要原因是Gabor相位信息随着空间位置呈周期性变化，而幅值的变化相对平滑而稳定，幅值反映了图像的能量谱，Gabor幅值特征通常称为Gabor
能量特征(Gabor Energy Features）。Gabor小波可像放大镜一样放大灰度的变化，人脸的一些关键功能区域(眼睛、鼻子、嘴、眉毛等)的局部特征被强化，从而有利于区分不同的人脸图像。

4. Gabor小波核函数具有与哺育动物大脑皮层简单细胞的二维反射区相同的特性，即具有较强的空间位置和方向选择性，并且能够捕捉对应于空间和频率的局部结构信息；Gabor滤波器对于图像的亮度和对比度变化以及人脸姿态变化具有较强的健壮性，并且它表达的是对人脸识别最为有用的局部特征。Gabor
小波是对高级脊椎动物视觉皮层中的神经元的良好逼近，是时域和频域精确度的一种折中。

三、Gabor滤波器的编程实现

这里我用opencv做的，因为对matlab不是很熟，二是工业应用一般都是C或C++做的，代码给出吧，我用mfc做的界面，之前有好几个版本，验证了几个算法，最后这个就是为了看效果的，比较简单。因为用的opencv3.0，所以显示图像时遇到了麻烦，查找相关资料，解决办法是自己建立CvvImage.h
CvvImage.cpp文件，添加到工程中即可。还有一个要说明的是，Gabor滤波器的实现，是用到的网上搜到的Gabor.h Gabor.cpp实现的，只是进行了简单的注释和文档整理。原作者。。。额，sorry，找的资料太多，不知道出处了，见谅啊！

Gabor.h头文件

#ifndef _GABOR_H
#define _GABOR_H
#include <stdio.h>
#include <iostream>
#include <cv.h>

#define PI 3.14159
#define GAMMA  0.5   //The default value of γ,which is the spatial aspect ratio (sigma_x/sigma_y)
#define RATIO_S2L 0.56  //The default value of σ/λ
#define THETA 45
class Gabor
{
public:
//@construct 构造函数
Gabor(float dLambda, float dTheta, float dRatio_S2L = RATIO_S2L, float dGamma = GAMMA, float dPhi = 0);
//@abolish 析构函数
~Gabor();

//@init 初始化函数
void init(float dLambda, float dTheta, float dPhi, float dGamma = GAMMA);
//@init 初始化函数
void init(float dSigma, float dTheta, float dPhi);
//@init 初始化函数
void init();

//判断Gabor内核是否创建成功--@To find whether the Gabor kernel is created
bool is_kernel(){ return bKernel; }
//判断是否初始化成功--@To find whether the parameters is inited
bool is_init() { return bInit; }
//判断初始化的参数是否足够--@To find whether the parameters inited is enough
bool is_param() { return bParam; }
//得到内核函数所在矩阵--@Get the kernel in matrix form
CvMat* get_Mat() { return pGaborfilter; }
//得到归一化图像--@Get the kernel in image form
IplImage* get_NormImage();
//使用Gabor核函数对输入图像进行处理--@Do the filtering operation to input image with Gabor kernel
IplImage* do_Filter(const IplImage *src);

protected:
bool bParam;        //初始化的参数--bool , if the parameters inited are enough
bool bKernel;        //bool
bool bInit;            //bool
float Lambda;      //余弦函数波长--Wavelength of the cosine factor, which represent the central frequency of Gabor filter
float Theta;          //核函数的方向--Orientation of the Gabor function, the axis x'
float Sigma;         // 标准差--The standard deviation of x, and for y , it is Sigma/Gamma;
float Gamma;       // 空间方向率，指定Gabor函数支持的椭圆率--The spatial aspect ratio
float Phi;              //Gabor的相位偏移--The phase offset of Gabor
CvSize GaborWindow;    //Gabor窗口的宽度--The width of  window
CvMat *pGaborfilter;     //The kernel of Gabor filter

private:
void create_kernel();
};
#endif

Gabor.cpp源文件

#include "stdafx.h"
#include <iostream>
#include <cv.h>
#include <highgui.h>
#include <cstdlib>
#include "Gabor.h"

Gabor::Gabor(float dLambda, float dTheta, float dRatio_S2L, float dGamma, float dPhi)
{
Lambda = dLambda;
Theta = dTheta;
Sigma = dLambda*dRatio_S2L;
Gamma = dGamma;
Phi = dPhi;
pGaborfilter = NULL;
bParam = 1;
}

Gabor::~Gabor()
{
cvReleaseMat(&pGaborfilter);
}
void Gabor::init()
{
float dtmp;
int itmp;
if (is_param() == 0)
{
AfxMessageBox("The parameters are not enough!");
return;
}

//没明白这里是什么意思？
dtmp = sqrt(48 * pow(Sigma, 2) + 1);//根号下( 48*Sigma^2 +1 )
itmp = cvRound(dtmp);//对一个double型的数进行四舍五入，并返回一个整型数！
if (itmp % 2 == 0)
itmp++;
GaborWindow.height = GaborWindow.width = itmp;//创建itmp*itmp的Gabor窗函数
bInit = 1;

create_kernel();
}

void Gabor::init(float dSigma, float dTheta, float dPhi)
{
float dtmp;
int itmp;

Sigma = dSigma;
Theta = dTheta;
Phi = dPhi;
Gamma = GAMMA;
Lambda = Sigma / RATIO_S2L;
bParam = 1;

dtmp = sqrt(24 * pow(Sigma, 2));
itmp = cvRound(dtmp);
if (itmp % 2 == 0)
itmp++;
GaborWindow.height = GaborWindow.width = itmp;
bInit = 1;

create_kernel();
}

void Gabor::init(float dLambda, float dTheta, float dPhi, float dGamma)
{

float dtmp;
int itmp;

Lambda = dLambda;
Theta = dTheta;
Phi = dPhi;
Gamma = dGamma;
Sigma = Lambda * RATIO_S2L;
bParam = 1;

dtmp = sqrt(24 * pow(Sigma, 2));
itmp = cvRound(dtmp);
if (itmp % 2 == 0)
itmp++;
GaborWindow.height = GaborWindow.width = itmp;
bInit = 1;

create_kernel();
}

void Gabor::create_kernel()
{
float tmp1, tmp2, xtmp, ytmp, re;
int i, j, x, y;

if (is_init() == 0)
{
AfxMessageBox("The paremeters haven't been initialed!");
}

pGaborfilter = cvCreateMat(GaborWindow.height, GaborWindow.width, CV_32FC1);
for (i = 0; i < GaborWindow.height; i++)
{
for (j = 0; j < GaborWindow.width; j++)
{
x = j - GaborWindow.width / 2;
y = i - GaborWindow.height / 2;

//源代码此处的计算公式有误
//xtmp = (float)x*cos(Theta) - (float)y*sin(Theta);
//ytmp = (float)x*sin(Theta) + (float)y*cos(Theta);
xtmp = (float)x*cos(Theta) + (float)y*sin(Theta);
ytmp = -(float)x*sin(Theta) + (float)y*cos(Theta);

tmp1 = exp(-(pow(xtmp, 2) + pow(ytmp*Gamma, 2)) / (2 * pow(Sigma, 2)));
tmp2 = cos(2 * PI*xtmp / Lambda + Phi);
// int p=sizeof(float);
re = tmp1*tmp2;
cvSetReal2D((CvMat*)pGaborfilter, i, j, re);
}
}

bKernel = 1;
}

IplImage* Gabor::get_NormImage()
{
if (is_kernel() == 0)
{
AfxMessageBox("The filter hasn't bee created!");
return NULL;
}

IplImage *pImg	= cvCreateImage(GaborWindow, IPL_DEPTH_32F, 1);
IplImage *pImgU8 = cvCreateImage(GaborWindow, IPL_DEPTH_8U, 1);
CvMat * pMat = cvCreateMat(GaborWindow.height, GaborWindow.width, CV_32FC1);

cvCopy(pGaborfilter, pImg);

//归一化，数组的数值被平移或缩放到一个指定的范围
cvNormalize((IplImage*)pImg, (IplImage*)pImg, 0, 255, CV_MINMAX, NULL);
//使用线性变换转换输入数组元素成8位无符号整型
cvConvertScaleAbs(pImg, pImgU8, 1, 0);

return pImgU8;
//return pImg;
}

IplImage * Gabor::do_Filter(const IplImage *src)
{
//判断核函数时候创建成功
if (is_kernel() == false)
{
printf("The Gabor Kernel has not been created!");
return NULL;
}

IplImage *pDestImage = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U, 1);//创建一个同原图像相同的目标图像
IplImage *tmpImg = cvCloneImage(src);//同cvCreateImage一样，只是不需要开辟内存空间，直接将原图像复制到目标图像
IplImage *tmpGrayImg = cvCreateImage(cvSize(src->width, src->height), IPL_DEPTH_8U, 1);

//判断是不是有色图像
if (tmpImg->nChannels != 1)
{
cvCvtColor(tmpImg, tmpGrayImg, CV_BGR2GRAY);//从有色图转到灰度图
}
else
{
cvReleaseImage(&tmpGrayImg);
tmpGrayImg = tmpImg;
}

CvMat *pGaborKernel = get_Mat();//其实就是CvMat *pGaborKernel=pGaborfilter;
//Gabor核函数与原图像进行卷积计算
cvFilter2D(tmpGrayImg, pDestImage, pGaborKernel, cvPoint((GaborWindow.width - 1) / 2, (GaborWindow.height - 1) / 2));

cvReleaseImage(&tmpImg);
cvReleaseImage(&tmpGrayImg);

return  pDestImage;
}

对话框Dlg.cpp，这里只给出了实现按钮的实现函数代码，完整的工程代码，看文章最后的附录：

void CTyreXDlg::OnBnClickedBtndstimg()
{
// TODO:  在此添加控件通知处理程序代码
UpdateData(TRUE);

CSliderCtrl   *pSlidCtrlTheta = (CSliderCtrl*)GetDlgItem(IDC_SLIDER_THETA);
int tmptheta = pSlidCtrlTheta->GetPos();//取得当前位置值

CString sValueTheta = "";
sValueTheta.Format("%d", tmptheta);
SetDlgItemText(IDC_StaticTheta, sValueTheta);

//显示lambda值
CSliderCtrl   *pSlidCtrlLambda = (CSliderCtrl*)GetDlgItem(IDC_SLIDER_LAMBDA);
int tmplambda = pSlidCtrlLambda->GetPos();//取得当前位置值

CString sValueLambda = "";
sValueLambda.Format("%d", tmplambda);
SetDlgItemText(IDC_StaticLambda, sValueLambda);

//图像是否加载成功
if (SrcImg == NULL)
{
AfxMessageBox("没有可处理图像！！！");
return;
}

/*
Gabor滤波
*/
//构造函数
Gabor GaborFilter(/*Lambda*/tmplambda *0.1, /*THETA*/PI*tmptheta / 180, RATIO_S2L, GAMMA, 0);
//初始化
GaborFilter.init();
//获取归一化图像
IplImage* pGaborNormImg = GaborFilter.get_NormImage();
//用Gabor核函数对输入图像处理，返回目标图像
IplImage* poutGaborimg;
poutGaborimg = GaborFilter.do_Filter(SrcImg);

/****************************************/

IplImage *pGrayImage;
// 转为灰度图
pGrayImage = cvCreateImage(cvGetSize(poutGaborimg), IPL_DEPTH_8U, 1);
//判断是不是彩色图像
if (poutGaborimg->nChannels != 1)
cvCvtColor(poutGaborimg, pGrayImage, CV_BGR2GRAY);//从彩色图转到灰度图
else
{
cvReleaseImage(&pGrayImage);
pGrayImage = poutGaborimg;
}

//二值化
IplImage *pBinaryImage;
pBinaryImage = cvCreateImage(cvGetSize(pGrayImage), IPL_DEPTH_8U, 1);
cvThreshold(pGrayImage, pBinaryImage, 0/*二值化阈值*/, 255, CV_THRESH_BINARY);

//显示到picture控件
ShowImgFunc(pBinaryImage, IDC_ShowDstImg);			//二值化图像
//ShowImgFunc(pGaborNormImg, IDC_ShowSrcImg);	//归一化图像
//ShowImgFunc(pGrayImage, IDC_ShowSrcImg);			//二值化之前的灰度图

cvReleaseImage(&DstImg);
cvReleaseImage(&DstGaborImg);
cvReleaseImage(&DstBinaryImg);
}

附录：

1. 上面代码的完整的工程文件，可以到这里下载，不需要积分。我运行成功，但是出现错误，或者有不合理的地方，希望各位看到能告诉我一声，共同进步。

工程文件完整下载地址：http://download.csdn.net/detail/jorg_zhao/8949247