OpenCV实战（10）——积分图像详解
0. 前言
我们知道直方图是通过遍历图像的所有像素并累积每个强度值在该图像中出现的频率来计算的。有时，我们只对计算图像某些区域的直方图感兴趣，在许多计算机视觉算法中，累积图像子区域内的像素之和是一项常见任务。现在，假设必须在图像内的多个感兴趣区域上计算直方图，此时计算代价将变得非常高昂，在这种情况下，我们可以使用积分图像 (integral image) 极大地提高对图像子区域的像素计数效率。积分图像作为在图像 (Region of Interest, ROI) 中求和像素的有效方法，已经广泛用于多尺度滑动窗口计算等应用中。
1. 积分图像计算
本节，我们将解释积分图像的基本原理，展示如何仅使用三个算术运算在矩形区域上对像素求和，并通过两个示例介绍如何使用积分图像。
(1) 当我们需要对多个图像区域的像素求和时，积分图像很有用。通常，我们使用以下代码获得 ROI 上所有像素的总和：
cv::Mat image = cv::imread("example.png", 0);
// 定义 ROI
int xo = 97, yo = 112;
int width = 25, height = 30;
cv::Mat roi(image, cv::Rect(xo, yo, width, height));
cv::Scalar sum = cv::sum(roi);
(2) cv::sum 函数遍历区域的所有像素并累加总和，使用积分图像，可以仅使用三个加法运算来实现。首先，需要计算积分图像：
// 计算积分图像
cv::Mat iimage;
cv::integral(image, iimage, CV_32S);
(3) 在计算出的积分图像上使用简单的算术表达式可以获得与 cv::sum 相同的结果：
int sumInt = integralImage.at<int>(yo+height, xo+width)
-integralImage.at<int>(yo+height, xo)
-integralImage.at<int>(yo, xo+width)
+integralImage.at<int>(yo, xo);
这两种方法可以得到相同的结果。但是，计算积分图像的成本很高，因为必须遍历所有图像像素。但一旦完成初始计算，无论该区域的大小如何，都只需使用区域的 4 个顶点像素即可获得 ROI 的总和；当必须在不同大小的多个区域上计算多个像素总和时，积分图像将变得非常有效。
通过对积分图像的简要演示，我们了解了积分图像的概念，即如何使用它们来高效地计算矩形区域内的像素总和。在积分图像中，每个像素值都是通过该像素限定的左上象限内的所有像素的总和获得的。积分图像可以通过扫描图像一次来计算，因为当前像素的积分值由上一行(同一列)像素的积分值加上当前行的累积和的值给出。因此，积分图像是包含像素和的新图像，为避免溢出，此图像通常是 int 值 (CV_32S) 或浮点值 (CV_32F) 图像。例如，在下图中，该积分图像中的像素 A 包含左上角区域内包含的像素的总和：
一旦计算了积分图像，就可以通过访问 4 个顶点像素获得矩形区域的像素值之和。在上图中，我们可以看到，通过读取像素 D 处的积分值，减去像素 A、B、C 和 D 所划定的区域内的像素之和，可以得到 B 上方和 C 左侧的像素。但是，这样相当于多次减去了位于 A 左上角的像素之和，因此必须重新添加 A 处的积分和，综上，矩形 ABCD 内的像素总和可以根据 D-B-C+A 得出。使用 cv::Mat 方法访问像素值，则可以转换为以下代码：
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)
);
因此，无论 ROI 的大小如何，使用积分图像计算区域像素和的复杂度都是恒定的。当我们需要求取多个像素区域之和时，可以使用积分图像，积分图像可用于高效的计算多个矩形区域的直方图。接下来，我们通过引入自适应阈值的概念来说明如何使用积分图像。
2. 自适应阈值
在图像上应用阈值算法创建二值图像是提取图像中目标元素的好方法。假设我们有以下文档图像：
2.1 固定阈值的缺陷
如果我们需要分析此图像中的文本，可以对该图像应用阈值算法：
// 使用固定阈值
cv::Mat binaryFixed;
cv::threshold(image, binaryFixed, 70, 255, cv::THRESH_BINARY);
应用阈值算法后可以获得以下结果：
事实上，无论使用任何阈值，都会丢失图像中的某些部分，具体到以上文档图像，文本会被强光所覆盖。为了克服这个问题，一种解决方案是使用从每个像素的邻域计算出的局部阈值，这种策略称为自适应阈值 (Adaptive Thresholding)，将每个像素与相邻像素的平均值进行比较，明显不同于其局部平均值的像素将被视为异常值，并被阈值处理截断。
2.2 使用自适应阈值
因此，自适应阈值需要计算每个像素周围的局部均值，这需要可以通过积分图像计算多个图像窗口之和。因此，第一步是计算以下积分图像：
// 计算积分图像
cv::Mat iimage;
cv::integral(image, iimage, CV_32S);
遍历所有像素并计算像素邻域的均值，我们可以使用 IntegralImage 类实现，但这个类使用低效的 at 方法进行像素访问，我们可以通过使用指针循环图像来提高效率：
// 使用自适应阈值
int blockSize = 21;
int threshold = 10; // 像素将与(平均阈值)进行比较
int halfSize = blockSize / 2;
for (int j=halfSize; j<nl-halfSize-1; 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] - idata2[i-halfSize] -
idata1[i+halfSize+1] + idata1[i-halfSize]) /
(blockSize*blockSize);
// 应用自适应阈值
if (data[i] < (sum-threshold)) data[i] = 0;
else data[i] = 255;
}
}
在以上代码中，使用大小为 21 x 21 的邻域。为了计算均值，我们需要访问界定正方形邻域的四个积分像素——两个位于 idata1 线上，两个位于 idata2 线上。当前像素与由此计算的平均值进行比较，为了确保被截断的像素与其局部平均值明显不同，我们从平均值中减去阈值(此处设置为 10)，可以得到以下二值图像：
2.3 其它自适应阈值计算方法
显然，这比使用固定阈值得到的结果要好得多。自适应阈值是一种常见的图像处理技术，因此，在 OpenCV 中也提供了自适应阈值函数的实现：
cv::adaptiveThreshold(
image, // 输入图像
binaryAdaptive, // 输出二值图像
255, // 最大值
cv::ADAPTIVE_THRESH_MEAN_C, // 自适应均值
cv::THRESH_BINARY, // 阈值类型
blockSize, // 块大小
threshold // 阈值
);
这个函数调用产生的结果与我们使用积分图像获得的结果完全相同。此外，该函数允许我们使用高斯加权和 (ADAPTIVE_THRESH_GAUSSIAN_C)，而不是使用局部均值进行阈值处理。
除此之外，我们还可以使用 OpenCV 图像算子编写自适应阈值处理程序：
// 基于图像算子的自适应阈值
time= cv::getTickCount();
cv::Mat filtered;
cv::Mat binaryFiltered;
// 使用滤波器计算区域上像素的平均值
cv::boxFilter(image, filtered, CV_8U, cv::Size(blockSize, blockSize));
// 检查像素值是否大于(均值+阈值)
binaryFiltered = image >= (filtered - threshold);
2.4 完整代码
头文件 integral.h 代码如下所示：
#if !defined IINTEGRAL
#define IINTEGRAL
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <vector>
template <typename T, int N>
class IntegralImage {
cv::Mat integralImage;
public:
IntegralImage(cv::Mat image) {
// 积分图像
cv::integral(image, integralImage, cv::DataType<T>::type);
}
// 访问 4 个像素计算任意大小的子区域的像素总和
cv::Vec<T, N> operator()(int xo, int yo, int width, int 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::Vec<T, N> operator()(int x, int y, int radius) {
return (
integralImage.at<cv::Vec<T, N> >(y+radius+1, x+radius+1)
-integralImage.at<cv::Vec<T, N> >(y+radius+1, x-radius)
-integralImage.at<cv::Vec<T, N> >(y-radius, x+radius+1)
+integralImage.at<cv::Vec<T, N> >(y-radius, x-radius)
);
}
};
// 转换为由二进制平面构成的多通道图像
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;
// 创建16个二进制图像的矢量
std::vector<cv::Mat> planes;
// 通过消除最低有效位将 bins 数量减少到nBins
cv::Mat reduced = input & mask;
// 计算每个二进制图像平面
for (int i=0; i<nPlanes; i++) {
planes.push_back((reduced==(i<<n))&0x1);
}
// 创建多通道图像
cv::merge(planes, output);
}
#endif
主函数文件 integral.cpp 代码如下所示：
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include "integral.h"
int main() {
cv::Mat image = cv::imread("test.jpeg", 0);
if (!image.data) return 0;
// 旋转图像
cv::transpose(image, image);
cv::flip(image, image, 0);
cv::namedWindow("Original Image");
cv::imshow("Original Image", image);
// 使用固定阈值
cv::Mat binaryFixed;
cv::Mat binaryAdaptive;
cv::threshold(image, binaryFixed, 70, 255, cv::THRESH_BINARY);
// 使用自适应阈值
int blockSize = 21;
int threshold = 10; // 像素将与(平均阈值)进行比较
int64 time;
time= cv::getTickCount();
cv::adaptiveThreshold(
image, // 输入图像
binaryAdaptive, // 输出二进制图像
255, // 最大值
cv::ADAPTIVE_THRESH_MEAN_C, // 自适应均值
cv::THRESH_BINARY, // 阈值类型
blockSize, // 块大小
threshold // 阈值
);
time = cv::getTickCount() - time;
std::cout << "time (adaptiveThreshold) = " << time << std::endl;
// 计算积分图像
IntegralImage<int, 1> integral(image);
std::cout << "sum = " << integral(18, 45, 30, 50) << std::endl;
cv::Mat test(image, cv::Rect(18, 45, 30, 50));
cv::Scalar t = cv::sum(test);
std::cout << "sum test = " << t[0] << std::endl;
cv::namedWindow("Fixed Threshold");
cv::imshow("Fixed Threshold",binaryFixed);
cv::namedWindow("Adaptive Threshold");
cv::imshow("Adaptive Threshold",binaryAdaptive);
cv::Mat binary = image.clone();
time = cv::getTickCount();
int nl = binary.rows;
int nc = binary.cols;
// 计算积分图像
cv::Mat iimage;
cv::integral(image, iimage, CV_32S);
int halfSize = blockSize / 2;
for (int j=halfSize; j<nl-halfSize-1; 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] - idata2[i-halfSize] -
idata1[i+halfSize+1] + idata1[i-halfSize]) /
(blockSize*blockSize);
// 应用自适应阈值
if (data[i] < (sum-threshold)) data[i] = 0;
else data[i] = 255;
}
}
// 添加白板
for (int j=0; j<halfSize; j++) {
uchar* data = binary.ptr<uchar>(j);
for (int i=0; i<binary.cols; i++) {
data[i] = 255;
}
}
for (int j=binary.rows-halfSize-1; j<binary.rows; j++) {
uchar* data = binary.ptr<uchar>(j);
for (int i=0; i<binary.cols; i++) {
data[i] = 255;
}
}
for (int j=halfSize; j<nl-halfSize-1; j++) {
uchar* data = binary.ptr<uchar>(j);
for (int i=0; i<halfSize; i++) {
data[i] = 255;
}
for (int i = binary.cols-halfSize-1; i<binary.cols; i++) {
data[i] = 255;
}
}
time= cv::getTickCount()-time;
std::cout << "time integral= " << time << std::endl;
cv::namedWindow("Adaptive Threshold (integral)");
cv::imshow("Adaptive Threshold (integral)",binary);
// 基于图像算子的自适应阈值
time= cv::getTickCount();
cv::Mat filtered;
cv::Mat binaryFiltered;
// 使用滤波器计算区域上像素的平均值
cv::boxFilter(image, filtered, CV_8U, cv::Size(blockSize, blockSize));
// 检查像素值是否大于(均值+阈值)
binaryFiltered = image >= (filtered - threshold);
time= cv::getTickCount()-time;
std::cout << "time filtered= " << time << std::endl;
cv::namedWindow("Adaptive Threshold (filtered)");
cv::imshow("Adaptive Threshold (filtered)",binaryFiltered);
cv::waitKey();
}
3. 使用直方图进行视觉跟踪
我们已经知道直方图构成了对象外观的可靠全局表示。在本节中，我们将介绍如何通过搜索与目标对象具有相似的直方图的图像区域来定位图像中的对象。我们已经学习了使用均值偏移算法通过使用直方图反投影和通过均值偏移进行局部搜索来查找目标对象。在本节中，我们将通过在整个图像上对搜索具有相似直方图的区域来查找目标对象。
3.1 查找目标对象
由 0 和 1 值组成的二值图像中，积分图像中的值本质上等于指定区域内值为 1 的像素数。本节中，我们将利用这一事实来计算灰度图像的直方图。
cv::integral 函数也适用于多通道图像，可以使用积分图像计算图像子区域的直方图，只需要将图像转换为由二元平面组成的多通道图像；每个平面中都与直方图中一个 bin 相关联，并显示哪些像素的值属于该 bin。使用灰度图像创建多平面图像如下：
// 转换为由二进制平面构成的多通道图像
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;
// 创建16个二值图像的矢量
std::vector<cv::Mat> planes;
// 通过消除最低有效位将 bins 数量减少到nBins
cv::Mat reduced = input & mask;
// 计算每个二值图像平面
for (int i=0; i<nPlanes; i++) {
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);
}
// 访问 4 个像素计算任意大小的子区域的像素总和
cv::Vec<T, N> operator()(int xo, int yo, int width, int 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 个 bin 的直方图：
// 直方图
Histogram1D h;
h.setNBins(16);
cv::Mat refHistogram = h.getHistogram(roi);
以上直方图将作为参考表示用于在后续图像中定位目标对象。由于我们将在不同位置计算多个直方图，因此我们需要首先计算积分图像：
// 创建具有16个平面的二值图像
cv::Mat planes;
convertToBinaryPlanes(image, planes, 16);
// 计算积分图像
IntegralImage<float, 16> intHisto(planes);
为了执行搜索，我们遍历一系列可能的位置并将当前直方图与参考直方图进行比较，目标是找到与参考直方图最相似的直方图位置：
double maxSimilarity = 0.0;
int xbest, ybest;
// 循环查找
for (int y=160; y<200; y++) {
for (int x=0; x<secondImage.cols-width; x++) {
// 计算直方图
histogram = intHistogram(x, y, width, height);
double distance = cv::compareHist(refHistogram, histogram, cv::HISTCMP_INTERSECT);
// 查找最相似直方图
if (distance>maxSimilarity) {
xbest = x;
ybest = y;
maxSimilarity = distance;
}
}
}
// 在目标对象周围绘制矩形框
cv::rectangle(image, cv::Rect(xo, yo, width, height), 0);
识别与参考直方图最相似直方图的位置：
白色矩形表示搜索区域(假设目标对象只会水平左右移动)，需要计算该区域内的所有窗口的直方图，以上代码中我们保持窗口大小不变(我们也可以通过搜索稍小或稍大的窗口以考虑目标对象最终的比例变化)。为了限制计算的复杂性，我们使用的直方图中包含的 bin 数量较低(仅使用 16 个 bin)。因此，此多平面图像中的的每个平面都是一个二值图像，平面 0 包含的图像显示值在 0 和 15 之间的所有像素，而平面 1 显示值在 16 和 31 之间的像素，依此类推。
对对象的搜索包括计算给定大小的所有窗口在预定像素范围内的直方图，这表示我们需要计算 3200 个不同直方图，这些直方图是从积分图像中得到的。IntegralImage 类返回的所有直方图都包含在 cv::Vec 对象中(使用 at 方法)，然后我们使用 cv::compareHist 函数来识别最相似的直方图，该函数与大多数 OpenCV 函数一样，可以通过 cv::InputArray 泛型参数类型接受 cv::Mat 或 cv::Vec 对象。
3.2 完整代码
关于头文件 histogram.h 和 integral.h 可以分别参考直方图一节和上一小节，主函数文件 tracking.cpp 代码如下：
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <vector>
#include "histogram.h"
#include "integral.h"
int main() {
// 打开图像
cv::Mat image = cv::imread("4.png", 0);
// 定义 ROI
int xo = 350, yo = 265;
int width = 25, height =35;
cv::Mat roi(image, cv::Rect(xo, yo, width, height));
// 计算像素和
cv::Scalar sum = cv::sum(roi);
std::cout << sum[0] << std::endl;
// 计算积分图像
cv::Mat integralImage;
cv::integral(image, integralImage, CV_32S);
int sumInt = integralImage.at<int>(yo+height, xo+width)
-integralImage.at<int>(yo+height, xo)
-integralImage.at<int>(yo, xo+width)
+integralImage.at<int>(yo, xo);
std::cout << sumInt << std::endl;
// 直方图
Histogram1D h;
h.setNBins(16);
cv::Mat refHistogram = h.getHistogram(roi);
cv::namedWindow("Reference Histogram");
cv::imshow("Reference Histogram", h.getHistogramImage(roi, 16));
std::cout << refHistogram << std::endl;
// 创建具有16个平面的二进制图像
cv::Mat planes;
convertToBinaryPlanes(image, planes, 16);
// 计算积分图像
IntegralImage<float, 16> intHisto(planes);
cv::Vec<float,16> histogram= intHisto(xo,yo,width,height);
std::cout<< histogram << std::endl;
cv::namedWindow("Reference Histogram (2)");
cv::Mat im= h.getImageOfHistogram(cv::Mat(histogram),16);
cv::imshow("Reference Histogram (2)",im);
// 在第二张图像中搜索
cv::Mat secondImage = cv::imread("5.png", 0);
if (!secondImage.data) return 0;
// 创建具有16个平面的二进制图像
convertToBinaryPlanes(secondImage, planes, 16);
// 计算积分图像
IntegralImage<float, 16> intHistogram(planes);
histogram = intHistogram(135, 114, width, height);
std::cout << histogram << std::endl;
cv::namedWindow("Current Histogram");
cv::Mat im2= h.getImageOfHistogram(cv::Mat(histogram),16);
cv::imshow("Current Histogram",im2);
std::cout << "Distance= " << cv::compareHist(refHistogram,histogram, cv::HISTCMP_INTERSECT) << std::endl;
double maxSimilarity = 0.0;
int xbest, ybest;
// 循环查找
for (int y=160; y<200; y++) {
for (int x=0; x<secondImage.cols-width; x++) {
// 计算直方图
histogram = intHistogram(x, y, width, height);
double distance = cv::compareHist(refHistogram, histogram, cv::HISTCMP_INTERSECT);
// 查找最相似直方图
if (distance>maxSimilarity) {
xbest = x;
ybest = y;
maxSimilarity = distance;
}
std::cout << "Distance(" << x << ", " << y << ") = " << distance << std::endl;
}
}
std::cout << "Best solution = (" << xbest << ", " << ybest << ") = " << maxSimilarity << std::endl;
// 在目标对象周围绘制矩形框
cv::rectangle(image, cv::Rect(xo, yo, width, height), 0);
cv::namedWindow("Initial Image");
cv::imshow("Initial Image", image);
cv::namedWindow("New Image");
cv::imshow("New Image", secondImage);
// 在检测到的最佳位置绘制矩形框
cv::rectangle(secondImage, cv::Rect(xbest, ybest, width, height), 0);
cv::rectangle(secondImage, cv::Rect(0, 200, secondImage.cols, height+10), 255);
cv::namedWindow("Object location");
cv::imshow("Object location", secondImage);
cv::waitKey();
return 0;
}
小结
积分图像作为在图像 (Region of Interest, ROI) 中求和像素的有效方法，已经广泛用于多尺度滑动窗口计算等应用中。本节中，介绍了积分图像的基本概念，并且通过多种不同方法实现了积分图像的计算，最后还通过两个实际应用讲解了积分图像在图像处理中的重要用途。
系列链接
OpenCV实战（1）——OpenCV与图像处理基础
OpenCV实战（2）——OpenCV核心数据结构
OpenCV实战（3）——图像感兴趣区域
OpenCV实战（4）——像素操作
OpenCV实战（5）——图像运算详解
OpenCV实战（6）——OpenCV策略设计模式
OpenCV实战（7）——OpenCV色彩空间转换
OpenCV实战（8）——直方图详解
OpenCV实战（9）——基于反向投影直方图检测图像内容