學習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的通道List

mask:不解釋

hist:生成的hist

dims: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;
}