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 (Gaussian blur - Wikipedia) 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;
}
}