最大熵阈值分割代码

理论在这里

http://blog.youkuaiyun.com/happygirlliu520/article/details/47026649

自己写的代码

//最大熵阈值分割
float opencv::calEntropy(const Mat& src)
{
ASSERT(CV_8UC1 == src.type());
double pix_probility[256] = { 0 };
int pix_counter[256] = { 0 };
size_t sum = static_cast<size_t>(src.cols*src.rows); //分辨率
size_t max_threshold = 255, min_threshold = 0;
uchar *idata = src.data;
for (size_t i = 0; i < sum; i++ )
{
pix_counter[(uchar)*idata] += 1;
idata++;
}

for (; max_threshold >= 0; max_threshold--)
{
if (0 != pix_counter[max_threshold]) break;
}


for (; min_threshold <256; min_threshold++)
{
if (0 != pix_counter[min_threshold]) break;
}

size_t threshold = 0;
double entropy = 0;
double pt =  0; //前景概率总和
double H = 0;
for (size_t i = min_threshold; i <= max_threshold; i++)
{
if (0<pix_counter[i])
{
pix_probility[i] = (double)pix_counter[i] / sum;
H -= pix_probility[i] * logf(pix_probility[i]);
}
}
double Ht = 0; 
for (size_t i = min_threshold; i <= max_threshold; i++) 
{
if (0 >= pix_probility[i]) continue;
pt += pix_probility[i];
Ht -= pix_probility[i] * logf(pix_probility[i]);
double p = pt*(1 - pt);
double l_entropy = logf(p) + (Ht + H*pt - 2*Ht*pt) / (p);//图像总熵，化简以后的式子
if (entropy < l_entropy)
{
entropy = l_entropy; 
threshold = i;
}
}
return (float)threshold;
}