## Jul 8, 2015

### Integral image and image blur

In this post I will show you an image processing technique that provides a great boost to algorithms that requires the summation of a large number of pixels, this technique is called integral image, also known as Summed Area Table (SAT). Additionally, I will show you an algorithm used to blur an image that take advantage of the integral image.
In the previous post when I talked about image convolution, I had shown a convolution kernel that produces the effect of blurring an image.

// Blur
int kw = 7;
int kh = 7;
qreal kernel[] = {1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49,
1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49,
1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49,
1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49,
1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49,
1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49,
1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49, 1. / 49};


Bigger the kernel, more noticeable is the blurring effect. The problem is that a convolution is very slow and expensive image processing technique, bigger the kernel, more time it takes to process a single image frame. The integral image allow to process a frame in constant time regardless of kernel size.
As said before, an image can be treated as a function of two variables, f(x, y), so the integral of f(x, y) is:

$$I(x,y)=\iint f(x,y)dxdy$$

and it's discrete form:

$$I(x, y)=\sum_{j=0}^{y} \sum_{i=0}^{x} f(i,j)$$

For example, if the pixels values of the image are:

1 2 3 4 5
6 7 8 9 10
11 12 13 14 15
16 17 18 19 20
21 22 23 24 25

Then, the resulting integral image is:

1 3 6 10 15
7 16 27 40 55
18 39 63 90 120
34 72 114 160 210
55 115 180 250 325

The formula for calculating sum of all pixels contained in the rectangle defined by (x, y, width, height) is:

\begin{align} sum & = integral(x-1,y-1) \\ & - integral(x+width-1,y-1) \\ & - integral(x-1,y+height-1) \\ & + integral(x+width-1,y+height-1) \end{align}

For example, if I want to calculate the sum of all pixels in the rectangle (1, 2, 4, 2), that is:

sum = 12 + 13 + 14 + 15
+ 17 + 18 + 19 + 20 = 128


using the above formula:

p0 = integral(x - 1, y - 1) = integral(1 - 1, 2 - 1) = integral(0, 1) = 7
p1 = integral(x + width - 1, y - 1) = integral(1 + 4 - 1, 2 - 1) = integral(4, 1) = 55
p2 = integral(x - 1, y + height - 1) = integral(1 - 1, 2 + 2 - 1) = integral(0, 3) = 34
p3 = integral(x + width - 1, y + height - 1) = integral(1 + 4 - 1, 2 + 2 - 1) = integral(4, 3) = 210

sum = p0 - p1 - p2 + p3 = 7 - 55 - 34 + 210 = 128


Now, what happens if I want to calculate the sum of the rectangle (0, 0, 5, 5) or any rectangle where it's coordinates x or y are zero? It will return a pixel that doesn't exists, but we can safely assume is equal to zero, so it will not affect the result of the formula. We can fix this problem padding the top and the left side of the SAT, making the values equal to zero.

0 0 0 0 0 0
0 1 3 6 10 15
0 7 16 27 40 55
0 18 39 63 90 120
0 34 72 114 160 210
0 55 115 180 250 325

In that case, the new formula is:

\begin{align} sum & = integral(x,y) \\ & - integral(x+width,y) \\ & - integral(x,y+height) \\ & + integral(x+width,y+height) \end{align}

And the code for the integral image is this:

#include <QImage>

class Pixel
{
public:
explicit Pixel():
r(0), g(0), b(0)
{
}
Pixel(quint32 r, quint32 g, quint32 b):
r(r), g(g), b(b)
{
}

Pixel operator +(const Pixel &other)
{
return Pixel(this->r + other.r,
this->g + other.g,
this->b + other.b);
}

Pixel &operator +=(QRgb pixel)
{
this->r += qRed(pixel);
this->g += qGreen(pixel);
this->b += qBlue(pixel);

return *this;
}

quint32 r;
quint32 g;
quint32 b;
};

void integralImage(const QImage &image, QVector<Pixel> &integral)
{
int oWidth = image.width() + 1;
int oHeight = image.height() + 1;
integral.resize(oWidth * oHeight);

for (int y = 1; y < oHeight; y++) {
const QRgb *line = (const QRgb *) image.constScanLine(y - 1);

// Reset current line summation.
Pixel sum;

for (int x = 1; x < oWidth; x++) {
QRgb pixel = line[x - 1];

// Accumulate pixels in current line.
sum += pixel;

// Offset to the current line.
int offset = x + y * oWidth;

// Offset to the previous line.
// equivalent to x + (y - 1) * oWidth;
int offsetPrevious = offset - oWidth;

// Accumulate current line and previous line.
integral[offset] = sum + integral[offsetPrevious];
}
}
}


The blur effect can be defined as the average of all pixels surrounding a pixel I(x, y) in a radius r:

$$Blur(x, y)=\frac{1}{(2r+1)^2}\sum_{j=-r}^{r}\sum_{i=-r}^{r}I(x+i, y+j)$$

And the final code for the image blur:

int main(int argc, char *argv[])
{
QCoreApplication a(argc, argv);
Q_UNUSED(a)

QImage inImage("lena.png");
inImage = inImage.convertToFormat(QImage::Format_RGB32);
QImage outImage(inImage.size(), inImage.format());

int oWidth = inImage.width() + 1;
QVector<Pixel> integral;
integralImage(inImage, integral);

// Here we configure how much strong the blur will be.

for (int y = 0; y < inImage.height(); y++) {
const QRgb *iLine = (const QRgb *) inImage.constScanLine(y);
QRgb *oLine = (QRgb *) outImage.scanLine(y);
int yp = qMax(y - radius, 0);
int kh = qMin(y + radius, inImage.height() - 1) - yp + 1;

for (int x = 0; x < inImage.width(); x++) {
int xp = qMax(x - radius, 0);
int kw = qMin(x + radius, inImage.width() - 1) - xp + 1;

Pixel *p0 = integral.data() + xp + yp * oWidth;
Pixel *p1 = p0 + kw;
Pixel *p2 = p0 + kh * oWidth;
Pixel *p3 = p2 + kw;

quint32 sumR = p0->r - p1->r - p2->r + p3->r;
quint32 sumG = p0->g - p1->g - p2->g + p3->g;
quint32 sumB = p0->b - p1->b - p2->b + p3->b;

qreal ks = kw * kh;

quint8 r = qRound(sumR / ks);
quint8 g = qRound(sumG / ks);
quint8 b = qRound(sumB / ks);

oLine[x] = qRgba(r, g, b, qAlpha(iLine[x]));
}
}

outImage.save("out.png");

return EXIT_SUCCESS;
}


And this is the result for r=50: