循序渐进之（五）空间域图像增强之自适应直方图均衡化（AHE）

                      循序渐进之（五）空间域图像增强之自适应直方图均衡化（AHE）

文字摘自：对比度受限的自适应直方图均衡化(CLAHE)

直方图均衡化（HE）是一种很常用的直方图类方法，基本思想是通过图像的灰度分布直方图确定一条映射曲线，用来对图像进行灰度变换，以达到提高图像 对比度的目的。该映射曲线其实就是图像的累计分布直方图（CDF）（严格来说是呈正比例关系）。然而HE是对图像全局进行调整的方法，不能有效地提高局部 对比度，而且某些场合效果会非常差。如：

上述原图和HE结果图的直方图分别为：

因为从原图的直方图中求取的映射函数（CDF）形状为：

将它作用于原图像会导致直方图被整体右移，没有充分利用整个灰度动态范围。

为了提高图像的局部对比度，有人提出将图像分成若干子块，对子块进行HE处理，这便是AHE（自适应直方图均衡化），使用AHE处理上图得到：

结果直方图：

可 以看出结果图像的灰度较好地分布在了全部动态范围上。从结果图像上也可以看出，局部对比度的确得到了提高，视觉效果要优于HE。但是仍然有个问题：AHE 对局部对比度提高过大，导致图像失真。看看背景区，本来的黑色背景现在已经变成白色了，原因是因为背景区中的局部子块统计得到的直方图在0灰度处幅值太高 （实际上全黑子图基本上就集中在0灰度处），这样导致映射曲线斜率过高，将所有灰度值都映射到整个灰度轴的右侧，所以结果图中背景偏白。

（上边的AHE算法实现的图像处理效果为添加了线性插值功能，这里代码实现的AHE,暂时没有加入线性插值功能，所以可以看出图像分块）

代码参考：CLAHE的实现和研究   中的AHE算法部分。

#include<iostream>
#include<opencv.hpp>
using namespace std;
using namespace cv;
int main()
{
	Mat originImage = imread("E://匹配.jpg", 0);
	imshow("原图", originImage);
	Mat src = originImage.clone();
	const int blockNumber = 8;//把图像分成block数量
	int width = src.cols;
	int height = src.rows;
	int singleBlockWidth = src.cols / 8;//每个block大小
	int singleBlockHeight = src.rows / 8;
	int pixelNumber[blockNumber *blockNumber][256] = { 0 };//存储不同block中各个不同灰度像素数量
	float total[blockNumber *blockNumber][256] = { 0.0 };//累计直方图
	for (int i = 0; i < blockNumber; i++)
	{
		for (int j = 0; j < blockNumber; j++)
		{
			int startPixelW = (i)*singleBlockWidth;
			int endPixelW = (i+1)*singleBlockWidth;
			int startPixelH = (j)*singleBlockHeight;
			int endPixelH = (j+1)*singleBlockHeight;
			int number = i + 8 * j;//统计运算到哪一个block了
			int singleBlockPixelNumber = singleBlockWidth*singleBlockHeight;
			for (int x = startPixelW; x < endPixelW; x++)//统计不同block中各个不同灰度像素数量
				for (int y = startPixelH; y < endPixelH; y++)
				{
					int pixelValue = src.at<uchar>(y, x);
					pixelNumber[number][pixelValue]++;
				}
			for (int k = 0; k < 256; k++)//计算累计直方图
			{
				if (k == 0)
					total[number][k] = 1.0*pixelNumber[number][k] / singleBlockPixelNumber;
				else
					total[number][k] = total[number][k - 1] + 1.0*pixelNumber[number][k] / singleBlockPixelNumber;
			}
		}
	}
	for (int i = 0; i < blockNumber; i++)//利用累计直方图对于原像素灰度在各自block中进行映射
	{
		for (int j = 0; j < blockNumber; j++)
		{
			int startPixelW = (i)*singleBlockWidth;
			int endPixelW = (i+1)*singleBlockWidth;
			int startPixelH = (j)*singleBlockHeight;
			int endPixelH = (j+1)*singleBlockHeight;
			int number = i + 8 * j;
			int singleBlockPixelNumber = singleBlockWidth*singleBlockHeight;
			for (int x = startPixelW; x < endPixelW; x++)
				for (int y = startPixelH; y < endPixelH; y++)
				{
					int pixelValue = src.at<uchar>(y, x);
					src.at<uchar>(y, x) = total[number][pixelValue] * 255;
				}
		}
	}
	imshow("均衡图", src);
	waitKey(0);
	return 0;
}

效果如下：

加入线性插值功能，代码参考：CLAHE的实现和研究 中的CLAHE算法部分（这里的线性插值代码好像较为复杂，后续看有没有更简便的实现手段）。。

#include<iostream>
#include<opencv.hpp>
using namespace std;
using namespace cv;
int main()
{
	Mat originImage = imread("E://匹配.jpg", 0);
	imshow("原图", originImage);
	Mat src = originImage.clone();
	Mat src1 = originImage.clone();
	const int blockNumber = 8;//把图像分成block数量
	int width = src.cols;
	int height = src.rows;
	int singleBlockWidth = src.cols / 8;//每个block大小
	int singleBlockHeight = src.rows / 8;
	int pixelNumber[blockNumber *blockNumber][256] = { 0 };//存储不同block中各个不同灰度像素数量
	float total[blockNumber *blockNumber][256] = { 0.0 };//累计直方图
	for (int i = 0; i < blockNumber; i++)
	{
		for (int j = 0; j < blockNumber; j++)
		{
			int startPixelW = (i)*singleBlockWidth;
			int endPixelW = (i + 1)*singleBlockWidth;
			int startPixelH = (j)*singleBlockHeight;
			int endPixelH = (j + 1)*singleBlockHeight;
			int number = i + 8 * j;//统计运算到哪一个block了
			int singleBlockPixelNumber = singleBlockWidth*singleBlockHeight;
			for (int x = startPixelW; x < endPixelW; x++)//统计不同block中各个不同灰度像素数量
				for (int y = startPixelH; y < endPixelH; y++)
				{
					int pixelValue = src.at<uchar>(y, x);
					pixelNumber[number][pixelValue]++;
				}
			for (int k = 0; k < 256; k++)//计算累计直方图
			{
				if (k == 0)
					total[number][k] = 1.0*pixelNumber[number][k] / singleBlockPixelNumber;
				else
					total[number][k] = total[number][k - 1] + 1.0*pixelNumber[number][k] / singleBlockPixelNumber;
			}
		}
	}
	for (int i = 0; i < blockNumber; i++)//利用累计直方图对于原像素灰度在各自block中进行映射
	{
		for (int j = 0; j < blockNumber; j++)
		{
			int startPixelW = (i)*singleBlockWidth;
			int endPixelW = (i + 1)*singleBlockWidth;
			int startPixelH = (j)*singleBlockHeight;
			int endPixelH = (j + 1)*singleBlockHeight;
			int number = i + 8 * j;
			int singleBlockPixelNumber = singleBlockWidth*singleBlockHeight;
			for (int x = startPixelW; x < endPixelW; x++)
				for (int y = startPixelH; y < endPixelH; y++)
				{
					int pixelValue = src1.at<uchar>(y, x);
					src1.at<uchar>(y, x) = total[number][pixelValue] * 255;
				}
		}
	}
	imshow("均衡图无线性差值", src1);
	for (int i = 0; i < width; i++)
	{
		for (int j = 0; j < height; j++)
		{
			//four coners  
			if (i <= singleBlockWidth / 2 && j <= singleBlockHeight / 2)
			{
				int num = 0;
				src.at<uchar>(j, i) = (int)(total[num][src.at<uchar>(j, i)] * 255);
			}
			else if (i <= singleBlockWidth / 2 && (j >= ((blockNumber - 1)*singleBlockHeight +singleBlockHeight / 2))) {
				int num = blockNumber*(blockNumber - 1);
				src.at<uchar>(j, i) = (int)(total[num][src.at<uchar>(j, i)] * 255);
			}
			else if (i >= ((blockNumber - 1)*singleBlockWidth + singleBlockHeight / 2) && j <=singleBlockHeight / 2) {
				int num = blockNumber - 1;
				src.at<uchar>(j, i) = (int)(total[num][src.at<uchar>(j, i)] * 255);
			}
			else if (i >= ((blockNumber - 1)*singleBlockWidth + singleBlockWidth / 2) && j >= ((blockNumber - 1)*singleBlockHeight + singleBlockHeight / 2)) {
				int num = blockNumber*blockNumber - 1;
				src.at<uchar>(j, i) = (int)(total[num][src.at<uchar>(j, i)] * 255);
			}
			//four edges except coners  
			else if (i <= singleBlockWidth/ 2)
			{
				//线性插值  
				int num_i = 0;
				int num_j = (j -singleBlockHeight / 2) /singleBlockHeight;
				int num1 = num_j*blockNumber + num_i;
				int num2 = num1 + blockNumber;
				float p = (j - (num_j*singleBlockHeight +singleBlockHeight / 2)) / (1.0f*singleBlockHeight);
				float q = 1 - p;
				src.at<uchar>(j, i) = (int)((q*total[num1][src.at<uchar>(j, i)] + p*total[num2][src.at<uchar>(j, i)]) * 255);
			}
			else if (i >= ((blockNumber - 1)*singleBlockWidth+ singleBlockWidth/ 2)) {
				//线性插值  
				int num_i = blockNumber - 1;
				int num_j = (j -singleBlockHeight / 2) /singleBlockHeight;
				int num1 = num_j*blockNumber + num_i;
				int num2 = num1 + blockNumber;
				float p = (j - (num_j*singleBlockHeight +singleBlockHeight / 2)) / (1.0f*singleBlockHeight);
				float q = 1 - p;
				src.at<uchar>(j, i) = (int)((q*total[num1][src.at<uchar>(j, i)] + p*total[num2][src.at<uchar>(j, i)]) * 255);
			}
			else if (j <=singleBlockHeight / 2) {
				//线性插值  
				int num_i = (i - singleBlockWidth/ 2) / singleBlockWidth;
				int num_j = 0;
				int num1 = num_j*blockNumber + num_i;
				int num2 = num1 + 1;
				float p = (i - (num_i*singleBlockWidth+ singleBlockWidth/ 2)) / (1.0f*singleBlockWidth);
				float q = 1 - p;
				src.at<uchar>(j, i) = (int)((q*total[num1][src.at<uchar>(j, i)] + p*total[num2][src.at<uchar>(j, i)]) * 255);
			}
			else if (j >= ((blockNumber - 1)*singleBlockHeight +singleBlockHeight / 2)) {
				//线性插值  
				int num_i = (i - singleBlockWidth/ 2) / singleBlockWidth;
				int num_j = blockNumber - 1;
				int num1 = num_j*blockNumber + num_i;
				int num2 = num1 + 1;
				float p = (i - (num_i*singleBlockWidth+ singleBlockWidth/ 2)) / (1.0f*singleBlockWidth);
				float q = 1 - p;
				src.at<uchar>(j, i) = (int)((q*total[num1][src.at<uchar>(j, i)] + p*total[num2][src.at<uchar>(j, i)]) * 255);
			}
			//双线性插值
			else {
				int num_i = (i - singleBlockWidth/ 2) / singleBlockWidth;
				int num_j = (j -singleBlockHeight / 2) /singleBlockHeight;
				int num1 = num_j*blockNumber + num_i;
				int num2 = num1 + 1;
				int num3 = num1 + blockNumber;
				int num4 = num2 + blockNumber;
				float u = (i - (num_i*singleBlockWidth+ singleBlockWidth/ 2)) / (1.0f*singleBlockWidth);
				float v = (j - (num_j*singleBlockHeight +singleBlockHeight / 2)) / (1.0f*singleBlockHeight);
				src.at<uchar>(j, i) = (int)((u*v*total[num4][src.at<uchar>(j, i)] +
					(1 - v)*(1 - u)*total[num1][src.at<uchar>(j, i)] +
					u*(1 - v)*total[num2][src.at<uchar>(j, i)] +
					v*(1 - u)*total[num3][src.at<uchar>(j, i)]) * 255);
			}
			//最后这步，类似高斯平滑
			src.at<uchar>(j, i) = src.at<uchar>(j, i) + (src.at<uchar>(j, i) << 8) + (src.at<uchar>(j, i) << 16);
		}
	}
	imshow("均衡图线性差值", src);
	waitKey(0);
	return 0;
}

效果如下：