对比度有限自适应直方图均衡化

前注：

ACE在图像处理方面可以有两种表示，一种是在说的:Automatic Color Equalization，即自动彩色均衡；还有一种是本篇要说的：Adaptive Contrast Enhancement，即自适应对比度增强 。不要搞混了~_~

医学图像由于本身及成像条件的限制，图像的对比度很低，因此对比度增强算法在医学图像方面显得很有必要。人眼对高频信号（边缘）比较敏感，但是如果高频信号嵌入在大量的低频背景+噪声信号中，其视觉可见性会降低，此时适当提高高频部分能够提高视觉效果。在这方面，传统的线性对比度拉升以及直方图均衡是最为广泛的全局图像增强方法。对比度拉升可以调整图像的动态范围，直方图均衡则是利用直方图分布概率重新映射图像的数据。但是，全局直方图均衡化可能加强部分噪声信号。为了避免低频背景的干扰，采用“局部”的增强方法可能会取得更好的效果。

在局部对比度增强方面，有两种方法较为常用：一是在中说的自适应直方图均衡，另一种是本篇要说的自适应对比度增强。

自适应对比度增强原理：

ACE原理其实挺简单的，即是将一幅图像分成两个部分：一是低频部分，可以通过图像的低通滤波（平滑模糊）获得；二是高频部分，可以由原图减去低频部分得到。而算法的目标是增强代表细节的高频部分，即对高频部分乘以某个增益值，然后重组得到增强的图像。所以ACE算法的核心就是高频部分增益系数的计算，一种方案是将增益设为一个固定值，另一种方案是将增益值表示为与方差相关的量，后面在附加的公式中再说明。

假设一幅图像中像素点表示为 x(i,j), 那么以(i，j)为中心，窗口大小为(2n+1)*(2n+1)的区域内，其局部均值和方差可以表示为：

，

，

均值mx可以近似认为是背景部分，此时 x-m 即是高频细节部分，对高频作增益乘积，有：

对于增益G,方案一是取大于1的常数，达到增强的效果，即为：

方案二是表示为与局部均方差成反比的变化值，即为：

在图像的高频区域，局部均方差较大，此时增益值就比较小，这样结果不会出现过亮情况。但是在图像平滑的区域，局部均方差很小，此时增益值比较大，从而可能会放大噪声信号，所以需要对增益最大值做一定的限制才能取得较好的效果。

自适应对比度增强实现：

根据上面的原理部分，及参照网上的代码，我用opencv实现了一下，代码如下：

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

using namespace cv;
using namespace std;

Mat matrixWiseMulti(Mat &m1, Mat &m2){
	Mat dst = m1.mul(m2);
	return dst;
}

 
//float MaxCG:对高频成分的最大增益值,int n：局部半径,int C:对高频的直接增益系数  
void ACE(Mat &src, int C = 3, int n = 3, float MaxCG = 7.5){
	int rows = src.rows;
	int cols = src.cols;

	Mat meanLocal; //图像局部均值  
	Mat varLocal;  //图像局部方差  
	Mat meanGlobal;//全局均值
	Mat varGlobal; //全局标准差  

	blur(src.clone(), meanLocal, Size(n, n)); 
	imshow("低通滤波", meanLocal);
	Mat highFreq = src - meanLocal;//高频成分 
	imshow("高频成分", highFreq);

	varLocal = matrixWiseMulti(highFreq, highFreq);
	blur(varLocal, varLocal, Size(n, n));      
	//换算成局部标准差  
	varLocal.convertTo(varLocal, CV_32F);
	for (int i = 0; i < rows; i++){
		for (int j = 0; j < cols; j++){
			varLocal.at<float>(i, j) = (float)sqrt(varLocal.at<float>(i, j));
		}
	}
	meanStdDev(src, meanGlobal, varGlobal);
	Mat gainArr = 0.5 * meanGlobal / varLocal;//增益系数矩阵  

	//对增益矩阵进行截止  
	for (int i = 0; i < rows; i++){
		for (int j = 0; j < cols; j++){
			if (gainArr.at<float>(i, j) > MaxCG){
				gainArr.at<float>(i, j) = MaxCG;
			}
		}
	}
	gainArr.convertTo(gainArr, CV_8U);
	gainArr = matrixWiseMulti(gainArr, highFreq);
	Mat dst1 = meanLocal + gainArr;
	imshow("变增益方法", dst1);
	Mat dst2 = meanLocal + C*highFreq;  
	imshow("恒增益方法", dst2);
}

int main()
{
	const char* img_path = "1.bmp"; 
	Mat src = imread(img_path, 0);
	imshow("src", src);
	int C = 5;
	int n = 7;
	float MaxCG = 8;
	ACE(src, C,n, MaxCG);
	waitKey();
	return  0;
}