OpenCV直方图与最大熵分割

原来一直觉得OpenCV里的直方图函数十分简单，今天临时需要用才发现原来OpenCV的calcHist功能如此强大，不仅能计算常见的1D Hist， calcHist理论上支持32维以下的Hist.(32维啊 有木有！)void calcHist(const Mat* images, int nimages, const int* channels, InputArray mask, OutputArray hist, int dims, const int* histSize, const float** ranges, bool uniform=true, bool accumulate=false )image: 输入图像序列nimages: 源图像数量channels: 用于计算hist的通道Listmask:不解释hist:生成的histdims:hist的维数（必须大于0，目前版本支持不大于CV_MAX_DIMS=32）histSzie:每一维的size（即bins）ranges: 直方图每一维的bins的边界的序列uniform:是否对齐accumlate:如果set为true，则直方图在使用之前不clear，用于在一个hist内保存多个图像集的直方图，或者及时更新hist。

实现代码：在此实现两个直方图绘制函数，分别绘制1D，2D直方图。
BGR三通道1D直方图：


附带的在程序里实现了一个窗口内显示多幅图像，也可以使用Mat1.put_back(Mat2)，但是push_back是向y轴方向push mat ，酱紫显示的就是细长细长窗口，所以还是使用ROI的用户体验比较好。

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;  
}  
  
void darwHistRGB(const Mat& src)  
{  
    Mat histB,histG,histR;  
  
    int bins=256;  
    int histSize[] = {bins};  
    float range[] = {0,256};  
    const float* ranges[] = {range};  
    int channelsB[] = {0};  
        int channelsG[] = {1};  
        int channelsR[] = {2};  
    calcHist(&src,1,channelsB,Mat(),histB,1,histSize,ranges,true,false);  
    calcHist(&src,1,channelsG,Mat(),histG,1,histSize,ranges,true,false);  
    calcHist(&src,1,channelsR,Mat(),histR,1,histSize,ranges,true,false);  
  
    Mat histBImg = drawHist(histB,bins,200,Scalar(255,0,0));  
    Mat histGImg = drawHist(histG,bins,200,Scalar(0,255,0));  
    Mat histRImg = drawHist(histR,bins,200,Scalar(0,0,255));  
      
    //在一个窗口中显示多幅图像  
    Mat display(200,bins*3,CV_8UC3);  
    Mat displayROI = display(Rect(0,0,bins,200));  
    resize(histBImg,displayROI,displayROI.size());  
    displayROI = display(Rect(bins,0,bins,200));  
    resize(histGImg,displayROI,displayROI.size());  
    displayROI = display(Rect(bins*2,0,bins,200));  
    resize(histRImg,displayROI,displayROI.size());  
  
    imshow("histRGB",display);  
    waitKey();  
}  
  
int main()  
{  
    Mat src = imread("D:/demo.jpg",1);  
    darwHistRGB(src);  
    return 1;  
}

HSV的H-S通道2D直方图：



灰度值的大小代表了直方图的高度。可以看做是一个从上向下（-Z轴方向）看的三维柱状图。

int main(  )  
{  
    Mat src, hsv;  
    src = imread("D:/demo.jpg", 1);  
  
  
    cvtColor(src, hsv, CV_BGR2HSV);  
  
    // Quantize the hue to 30 levels  
    // and the saturation to 32 levels  
    int hbins = 30, sbins = 32;  
    int histSize[] = {hbins, sbins};  
    // hue varies from 0 to 179, see cvtColor  
    float hranges[] = { 0, 180 };  
    // saturation varies from 0 (black-gray-white) to  
    // 255 (pure spectrum color)  
    float sranges[] = { 0, 256 };  
    const float* ranges[] = { hranges, sranges };  
    MatND hist;  
    // we compute the histogram from the 0-th and 1-st channels  
    int channels[] = {0, 1};  
  
    calcHist( &hsv, 1, channels, Mat(), // do not use mask  
        hist, 2, histSize, ranges,  
        true, // the histogram is uniform  
        false );  
    double maxVal=0;  
    minMaxLoc(hist, 0, &maxVal, 0, 0);  
  
    int scale = 10;  
    Mat histImg = Mat::zeros(sbins*scale, hbins*10, CV_8UC3);  
  
    for( int h = 0; h < hbins; h++ )  
        for( int s = 0; s < sbins; s++ )  
        {  
            float binVal = hist.at<float>(h, s);  
            int intensity = cvRound(binVal*255/maxVal);  
            rectangle( histImg, Point(h*scale, s*scale),  
                Point( (h+1)*scale - 1, (s+1)*scale - 1),  
                Scalar::all(intensity),  
                CV_FILLED );  
        }  
  
        namedWindow( "Source", 1 );  
        imshow( "Source", src );  
  
        namedWindow( "H-S Histogram", 1 );  
        imshow( "H-S Histogram", histImg );  
        waitKey();  
}

最大熵分割
利用Hist实现最大熵模型信息熵：

最大熵分割：

#include <iostream>  
#include <opencv/cv.h>  
#include <opencv/highgui.h>  
  
  
using namespace std;  
using namespace cv;  
  
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);  
    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;  
        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;  
        entropy_object += -percentage * logf(percentage); // 能量的定义公式；  
    }  
  
    entropy = entropy_object+entropy_back;  
    return entropy;  
}  
  
void MaxEntropy(Mat img, Mat hist)  
{  
    total = sum(hist)[0];  
    float MaxEntropyValue = 0.0, MaxEntropyThreshold=0.0;  
    float tmp;  
    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;  
}  
  
int main()  
{  
    Mat src = imread("D:/test1.jpg",0);  
    imshow("SRC",src);  
    Mat hist = Hist(src).t();  
    MaxEntropy(src, hist);  
    waitKey();  
    return 1;  
}