Click here to Skip to main content
Click here to Skip to main content

Fast statistical calculations of sub matrices for image processing

, 1 Oct 2005
Rate this:
Please Sign up or sign in to vote.
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.

License

This article has no explicit license attached to it but may contain usage terms in the article text or the download files themselves. If in doubt please contact the author via the discussion board below.

A list of licenses authors might use can be found here

About the Author

ediazc
Web Developer
Chile Chile
Eduardo Diaz
personal blog

Comments and Discussions

 
JokeWhy don't use managed DirectX 9... Pinmemberhakervytas29-Apr-06 2:07 
GeneralSmall error in VarSubMatrix(), and efficiency Pinmembermi5ke27-Jan-06 1:03 
Hi,
 
In your VarSubMatrix code, I noticed a small error:
 
   double N = (bottom - top)+(right - left);
 
should be:
 
   double N = (bottom - top)*(right - left);
 
since N is the number of elements in the region.
 
By the way, another optimisation often used is to linearise the ii array (as others have mentioned) and then to simply autoincrement an index through the array, avoiding any multiplication:
 
   int p = 0;
   for(y = ....
         for(x = ....
            ii[p] = ....
            p++;
 
Lastly, a common trick for linearising two-dimensional arrays is to use a width which is rounded up to the next power of two, referred to as the 'pitch' of the data (rather than the 'width'). Then, the conventional linearised array access:
 
   myArray[x + y * width]
 
can make use of the fact that width is known to be a power of two:
 
   myArray[x + y << pitchShift]
 
To initialise the pitch and pitchShift, just use:
 
   int pitchShift = 0;
   int pitch = 1;
   while(pitch < width) {
         pitchShift++;
         pitch <<= 1;
   }
 
Hope that helps.
 
Mi5ke
GeneralRe: Small error in VarSubMatrix(), and efficiency Pinmemberediazc27-Jan-06 5:35 
Generalone less array Pinmemberregnwald20-Dec-05 14:32 
Generalnumber plate recognition--- Need Help!! Pinmemberharshi_mm214-Dec-05 0:05 
GeneralA fast algorithm using slow constructs PinmemberJeffrey Sax13-Oct-05 16:33 
GeneralRe: A fast algorithm using slow constructs Pinmemberediazc13-Oct-05 17:53 
GeneralRe: A fast algorithm using slow constructs PinmemberJeffrey Sax13-Oct-05 18:48 
GeneralRe: A fast algorithm using slow constructs Pinmemberediazc14-Oct-05 3:22 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Mobile
Web01 | 2.8.140718.1 | Last Updated 1 Oct 2005
Article Copyright 2005 by ediazc
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid