Fast Image Convolution Using Radial Symmetry (with source)


#1

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

#2

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.


#3

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


#4

Very cool


#5

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


#6

cough DropShadow…

Everything discussed in this post is implemented here:

It works on single channel images…