色調映射

色調映射是在有限動態範圍媒介上近似顯示高動態範圍圖像的一項電腦圖形學技術。列印結果、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 圖像處理工具集合 （包括一些色調映射算子）