Gabor及其编程

Gabor变换属于加窗傅立叶变换，Gabor函数可以在频域不同尺度、不同方向上提取相关的特征。另外Gabor函数与人眼的生物作用相仿，所以经常用作纹理识别上，并取得了较好的效果。Gabor变换是短时Fourier变换中当窗函数取为高斯函数时的一种特殊情况.

Gabor变换的本质实际上还是对二维图像求卷积。因此二维卷积运算的效率就直接决定了Gabor变换的效率。在这里我先说说二维卷积运算以及如何通过二维傅立叶变换提高卷积运算效率.关于离散二维叠加和卷积的运算介绍的书籍比较多，我这里推荐William K. Pratt著，邓鲁华 张延恒 等译的《数字图像处理（第3版）》，其中第7章介绍的就是这方面的运算.

A可以理解成是待处理的笔迹纹理，B可以理解成Gabor变换的核函数，现在要求A与B的离散二维叠加卷积，我们首先对A的右（？）边界和下（？）边界填充0（zero padding)，然后将B进行水平翻转和垂直翻转，如下图：

然后用B中的每个值依次乘以A中相对位置处的值并进行累加，结果填入相应位置处（注意红圈位置）。通常二维卷积的结果比A、B的尺寸要大。如下图所示：

2、快速傅立叶变换卷积

根据傅立叶变换理论，对图像进行二维卷积等价于对图像的二维傅立叶变换以及核函数的二维傅立叶变换在频域求乘法。通过二维傅立叶变换可以有效提高卷积的运算效率。但在进行傅立叶变换时一定要注意“卷绕误差效应”，只有正确对原有图像以及卷积核填补零后，才能得到正确的卷积结果。

二维Gabor函数可以表示为：

其中:

v的取值决定了Gabor滤波的波长，u的取值表示Gabor核函数的方向，K表示总的方向数。参数决定了高斯窗口的大小，这里取. 程序中取4个频率（v=0, 1, ..., 3），8个方向（即K=8，u＝0, 1, ... ,7），共32个Gabor核函数。不同频率不同方向的Gabor函数可通过下图表示：

 图片来源：http://www.neuroinformatik.ruhr-uni-bochum.de/ini/VDM/research/computerVision/imageProcessing/wavelets/gabor/gaborFilter.html

图片来源：http://www.bmva.ac.uk/bmvc/1997/papers/033/node2.html

所谓目标识别，从某种意义上说就是特征识别的问题。而红外图特征提取的角度一般可以从几何形状和上下文的判断中得到，比如说，当你在一幅图像中去搜索桥梁和机场的跑道的时候，我们可以从形状上发现他们都是一组平行线；然后基于上下文去判断，桥梁的长度肯定比跑道短，桥梁的两边一般是水域，而跑道的中间有一些联络道，而这些这些联络道本身也是平行线。
  
通常用这种仅仅利用图像本身的知识去判断的时候，时常由于经验的限制，识别率总是没有办法做到很好！虽然图像所一些早期的论文里都声称，识别率能够达到90％以上，但是那只是用自己搜集的有利于自己算法的图片得出来的结果而已！
 
言规正传，由于种种努力，我们发现了频率的概念，频率的本质就是反映变化快慢的物理量。因此当我们把一幅图像f(x,y)通过傅立叶变换变成F(u,v)进入频率域的时候，w＝sqrt(u2+v2) 的大小，则直接体现了原图像中灰度值变化速率的大小。于是人们发现了滤波的方法，如通过低通滤波器将图像中的噪声去掉，同时把一些突变的地方变得平滑，而高通滤波器则可以把图像的边缘特征加以锐化。
 
但是很快大家就发现了傅立叶变换的一个缺陷，傅立叶变换在很大程度上只是一个体现大局观的概念，而无法反映图像局部的特点，因此我们无法通过傅立叶变换来对图像进行细节上更准确的操作。幸运的是，人们之后的研究中，发现了Gabor于1946年发表的论文《the Theory of Communication》,于是在识别的方法上，尤其在纹理识别上达到了一个新的高度。当然小波变换的应用也是差不多这个时间段开始的，由此我们也可以理解数学在图像领域的价值！
 
Gabor滤波的本质是提取图像的特征分量，通常应用的领域有指纹识别，虹膜识别，人脸识别等。
 
所谓gabor的价值由来是由于人们对傅立叶频率分析的需求来的，因此，当我们去使用gabor滤波器的时候，就一定要结合对图象频率的分析。这一点表面上简单而直接，实际却被大部分的研究者忽略。关于这部分的研究，可以在各大搜索引擎上检索关键："fourier power spectrum"。如果使用matlab，可以使用下面的语句将图象转换到频率域：

function fps(I)
g=log(abs(fftshift(fft2(double(I)))));  %通过一定的技巧我们也可以将频率精度扩大
imview(g,[]); % 注意不能用imshow;

      传统的gabor滤波器的思想，就是把图象不同的频率范围分别滤出来，然后提取滤波图象的特征分量加以分类。而一些对gabor滤波器的设计正是从这个角度来考虑。至于 相邻的滤波器在频率区域的有效部分是相切还是重叠，则一直是争论的话题。甚至有些人因此从gabor小波的角度去设计。
       关于gabor滤波器的参数本身有4个（wf:center frequency  theta：angle  sigmax,sigmay：extend of gaussian function），90年代的时候，很多文献将sigmax=sigmay，于是参数变成3个；本世纪的时候，人们由于我在上面一段阐述的道理，将sigma与wf组成反比关系。参数又变成2个。这就是现在通用的关于尺度和角度的设计。

      在这里我仍然不想去针对特定的纹理，比如指纹，虹膜，文字之类的具体实现来讲更多的东西；因为那实在是已经成熟的技术。只要你慢慢照着做，就一定能达到一定程度的成功率。

我想说的是一些更有趣的东西，这就是频率域本身的矛盾性；即大部分信息包含在低频部分，但低频部分的范围却很小，难以区分，因此我们其实浪费了最大量的信息。关于角度的信息，我们只要侦察到特定的角度就一定可以从相当广的频率范围下找到；但是尺度信息的话，我们只能徘徊在一个很小的区间。并且有趣的是这里很多现象象极了物理上的衍射现象。而我们一旦想深入到低频部分的话，对滤波器的设计简直就是苛刻。滤波器的鲁棒性在这里没有价值，我们需要象狐狸一样敏感。但是我相信，如果能够解决这一问题，那么我们甚至可以摈弃纹理识别的概念，而是可以说对图象识别的成功率将到达一个新的高度。

３.编程

目前网上可以找到开源C#版的快速傅立叶变换代码（Exocortex.DSP），我使用的是1.2版，2.0版似乎只能通过CVS从SourceForge上签出， 并且功能没有什么太大改变。将Exocortex.DSP下载下来后，将源代码包含在自己的项目中，然后就可以利用它里面提供的复数运算以及傅立叶变换功能了。为了测试通过傅立叶变换求卷积的有效性，特编写以下代码：

using System;
using Exocortex.DSP;
{
	static void Main()
	{
		fftConv2 c = new fftConv2();
		c.DoFFTConv2();
	}
}
public class fftConv2
{
	double[,] kernel = {	
				{-1, 1}, 
				{0, 1}
			};
	double[,] data = {
					{10,5,20,20,20},
					{10,5,20,20,20},
					{10,5,20,20,20},
					{10,5,20,20,20},
					{10,5,20,20,20}
			};
	Complex[] Kernel = new Complex[8*8];
	Complex[] Data = new Complex[8*8];
	Complex[] Result = new Complex[8*8];
	
	private void Init()
	{
		for(int y=0; y<2; y++)
			for(int x=0; x<2; x++)
				Kernel[y*8+x].Re = kernel[y,x];
		for(int y=0; y<5; y++)
			for(int x=0; x<5; x++)
				Data[y*8+x].Re = data[y,x];
	}
	public void DoFFTConv2()
	{
		Init();
		Fourier.FFT2(Data, 8, 8, FourierDirection.Forward);
		Fourier.FFT2(Kernel, 8, 8, FourierDirection.Forward);
		for(int i=0; i<8*8; i++)
			Result[i] = Data[i] * Kernel[i] / (8*8);
		Fourier.FFT2(Result, 8, 8, FourierDirection.Backward);
		for(int y=0; y<6; y++)
		{
			for(int x=0; x<6; x++)
				Console.Write("{0,8:F2}", Result[y*8+x].Re);
			Console.WriteLine();
		}
	}
}
class MainEntry

程序的运行结果与离散二维叠加和卷积的运算结果完全相同。

由于卷积结果与原始输入图片的大小是不一样的，存在着所谓“边界”，在我的实际应用程序中，为了避免这些“边界”对结果过多的影响，我采用的是居中阵列定义，并且从卷积结果中只截取需要的那部分内容，确保和原始图片的大小完全一致，如下图：


这就需要对卷积的傅立叶求法做些微小的调整，具体调整办法就不说了，主要是坐标的变换，将示例代码贴上来供大家参考：

using System;
using Exocortex.DSP;
class MainEntry
{
	static void Main()
	{
		CenterfftConv2 s = new CenterfftConv2();
		s.CommonMethod();
		s.DoFFTConv2();
	}
}
public class CenterfftConv2
{
	double[,] kernel = {
				{0, 1, 0}, 
				{1, 2, 0},
				{0, 0, 3}
				};
	double[,] data = new double[12,12];
	Complex[] Kernel = new Complex[16*16];
	Complex[] Data = new Complex[16*16];
	Complex[] Result = new Complex[16*16];
	public CenterfftConv2()
	{
		Random r = new Random();
		for(int y=0; y<12; y++)
			for(int x=0; x<12; x++)
				data[y,x] = r.NextDouble();
		for(int y=0; y<3; y++)
			for(int x=0; x<3; x++)
				Kernel[y*16+x].Re = kernel[y,x];
		for(int y=1; y<13; y++)
			for(int x=1; x<13; x++)
				Data[y*16+x].Re = data[y-1,x-1];
	}
public void DoFFTConv2()
{
	Console.WriteLine(" ========= By FFT2Conv2 Method =========");
	Fourier.FFT2(Data, 16, 16, FourierDirection.Forward);
	Fourier.FFT2(Kernel, 16, 16, FourierDirection.Forward);
	for(int i=0; i<16*16; i++)
		Result[i] = Data[i] * Kernel[i] / (16*16);
	Fourier.FFT2(Result, 16, 16, FourierDirection.Backward);
	for(int y=2; y<14; y++)
	{
		for(int x=2; x<14; x++)
			Console.Write("{0,5:F2}", Result[y*16+x].GetModulus());
		Console.WriteLine();
	}
}
public void CommonMethod()
{
	double real = 0;
	Console.WriteLine(" ========== Direct Transform ===========");
	for(int y=0; y < 12; y++)
	{
		for(int x=0; x < 12; x++)
		{
			for(int y1=0; y1 < 3; y1++)
				for(int x1=0; x1 < 3; x1++)
				{
					// (kernel.Length-1)/2 = 1
					if(((y - 1 + y1)>=0) &&((y - 1 + y1)<12) &&((x - 1 + x1)>=0) &&((x - 1 + x1)<12))
					{
						real += data[y - 1 + y1, x - 1 + x1] * kernel[2 - x1, 2 - y1]; 
					}
				}
			Console.Write("{0,5:F2}", real);
			real=0;
		}
		Console.WriteLine();
	}
	Console.WriteLine("n");
}

有了此部分的基础知识后，我们就要步入笔迹识别中最核心的部分Gabor变换，提取笔迹的特征了。

４、代码实现
Gabor函数是复值函数，因此在运算过程中要分别计算其实部和虚部。代码如下：

private void CalculateKernel(int Orientation, int Frequency)
{
	double real, img;
	for(int x = -(GaborWidth-1)/2; x<(GaborWidth-1)/2+1; x++)
		for(int y = -(GaborHeight-1)/2; y<(GaborHeight-1)/2+1; y++)
		{
			real = KernelRealPart(x, y, Orientation, Frequency);
			img = KernelImgPart(x, y, Orientation, Frequency);
			KernelFFT2[(x+(GaborWidth-1)/2) + 256 * (y+(GaborHeight-1)/2)].Re = real;
			KernelFFT2[(x+(GaborWidth-1)/2) + 256 * (y+(GaborHeight-1)/2)].Im = img;
		}
}
private double KernelRealPart(int x, int y, int Orientation, int Frequency)
{
	double U, V;
	double Sigma, Kv, Qu;
	double tmp1, tmp2;
	U = Orientation;
	V = Frequency;
	Sigma = 2 * Math.PI * Math.PI;
	Kv = Math.PI * Math.Exp((-(V+2)/2)*Math.Log(2, Math.E));
	Qu = U * Math.PI / 8;
	tmp1 = Math.Exp(-(Kv * Kv * ( x*x + y*y)/(2 * Sigma)));
	tmp2 = Math.Cos(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(Sigma/2));
	return tmp1 * tmp2 * Kv * Kv / Sigma; 
}
private double KernelImgPart(int x, int y, int Orientation, int Frequency)
{
	double U, V;
	double Sigma, Kv, Qu;
	double tmp1, tmp2;
	U = Orientation;
	V = Frequency;
	Sigma = 2 * Math.PI * Math.PI;
	Kv = Math.PI * Math.Exp((-(V+2)/2)*Math.Log(2, Math.E));
	Qu = U * Math.PI / 8;
	tmp1 = Math.Exp(-(Kv * Kv * ( x*x + y*y)/(2 * Sigma)));
	tmp2 = Math.Sin(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(Sigma/2));
	return tmp1 * tmp2 * Kv * Kv / Sigma; 
}

有了Gabor核函数后就可以采用前文中提到的“离散二维叠加和卷积”或“快速傅立叶变换卷积”的方法求解Gabor变换，并对变换结果求均值和方差作为提取的特征。32个Gabor核函数对应32次变换可以提取64个特征（包括均值和方差）。由于整个变换过程代码比较复杂，这里仅提供测试代码供下载。该代码仅计算了一个101×101尺寸的Gabor函数变换，得到均值和方差。代码采用两种卷积计算方式，从结果中可以看出 ，快速傅立叶变换卷积的效率是离散二维叠加和卷积的近50倍。