频域中的同态滤波器实现

Question

我找到了 Frequency Domain representation of Homomorphic Filter 的各种公式。我正在使用以下一个：

我按照与 .

相同的源代码模式实现了它

源代码

    private Array2d<Complex> HomoMorphicFilterFft(double M, double N, double yH, double yL, double c, double D0)
    {
        Array2d<Complex> kernel = new Array2d<Complex>((int)M, (int)N);

        for (double y = 0; y < N; y++)
        {
            double v = y / N;

            for (double x = 0; x < M; x++)
            {
                double u = x / M;

                double kw = HMFft(u, v, M, N, yH, yL, c, D0);

                kernel[(int)x, (int)y] = new Complex(kw, 0);
            }
        }

        return kernel;
    }

    private double HMFft(double u, double v, double M, double N, double yH, double yL, double c, double D0)
    {
        double p = u - M / 2;
        double q = v - N / 2;

        double Duv = Math.Sqrt(p * p - q * q);

        double d = (Duv / D0) * (Duv / D0);
        double e = Math.Exp((-1) * c * d);

        double homo = (yH - yL) * (1-e) + yL;

        return homo;
    }
}

内核公式正在生成NaN.

这次我做错了什么？

更新： 我遵循了 Duurt 的回答，但没有输出：

然后我在源码中做了一些修改：

Array2d<double> dOutput = Rescale2d.Rescale(DataConverter2d.ToDouble(cOutput));

替换为

Array2d<double> dOutput = Rescale2d.Limit(DataConverter2d.ToDouble(cOutput));

并且，

Array2d<double> dLimitedKernel = Rescale2d.Limit(dKernel);

被

取代

Array2d<double> dLimitedKernel = Rescale2d.Rescale(dKernel);

但是，输出仍然不是预期的：

预期输出类似于以下内容（或者，是吗？）：

^{Limit()和Rescale()的区别是：Limit() 只修剪那些超出 0-1 范围的值。 Rescale() 通过将数组中的所有值除以数组中的最大值来重新缩放数组中的所有值。} .

源代码

下面是更详细的源码：

public partial class Form1 : Form
{
    public Form1()
    {
        InitializeComponent();

        Bitmap image = DataConverter2d.ReadGray(StandardImage.LenaGray);
        Array2d<double> dImage = DataConverter2d.ToDouble(image);

        int newWidth = Tools.ToNextPowerOfTwo(dImage.Width);
        int newHeight = Tools.ToNextPowerOfTwo(dImage.Height);

        double yH = 2;//2;
        double yL = 0.5;//0.5;
        double c = 0.5;
        double D0 = 1;//0.5;

        Array2d<Complex> kernel2d = HomoMorphicFilterFft(newWidth, newHeight, yH, yL, c, D0);

        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);
        // trims the values to keep them between 0 and 1.
        Array2d<double> dOutput = Rescale2d.Limit(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);
        // Rescales the values to keep them between 0 and 1.
        Array2d<double> dLimitedKernel = Rescale2d.Rescale(dKernel);

        Bitmap kernel = DataConverter2d.ToBitmap(dLimitedKernel, image.PixelFormat);

        pictureBoxExt1.Image = image;
        pictureBoxExt2.Image = kernel;
        pictureBoxExt3.Image = output;
    }

    private Array2d<Complex> HomoMorphicFilterFft(double M, double N, double yH, double yL, double c, double D0)
    {
        Array2d<Complex> kernel = new Array2d<Complex>((int)M, (int)N);

        for (double y = 0; y < N; y++)
        {
            double v = y / N;

            for (double x = 0; x < M; x++)
            {
                double u = x / M;

                double kw = HMFft(u, v, M, N, yH, yL, c, D0);

                kernel[(int)x, (int)y] = new Complex(kw, 0);
            }
        }

        return kernel;
    }

    private double HMFft(double u, double v, double M, double N, double yH, double yL, double c, double D0)
    {
        double p = u - M / 2;
        double q = v - N / 2;

        double Duv = Math.Sqrt(p * p + q * q);

        double d = (Duv / D0) * (Duv / D0);
        double e = Math.Exp((-1) * c * d);

        double homo = (yH - yL) * (1-e) + yL;

        return homo;
    }
}

^{此时只专注于算法}

Answer 1

Duv 在 Sqrt 中有一个负号，而在公式中它有一个加号。取负数的平方可以解释你的问题。

频域中的同态滤波器实现

Homomorphic Filter implementation in Frequency domain

c#

image-processing

imagefilter