C++实现Retinex算法（还有伽马矫正、直方图均衡化）

        之前的代码是学习MC庆仔的http://t.csdnimg.cn/0B5Lr，后来自己为了加快运算时间，做了一些改动。同时也附上了伽马矫正、直方图均衡化的结果做对比。本次未写带颜色恢复的多尺度Retinex算法，如果后面需要的人多我再改一改附上。

一、引言

        一般用于低照度环境下图像增强算法，可以分为传统的图像增强算法和基于深度学习的图像增强算法。目前个人未实践过基于深度学习的图像增强算法，所以本次就着重些传统的图像增强算法。

        传统的图像增强算法可以分为两类：基于底层处理的图像增强算法、基于Retinex理论的图像增强算法。基于底层处理的图像增强算法最常见的就是直方图均衡化与伽马矫正，因为方法简单，计算速度快，所以一般的增强最常用这两种方法。但是为了追求更好的图像增强效果，应用比较多的是基于Retinex理论的图像增强算法。

二、基于底层处理的图像增强

2.1  伽马矫正

伽马矫正的公式如下。计算时要将像素值从0-255映射到0-1，计算后再转回0-255

伽马矫正的系数为，当大于1时图像变暗，当小于1时，图像变亮。函数图像如下：

 

附上C++代码

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

using namespace cv;

Mat gama_enhance(Mat input,double gama)
{
	//建立查找表，减少整体计算时间
	uchar lut[256];
	for (int i = 0; i < 256; i++)
	{
		lut[i] = saturate_cast<uchar>(pow((float)(i / 255.0f), gama) * 255.0f);
	}

	//用指针访问每个像素值，并且进行更改
	for (int j = 0; j < input.rows; j++)
	{
		uchar* data = input.ptr<uchar>(j);
		for (int i = 0; i < input.cols * input.channels(); i++)
		{
			data[i] = lut[data[i]];
		}
	}
	return input;
}

2.2 直方图均衡化

 直方图均衡技术将原始图像的灰度直方图从比较集中的某个灰度区间变成在全部灰度范围内的均匀分布。（来源于百度）

OpenCV中有直方图均衡化的函数，但是只能处理灰度图像。我这里为了彩色图直方图均衡化的颜色不失真，我先把图像转为HSV空间，再对V值直方图均衡化。附C++代码

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

using namespace cv;

//直方图均衡化代码
Mat HE_enhance(Mat& input)
{
	Mat output;
	//如果是彩色图像，为了防止颜色失真，先转换到HSV空间，对V值进行直方图均衡化增强
	if (input.channels() > 1)
	{
		cvtColor(input, input, COLOR_BGR2HSV);
		vector<Mat> input_HSV;
		split(input, input_HSV);
		equalizeHist(input_HSV[2], input_HSV[2]);
		merge(input_HSV, output);
		cvtColor(output, output, COLOR_HSV2BGR);
	}
	else
	{
		equalizeHist(input, output);
	}
	return output;
}

三、基于Retinex理论的图像增强 

 3.1 Retinex理论

Retinex理论的基本内容：物体的颜色是由物体对较长波段（红色）、中等波段（绿色）、较短波段（蓝色）光线的反射能力来决定的，而不是由反射光的强度来决定的，物体的色彩不受光照非均匀性的影响。

 L(x,y)表示入射光场分布，R(x,y)表示物体反射特性，I(x,y)表示相机接收的图像信息。

 3.2 Retinex算法

一般来说，常用的Retinex算法如下：

流程就是如下：

照度分量估计这里，最开始用这个算法的人用的是高斯滤波，后来有的人用引导滤波或者双边滤波，或者更难的处理方式。最常用的流程就是以上，后续的大部分算法只是在中增加一些步骤。最开始提出该算法的人Jobson，把该算法分为了单尺度Retinex算法（SSR)以及多尺度Retinex算法 (MSR)。相关文献：

1. Jobson D J, Rahman Z, Woodell G A, et al Properties and performance of a center/surround retinex[J] IEEE Transactions on Image Processing A Publication of the IEEE Signal Processing Society, 1997, 6(3):451-462
2. Jobson D J , Rahman Z , Woodell G A . A multiscale retinex for bridging the gap between color images and the human observation of scenes[J]. IEEE Transactions on Image Processing, 1997, 6(7):965-976.

3.3 Retinex算法的C++代码

我为了让计算速度加快，log值计算用了查找表，像素值访问用了指针，然后就是将高斯滤波换为了引导滤波。引导滤波代码来自于小武~~，原帖http://t.csdnimg.cn/wgZ0y。如果不能直接附代码，我就删了。

先附灰度图像的引导滤波代码 

#include<opencv2/opencv.hpp>
#include<vector>
//引导滤波函数
//输入引导图像I
//输入滤波图像p
//输入r
//输入eps，eps越大，滤波效果约明显，eps=0时与原图一致

cv::Mat GuidedFilter_gray(cv::Mat& I, cv::Mat& p, int r, double eps) {
	int wsize = 2 * r + 1;
	//数据类型转换
	I.convertTo(I, CV_32F, 1.0 / 255.0);
	p.convertTo(p, CV_32F, 1.0 / 255.0);

	//meanI=fmean(I)
	cv::Mat mean_I;
	cv::boxFilter(I, mean_I, -1, cv::Size(wsize, wsize), cv::Point(-1, -1), true, cv::BORDER_REFLECT);//盒子滤波

	//meanP=fmean(P)
	cv::Mat mean_p;
	cv::boxFilter(p, mean_p, -1, cv::Size(wsize, wsize), cv::Point(-1, -1), true, cv::BORDER_REFLECT);//盒子滤波

	//corrI=fmean(I.*I)
	cv::Mat mean_II;
	mean_II = I.mul(I);
	cv::boxFilter(mean_II, mean_II, -1, cv::Size(wsize, wsize), cv::Point(-1, -1), true, cv::BORDER_REFLECT);//盒子滤波

	//corrIp=fmean(I.*p)
	cv::Mat mean_Ip;
	mean_Ip = I.mul(p);
	cv::boxFilter(mean_Ip, mean_Ip, -1, cv::Size(wsize, wsize), cv::Point(-1, -1), true, cv::BORDER_REFLECT);//盒子滤波

	//varI=corrI-meanI.*meanI
	cv::Mat var_I, mean_mul_I;
	mean_mul_I = mean_I.mul(mean_I);
	cv::subtract(mean_II, mean_mul_I, var_I);

	//covIp=corrIp-meanI.*meanp
	cv::Mat cov_Ip;
	cv::subtract(mean_Ip, mean_I.mul(mean_p), cov_Ip);

	//a=conIp./(varI+eps)
	//b=meanp-a.*meanI
	cv::Mat a, b;
	cv::divide(cov_Ip, (var_I + eps), a);
	cv::subtract(mean_p, a.mul(mean_I), b);

	//meana=fmean(a)
	//meanb=fmean(b)
	cv::Mat mean_a, mean_b;
	cv::boxFilter(a, mean_a, -1, cv::Size(wsize, wsize), cv::Point(-1, -1), true, cv::BORDER_REFLECT);//盒子滤波
	cv::boxFilter(b, mean_b, -1, cv::Size(wsize, wsize), cv::Point(-1, -1), true, cv::BORDER_REFLECT);//盒子滤波

	//q=meana.*I+meanb
	cv::Mat q;
	q = mean_a.mul(I) + mean_b;

	//数据类型转换
	I.convertTo(I, CV_8U, 255);
	p.convertTo(p, CV_8U, 255);
	q.convertTo(q, CV_8U, 255);
	return q;
}

 注意把该灰度图像的引导滤波的函数保存后，再写彩色图像的引导滤波。

#include<opencv2/opencv.hpp>
#include<vector>
//引导滤波函数
//输入引导图像I
//输入滤波图像input
//输入r
//输入eps，eps越大，滤波效果约明显，eps=0时与原图一致

cv::Mat GuidedFilter_bgr(cv::Mat& I, cv::Mat& input, int r, double eps) {
	//判断引导图像是否为灰度图
	if (I.channels() > 1)
	{
		cvtColor(I, I, cv::COLOR_BGR2GRAY);
	}

	cv::Mat dst1, src_input;
	std::vector<cv::Mat> p, q;
	src_input = input;

	if (input.channels() > 1) {             //输入为彩色图
		cv::split(src_input, p);
		for (int i = 0; i < input.channels(); ++i) {
			dst1 = GuidedFilter_gray(I, p[i], r, eps * eps);
			q.push_back(dst1);
		}
		cv::merge(q, dst1);
	}
	else {                               //输入为灰度图
		input.copyTo(I);
		dst1 = GuidedFilter_gray(I, src_input, r, eps * eps);
	}
	return dst1;
}

 把引导滤波的函数都写好后，再写一下log值计算代码

#include<opencv2/opencv.hpp>
#include<cmath>
using namespace cv;

Mat fast_logmat(Mat& input)
{
	Mat output=Mat(input.size(), CV_32FC3);
	//建立查找表，减少整体计算时间
	float lut[256];
	lut[0] = 0;
	for (int i = 1; i < 256; i++)
	{
		lut[i] = log10(i);
	}

	//用指针访问每个像素值，并且进行更改
	for (int j = 0; j < input.rows; j++)
	{
		uchar* in = input.ptr<uchar>(j);
		float* out = output.ptr<float>(j);
		for (int i = 0; i < input.cols * input.channels(); i++)
		{
			out[i] = lut[in[i]];
		}
	}
	return output;
}

现在可以写Retinex算法的代码了，先附SSR代码

#include<opencv2/opencv.hpp>
using namespace cv;

//引导函数做滤波的retinex函数
//输入图片inputimage
//引导函数的引导图片guideimage，半径r，正则系数eps
//半径r为滤波尺寸，eps为正则化系数，eps越大，滤波效果越明显
//返回retinex增强后图像
Mat retinexfunc_guidefilter(Mat& inputimage, int r, double eps)
{
	Mat outputimage, gfimage;
	Mat image_log, guidefimage_log;
	Mat guideimage = inputimage;

	//对图片进行引导滤波
	gfimage = GuidedFilter_bgr(guideimage, inputimage, r, eps);

	outputimage = Mat(inputimage.size(), CV_32FC3);

	//求输入图片的log值
	image_log = fast_logmat(inputimage);

	//求滤波后的入射光图像log值
	guidefimage_log = fast_logmat(gfimage);

	outputimage = image_log - guidefimage_log;

	//归一化处理
	normalize(outputimage, outputimage, 0, 255, NORM_MINMAX);
	outputimage.convertTo(outputimage, CV_8UC3);
	return outputimage;
}

然后是MSR。解释一些MSR就是，输入图像进行多次照度分量估计，然后得到多个logR，然后把这些logR加权相加，得到最后的增强图像。

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

using namespace cv;

//多尺度的基于引导滤波的retinex增强
//输入需要增强的图像inputimage
//输入每个尺度的引导滤波大小r（一般为15-250的整数，如果图像过暗，可以增大）
//输入每个尺度的引导滤波大小eps（我测试对结果影响不大，0.6-0.9都可）
//输入每个尺度的权值weight（一般weight里数值相加后为1）
//输出增强图像output
Mat MSR_guidefilter(Mat& inputimage, vector<int> r, vector<double> eps, vector<double> weight)
{
	int escale = eps.size();
	int rsacle = r.size();
	int wsacle = weight.size();
	Mat output = Mat::zeros(inputimage.size(), CV_32FC3);
	Mat guideimage = inputimage;

	if ((escale == rsacle && escale == wsacle) == false)
	{
		cout << "输出错误！" << endl;
		return output;
	}

	//求输入图片的log值
	Mat input_log = fast_logmat(inputimage);

	//多尺度求解
	for (int i = 0; i < wsacle; i++)
	{
		int rvalue = r[i];
		double epsvalue = eps[i];
		double weightvalue = weight[i];
		//计算每个尺度的滤波图片
		Mat filterimage = GuidedFilter_bgr(guideimage, inputimage, rvalue, epsvalue);
		filterimage = fast_logmat(filterimage);

		for (int y = 0; y < input_log.rows; y++)
		{
			float* in = input_log.ptr<float>(y);
			float* f_middle = filterimage.ptr<float>(y);
			float* out = output.ptr<float>(y);
			for (int x = 0; x < input_log.cols * input_log.channels(); x++)
			{
				float value;
				value = weight[i] * (in[x] - f_middle[x]);
				out[x] += value;
			}
		}
	}

	normalize(output, output, 0, 255, NORM_MINMAX);
	output.convertTo(output, CV_8UC3);
	return output;
}

四、结果

原图

增强结果 

结束了，如果哪里不对希望指正。