一种Tone Mapping 算法---Adaptive Logarithmic Mapping自适应对数映射

一种Tone Mapping 算法—Adaptive Logarithmic Mapping自适应对数映射微信公众号：幼儿园的学霸在查找图像增强方面的相关算法时,看到一篇博客,而博文中的核心是实现2003年的一篇论文,论文的原意是用在高对比场景下的自适应对数映射算法,但是查看博文中展示的图片效果,发现该算法在低光照下对色彩恢复效果也是很好，论文的名字是《Adaptive Logarithmic M

leonardohaig

2322人浏览 · 2020-12-30 22:34:37

leonardohaig · 2020-12-30 22:34:37 发布

一种Tone Mapping 算法—Adaptive Logarithmic Mapping自适应对数映射

微信公众号：幼儿园的学霸

在查找图像增强方面的相关算法时,看到一篇博客,而博文中的核心是实现2003年的一篇论文,论文的原意是用在高对比场景下的自适应对数映射算法,但是查看博文中展示的图片效果,发现该算法在低光照下对色彩恢复效果也是很好，论文的名字是《Adaptive Logarithmic Mapping For Displaying High Contrast Scenes》。

但是那篇博客中仅仅展示了一小部分实现代码,对于论文中的完整算法实现并未展示.因此特意将论文下载下来,参考博客对算法进行理解与完整复现.

引入

we assumed the logarithmic relation in our tone mapping solution following Stockham who recommended such a relation for image processing purposes:

首先,论文提出一个假设,认为显示器给人眼看到的亮度值 $L_d$ 和场景亮度值/原图像亮度值 $L_w$ 之前存在如下对数映射关系:
公式1
其中 $L_{max}$ 为场景亮度的最大值.

该映射公式能够保证无论场景的动态范围是什么，场景亮度的最大值被重新映射为1（白色），而其他亮度值能够平滑递增.虽然这个公式对一些图像能够获得较好的效果,但是对于有些图的亮度过于压缩了,高对比度的内容也丢失了.

This mapping ensures that whatever the dynamic range of the scene is, the maximum value is remapped to one (white) and other luminance values are smoothly incremented.While this formula leads to pleasant images, we found that the luminance compression is excessive and the feeling of high contrast content is lost.

自适应对数映射

基于上面公式缺点,论文提出了自己的色调映射方案,其设计的方案遵循以下原则：(1)不管外界场景颜色如何丰富多变,算法必须都能输出连续的结果;(2)它应该具有自适应性和可扩展性，能显示出场景的物理本质同时不得引入对比度反转和光晕。总的亮度也必须忠实于实际的内容;(3)算法也需要对用户友好，也就是说在大部分情况下需要能自动实现，少数情况需要调节一些比较直观的参数。

以上原则对应论文中第3部分,接下来论文介绍了: 3.1场景亮度到图像亮度的映射; 3.2对比度调整; 3.3算法. 该部分内容感兴趣的可以看一下,核心在于论文中提到的公式.
论文中的第2个公式是对数的换底公式:

然后是第3个公式:偏置指数函数公式:

这个公的来源在于实验表明人眼对亮度的感受值和实际值符合对数曲线,为了使得对数变换变得"平滑",使用了上面的偏置函数,该函数就是将一个数值t做一个指数变换,来达到调节的目的,当b=0.5时, $bias_{b}(t)=t$ ,当b=0.73时,得到的调整函数最接近 $\gamma=2.2$ 的gamma矫正结果.论文尝试了一系列的b值,最后认为b=0.85时,效果是最好的

接下来,就是论文算法中的核心公式了:

式中多了个参数Ldmax,这个值表示显示设备的最大显示能力，对于普通的CRT显示器，论文直接取值为100。单独观察公式(4)可能看不到什么内容,我们将公式(4)利用换底公式进行一点变换:
公式4展示表示
观察上面展开后的公式:
1)观察标红的对数基(底数),显然,底数部分是变化的,
2)再次观察标红的对数基,很明显,其取值范围是[log(2),log(10)],该取值范围和论文中3.2部分描述一致;

综上,公式中体现了自适应的关键所在:对数基的自适应调整,根据每个像素的信息来自适应的调整函数中的对数基

The principal characteristic of our tone mapping function is an adaptive adjustment of logarithmic base depending on each pixel’s radiance.
对数基自适应:自适应体现在对数基根据亮度值的自适应

最后部分,论文有一点关于gamm矫正内容,其采的公式如下:
gamma矫正公式
其中L表示RGB值, $E^{'}$ 表示矫正后用于显示的像素值;

Where slope is the elevation ratio of the line passing by the origin and tangent to the curve, and start is the abscissa at the point of tangency.

但是具体slope和start的取值是多少,在论文中并没有看到描述.

根据以上介绍和公式(4)就可以写代码以实现该算法了,需要注意的是,作者在映射过程中首先将RGB转换到CIE XYZ色彩空间中进行实现,转换完毕后再次转回到RGB色彩空间.

实现及优化

实现过程中,发现,OpenCV中已经封装有该算法,不过还是自己复现了一把,观察和自带实现效果和效率上的差异.

实现代码如下.在代码中已经有详细的注释,以方便理解

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

#include <chrono>

//Tone Mapping,自适应对数映射
//Input Param:src--输入图像,CV_8UC3,BGR
//            gamma--positive value for gamma correction. Gamma value of 1.0 implies no correction, gamma
//                      equal to 2.2f is suitable for most displays.Generally gamma > 1 brightens the image
//                      and gamma < 1 darkens it.
//            bias--偏置值,value for bias function in [0, 1] range. Values from 0.7 to 0.9
//                          usually give best results, default value is 0.85.
//Output Param:dst--结果图像,type as src
//Return:   null
//Reference:Adaptive Logarithmic Mapping For Displaying High Contrast Scenes
void TMO_Drago03(const cv::Mat& src,cv::Mat& dst,float gamma=2.2f,float bias=0.85f);
//调用OpenCV自带的TMO 算法
void CV_TMO_Drag03(const cv::Mat& src,cv::Mat& dst,float gamma=2.2f,float bias=0.85f);

void TMO_Drago03(const cv::Mat &src, cv::Mat &dst, float gamma, float bias)
{
    assert( CV_8UC3==src.type() );

    const int rows = src.rows;
    const int cols = src.cols;

    //1°BGR色彩空间转换到CIE XYZ色彩空间,并转换到[0,1]区间
    //or 自定义实现,可以采用 位运算符&&并行 以加速转换过程
    std::vector<cv::Mat> xyzChannels;
    cv::cvtColor(src,dst,cv::COLOR_BGR2XYZ);//opencv内部已加速
    dst.convertTo(dst,CV_32FC3,1.0/255);
    cv::split(dst,xyzChannels);

    //1.2° 计算原图像中各点X,Y,Z分量所占比重,以方便后面恢复图像
    cv::Mat xyzSum,xRatio,yRatio,zRatio;
    xyzSum = xyzChannels[0]+xyzChannels[1]+xyzChannels[2]+0.00001f;
    xyzSum = 1.0/xyzSum;
    xRatio = xyzChannels[0].mul(xyzSum);//CV_32FC1
    yRatio = xyzChannels[1].mul(xyzSum)+0.00001f;//CV_32FC1
    zRatio = 1-xRatio-yRatio;//CV_32FC1



    //2°求Lwmax,以及对数和Lwsum,亮度的对数平均值,Lwav;
    double Lwmax(0.0);
    cv::minMaxLoc(xyzChannels[1],nullptr,&Lwmax);
    //cv::Mat _logedY;
    //cv::log(xyzChannels[1]+0.001,_logedY);//对亮度求对数
    //auto Lwsum = cv::sum(_logedY)[0];//求对数和
    //auto Lwav = cv::exp(Lwsum*1.0/(rows*cols));//对数平均值
    //Lwmax /= Lwav;

    //3°自适应对数变换
    //3.1°计算映射表/查找表,以加速求解过程
    float lut_YTable[256]={0};
    const float Ldmax = 100;//一般显示设备该值取100
    auto section1 = Ldmax*0.01/std::log10(Lwmax + 1);//论文公式4第1部分
    auto biase_b = std::log(bias)/std::log(0.5);//论文公式4第2部分分母中的指数
    for(int Y=0; Y!=256;++Y)//tone mapping on each Lum(or Y) pixel's value
    {
        //将[0,255]区间的Y值缩放到[0,1]
        auto oriY = Y*0.0039215;//idx/255.0

        //核心公式:自适应对数变换.对数基的自适应
        //auto _newY = 1.0 * log(oriY + 1) / log(2 + 8.0*pow((oriY / Lwmax), log(bias) / log(0.5))) / log10(Lwmax + 1);
        //将上面公式中重复计算部分提取出来,改写为下面形式
        auto denominator2 = std::log(2 + 8.0*pow((oriY / Lwmax), biase_b));//论文公式4第2部分的分母
        lut_YTable[Y] = section1 * (std::log(oriY + 1) / denominator2);
        //if( lut_YTable[Y]>1)
        //    lut_YTable[Y] = 1.0f;

        //调整饱和度
        //float saturation = 1.0f;
        //lut_YTable[Y] = cv::pow(lut_YTable[Y]/(oriY+0.0001),saturation);
    }


    //3.2°根据查找表计算新的Y;同时根据X,Y,Z的比重恢复X,Z分量
    for(auto iter=xyzChannels[1].begin<float>(); iter != xyzChannels[1].end<float>();++iter)
        *iter = lut_YTable[int(*iter*255)];

    //xyzChannels[0] = xyzChannels[1].mul(xRatio/yRatio);//X X:Y = xratio:yratio----> X = Y*(xratio/yRatio)
    //xyzChannels[2] = xyzChannels[1].mul(zRatio/yRatio);//Z
    //将上面的重复计算部分提取出,更改为下面形式 X = xratio*(Y/yratio)
    auto _yRatio = xyzChannels[1]/yRatio;
    xyzChannels[0] = xRatio.mul(_yRatio);//X
    xyzChannels[2] = zRatio.mul(_yRatio);//Z

    //4° XYZ转换为BGR,并进行gamma校正
    cv::merge(xyzChannels,dst);
    cv::cvtColor(dst,dst,cv::COLOR_XYZ2BGR);

    //gamma 校正函数1
    const float fgamma = (float)((0.45f/gamma) * 2.0f );//即:0.9f/gamma
    auto gamma_correct = [fgamma](float& val) {
        if( val<=0.05)
            val *= 2.64;
        else
            val = 1.099f*cv::pow(val,fgamma)-0.099f;
    };
    //或者gamma校正函数2
    float slope = 4.5f;
    float start = 0.018f;
    if( gamma >= 2.1f )
    {
        start = (float)( 0.018f / ( ( gamma - 2.0f ) * 7.5f ) );
        slope = (float)( 4.5f   * ( ( gamma - 2.0f ) * 7.5f ) );
    }
    else if( gamma <= 1.9f )
    {
        start = (float)( 0.018f * ( ( 2.0f - gamma ) * 7.5f ) );
        slope = (float)( 4.5f   / ( ( 2.0f - gamma ) * 7.5f ) );
    }

    auto rec709_gamma_correct = [fgamma,start,slope](float& val){
        if( val<=start) val *= slope;
        else val = 1.099f*cv::pow(val,fgamma)-0.099f;
    };

    //对各点B,G,R分量分别进行gamma校正//考虑采用查找表的方式,以减少计算量
    //同时,当gamma值确定后,查找表仅用计算一次即可,同样能够减少时间消耗
    for(auto iter=dst.begin<cv::Vec3f>(); iter != dst.end<cv::Vec3f>(); ++iter)
    {
        for(int i=0; i<3; ++i)
        {
            auto& val = (*iter)[i];
            if( val<0) val = 0.0f;
            else if( val>1.0) val=1.0f;

            //gamma_correct(val);
            rec709_gamma_correct(val);
        }
    }

    dst.convertTo(dst,CV_8UC3,255);
}

void CV_TMO_Drag03(const cv::Mat &src, cv::Mat &dst, float gamma, float bias)
{
    cv::Mat ldrDrago;
    src.convertTo(ldrDrago, CV_32FC3, 1.0f/255);
    cvtColor(ldrDrago, ldrDrago, cv::COLOR_BGR2XYZ);

    cv::Ptr<cv::TonemapDrago> ptrTonemapDrago = cv::createTonemapDrago(gamma, 1.f, bias);
    ptrTonemapDrago->process(ldrDrago, dst);

    cvtColor(dst, dst, cv::COLOR_XYZ2BGR);
    dst.convertTo(dst, CV_8UC3, 255);
}

int main()
{
    cv::Mat src = cv::imread("./dark.jpg",cv::IMREAD_COLOR);
    
    const int ITER_TIMES = 500;

    cv::Mat dst1;
    cv::Mat dst2;

    auto time1 = std::chrono::high_resolution_clock::now();

    //自定义实现TMO
    for(int i=0; i<ITER_TIMES; ++i)
        TMO_Drago03(src,dst1);

    auto time2 = std::chrono::high_resolution_clock::now();

    //OpenCV自带TMO
    for(int i=0; i<ITER_TIMES; ++i)
        CV_TMO_Drag03(src,dst2);

    auto time3 = std::chrono::high_resolution_clock::now();


    auto duration1 = std::chrono::duration_cast<std::chrono::milliseconds>(time2-time1).count();
    auto duration2 = std::chrono::duration_cast<std::chrono::milliseconds>(time3-time2).count();
    printf("duration1:%ld ms,duration2:%ld ms.\n",duration1,duration2);


    cv::imshow("src",src);
    cv::imshow("tmo_dst",dst1);
    cv::imshow("cv_tmo_dst",dst2);
    cv::imshow("auto_gamma",dst3);
    
    cv::waitKey(0);
    
    return 0;

}

原图	自定义实现	OpenCV封装实现

附:
对TMO感兴趣的可以搜索以下关键词:
Tone Mapping Operator(TMO)
github搜索topic:tone-mapping

观察发现自定义实现的效果要比OpenCV自带的封装效果好,但是在效率上比自带封装要低,因此原实现的基础上采用并行计算以及gamma校正查找表方式以提高运算效率,对代码稍作修改,修改后的代码如下:

实测发现,修改后的代码效率明显优于自带封装实现,

实际上,该代码还可以继续优化加速,优化建议已写在注释中,不在赘述

void TMO_Drago03_A(const cv::Mat &src, cv::Mat &dst, float gamma, float bias)
{
    assert( CV_8UC3==src.type() );

    const int rows = src.rows;
    const int cols = src.cols;

    //1°BGR色彩空间转换到CIE XYZ色彩空间,并转换到[0,1]区间
    //or 自定义实现,可以采用 位运算符&&并行 以加速转换过程
    std::vector<cv::Mat> xyzChannels;
    cv::cvtColor(src,dst,cv::COLOR_BGR2XYZ);//opencv内部已加速
    dst.convertTo(dst,CV_32FC3,1.0/255);
    cv::split(dst,xyzChannels);

    //1.2° 计算原图像中各点X,Y,Z分量所占比重,以方便后面恢复图像
    cv::Mat xyzSum,xRatio,yRatio,zRatio;
    xyzSum = xyzChannels[0]+xyzChannels[1]+xyzChannels[2]+0.00001f;
    xyzSum = 1.0/xyzSum;
    xRatio = xyzChannels[0].mul(xyzSum);//CV_32FC1
    yRatio = xyzChannels[1].mul(xyzSum)+0.00001f;//CV_32FC1
    zRatio = 1-xRatio-yRatio;//CV_32FC1



    //2°求Lwmax,以及对数和Lwsum,亮度的对数平均值,Lwav;
    double Lwmax(0.0);
    cv::minMaxLoc(xyzChannels[1],nullptr,&Lwmax);
    //cv::Mat _logedY;
    //cv::log(xyzChannels[1]+0.001,_logedY);//对亮度求对数
    //auto Lwsum = cv::sum(_logedY)[0];//求对数和
    //auto Lwav = cv::exp(Lwsum*1.0/(rows*cols));//对数平均值
    //Lwmax /= Lwav;

    //3°自适应对数变换
    //3.1°计算映射表/查找表,以加速求解过程
    float lut_YTable[256]={0};
    const float Ldmax = 100;//一般显示设备该值取100
    auto section1 = Ldmax*0.01/std::log10(Lwmax + 1);//论文公式4第1部分
    auto biase_b = std::log(bias)/std::log(0.5);//论文公式4第2部分分母中的指数
#pragma omp parallel for
    for(int Y=0; Y<256;++Y)//tone mapping on each Lum(or Y) pixel's value
    {
        //将[0,255]区间的Y值缩放到[0,1]
        auto oriY = Y*0.0039215f;//idx/255.0

        //核心公式:自适应对数变换.对数基的自适应
        //auto _newY = 1.0 * log(oriY + 1) / log(2 + 8.0*pow((oriY / Lwmax), log(bias) / log(0.5))) / log10(Lwmax + 1);
        //将上面公式中重复计算部分提取出来,改写为下面形式
        auto denominator2 = std::log(2 + 8.0*pow((oriY / Lwmax), biase_b));//论文公式4第2部分的分母
        lut_YTable[Y] = section1 * (std::log(oriY + 1) / denominator2);
        //if( lut_YTable[Y]>1)
        //    lut_YTable[Y] = 1.0f;

        //调整饱和度
        //float saturation = 1.0f;
        //lut_YTable[Y] = cv::pow(lut_YTable[Y]/(oriY+0.0001),saturation);
    }


    //3.2°根据查找表计算新的Y;同时根据X,Y,Z的比重恢复X,Z分量
#pragma omp parallel for
    for(auto iter=xyzChannels[1].begin<float>(); iter < xyzChannels[1].end<float>();++iter)
        *iter = lut_YTable[int(*iter*255)];

    //xyzChannels[0] = xyzChannels[1].mul(xRatio/yRatio);//X X:Y = xratio:yratio----> X = Y*(xratio/yRatio)
    //xyzChannels[2] = xyzChannels[1].mul(zRatio/yRatio);//Z
    //将上面的重复计算部分提取出,更改为下面形式 X = xratio*(Y/yratio)
    auto _yRatio = xyzChannels[1]/yRatio;
    xyzChannels[0] = xRatio.mul(_yRatio);//X
    xyzChannels[2] = zRatio.mul(_yRatio);//Z


    //4° XYZ转换为BGR,并进行gamma校正
    //4.1°转换为BGR
    cv::merge(xyzChannels,dst);
    cv::cvtColor(dst,dst,cv::COLOR_XYZ2BGR);

    //4.2°计算gamma校正表,以加速求解过程
    //gamma 校正函数1
    const float fgamma = (float)((0.45f/gamma) * 2.0f );//即:0.9f/gamma
    auto gamma_correct = [fgamma](float& val) {
        if( val<=0.05)
            val *= 2.64;
        else
            val = 1.099f*cv::pow(val,fgamma)-0.099f;
    };
    //或者gamma校正函数2
    float slope = 4.5f;
    float start = 0.018f;
    if( gamma >= 2.1f )
    {
        start = (float)( 0.018f / ( ( gamma - 2.0f ) * 7.5f ) );
        slope = (float)( 4.5f   * ( ( gamma - 2.0f ) * 7.5f ) );
    }
    else if( gamma <= 1.9f )
    {
        start = (float)( 0.018f * ( ( 2.0f - gamma ) * 7.5f ) );
        slope = (float)( 4.5f   / ( ( 2.0f - gamma ) * 7.5f ) );
    }

    auto rec709_gamma_correct = [fgamma,start,slope](float& val){
        if( val<=start) val *= slope;
        else val = 1.099f*cv::pow(val,fgamma)-0.099f;
    };

    //对各点B,G,R分量分别进行gamma校正//考虑采用查找表的方式,以减少计算量
    //同时,当gamma值确定后,查找表仅用计算一次即可,同样能够减少时间消耗
    float gamma_Table[256] = {0};
#pragma omp parallel for
    for(int Y=0; Y<256;++Y)//calculate the pixel's gamm corrected value
    {
        //将[0,255]区间的Y值缩放到[0,1]
        float oriY = Y * 0.0039215f;//idx/255.0
        //gamma_correct(oriY);
        rec709_gamma_correct(oriY);
        gamma_Table[Y] = oriY;
    }

    //4.3°利用gamma表进行gamma校正
#pragma omp parallel for
    for(auto iter=dst.begin<cv::Vec3f>(); iter < dst.end<cv::Vec3f>(); ++iter)
    {
        for(int i=0; i<3; ++i)
        {
            auto& val = (*iter)[i];
            if( val<0) val = 0.0f;
            else if( val>1.0) val=1.0f;

            val = gamma_Table[int(val*255)];
        }
    }

    dst.convertTo(dst,CV_8UC3,255);
}