频域同态滤波器的实现

我找到了各种计算公式 Frequency Domain representation of Homomorphic Filter . 我正在使用以下方法:

哪里, D(u, v) is

我已经按照与 FFT Gabor 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 .

在这种情况下我做错了什么?

我听从了杜尔特的回答,结果却不见了:

然后我对源代码做了一些修改:

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() 是: 仅修剪超出0-1范围的值。 重缩放()} .

源代码

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;
    }
}

^{现在只关注算法。}