Gabor滤波器在频域中的实现 [英] Gabor Filter implementation in Frequency domain
问题描述
此处 我们有Gabor过滤器的Spatial域实现.但是，出于性能原因，我需要在频域中实现Gabor过滤器.
我找到了 Gabor滤波器的频域方程:
我实际上对该公式的正确性和/或适用性有疑问.
源代码
Source Code
因此，我实现了以下内容:
So, I have implemented the following :
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;
}
}
这时仅关注算法步骤.
我在频域中生成了一个Gabor内核.由于内核已经在频域中，因此我没有对其应用FFT，而对图像进行了FFT处理.然后，我将内核和图像相乘以实现FFT卷积.然后将它们进行逆FFT并照常转换回Bitmap.
I have generated a Gabor kernel in the frequency domain. Since, the kernel is already in Frequency domain, I didn't apply FFT to it, whereas image is FFT-ed. Then, I multiplied the kernel and the image to achieve FFT-Convolution. Then they are inverse-FFTed and converted back to Bitmap as usual.
输出
Output
- 内核看起来还不错.但是，过滤器的输出看起来并不十分有前景(或者，是吗?).
- 方向( theta )对内核没有任何影响.
- 计算/公式经常会因数值变化而被零除异常困扰.
- The kernel looks okay. But, The filter-output doesn't look very promising (or, does it?).
- The orientation (theta) doesn't have any effect on the kernel.
- The calculation/formula is frequently suffering from divide-by-zero exception up on changing values.
如何解决这些问题?
哦，还有
- 参数α，β代表什么?
- f 0 的合适值是什么?
- what do the parameters α, β, represent?
- what should be the appropriate value of f0?
更新:
我已根据 @Cris Luoengo 的答案修改了代码.
I have modified my code according to @Cris Luoengo's answer.
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;
}
哪里
double u0 = 0.2;
double v0 = 0.2;
double alpha = 10;//1.5;
double beta = alpha;
我不确定这是否是一个好输出.
I am not sure whether this is a good output.
推荐答案
您发现的Gabor过滤器的方程式中似乎有一个错字.伽柏滤波器在频域中是转换后的高斯.因此，它需要在指数中包含u²和v².
There seems to be a typo in the equation for the Gabor filter that you found. The Gabor filter is a translated Gaussian in the frequency domain. Hence, it needs to have u² and v² in the exponent.
您链接中的方程式(2)似乎更明智，但仍然缺少2 :
Equation (2) in your link seems more sensible, but still misses a 2:
exp( -2(πσ)² (u-f₀)² )
这是一维情况，这是我们要在θ方向上使用的滤波器.现在，我们在垂直方向v上乘以不变的高斯.我将α和β设置为两个西格玛的倒数:
It is the 1D case, this is the filter we want to use in the direction θ. We now multiply in the perpendicular direction, v, with a non-shifted Gaussian. I set α and β to be the inverse of the two sigmas:
exp( -2(π/α)² (u-f₀)² ) exp( -2(π/β)² v² ) = exp( -2π²((u-f₀)/α)² + -2π²(v/β)² )
您应该像以前一样在u和v旋转了θ的情况下实现上述方程式.
You should implement the above equation with u and v rotated over θ, as you already do.
此外，u和v应该从-0.5到0.5，而不是-sizeX/2到sizeX/2.这就是假设您的FFT将原点设置在图像的中间，这并不常见.通常，FFT算法将原点设置在图像的一角.因此，您应该将u和v的值从0改为(sizeX-1)/sizeX.
Also, u and v should run from -0.5 to 0.5, not from -sizeX/2 to sizeX/2. And that is assuming your FFT sets the origin in the middle of the image, which is not common. Typically the FFT algorithms set the origin in a corner of the image. So you should probably have your u and v run from 0 to (sizeX-1)/sizeX instead.
使用上述更正的实现，您应该将f₀设置为0到0.5之间(尝试从0.2开始)，并且α和β应该足够小，以使高斯不能达到0频率(您希望滤波器在那里为0)
With a corrected implementation as above, you should set f₀ to be between 0 and 0.5 (try 0.2 to start with), and α and β should be small enough such that the Gaussian doesn't reach the 0 frequency (you want the filter to be 0 there)
在频域中，您的滤波器看起来像是远离原点的旋转高斯.
In the frequency domain, your filter will look like a rotated Gaussian away from the origin.
在空间域中，滤波器的幅度应再次看起来像高斯.虚构的组件应如下所示(图片链接指向我在其上找到的Wikipedia页面):
In the spatial domain, the magnitude of your filter should look again like a Gaussian. The imaginary component should look like this (picture links to Wikipedia page I found it on):
(即它在&θ;方向上是反对称的(奇数))，可能有更多的波瓣，具体取决于α，β和f₀.实数部分应相似但对称(偶数)，中间最大.请注意，在进行IFFT之后，您可能需要将原点从图像的左上角移到图像的中间(Google"fftshift").
(i.e. it is anti-symmetric (odd) in the θ direction), possibly with more lobes depending on α, β and f₀. The real component should be similar but symmetric (even), with a maximum in the middle. Note that after IFFT, you might need to shift the origin from the top-left corner to the middle of the image (Google "fftshift").
请注意，如果将α和β设置为相等，则u-v平面的旋转无关紧要.在这种情况下，可以使用直角坐标而不是极坐标来定义频率.也就是说，您无需定义f₀和θ作为参数，而是定义u₀和v₀.然后在指数中将u-f₀替换为u-u₀，将v替换为v-v₀.
Note that if you set α and β to be equal, the rotation of the u-v plane is irrelevant. In this case, you can use cartesian coordinates instead of polar coordinates to define the frequency. That is, instead of defining f₀ and θ as parameters, you define u₀ and v₀. In the exponent you then replace u-f₀ with u-u₀, and v with v-v₀.
问题编辑后的代码再次丢失正方形.我将编写如下代码:
The code after the edit of the question misses the square again. I would write the code as follows:
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;
}
这篇关于Gabor滤波器在频域中的实现的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
