图像自适应阈值（最大类间方差法、迭代法、矩形区域块自适应阈值（使用积分处理的自适应处理)）

一、最大类间方差法 或者 大津法 或者 OTSU

1、简介

最大类间方差法是由日本学者大津(Nobuyuki Otsu)于1979年提出的,是一种自适应的阈值确定的方法,又叫大津法,简称OTSU。它是按图像的灰度特性,将图像分成背景和目标两部分，或者说，是寻找一个阈值为K，将图像的颜色分为1,2.....K和K+1.....256两部分。

如何确定这个阈值K？算法分类的原理是让背景和目标之间的类间方差最大，因为背景和目标之间的类间方差越大,说明构成图像的2部分的差别越大,错分的可能性越小。下面进行公式推导：

2、最大类间方差法推到公式：

M*N 为一张图像所有的像素点。N1为小于等于当前像素的所有点数目，w1为 小于等于当前像素的所有点数占所有像素点的比。

N2为大于当前像素点的所有点的数目，w2为 大于当前像素点的所有点的数目的占比。u1 为所有小于等于当前像素点值的和 除上N1，也就是小于等于当前像素值所有数的平均值，u2为大于当前像素所有数的平均值，u为以当前像素为分割点时的期望（平均值），g为方差。

代码：

int ThresholdByOtsu( stImage* pImage)           // 最大类间方差法求阈值。
{
	int iWidth = pImage->cols;
	int iHeight = pImage->rows;
	unsigned char* pGrayImg = pImage->buffer;

	// verify the input parameters
	if ((pGrayImg == 0) || (iWidth <= 0) || (iHeight <= 0))
	{
		return -1;
	}
	int thresholdValue = 0; // Threshold
	int n, n1, n2;
	double m1, m2, sum, csum, fmax, sb;
	memset(ihist, 0, sizeof(ihist)); // set ihist[] to zero 
	n = iHeight*iWidth;
	sum = csum = 0.0;
	fmax = -1.0;
	n1 = 0;
	// hist generation
	for (int i = 0; i < iHeight; i++)
	{
		for (int j = 0; j < iWidth; j++)
		{
			ihist[*pGrayImg]++;         //统计每个像素值的个数。
			pGrayImg++;
		}
	}
	pGrayImg -= n;                           // 等于 pGrayImg = pImage->buffer;
	for (int k = 0; k <= 255; k++)
	{
		sum += (double)k * (double)ihist[k];         // 像素值乘以当前个数。
	}
	// Otsu Algorithm
	for (int k = 0; k <= 255; k++)
	{
		n1 += ihist[k];                                          //n1为小于等于当前像素的个数
		if (n1 == 0)continue;
		n2 = n - n1;                                             // n2 为当前当前像素的个数。
		if (n2 == 0)break;

		csum += (double)k *ihist[k];
		m1 = csum / n1;         // average level of class 1       // m1 遍历到当前像素值得，平均值水平。

		m2 = (sum - csum) / n2; // average level of class 2        //m2 剩下的（大于当前）像素值的平均值水平。

		sb = (double)n1 *(double)n2 *(m1 - m2) * (m1 - m2);        // 这里虽然是用的n1和n2不是公式中的w1、w2(比重).其实最终相求的，等同于w1、w2都扩大了iHeight*iWidth这么多倍。
		if (sb > fmax)
		{
			fmax = sb;
			thresholdValue = k;
		}
	}
	return(thresholdValue);
}

二、迭代法求图像阈值

迭代法阈值选择算法是对双峰法的改进，他首先选择一个近似的阈值T，将图像分割成两个部分，R1和R2，计算出区域R1和R2的均值u1和u2，再选择新的

阈值T=(u1+u2)/2;

重复上面的过程，知道u1和u2不在变化为止。

由于迭代法比较易实现，在这就不写代码了。

三，矩阵区域块自适应阈值（使用积分阈值自适应处理）

 摘要：  图像阈值处理是许多计算机视觉和图形应用中的常见任务。 对图像进行阈值处理的目的是将像素分类为“暗”或“亮”。 自适应阈值处理是一种阈值处理，它考虑了照明的空间变化。 我们提出了一种使用输入的积分图像进行实时自适应阈值处理的技术。 我们的技术是以前方法的扩展。 但是，我们的解决方案对图像中的照明变化更加稳健。 此外，我们的方法简单易行。我们的技术适用于以实时帧速处理实时视频流，使其成为增强现实等交互式应用的有用工具。

1、简介（introduction）
图像阈值化基于像素的特定特征（例如强度值）来分割数字图像。目标是创建二值图像，将每个像素分类为两个类别之一，例如“暗”或“亮”。这是许多图像处理应用程序和一些计算机图形应用程序中的常见任务。例如，它通常是基于标记的增强现实系统的第一步[Billinghurstet al。 2001;布拉德利和罗斯2004; Fiala 2005]，它已被用于高动态范围摄影[Ward2003]。最基本的阈值处理方法是选择固定的阈值并将每个像素与该值进行比较。这些技术已在许多调查论文中进行了广泛的描述和评估[Weszka andRosenfeld 1978; Palumbo等人。 1986; Sahoo等。 1988;李等人。 1990; Glasbey 1993;特里尔和耆那教1995; Sezginand Sankur 2004]。然而，如果照明在图像中在空间上变化或在视频流中超时，则固定阈值经常失败。

为了解释照明的变化，常见的解决方案是自适应阈值处理。主要区别在于为图像中的每个像素计算不同的阈值。该技术为照明变化提供了更强大的功能。存在许多自适应阈值方法[White和Rohrer 1983; Bernsen1986;帕克1991; Wellner 1993;杨等人。 1994; Shen和Ip 1997;陈等人。 1998; Savakis 1998; Sauvolaand Pietikainen 2000;杨和燕2000]。进一步的例子和比较可以在[Venkateswarlu andBoyle 1995; Sezgin和Sankur 2004]。我们使用积分图像提出了一种非常简单明了的技术。我们的方法易于实现实时视频流的实时性能。虽然我们的技术是以前方法的扩展[Wellner 1993]，但我们增强了强光照变化的稳健性。此外，我们提供清晰整洁的解决方案，而不会增加实施的复杂性。我们的技术也类似于White和Rohrer用于光学字符识别的阈值处理方法[White and Rohrer 1983]，但我们提出了一种针对实时视频设计的实现方法。这项工作的动机是找到真正的现实应用程序。 Pintaric还提出了一种适用于增强现实标记的自适应阈值算法[Pintaric 2003]，但是他的方法要求在前一帧中定位一个元素，以便技术正确地进行阈值处理。我们的算法不做任何假设，更通用，适合在任何应用中使用。源代码可在本文末尾列出的地址在线获取。

2、背景 

2.1实时自适应阈值处理
在本文中，我们专注于从实时视频流中自适应地阈值化图像。 为了保持实时性能，阈值算法必须限制通过每个图像的迭代次数。 阈值处理通常是一个子任务，它构成了一个更大的过程的一部分。 例如，在增强现实中，必须对输入图像进行分段以定位场景中用于动态建立相机设置的已知标记。 因此，简单快速的自适应阈值技术是一个重要的工具

2.2积分图像
积分图像（也称为summed-area table）是一种工具，当我们有一个从像素到数值（比如灰度值）的映射函数f(x, y)，并且希望计算一个矩形区域内的函数值的和时，积分图是一个非常高效的工具。已应用积分图像的示例包括纹理映射[Crow 1984]，图像中的人脸检测[Viola and Jones 2004]和立体对应[Veksler 2003]。在没有积分图像的情况下，我们计算每个矩形窗口中像素点的映射函数的和，需要通过分别计算每个像素的函数值最后叠加实现。但是，如果我们需要计算多个重叠矩形窗口的总和，为了降低时间复杂度和操作次数，我们可以使用积分图像。

为了计算积分图像，我们在每个位置I（x，y）存储左边和上面像素（x，y）的所有f（x，y）项的总和。 对于每个像素，使用以下等式在线性时间内完成（考虑边界情况），

图2（左和中）说明了积分图像的计算。 一旦我们得到积分图像，任何具有左上角（x1，y1）和右下角（x2，y2）的矩形的函数的总和可以使用以下等式计算，计算时间不稳定。

图2（右）说明使用等式2计算矩形D上的f（x，y）之和相当于计算矩形上的和（A + B + C + D） - （A + B） - （A+ C）+ A。

图2：积分图像。 左：图像值的简单输入。 中心：计算的积分图像。 右：使用积分图像计算矩形D的总和。

3、技术                         ----> 这个是最重要的。

    我们的自适应阈值技术是Wellner方法的简单扩展[WELNER 1993 ]。Welnne算法的主要思想是将每个像素与周围像素的平均值进行比较。特别地，在遍历图像时计算了最后一个像素的近似移动平均值。如果当前像素的值比平均值低百分之t，那么它被设置为黑色，否则它被设置为白色。这种方法的工作是将像素与附近像素的平均值进行比较，从而保持硬对比线，忽略软梯度变化。这种方法的优点是只需要遍历一次图像。WELNER使用1／8倍的图像宽度来表示s的值，t=15。然而，这种方法的一个问题是它依赖于像素的扫描顺序。此外，由于邻域样本在各个方向上不均匀分布，所以在使用平均步长并不是一个很好的表现形式。通过使用积分图像（和骶一个额外的迭代通过图像），我们提出了一个解决方案，不遭受这些问题。我们的技术是干净的，直截了当的，易于编码，并产针对不同的处理方式产生生相同的输出。我们不是计算最后看到的s个像素的平均值，而是计算以每个像素为中心的s x s像素窗口的平均值。 这是比较好的平均值，因为它考虑了所有边上的相邻像素。利用积分图像进行平均计算的时间是线性的。我们在第一次遍历图像的时候计算积分图像。在第二遍中，我们利用积分图像计算每个像素的s x s的平均值，使用时间为常量，然后进行比较。如果当前像素的值小于这个平均值，则将其设置为黑色，否则它被设置为白色。

代码：

void AdaptiveThresholdPartOptUsingSram(unsigned char * org_img, unsigned char * new_img, int width, int height, int block, int ratio)
{
	//org_img 原图   new_img 二值化后的图   block块略等于1/8宽度 ratio？
	int S = (block - 1) / 2;//101
	int T = 20; //15
	T = (ratio >= 100 || ratio <= 0) ? 20 : ratio;

	int i, j;
	long sum = 0;
	int count = 0;
	int index = 0;
	int x1, y1, x2, y2;
	// 计算积分图像    别说积分好不好，挺吓人的，不就是求dp[i][j]吗，dp[i][j] 为从（0,0）为起点到（i,j）为终点的矩形内，所有点之和。
	Lx_2Row[0][0] = Lx_1Col[0] = org_img[0];
	// 原来这就是图像积分。
	for (i = 1; i<width; i++)
	{
		Lx_2Row[0][i] = Lx_2Row[0][i - 1] + org_img[i];      // 横向累加  第一行
	}

	for (j = 1; j<height; j++)
	{
		Lx_1Col[j] = Lx_1Col[j - 1] + org_img[j*width];      //  第一列 横向累加。        
	}
	
	memcpy(&integralImg[0][0], &Lx_2Row[0][0], width * sizeof(unsigned int));

	for (j = 1; j<height; j++)
	{
		int row0 = (j) & 0x01, row1 = (j + 1) & 0x01;
		Lx_2Row[row0][0] = Lx_1Col[j];
		index = j*width + 1;
		for (i = 1; i<width; i++)
		{
			Lx_2Row[row0][i] = Lx_2Row[row0][i - 1] + Lx_2Row[row1][i] - Lx_2Row[row1][i - 1] + org_img[index++];
			// 其实这句话的推到公式可以简化为 dp[i][j] = dp[i][j-1] + dp[i-1][j] - dp[i-1][j] + map[i][j];
			// dp[i][j] 为 以0,0为点，以i、j终点的矩形中所有点的和。
			// 在第一行和一列已经累加出来的前提下。
		}
		memcpy(&integralImg[j][0], &Lx_2Row[row0][0], width * sizeof(unsigned int));
	}
	// integralImg 积分后的图像。

	memset(new_img, 255, (height*width));    //放这输入输出图可以用同一帧。全是255为了中值后图片与固定阈值分割的图逻辑与操作使用

	count = block*block * 100;
	for (i = (S + 1); i< (height - S); i++)
	{
		for (j = (S + 1); j< (width - S); j++)
		{
			index = i*width + j;

			// set the SxS region
			x1 = j - S;
			x2 = j + S;
			y1 = i - S;
			y2 = i + S;

			sum = integralImg[y2][x2] - integralImg[y1][x2] - integralImg[y2][x1] + integralImg[y1][x1];
			// 求以(y1,x1)为起点，（y2,x2）为终点的 矩阵所有点之和。

			//bin[index] = (((unsigned int)( input[index] * count- sum*(100 - T) ))&&(1)) * 255;

			if ((org_img[index] * count) < (sum*(100 - T)))  //不懂这句话的，请看上文技术，如果当前像素的值比平均值低百分之t，那么它被设置为黑色，否则它被设置为白色       
				new_img[index] = 0;
			//else
			//	new_img[index] = 255;

		}
	}
}

四、就是给图像加边框后，再进行块状区域积分求阈值

static int SurroundImage(unsigned char* src, unsigned char* dst, int width, int height, int exWidth, int exHeight)
{
	//原图像是data；扩展后是DataDst width、height为原图的宽高，exwidth、exheight为块状区域的宽高
	int i, j;
	unsigned char* pSrc, *pDst;
	int dstWidth = width + exWidth;
	int dstHeight = height + exHeight;

	if ((0 != width % 2) || (0 != height % 2))
	{
		memset(src, 0, (height*width));
		return -1;
	}

	//center
	pSrc = src;
	pDst = dst + (exHeight / 2)*(width + exWidth) + exWidth / 2;
	for (j = exHeight / 2; j < (height + exHeight / 2); j++)   // 把原图复制过来。
	{
		memcpy(pDst, pSrc, (width));             //
		pSrc += (width);
		pDst += (width + exWidth);
	}

	//  9月17号 打卡
	//  9月18号 继续？ 继续
	//4 corners /*加4个边角*/  
	pSrc = src;
	pDst = dst;
	/*左上角的方框，全部赋值为原图的左上角点*/
	for (i = 0; i < (exWidth / 2); i++)
		for (j = 0; j < (exHeight / 2); j++)
			*(pDst + i + j*(dstWidth)) = *(pSrc + 0 + 0 * width);//pSrc[0][0];   // dstWidth  = width + exWidth
	/*右上角的方框，全部赋值为原图的右上角点*/
	for (i = (width + exWidth / 2); i < (width + exWidth); i++)
		for (j = 0; j < (exHeight / 2); j++)
			*(pDst + i + j*(dstWidth)) = *(pSrc + width - 1 + 0 * width); //pSrc[width-1][0];

	/*右下角的方框，全部赋值为原图的右下角点*/
	for (i = (width + exWidth / 2); i < (width + exWidth); i++)
		for (j = (height + exHeight / 2); j < (height + exHeight); j++)    // 行数
			*(pDst + i + j*(dstWidth)) = *(pSrc + width - 1 + (height - 1) * width); ;// pSrc[width - 1][height - 1];

	/*左下角的方框，全部赋值为原图的左下角点*/
	for (i = 0; i < (exWidth / 2); i++)
		for (j = (height + exHeight / 2); j < (height + exHeight); j++)
			*(pDst + i + j*(dstWidth)) = *(pSrc + 0 + (height - 1) * width); //pSrc[0][height-1];

																			 //4line
	pSrc = src;
	pDst = dst;
	//up
	/* 上方的线，都用原图中的第一行赋值*/
	for (i = 0; i < (width); i++)
		for (j = 0; j < (exHeight / 2); j++)
		{
			*(pDst + i + exWidth / 2 + j*(dstWidth)) = *(pSrc + i + 0 * width); //pSrc[i][0];
		}
	//right
	/*右边的线exWidth/2 * height 这么多点，有height行，每行exWidth/2个点，全部用原图最右边的线上对应行的那个点赋值*/
	for (j = 0; j < (height); j++)
	{
		memset((pDst + width + exWidth / 2 + (j + exHeight / 2)*(dstWidth)), *(pSrc + (j + 1)*width - 1), (exWidth / 2)); //pSrc[][];
	}
	//down

	for (i = 0; i < (width); i++)
		for (j = (exHeight / 2 + height); j < (exHeight + height); j++)
		{
			*(pDst + i + exWidth / 2 + j*(dstWidth)) = *(pSrc + width*(height - 1) + i); //pSrc[i][0];
		}
	//left
	for (j = 0; j < (height); j++)
	{
		memset((pDst + 0 + (j + exHeight / 2)*(dstWidth)), *(pSrc + j*width), (exWidth / 2)); //pSrc[][];
	}
	/*关于为什么这么赋值，我想应该经过验证的。*/
	return 0;
}

void BinAdaptiveThresholdFull(unsigned char* input, int height, int width, unsigned char* bin, int block, int ratio)
{
	//input 原图，width、height原图的宽高，bin要输出的图像，block 块状区域大小，ratio根据这个值，调整T的。
	//ratio can be modified : 20 means 20%
	int T = 20; //15
	int halfBlock = (block - 1) / 2;
 
	long sum = 0;
	int count = 0;
	int index = 0;
	int x1, y1, x2, y2;

	if (0 == block % 2)//block must be odd
		return;

	T = (ratio >= 100 || ratio <= 0) ? 20: ratio;
	count = block*block * 100;

	SurroundImage(input, (unsigned char*)pSurroundImage, width, height, block - 1, block - 1);
#if (!DSP_HARDWARE_) &&( !UNIT_TEST_VERSION)
//IplSurroundImage = cvCreateImage(cvSize(width + block - 1 , height + block - 1), IPL_DEPTH_8U, 1);//debug
//IplSurroundImage->imageData = pSurroundImage;
#endif
	// 计算积分图像	
	integralFullImg[0] = pSurroundImage[0];
	//Y=0
	for (int i = 1; i<(width + block - 1); i++)
	{
		int j = 0;
		integralFullImg[i + j*(width + block - 1)] = integralFullImg[i - 1 + j*(width + block - 1)] + pSurroundImage[i + j*(width + block - 1)];
	}
	//X=0
	for (int j = 1; j<(height + block - 1); j++)
	{
		int i = 0;
		integralFullImg[i + j*(width + block - 1)] = integralFullImg[i + (j - 1)*(width + block - 1)] + pSurroundImage[i + j*(width + block - 1)];
	}

	for (int j = 1; j<(height + block - 1); j++)
	{
		for (int i = 1; i<(width + block - 1); i++)
		{
			integralFullImg[i + j*(width + block - 1)] = integralFullImg[i - 1 + j*(width + block - 1)] + integralFullImg[i + (j - 1)*(width + block - 1)] - integralFullImg[i - 1 + (j - 1)*(width + block - 1)] + pSurroundImage[i + j*(width + block - 1)];
		}
	}

	for (int i = 0; i< (width); i++)
	{
		for (int j = 0; j< (height); j++)
		{
			// set the (block-1)*(block-1) region
			x1 = i;
			x2 = i + block - 1;
			y1 = j;
			y2 = j + block - 1;

			sum = integralFullImg[x2 + y2*(width + block - 1)] - integralFullImg[x2 + y1*(width + block - 1)] - integralFullImg[x1 + y2*(width + block - 1)] + integralFullImg[x1 + y1*(width + block - 1)];

			//bin[index] = (((unsigned int)( input[index] * count- sum*(100 - T) ))&&(1)) * 255;
			if ((input[i + width*j] * count) < (sum*(100 - T)))
				bin[i + width*j] = 0;
			else
				bin[i + width*j] = 255;
		}
	}
}