11,716,249 members (88,124 online)

# Fast statistical calculations of sub matrices for image processing

, 1 Oct 2005 46.4K 34
 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

## You may also be interested in...

 First Prev Next
 Why don't use managed DirectX 9... hakervytas29-Apr-06 2:07 hakervytas 29-Apr-06 2:07
 Small error in VarSubMatrix(), and efficiency mi5ke27-Jan-06 1:03 mi5ke 27-Jan-06 1:03
 Re: Small error in VarSubMatrix(), and efficiency ediazc27-Jan-06 5:35 ediazc 27-Jan-06 5:35
 one less array regnwald20-Dec-05 14:32 regnwald 20-Dec-05 14:32
 number plate recognition--- Need Help!! harshi_mm214-Dec-05 0:05 harshi_mm2 14-Dec-05 0:05
 A fast algorithm using slow constructs Jeffrey Sax13-Oct-05 16:33 Jeffrey Sax 13-Oct-05 16:33
 Re: A fast algorithm using slow constructs ediazc13-Oct-05 17:53 ediazc 13-Oct-05 17:53
 Re: A fast algorithm using slow constructs Jeffrey Sax13-Oct-05 18:48 Jeffrey Sax 13-Oct-05 18:48
 ediazc wrote:is this is a compiler flaw, or a restriction of CLI? I would say it's mostly a CLI issue. There are specialized IL (Intermediate Language) instructions for handling one-dimensional, zero-based arrays or 'vectors' (stelem, ldelem and ldelema to store, load, and load the address of an array element). The JIT compiler can convert these to near-optimal code. If bounds checking can be eliminated the resulting code's speed is of the same order as C code. In fact, the optimizations are so good that the 'standard' implementation of a nested loop (using a[x+width*y]) outperforms any attempt at manual optimization (for example, by keeping the index in a helper variable that is incremented rather than recalculating it each time). There are no such specialized instructions for multi-dimensional (MD) arrays. It has to go through generic methods that treat each element as an object. So every access causes the value to be boxed and unboxed. With the introduction of generics in .NET 2.0, it would be fairly straightforward to write a multi-dimensional array class that easily outperforms the built-in multidimensional arrays. For an excellent backgrounder on arrays, see Wesner Moise's ARRAYS Undocumented[^] article here on the Code Project. Jeffrey Everything should be as simple as possible, but not simpler.    -- Albert Einstein Numerical components for C# and VB.NET
 Re: A fast algorithm using slow constructs ediazc14-Oct-05 3:22 ediazc 14-Oct-05 3:22
 Last Visit: 31-Dec-99 18:00     Last Update: 2-Sep-15 23:36 Refresh 1