【OpenCV】对比度增强之直方图均衡化(全局)

直方图均衡化属于数字图像处理中灰度变换(intensity transformation)的内容，灰度变换的目的就是找到一个合适的映射函数s=T$(r)$(r).将原图像的灰度值映射到新的图像中，已达到优化图像的目的。
直方图均衡化是通过调整图像的灰阶分布，使得在0~255灰阶上的分布更加均衡，提高了图像的对比度，达到改善图像主观视觉效果的目的。对比度较低的图像适合使用直方图均衡化方法来增强图像细节。
直方图均衡数学背景是将一个分布（强度值给定的直方图）映射到另一个分布（强度值更宽和理想的均匀分布）。也就是说，我们希望在新分配中尽可能均匀分布原始分布的y值。事实证明，解决扩展分布值的问题的一个好方法是：重映射函数应该是累积分布函数。

公式的连续化
假设原图像的灰度统计直方图标准化后为$P_r(r)$Pr​(r).原图像灰度范围为(0~L-1).那么直方图均衡化找到的就是这样一个映射函数：
S=(L-1)$\int^{r}_0p_{r}(w){\rm d}x$∫0r​pr​(w)dx
设映射后的图像的灰度分布为$p_s(s)$ps​(s),再由概率论相关理论（随机变量函数的概率密度与随机变量概率密度的关系）可知：
$p_s(s)$ps​(s)=$p_r(r)$pr​(r)|$\frac{{\rm d}r}{{\rm d}s}$dsdr​|
对映射函数两边进行求导
$\frac{{\rm d}s}{{\rm d}r}$drds​=(L-1)$p_r(r)$pr​(r)
所以我们可以得到变换后的图像直方图分布为
$p_s(s)$ps​(s)=$p_r(r)$pr​(r)|$\frac{1}{\frac{{\rm d}s}{{\rm d}r}}$drds​1​|=$p_r(r)\frac{1}{(L-1)p_r}$pr​(r)(L−1)pr​1​=$\frac{1}{L-1}$L−11​
我们可以看到变换后的图像灰度直方图分布恒为1/(L-1),这就达到了上面的目的,使得图像的灰度分布更均匀,层次感更强.
注:在灰度变换中,变化函数T®需要满足下面的两点要求，
1.当0$\leq$≤r$\leq$≤L-1时, $T(r)$T(r) 是一个严格递增函数。
2.当0$\leq$≤r$\leq$≤L-1时，0$\leq$≤$T(r)$T(r)$\leq$≤L-1
第一点要求的原因是，对于变化前像素和变换后像素灰度的明暗顺序不能改变，之所以要严格递增，是为了确保变化前和变换后像素可以一一对应。
第二点要求的原因是，变换后的图像不能超过原先的灰度级数。
不难发现，其实直方图均衡化的过程并不一定满足条件1。所以该变换式不可逆的。

公式的离散化
设原图像灰度等级为0、1、2…L-1,离散化后的映射公式就是
$s(r)$s(r)=(L-1)$\sum\limits_{i=0}^rp_r(i)$i=0∑r​pr​(i)
在利用上面的公式进行计算的时候，需要把计算的结果$s(r)$s(r),近似为最近的整数。
1)计算输入图像的直方图；
2）进行直方图归一化，直方图的组距和为255；
3）计算直方图积分：
$\qquad$ $H^{'}(i)$H′(i)= $\sum_{0 $\leq$≤ j $\leq$≤i }$H(j)$H(j)$
4）以$H^{'}$H′作为查询表进行图像变换：
$\qquad$ dst(x,y)=$H^{'}(src(x,y))$H′(src(x,y))
简而言之，由equalizeHist()函数实现的灰度直方图均衡化算法，就是把直方图的每个灰度级进行归一化处理，求每种灰度的累积分布，得到一个映射的灰度映射表，然后根据相应的灰度值来修正原图中的每个像素。

原理详解：
假设输入图像为I，高为H、宽为W, $hist_I$histI​代表灰度直方图，$hist_I(k)$histI​(k)代表灰度值等于k的像素点个数，其中k$\epsilon$ϵ[0,255]。全局直方图均衡化操作是对图像I进行改变，使得输出图像O的灰度直方图$hist_O$histO​是"平“的，即每一个灰度级的像素点个数是”相等“的。注意，其实这里的”相等“不是严格意义上的等于，而是约等于，比如高为137、宽为255的图像矩阵不可能出现每一个灰度级的像素点个数是严格相等的，即$hist_O(k)$histO​(k)$\approx$≈$\frac{H*W}{256}$256H∗W​,k$\epsilon$ϵ[0,255],那么对于任意的灰度级p,0$\leq$≤q$\leq$≤ 255,总能找到q,0$\leq$≤q$\leq$≤ 255,使得：

$\sum\limits_{k=0}^p{hist_I(k)}$k=0∑p​histI​(k)=$\sum\limits_{k=0}^q{hist_O(k)}$k=0∑q​histO​(k)
其中$\sum\limits_{k=0}^p{hist_I(k)}$k=0∑p​histI​(k)和$\sum\limits_{k=0}^q{hist_O(k)}$k=0∑q​histO​(k)称为I和O的累加直方图。又因为$hist_O(k)$histO​(k)$\approx$≈$\frac{H*W}{256}$256H∗W​,所以

$\sum\limits_{k=0}^p{hist_I(k)}$k=0∑p​histI​(k)$\approx$≈(q+1)$\frac{H*W}{256}$256H∗W​

化简上式可得
q$\approx$≈ $\frac{\sum\limits_{k=0}^{p}{hist_I(k)}}{H*W}$H∗Wk=0∑p​histI​(k)​*256-1

上式给出了一个从亮度级为p的输入像素到亮度级为q的输出像素的映射，那么令

O(r,c)$\approx$≈ $\frac{\sum\limits_{k=0}^{I(r,c)}{hist_I(k)}}{H*W}$H∗Wk=0∑I(r,c)​histI​(k)​*256-1

其中I(r,c)是I的第r行第c列的灰度值，O(r,c)是对应位置输出的灰度值，其中0$\leq$≤r<H,0$\leq$≤c<W,这样就计算出了输出图像O的每一个位置的灰度值。

C++实现
对于直方图均衡化的C++实现，通过定义函数equalHist依次按照四个步骤来实现，输入参数为8位的灰度图，代码如下

Mat equalHist(Mat image)
{
	CV_Assert(image.type()==CV_8UC1);
	//灰度图像的高、宽
	int rows = image.rows;
	int cols = image.cols;
	//第一步：计算图像的灰度直方图
	Mat grayHist = calcGrayHist(image);
	//第二步：计算累加灰度直方图
	Mat zeroCumuMoment = Mat::zeros(Size(256,1),CV_32SC1);
	for(int p =0;p<256;p++)
	{
		if(p==0)
			zeroCumuMoment.at<int>(0,p)=grayHist.at<int>(0,0);
		else
			zeroCumuMoment.at<int>(0,p)=zeroCumuMoment.at<int>(0,p-1)+grayHist.at<int>(0,p);
	}
	// 第三步，根据累加直方图得到输入灰度级和输出灰度级之间的映射关系
	Mat outPut_q=Mat::zeros(Size(256,1),CV_8UC1);
	float cofficient = 256.0/(rows*cols);
	for(int p=0; p<256;p++)
	{
		float q = cofficient*zeroCumuMoment.at<int>(0,p)-1;
		if(q>=0)
			outPut_q.at<uchar>(0,p)=uchar(floor(q));
		else
			outPut_q.at<uchar>(0,p)=0;
	}
	//第四步：得到直方图均衡化后的图像
	Mat equalHistImage=Mat::zeros(image.size(),CV_8UC1);
	for(int r =0;r<rows;r++)
	{
		for(int c=0;c>cols;c++)
		{
			int p = image.at<uchar>(r,c);
			equalHistImage.at<uchar>(r,c)=outPut_q.at<uchar>(0,p);
		}
	}

}