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

Share

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 3:07 
GeneralSmall error in VarSubMatrix(), and efficiency Pinmembermi5ke27-Jan-06 2:03 
GeneralRe: Small error in VarSubMatrix(), and efficiency Pinmemberediazc27-Jan-06 6:35 
Generalone less array Pinmemberregnwald20-Dec-05 15:32 
Generalnumber plate recognition--- Need Help!! Pinmemberharshi_mm214-Dec-05 1:05 
GeneralA fast algorithm using slow constructs PinmemberJeffrey Sax13-Oct-05 17:33 
GeneralRe: A fast algorithm using slow constructs Pinmemberediazc13-Oct-05 18:53 
GeneralRe: A fast algorithm using slow constructs PinmemberJeffrey Sax13-Oct-05 19:48 
GeneralRe: A fast algorithm using slow constructs Pinmemberediazc14-Oct-05 4:22 
Very insightful, thanks
Smile | :)
 
Eduardo Diaz
 
site | english blog | spanish blog

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 | Terms of Use | Mobile
Web04 | 2.8.141223.1 | Last Updated 1 Oct 2005
Article Copyright 2005 by ediazc
Everything else Copyright © CodeProject, 1999-2014
Layout: fixed | fluid