频域中的Gabor滤波器实现

Question

这里我们有Gabor滤波器的空间域实现。但是，出于性能原因，我需要在频域中实现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<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并像往常一样转换回位图。

输出

1. 核看起来还不错。但是，滤波器输出看起来并不十分理想（或者说它是吗？）。
2. 方向（theta）对核没有任何影响。
3. 在更改值时，计算/公式经常遭受除以零的异常情况。

我该如何解决这些问题？

哦，还有,

• 参数α、β表示什么？
• 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<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;

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

- user366312

---

中间的图是您在频域中的 Gabor 滤波器吗？它应该看起来像一个平移的高斯函数（即远离原点，即左上角像素）。在空间域中，它将具有复杂值。过滤后的结果也是如此。您是否显示了幅度？ - Cris Luengo

另外，出于性能考虑，您应该在空间域中实现递归算法：https://www.researchgate.net/publication/3318442_Recursive_Gabor_filtering--这比FFT快得多。 - Cris Luengo

@CrisLuengo，中间的图是您在频域中的Gabor滤波器吗？ --- 是的。 您显示了幅度吗？ --- 是的。 - user366312

从 f0= sizeX/4 开始。这应该给你一个介于0和奈奎斯特频率之间的高斯函数。Alpha和beta是高斯函数的标准差，它们越小，你对频率就越有选择性，对空间的选择性就越小。如果 sizeX 和 sizeY 是256，那么你的10和15相当大。这没关系，但如果结果仍然不好，请尝试使用较小的值。 - Cris Luengo

Gabor是一种可适应的滤波器，就像高斯、FFT和许多其他滤波器一样。您可以将其应用于每个图像行，然后在结果上，在每个图像列上进行操作。在这种情况下，可分离的含义是您可以将空间域滤波器编写为两个1D滤波器的外积。 - Cris Luengo

Answer 1

你找到的Gabor滤波器方程似乎有一个打字错误。Gabor滤波器是频域中的平移高斯滤波器。因此，它需要在指数中具有u²和v²。

你链接中的方程（2）似乎更合理，但仍然缺少一个2：

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

这是1D情况，我们想要在方向θ上使用这个滤波器。现在我们在垂直方向v上用一个未移位的高斯函数进行乘法。我将α和β设置为两个标准差的倒数。

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.5 之间（尝试从0.2开始），并且 α 和 β 应该足够小，以使高斯函数不到达 0 频率（你希望滤波器在那里为 0）。

在频域中，你的滤波器看起来像是一个远离原点的旋转高斯函数。

在空间域中，您的滤波器的幅度应该再次类似于高斯分布。虚部应该像这样（图片链接到我找到的维基百科页面）：

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

---

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

- Cris Luengo

---

你完全移除了sin(θ)和cos(θ)吗？你的公式中没有三角函数项吗？ - user366312

@匿名用户：我简化了。将此答案中的 u 和 v 视为您代码中的 uDash 和 vDash。也就是说，在使用我提供的方程之前，您需要旋转 u-v 坐标轴。 - Cris Luengo

尽管如此，我总是使用α==β，然后旋转就不是必要的了，一切都变得更加容易。然后，您可以将f₀和θ定义为参数，而将u₀和v₀定义为参数。在指数中，您可以用u-u₀代替u-f₀，用v-v₀代替v。也就是说，使用笛卡尔坐标系来定义频率，而不是使用极坐标系。 - Cris Luengo

这个内核似乎无法旋转。 - user366312

旋转来自于选择 u₀ 和 v₀。 - Cris Luengo