# Fast statistical calculations of sub matrices for image processing

, 1 Oct 2005
 Rate this:
Describes the summed area table algorithm of Franklin Crow.

## Introduction

When you are involved in image processing or graphic programming, pattern recognition is usually done with matrices that represent the images. This article shows an algorithm described first by Franklin Crow on SIGGRAPH '83. The algorithm is used to do efficient sums of sub matrices. These techniques are used, for example, for finding patterns inside pictures for face recognition.

## The summed area table

We have a matrix i, of N rows x M columns, usually this matrix is denoted as i(x,y), where 0 <= x < N, and 0 <= y < M. We will build a summed area table, or integral image ii(x,y) by the following recurrence:

```s(x, y) = s(x, y-1) + i(x,y)
ii(x, y) = ii(x-1, y) + s(x,y)```

where s(x,y) is the cumulative row sum, s(x, -1) = 0, and ii(-1,y) = 0. ii(x,y) is the sum over x' < x, y' < y of i(x',y') values.

Given the following matrix i:

 1 2 3 4 5 6 7 8 9 10 11 12

the ii matrix is:

 1 3 6 5 12 21 12 27 45 22 48 78

Now, if we want to calculate the sum of the sub matrix (1,1)-(2,2) = 5+6+8+9 = 28. This value is equal to (45 + 1) - (12+6) = 28.

We have the following pattern:

 A 1 B 2 C 3 D 4

Here 1, 2, 3 and 4 are the cells at the right bottom corner on each sub matrix A, B, C and D.

If we choose x, y, a and b such as (x+a)*(y+b) is the sum of the internal rectangles, x*y is the sum of rectangle A, and a*b is the sum of the internal rectangle D, (x+a)*y is the sum of the rectangles A and B, (y+b)*x is the sum of the rectangles A and C.

We write |D| = a*b, |A+B+C+D| = (x+a)*(y+b), |A+B| = y*(x+a), |A+C| = x*(y+b).

We have (x+a)*(y+b) = x*y+a(x+y)+b(x+y)+a*b, this implies a*b = (x+a)*(y+b)-a(x+y)+b(x+y) = (x+a)*(y+b)-y(x+a)-x(y+b)+xy = |A+B+C+D|-|A+B|-|A+C| + |A| = |D|.

|A+B+C+D| is the value of the cell 4. |A| is the value on cell 1, |A+B| is the value on cell 2 and |A+C| is the value on cell 3. So the sum of D is given by (4+1)-(2+3).

## Implementation

One of the things I like about C# is the multi-dimensional array. You can have an array indexed like in Pascal. For example, you can create a matrix of 4x3 this way:

```double[,] matrix = new double[2,3];
for (int y = 0; y < 3; y++)
for (int x = 0; x < 2; x++)
matrix[x,y] = x*y;```

Given this, I wrote a simple console program to demonstrate this algorithm, the program generates a matrix with random values and calculates the sum of a given sub matrix. You must provide 6 parameters to the program, the dimensions of matrix: N, M, and the "coordinates" of the sub matrix: left, top, right, bottom.

Here's the source code:

```using System;
using System.IO;
public class FastMatrixSum
{
public static void Main(string[] args)
{
if (args.Length != 6) {
Console.WriteLine("usage: matrixsum N M " +
"left top right bottom");
Console.WriteLine("This program build a random " +
"matrix of N rows with M columns,");
Console.WriteLine("and calculate the sum of the " +
"sub matrix (left,top)-(right,bottom)");
Console.WriteLine("This are the conditions for the parameters");
Console.WriteLine("N > 0 and M > 0");
Console.WriteLine("(0 <= left < M) and (0 <= right < M)");
Console.WriteLine("(0 <= top < N) and (0 <= bottom < M)");
Console.WriteLine("left <= right and top <= bottom");

}

int n = Convert.ToInt32(args[0]);
int m = Convert.ToInt32(args[1]);
int left = Convert.ToInt32(args[2]);
int top  = Convert.ToInt32(args[3]);
int right = Convert.ToInt32(args[4]);
int bottom = Convert.ToInt32(args[5]);
double[,] matrix = FillMatrix(n,m);
double[,] sum = SumMatrix(matrix);

Console.WriteLine("Matrix");
PrintMatrix(matrix);
Console.WriteLine("Sum(Matrix)");
PrintMatrix(sum);

Console.WriteLine();
Console.WriteLine("The sum fo the sub matrix " +
"({0},{1})-({2},{3}) is = {4}",
left,top,right,bottom,SumSubMatrix(sum,
left,top,right,bottom));
}
static double SumSubMatrix(double[,] sum, int left,
int top, int right, int bottom)
{
double v1, v2, v3, v4;
v1 = (left == 0 || top == 0) ? 0 : sum[left-1, top-1];
v2 = (top == 0) ? 0 : sum[right,top-1];
v3 = (left == 0) ? 0 : sum[left-1,bottom];
v4 = sum[right,bottom];
return (v4+v1)-(v2+v3);

}
static double[,] FillMatrix(int n, int m)
{
Random rnd = new Random();
double[,] result = new double[m,n];
for (int i = 0; i < n; i++)
for (int j = 0; j < m; j++)
result[j,i] = rnd.NextDouble() * 100.0;
return result;
}

static double[,] SumMatrix(double[,] matrix)
{
int n = matrix.GetLength(0);
int m = matrix.GetLength(1);
double[,] s = new double[n,m];
double[,] ii = new double[n,m];

s[0,0] = matrix[0,0];
ii[0,0] = s[0,0];
for (int x = 1; x < n; x++)
{
s[x,0] = matrix[x,0];
ii[x,0] = ii[x-1,0] + s[x,0];
}
for (int y = 1; y < m; y++) {
ii[0,y] = s[0,y] = s[0,y-1] + matrix[0, y];

for (int x = 1; x < n; x++) {
s[x,y] = s[x,y-1] + matrix[x,y];
ii[x,y] = ii[x-1,y] + s[x,y];
}
}
return ii;
}

static void PrintMatrix(double[,] matrix)
{
int n = matrix.GetLength(0);
int m = matrix.GetLength(1);
for (int y = 0; y < m; y++)
{
for (int x = 0; x < n; x++)
Console.Write(" {0,8:00.00} ", matrix[x,y]);
Console.WriteLine();
}
}
}```

## Extensions

You can build an integration image ii of the square values for the matrix, this allows to calculate the variance of a sub matrix.

Remember that the variance can be calculated by the following equation:

`Var = m<SUP>2</SUP> - 1/N* SUM(x<SUP>2</SUP>)`

where m is the mean and x is the value at each cell, and N is the total number of cells on the sub matrix.

You can calculate mean with the following method:

```static double MeanSubMatrix(double[,] sum, int left,
int top, int right, int bottom)
{
double v1, v2, v3, v4;
v1 = (left == 0 || top == 0) ? 0 : sum[left-1, top-1];
v2 = (top == 0) ? 0 : sum[right,top-1];
v3 = (left == 0) ? 0 : sum[left-1,bottom];
v4 = sum[right,bottom];
return ((v4+v1)-(v2+v3)) / ((bottom-top)+(right-left));
}```

and the variance as:

``` static double[,] SumMatrixOfSquare(double[,] matrix)
{
int n = matrix.GetLength(0);
int m = matrix.GetLength(1);
double[,] s = new double[n,m];
double[,] ii = new double[n,m];

s[0,0] = matrix[0,0];
ii[0,0] = s[0,0];
for (int x = 1; x < n; x++)
{
s[x,0] = matrix[x,0]*matrix[x,0]; // x*x
ii[x,0] = ii[x-1,0] + s[x,0];
}
for (int y = 1; y < m; y++) {
ii[0,y] = s[0,y] = s[0,y-1] + matrix[0, y]*matrix[0,y];

for (int x = 1; x < n; x++) {
s[x,y] = s[x,y-1] + matrix[x,y]*matrix[x,y];
ii[x,y] = ii[x-1,y] + s[x,y];
}
}
return ii;
}
static double VarSubMatrix(double[,] sum, double[,] sumxx,
int left, int top, int right, int bottom)
{
double N = (bottom - top)+(right - left);
double   m = MeanSubMatrix(sum, left, top, right, bottom);
double v1, v2, v3, v4;
v1 = (left == 0 || top == 0) ? 0 : sumxx[left-1, top-1];
v2 = (top == 0) ? 0 : sumxx[right,top-1];
v3 = (left == 0) ? 0 : sumxx[left-1,bottom];
v4 = sumxx[right,bottom];
double sumS2 = (v4+v1)-(v2+v3);
return m*m - sumS2/N;
}```

## Conclusion

The algorithm is useful for any kind of two dimensional matrix. This algorithm guarantees fewer operations when you have to calculate sums, other statistics values, like mean, and variance on sub matrices.

A list of licenses authors might use can be found here

## Share

Web Developer
Chile
Eduardo Diaz
personal blog

 First Prev Next
 Why don't use managed DirectX 9... hakervytas 29-Apr-06 3:07
 Small error in VarSubMatrix(), and efficiency mi5ke 27-Jan-06 2:03
 Re: Small error in VarSubMatrix(), and efficiency ediazc 27-Jan-06 6:35
 one less array regnwald 20-Dec-05 15:32
 number plate recognition--- Need Help!! harshi_mm2 14-Dec-05 1:05
 A fast algorithm using slow constructs Jeffrey Sax 13-Oct-05 17:33
 Re: A fast algorithm using slow constructs ediazc 13-Oct-05 18:53
 Yes, you are right, I tried the linearized array implementation and performance impact is big. I will add a section on the article about this issue, in fact, I pretend to test an unmanaged version for my particular project.   Is sad that multidimensional array were not optimized, because syntax is more clear, like in Delphi (Object Pascal).   Sure, you are an expert on optimization, is this is a compiler flaw, or a restriction of CLI?   Eduardo Diaz   site | english blog | spanish blog
 Re: A fast algorithm using slow constructs Jeffrey Sax 13-Oct-05 19:48
 Re: A fast algorithm using slow constructs ediazc 14-Oct-05 4:22
 Last Visit: 31-Dec-99 19:00     Last Update: 25-Dec-14 8:50 Refresh 1