阅读opencv计算机视觉编程三（积分图像统计像素）

使用积分图像统计图像感兴趣区域的像素是一种高效的方法。它在程序中的应用非常广泛，
例如用于计算基于不同大小的滑动窗口。
本节将讲解积分图像背后的原理。这里的目标是说明如何只用三次算术运算，就能累加一个
矩形区域的像素
实现原理
为了理解积分图像的实现原理，我们先对它下一个定义：
取图像左上方的全部像素计算累加和，并用这个累加和替换图像中的每一个像素，用这种方式得
到的图像称为积分图像。计算积分图像时，只需对图像扫描一次。实际上，当前像素的积分值等
于上方像素的积分值加上当前行的累计值。因此积分图像就是一个包含像素累加和的新图像。为
防止溢出，积分图像的值通常采用int 类型（CV_32S）或float 类型（CV_32F）。例如在下图
中，积分图像的像素A 包含左上角区域，即双阴影线图案标识的区域的像素的累加和。

计算完积分图像后，只需要访问四个像素就可以得到任何矩形区域的像素累加和。这里解释
一下原因。再来看上面的图片，计算由A、B、C、D 四个像素表示区域的像素累加和，先读取D
的积分值，然后减去B 的像素值和C 的左手边区域的像素值。但是这样就把A 左上角的像素累
加和减了两次，因此需要重新加上A 的积分值。所以计算A、B、C、D 区域内的像素累加的正
式公式为：ABC + D。如果用cv::Mat 方法访问像素值，公式可转换成以下代码：
// 窗口的位置是(xo,yo)，尺寸是width×height

return (integralImage.at<cv::Vec<T,N>>(yo+height,xo+width)-
integralImage.at<cv::Vec<T,N>>(yo+height,xo)-
integralImage.at<cv::Vec<T,N>>(yo,xo+width)+
integralImage.at<cv::Vec<T,N>>(yo,xo));

不管感兴趣区域的尺寸有多大，使用这种方法计算的复杂度是恒定不变的。注意，为了简化，
这里使用了cv::Mat 类的at 方法，它访问像素值的效率并不是最高的。

积分图像适合用来执行多次像素累计值的统计。本段将通过介绍自适应阈值化的概念，说明
积分图像的使用方法。在需要快速计算多个窗口的直方图时，积分图像非常有用。本节也将对此
进行解释。

1. 自适应的阈值化
   通过对图像应用阈值来创建二值图像是从图像中提取有意义元素的好方法。假设有下面这个
   关于本书的图像。
   

为了分析图像中的文字，对该图像应用一个阈值，代码如下所示：
// 使用固定的阈值

cv::Mat binaryFixed;
cv::threshold(image,binaryFixed,70,255,cv::THRESH_BINARY);

得到如下结果。

实际上，不管选用什么阈值，图像都会丢失一部分文本，还有部分文本会消失在阴影下。要
解决这个问题，有一个办法就是采用局部阈值，即根据每个像素的邻域计算阈值。这种策略称为
自适应阈值化，将每个像素的值与邻域的平均值进行比较。如果某像素的值与它的局部平均值差
别很大，就会被当作异常值在阈值化过程中剔除。
因此自适应阈值化需要计算每个像素周围的局部平均值。这需要多次计算图像窗口的累计
值，可以通过积分图像提高计算效率。正因为如此，方法的第一步就是计算积分图像：
// 计算积分图像

cv::Mat iimage;
cv::integral(image,iimage,CV_32S);

现在就可以遍历全部像素，并计算方形邻域的平均值了。我们也可以使用IntegralImage
类来实现这个功能，但是这个类在访问像素时使用了效率很低的at 方法。根据第2 章学过的方法，
我们可以使用指针遍历图像以提高效率，循环代码如下所示：

int blockSize= 21; // 邻域的尺寸
int threshold=10; // 像素将与(mean-threshold)进行比较
// 逐行
int halfSize= blockSize/2;
for (int j=halfSize; j<nl-halfSize-1; j++) {
// 得到第j 行的地址
uchar* data= binary.ptr<uchar>(j);
int* idata1= iimage.ptr<int>(j-halfSize);
int* idata2= iimage.ptr<int>(j+halfSize+1);
// 一个线条的每个像素
for (int i=halfSize; i<nc-halfSize-1; i++) {
// 计算累加值
int sum= (idata2[i+halfSize+1]-data2[i-halfSize]-
idata1[i+halfSize+1]+idata1[i-halfSize])
/(blockSize*blockSize);
// 应用自适应阈值
if (data[i]<(sum-threshold))
data[i]= 0;
else
data[i]=255;
}
}

本例使用了21×21 的邻域。为计算每个平均值，我们需要访问界定正方形邻域的四个积分像
素：两个在标有idata1 的线条上，另两个在标有idata2 的线条上。将当前像素与计算得到的
平均值进行比较。为了确保被剔除的像素与局部平均值有明显的差距，这个平均值要减去阈值（这
里设为10）。由此得到下面的二值图像。

很明显，这比用固定阈值得到的结果好得多。自适应阈值化是一种常用的图像处理技术。
OpenCV 中也实现了这种方法：

cv::adaptiveThreshold(image, // 输入图像
binaryAdaptive, // 输出二值图像
255, // 输出的最大值
cv::ADAPTIVE_THRESH_MEAN_C, // 方法
cv::THRESH_BINARY, // 阈值类型
blockSize, // 块的大小
threshold); // 使用的阈值
调用这个函数得到的结果与使用积分图像的结果完全相同。另外，除了在阈值化中使用局部
平均值，本例中的函数还可以使用高斯（Gaussian）加权累计值（该方法的标志为ADAPTIVE_
THRESH_GAUSSIAN_C）。有意思的是，这种实现方式要比调用cv::adaptiveThreshold 稍微
快一些。
最后需要注意，我们也可以用OpenCV 的图像运算符来编写自适应阈值化过程。具体方法如
下所示：
cv::Mat filtered;
cv::Mat binaryFiltered;
// boxFilter 计算矩形区域内像素的平均值
cv::boxFilter(image,filtered,CV_8U,cv::Size(blockSize,blockSize));
// 检查像素是否大于(mean + threshold)
binaryFiltered= image>= (filtered-threshold);

2用直方图实现视觉追踪
通过前面几节的学习，我们知道可用直方图表示物体外观的全局特征。本节将搜寻一个所呈
现直方图与目标物体相似的图像区域，演示如何在图像中定位物体，以此说明积分图像的用途。
我们在4.6 节实现了这个功能，用的是直方图反向投影概念和通过均值偏移局部搜索的方法。这
次我们在整幅图像上显式地搜索具有类似直方图的区域，以此找到物体。
由0 和1 组成的二值图像生成积分图像是一种特殊情况，这时的积分累计值就是指定区域内
值为1 的像素总数。本节将利用这一现象计算灰度图像的直方图。
cv::integral 函数也可用于多通道图像。你可以充分利用这点，用积分图像计算图像子
区域的直方图。只需简单地把图像转换成由二值平面组成的多通道图像，每个平面关联直方图
的一个箱子，并显示哪些像素的值会进入该箱子。下面的函数将从一个灰度图像创建这样的多
图层图像：

// 转换成二值图层组成的多通道图像
// nPlanes 必须是2 的幂
void convertToBinaryPlanes(const cv::Mat& input,
cv::Mat& output, int nPlanes) {
// 需要屏蔽的位数
int n= 8-static_cast<int>(
log(static_cast<double>(nPlanes))/log(2.0));
// 用来消除最低有效位的掩码
uchar mask= 0xFF<<n;
// 创建二值图像的向量
std::vector<cv::Mat> planes;
// 消除最低有效位，箱子数减为nBins
cv::Mat reduced= input&mask;
// 计算每个二值图像平面
for (int i=0; i<nPlanes; i++) {
// 将每个等于i<<shift 的像素设为1
planes.push_back((reduced==(i<<n))&0x1);
}
// 创建多通道图像
cv::merge(planes,output);
你也可以把积分图像的计算过程封装进模板类中：
template <typename T, int N>
class IntegralImage {
cv::Mat integralImage;
public:
IntegralImage(cv::Mat image) {
// 计算积分图像（很耗时）
cv::integral(image,integralImage,
cv::DataType<T>::type);
}
// 通过访问四个像素，计算任何尺寸子区域的累计值
cv::Vec<T,N> operator()(int xo, int yo, int width, int height) {
// （xo,yo）处的窗口，尺寸为width×height
return (integralImage.at<cv::Vec<T,N>>(yo+height,xo+width)-
integralImage.at<cv::Vec<T,N>>(yo+height,xo)-
integralImage.at<cv::Vec<T,N>>(yo,xo+width)+
integralImage.at<cv::Vec<T,N>>(yo,xo));
}
};

用骑车女孩的图片举例

我们在前面的图像中识别出了骑车的女孩，现在要在后面的图像中找到她。首先计算原始图
像中女孩的直方图，这可通过Histogram1D 类实现。以下代码将生成16 个箱子的
直方图：

// 16 个箱子的直方图
Histogram1D h;
h.setNBins(16);
// 计算图像中ROI 的直方图
cv::Mat refHistogram= h.getHistogram(roi);
这个直方图将作为基准，在下面的图像中定位目标（即骑车的女孩）。
假设我们仅知道图像中女孩在水平方向移动。因为需要对不同的位置计算很多直方图，我们
先做准备工作，即计算积分图像。参见以下代码：
// 首先创建16 个平面的二值图像
cv::Mat planes;
convertToBinaryPlanes(secondIimage,planes,16);
// 然后计算积分图像
IntegralImage<float,16> intHistogram(planes);
执行搜索时，循环遍历可能出现目标的位置，并将它的直方图与基准直方图做比较，目的是
找到与直方图最相似的位置，参见以下代码：
double maxSimilarity=0.0;
int xbest, ybest;
// 遍历原始图像中女孩位置周围的水平长条
for (int y=110; y<120; y++) {
for (int x=0; x<secondImage.cols-width; x++) {
// 用积分图像计算16 个箱子的直方图
histogram= intHistogram(x,y,width,height);
// 计算与基准直方图的差距
double distance= cv::compareHist(refHistogram,
histogram,
CV_COMP_INTERSECT);
// 找到最相似直方图的位置
if (distance>maxSimilarity) {
xbest= x;
ybest= y;
maxSimilarity= distance;
}
}
}

// 在最准确的位置画矩形

cv::rectangle(secondImage, cv::Rect(xbest,ybest,width,height),0));

然后就可确定直方图最相似的位置，如下图所示。

白色矩形表示搜索的区域。计算区域内部所有窗口的直方图。这里的窗口尺寸是固定的，但
是更好的做法是也搜索稍小或稍大的窗口，以便应对缩放比例可能带来的变动。有一点需要注意，
为了降低计算复杂度，要减少直方图中要计算的箱子数量。本例减少到16 个箱子。因此，在这
个多平面图像中，平面0 包含一个二值图像，表示值从0 到15 的所有像素；平面1 表示值从16
到31 的全部像素，等等。
对物体的搜索过程包含了用预定范围的像素，计算指定尺寸的所有窗口的直方图的计算过
程，这意味着从积分图像对3200 个直方图进行了高效计算。IntegralImage 类返回的直方图
都存储在cv::Vec 对象中（因为用了at 方法）。然后用cv::compareHist 函数找到最相似的
直方图（和大多数OpenCV 函数一样，这个函数可以利用实用的通用参数类型cv::InputArray
获得cv::Mat 或cv::Vec）。