《数字图像处理》直方图匹配（规定化）- C++源码

原理见《数字图像处理》第三版3.3.2

下面直接上源码（VS2013+OpenCV3.0）

/*
//***************************************************************************************************************************************
//-----------------------------------------------作者：我三食堂不服-------------------------------------------------------------------
//***************************************************************************************************************************************
*/


#include <opencv2/opencv.hpp>  
#include <math.h>
#include <vector>

using namespace std;
using namespace cv;

//计算直方图
int CalHist(const Mat &input, vector<long> &histogram);

//绘制直方图，dep - 位深度
Mat DrawHist(const vector<long> &hist);

//获取最大值
template<class T> T Max(const vector<T> &inputArray);

//直方图规定化，z - 规定的各灰度值的概率
int MatchHist(const Mat &input, double *z, Mat &out);

int main(int argc, char** argv)
{
	Mat im_source1;				//原图
	Mat im_pro1;				//处理后的图片
	Mat im_Hist1;				//灰度直方图
	vector<long> histogram1;	//灰度直方图数据
	double *z;					//规定的直方图概率分布

	z = new double[256]{};

	//输入概率分布规定
	for (int i = 0; i < 5; i++)
		z[i] = i * 7 / 4;
	for (int i = 5; i < 21; i++)
		z[i] = 7 - (i - 4)*6.3 / 16;
	for (int i = 21; i < 180; i++)
		z[i] = 0.7 - (i - 20)*0.7 / 159;
	for (int i = 180; i < 201; i++)
		z[i] = (i - 179)*0.48 / 21;
	for (int i = 201; i < 256; i++)
		z[i] = 0.48 - (i - 200)*0.48 / 55;
	
	//载入原图（灰度图方式）
	im_source1 = imread(\"3-23_a.tif\", IMREAD_GRAYSCALE);

	//图片处理

	MatchHist(im_source1, z, im_pro1);
	CalHist(im_pro1, histogram1);
	im_Hist1 = DrawHist(histogram1);
	//创建窗口
	namedWindow(\"【原图1】\");
	namedWindow(\"【处理后b1】\");
	namedWindow(\"【处理后c1】\");
	//显示图片
	imshow(\"【原图1】\", im_source1);
	imshow(\"【处理后b1】\", im_Hist1);
	imshow(\"【处理后c1】\", im_pro1);

	//保存图片
	imwrite(\"3-25_d.tif\", im_Hist1);
	imwrite(\"3-25_c.tif\", im_pro1);
	
	//等待键盘操作
	cvWaitKey(0);
	
	return 0;
}


//直方图规定化，z - 规定的各灰度值的概率
int MatchHist(const Mat &input, double *z, Mat &out)
{
	int h{ input.rows };				// 图片高
	int w{ input.cols };				// 图片宽
	uchar* data_in{ input.data };		// 原图图片像素数据
	uchar* data_out;					// 输出图片像素数据
	int data;							// 灰度值
	int index;							// 映射索引
	int index1;							// 映射索引
	int index2;							// 映射索引
	double *s;							// 输入图像r到均衡化s的灰度映射s[r]
	double *g;							// 规定的概率分布z到均衡化s\'的映射g[z]
	int *_g;							// 反映射，即z均衡化后s\'到z的映射_g[s\']
	vector<long> hist;					// 灰度直方图数据
	int bit_dep;						// 位比特数
	double k{ 0 };						// 概率分布修正系数

	//输出图像初始化
	out.release();
	out.create(h, w, input.type());
	data_out = out.data;

	//计算直方图
	CalHist(input, hist);

	//获取位深度
	bit_dep = hist.size();

	s = new double[bit_dep]{};
	g = new double[bit_dep]{};
	_g = new int[bit_dep]{};
	for (int i = 0; i < bit_dep; i++)
	{
		s[i] = 0;
		g[i] = 0;
		_g[i] = -1;
	}

	//修正规定的概率分布
	for (int i = 0; i < bit_dep; i++)
		k += z[i];
	for (int i = 0; i < bit_dep; i++)
		z[i] /= k;

	//计算r到s的映射s[r]
	for (int i = 0; i < bit_dep; i++)
	for (int j = 0; j <= i; j++)
	{
		s[i] += (double)(bit_dep - 1) / (h*w)*hist[j];
	}

	//计算规定的概率分布z到均衡化s\'的映射g[z]，及其反映射_g[s\']
	for (int i = 0; i < bit_dep; i++)
	{
		for (int j = 0; j <= i; j++)
			g[i] += (bit_dep - 1) * z[j];

		if (_g[cvRound(g[i])]>i || _g[cvRound(g[i])] < 0)
			_g[cvRound(g[i])] = i;
	}

	//映射灰度值
	for (int y = 0; y < h; y++)
	for (int x = 0; x < w; x++)
	{
		index = cvRound(s[data_in[y*w + x]]);
		data = _g[index];

		if (data < 0)
		for (int i = 1;;i++)
		{
			index1 = index - i;
			index2 = index + i;
			if (index1 >= 0 && index1 < bit_dep)
			if (_g[index1] >= 0)
			{
				data = _g[index1];
				break;
			}

			if (index2 >= 0 && index2 < bit_dep)
			if (_g[index2] >= 0)
			{
				data = _g[index2];
				break;
			}
		}
		
		data_out[y*w + x] = (uchar)data;
	}

	delete s;
	s = nullptr;
	delete g;
	g = nullptr;
	delete _g;
	_g = nullptr;

	return 1;
}


//计算直方图
int CalHist(const Mat &input, vector<long> &histogram)
{
	int h{ input.rows };				// 图片高
	int w{ input.cols };				// 图片宽
	uchar* data_in{ input.data };		// 原图图片像素数据
	long *hist;							// 直方图数据
	int bit_dep;						// 位比特数

	histogram.clear();

	// 获取比特数
	switch (input.depth())
	{
	case CV_8U:
		bit_dep = 256;
		break;
	case CV_16U:
		bit_dep = 65536;
		break;
	default:
		break;
	}

	hist = new long[bit_dep]{};

	//计算直方图
	for (int y = 0; y < h; y++)
	for (int x = 0; x < w; x++)
	{
		hist[data_in[y*w + x]]++;
	}

	//输出
	for (int i = 0; i < bit_dep; i++)
	{
		histogram.push_back(hist[i]);
	}
	
	delete hist;
	hist = nullptr;
	return 1;
}

//绘制直方图，dep - 位深度
Mat DrawHist(const vector<long> &hist)
{

	double kh, kw;			//直方图横、纵坐标放大系数
	kh = 590.0 / Max(hist);
	kw = 1000.0 / hist.size();

	Mat out = Mat::zeros(600, 1000, CV_8UC1);

	for (int i = 0; i < (int)hist.size(); i++)
	{
		line(out, Point((int)(i*kw), (int)(out.rows - hist[i] * kh)), Point((int)(i*kw), out.rows), Scalar(255), (int)(0.7*kw));
	}

	return out;
}

//获取最大值
template<class T> T Max(const vector<T> &inputArray)
{
	T re{ inputArray[0] };
	for (size_t i = 0; i < inputArray.size(); i++)
	{
		if (re < inputArray[i])
			re = inputArray[i];
	}
	
	return re;
}

效果如下：