Learning Opencv 3 —— 十三章 直方图 Histogram

Opencv3 直方图 Histogram

直方图广泛应用于很多计算机视觉应用中，是计算机视觉中最经典的工具之一。简单的说，直方图就是对数据进行统计，将统计值组织到一系列事先定义好的 bin 中，bin 中的数值就是从数据中得到的特征统计量。边缘，色彩，角等直方图构成了可以被传递给目标识别分类器的一个通用特征类型。

直方图表示

由于直方图就是每个 bin 中记录的统计数值，因此直方图能够直接使用 Mat，vector<> 以及 sparse mat进行表示。事实上，就是由于高维直方图绝大多数的数值都为零，为更好地对其进行描述，sparse mat 才被引入。

直方图计算：calcHist

void cv::calcHist(
	// 源输入（图像）数组，必须是相同深度的CV_8U或者CV_32F(即uchar或者float),相同大小，每一个可以是任意通道的;
	const cv::Mat* images, 
	// 源输入数组中的元素个数;
	int nimages, 
	// 用来计算直方图的通道维数数组,第一个数组的通道由0到arrays[0].channels()-1列出，
	// 第二个数组的通道从arrays[0].channels()到arrays[0].channels()+arrays[1].channels()-1以此类推;
	const int* channels, 
	// 可选的掩膜，如果该矩阵不是空的，则必须是8位的并且与arrays[i]的大小相等，
	// 掩膜的非零值标记需要在直方图中统计的数组元素;
	cv::InputArray mask, 
	// 输出直方图，是一个稠密或者稀疏的dims维的数组;
	cv::OutputArray hist, 
	// 直方图的维数，必须为正，并且不大于CV_MAX_DIMS(当前的OpenCV版本中为32，即最大可以统计32维的直方图);
	int dims, 
	// 用于指出直方图数组每一维的大小的数组，即指出每一维的bin的个数的数组;
	const int* histSize, 
	// 用于指出直方图每一维的每个bin的上下界范围数组的数组，
	// 当直方图是均匀的(uniform =true)时，对每一维i指定直方图的第0个bin的下界(包含即[)L0
	// 和最后一个即第histSize[i]-1个bin的上界(不包含的即))U_histSize[i]-1，
	// 也就是说对均匀直方图来说，每一个ranges[i]都是一个两个元素的数组【指出该维的上下界】。
	// 当直方图不是均匀的时，每一个ranges[i]数组都包含histSize[i]+1个元素:
	// L0,U0=L1,U1=L1,...,U_histSize[i]-2 = L_histSize[i] - 1, U_histSize[i] - 1.
	// 不在L0到U_histSize[i] - 1之间的数组元素将不会统计进直方图中;
	const float** ranges,
	// 直方图是否均匀的标志;【指定直方图每个bin统计的是否是相同数量的灰度级】
	bool uniform = true, 
	// 累加标志；为true，将在现有统计数据基础上继续计算
	bool accumulate = false 
);

opencv文件中的范例代码

#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>

using namespace cv;

int main(int argc, char** argv)
{
	Mat src, hsv;
	if (argc != 2 || !(src = imread(argv[1], 1)).data)
		return -1;

	cvtColor(src, hsv, COLOR_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();
}

直方图其他操作

直方图归一化

cv::Mat normalized = my_hist / sum(my_hist)[0];

// Or
cv::normalize( my_hist, my_hist, 1, 0, NORM_L1 );

直方图取阈值

cv::threshold(
	my_hist, // input histogram
	my_thresholded_hist, // result, all values<threshold set to zero
	threshold, // cutoff value
	0, // value does not matter in this case
	cv::THRESH_TOZERO // threshold type
);

直方图取最大最小值

void cv::minMaxLoc(
	cv::InputArray src, // Input array
	double* minVal, // put minimum value (if not NULL)
	double* maxVal = 0, // put maximum value (if not NULL)
	cv::Point* minLoc = 0, // put minimum location (if not NULL)
	cv::Point* maxLoc = 0, // put maximum location (if not NULL)
	cv::InputArray mask = cv::noArray() // ignore points for which mask is zero
);

double max_val;
cv::Point max_pt; // Mat 水平向右为x轴正向，垂直向下为y轴正向
cv::minMaxLoc(
	my_hist, // input histogram
	NULL, // don't care about the min value
	&max_val, // place to put the maximum value
	NULL, // don't care about the location of the min value
	&max_pt // place to put the maximum value location (a cv::Point)
);

对于稀疏高维直方图，可以使用整数向量来保存位置

double maxval;
int max_pt[CV_MAX_DIM];
cv::minMaxLoc(
	my_hist, // input sparse histogram
	NULL, // don't care about the min value
	&max_val, // place to put the maximum value
	NULL, // don't care about the location of the min value
	max_pt // place to put the maximum value location (a cv::Point)
);

对于常规高维直方图，使用 minMaxIdx。需要注意的是，如果直方图是一维的，位置数组也应该是二维的。

void cv::minMaxIdx(
	cv::InputArray src,
	double* minVal, // put min value (if not NULL)
	double* maxVal = 0, // put max value (if not NULL)
	int* minLoc = 0, // put min location indices (if not NULL)
	int* maxLoc = 0, // put max location indices (if not NULL)
	cv::InputArray mask = cv::noArray() // ignore points if mask is zero
);

直方图对比

double cv::compareHist(
	cv::InputArray H1, // First histogram to be compared
	cv::InputArray H2, // Second histogram to be compared
	int method // comparison method (see options below)
);

对比方法

1. Correlation method (cv::COMP_CORREL)，最佳匹配为 1，最不匹配为 -1，无关为 0
2. Chi-square method (cv::COMP_CHISQR_ALT)，最大匹配为 0，越不匹配值越大
3. Intersection method (cv::COMP_INTERSECT)，越匹配值越大，归一化数据最佳匹配为 1，最不匹配为 0
4. Bhattacharyya distance method (cv::COMP_BHATTACHARYYA)，最佳匹配为 0，最不匹配为 1
5. Earth Mover’s Distance，其为独立的函数，不能在 compareHist 中作为方法使用

最后，虽然直方图对比允许使用未归一化的数据，但是对于比如 Intersection method 的最后结果还是会有一定的影响。因此，在进行直方图对比之前，最后对数据进行归一化。速度和性能方便，3 最快，性能也最差；2 和 4 慢一些但性能更好。5 最慢但效果最佳。

陆地移动距离 Earth Mover’s Distance

之前例举直方图对比方法，都是基于每个 bin 中的数值进行计算，因此如果两个直方图形状相似，只是相对平移了一段距离，那么之前的方法都将得出不匹配的结果。而 EMD 本质上度量了怎么将一个直方图的形状转变为另一个直方图的形状。同时，EMD 允许自定义距离度量标准和移动代价。

float cv::EMD(
	cv::InputArray signature1, // sz1-by-(dims+1) float array
	cv::InputArray signature2, // sz2-by-(dims+1) float array
	int distType, // distance type (e.g., 'cv::DIST_L1')
	cv::InputArray cost = noArray(), // sz1-by-sz2 array (if cv::DIST_USER)
	float* lowerBound = 0, // input/output low bound on distance
	cv::OutputArray flow = noArray() // output, sz1-by-sz2
);

• signature：对于三维直方图，如果在（7，43，11）位置上数量为537，则对应数据为 [ 537, 7, 43, 11 ]
• distType：汉密尔顿距离 (cv::DIST_L1), 欧式距离 (cv::DIST_L2), 棋盘距离 (cv::DIST_C), 或者用户自定义 (cv::DIST_USER)。
• cost：如果为用户自定义距离，在此指定
• lowerBound：只有当不使用自定义距离，直方图权重相同，且 signature 不只包含权重（有一列以上）时，才被使用。同时，质心（mass center）的距离大于或等于 lowerBound 时，将不计算 EMD 而只返回 lowerBound。
• flow：flow(i,j) = 从 signature1 的第 i 个点到 signature2 的第 j 个点的流

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;

void help(char** argv) {
	cout << "\nCall is:\n"
		<< argv[0] << " modelImage0 testImage1 testImage2 badImage3\n\n"
		<< "for example: "
		<< " ./ch7_ex7_3_expanded HandIndoorColor.jpg HandOutdoorColor.jpg "
		<< "HandOutdoorSunColor.jpg fruits.jpg\n"
		<< "\n";
}

// Compare 3 images' histograms
int main(int argc, char** argv) {
	if (argc != 5) { help(argv); return -1; }
	vector<cv::Mat> src(5);
	cv::Mat tmp;
	int i;
	tmp = cv::imread(argv[1], 1);
	if (tmp.empty()) {
		cerr << "Error on reading image 1," << argv[1] << "\n" << endl;
		help(argv);
		return(-1);
	}
	// Parse the first image into two image halves divided halfway on y
	//
	cv::Size size = tmp.size();
	int width = size.width;
	int height = size.height;
	int halfheight = height >> 1;
	cout << "Getting size [[" << tmp.cols << "] [" << tmp.rows << "]]\n" << endl;
	cout << "Got size (w,h): (" << size.width << "," << size.height << ")" << endl;
	src[0] = cv::Mat(cv::Size(width, halfheight), CV_8UC3);
	src[1] = cv::Mat(cv::Size(width, halfheight), CV_8UC3);
	// Divide the first image into top and bottom halves into src[0] and src[1]
	//
	cv::Mat_<cv::Vec3b>::iterator tmpit = tmp.begin<cv::Vec3b>();
	// top half
	//
	cv::Mat_<cv::Vec3b>::iterator s0it = src[0].begin<cv::Vec3b>();
	for (i = 0; i < width*halfheight; ++i, ++tmpit, ++s0it) *s0it = *tmpit;
	// Bottom half
	//
	cv::Mat_<cv::Vec3b>::iterator s1it = src[1].begin<cv::Vec3b>();
	for (i = 0; i < width*halfheight; ++i, ++tmpit, ++s1it) *s1it = *tmpit;
	// Load the other three images
	//
	for (i = 2; i<5; ++i) {
		src[i] = cv::imread(argv[i], 1);
		if (src[i].empty()) {
			cerr << "Error on reading image " << i << ": " << argv[i] << "\n" << endl;
			help(argv);
			return(-1);
		}
	}
	// Compute the HSV image, and decompose it into separate planes.
	//
	vector<cv::Mat> hsv(5), hist(5), hist_img(5);
	int h_bins = 8;
	int s_bins = 8;
	int hist_size[] = { h_bins, s_bins }, ch[] = { 0, 1 };
	float h_ranges[] = { 0, 180 }; // hue range is [0,180]
	float s_ranges[] = { 0, 255 };
	const float* ranges[] = { h_ranges, s_ranges };
	int scale = 10;
	for (i = 0; i<5; ++i) {
		cv::cvtColor(src[i], hsv[i], cv::COLOR_BGR2HSV);
		cv::calcHist(&hsv[i], 1, ch, cv::Mat(), hist[i], 2, hist_size, ranges, true);
		cv::normalize(hist[i], hist[i], 0, 255, cv::NORM_MINMAX);
		hist_img[i] = cv::Mat::zeros(hist_size[0] * scale, hist_size[1] * scale, CV_8UC3);
		// Draw our histogram For the 5 images
		//
		for (int h = 0; h < hist_size[0]; h++)
			for (int s = 0; s < hist_size[1]; s++) {
				float hval = hist[i].at<float>(h, s);
				cv::rectangle(
					hist_img[i],
					cv::Rect(h*scale, s*scale, scale, scale),
					cv::Scalar::all(hval),
					-1
				);
			}
	}
	// Display
	//
	cv::namedWindow("Source0", 1); cv::imshow("Source0", src[0]);
	cv::namedWindow("HS Histogram0", 1); cv::imshow("HS Histogram0", hist_img[0]);
	cv::namedWindow("Source1", 1); cv::imshow("Source1", src[1]);
	cv::namedWindow("HS Histogram1", 1); cv::imshow("HS Histogram1", hist_img[1]);
	cv::namedWindow("Source2", 1); cv::imshow("Source2", src[2]);
	cv::namedWindow("HS Histogram2", 1); cv::imshow("HS Histogram2", hist_img[2]);
	cv::namedWindow("Source3", 1); cv::imshow("Source3", src[3]);
	cv::namedWindow("HS Histogram3", 1); cv::imshow("HS Histogram3", hist_img[3]);
	cv::namedWindow("Source4", 1); cv::imshow("Source4", src[4]);
	cv::namedWindow("HS Histogram4", 1); cv::imshow("HS Histogram4", hist_img[4]);
	// Compare the histogram src0 vs 1, vs 2, vs 3, vs 4
	cout << "Comparison:\n"
		<< "Corr Chi Intersect Bhat\n"
		<< endl;
	for (i = 1; i<5; ++i) { // For each histogram
		cout << "Hist[0] vs Hist[" << i << "]: " << endl;;
		for (int j = 0; j<4; ++j) { // For each comparison type
			cout << "method[" << j << "]: " << cv::compareHist(hist[0], hist[i], j) << " ";
		}
		cout << endl;
	}
	//Do EMD and report
	//
	vector<cv::Mat> sig(5);
	cout << "\nEMD: " << endl;
	// Oi Vey, parse histograms to earth movers signatures
	//
	for (i = 0; i<5; ++i) {
		vector<cv::Vec3f> sigv;
		// (re)normalize histogram to make the bin weights sum to 1.
		//
		cv::normalize(hist[i], hist[i], 1, 0, cv::NORM_L1);
		for (int h = 0; h < h_bins; h++)
			for (int s = 0; s < s_bins; s++) {
				float bin_val = hist[i].at<float>(h, s);
				if (bin_val != 0)
					sigv.push_back(cv::Vec3f(bin_val, (float)h, (float)s));
			}
		// make Nx3 32fC1 matrix, where N is the number of nonzero histogram bins
		//
		sig[i] = cv::Mat(sigv).clone().reshape(1);
		if (i > 0)
			cout << "Hist[0] vs Hist[" << i << "]: "
			<< EMD(sig[0], sig[i], cv::DIST_L2) << endl;
	}
	cv::waitKey(0);
}

反向投影 Back Projection

反向投影度量了像素块与已知直方图的匹配性，其能够被用来从一副图像中查找感兴趣的区域（与模式直方图匹配）。从概率统计的角度，反向投影可以理解为计算图像中各像素来自给定分布（直方图）的概率。

void cv::calcBackProject(
	const cv::Mat* images, // C-style array of images, 8U or 32F
	int nimages, // number of images in 'images' array
	const int* channels, // C-style list, ints identifying channels
	cv::InputArray hist, // input histogram array
	cv::OutputArray backProject, // output single channel array
	const float** ranges, // C-style array, 'dims' pairs of bin sizes
	double scale = 1, // Optional scale factor for output
	bool uniform = true // true for uniform binning
);

其参数意义与 calcHist 基本一致。