最大熵阈值分割法

1.频率和概率

    直方图每个矩形框的数值描述的是图像中相应灰度值的频率。因此，可以说直方图是一种离散的频率分布。给定一个大小为M*N的图像I，直方图中所有矩形框所代表的数值之和，即为图像中的像素数量，即：

    相对应的归一化直方图表示为：

   0<=i<K 通常被解释为一个随机过程的概率分布或概率密度函数，表示的是图像中像素灰度值为i所出现的概率。i的累积概率值为1，即概率分布p必须满足以下关系:

    与累积概率所所对应的累积直方图H是一个离散的分布函数P(),(通常也称为累积分布函数或cdf)：

2.最大熵阈值分割

    熵是信息理论中一个重要的概念，这种方法常用于数据压缩领域。熵是一种统计测量方法，用以确定随机数据源中所包含的信息数量。例如，包含有N个像素的图像I，可以解释为包含有N个符号的信息，每一个符号的值都独立获取于有限范围K（e.g.,256）中的不同灰度值。

   将数字图像建模为一种随机信号处理，意味着必须知道图像灰度中每一个灰度g所发生的概率，即：

  因为所有概率应该事先知道，所以这些概率也称为先验概率。对于K个不同灰度值g=0,…,K-1的概率向量可以表示为：

 上述概率向量也称为概率分布或称为概率密度函数(pdf)。实际数字图像处理应用当中，先验的概率通常是不知道的，但是这些概率可以通过在一幅图像或多幅图像中观察对应灰度值所发生的频率，从而估算出其概率。图像概率密度函数p(g)可以通过归一化其对应的直方图获得其概率，即：

3.最大熵

    数字图像中给定一个估算的概率密度函数p(g)，数字图像中的熵定义为：

4.用图像熵进行图像分割

    利用图像熵为准则进行图像分割有一定历史了，学者们提出了许多以图像熵为基础进行图像分割的方法。我们介绍一种由Kapuret al提出来，现在仍然使用较广的一种图像熵分割方法。

给定一个特定的阈值q(0<=q<K-1),对于该阈值所分割的两个图像区域C0,C1，其估算的概率密度函数可表示为：

 同样的，背景熵可以改写为：

5.代码实现及结果

int MakeTable(double *S1,double *S2,double *NormalizeHist,int HistSize)
{	int i;

	double s0 = 0,s1;
	for (i=0;i < HistSize;i++)
	{
		if (NormalizeHist[i]>0)
		{
			 s0 = s0+NormalizeHist[i]*log(NormalizeHist[i]);
		}
		S1[i] = s0;
	}

	s1 = 0;
	for (i = HistSize-1;i >= 0;i--)
	{
		S2[i] = s1;
		if (NormalizeHist[i]>0)
		{
			s1 = s1+NormalizeHist[i]*log(NormalizeHist[i]);
		}
	}
	return 0;
}

int Get_Max_Entropy_Threshold(T_U32 *Hist,DWORD width,DWORD height)
{
	int i,q,HistSize = 256,qMax = -1;
	double NormalizeHist[256] = {0},H01 = 0,H0 = 0,H1 = 0,P1 = 0,P0 = 0, Hmax = -9999 ;
	double S1[256] = {0},S2[256] = {0};

	for (i = 0;i < HistSize;i++)
		NormalizeHist[i] = Hist[i]/(double)(width*height);

	MakeTable(S1,S2,NormalizeHist,HistSize);

	for (q = 0; q < HistSize-1;q++)
	{
		P0 += NormalizeHist[q];
		P1 = 1 - P0;

		if(P0 > 0)
			H0 = -(1*S1[q])/P0+log(P0);
		else
			P0 = 0;

		if(P1 > 0)
			H1 = -(1*S2[q])/P1+log(P1);
		else
			H1 = 0;

		H01 = H0+H1;

		if (H01 > Hmax)
		{
			Hmax = H01;
			qMax = q;
		}

	}
	return qMax;
}
int Max_Entropy_Threshold(IMAGE_TYPE *BMP8_img)
{
	DWORD width,height,i,j,bfsize;
	WORD  biBitCount;
	T_U8 *dst,*bmp,k;
	T_U32 line_byte,Threshold,index,Hist[256] = {0};
	
	BITMAPFILEHEADER bf;
	BITMAPINFOHEADER bi;
	
	FILE *Max_Entropy_ThresholdBMP_fp = fopen("Max_Entropy_Threshold.bmp","wb");
	
	if(NULL == Max_Entropy_ThresholdBMP_fp)
	{
		printf("Can't open IsoData_Threshold.bmp\n");
		return -1;
	}
	bmp = BMP8_img;
	memset(&bf, 0, sizeof(bf));
	memset(&bi, 0, sizeof(bi));
	memcpy(&bf,bmp,14);
	memcpy(&bi,&bmp[14],40);
	
	height = bi.biHeight;
	width  = bi.biWidth;
	bfsize = bf.bfSize;
	biBitCount = bi.biBitCount;
	line_byte = (width * 8 / 8 + 3) / 4 * 4; 
	dst = BMP8_img+54+1024;

	for (i = 0; i < height;i++)
	{
		for (j = 0;j < width;j++)
		{
			Hist[dst[i*line_byte+j]]++;
		}
	}
	

	dst = BMP8_img+54+1024;
	Threshold = Get_Max_Entropy_Threshold(Hist,width,height);
	for(i = 0;i < height;i++)
	{
		for(j = 0;j < width;j++)
		{

			index = i*line_byte+j;
			k = dst[index];
			
			if(k >= Threshold)
				dst[index] = 255;
			else
				dst[index] = 0;
			
		}
	}
	
	fwrite(BMP8_img,bfsize,1,Max_Entropy_ThresholdBMP_fp);
	fclose(Max_Entropy_ThresholdBMP_fp);
	Max_Entropy_ThresholdBMP_fp = NULL;
	
	return 0;

}

结果：