Here我们有 Gabor 滤波器的空间域实现。但是，我需要在频域中实现 Gabor 滤波器出于性能原因。
我已经找到了Gabor滤波器的频域方程:
I am actually in doubt about the correctness and/or applicability of this formula.
源代码
所以，我实施了以下措施：
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<Complex> GaborKernelFft(int sizeX, int sizeY, double f0, double theta, double a, double b)
{
int halfX = sizeX / 2;
int halfY = sizeY / 2;
Array2d<Complex> kernel = new Array2d<Complex>(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<double> 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<Complex> kernel2d = GaborKernelFft(newWidth, newHeight, u0, v0, alpha, beta);
dImage.PadTo(newWidth, newHeight);
Array2d<Complex> cImage = DataConverter2d.ToComplex(dImage);
Array2d<Complex> fImage = FourierTransform.ForwardFft(cImage);
// FFT convolution .................................................
Array2d<Complex> fOutput = new Array2d<Complex>(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<Complex> cOutput = FourierTransform.InverseFft(fOutput);
Array2d<double> dOutput = Rescale2d.Rescale(DataConverter2d.ToDouble(cOutput));
//dOutput.CropBy((newWidth-image.Width)/2, (newHeight - image.Height)/2);
Bitmap output = DataConverter2d.ToBitmap(dOutput, image.PixelFormat);
Array2d<Complex> cKernel = FourierTransform.InverseFft(kernel2d);
cKernel = FourierTransform.RemoveFFTShift(cKernel);
Array2d<double> dKernel = DataConverter2d.ToDouble(cKernel);
Array2d<double> dRescaledKernel = Rescale2d.Rescale(dKernel);
Bitmap kernel = DataConverter2d.ToBitmap(dRescaledKernel, image.PixelFormat);
pictureBox1.Image = image;
pictureBox2.Image = kernel;
pictureBox3.Image = output;
}
}
Just concentrate on the algorithmic steps at this time.
我已经在频域中生成了 Gabor 核。由于内核已经在频域中，因此我没有对其应用 FFT，而图像是经过 FFT 处理的。然后，我将内核与图像相乘以实现 FFT 卷积。然后它们像往常一样被逆 FFT 并转换回位图。
Output
我该如何解决这些问题？
哦，还有，
Update:
我已经根据修改了我的代码@克里斯·罗恩戈的回答。
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<Complex> 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<Complex> kernel = new Array2d<Complex>(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;
}
where,
double u0 = 0.2;
double v0 = 0.2;
double alpha = 10;//1.5;
double beta = alpha;
我不确定这是否是一个好的输出。
您发现的 Gabor 滤波器的方程似乎存在拼写错误。 Gabor 滤波器是频域中的平移高斯滤波器。因此，它需要有u² and v²在指数中。
您链接中的方程（2）似乎更合理，但还是错过了2:
exp( -2(πσ)² (u-f₀)² )
这是一维情况，这是我们要在 θ 方向使用的滤波器。我们现在在垂直方向相乘，v，具有非平移高斯分布。我设置α and β是两个西格玛的倒数：
exp( -2(π/α)² (u-f₀)² ) exp( -2(π/β)² v² ) = exp( -2π²((u-f₀)/α)² + -2π²(v/β)² )
您应该使用以下命令来实现上面的等式u and v旋转 θ，就像你已经做的那样。
Also, u and v应该从 -0.5 到 0.5 运行，而不是从-sizeX/2 to sizeX/2。假设您的 FFT 将原点设置在图像的中间，这并不常见。通常，FFT 算法将原点设置在图像的一角。所以你可能应该有你的u and v从 0 运行到(sizeX-1)/sizeX反而。
通过上述更正的实现，您应该设置f₀介于 0 和 0.5 之间（尝试从 0.2 开始），并且α and β应足够小，以便高斯不会达到 0 频率（您希望滤波器在那里为 0）
在频域中，您的滤波器看起来就像远离原点的旋转高斯函数。
在空间域中，滤波器的幅度应该再次看起来像高斯分布。假想部分应该如下所示（图片链接到我在其中找到的维基百科页面）：
（即，它在 θ 方向上是反对称（奇数）），可能有更多的波瓣，具体取决于α, β and f₀。实数分量应该相似但对称（偶数），最大值位于中间。请注意，在 IFFT 之后，您可能需要将原点从左上角移动到图像的中间（Google“fftshift”）。
请注意，如果您设置α and β为相等，则旋转u-v飞机无关。在这种情况下，您可以使用笛卡尔坐标而不是极坐标来定义频率。也就是说，您不是将 f₀ 和 θ 定义为参数，而是定义 u₀ 和 v₀。然后在指数中将 u-f₀ 替换为 u-u₀，将 v 替换为 v-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<Complex> GaborKernelFft(int sizeX, int sizeY, double u0, double v0, double a, double b)
{
double halfX = sizeX / 2;
double halfY = sizeY / 2;
Array2d<Complex> kernel = new Array2d<Complex>(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;
}