/*****************************************************************************
* 函数名称: 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;
}