毕设—高动态范围图像渲染算法之基于对数变换的色调映射技术（三）人工智能AACM的博客-

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

这个算法来源于论文《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了。

空间转化矩阵如下：


运行效果如下：