Search code examples
c#emgucvgaussianaforgefrequency-domain

Applying Gaussian blur to image in frequency domain


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:

  1. Load the image.
  2. 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();
    }
    
  3. 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;
    }
    
  4. Apply Gaussian blur.

    1. 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))));
      }
      
    2. Apply FFT to every channel and kernel.

    3. 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]);
              }
          }
      }
      
  5. Apply IFFT to every channel.
  6. 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;
    }
    
  7. 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.


Solution

  • 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:

    1. I've found how kernel should look like after FFT (see below).
    2. Looked up colors of that image.
    3. Calculated gauss2d for x = -2; y = -2; std = 1.
    4. Calculated the prescaler to receive color value from value calculated in pt. 3 (see wolfram).
    5. 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):

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