Fast Image Convolution Using Radial Symmetry (with source)

6 posts / 0 new
Last post
TheVinn
Offline
Last seen: 6 months 2 weeks ago
Joined: 29 Aug 2009 - 11:31
Fast Image Convolution Using Radial Symmetry (with source)

Okay so I was bored and burned out so I took a break and whipped out this class I've been meaning to create.

Taking advantage of the separable nature of radially symmetric filters (http://en.wikipedia.org/wiki/Gaussian_blur) I wrote my own version of an image convolver that can apply a Gaussian blur to an image. This one correctly performs edge replication of color channel information to give values to the undefined areas outside the image that must necessarily be sampled in order to properly apply the filter. There are two routines of interest. createConvolvedImageFull() will blur the source image, producing a new image that is larger by 2*radius-1 pixels in each dimension. createConvolvedImage() will return a convolved image that with the same size and origin as the source image. Note that in this case, some data is lost.

These routines work with any type of image. If you provide an image with an associated alpha channel, it will correctly blur the alpha channel as well (treating undefined areas in the source as if they were fully transparent).

// Takes advantage of radial symmetry to implement an efficient image convolution
class RadialImageConvolutionKernel
{
public:
  // This includes the center sample
  explicit RadialImageConvolutionKernel (int radiusInSamples);
  ~RadialImageConvolutionKernel ();

  int getRadiusInSamples () { return m_radius; }

  // the gaussian radius will be (radiusInSamples-1)/2
  void createGaussianBlur ();

  void setOverallSum (float desiredTotalSum);

  void rescaleAllValues (float multiplier);
  
  // creates a convolved image with the same dimensions as the source
  Image createConvolvedImage (const Image& sourceImage) const;

  // creates a larger image that fully contains the result
  // of applying the convolution kernel to the source image.
  Image createConvolvedImageFull (const Image& sourceImage) const;

  //
  // copy a line of color components, performing edge
  // replication on either side. dest must have 2 * replicas
  // more components than src.
  //
  static void copy (int pixels,
                    uint8* dest,
                    int destSkip,
                    const uint8* src,
                    int srcSkip,
                    int replicas);

  //
  // copy a line of alpha values, inserting fully transparent
  // values (0) as replicas on either side
  //
  static void copy_alpha (int pixels,
                          uint8* dest,
                          int destSkip,
                          const uint8* src,
                          int srcSkip,
                          int replicas);

  //
  // convolve a line of color components.
  // src must have 2 * (radius - 1) more components than dest
  //
  static void convolve (int pixels,
                        uint8* dest,
                        int destSkip,
                        const uint8* src,
                        int srcSkip,
                        int radius,
                        const float* kernel);

private:
  const int m_radius;
  HeapBlock<float> m_kernel;
};

RadialImageConvolutionKernel::RadialImageConvolutionKernel (int radiusInSamples)
  : m_radius (radiusInSamples)
{
  jassert (radiusInSamples >= 2);

  m_kernel.malloc (radiusInSamples);
}

RadialImageConvolutionKernel::~RadialImageConvolutionKernel()
{
}

void RadialImageConvolutionKernel::createGaussianBlur ()
{
  const double blurRadius = double(m_radius-1) / 2;
  const double radiusFactor = -1 / (2 * blurRadius * blurRadius);

  for (int r = 0; r < m_radius; ++r)
    m_kernel[r] = float( exp( radiusFactor * (r * r)));

  setOverallSum (1.0f);
}

void RadialImageConvolutionKernel::setOverallSum (float desiredTotalSum)
{
  double currentTotal = 0.0;

  for (int r = m_radius; --r >= 1; )
    currentTotal += m_kernel[r];
  currentTotal *= 2; // for symmetry
  currentTotal += m_kernel[0];

  rescaleAllValues (float (desiredTotalSum / currentTotal));
}

void RadialImageConvolutionKernel::rescaleAllValues (float multiplier)
{
  for (int r = m_radius; --r >= 0; )
    m_kernel[r] *= multiplier;
}

Image RadialImageConvolutionKernel::createConvolvedImage (const Image& sourceImage) const
{
  Image result (sourceImage.getFormat(), sourceImage.getWidth(), sourceImage.getHeight(), true);
  Image full = createConvolvedImageFull (sourceImage);

  Graphics g (result);
  g.drawImageAt (full, -(m_radius - 1), -(m_radius - 1));

  return result;
}

Image RadialImageConvolutionKernel::createConvolvedImageFull (const Image& sourceImage) const
{
  // calc destination size based on kernel radius
  int dw = sourceImage.getWidth() + 2 * m_radius - 1;
  int dh = sourceImage.getHeight() + 2 * m_radius - 1;
  Image result (sourceImage.getFormat(), dw, dh, false);

  // calc size of edge-replicated source dimensions
  int sw = dw + 2 * m_radius - 1;
  int sh = dh + 2 * m_radius - 1;

  // temp buffer is big enough for the largest edge-replicated line
  HeapBlock<uint8> temp;
  temp.malloc (jmax (sw, sh));

  const Image::BitmapData srcData (sourceImage, false);
  const Image::BitmapData destData (result, 0, 0, dw, dh, true);

  int ci[4];
  int nc = srcData.pixelStride;
  bool alpha = false;

  switch (srcData.pixelFormat)
  {
  case Image::RGB:
    ci[0] = PixelRGB::indexR;
    ci[1] = PixelRGB::indexG;
    ci[2] = PixelRGB::indexB;
    nc = 3;
    alpha = false;
    break;

  case Image::ARGB:
    ci[0] = PixelARGB::indexR;
    ci[1] = PixelARGB::indexG;
    ci[2] = PixelARGB::indexB;
    ci[3] = PixelARGB::indexA;
    nc = 3;
    alpha = true;
    break;

  case Image::SingleChannel:
    ci[0] = 0;
    nc = 0;
    alpha = true;
  }

  // edge-replicate each row in source to temp, and convolve into dest
  for (int y = 0; y < srcData.height; ++y)
  {
    for (int c = 0; c < nc; ++c)
    {
      copy (srcData.width,
            temp,
            1,
            srcData.getLinePointer (y) + ci[c],
            srcData.pixelStride,
            2 * (m_radius - 1));

      convolve (destData.width,
                destData.getPixelPointer (0, y + m_radius - 1) + ci[c],
                destData.pixelStride,
                temp,
                1,
                m_radius,
                m_kernel);
    }

    if (alpha)
    {
      copy_alpha (srcData.width,
                  temp,
                  1,
                  srcData.getLinePointer (y) + ci[nc],
                  srcData.pixelStride,
                  2 * (m_radius - 1));

      convolve (destData.width,
                destData.getPixelPointer (0, y + m_radius - 1) + ci[nc],
                destData.pixelStride,
                temp,
                1,
                m_radius,
                m_kernel);
    }
  }

  // edge-replicate each intermediate column from dest to temp, and convolve into dest
  for (int x = 0; x < destData.width; ++x)
  {
    for (int c = 0; c < nc; ++c)
    {
      copy (srcData.height,
            temp,
            1,
            destData.getPixelPointer (x, m_radius - 1) + ci[c],
            destData.lineStride,
            2 * (m_radius - 1));

      convolve (destData.height,
                destData.getPixelPointer (x, 0) + ci[c],
                destData.lineStride,
                temp,
                1,
                m_radius,
                m_kernel);
    }

    if (alpha)
    {
      copy_alpha (srcData.height,
                  temp,
                  1,
                  destData.getPixelPointer (x, m_radius - 1) + ci[nc],
                  destData.lineStride,
                  2 * (m_radius - 1));

      convolve (destData.height,
                destData.getPixelPointer (x, 0) + ci[nc],
                destData.lineStride,
                temp,
                1,
                m_radius,
                m_kernel);
    }
  }

  return result;
}

void RadialImageConvolutionKernel::copy (int pixels,
                                         uint8* dest,
                                         int destSkip,
                                         const uint8* src,
                                         int srcSkip,
                                         int replicas)
{
  jassert (pixels > 0);

  for (int i = replicas; --i >= 0;)
  {
    *dest = *src;
    dest += destSkip;
  }

  for (int i = pixels; --i > 0;) // pixels-1 iterations
  {
    *dest = *src;
    src += srcSkip;
    dest += destSkip;
  }

  for (int i = replicas + 1; --i >= 0;) // extra iteration
  {
    *dest = *src;
    dest += destSkip;
  }
}

void RadialImageConvolutionKernel::copy_alpha (int pixels,
                                               uint8* dest,
                                               int destSkip,
                                               const uint8* src,
                                               int srcSkip,
                                               int replicas)
{
  jassert (pixels > 0);

  for (int i = replicas; --i >= 0;)
  {
    *dest = 0;
    dest += destSkip;
  }

  for (int i = pixels; --i >= 0;)
  {
    *dest = *src;
    src += srcSkip;
    dest += destSkip;
  }

  for (int i = replicas; --i >= 0;)
  {
    *dest = 0;
    dest += destSkip;
  }
}

void RadialImageConvolutionKernel::convolve (int pixels,
                                             uint8* dest,
                                             int destSkip,
                                             const uint8* src,
                                             int srcSkip,
                                             int radius,
                                             const float* kernel)
{
  for (int n = pixels; --n >= 0;)
  {
    const uint8* in = src;
    const float* k = kernel + radius;
    float tot = 0;

    for (int i = radius; --i >= 0;)
    {
      tot += *--k * *in;
      in += srcSkip;
    }

    for (int i = radius; --i > 0;)
    {
      tot += *++k * *in;
      in += srcSkip;
    }

    *dest = tot;//uint8 (jmin (0xff, roundToInt (tot)));
    dest += destSkip;

    src += srcSkip;
  }
}

Quote:
In addition to being circularly symmetric, the Gaussian blur can be applied to a two-dimensional image as two independent one-dimensional calculations, and so is termed linearly separable. That is, the effect of applying the two-dimensional matrix can also be achieved by applying a series of single-dimensional Gaussian matrices in the horizontal direction, then repeating the process in the vertical direction.
TheVinn
Offline
Last seen: 6 months 2 weeks ago
Joined: 29 Aug 2009 - 11:31
Re: Fast Image Convolution Using Radial Symmetry (with sourc

Judging from the lack of responses in a little over a year I think, the significance of this class is not well understood. Perhaps a diagram will help:

[attachment=0]RadialBlur.png[/attachment]

The starting image is of type ARGB. Three color channels (red, green, blue) and an associated alpha channel. Four channels in total. The code I posted detects that the image has an alpha channel, and applies the proper algorithm to blur not only the color channels but also the transparency information.

attachment: 
jules
Offline
Last seen: 32 min 3 sec ago
Joined: 29 Apr 2013 - 18:37
Re: Fast Image Convolution Using Radial Symmetry (with sourc

Never noticed this post first time, but thanks for bumping it - very interesting!

Anima
Offline
Last seen: 10 hours 34 min ago
Joined: 16 Dec 2010 - 14:32
Re: Fast Image Convolution Using Radial Symmetry (with sourc

Very cool

X-Ryl669
Offline
Last seen: 4 months 1 week ago
Joined: 24 Apr 2005 - 17:30
Re: Fast Image Convolution Using Radial Symmetry (with sourc

Yep, it's cool. And special upvote for not using that boost-o-matic creature.

X-Ryl669

TheVinn
Offline
Last seen: 6 months 2 weeks ago
Joined: 29 Aug 2009 - 11:31
Re: Fast Image Convolution Using Radial Symmetry (with sourc

*cough* DropShadow...

Everything discussed in this post is implemented here:

https://github.com/vinniefalco/VFLib/blob/master/modules/vf_gui/graphics/vf_RadialImageConvolutionKernel.cpp

It works on single channel images...