Fast statistical calculations of sub matrices for image processing

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:

the ii matrix is:

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:

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.