Camera response curve calibration

公式參考引用數驚人的[1]”Paul E. Debevec and Jitendra Malik. Recovering High Dynamic Range Radiance Maps from Photographs. In SIGGRAPH 97, August 1997.”

1. Introduction

在影像處理中，有幾個主題如HDR, Lens shading correction, shape from shading等，都需精準的由intensity判斷一些重要訊息，然而在真實的相機中，input irradiance \(E\)和輸出的intensity \(Z\)未必是與曝光時間\(\Delta{t}\)呈線性關係，而是一個類似gamma curve的非線性的對應曲線。這會造成計算上的誤差，甚至是許多artifact的產生，因此校正這個關係使得整張影像的intensity \(Z\)與input irradiance \(E\)呈線性變成一個研究主題。

• Remark: 這部分是我自己的碎碎念

 

2. Response Curve

如下圖，橫座標代表進入感光元件的能量，縱坐標則表示顯示出的intensity亮度；圖中不同曲線則代表不同相機的response function (柯達的眼淚)，由於感光元件多以Bayer pattern排列，圖中多是取G channel做的結果。 我們發現大多intensity的反應比我們預期的更亮許多。

[2] “Modeling the Space of Camera Response Functions” IEEE Trans Pattern Anal Mach Intell. 2004

3. Model

我們令這個對應關係為一個response function \(f\)

\[\begin{equation} \label{eq:model}
Z=f(E\Delta{t})
\end{equation}\]

其中\(E\Delta{t}\)代表一個場景irradiance \(E\)經過\(\Delta{t}\)時間取樣的總能量，經過相機的response function轉變為我們讀取到的intensity \(Z\)

並且\(f\)被定義為有幾個重要特性：

• monotonic: 亮的輸入就亮的輸出，關係不會錯置
• smooth: 畢竟取樣總是不完美的，讀者可以嘗試不加上smooth regularization看看\(f\)是不是很noisy
• [2]中利用了PCA定義出\(f\)幾個重要的basis function如下圖

log basis function[2]

4. Optimization

我們的目標就是反求出\(f\)和真實場景的irradiance \(E\)。

\[\begin{equation}
f^{-1}{(Z)}=E\Delta{t}
\end{equation}\]

• Remark: 其中\(E\)就是所謂的High Dynamic Range (HDR) image，經過normalize和tone mapping之後成我們在螢幕看到值域為0~255的影像。這裡和攝影常用的HDR名詞與CG中的HDR名詞有點混用，。

左右各取\(\ln\)將乘法拿掉方便最佳化

\[\ln{f^{-1}{(Z)}}=\ln{E}+\ln{\Delta{t}}\]

令等號左邊\(g(x)=\ln{f^{-1}{(x)}}\)

\[\begin{equation} \label{eq:1}
g(Z)=\ln{E}+\ln{\Delta{t}}
\end{equation}\]

4-1. Quadratic form

在校正時，\(Z\)和曝光時間\(\Delta{t}\)為已知，\(g\)和\(E\)未知。我們對同一個場景取多張不同曝光時間的影像來校正。令\(i\)為pixel位置，\(j\)表不同曝光的影像，\(Z_{ij}\)為在\(j\)影像下pixel \(i\)的intensity，\(E_{i}\)為場景在pixel \(i\)位置的irradiance，\(\Delta{t_{j}}\)表第\(j\)影像的曝光時間。

16張不同曝光時間的影像[1]

我們嘗試最佳化所有取樣pixel符合上式\(\eqref{eq:1}\)

\[\begin{equation}
\min_{g, E}{\sum_{i=1}^{N}{\sum_{j=1}^{P}{{\left[g({Z}_{ij})-\ln{E_{i}}-\ln{\Delta{t_{j}}}\right]}^{2}}}}
\end{equation}\]

再加上smooth regularization，\(Z_{min}, Z_{max}\)通常是0~255 (8-bit intensity)，\(\lambda\)是一個使用者調整的參數，控制smooth的程度，實驗中帶入100

\[\begin{equation}
\min_{g, E}{\sum_{i=1}^{N}{\sum_{j=1}^{P}{{\left[g({Z}_{ij})-\ln{E_{i}}-\ln{\Delta{t_{j}}}\right]}^{2}}+\lambda\sum_{z=Z_{min}}^{Z_{max}}{g^{”}{(z)^2}}}}
\end{equation}\]

由於intensity是離散取樣的，所以該論文也定義\(f\)是一個離散的對應table，有些論文則是會將其參數化來求解(求已知basis function的線性組合)。其中\(g^{”}{(z)}\)的部分以離散Laplacian實作

\[g^{”}{(z)}=g(z-1)-2g(z)+g(z+1)\]

• Remark: 注意這個\(g^{”}{(z)^2}\)的L-2 regularization會傾向把\(g(z)\)變成接近線性，也就是使\(g^{‘}{(z)}\)接近一個常數(Gaussian prior的均質特性，會想把數字變成接近0且每個數字都很接近)

再加上intensity weighting \(w(z)\)，基本想法是sensor取樣在極值的intensity是不穩定的(容易被clamping)，所以令intensity越接近sensor中間灰的權重提高。\(w(z)\)為一個三角波，中央位置為\(Z_{mid}=(Z_{max}+Z_{min})/2\)，通常在8-bit影像是128

\[\begin{equation} \label{eq:2}
\min_{g,E}{\sum_{i=1}^{N}{\sum_{j=1}^{P}{{\left\{w(Z_{ij})\left[g({Z}_{ij})-\ln{E_{i}}-\ln{\Delta{t_{j}}}\right]\right\}}^{2}}+\lambda\sum_{z=Z_{min}}^{Z_{max}}{\left[w(z)g^{”}{(z)}\right]^{2}}}}
\end{equation}\]

\[
w(z)=\Big\{\begin{matrix}z-Z_{min},z{\le}Z_{mid}\\Z_{max}-z,z>Z_{mid}\end{matrix}
\]

• Remark: 在\(g^{”}{(z)}\)前面乘上\(w(z)\)我覺得很怪，因為這樣接近中間灰的地方smooth很強，而讓接近極值的地方smooth效果更弱，正好和論文前段的假設相左；不是該相信中間灰的數值是正確的嗎?為什麼反而要加強smooth中間灰?

4-2. Matrix form

這裡可謂本論文次困難的部分，將\(\eqref{eq:2}\)改寫成\(Ax=b\)的形式，並且求得\(x=(A^{T}A)^{-1}A^{T}b\)，將所有未知數\(g, E\)巧妙的安排到\(x\)中。將\(\eqref{eq:2}\)改寫

\[
\min_{g,E}{\sum_{i=1}^{N}{\sum_{j=1}^{P}{\left[w(Z_{ij})g({Z}_{ij})-w(Z_{ij})\ln{E_{i}}-w(Z_{ij})\ln{\Delta{t_{j}}}\right]^{2}}+\lambda\sum_{z=Z_{min}}^{Z_{max}}{\left[w(z)g^{”}{(z)}\right]^{2}}}}
\]

我們把\(Ax=b\)拆成數個部分

\[
Ax=b\Rightarrow\begin{bmatrix}a_{11}&a_{12}\\a_{21}&0\end{bmatrix}\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}=\begin{bmatrix}b_{1}\\b_{2}\end{bmatrix}
\]

其中(詳細請見Appendix)

• \(a_{11}\)描述\(w(Z_{ij})g({Z}_{ij})\)中的\(w(Z_{ij})\)，而\(a_{11}\)的column index就代表\(Z_{ij}\)
• \(a_{12}\)描述\(-w(Z_{ij})\ln{E_{i}}\)中的\(-w(Z_{ij})\)，而\(a_{12}\)的column index就代表第\(i\)個pixel 取樣點
• \(a_{21}\)描述\(w(z)g^{”}{(z)}\)

\(A\)是一個sparse matrix，在\(a_{11}\)與\(a_{12}\)的部分僅有取樣\(Z_{ij}\)的column有值，而\(a_{21}\)也僅在\(z\)的column前後有值

4-3. Unit exposure

回到\(\eqref{eq:model}\)，等式左右各乘上一個常數\(\alpha\)皆成立。

\[
\begin{split}
\alpha{Z}&=\alpha{f(E\Delta{t})}\\
\alpha{f^{-1}{(Z)}}&=\alpha{E\Delta{t}}\\
\ln{\alpha}+\ln{f^{-1}{(Z)}}&=\ln{\alpha}+\ln{E}+\ln{\Delta{t}}\\
\end{split}
\]

\[
\begin{equation} \label{eq:alphaShift}
\ln{\alpha}+g(Z)=\ln{\alpha}+\ln{E}+\ln{\Delta{t}}
\end{equation}
\]

在式\(\eqref{eq:alphaShift}\)中，\(g(Z)\)的值可能會有一個\(\ln{\alpha}\)的平移，為了方便顯示，我們定義\(g(Z_{mid})=0\)，也就是在8-bit影像中，將intensity=128對齊總照射能量\(E\Delta{t}=1\)的位置，對齊unit exposure

\[
\begin{split}
g(Z_{mid})&=0\\
\ln{f^{-1}(Z_{mid})}&=\ln{E}+\ln{\Delta{t}}&=0\\
f^{-1}(Z_{mid})&=E\Delta{t}&=1
\end{split}
\]

• Remark: 這部分是全論文最詭異的地方，附圖和說明也看不太懂，我是這樣理解的；如果你有圖組5個曝光，選3個pixel位置”各別”去做optimization，得到3組\(g(z)\)，因為這三組\(g(z)\)可能不是在同一條curve上，所以需要做對齊不同線性系統的scale；但對單一\(g(z)\)自身來看，所有取樣都是在一條curve上的。問題在什麼時候我們會針對單一pixel位置的不同曝光去做calibration? 一般都是廣泛的取樣吧? (抱頭…)
• 原文: plot of \(g(Z_{ij})\) from three pixels observed in five images, assuming unit radiance at each pixel

5. Usage

當我們求解出\(g(Z)\)和之後要如何使用? 求解任意pixel的irradiance\(E\)

\[\ln{E_{ij}}=g(Z_{ij})-\ln{\Delta{t_{j}}}\]

取自然對數\(e\)

\[E=exp(\ln{E})\]

對\(P\)張曝光的影像，簡單取平均

\[\ln{E_{i}}=\frac{1}{P}\sum_{j=1}^{P}{\ln{E_{ij}}}\]

如果考慮上wieghting \(w(Z)\)，據說可以消除極值的誤差，減低noise

\[\ln{E_{i}}=\frac{1}{W_{i}}\sum_{j=1}^{P}{w(Z_{ij})\ln{E_{ij}}}\\W_{i}=\sum_{j=1}^{P}{w(Z_{ij})}\]

 

---

 Appendix

\[
\begin{matrix}b_{1}=\begin{bmatrix}w(Z_{1,1})\ln{\Delta{t_{1}}}\\w(Z_{2,1})\ln{\Delta{t_{1}}}\\w(Z_{3,1})\ln{\Delta{t_{1}}}\\\vdots\\w(Z_{1,2})\ln{\Delta{t_{2}}}\\w(Z_{2,2})\ln{\Delta{t_{2}}}\\w(Z_{3,2})\ln{\Delta{t_{2}}}\\\vdots\\w(Z_{N,P})\ln{\Delta{t_{P}}}\end{bmatrix}_{NP\times{1}}b_{2}=\begin{bmatrix}0\\0\\\vdots\\0\end{bmatrix}_{\left|Z\right|\times{1}}\end{matrix}
\]

\[
\begin{matrix}x_{1}=\begin{bmatrix}g_{z=Z_{min}}\\g_{z=Z_{min}+1}\\\vdots\\g_{z=Z_{max}}\end{bmatrix}x_{2}=\begin{bmatrix}{\ln{E_{1}}}\\\ln{E_{2}}\\\vdots\\\ln{E_{N}}\end{bmatrix}\end{matrix}
\]

\[
a_{21}=\lambda\begin{bmatrix}w(1)&-2w(1)&w(1)&0&0&0&0&\cdots\\0&w(2)&-2w(2)&w(2)&0&0&0&\cdots\\0&0&w(3)&-2w(3)&w(3)&0&0&\cdots\\0&0&0&w(4)&-2w(4)&w(4)&0&\cdots\\0&0&0&0&w(5)&-2w(5)&w(5)&\cdots\\\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{bmatrix}
\]