毕设—高动态范围图像渲染算法之基于对数变换的色调映射技术（三）

这个算法是我第一个实现的色调映射算法，也是五种算法中原理最清晰，实现最简单的。（原理见：一种自适应对数映射的高对比度图像显示技术）

这个算法来源于论文《Adaptive Logarithmic Mapping For Displaying High Contrast Scenes》是一篇2003年的古老算法。其实这个算法OpenCV 3.0已经实现了，但是我还是根据自己的理解又实现了一下。

实现该算法共分为四步：

第一步：把场景亮度映射到图像亮度

其实完成这步有很多方法，上述论文中提到了两种方法：一、对于静态图像，使用场景的对数平均值作为初始缩放因子，进行压缩。二、对于需要进行交互式的图像处理，可以使用亮度的对数信息的高斯模糊后的结果，作为压缩结果，范围用户自己调，默认是15%。

但是本项目使用了另一种压缩方法，即归一化。具体步骤如下：首先我们先把提取出来的场景信息进行色彩空间转化，得到亮度值；随后我们对亮度值进行归一化处理；

	for(int i = 0; i < rows; i++){
		for(int j = 0; j < cols; j++){

			//通过这么写，不仅可以实现色调映射，还可以实现图像增强
			if(strstr(fileName,".hdr")!=NULL)
			{
				//提取场景的RGB数据
				RGB[0]=result.cols[i*cols*3+j*3+0];
				RGB[1]=result.cols[i*cols*3+j*3+1];
				RGB[2]=result.cols[i*cols*3+j*3+2];
			}
			else
			{
				RGB[2]=(double)re.at<Vec3b>(i,j)[0];
				RGB[1]=(double)re.at<Vec3b>(i,j)[1];
				RGB[0]=(double)re.at<Vec3b>(i,j)[2];
			}
			
			//色彩空间转化
			//转化成XYZ表示，虽然缩放只缩放Y，但是为了后面进行转化，这里把x，z也求出来;
			xyz[i][j].x=Mrgb[0][0]*RGB[0]+Mrgb[1][0]*RGB[1]+Mrgb[2][0]*RGB[2];
			xyz[i][j].y=Mrgb[0][1]*RGB[0]+Mrgb[1][1]*RGB[1]+Mrgb[2][1]*RGB[2];
			xyz[i][j].z=Mrgb[0][2]*RGB[0]+Mrgb[1][2]*RGB[1]+Mrgb[2][2]*RGB[2];

			//获取XYZ的比例
			double rateX=xyz[i][j].x/(xyz[i][j].x+xyz[i][j].y+xyz[i][j].z);
			double rateY=xyz[i][j].y/(xyz[i][j].x+xyz[i][j].y+xyz[i][j].z);

			//归一化
			if(strstr(fileName,".hdr")!=NULL)
				xyz[i][j].y=atan(xyz[i][j].y)*2/PI;
			else
				xyz[i][j].y/=255.0;

			//还原x，y，z的比率
			xyz[i][j].x=xyz[i][j].y/rateY*rateX;
			xyz[i][j].z=xyz[i][j].y/rateY*(1-rateX-rateY);

			Lmax = max(xyz[i][j].y,Lmax);//取出场景最亮
			Lsum += log(xyz[i][j].y+a);//对场景亮度进行对数求和
		}
	}

接下来是第二步，就是运用传说中的自适应函数：

这个公式看起来很唬人，其实并不难：

//对数变换，本算法的核心
xyz[i][j].y = 1.0 * log(xyz[i][j].y/Lav + 1) / log(2.0 + 8.0*pow((xyz[i][j].y/Lav / Lmax), log(bias) / log(0.5))) / log10(Lmax + 1);

其中Lav是对数平均值，简单来说除以对数平均值，映射效果会更好。为什么会更好，原因在讲基于Retinex的色调映射算法时会讲。

第三步：伽马校正

//伽马校正
RGB[0]=C*pow(RGB[0],r);
RGB[1]=C*pow(RGB[1],r);
RGB[2]=C*pow(RGB[2],r);

第四步：色彩空间还原

//转化成RGB表示
RGB[0] = Mxyz[0][0]*xyz[i][j].x + Mxyz[1][0]*xyz[i][j].y + Mxyz[2][0]*xyz[i][j].z;
RGB[1] = Mxyz[0][1]*xyz[i][j].x + Mxyz[1][1]*xyz[i][j].y + Mxyz[2][1]*xyz[i][j].z;
RGB[2] = Mxyz[0][2]*xyz[i][j].x + Mxyz[1][2]*xyz[i][j].y + Mxyz[2][2]*xyz[i][j].z;

完整代码如下：

Mat Tone_Mapping(char* fileName, double bias, double C, double r)//基于对数变换的HDR
{
	int rows,cols;
	HDRLoaderResult result;
	Mat re;
	//判断图像类型
	if(strstr(fileName,".hdr")!=NULL)
	{
		//获取图像信息
		int ret = HDRLoader::load(fileName, result);
		rows=result.height;
		cols=result.width;
	}
	else
	{
		//加载图像
		re = imread(fileName, 1);
		rows=re.rows;
		cols=re.cols;
	}

	//定义矩阵
	Mat out_image (rows,cols,CV_8UC3);//结果图像
	Mat img(rows,cols,CV_8UC3);//过程图像

	//定义变量
	double PI=acos(-1.0);
	double Ymin=INF,Ymax=-INF,Ysum=0,Yav=0;
	double Lmin=INF,Lmax=-INF,Lsum=0,Lav=0;
	for(int i = 0; i < rows; i++){
		for(int j = 0; j < cols; j++){

			//通过这么写，不仅可以实现色调映射，还可以实现图像增强
			if(strstr(fileName,".hdr")!=NULL)
			{
				RGB[0]=result.cols[i*cols*3+j*3+0];
				RGB[1]=result.cols[i*cols*3+j*3+1];
				RGB[2]=result.cols[i*cols*3+j*3+2];
			}
			else
			{
				RGB[2]=(double)re.at<Vec3b>(i,j)[0];
				RGB[1]=(double)re.at<Vec3b>(i,j)[1];
				RGB[0]=(double)re.at<Vec3b>(i,j)[2];
			}
			

			//转化成XYZ表示，虽然缩放只缩放Y，但是为了后面进行转化，这里把x，z也求出来;
			xyz[i][j].x=Mrgb[0][0]*RGB[0]+Mrgb[1][0]*RGB[1]+Mrgb[2][0]*RGB[2];
			xyz[i][j].y=Mrgb[0][1]*RGB[0]+Mrgb[1][1]*RGB[1]+Mrgb[2][1]*RGB[2];
			xyz[i][j].z=Mrgb[0][2]*RGB[0]+Mrgb[1][2]*RGB[1]+Mrgb[2][2]*RGB[2];

			//获取XYZ的比例
			double rateX=xyz[i][j].x/(xyz[i][j].x+xyz[i][j].y+xyz[i][j].z);
			double rateY=xyz[i][j].y/(xyz[i][j].x+xyz[i][j].y+xyz[i][j].z);

			//归一化
			if(strstr(fileName,".hdr")!=NULL)
				xyz[i][j].y=atan(xyz[i][j].y)*2/PI;
			else
				xyz[i][j].y/=255.0;

			//还原x，y，z的比率
			xyz[i][j].x=xyz[i][j].y/rateY*rateX;
			xyz[i][j].z=xyz[i][j].y/rateY*(1-rateX-rateY);

			Lmax = max(xyz[i][j].y,Lmax);//取出场景最亮
			Lsum += log(xyz[i][j].y+a);//对场景亮度进行对数求和
		}
	}

	Lav = exp(Lsum/(1.0*rows*cols));//求取对数平均值

	Lmax=Lmax/Lav;

	for(int i = 0; i < rows; i++){
		for(int j = 0; j < cols; j++){

			//记录XYZ空间下x，y，z的比率
			double rateX=xyz[i][j].x/(xyz[i][j].x+xyz[i][j].y+xyz[i][j].z);
			double rateY=xyz[i][j].y/(xyz[i][j].x+xyz[i][j].y+xyz[i][j].z);

			//对数变换，本算法的核心
			xyz[i][j].y = 1.0 * log(xyz[i][j].y/Lav + 1) / log(2.0 + 8.0*pow((xyz[i][j].y/Lav / Lmax), log(bias) / log(0.5))) / log10(Lmax + 1);

			//还原x，y，z的比率
			xyz[i][j].x=xyz[i][j].y/rateY*rateX;
			xyz[i][j].z=xyz[i][j].y/rateY*(1-rateX-rateY);

			//转化成RGB表示，色彩空间还原
			RGB[0] = Mxyz[0][0]*xyz[i][j].x + Mxyz[1][0]*xyz[i][j].y + Mxyz[2][0]*xyz[i][j].z;
			RGB[1] = Mxyz[0][1]*xyz[i][j].x + Mxyz[1][1]*xyz[i][j].y + Mxyz[2][1]*xyz[i][j].z;
			RGB[2] = Mxyz[0][2]*xyz[i][j].x + Mxyz[1][2]*xyz[i][j].y + Mxyz[2][2]*xyz[i][j].z;

			if (RGB[0] < 0)RGB[0] = 0; if (RGB[0]>1)RGB[0] = 1;
			if (RGB[1] < 0)RGB[1] = 0; if (RGB[1]>1)RGB[1] = 1;
			if (RGB[2] < 0)RGB[2] = 0; if (RGB[2]>1)RGB[2] = 1;

			//伽马校正
			RGB[0]=C*pow(RGB[0],r);
			RGB[1]=C*pow(RGB[1],r);
			RGB[2]=C*pow(RGB[2],r);

			if (RGB[0] < 0)RGB[0] = 0; if (RGB[0]>1)RGB[0] = 1;
			if (RGB[1] < 0)RGB[1] = 0; if (RGB[1]>1)RGB[1] = 1;
			if (RGB[2] < 0)RGB[2] = 0; if (RGB[2]>1)RGB[2] = 1;

			//修正补偿
			RGB[0] = Transform(RGB[0]);
			RGB[1] = Transform(RGB[1]);
			RGB[2] = Transform(RGB[2]);

			out_image.at<Vec3b>(i, j)[0] = (int)(RGB[2]*255);
			out_image.at<Vec3b>(i, j)[1] = (int)(RGB[1]*255);
			out_image.at<Vec3b>(i, j)[2] = (int)(RGB[0]*255);
		}
	}
	return out_image;
}

这些步骤没有很明显的前后之分，比如在我的项目中，色彩空间还原就放在了伽马校正的前面。

上面用到的色彩空间转化，是RGB到XYZ，其中Y代表亮度，在映射过程中就使用了亮度这一个变量。色彩空间还原那肯定是XYZ到RGB了。

空间转化矩阵如下：


运行效果如下：