Помогите с реализацией виннеровского фильтра

Підписуйтеся на Telegram-канал «DOU #tech», щоб не пропустити нові технічні статті

Собственно сама формула фильтра F(u,v)=conj(H(u,v)) / ( |H(u,v)|^2 + K ) * G(u,v)
где G — двухмерное преобразование Фурье исходной картинки
H — искажающая функция.
en.wikipedia.org/...​wiki/Wiener_deconvolution

        public static Complex[,] WienerFilter(Complex[,] Gx, Complex[,] Hx)
        {
            int w = Gx.GetLength(0);
            int h = Gx.GetLength(1);
            Complex[,] lRes = new Complex[w, h];
            //F(u,v)=conj(H(u,v)) / ( |H(u,v)|^2 + K ) * G(u,v)            
            for (int ly = 0; ly < h; ly++)
                for (int lx = 0; lx < w; lx++)
                {
                    Complex lH = Hx[lx, ly];
                    Complex lG = Gx[lx, ly];
                    double lHmod = lH.Re * lH.Re + lH.Im * lH.Im;
                    Complex lHConj = lH;
                    lHConj.Im = -lHConj.Im;
                    double lDiv = (lHmod + 0.003f);

                    Complex lFilt = lHConj;
                    lFilt.Im = lFilt.Im / lDiv;
                    lFilt.Re = lFilt.Re / lDiv;
                    lFilt.Re = Math.Abs(lFilt.Re);
                    lFilt.Im = Math.Abs(lFilt.Im);

                    Complex lFiltRes = lG * lFilt;

                    lRes[lx, ly] = lFiltRes;                
                }

            return lRes;
        }

0.003f — это параметр K отношение шума к сигналу
Gx — оригинальная картинка после преобразования фурье.
Hx — двухмерная гаусова свертка после преобразования Фурье
Нормализация на количество отсчетов 1/N производится при прямом преобразовании Фурье.

Если не делать
lFilt.Re = Math.Abs(lFilt.Re);
lFilt.Im = Math.Abs(lFilt.Im);
То исходная картинка отражается по 4 квадрантам изображения.

Вспомогательные функции

        public static ByteBitmap GetGausianByteBitmapImage(Size OutImageSize, int HalfWindow, float aPower, float aSigma)
        {
            ByteBitmap lRes = null;

            int lMinHalf = Math.Min(OutImageSize.Width, OutImageSize.Height) / 2;
            HalfWindow = Math.Min(HalfWindow, lMinHalf);

            ByteBitmap Gausian = new ByteBitmap(2 * HalfWindow, 2 * HalfWindow);
            float lCenterX = (float)HalfWindow - 0.5f;
            float lCenterY = (float)HalfWindow - 0.5f;
            //
            //-1/(2 * aSigma ^ 2);
            float larg = -1 / (2 * aSigma * aSigma);
            for (int ly = 0; ly < Gausian.Height; ly++)
                for (int lx = 0; lx < Gausian.Width; lx++)
                {
                    float xdiff = (float)lx - lCenterX;
                    float ydiff = (float)ly - lCenterY;
                    Gausian[lx, ly] = (byte)(aPower * Math.Exp(larg * (xdiff * xdiff + ydiff * ydiff)));
                }

            //fit for res
            lRes = new ByteBitmap(OutImageSize.Width, OutImageSize.Height);
            int lNewShiftX = (OutImageSize.Width - Gausian.Width) / 2;
            int lNewShiftY = (OutImageSize.Height - Gausian.Height) / 2;
            for (int ly = 0; ly < Gausian.Height; ly++)
                for (int lx = 0; lx < Gausian.Width; lx++)
                {
                    lRes[lx + lNewShiftX, ly + lNewShiftY] = Gausian[lx, ly];

                }

            return lRes;
        }

        public static ByteBitmap FitByteBitmapForFFT2D(ByteBitmap aBitmap, int Shift, out Point StartPoint, byte DefVal)
        {
            StartPoint = new Point(0, 0);
            ByteBitmap lRes = null;
            //calculate new shift
            int lNewWidth = aBitmap.Width + 2 * Shift;
            int lNewHeight = aBitmap.Height + 2 * Shift;
            lNewWidth += AdditionToPowerOfTwo(lNewWidth);
            lNewHeight += AdditionToPowerOfTwo(lNewHeight);
            int lNewShiftX = (lNewWidth - aBitmap.Width) / 2;
            int lNewShiftY = (lNewHeight - aBitmap.Height) / 2;


            if (0 != DefVal)
                lRes = new ByteBitmap(lNewWidth, lNewHeight, DefVal);
            else
                lRes = new ByteBitmap(lNewWidth, lNewHeight);

            for (int ly = 0; ly < aBitmap.Height; ly++)
                for (int lx = 0; lx < aBitmap.Width; lx++)
                    lRes[lx + lNewShiftX, ly + lNewShiftY] = aBitmap[lx, ly];
            //
            StartPoint.X = lNewShiftX;
            StartPoint.Y = lNewShiftY;

            return lRes;
        }

        public static ByteBitmap TestWienerFilter(ByteBitmap aBitmap)
        {
            Point lStart;
            ByteBitmap lLoc = FFTHelpers.FitByteBitmapForFFT2D(aBitmap, 10, out lStart, 255);
            Complex[,] lGx = FFTHelpers.ByteBitmapToComplex(lLoc);
            FourierTransform.FFT2(lGx, FourierTransform.Direction.Forward);

            //convolution function
            //gauss distribution
            ByteBitmap lBt = GetGausianByteBitmapImage(new Size(lGx.GetLength(0), lGx.GetLength(1)), 5, 255, 1.0f);
            Complex[,] lHx = FFTHelpers.ByteBitmapToComplex(lBt);
            FourierTransform.FFT2(lHx, FourierTransform.Direction.Forward);

            //make wiener
            Complex[,] lWR = WienerFilter(lGx, lHx);
            
            //backward result
            FourierTransform.FFT2(lWR, FourierTransform.Direction.Backward);


            //convert to bytebitmap
            ByteBitmap lRes = FFTHelpers.ComplexToByteBitmap(lWR);
            lRes = FFTHelpers.GetFittedByteBitmapFromFFT2D(lRes, lStart, new Size(aBitmap.Width, aBitmap.Height));
            return lRes;
        }

Так вот в результате TestWienerFilter фильтр почему-то делает не восстановление сигнала после применения функции, а наоборот размазывает сигнал Hx — гаусовым блюром.

Проект для сборки(VS2008 но легко конвертируется и в 2010) можно найти здесь:
yadi.sk/d/lJ_sQbss4clws

PS
Оригинальный фильтр F = G (1/H) * (|H|^2/(|H|^2 + SNR)) вообще даёт полностью черную картинку.