Otsu2D

/***************************************************************************** 
* 函数名称： cvOtsu2D() 
* 函数参数:  CvMat* pGrayMat:灰度图形相对应的矩阵 
* 返回值: int nThreshold
* 函数说明:实现灰度图的二值化分割——最大类间方差法（二维Otsu算法） 
* 备注：在构建二维直方图的时候，采用灰度点的3*3邻域均值 
******************************************************************************/  
int cvOtsu2D(CvMat *pGrayMat)  
{  
	double dHistogram[256][256];        //建立二维灰度直方图  
	double dTrMatrix = 0.0;				//离散矩阵的迹
	int height = pGrayMat->height;
	int width = pGrayMat->width;
	int N = height*width;               //总像素数
	int i, j;  

	for(i = 0; i < 256; i++)  
	{  
		for(j = 0; j < 256; j++)  
			dHistogram[i][j] = 0.0;      //初始化变量  
	}  

	for(i = 0; i < height; i++)  
	{  
		for(j = 0; j < width; j++)  
		{  
			unsigned char nData1 = (unsigned char)cvGetReal2D(pGrayMat, i, j);//当前的灰度值  
			unsigned char nData2 = 0;                                                                                                                                
			int nData3 = 0;//注意9个值相加可能超过一个字节                                                                   
			for(int m = i-1; m <= i+1; m++)  
			{  
				for(int n = j-1; n <= j+1; n++)  
				{  
					if((m >= 0) && (m < height) && (n >= 0) && (n < width))  
						nData3 += (unsigned char)cvGetReal2D(pGrayMat, m, n); //当前的灰度值                                                                        
				}  
			}  
			nData2 = (unsigned char)(nData3 / 9);    //对于越界的索引值进行补零,邻域均值  
			dHistogram[nData1][nData2]++;  
		}  
	}  
	for(i = 0; i < 256; i++)  
		for(j = 0; j < 256; j++)  
			dHistogram[i][j] /= N;  //得到归一化的概率分布  

	double Pai = 0.0;      //目标区均值矢量i分量  
	double Paj = 0.0;      //目标区均值矢量j分量  
	double Pbi = 0.0;      //背景区均值矢量i分量  
	double Pbj = 0.0;      //背景区均值矢量j分量  
	double Pti = 0.0;      //全局均值矢量i分量  
	double Ptj = 0.0;      //全局均值矢量j分量  
	double W0 = 0.0;       //目标区的联合概率密度  
	double W1 = 0.0;       //背景区的联合概率密度  
	double dData1 = 0.0;  
	double dData2 = 0.0;  
	double dData3 = 0.0;  
	double dData4 = 0.0;   //中间变量  
	int nThreshold_s = 0;  
	int nThreshold_t = 0;  
	double temp = 0.0;     //寻求最大值  
	for(i = 0; i < 256; i++)  
	{  
		for(j = 0; j < 256; j++)  
		{  
			Pti += i*dHistogram[i][j];  
			Ptj += j*dHistogram[i][j];  
		}  
	}  
	for(i = 0; i < 256; i++)  
	{  
		for(j = 0; j < 256; j++)  
		{  
			W0 += dHistogram[i][j];  
			dData1 += i*dHistogram[i][j];  
			dData2 += j*dHistogram[i][j];  

			W1 = 1-W0;  
			dData3 = Pti-dData1;  
			dData4 = Ptj-dData2;  

			Pai = dData1 / W0;  
			Paj = dData2 / W0;  
			Pbi = dData3 / W1;  
			Pbj = dData4 / W1;   // 得到两个均值向量，用4个分量表示  
			dTrMatrix = ((W0 * Pti - dData1) * (W0 * Pti - dData1) + (W0 * Ptj - dData2) * (W0 * Ptj- dData2)) / (W0 * W1);  
			if(dTrMatrix > temp)  
			{  
				temp = dTrMatrix;  
				nThreshold_s = i;  
				nThreshold_t = j;  
			}  
		}  
	}

	int nThreshold = (nThreshold_s + nThreshold_t) / 2;//返回阈值，有多种形式，可以单独某一个，也可是两者的均值 
	return nThreshold; 
}