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;
}