Gabor滤波器在频域中的实现

这里我们有Gaborfilter的Spatial域实现。 但是， 出于性能原因 ，我需要在频域中实现Gabor滤波器。

我找到了Gabor滤波器的频域方程 ：

^{我实际上对这个公式的正确性和/或适用性存有疑问。}

源代码

所以，我已经实现了以下内容：

public partial class GaborFfftForm : Form { private double Gabor(double u, double v, double f0, double theta, double a, double b) { double rad = Math.PI / 180 * theta; double uDash = u * Math.Cos(rad) + v * Math.Sin(rad); double vDash = (-1) * u * Math.Sin(rad) + v * Math.Cos(rad); return Math.Exp((-1) * Math.PI * Math.PI * ((uDash - f0) / (a * a)) + (vDash / (b * b))); } public Array2d GaborKernelFft(int sizeX, int sizeY, double f0, double theta, double a, double b) { int halfX = sizeX / 2; int halfY = sizeY / 2; Array2d kernel = new Array2d(sizeX, sizeY); for (int u = -halfX; u < halfX; u++) { for (int v = -halfY; v < halfY; v++) { double g = Gabor(u, v, f0, theta, a, b); kernel[u + halfX, v + halfY] = new Complex(g, 0); } } return kernel; } public GaborFfftForm() { InitializeComponent(); Bitmap image = DataConverter2d.ReadGray(StandardImage.LenaGray); Array2d dImage = DataConverter2d.ToDouble(image); int newWidth = Tools.ToNextPowerOfTwo(dImage.Width) * 2; int newHeight = Tools.ToNextPowerOfTwo(dImage.Height) * 2; double u0 = 0.2; double v0 = 0.2; double alpha = 10;//1.5; double beta = alpha; Array2d kernel2d = GaborKernelFft(newWidth, newHeight, u0, v0, alpha, beta); dImage.PadTo(newWidth, newHeight); Array2d cImage = DataConverter2d.ToComplex(dImage); Array2d fImage = FourierTransform.ForwardFft(cImage); // FFT convolution ................................................. Array2d fOutput = new Array2d(newWidth, newHeight); for (int x = 0; x < newWidth; x++) { for (int y = 0; y < newHeight; y++) { fOutput[x, y] = fImage[x, y] * kernel2d[x, y]; } } Array2d cOutput = FourierTransform.InverseFft(fOutput); Array2d dOutput = Rescale2d.Rescale(DataConverter2d.ToDouble(cOutput)); //dOutput.CropBy((newWidth-image.Width)/2, (newHeight - image.Height)/2); Bitmap output = DataConverter2d.ToBitmap(dOutput, image.PixelFormat); Array2d cKernel = FourierTransform.InverseFft(kernel2d); cKernel = FourierTransform.RemoveFFTShift(cKernel); Array2d dKernel = DataConverter2d.ToDouble(cKernel); Array2d dRescaledKernel = Rescale2d.Rescale(dKernel); Bitmap kernel = DataConverter2d.ToBitmap(dRescaledKernel, image.PixelFormat); pictureBox1.Image = image; pictureBox2.Image = kernel; pictureBox3.Image = output; } } 

^{现在只关注算法步骤。}

我在频域中生成了一个Gabor内核。 因为内核已经在频域中，所以我没有对它应用FFT，而图像是FFT编辑的。 然后，我将内核与图像相乘以实现FFT-Convolution。 然后它们被反FFT并像往常一样转换回位图。

产量

1. 内核看起来没问题。 但是，filter输出看起来不是很有希望（或者，它呢？）。
2. 方向（ theta ）对内核没有任何影响。
3. 计算/公式经常在变化值的情况下遭受除零exception。

我该如何解决这些问题？

哦，而且，

• 参数α ， β代表什么？
• 什么应该是f ₀的适当值？

---

更新：

我根据@Cris Luoengo的回答修改了我的代码。

  private double Gabor(double u, double v, double u0, double v0, double a, double b) { double p = (-2) * Math.PI * Math.PI; double q = (u-u0)/(a*a); double r = (v - v0) / (b * b); return Math.Exp(p * (q + r)); } public Array2d GaborKernelFft(int sizeX, int sizeY, double u0, double v0, double a, double b) { double xx = sizeX; double yy = sizeY; double halfX = (xx - 1) / xx; double halfY = (yy - 1) / yy; Array2d kernel = new Array2d(sizeX, sizeY); for (double u = 0; u <= halfX; u += 0.1) { for (double v = 0; v <= halfY; v += 0.1) { double g = Gabor(u, v, u0, v0, a, b); int x = (int)(u * 10); int y = (int)(v * 10); kernel[x,y] = new Complex(g, 0); } } return kernel; } 

哪里，

  double u0 = 0.2; double v0 = 0.2; double alpha = 10;//1.5; double beta = alpha; 

我不确定这是否是一个好的输出。

Answer 1

您找到的Gabor滤波器的等式中似乎存在拼写错误。 Gabor滤波器是频域中的平移高斯滤波器。 因此，它需要在指数中包含u²和v² 。

链接中的等式（2）似乎更明智， 但仍然错过了2 ：

 exp( -2(πσ)² (u-f₀)² ) 

这是1D情况，这是我们想要在θ方向使用的滤波器。 我们现在在垂直方向v上乘以非移位高斯。 我将α和β设置为两个sigma的倒数：

 exp( -2(π/α)² (u-f₀)² ) exp( -2(π/β)² v² ) = exp( -2π²((u-f₀)/α)² + -2π²(v/β)² ) 

你应该用u和v旋转θ来实现上面的等式，就像你已经做的那样。

此外， u和v应该从-0.5到0.5，而不是从-sizeX/2到sizeX/2 。 这就是假设您的FFT在图像中间设置原点，这并不常见。 通常，FFT算法将原点设置在图像的一角。 所以你应该让你的u和v从0运行到(sizeX-1)/sizeX 。

使用上面的校正实现，你应该将f₀ 0设置在0和0.5之间（尝试0.2开始），并且α和β应该足够小，使得高斯不会达到0频率（你想要滤波器到那里是0）

在频域中，您的滤镜看起来像是远离原点的旋转高斯。

在空间域中，滤镜的幅度应该再次像高斯一样。 虚构组件应如下所示（图片链接到我发现的维基百科页面）：

（即它在θ方向上是反对称的（奇数）），可能有更多的波瓣取决于α ， β和f₀ 。 真实组件应该相似但对称（偶数），中间最大。 请注意，在IFFT之后，您可能需要将原点从左上角移动到图像的中间（Google“fftshift”）。

---

请注意，如果将α和β设置为相等，则u – v平面的旋转无关紧要。 在这种情况下，您可以使用笛卡尔坐标而不是极坐标来定义频率。 也就是说，不是将f 0和θ定义为参数，而是定义u 0和v 0。 在指数中，然后用u-u 0代替u-f 0，用v-v 0代替v。

---

编辑问题后的代码再次错过了正方形。 我会写代码如下：

 private double Gabor(double u, double v, double u0, double v0, double a, double b) { double p = (-2) * Math.PI * Math.PI; double q = (u-u0)/a; double r = (v - v0)/b; return Math.Exp(p * (q*q + r*r)); } public Array2d GaborKernelFft(int sizeX, int sizeY, double u0, double v0, double a, double b) { double halfX = sizeX / 2; double halfY = sizeY / 2; Array2d kernel = new Array2d(sizeX, sizeY); for (double y = 0; y < sizeY; y++) { double v = y / sizeY; // v -= HalfY; // whether this is necessary or not depends on your FFT for (double x = 0; x < sizeX; x++) { double u = x / sizeX; // u -= HalfX; // whether this is necessary or not depends on your FFT double g = Gabor(u, v, u0, v0, a, b); kernel[(int)x, (int)y] = new Complex(g, 0); } } return kernel; }