Overview

As both the human eye and a camera sensor capture light, it is hard to represent the full range of possible brightnesses, from the sun itself down to the depths of an indoor shadow. However, this range can be represented across several images with different exposures (long exposure capturing the detail in dark objects, short exposure for bright objects), and the so called High Dynamic Range (HDR) for an image can be derived by finding the "true" radiance values for each pixel by analyzing multiple exposures. Once this HDR radiance map is obtained, the values for each channel can be re-formed into a new, single image, which aesthetically represents the full dynamic range of the original scene.

Algorithm

Radiance Map Construction

This is discussed in detail in Debevec and Malik 1997, but I'll summarize the process here.

Each pixel value across all exposures can be thought of as a function of scene radiance at that point and the duration of the exposure. Because this function could a be a fairly complicated response curve, it's easier to solve for the function g: the log of f's inverse. Therefore, g maps pixel values to the log of exposure values: g(Zij)= ln(Ei) + ln(tj), for pixel i in image j, which is equation 2 in the Debevec paper. Though this is infinitely under-constrained if g is considered to be a continuous function, we must observe that we only want the output for values 0-255, which is essentially 256 different equations in a linear system. Adding a smoothness constraint to this (second derivative should be zero) and an anchor constraint (directly defining one value for g), we have an over-constrained linear system of equations which MATLAB can easily solve for us. (The formulation of this matrix equation is illustrated here.) Once we have solved for the 256 pertinent values of g, we can easily reformulate the equation above (equation 2) to determine a mapping from pixel to radiance value. Because we want to incorporate the values observed for each exposure, we take the sum of g - ln(exposure) across each different exposure. This sum is weighted to give higher relevance to exposures where a pixel's value is closer to the middle of the response function; this is equation 6 in Debevec. Solving this gives us the log of the radiance for each pixel in the scene, so after exponentiating these values the radiance map is complete.

Tone Mapping

Given input from the above (i.e. linear values for radiance in each channel), we
  1. Compute the intensity by finding a weighted average from each channel (.299 * R + .587 * G + .114 * B)
  2. Run a bilateral filter on the log of this intensity: for each pixel of the input, solve this equation:
    where J(x) is the pixel value, we're summing over all values of nearby pixels to x, f is the spatial gaussian of x and a nearby pixel E, g is a gaussian of the intensity difference between pixels x and E, I(x) is the intensity value of x, and k(x) is a normalization factor (simply the sum of the products of the above gaussians). This is calculated once for each pixel, and in my implementation, I save time by using element-wise matrix multiplies at each pixel instead of using another pair of nested for-loops to accumulate the sums.
  3. Compute the detail layer by subtracting the bilateral filtered L from L itself.
  4. Apply an offset (the maximum value in the bilateral filtered L) and scale the result to the base image.
  5. Reconstruct the log intensity of the recombined offset/scaled result and detail layer (by exponentiation).
  6. Reincorporate the color channels.
  7. Apply gamma compression to brighten the result.
For more detail on the above process, consult the assignment handout and bilateral filter slides, or read through the Durand 2002 paper.

Results

Following is the breakdown for one HDR image generation result, followed by more abbreviated results for other sets of images.
Original Images (Mouse Over)
Radiance Map
Exposure vs. Pixel Values
Simple Global Tone-Mapping
Bilateral Detail
Bilateral Large Scale
Durand Tone Mapping

Original Images (Mouse Over)
Simple Global Tone-Mapping
Durand Tone Mapping
November 12th, 2012