OpenCV图像处理——积分图计算与积分二值化

1.积分图原理

积分图（Integral Image）是计算机视觉中一个非常有用的数据结构，特别是在图像处理和特征提取领域。它主要用于快速计算图像中任意矩形区域的像素和。这个特性使得积分图在诸如目标检测、运动估计和图像分割等任务中非常有用。
积分图（Integral Image）的定义：取图像左上侧的全部像素计算累加和，并用这个累加和替换图像中的每一个像素，使用这种方式得到的图像称为积分图像。

$$I(x,y)=\sum _{\begin{array}{c}{x}^{\prime }⩽x\\ y^{\prime }⩽y\\ \end{array}}i(x^{\prime },y^{\prime })$$

积分图又称总和面积表（summed area table，简称SAT），是一个快速且有效的对一个网格的矩形子区域中计算和的数据结构和算法。

积分图可以只遍历一次图像即可有效的计算出来，其通常被用来加速计算过程。一旦积分图计算完毕，对任意矩形区域的和的计算就可以在常数时间内（一次加法，两次减法）完成。

2. cv::integral

cv::integral 有三种重载形式：

//第一种重载形式
void integral(
    InputArray src,
    OutputArray sum,
    int sdepth = -1
);

//第二种重载形式
void integral(
    InputArray src,
    OutputArray sum,
    OutputArray sqsum,
    int sdepth = -1,
    int sqdepth = -1
);

//第三种重载形式
void integral(
    InputArray src,
    OutputArray sum,
    OutputArray sqsum,
    OutputArray tilted,
    int sdepth = -1,
    int sqdepth = -1
);

参数详解：

1. InputArray src：大小为$W\times H$的输入图像。类型为8-bit或浮点型（32f或64f）。

2. OutputArray sum：输出大小为$(W+1)\times (H+1)$的积分图，类型为32-bit整型或浮点型（32f或64f）。
   

3. OutputArray sqsum：对像素值进行平方之后再计算得到的积分图，积分图大小是$(W+1)\times (H+1)$，类型为双精度浮点型（64f）。
   

4. OutputArray tilted：对原始图像旋转45度后再计算得到的积分图，积分图的大小依然是$(W+1)\times (H+1)$，类型和 sum 一样。
   

5. int sdepth = -1： sum 和 tilted 的位图深度，可以是CV_32S、CV_32F或CV_64F。

6. int sqdepth = -1： sqsum 的位图深度，可以是CV_32F或CV_64F。
   
   上图是一个例子。上图左是Rect(4,4,3,2)矩形内的积分计算，上图右是一个倾斜矩形Rect(5,1,2,3)内的积分计算。

3.代码测试：

Mat src = (Mat_<uchar>(2,3)<<1,2,3,4,5,6);
for(int r =0 ; r<src.rows; r++)
{
    for(int c = 0 ; c<src.cols; c++)
    {
        int i = src.ptr<uchar>(r)[c];
        std::cout<<i<<" ";
    }
    std::cout<<std::endl;
}
/*
1 2 3
4 5 6
*/
for(int r =0 ; r<src.rows; r++)
{
    for(int c = 0 ; c<src.cols; c++)
    {
        int i = src.ptr<uchar>(r)[c];
        std::cout<<i*i<<" ";
    }
    std::cout<<std::endl;
}
/*
1 4 9
16 25 36
*/

Mat sumii = Mat::zeros(src.rows + 1, src.cols + 1, CV_32FC1);
Mat sqsumii = Mat::zeros(src.rows + 1, src.cols + 1, CV_64FC1);
Mat tilt = Mat::zeros(src.rows + 1, src.cols + 1, CV_32FC1);
integral(src, sumii, sqsumii, tilt);

for(int r =0 ; r<sumii.rows; r++)
{
    for(int c = 0 ; c<sumii.cols; c++)
    {
        int i = sumii.ptr<int>(r)[c];//用float或double输出的数值不对
        std::cout<<i<<" ";
    }
    std::cout<<std::endl;
}
/*
0 0 0 0
0 1 3 6
0 5 12 21
*/
for(int r =0 ; r<sqsumii.rows; r++)
{
    for(int c = 0 ; c<sqsumii.cols; c++)
    {
        double i = sqsumii.ptr<double>(r)[c];
        std::cout<<i<<" ";
    }
    std::cout<<std::endl;
}
/*
0 0 0 0
0 1 5 14
0 17 46 91
*/
for(int r =0 ; r<tilt.rows; r++)
{
    for(int c = 0 ; c<tilt.cols; c++)
    {
        int i = tilt.ptr<int>(r)[c];
        std::cout<<i<<" ";
    }
    std::cout<<std::endl;
}
/*
0 0 0 0
0 1 2 3
1 7 11 11
*/

tilt 就是计算旋转45度后矩形的面积，比如：

积分二值化

主要步骤如下：

1. 检查输入图像：首先，确保输入的图像不为空且是 8 位深度（即灰度图像或彩色图像）。
2. 灰度转换：如果输入图像是彩色的，则转换为灰度图像，否则直接使用原图。
3. 积分图像计算：通过 cv::integral 计算输入图像的积分图像，这在之后的局部区域求和时使用。
4. 设置参数：定义了局部窗口大小 S 和阈值系数 T，其中 S 是图像宽高的 1/8，T 用于控制阈值的敏感度。
5. 区域扫描与阈值计算：
   • 对每个像素点，基于局部的 SxS 窗口计算图像在该区域内的积分和（求和）。
   • 然后通过对比像素值与该区域内像素值的均值（结合系数 T 和输入的全局阈值 thre）决定该像素是否置为 0 或 255（即黑或白）。
6. 输出结果：生成二值化图像 cv_dst。

这个自适应阈值方法比全局阈值更加灵活，特别适合处理光照不均匀的图像。

void threshold_integral(cv::Mat& cv_src, double thre, cv::Mat& cv_dst)
{
	// accept only char type matrices
	CV_Assert(!cv_src.empty());
	CV_Assert(cv_src.depth() == CV_8U);

	cv::Mat cv_bin;

	if (cv_src.channels() != 1)
	{
		cv::cvtColor(cv_src, cv_bin, cv::COLOR_BGR2GRAY);
	}
	else
	{
		cv_bin = cv_src.clone();
	}

	cv_dst = cv::Mat(cv_bin.size(), CV_8UC1, 1);

	// rows -> height -> y
	int nRows = cv_bin.rows;
	// cols -> width -> x
	int nCols = cv_bin.cols;

	// create the integral image
	cv::Mat sumMat;
	cv::integral(cv_bin, sumMat);

	CV_Assert(sumMat.depth() == CV_32S);
	CV_Assert(sizeof(int) == 4);

	int S = MAX(nRows, nCols) / 8;
	double T = 0.15;

	// perform thresholding
	int s2 = S / 2;
	int x1, y1, x2, y2, count, sum;

	// CV_Assert(sizeof(int) == 4);
	int* p_y1, * p_y2;
	uchar* p_inputMat, * p_outputMat;

	for (int i = 0; i < nRows; ++i)
	{
		y1 = i - s2;
		y2 = i + s2;

		if (y1 < 0)
		{
			y1 = 0;
		}
		if (y2 >= nRows)
		{
			y2 = nRows - 1;
		}

		p_y1 = sumMat.ptr<int>(y1);
		p_y2 = sumMat.ptr<int>(y2);
		p_inputMat = cv_bin.ptr<uchar>(i);
		p_outputMat = cv_dst.ptr<uchar>(i);

		for (int j = 0; j < nCols; ++j)
		{
			// set the SxS region
			x1 = j - s2;
			x2 = j + s2;

			if (x1 < 0)
			{
				x1 = 0;
			}
			if (x2 >= nCols)
			{
				x2 = nCols - 1;
			}

			count = (x2 - x1) * (y2 - y1);

			// I(x,y)=s(x2,y2)-s(x1,y2)-s(x2,y1)+s(x1,x1)
			sum = p_y2[x2] - p_y1[x2] - p_y2[x1] + p_y1[x1];

			if ((int)(p_inputMat[j] * count) < (int)(sum * (1.0 - T) * thre))
				p_outputMat[j] = 0;
			else
				p_outputMat[j] = 255;
		}
	}
}