OpenCv，局部自适应图像增强(Local Adaptive Contrast Enhancement)

一、理论

         图像增强算法的基本原则是“降低低频区域，突出高频区域”，以此强化边缘，达到增强的目的。最简单的例子就是通过原始图像减去高斯模糊处理后的图像，就能够将边缘强化出来。

         直方图均衡化也是一种非常常见的增强方法。但是为了避免背景的干扰，更倾向于采用“局部”方法进行处理。我们这里着重研究自适应对比度增强(ACE)的相关内容。

        ACE的定义和原理看上去还是比较简单的。这里的和都可以根据图像本身计算出来。而则需要单独计算。

          可以为单独的常量，或者通过来代替。这里的D是一个全局的值，比如平均值。

二、实现

        涉及到局部的运算，自然而然会想到使用卷积的方法。更好的是Opencv提供了专门的函数用来做这个工作—BLUR

文档中写到：

bool ACE_Enhance(cv::Mat& src_img, cv::Mat& dst_img, unsigned int half_winSize, double Max_Q);
double GetMeanValue(cv::Mat& src_img);
double GetVarianceValue(cv::Mat& src_img, double MeanVlaue);

//************************************************************************  
// 函数名称:    ACE_Enhance  
// 访问权限:    public   
// 创建日期:    2016/11/23  
// 创 建 人:         
// 函数说明:    单通道增强  
// 函数参数:    cv::Mat & src_img   输入图像  
// 函数参数:    cv::Mat & dst_img   输出图像  
// 函数参数:    unsigned int half_winSize   增强窗口的半窗大小  
// 函数参数:    double Max_Q            最大Q值  
// 返 回 值:    bool  
//************************************************************************  
bool ACE_Enhance(Mat& src_img, Mat& dst_img, unsigned int half_winSize, double Max_Q)
{
	if (!src_img.data)
	{
		cout << "没有输入图像" << endl;
		return false;
	}

	int rows(src_img.rows);
	int cols(src_img.cols);
	unsigned char* pdata = nullptr;
	unsigned char* pdata1 = nullptr;
	cv::Mat tmpDstImg(rows, cols, CV_8UC1, cv::Scalar::all(0));

	for (int i = half_winSize; i < (rows - half_winSize); i++)
	{
		pdata = tmpDstImg.ptr<unsigned char>(i);
		pdata1 = src_img.ptr<unsigned char>(i);
		for (int j = half_winSize; j < (cols - half_winSize); j++)
		{
			cv::Mat tempArea = src_img(cv::Rect(j - half_winSize, i - half_winSize, half_winSize * 2 + 1, half_winSize * 2 + 1));   //截取窗口图像  
			double MeanVlaue = GetMeanValue(tempArea);//平均值
			double varian = GetVarianceValue(tempArea, MeanVlaue);//均方差
			if (varian > 0.1)
			{
				double cg = 100.0 / sqrt(varian);
				cg = cg > Max_Q ? Max_Q : cg;
				double pixelvalue = cg*((double)pdata1[j] - MeanVlaue);

				int tempValue = MeanVlaue + pixelvalue;
				tempValue = tempValue > 255 ? 255 : tempValue;
				tempValue = tempValue < 0 ? 0 : tempValue;
				pdata[j] = tempValue;
			}
		}
	}
	dst_img = tmpDstImg;

	return true;
}

//************************************************************************  
// 函数名称:    GetMeanValue  
// 访问权限:    public   
// 创建日期:    2016/11/18  
// 创 建 人:         
// 函数说明:    获取图像的平均灰度值  
// 函数参数:    cv::Mat & src_img   输入图像  
// 返 回 值:    double  
//************************************************************************  
double GetMeanValue(Mat& src_img)
{
	if (CV_8UC1 != src_img.type())
	{
		return -1.0;
	}

	int rows(src_img.rows); //height  
	int cols(src_img.cols); //width  
	unsigned char* data = nullptr;
	double PixelValueSum(0.0);  //总共的像素值  

	for (int i = 0; i < rows; i++)
	{
		data = src_img.ptr<unsigned char>(i);
		for (int j = 0; j < cols; j++)
		{
			PixelValueSum += (double)data[j];
		}   //计算图像的总共像素值  
	}

	double result(PixelValueSum / static_cast<double>(rows*cols));    //计算图像的均值  

	return result;
}

// 
// 函数名称:    GetVarianceValue  
// 访问权限:    public   
// 创建日期:    2016/11/18  
// 创 建 人:         
// 函数说明:        计算图像的均方差  
// 函数参数:    cv::Mat & src_img   输入图像  
// 函数参数:    double MeanVlaue        图像的均值  
// 返 回 值:    double  
//  
double GetVarianceValue(Mat& src_img, double MeanVlaue)
{
	if (CV_8UC1 != src_img.type())
	{
		return -1.0;
	}

	int rows(src_img.rows); //height  
	int cols(src_img.cols); //width  
	unsigned char* data = nullptr;
	double PixelValueSum(0.0);  //总共的像素值  

	for (int i = 0; i < rows; i++)
	{
		data = src_img.ptr<unsigned char>(i);
		for (int j = 0; j < cols; j++)
		{
			PixelValueSum += pow((double)(data[j] - MeanVlaue), 2);
		}   //计算图像方差  
	}

	double result(PixelValueSum / static_cast<double>(rows*cols));    //计算图像的均方差  

	return result;
}

程序保留：

#include <stdio.h>
#include <iostream>
#include "fftw3.h"
#include <stdlib.h>
#include <string.h>
#include <limits.h>

#include "string"
#include "vector"
#include <windows.h>
#include <opencv2/legacy/legacy.hpp>
#include <opencv2/nonfree/nonfree.hpp>//opencv_nonfree模块：包含一些拥有专利的算法，如SIFT、SURF函数源码。 
#include "opencv2/core/core.hpp"
#include "opencv2/features2d/features2d.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <opencv2/nonfree/features2d.hpp>

//参考1：基于“局部标准差”的图像增强（原理、算法、代码）
//http://www.cnblogs.com/jsxyhelu/p/4857721.html
//参考2：使用局部标准差实现图像的局部对比度增强算法。
//http://www.cnblogs.com/Imageshop/p/3324282.html
using namespace cv;
using namespace std;


//点乘法,elementWiseMultiplication
Mat matrixWiseMulti(Mat &m1, Mat &m2){
	Mat dst = m1.mul(m2);//注意是对应矩阵位置的元素相乘
	return dst;
}
//ACE算法原理：
//ACE算法原理表达式：f(i,j)=Mx(i,j)+G(i,j)*[x(i,j)-Mx(i,j)]
//高频成分：x(i,j)-Mx(i,j)，x(i,j)表示当前中心像素，Mx(i,j)表示局部平均值
//增益系数：G(i,j)，可以直接令其为系数C（一般总是大于1）
//

//图像局部对比度增强算法
//float MaxCG:对高频成分的最大增益值
//int n：局部半径
//int C；对高频的直接增益系数
//Mat src：原图像
Mat 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));//meanMask为图像局部均值 
	imshow("低通滤波", meanLocal);
	Mat highFreq = src - meanLocal;//高频成分：x(i,j)-Mx(i,j)
	imshow("高频成分", highFreq);

	varLocal = matrixWiseMulti(highFreq, highFreq);
	blur(varLocal, varLocal, Size(n, n));    //varMask为此时为图像局部方差   
	//换算成局部标准差（开根号）
	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); //meanGlobal为全局均值 varGlobal为全局标准差，实际均是一个数
	Mat gainArr =0.5 * meanGlobal / varLocal;//增益系数矩阵：G(i,j)，可以直接令其为系数C（一般总是大于1）
	/*
	for (int i = 0; i < rows; i++){
		for (int j = 0; j < cols; j++)
				cout<<gainArr.at<float>(i, j)<<" " ;
		cout << endl;
		if (i == 1)
			break;
	}*/
	
	//对增益矩阵进行截止
	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("Lee改进的D方法", dst1);
	Mat dst2 = meanLocal + C*highFreq;//直接利用系数C进行高频成分放大
	imshow("直接系数C方法", dst2);
	return dst2;
}


Mat myACE(Mat &src, int n = 7, float MaxCG = 7.5){
	int rows = src.rows;
	int cols = src.cols;
	Mat dst(src.rows, src.cols, CV_8UC1, Scalar::all(0));
	
	if (src.type() == CV_8UC1)
		int aa = src.type();
	Mat meanLocal;
	Mat varLocal;
	Mat meanGlobal;
	Mat varGlobal;

	blur(src.clone(), meanLocal, Size(n, n));//meanMask为图像局部均值 
	Mat highFreq = src - meanLocal;//高频成分：x(i,j)-Mx(i,j)
	
	varLocal = matrixWiseMulti(highFreq, highFreq);
	blur(varLocal, varLocal, Size(n, n));    //varMask为此时为图像局部方差   
	//换算成局部标准差（开根号）
	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); //meanGlobal为全局均值 varGlobal为全局标准差，实际均是一个数
	Mat gainArr = 0.5 * meanGlobal / varLocal;//增益系数矩阵：G(i,j)，可以直接令其为系数C（一般总是大于1）
	//对增益矩阵进行截止
	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);
	dst = meanLocal + gainArr;
	//imshow("Lee改进的D方法", dst);
	return dst;
}

int main()
{
	const char* img_path = "oct.bmp";
	//const char* img_path = "yanjing.jpg";
	//const char* img_path = "luowen.tif";
	//const char* img_path = "colors_large.png";  
	//const char* img_path = "input_0.png";  
	//const char* img_path = "a1.jpg";  
	//const char* img_path = "seed.tif";  
	//const char* img_path = "flyman_gray.bmp";
	//const char* img_path = "CT.bmp";
	//const char* img_path = "rcc.jpg";
	Mat src = imread(img_path, 0);
	imshow("src", src);
	int n = 50;
	float MaxCG = 10.5;
	Mat dst=myACE(src,n,MaxCG);
	imshow("myACE",dst);
	waitKey();
	return  0;
}