【ISP图像处理】Tone Mapping基础知识及相关算法（附代码）

蘑菇桑巴

已于 2023-10-26 11:25:05 修改

阅读量2.8k

点赞数 8

分类专栏：图像处理文章标签： 1024程序员节图像处理 isp

于 2023-10-24 22:07:53 首次发布

本文链接：https://blog.csdn.net/royole98/article/details/134023035

版权

图像处理专栏收录该内容

8 篇文章 9 订阅

订阅专栏

1、基本定义

一般LDR（Low Dynamic Range）图像的颜色显示范围通常只有8位，即每个颜色通道的颜色数值有2^8=256个等级。这个量级用于描述现实场景中的景象往往十分有限，以LDR储存图像往往需要对颜色进行压缩。为了更加真实还原真实场景的颜色，HDR图像应运而生，一般通道位数超过8位，便可称为HDR，常见有12位和16位。

虽然存储图像的信息量提升了，但是现在使用的大部分显示设备宽动态范围只有100:1甚至更低。为了使得HDR图像能够在低动态范围的显示设备上显示，Tone Mapping技术便十分重要。它可以将HDR图像颜色范围进行压缩，且这种压缩并不是简单的线性压缩，它需要保存原图的颜色，并且把高亮区域和阴影区域向中等亮度方向压缩。这一变化如图所示。

2、相关算法

2.1 Global Tone Mapping(全局色调映射，GTM)

GTM将每个像素值映射到另一个值，而不考虑像素的位置。这种方式因为简单有效，在早期被大量使用。

（1）Photographic Tone Reproduction for Digital Images [PDF]

首先，文章作者对整个图像做了一个光亮度的映射，其作用类似于设置相机的曝光。这一操作实际就是对每一个像素做固定缩放。根据前人对Tone Mapping的研究结论，论文作者认为光亮度的log平均值能够反映图片中像素光亮度的特征。因此，作者用该值对每个像素作缩放。

缩放后的光亮度 L(x，y) 可用如下公式表示：

其中，α 是一个缩放参数，被称为Key Value，不同的 α 值对应了不同的缩放程度，如下图所示，α 值越小，画面整体色调越偏灰偏暗：

但是，仅仅做简单线性缩放是不够的。对于光亮度变化不是很大的图片，这种方法可以将像素的光亮度很好地压缩到一定范围，但是对于大多数图片，绝大部分像素光亮度是在某一个范围之内，而少数高亮的像素比平均值高太多，很容易产生过曝，如：光源、高光反射等。因此，通常在经过线性缩放之后，还需要利用非线性的算子对图像进行处理。其中最简单的算子是：

可以看到，当亮度越小，缩放因子接近于1；而亮度越大，则缩放因子趋近于1/L。这样对高亮区域有着更强的压缩效果。

公式（3）的算子可以进一步被拓展为：

其中，Ld（x，y）表示经过非线性算子处理后的像素，Lwhite表示图片中被映射到白色亮度的像素中的最小值。这样可以使得最亮像素点被映射到1，以减少信息的损失。

以下是算法对应C++代码，注意算子是应用在图像的亮度通道。

float luminance(vec3 v)
{
    return dot(v, vec3(0.2126f, 0.7152f, 0.0722f));
}

vec3 change_luminance(vec3 c_in, float l_out)
{
    float l_in = luminance(c_in);
    return c_in * (l_out / l_in);
}

vec3 reinhard_extended_luminance(vec3 v, float max_white_l)
{
    float l_old = luminance(v);
    float numerator = l_old * (1.0f + (l_old / (max_white_l * max_white_l)));
    float l_new = numerator / (1.0f + l_old);
    return change_luminance(v, l_new);
}

（2）Filmic Tonemapping with Piecewise Power Curves [Blog]

在Tone Mapping中，还有一种最直观的映射曲线，即“S”曲线。如图所示：x曲线是输入高动态范围像素（一般会先转换到log域），y是映射后输出的低动态范围像素，高亮暗部像素都向中间压缩，即实现抑制过曝提亮暗部的目的。

Uncharted 2人工计算一条S曲线进行拟合，得到一个多项式曲线的表达式。这样的表达式应用到渲染结果后，就能在很大程度上自动接近人工调整的结果。

代码如下：

vec3 uncharted2_tonemap_partial(vec3 x)
{
    float A = 0.15f;
    float B = 0.50f;
    float C = 0.10f;
    float D = 0.20f;
    float E = 0.02f;
    float F = 0.30f;
    return ((x*(A*x+C*B)+D*E)/(x*(A*x+B)+D*F))-E/F;
}

vec3 uncharted2_filmic(vec3 v)
{
    float exposure_bias = 2.0f;
    vec3 curr = uncharted2_tonemap_partial(v * exposure_bias);

    vec3 W = vec3(11.2f);
    vec3 white_scale = vec3(1.0f) / uncharted2_tonemap_partial(W);
    return curr * white_scale;
}

（3）ACES Filmic Tone Mapping Curve [Blog]

类似Uncharted 2使用一个多项式进行拟合，ACES提出的多项式则更加简单，S曲线如图所示：

可以看到，亮部和暗部区域像素被压缩范围更加均衡，应用中色彩对比度还原会更加优异。

代码如下：

float3 ACESToneMapping(float3 color, float adapted_lum)
{
	const float A = 2.51f;
	const float B = 0.03f;
	const float C = 2.43f;
	const float D = 0.59f;
	const float E = 0.14f;

	color *= adapted_lum;
	return (color * (A * color + B)) / (color * (C * color + D) + E);
}

方法效果可以参考另一篇博文。

（4）直方图均衡 [PDF]

图像直方图描述了一张图片像素值分布的情况，我们知道，例如对于一张sRGB图像，我们可以看作是一个H*W*1的矩阵，矩阵中每个元素的取值范围为[0,255],图像直方图h(g)定义为像素值为g的元素的个数。如图所示，原始HDR图像大部分像素都是偏暗，因此其对应直方图像素值也集中分布在值较小的区间。经过直方图均衡后，图像中个数多的灰度级进行展宽，而对像素个数少的灰度级进行缩减，从而达到提升暗部，抑制高亮区域，清晰图像的目的。

实现代码如下：

using namespace cv;
using namespace std;

int main()
{ 
	Mat image, image_gray, image_enhanced;   //定义输入图像，灰度图像, 直方图
	image = imread("lena.png");  //读取图像；
	if (image.empty())
	{
		cout << "读取错误" << endl;
		return -1;
	}

	cvtColor(image, image_gray, COLOR_BGR2GRAY);  //灰度化
	imshow(" image_gray", image_gray);   //显示灰度图像

	equalizeHist(image_gray, image_enhanced);//直方图均衡化
	imshow(" image_enhanced", image_enhanced);   //显示增强图像

	waitKey(0);  //暂停，保持图像显示，等待按键结束
	return 0;
}

2.2 Local Tone Mapping(局部色调映射，LTM)

在真实场景中，仅靠GTM增强图像会产生局部区域过度/不足增强。LTM提供更细粒度的控制，并有助于突出重点区域细节。

（1）Photographic Tone Reproduction for Digital Images [PDF]

这篇论文除了提出了一个GTM算子，同时还提出了一个LTM算子。作者认为仅仅依赖GTM算子并不足以应对动态范围极广的图像，比如，强光反射区域中的小面积纹理细节容易在被压缩过程中丢失。因此，作者进一步结合dodging-and-burning方法，提出了一个LTM算子，希望做颜色映射时进一步保持对比度。

首先介绍dodging-and-burning方法，该方法的核心思想是对于每一个像素点都考察不同尺度的邻域，以找到使得对比度没有明显突变的最大邻域。如下面示意图所示，尺度过小和过大时，中心区域和环绕区域间的对比度都比较低，而在正确的尺度上，中心区域与环绕区域的对比度就比较大。不同尺度的运用其实是通过使用不同尺度的高斯滤波器对图像滤波来实现的，中心区域与环绕区域也是半径具有一定比例的两个邻域。记中心区域的尺度为s1,环绕区域的对比度为s2，对于每个像素点，比较两个尺度下的滤波结果V1和V2，设置一个可以利用正确尺度下V1和V2的高对比度性质的公式进行尺度甄选即可。

具体计算过程，首先采用高斯滤波器对图像进行卷积：

其中， $\alpha_i$ 的 $i$ 有1和2两个取值，分别对应中心区域与环绕区域的尺度缩放系数，论文中取值为 $\alpha_1=1/2\sqrt{2}$ , $\alpha_2=1.6\alpha_1$ 。在获得了每个尺度的 $V_1$ 和 $V_2$ 后，根据中心环绕函数计算一个对比度指标：

其中 $\phi$ 默认取值为8.0。分母中前面这一项是为了防止V1趋于0时V取值过大。

如果V1和V2的值接近，则说明这一片区域没有变化强烈的的内容，即对比度低，计算出的V值也趋近于0；反之则V越大。考虑局部区域原则上是给定像素周围没有大对比度变化的最大区域，所以我们需要找到一个大下适中的滤波核，不会因半径过小没有领域信息卷积后变化不明显，也不会因半径过大在高亮区域周围形成暗环。

通过设置阈值，第一个对比度指标大于阈值的小半径滤波核尺度即为所需的卷积半径：

$\left | V(x,y,s_m) \right |> \varepsilon$

确定该高斯滤波核之后，映射公式可表达为：

由公式可知，在相对明亮的区域，暗像素的亮度会满足L>V1，因此该算子会降低显示亮度Ld，从而增加该像素处的对比度，同样的，一个像素在一个相对较暗的区域将被压缩较少，像素相对于周围区域的对比度增加。

python代码如下：

def reinhard_local(img,eps,refV):
    gamma = 1 / 1.6
    picNum = 2
    ratio = 2.7
    # refV = 1.0
    phi = 8.0   # 8.0
    # eps = 0.05
    alpha = 1 / np.sqrt(8)
    epsilon = 0.00001
    a = 0.54

    # L = 0.27 * img[:, :, 0] + 0.67*  img[:, :, 1] + 0.06 * img[:, :, 2] + epsilon
    L = 0.299 * img[:, :, 0] + 0.587 * img[:, :, 1] + 0.114 * img[:, :, 2] + epsilon
    R = img[:, :, 0] / L
    G = img[:, :, 1] / L
    B = img[:, :, 2] / L

    L_log = np.log(L + 1e-9)
    temp = np.exp(L_log.mean())
    L = L * a / temp
    H, W = np.shape(L)

    V1 = np.zeros((W, H, picNum + 1))

    for i in range(picNum + 1):
        s = ratio**i
        sigma = alpha * s
        kernelRadius = int(np.ceil(2 * sigma))
        kernelSize = 2 * kernelRadius + 1
        kernel_H = np.multiply(cv2.getGaussianKernel(kernelSize, sigma), (cv2.getGaussianKernel(1, sigma)).T) 
        temp = tmu.conv2(L, kernel_H, 'same')
        kernel_V = np.multiply(cv2.getGaussianKernel(1, sigma), (cv2.getGaussianKernel(kernelSize, sigma)).T)
        V1[:,:,i] = tmu.conv2(temp, kernel_V, 'same')
    
    V = np.zeros((W, H, picNum))
    for i in range(picNum):
        s = ratio ** i
        V[:,:,i]=abs(V1[:,:,i+1]-V1[:,:,i]) / (a * (2 ** phi)/(s**2) + V1[:,:,i])
    
    Vsm = V1[:,:,picNum]
    for i in range(W):
        for j in range(H):
            for k in range(picNum):
                if V[i,j,k] > eps:
                    Vsm[i,j]=V1[i,j,k] 
                    break
    
    Ld = 1 * L / (refV + 0.6*Vsm)
    Ld[Ld > 1] = 1  #亮度通道重映射结果

    out = np.zeros(img.shape)

    out[:,:,0] = R * Ld
    out[:,:,1] = G * Ld
    out[:,:,2] = B * Ld

    return out

matlab代码可以参考这篇博文。

（2）基于快速双边滤波的方法 [PDF]

前面说到，如果直接整张图进行压缩，容易使得图像细节丢失。为了保持图像细节，提出基于快速双边滤波方法。该方法的核心思想是对图片进行分层：主要为亮度变化信息的基础层（base），和包含高频纹理的细节层（detail）。所以在压缩的时候，只需要对基础层做对比度压缩，压缩之后的基础层再和原来的细节层相加，得到保留了细节信息的低动态范围图像。

为了很好地将低频基础层和高频细节层分离，文章采用双边滤波算子与原图做卷积，卷积后的图像作为基础层，原图再减去基础层作为细节层。双边滤波的原理可以参看这篇博文。

基础层的压缩则通过一个压缩因子实现。首先计算压缩因子：

然后将该压缩因子和基础层相乘则实现压缩过程。需要注意的是，整个过程都是在log域完成的。最终图像获取见如下公式：

完整的流程可总结为：

实现的python代码如下：

def durandanddorsy(img, c):
    height, width = np.shape(img)[:2]
    img = img/img.max()
    epsilon = 0.000001
    # L = 0.299 * img[:, :, 0] + 0.587 * img[:, :, 1] + 0.114 * img[:, :, 2] + epsilon #Y - grey
    L = 1 / 61 * (20 * img[:, :, 0] + 40 * img[:, :, 1] + img[:, :, 2] + epsilon) 
    R = img[:, :, 0] / L
    G = img[:, :, 1] / L
    B = img[:, :, 2] / L

    L_log = np.log(L) 

    base_log = cv2.bilateralFilter(L_log, 9, 0.3, 0.3)
    detail_log = L_log - base_log
    compressionFactor = np.log(c) / (base_log.max() - base_log.min())
    log_abs_scale = base_log.max() * compressionFactor
    # Ld_log = base_log * compressionFactor + detail_log
    Ld_log = base_log * compressionFactor + detail_log - log_abs_scale
    
    out = np.zeros(img.shape)
    out[:,:,0] = R * np.exp(Ld_log)
    out[:,:,1] = G * np.exp(Ld_log)
    out[:,:,2] = B * np.exp(Ld_log)
    
    outt = np.zeros(out.shape, dtype=np.float32)
    cv2.normalize(out, outt, alpha=0, beta=1, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)

    return outt

c++代码可以参考这篇博文。

（3）Adaptive Local Tone Mapping Based on Retinex for High Dynamic Range Images [PDF]

这篇文章做LTM的思路和Reinhard很像，也是先做全局的色调映射，然后再做局部色调映射。首先介绍GTM算子：

其中N是图像的像素总数， $\delta$ 是一个很小的值，为了防止遇到图像中亮度为0的情况；

通过GTM将HDR图像颜色范围压缩到[0, 1]之间后，再进一步应用LTM算子。这篇论文也提到了，在做完GTM后，亮部区域周围可能出现黑边情况，这一问题可以在局部映射中引入保边滤波器解决。

在介绍LTM之前，先介绍Reninex原理。其核心就是类似于人视网膜的视觉观察，我们看到的物体有一部分信息是入射光影响造成的，我们要通过我们看到的图像推算出入射光的影响所造成的图像，然后求二者的差，就是真实的反射物体该有的样子。

基于该原理，Retinex输出可以通过如下公式计算获得：

其中，这篇文章使用引导滤波（原理介绍可见博文）作为F计算获得观测后的图像。我们需要的是实际的反射图像R。

为了进一步防止滤波器造成的平面外观，文章引入了两个因子： $\alpha$ 用于提升图像对比度， $\beta$ 作为自适应非线性偏移，通过调整对数起点来控制非线性强度。因子计算和LTM过程可用如下公式表达（注意公式（13）把log相减变为了除法，实际还是Retinex公式）：

实现的python代码如下：

def maxfilt2(X, size=(3, 3), mode='constant'):
    # filtering
    Y = maximum_filter(X, size)

    return Y

def guidedfilter(I, p, r, eps):
    """
    I: guided image
    p: input image
    r: radius
    eps: regularization
    """

    #均值滤波
    def boxfilter(imSrc, r):
        h, w = np.shape(imSrc)
        imDst = np.zeros((h, w))

        imCum = np.cumsum(imSrc, 0)
        imDst[:r+1, :] = imCum[r:2*r+1, :]
        imDst[r+1:(h-r), :] = imCum[2*r+1:h, :] - imCum[:(h-2*r-1), :]
        imDst[(h-r):h, :] = np.tile(imCum[h-1, :], [r, 1]) - imCum[h-2*r-1:h-r-1, :]

        imCum = np.cumsum(imDst, 1)
        imDst[:, :r+1] = imCum[:, r:2*r+1]
        imDst[:, r+1:w-r] = imCum[:, 2*r+1:w] - imCum[:, :w-2*r-1]
        imDst[:, w-r:w] = np.tile(imCum[:, w-1].reshape(len(imCum[:, w-1]),1), (r)) - imCum[:, w-2*r-1:w-r-1]
        return imDst

    h, w = np.shape(I)
    N = boxfilter(np.ones((h, w)), r)

    mean_I = boxfilter(I, r) / N
    mean_p = boxfilter(p, r) / N
    mean_Ip = boxfilter(I*p, r) / N
    cov_Ip = mean_Ip - mean_I * mean_p

    mean_II = boxfilter(I*I, r) / N
    var_I = mean_II - mean_I * mean_I

    a = cov_Ip / (var_I + eps)
    b = mean_p - a * mean_I

    mean_a = boxfilter(a, r) / N
    mean_b = boxfilter(b, r) / N

    q = mean_a * I + mean_b
    return q

def ALTM_local(img, krnlsz, lamda):
    h, w = img.shape[:2]
    img = np.float32(img)  # res = np.array(img, dtype=np.float32)  # 转换为32位图像
    Lw = 0.299 * img[:, :, 0] + 0.587 * img[:, :, 1] + 0.114 * img[:, :, 2]
    R = img[:, :, 0] / Lw
    G = img[:, :, 1] / Lw
    B = img[:, :, 2] / Lw
    
    Lwmax = Lw.max()
    Lwaver= np.exp(np.log(0.00001 + Lw).sum() / (h * w))
    Lg = np.log(Lw / Lwaver + 1) / np.log(Lwmax / Lwaver + 1)

    #local adaptation
    kenlRatio = 0.1
    # krnlsz = 10
    # krnlsz = np.floor(max([3, h * kenlRatio, w * kenlRatio]))
    Lg1 = maxfilt2(Lg, (krnlsz, krnlsz))
    Hg = guidedfilter(Lg, Lg1, 10, 0.0001)

    eta = 36
    alpha = 1 + eta * Lg / Lg.max()
    Lgaver = np.exp( np.log(0.00001 + Lg).sum() / (h * w))
    beta = 10 * Lgaver

    Lout = alpha * np.log(Lg / Hg + beta)
    Lout = (Lout - Lout.min()) / (Lout.max() - Lout.min())

    Ld = Lout
    r_altm = R * Ld
    g_altm = G * Ld
    b_altm = B * Ld
    outt = cv2.merge([r_altm,g_altm,b_altm])
    
    return outt

2.3 其他

由于笔者工作场景中硬件设备限制，能使用的tone mapping方法计算复杂度不能太高，所以以上整理基本为传统方法。在收集资料过程中，也有许多基于深度学习的方法，虽不能应用，但是一并整理在此，也希望可以给做相关任务的朋友一些启发。

（1）A comparative review of tone-mapping algorithms for high dynamic range video

（2）Real-Time Tone Mapping: A Survey and Cross-Implementation Hardware Benchmark

（3）Which Tone-Mapping Operator Is the Best? A Comparative Study of Perceptual Quality

（4）Comprehensive Fast Tone Mapping for High Dynamic Range Image Visualization

（5）Globally Optimized Linear Windowed Tone-Mapping

（6）Joint Multi-Scale Tone Mapping and Denoising for HDR Image Enhancement

（7）High Dynamic Range Image Tone Mapping: Literature review and performance benchmark

参考文献：

（1）基于物理的渲染—HDR Tone Mapping - 知乎

（2）Tone Mapping | δelta

（3）Tone mapping进化论 - 知乎

（4）HDR技术讲解之Photographic Tone Reproduction for Digital Images论文算法复现及改进-CSDN博客

（5）图像处理基础（十五）基于双边滤波的 tone-mapping - 知乎

蘑菇桑巴

关注

8
点赞
踩
26

收藏

觉得还不错? 一键收藏
4
评论
【ISP图像处理】Tone Mapping基础知识及相关算法（附代码）

一般LDR（Low Dynamic Range）图像的颜色显示范围通常只有8位，即每个颜色通道的颜色数值有2^8=256个等级。这个量级用于描述现实场景中的景象往往十分有限，以LDR储存图像往往需要对颜色进行压缩。为了更加真实还原真实场景的颜色，HDR图像应运而生，一般通道位数超过8位，便可称为HDR，常见有12位和16位。虽然存储图像的信息量提升了，但是现在使用的大部分显示设备宽动态范围只有100:1甚至更低。为了使得HDR图像能够在低动态范围的显示设备上显示，Tone Mapping技术便十分重要。
复制链接

扫一扫