色调映射

色调映射是在有限动态范围媒介上近似显示高动态范围图像的一项计算机图形学技术。打印结果、CRT 或者 LCD 显示器以及投影仪等都只有有限的动态范围。

本质上来讲，色调映射是要解决的问题是进行大幅度的对比度衰减将场景亮度变换到可以显示的范围，同时要保持图像细节与颜色等对于表现原始场景非常重要的讯息。

根据应用的不同，色调映射的目标可以有不同的表述。在有些场合，生成“好看”的图像是主要目的，而在其它一些场合可能会强调生成尽可能多的细节或者最大的图像对比度。在实际的渲染应用中可能是要在真实场景与显示图像中达到匹配，尽管显示设备可能并不能够显示整个的亮度范围。

<image>

在最近几年中已经开发了各种各样的色调映射算法，其中一个简单的色调映射滤波器例子是 $L=Y/(Y+1)$，这个函数将场景在 $[0,\infty)$ 域内的 radiance 值 $Y$ 映射到显示输出范围 $[0,1)$。

<image>

另外一组更加复杂的算法是基于对比度或者梯度域的方法，这些算法的侧重点在于对比度的保持而不是亮度的映射，这种方法的思路起源于人眼对于对比度或者不同亮度区域的亮度比例最为敏感。这种色调映射由于较好地保存了对比度细节，所以通常会产生非常锐利的图像，但是这样做的代价是使得整体的图像对比度变得平缓。这种色调映射方法的例子包括：梯度域高动态范围压缩以及高动态范围图像感知框架（有一种色调映射就是这种框架的一个应用）。

一个有趣的高动态范围图像色调映射实现方法是从anchoring theory of lightness perception 得到的灵感。这个理论全面地解释了人类视觉系统中如 lightness constancy 以及 its spectacular failures 这些对于感知图像非常重要部分的特点。这个色调映射方法（色调重建中的亮度感知）的核心概念是将高动态范围图像分解成照明一致的区域或结构以及计算局部亮度值。根据各个区域亮度的比例进行合并，从而计算图像的纯亮度值。其中尤其重要的是亮度定位（anchoring）——将照明的亮度对应于已知的亮度值，也就是说哪个亮度值在场景中感知为白色。这种色调映射实现方法不影响局部对比度，并且由于对亮度进行线性处理所以也保留了高动态图像的自然颜色。

## Photographic 算法

在此算法中，主要分为两个步骤来改变场景中的亮度，第一个步骤为整个场景的亮度改变(global)，就如同在拍摄景物的时候需要调整相机的曝光值；第2个步骤为局部区块的亮度改变(local)，这个步骤使用到了局部加亮和加深(Dodging & Burning)的概念

### 步骤一 整个场景的亮度改变

首先，根据照片红、绿、蓝3个颜色可以计算出照片中每个像素的亮度$L_w(x,y)$，再来需要算出场景中的关键值(key value)，在这里关键值可以代表一个场景整体来说是亮的还是暗的，如果照片是白天拍摄的，就可以说其关键值是高的，若是在晚上拍摄的，则关键值会较低。在这里我们可以假设关键值会与场景中亮度的对数平均值相关。

我们以$\overline{L_w}$来估计场景的关键值。

$\overline{L_w} =\exp{(\frac{(\sum_{x,y}\log{(\delta+L_w(x,y))}}{N})}$

$\delta$: 一个极小值，避免在算黑色像素的时候会遇到$\log = 0$的情况算出场景中的关键值后，我们要将此关键值重新调整到中间值$\alpha$，调整方法为$L(x,y) = \frac{\alpha}{\overline{L_w}}L_w(x,y)$

其中$\alpha$为调整后新的关键值。

调整完后，会发现此方法尚无法解决场景过亮的地方会超出所能显示的颜色范围的问题。所以下一步是要将亮度高的地方压缩。

$L_d(x,y) = \frac{L(x,y)(1+\frac{L(x,y)}{{L_{white}}^2})}{1+L(x,y)}$

其中$L_{white}$场景中亮度最高的值。

由上述式子可以发现亮度接近最大亮度的像素会被压缩成$\frac{1}{L}$，而亮度较低的地方压缩倍数则较少，如此即可保留低对比区域的细节，并且将高亮度的区域压缩至合理范围。

### 步骤二 局部区块的亮度

在使用完simple operator后，对于动态范围很高的影像，可以发现有些景物的细节还是无法被保留，所以此算法又提出了第二个步骤- 局部加亮与加深(Dodging & Burning)的作法。此作法是参考早期冲洗相片时的一个常用技巧。

早期在冲洗底片的时候，为了能让照片中的某些区域曝光高一点，有些区域低一点，所以使用加亮和加深(Dodging & Burning)的技巧，在冲洗照片的过程中，将一些照片中的区域用纸张或是木头遮挡住，来调整不同区域的亮暗，此作法可以提升原本暗处的亮度或是让过曝的地方能够变暗一点，以看清楚照片中的细节。

通常加亮与加深是使用在一个被高对比度包围的影像区域上，所以此方法首先要先找到对于每个像素而言，不存在强烈对比的最大区域。

找寻方法为，将两个不同变异数的高斯函数套在影像上，得出两个结果$V_1(x,y,s)$ 和 $V_2(x,y,s)$，$V_1, V_2$代表了给定某一个像素下，以此像素为圆心，在一个半径内的亮度平均值，而计算$V_2$所使用的高斯函数的变异数较大，所以得出的结果为在一个较大半径内的平均值$V_1, V_2$算法如下：

$V_i(x,y,s) = L(x,y,)\otimes R_i(x,y,s)$ $R_i(x,y,s) = \frac{1}{\pi(\alpha_i s)^2}\exp{-\frac{x^2+y^2}{(\alpha_i s)^2}}$ 并以函数$V(x,y,s)$来凭断套用这两个不同高斯后的差别。当所计算区域内的亮度变化很小的时候，$V_1(x,y,s), V_2(x,y,s)$的差异很小，$V(x,y,s)$也就很小；当所计算区域内的亮度有很明显差异的时候，通常也就是位于物体边界的时候，$V(x,y,s)$会很大。

$V(x,y,s) = \frac{V_1(x,y,s)-V_2(x,y,s)}{2^\phi a/s^2 + V_1(x,y,s)}$ 其中$s$代表scale，此步骤会重复多次在不同的scale上，每次都会选用更大的变异数值。从最小的scale开始找起，目标是要找到$V(x,y,s)$最接近临界值，但是又小于临界值的scale。

最后的亮度可以表示为$L_d(x,y) = \frac{L(x,y)}{1+V_1(x,y,s_m(x,y))}$ 其中$V_1(x,y,s_m(x,y))$为满足$V(x,y,s)$小于临界值下的$V_1$。

值得注意的是，在使用此算法的时候，如果选择的scale太大，也就是说每次变异数都变化很大的话或是临界值选择太大，都有可能会造成无法选择到正确的亮度无变化的最大区域，导致算法估计错误，最后影像会有明显黑边。

## 实作Photographic 算法

以下为使用Matlab所实作的Photographic 算法之程式码

```
function rgb = tonemap_p( hdr )
    % HDR need be double format 
    %步驟1,調整整個場景的亮度  
    a = 0.6;%關鍵值設定 
    epsilon = 0.05; 
 
    %利用rgb計算照片亮度 
    lum_w = 0.27 * hdr(:,:,1) + 0.67 * hdr(:,:,2) + 0.06 * hdr(:,:,3); 
    image_size = size( hdr );
    [height, width, channel] = size(hdr)
    %計算L_w_bar(x,y) 
    sum_all = sum(sum(log(0.00000001+ lum_w)));

    N = height * width;
    lum_white = 1e20; 
    lum_w_bar = double( exp (sum_all / N) ); 
 
    %計算L_d(x,y) 
    lum = a / lum_w_bar * lum_w; 
    
    lum_add1 = 1 + lum;
    luma_d_1 = lum .* (1 + lum ./ lum_white ./ lum_white) ./lum_add1;
    
    %步驟二：加亮與加深(Dodging & Burning) 
    phi = 8; 
    %計算兩個不同變異數大小的高斯函數 
    for i = 1: height 
        for j = 1: width 
            for k = 1: 8 
                s = 1 * 1.6 ^ k; 
                alpha1 = 0.35; 
                R1(i, j, k) = exp( -( i ^ 2 + j ^ 2) / ( alpha1 ) ^ 2 ) / ( pi * ( alpha1 * s ) ^ 2 ); 
                alpha2 = 0.35 * 1.6; 
                R2(i, j, k) = exp( -( i ^ 2 + j ^ 2) / ( alpha2 ) ^ 2 ) / ( pi * ( alpha2 * s ) ^ 2 ); 
            end 
        end 
    end 
    %將亮度與高斯函數做卷積 
    for k = 1: 8 
        V1(:,:,k) = conv2( lum_w, R1(:,:,k)); 
        V2(:,:,k) = conv2( lum_w, R2(:,:,k)); 
    end 
 
    %計算函數V(x, y, s) 
    for i = 1: height 
        for j = 1: width 
            for k = 1:8 
                V(i, j, k) = ( V1(i, j, k) - V2(i, j, k)) / (2 ^ phi * a / (k ^ 2) + V1(i, j, k)); 
            end 
        end 
    end 
 
    %找尋符合V小於臨界值epsilon的scale 
    for i = 1: height 
        for j = 1: width 
            for k = 1:8 
                if abs(V(i, j, k)) < epsilon 
                    V1_new(i, j) = V(i, j, k); 
                end 
            end 
        end 
    end 
 
    %計算新的亮度 
    for i = 1: height 
        for j = 1: width 
            lum_d2(i, j) = lum_w(i, j) / (1 + V1_new(i, j)); 
        end 
    end 
 
    %將rgb影像轉到hsv的系統後,改變亮度大小,最後輸出成色調映射後的結果 
    hsv = rgb2hsv( hdr ); 
    new_img = cat(3,hsv(:,:,1),hsv(:,:,2),lum_d2);    
    rgb = hsv2rgb( new_img); 
end
 
```

## 图库

<image>

- 新西兰惠灵顿 圣保罗教堂 大堂，从 1/20 秒到 30 秒的曝光的八幅照片

<image>

- 两幅四级曝光差别的照片组合在一起，然后动态范围进行压缩以在标准显示器上显示

<image>

- 高动态范围落日图像。Taken as an exposure bracket of three images 2 stops apart.

## 参考文献

- devlin Kate Devlin, Alan Chalmers, Alexander Wilkie, Werner Purgathofer. "STAR Report on Tone Reproduction and Physically Based Spectral Rendering" in Eurographics 2002. DOI: 10.1145/1073204.1073242
- fattal Raanan Fattal, Dani Lischinski, Michael Werman. "梯度域高动态范围压缩"
- mantiuk Rafal Mantiuk, Karol Myszkowski, Hans-Peter Seidel. "高动态范围图像对比度处理的感知架构"
- gilchrist Alan Gilchrist. "An Anchoring Theory of Lightness Perception".
- krawczyk Grzegorz Krawczyk, Karol Myszkowski, Hans-Peter Seidel. "高动态范围图像色调重现中的亮度感知"
- reinhard Reinhard, Erik, et al. "Photographic tone reproduction for digital images." ACM transactions on graphics (TOG) 21.3 (2002): 267-276.

## 外部链接

- pfstmo: 色调映射算子的实现
- exrtools: OpenEXR 图像处理工具集合 （包括一些色调映射算子）