最大熵阈值分割——opencv与matlab实现


这篇文章的opencv实现代码源于http://blog.csdn.net/yangtrees/article/details/8785377这篇博客,自己稍加改动,再此对博主表示感谢。

1.最大熵原理

原理自己还是不是很清楚,-p*log(p)当p增大的时候log(p)减小,那么随着阈值的增加,前景(或者背景)的熵值的变化应该不是单调的(没有证明,只是自己感觉)。但我试着求出随着阈值变化的背景熵时,发现基本是单调变化的。额,写到这里突然发现高数都白学了,求导结果为-(1+log(p)),这个是大于0的,也就是会一直单调增加,无极值。

2.opencv代码

  下面给出完整的(有点冗余)。
#include <opencv2\opencv.hpp>
#include<opencv2\imgproc\imgproc.hpp>
#include <iostream>
#include <string>
#include <windows.h>
#include "Plate.h"
#include <stdlib.h>
#include <ctime>
#define window_name "均衡化后直接二值化,呵呵"
using namespace cv;
using namespace std;


typedef enum {back,object} entropy_state;  
float total;  
//绘制hist;  
Mat drawHist(Mat hist,int bins,int height,Scalar rgb)  
{  
    double maxVal=0;  
    minMaxLoc(hist,0,&maxVal,0,0);  
    int scale=1;  
    Mat histImg = Mat::zeros(height, bins, CV_8UC3);  
    float *binVal = hist.ptr<float>(0);  
    for (int i=0; i<bins; i++)  
    {  
        int intensity = cvRound(binVal[i]*height/maxVal);  
        rectangle(histImg,Point(i*scale,0),  
            Point((i+1)*scale,(intensity)),rgb,CV_FILLED);  
    }  
    flip(histImg,histImg,0);  
    return histImg;  
}  
//计算直方图;  
Mat Hist(const Mat& src)  
{  
    Mat hist;  
    int bins=256;  
    int histSize[] = {bins};  
    float range[] = {0,256};  
    const float* ranges[] = {range};  
    int channels[] = {0};  
    calcHist(&src,1,channels,Mat(),hist,1,histSize,ranges,true,false);  
    Mat histImg = drawHist(hist,bins,200,Scalar(255,0,0));  
    imshow("histRGB",histImg);  
    return hist;  
}  
//计算当前熵;  
float calEntropy(const Mat& hist,int threshold)  
{  
    float total_back=0,total_object=0;  
    float entropy_back=0,entropy_object=0;  
    float entropy = 0;  
    int i=0;  
    const float* hist_p = (float*) hist.ptr<float>(0);  
	total=0;
	for (i=0; i<hist.cols; i++)  //total是总的点数
    { 
		total+=hist_p[i]; 
	}
    for (i=0; i<threshold; i++)  
    {  
        total_back += hist_p[i];  
    }  
    total_object=total-total_back;  
  
    //背景熵;  
    for (i=0; i<threshold; i++)  
    {  
//      if(hist_p[i]==0)  
//          continue;  
        float percentage = hist_p[i]/total_back;  
		if(percentage >0)
		{
        entropy_back += -percentage * logf(percentage); // 能量的定义公式 
		}
    }  
    //前景熵;  
    for (i=threshold; i<hist.cols; i++)  
    {  
//      if(hist_p[i]==0)  
//      {  
//          continue;  
//      }  
        float percentage = hist_p[i]/total_object;  
		if(percentage >0)
		{
        entropy_object += -percentage * logf(percentage); // 能量的定义公式; 
		}
    }  
  
    entropy = entropy_object+entropy_back;  
	//entropy =entropy_object;  
    return entropy;  
}  

float LeftBackEntropy(const Mat& hist,int threshold)  //这个函数是测试随着阈值不断增加,左侧也就是背景熵的变化
{  
    float total_back=0,total_object=0;  
    float entropy_back=0,entropy_object=0;  
    float entropy = 0;  
    int i=0;  
    const float* hist_p = (float*) hist.ptr<float>(0);  
	total=0;
	for (i=0; i<hist.cols; i++)  //total是总的点数
    { 
		total+=hist_p[i]; 
	}
    for (i=0; i<threshold; i++)  
    {  
        total_back += hist_p[i];  
    }  
    total_object=total-total_back;  
  
    //背景熵;  
    for (i=0; i<threshold; i++)  
    {  
//      if(hist_p[i]==0)  
//          continue;  
        float percentage = hist_p[i]/total_back;  
		if(percentage >0)
		{
        entropy_back += -percentage * logf(percentage); // 能量的定义公式 
		}
    }   
  
    entropy = entropy_back;  
	//entropy =entropy_object;  
    return entropy;  
}  


  
void MaxEntropy(Mat img, Mat hist)  
{  
    total = sum(hist)[0];  
    float MaxEntropyValue = 0.0, MaxEntropyThreshold=0.0;  
    float tmp;  


	cout<<hist.size()<<endl;

    for (int i=0; i<hist.cols; i++)  
    {  
        tmp = calEntropy(hist,i); 
		
        if(tmp>MaxEntropyValue)  
        {  
            MaxEntropyValue = tmp;  
            MaxEntropyThreshold = i;  
        }  
    }  
    threshold(img,img,MaxEntropyThreshold,255,CV_THRESH_BINARY);  
    imshow("thresholdImg",img);  
    //imwrite("D:/thresholdImg.png",img);  
    cout<<MaxEntropyThreshold<<endl;  
    cout<<MaxEntropyValue<<endl;  
} 

void LeftEntropy(Mat img, Mat hist)   //这个函数是测试随着阈值不断增加,左侧也就是背景熵的变化
{
	
	total = sum(hist)[0];  
    float MaxEntropyValue = 0.0, MaxEntropyThreshold=0.0;  
    float tmp;  
	Mat SingleHist(hist.size(),CV_32FC1);//测试左边图像的熵值
	

    for (int i=0; i<hist.cols; i++)  
    {  
        tmp = LeftBackEntropy(hist,i); 		
		SingleHist.at<float>(0,i)= tmp;//这是测试左边图像的熵值

    }  
 
	Mat histImg = drawHist(SingleHist,256,200,Scalar(255,0,0));  
    imshow("SingleHist",histImg); 

}

int main(int argc,char *argv[])
{
  //Mat img = imread("smallpicture.jpg");
      Mat src = imread("smallpicture.jpg",0);  
    imshow("SRC",src);  
    Mat hist = Hist(src).t();
	MaxEntropy(src, hist);  
	LeftEntropy(src, hist);
/*	
	cout<<hist.size()<<endl;
	const float* hist_p = (float*) hist.ptr<float>(0);  

	int threshold=50;

	float total_back=0,total_object=0;  
    float entropy_back=0,entropy_object=0;  
    float entropy = 0;  
    int i=0;  
    //const float* hist_p = (float*) hist.ptr<float>(0);  
	total=0;
	for (i=0; i<hist.cols; i++)  //total是总的点数
    { 
		total+=hist_p[i]; 
	}
    for (i=0; i<threshold; i++)  
    {  
        total_back += hist_p[i];  
    }  
    total_object=total-total_back;  
  
    //背景熵;  
    for (i=0; i<threshold; i++)  
    {  
//      if(hist_p[i]==0)  
//          continue;  
        float percentage = hist_p[i]/total_back; 
		if(percentage>0)
		{
        entropy_back += -percentage * logf(percentage); // 能量的定义公式  
		}
    }  
    //前景熵;  
    for (i=threshold; i<hist.cols; i++)  
    {  
//      if(hist_p[i]==0)  
//      {  
//          continue;  
//      }  
        float percentage = hist_p[i]/total_object;  
		if(percentage>0)
		{
			entropy_object += -percentage * logf(percentage); // 能量的定义公式;  
		}
    }  
  
    entropy = entropy_object+entropy_back;  


	float percentage = hist_p[0]/total_back; 
	if(percentage)
	{
    entropy_back = -percentage * logf(percentage); // 能量的定义公式 
	cout<<entropy_back<<endl;
	}
	cout<<total_object<<"  "<<total_back<<"  "<<total<<endl;
	cout<<entropy<<endl;
	/*
	int i;
	for(i=0;i<50;i++)
	{
	cout<<hist_p[i]<< "  ";
	}
	for(i=0;i++;i<50)
	{
	cout<<endl;
	float sum=calEntropy(hist,50);
	cout<<sum<<"  "<<endl;
	}
	*/
	
   
    waitKey();  
    return 1;  

}
效果图如下:

3.MATLAB最大熵阈值分割代码

%一维最大熵法
E = imread('picture.bmp');
E = rgb2gray(E);
vHist=imhist(E);
[m,n]=size(E);
p=vHist(find(vHist>0))/(m*n);%求每一不为零的灰度值的概率
Pt=cumsum(p);%计算出选择不同t值时,A区域的概率
Ht=-cumsum(p.*log(p));%计算出选择不同t值时,A区域的熵
HL=-sum(p.*log(p));%计算出全图的熵
Yt=log(Pt.*(1-Pt)+eps)+Ht./(Pt+eps)+(HL-Ht)./(1-Pt+eps);%计算出选择不同t值时,判别函数的值
[a,th]=max(Yt);%th即为最佳阈值
segImg=(E>th);
figure,imshow(segImg)


与opencv得到的结果还是有些差别的,不知道为什么。


阅读更多
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

加入CSDN,享受更精准的内容推荐,与500万程序员共同成长!
关闭
关闭