Реализация гомоморфного фильтра в частотной области

Я нашел различные формулы для представления гомоморфного фильтра в частотной области. Я использую следующий:

Где D(u, v) равно:

Я реализовал его по тому же шаблону исходного кода, что и 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() заключается в следующем: 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;
    }
}

^{Сейчас сосредоточьтесь только на алгоритме.}