图像处理子矩阵的快速统计计算
4.47/5 (8投票s)
2005 年 10 月 1 日
3分钟阅读
62930
描述了 Franklin Crow 的面积表算法。
引言
当您参与图像处理或图形编程时，模式识别通常使用代表图像的矩阵完成。本文展示了一个最初由 Franklin Crow 在 SIGGRAPH '83 上描述的 算法。该算法用于高效地计算子矩阵的和。这些技术用于，例如，在图片中查找模式，例如用于 面部识别。
面积表
我们有一个矩阵 i，由 N 行 x M 列组成，通常这个矩阵表示为 i(x,y)，其中 0 <= x < N，并且 0 <= y < M。我们将通过以下递归构建一个面积表，或积分图像 ii(x,y)
s(x, y) = s(x, y-1) + i(x,y)
ii(x, y) = ii(x-1, y) + s(x,y)
其中 s(x,y) 是累积行和，s(x, -1) = 0，并且 ii(-1,y) = 0。ii(x,y) 是 i(x',y') 值在 x' < x, y' < y 上的总和。
给定以下矩阵 i
|
1 |
2 |
3 |
|
4 |
5 |
6 |
|
7 |
8 |
9 |
|
10 |
11 |
12 |
ii 矩阵是
|
1 |
3 |
6 |
|
5 |
12 |
21 |
|
12 |
27 |
45 |
|
22 |
48 |
78 |
现在，如果我们想计算子矩阵 (1,1)-(2,2) 的和 = 5+6+8+9 = 28。这个值等于 (45 + 1) - (12+6) = 28。
我们有以下模式
| A 1 | B 2 | |
| C 3 | D 4 | |
这里 1、2、3 和 4 是每个子矩阵 A、B、C 和 D 右下角的单元格。
如果我们选择 x、y、a 和 b，使得 (x+a)*(y+b) 是内部矩形的和，x*y 是矩形 A 的和，而 a*b 是内部矩形 D 的和，(x+a)*y 是矩形 A 和 B 的和，(y+b)*x 是矩形 A 和 C 的和。
我们写 |D| = a*b, |A+B+C+D| = (x+a)*(y+b), |A+B| = y*(x+a), |A+C| = x*(y+b)。
我们有 (x+a)*(y+b) = x*y+a(x+y)+b(x+y)+a*b，这意味着 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| 是单元格 4 的值。|A| 是单元格 1 的值，|A+B| 是单元格 2 的值，而 |A+C| 是单元格 3 的值。因此，D 的总和由 (4+1)-(2+3) 给出。
实现
我喜欢 C# 的一个原因是多维数组。您可以使用类似于 Pascal 的索引创建数组。例如，您可以这样创建一个 4x3 的矩阵
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;
鉴于此，我编写了一个简单的控制台程序来演示此算法，该程序生成一个具有随机值的矩阵，并计算给定子矩阵的和。您必须向程序提供 6 个参数，即矩阵的维度：N、M，以及子矩阵的“坐标”：左、上、右、下。
这是源代码
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();
}
}
}
扩展
您可以构建矩阵的平方值的积分图像 ii，这允许计算子矩阵的方差。
请记住，方差可以通过以下等式计算
Var = m2 - 1/N* SUM(x2)
其中 m 是平均值，x 是每个单元格的值，N 是子矩阵中单元格的总数。
您可以使用以下方法计算平均值
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));
}
并将方差计算为
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;
}
结论
该算法适用于任何类型的二维矩阵。当您必须计算子矩阵上的总和、其他统计值（如均值和方差）时，该算法保证了更少的操作。