Applying Gaussian blur to image in frequency domain

Question

I've got torubles with appling gaussian blur to image in frequency domain. For unknown reasons (probably I've dont something wrong) I recieve wired image instead of blurred one.

There's what i do step by step:

Load the image.

Split image into separate channels.

private static Bitmap[] separateColorChannels(Bitmap source, int channelCount)
{
    if (channelCount != 3 && channelCount != 4)
    {
        throw new NotSupportedException("Bitmap[] FFTServices.separateColorChannels(Bitmap, int): Only 3 and 4 channels are supported.");
    }

    Bitmap[] result = new Bitmap[channelCount];
    LockBitmap[] locks = new LockBitmap[channelCount];
    LockBitmap sourceLock = new LockBitmap(source);
    sourceLock.LockBits();

    for (int i = 0; i < channelCount; ++i)
    {
        result[i] = new Bitmap(source.Width, source.Height, PixelFormat.Format8bppIndexed);
        locks[i] = new LockBitmap(result[i]);
        locks[i].LockBits();
    }

    for (int x = 0; x < source.Width; x++)
    {
        for (int y = 0; y < source.Height; y++)
        {
            switch (channelCount)
            {
                case 3:
                    locks[0].SetPixel(x, y, Color.FromArgb(sourceLock.GetPixel(x, y).R));
                    locks[1].SetPixel(x, y, Color.FromArgb(sourceLock.GetPixel(x, y).G));
                    locks[2].SetPixel(x, y, Color.FromArgb(sourceLock.GetPixel(x, y).B));

                    break;
                case 4:
                    locks[0].SetPixel(x, y, Color.FromArgb(sourceLock.GetPixel(x, y).A));
                    locks[1].SetPixel(x, y, Color.FromArgb(sourceLock.GetPixel(x, y).R));
                    locks[2].SetPixel(x, y, Color.FromArgb(sourceLock.GetPixel(x, y).G));
                    locks[3].SetPixel(x, y, Color.FromArgb(sourceLock.GetPixel(x, y).B));

                    break;
                default:
                    break;
            }
        }
    }

    for (int i = 0; i < channelCount; ++i)
    {
        locks[i].UnlockBits();
    }

    sourceLock.UnlockBits();
}

Convert every channel into complex images (with AForge.NET).

public static AForge.Imaging.ComplexImage[] convertColorChannelsToComplex(Emgu.CV.Image<Emgu.CV.Structure.Gray, Byte>[] channels)
{
    AForge.Imaging.ComplexImage[] result = new AForge.Imaging.ComplexImage[channels.Length];

    for (int i = 0; i < channels.Length; ++i)
    {
        result[i] = AForge.Imaging.ComplexImage.FromBitmap(channels[i].Bitmap);
    }

    return result;
}

Apply Gaussian blur.

First i create the kernel (For testing purposes kernel size is equal to image size, tho only center part of it is calculated with gaussian function, rest of kernel is equal to re=1 im=0).

private ComplexImage makeGaussKernel(int side, double min, double max, double step, double std)
{
    // get value at top left corner
    double _0x0 = gauss2d(min, min, std);

    // top left corner should be 1, so making scaler for rest of the values
    double scaler = 1 / _0x0;

    int pow2 = SizeServices.getNextNearestPowerOf2(side);

    Bitmap bitmap = new Bitmap(pow2, pow2, PixelFormat.Format8bppIndexed);

    var result = AForge.Imaging.ComplexImage.FromBitmap(bitmap);

    // For test purposes my kernel is size of image, so first, filling with 1 only.
    for (int i = 0; i < result.Data.GetLength(0); ++i)
    {
        for (int j = 0; j < result.Data.GetLength(0); ++j)
        {
            result.Data[i, j].Re = 1;
            result.Data[i, j].Im = 0;
        }
    }

    // The real kernel's size.
    int count = (int)((Math.Abs(max) + Math.Abs(min)) / step);

    double h = min;
    // Calculating kernel's values and storing them somewhere in the center of kernel.
    for (int i = result.Data.GetLength(0) / 2 - count / 2; i < result.Data.GetLength(0) / 2 + count / 2; ++i)
    {
        double w = min;
        for (int j = result.Data.GetLength(1) / 2 - count / 2; j < result.Data.GetLength(1) / 2 + count / 2; ++j)
        {
            result.Data[i, j].Re = (scaler * gauss2d(w, h, std)) * 255;
            w += step;
        }
        h += step;
    }

    return result;
}

// The gauss function
private double gauss2d(double x, double y, double std)
{
    return ((1.0 / (2 * Math.PI * std * std)) * Math.Exp(-((x * x + y * y) / (2 * std * std))));
}

Apply FFT to every channel and kernel.

Multiply center part of every channel by kernel.

void applyFilter(/*shortened*/)
{
    // Image's size is 512x512 that's why 512 is hardcoded here
    // min = -2.0; max = 2.0; step = 0.33; std = 11
    ComplexImage filter = makeGaussKernel(512, min, max, step, std);

    // Applies FFT (with AForge.NET) to every channel and filter
    applyFFT(complexImage);
    applyFFT(filter);

    for (int i = 0; i < 3; ++i)
    {
        applyGauss(complexImage[i], filter, side);
    }

    // Applies IFFT to every channel
    applyIFFT(complexImage);
}

private void applyGauss(ComplexImage complexImage, ComplexImage filter, int side)
{
    int width = complexImage.Data.GetLength(1);
    int height = complexImage.Data.GetLength(0);

    for(int i = 0; i < height; ++i)
    {
        for(int j = 0; j < width; ++j)
        {
            complexImage.Data[i, j] = AForge.Math.Complex.Multiply(complexImage.Data[i, j], filter.Data[i, j]);
        }
    }
}

Apply IFFT to every channel.

Convert every channel back to bitmaps (with AForge.NET).

public static System.Drawing.Bitmap[] convertComplexColorChannelsToBitmap(AForge.Imaging.ComplexImage[] channels)
{
    System.Drawing.Bitmap[] result = new System.Drawing.Bitmap[channels.Length];

    for (int i = 0; i < channels.Length; ++i)
    {
        result[i] = channels[i].ToBitmap();
    }

    return result;
}

Merge bitmaps into single bitmap

public static Bitmap mergeColorChannels(Bitmap[] channels)
{
    Bitmap result = null;

    switch (channels.Length)
    {
        case 1:
            return channels[0];
        case 3:
            result = new Bitmap(channels[0].Width, channels[0].Height, PixelFormat.Format24bppRgb);
            break;
        case 4:
            result = new Bitmap(channels[0].Width, channels[0].Height, PixelFormat.Format32bppArgb);
            break;
        default:
            throw new NotSupportedException("Bitmap FFTServices.mergeColorChannels(Bitmap[]): Only 1, 3 and 4 channels are supported.");
    }

    LockBitmap resultLock = new LockBitmap(result);
    resultLock.LockBits();

    LockBitmap red = new LockBitmap(channels[0]);
    LockBitmap green = new LockBitmap(channels[1]);
    LockBitmap blue = new LockBitmap(channels[2]);

    red.LockBits();
    green.LockBits();
    blue.LockBits();

    for (int y = 0; y < result.Height; y++)
    {
        for (int x = 0; x < result.Width; x++)
        {
            resultLock.SetPixel(x, y, Color.FromArgb((int)red.GetPixel(x, y).R, (int)green.GetPixel(x, y).G, (int)blue.GetPixel(x, y).B));
        }
    }

    red.UnlockBits();
    green.UnlockBits();
    blue.UnlockBits();

    resultLock.UnlockBits();

    return result;
}

As a result I've got shifted, red-colored blurred version of image: link.

@edit - Updated the question with several changes to the code.

Is the result identical to the original image, or identical to the weird image? If the kernel of all 1s gives the weird image, then you know the error is not in your blurring code. So it must be in the FFT/IFFT or in the process of splitting the channels and recombining them, or in the process of generating the complex image with aforge.net. On the other hand, if the kernel of all 1s gives the original image back, then the problem lies in the blurring/kernel. Does this line of thinking seem to help? — Moby Disk
the result is weird image, sorry for inaccurate reply. It helps a bit, thanks. However question is still open if someone wants to add anything. — user2475983
I've updated the question with several changes in the code and changes in resulted image. — user2475983

user2475983 user2475983 · Accepted Answer · 2015-04-12T16:46:55

I figured it out with some help at DSP stackexchange... and some cheating but it works. The main problem was kernel generation and applying FFT to it. Also important thing is that AForge.NET divides image pixels by 255 during conversion to ComplexImage and multiplies by 255 during conversion from ComplexImage to Bitmap (thanks Olli Niemitalo @ DSP SE).

How I solved this:

I've found how kernel should look like after FFT (see below).
Looked up colors of that image.
Calculated gauss2d for x = -2; y = -2; std = 1.
Calculated the prescaler to receive color value from value calculated in pt. 3 (see wolfram).
Generated kernel with scaled values with perscaler from pt. 4.

However I cant use FFT on generated filter, because generated filter looks like filter after FFT already. It works - the output image is blurred without artifacts so I think that's not too bad.

The images (I cant post more than 2 links, and images are farily big):

Input image
Generated filter (without FFT!)
- Parameters for below function:
- std = 1.0
- size = 8.0
- width = height = 512
Result image

The final code:

private ComplexImage makeGaussKernel(double size, double std, int imgWidth, int imgHeight)
{
    double scale = 2000.0;
    double hsize = size / 2.0;

    Bitmap bmp = new Bitmap(imgWidth, imgHeight, PixelFormat.Format8bppIndexed);
    LockBitmap lbmp = new LockBitmap(bmp);

    lbmp.LockBits();

    double y = -hsize;
    double yStep = hsize / (lbmp.Height / 2.0);
    double xStep = hsize / (lbmp.Width / 2.0);

    for (int i = 0; i < lbmp.Height; ++i)
    {
        double x = -hsize;

        for (int j = 0; j < lbmp.Width; ++j)
        {
            double g = gauss2d(x, y, std) * scale;

            g = g < 0.0 ? 0.0 : g;
            g = g > 255.0 ? 255.0 : g;

            lbmp.SetPixel(j, i, Color.FromArgb((int)g));

            x += xStep;
        }

        y += yStep;
    }

    lbmp.UnlockBits();

    return ComplexImage.FromBitmap(bmp);
}

private double gauss2d(double x, double y, double std)
{
    return (1.0 / (2 * Math.PI * std * std)) * Math.Exp(-(((x * x) + (y * y)) / (2 * std * std)));
}

private void applyGaussToImage(ComplexImage complexImage, ComplexImage filter)
{
    for (int i = 0; i < complexImage.Height; ++i)
    {
        for (int j = 0; j < complexImage.Width; ++j)
        {
            complexImage.Data[i, j] = AForge.Math.Complex.Multiply(complexImage.Data[i, j], filter.Data[i, j]);
        }
    }
}

Applying Gaussian blur to image in frequency domain

1 Answers