HDR Imaging and Tonemapping

§ The Problem

The real world is in high dynamic range. A camera sensor is unable to capture in HDR; as bright spots may oversaturate a photosensor and force it to white -- losing information in that area. We can modify our exposure time to adjust the amount of light that enters the sensor -- but there is no one exposure time that will give us a good picture, so we must take many.

§ Many-to-one

Camera Pipeline:

1. 3D Radiance ($L = \frac{W}{srm^2}$)
2. Lens
3. 2D Sensor (Irradiance, $E = \frac{W}{m^2}$)
4. Shutter (Exposure time, $\Delta t$)
5. Sensor Exposure $H = E\Delta t$
6. ADC (Analogue to Digital Converter)
7. Raw Image
8. Remapping
9. Pixel Values ($Z$)

Our goal is to go from $Z$ to $E$, however, to fully go to $E$ we must reverse the ADC -- we have no interest or intention of doing so. We can note that $E$ is invariant to exposure time, so if we know $\Delta t$ we can find $E$ -- typically $\Delta t$ is given as metadata on an image.

We need multiple images to recover the $Z$ vs $H$ curve -- typically with varying $\Delta t$.

§ Camera Response Curve

Then, our goal is to find $E_i$ or the exposure sensitivity for the $i$^th pixel.

Objective: Solve for $E_i$ in $Z_{i,j} = f(E_i\Delta t_j)$, where $Z_{i,j}$ represents the value of pixel $i$ on image $j$, and $\Delta t_j$ is the exposure time on the $j$^th image.

Using algebra, we get

$$\begin{align*} f^{-1}(Z_{i,j}) &= E_i\Delta t_j\\ \log f^{-1}(Z_{i, j}) &= \log E_i + \log \Delta t_j\\ g(Z) &= \log \circ f^{-1} (z)\\ g(Z_{i, j}) &= \log E_i + \log \Delta t_j. \end{align*}$$

In the above expression, the values $Z_{i,j}$ and $\Delta t_j$ are known, and $g(Z)$ and $E_i$ are unknown.

Through optimization, we may solve for the unknown values by via the following

$$[g, E_i] = \argmin_{g,\ E_i} \sum_{i=1}^N\sum_{j=1}^P \left[g(z_{i,j}) - \log E_i - \log \Delta t_j\right]^{2},$$

where $N$ is the number of pixels, $P$ the number of images.

This makes the assumption that the images are reliable -- however, this is not the case, we modify our equation to

$$[g, E_i] = \argmin_{g,\ E_i} \sum_{i=1}^N\sum_{j=1}^P \Big[w(z_{i,j})\left(g(z_{i,j}) - \log E_i - \log \Delta t_j\right)\Big]^{2} + \lambda \sum_{z=0}^{2^b} \left[w(z)\frac{d^2}{dz^2}g(z)\right]^2,$$

1. $w(z)$ is a function that assigns a weight to a pixel value; we want large $z$ values to map to smaller values, as these are generally unreliable. One possible function is the triangle function given by

$$w(z) = \begin{cases} z, &\text{ if } z \leq 2^b - 1\\ 2^b - (z + 1), &\text{otherwise} \end{cases}$$

where $b$ is the bit-depth of the image. 2. There is no reason for a manufacturer to use a noisy curve; therefore, by making the assumption that the camera's response curve is smooth,

$$\frac{d^2}{dz^2}g(z) \approxeq g(z - 1) - 2g(z) + g(z + 1) = 0.$$

3. $\lambda$ is a factor used to balance the two terms -- a larger $\lambda$ gives a smoother response curve.

We only need to perform this optimization once per camera, then we can reuse it for other images; sometimes this method is called "irradiance calibration".

§ Tonemapping

Previously we solved for $g: Z \to E,$ $g(Z_{i,j}) = \log E_i + \log \Delta t_j,$ or for every pixel we get $P$ estimates for $\log E_i,$ we cannot just use the arithmetic mean as the highly irradiant areas in longer exposures will bias the mean.

We can then reuse the weight function $w(Z)$ to get a weighted mean,

$$\log E_i \approxeq \frac{\sum_{j=1}^P w(Z_{i, j})\left[g(Z_{i,j}) - \log \Delta t_j\right]}{\sum_{j=1}^P w(Z_{i, j})}.$$

Having recovered the irradiance of each pixel in an image, we can now use it in various applications such as reusing the irradiance in a computer rendering or in viewing an image -- however, we must compress the high dynamic range into the range displayable by a device.

• Clipping
  • $I = \min (E, I)$
  • Causes the loss of detail in highly irradiant areas.
• Scaling
  • $I = \frac{E}{\max E}$
  • Causes the loss of detail in areas of low irradiance.
• Reinhard Global Operator
  • $I = \frac{E}{1 + E}$
  • Better than scaling and clipping as it increasly sharply at small values (boosting), and compresses the high irradiance areas.
• Gamma Compression
  • $I'' = (I')^\gamma, I' = \frac{E}{\max E}$
  • We first perform scaling, then a gamma correction to boost the smaller values like Reinhard.
  • Colour gets washed out.
• Gamma Compression on Intensity
  • Intensity is given by $L = \frac{1}{3}(R + G + B)$
  • Then, gamma compression on intensity $L'' = (L')^\gamma, L' = \frac{L}{\max L}$
  • Colour is given by ParseError: KaTeX parse error: Got function '\\' with no arguments as argument to '\left' at position 21: …frac{1}{L}\left\̲\̲{R, G, B\right\…
  • Then, the final image is computed via $I'' = L'' C$
  • This preserves colour, but instead washes out intensity -- or lowers the constrast of the image.

§ Chiu et. al (1993)

In Chiu^[1], the authors propose the following method:

$$I' = I - (I \star f(x)) = I - J,$$

where $f(x)$ is the Gaussian function,

$$f(x) = \frac{1}{\sqrt{2\pi}\sigma^2}e^{-\frac{x^2}{2\sigma^2}}.$$

Their process can be visualised by the following image,

Their motivation was to use a lowpass filter (here, we used the Gaussian as an illustration) to isolate ranges of radiance that are perceptually important.

Unfortunately, this approach leads to a "halo" effect around sharp irradiance transitions due to undesirable easing of a transition in a step-like function.

To understand why, we turn towards the definition of the convolution,

$$J(p) = \sum_{q} f(p - q) I(q).$$

We may notice that $f(p - q)$ is a weighting on $I(q)$ that is dependent on the spatial distance between the two pixels $p$ and $q.$

Visually, we can see the effect that this spatial weighting has on the step-like result that we desire

§ Tomasi and Manduchi (1998)

In Tomasi^[2], they explore a method called "bilateral filtering" to reduce the halo effect.

The method is given by the following adjustment the convolution calculation

$$J(p) = \frac{1}{k(p)} \sum_q f_s(p - q) f_c (I(p) - I(q)) I(q), \tag{1a}$$

where

$$k(p) = \sum_q f_s(p - q) f_c(I(p) - I(q)). \tag{1b}$$

The $f_s$ and $f_c$ functions represent the spatial and colour filters -- now instead of weighting just the spatial distance of two pixels, we also use the "colour distance".

§ Durand and Dorsey (2002); Paris and Durand (2006)

In Durand^[3], they explore a method of computing an approximation of the bilateral filter that is $\mathcal{O}(1)$ per pixel.

A further improvement upon Durand^[3:1] was made in Paris^[4]. The authors note that the product of the spatial and range Gaussian defines a higher dimensional Gaussian in the 3D product space between the domain and range of the image.

However, the definition of the bilateral filter in eq. 1 is not a convolution as it is a summation over the 2D spatial domain. The authors introduce a new dimension $\zeta$ for the intensity for each point in the product space in order to define a summation of the 3D space.

We now define for each point $(x, y, \zeta)$ an intensity $I$, and we let $\mathcal{R}$ be the interval for which intensites are defined. Also, we define the Kronecker symbol $\delta(\zeta)$ by $\delta(0) = 1, \delta(\zeta) = 0.$

$$J(p) = \frac{1}{k(p)} \sum_{q \in \mathcal{S}}\sum_{\zeta \in \mathcal{R}} f_s(||p - q||) f_c(|I(p) - \zeta|) \delta(\zeta - I(q)) I(q),$$

where

$$k(p) = \sum_{q \in \mathcal{S}}\sum_{\zeta \in \mathcal{R}} f_s(||p - q||) f_c(|I(p) - \zeta|) \delta(\zeta - I(q)).$$

We can use homogeneous coordinates to express the same equation,

$$\begin{pmatrix} W(p)J(p)\\\\ W(p) \end{pmatrix} = \sum_{q \in \mathcal{S}}\sum_{\zeta \in \mathcal{R}} f_s(||p - q||) f_c(I(p) - \zeta) \delta(\zeta - I(q)) \begin{pmatrix} W(q)I(q)\\\\ W(q) \end{pmatrix},$$

where $W(q) = 1.$

Note that the $\delta(\zeta - I(q))$ terms are cancelled when $\zeta \neq I(q)$.

The summation is now over the space $\mathcal{S} \times \mathcal{R},$ the product $f_s f_c$ defines a separable Gaussian kernel $f_{s, c}$ on $\mathcal{S} \times \mathcal{R},$

$$f_{s, c}: (x \in \mathcal{S}, \zeta \in \mathcal{R}) \to f_s(||x||) f_c(|\zeta|).$$

Then, we construct the functions $i$ and $w$,

$$i: (x \in \mathcal{S}, \zeta \in \mathcal{R}) \to I(x),$$

$$w: (x \in \mathcal{S}, \zeta \in \mathcal{R}) \to \delta(\zeta - I(x)).$$

Then, we rewrite the 3D summation

$$\delta(\zeta - I(q)) \begin{pmatrix} W(q)I(q)\\\\ W(q) \end{pmatrix}= \begin{pmatrix} \delta(\zeta - I(q))W(q)I(q)\\\\ \delta(\zeta - I(q))W(q) \end{pmatrix}= \begin{pmatrix} w(q, \zeta)i(q, \zeta)\\\\ w(q, \zeta) \end{pmatrix}.$$

Using separability of the 3D Gaussian,

$$\begin{pmatrix} W(p)J(p)\\\\ W(p) \end{pmatrix}= \sum_{(q, \zeta) \in \mathcal{S}\times\mathcal{R}} f_{s, c}(p - q, I(p) - \zeta) \begin{pmatrix} w(q, \zeta)i(q, \zeta)\\\\ w(q, \zeta) \end{pmatrix}.$$

This formula corresponds to the value of the convolution between $f_{s, c}$ and the two-dimensional function $(wi, w)$ at $(p, I(p))$,

$$\begin{pmatrix} W(p)J(p)\\\\ W(p) \end{pmatrix}= \left[ g_{s, c} \otimes \begin{pmatrix} wi\\\\ w \end{pmatrix} \right] (p, I(p)).$$

The functions $i^b$ and $w^b$ are given by

$$(w^b i^b, w^b) = g_{s, c} \otimes (wi, w).$$

Finally, the bilateral filter can be given by

$$J(p) = \frac{w^b(p, I(p))i^b(p, I(p))}{w^b(p, I(p))}.$$

First, we perform slicing or evaluating $w^bi^b$ and $w^b$ at $(p, I(p)),$ then division in order to normalize the result.

---

1. K. Chiu, M. Herf, P Shirley, S. Swamy, C. Wang, K. Zimmerman. "Spatially Nonuniform Scaling Functions for High Contrast Images" (1993). ↩︎

2. C. Tomasi, R. Manduchi. "Bilateral Filtering for Gray and Color Images" (1998). ↩︎

3. F. Durand, J. Dorsey. "Fast Bilateral Filtering for the Display of High-Dynamic-Range Images" (2002). ↩︎ ↩︎

4. S. Paris, F. Durand. "A Fast Approximation of the Bilateral Filter using a Signal Processing Approach" (2006). ↩︎