毕设—高动态范围图像渲染算法之基于直方图均衡化的色调映射技术（四）

本算法是我第二个实现的色调映射算法，因为我参考的这篇论文是中文论文，容易看懂。。。所以先实现了它。在进行介绍算法之前我先给大家推荐一篇博客，他把直方图均衡化介绍了很详细（链接https://blog.csdn.net/Enjoying_Science/article/details/51165295?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.nonecase&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.nonecase）。

参考文献是《一种基于直方图均衡化的色调映射方法 》（我有论文，但是不知道如何点击下载，希望有人可以教教我），这个文献详细介绍了算法流程，与操作步骤。但是他的过程有些繁琐，而且算法似乎针对的是RGB图像，所以我对其进行了一点点修改。

本算法分为五步：

第一步：色彩空间转化+归一化

			//转化成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];

			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=pow(xyz[i][j].y,r);
			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);

第二步：统计直方图

			//统计对应色阶的像素点个数
			NumPixel[(int)(xyz[i][j].y*255+0.5)]++;

第三步：计算分布函数

	for(int i=0;i<256;i++)
	{
		//计算分布函数
		if(i==0)RatePixel[i]=1.0*NumPixel[i]/(rows*cols);
		else RatePixel[i]=1.0*NumPixel[i]/(rows*cols)+RatePixel[i-1];
		if(RatePixel[i]!=0&&flag==0)
		{
			flag=1;
			PixelMin = 255.0 * RatePixel[i];
		}
	}

第四步：拉伸对比度

		//拉伸对比度
		map[i]=(255.0*RatePixel[i]-PixelMin)*255.0/(255-PixelMin);

第五步：映射+色彩空间还原

			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(xyz[i][j].y>1)xyz[i][j].y=1;
			xyz[i][j].y = map[(int)(255*xyz[i][j].y+0.5)]/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);

			//转化成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 EqualizeHist(char* fileName, double r)//基于直方图均衡化的色调映射技术
{
	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 NumPixel[256];
	memset(NumPixel,0,sizeof(NumPixel));
	double RatePixel[256];
	memset(RatePixel,0,sizeof(RatePixel));

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

			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=pow(xyz[i][j].y,r);
			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);

			//统计对应色阶的像素点个数
			NumPixel[(int)(xyz[i][j].y*255+0.5)]++;
		}
	}


	double PixelMin;
	int flag=0;
	for(int i=0;i<256;i++)
	{
		//计算分布函数
		if(i==0)RatePixel[i]=1.0*NumPixel[i]/(rows*cols);
		else RatePixel[i]=1.0*NumPixel[i]/(rows*cols)+RatePixel[i-1];
		if(RatePixel[i]!=0&&flag==0)
		{
			flag=1;
			PixelMin = 255.0 * RatePixel[i];
		}
	}

	int map[256]={0};
	for(int i=0;i<256;i++)
	{
		//拉伸对比度
		map[i]=(255.0*RatePixel[i]-PixelMin)*255.0/(255-PixelMin);
	}

	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);
			//if(xyz[i][j].y>1)xyz[i][j].y=1;
			xyz[i][j].y = map[(int)(255*xyz[i][j].y+0.5)]/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);

			//转化成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] = 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;
}