笔迹鉴别(3) —— Gabor变换(1)

一、二维卷积运算

Gabor变换的本质实际上还是对二维图像求卷积。因此二维卷积运算的效率就直接决定了Gabor变换的效率。在这里我先说说二维卷积运算以及如何通过二维傅立叶变换提高卷积运算效率。在下一步分内容中我们将此应用到Gabor变换上，抽取笔迹纹理的特征。

1、离散二维叠加和卷积

关于离散二维叠加和卷积的运算介绍的书籍比较多，我这里推荐William K. Pratt著，邓鲁华 张延恒 等译的《数字图像处理（第3版）》，其中第7章介绍的就是这方面的运算。为了便于理解，我用下面几个图来说明离散二维叠加和卷积的求解过程。

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

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

2、快速傅立叶变换卷积

根据傅立叶变换理论，对图像进行二维卷积等价于对图像的二维傅立叶变换以及核函数的二维傅立叶变换在频域求乘法。通过二维傅立叶变换可以有效提高卷积的运算效率。但在进行傅立叶变换时一定要注意“卷绕误差效应”，只有正确对原有图像以及卷积核填补零后，才能得到正确的卷积结果。关于这部分内容可以参考William K. Pratt著，邓鲁华 张延恒 等译的《数字图像处理（第3版）》第9章的相关内容，此处就不再赘述。

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

using System;
using Exocortex.DSP;
class MainEntry
{
   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();
      }
   }
}

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

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

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

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变换，提取笔迹的特征了。