CS129 / Project 5 / HDR

1.Radiance map construction

The goal of this part is to build an HDR radiance map from several LDR exposures. To achieve this, we are supposed to recover the response function g which shows the relationship between the pixel values and exposures in log domain (g(Z) = ln(exposure)), and then reconstruct the radiance map E according to ln(Ei) = g(Zij)-ln(tj).

To obtain g(Zij), where Zij ranges from 0 to 255, we could transfer the issue to linear least squares problem(as the equation below) and solve it using singular value decomposition (SVD) method(x = A\b, x is combined with the values of g(0)--g(255) and lnEi for the sampled pixels).

1.1 Results visualization

(1) Plots of obtained response function g for R,G,B channels:

Red
Green
Blue
Combined

(2) Radiance map for R,G,B channels:

Red
Green
Blue

2.Tone mapping

2.1 Global tone mapping

If we scale the radiance map we obtained from last session to [0,1] range directly, the result will not be reasonable. Because it compressed the radiance value uniformly. To improve this, global tone mapping, for instance log method, compress the extreme values of radiance map more heavily than the middle values, which makes the result better. In this session, I used L/(L+1) method to do global tone mapping and scale the results to [0,1] range.

2.2 Local tone mapping

The main idea of this local tone mapping method is that first we obtain the intensity image of the input HDR image, and then get the color information by subtracting intensity from input image. Next we use bilateral filter to blur the intensity image to get the low frequency and also edge information, which we call large-scale image. And then we can get the detail information by subtracting large-scale from intensity image. After reducing contrast of large-scale image, we add the detail information and color information back to get a reasonable tone mapping result.

Here is the basic implementation of the bilateral filter:


%Pseudocode
1. Set sigma_space and sigma_range.

2. Make spacial gaussian kernel in advance, where kernel radius = 3*sigma_spacial.

3. For each pixel p in the input_image
      p_value = 0;
      weight = 0;

      for each pixel q inside the range of spacial kernel
         p_value += Gaussian_spacial(|p-q|)*Gaussain_range(|Ip - Iq|)*Iq;
         weight += Gaussian_spacial(|p-q|)*Gaussain_range(|Ip - Iq|);
      end
      
      //normalize
      p_value /= weight;

   end

One thing I have noticed is that if I set the sigma_spacial and sigma_range to a relative small value(smaller kernel size), the final result will not be sharp enough compared with what it should have been. The reason I think is that when the sigma is small, the output of bilateral filter(large scale) will not be very blurred which results in less detail information because detail = intensity image - large scale. After reducing contrast on the large scale image and adding detail information back it loses some details. On the other hand, if we choose relative larger sigma, enough detail information will be retained, but it is time consuming based on basic implementation.

sigma_spacial=9 sigma_range=0.1
sigma_spacial=32 sigma_range=1.76

2.3 Result images

Scale
Global
Local
dR = 6, Gamma = 0.7
dR = 7, Gamma = 0.8
dR = 7, Gamma = 0.7
dR = 12, Gamma = 0.9
dR = 14, Gamma = 0.8
dR = 10, Gamma = 1
dR = 8, Gamma = 0.7
dR = 9, Gamma = 0.5

3.Fast bilateral filter

When doing bilateral filter, if we use relatively large gaussian kernel, it will become pretty time consuming. As we can see from the basic implementation pseudocode, it has to traverse every pixel of input image and do gaussian filter on the spacial and range domain on the fly, so the time complexity is O(n^4) in this part. Because the shape of bilateral filter is different every time, so it is impossible to precompute or use FFT to speed up. 'A Fast Approximation of the Bilateral Filter using a Signal Processing Approach' raises a effective method, in which the input image is downsampled to a much smaller size and filtered by 3D gaussian filter(2D-spacial and 1D-intensity range) with also very small size.

Here is the method:


1. Set sigma_space and sigma_range and make samplingRateSpatial = sigmaSpatial and samplingRateRange = sigmaRange.

2. Set sampledSigmaSpatial = sigmaSpatial / samplingRateSpatial and sampledSigmaRange = sigmaRange / samplingRateRange.

3. Make 3D gaussian kernel in advance with sampledSigmaSpacial and sampledSigmaRange.

4. For each pixel p in the input_image
      //downsample the spacial*range space on the fly by box filter
      (round(p.x/samplingRateSpacial), round(p.y/samplingRateSpacial), round((Ip-Imin)/samplingRateRange))
   end

5. Convoluting downsampled 3D space and according w by subsampled 3D gaussian kernel made by step 3.

6. Slicing and normalizing.

7. Interpolate the result of step 6.

As we can see above, downsampling session is the only place we have to traverse the input image and convolution part will not cost much time because the downsampled image and the kernel size are both very small. It only take 2-3 seconds when filtering the 'Chapel' HDR intensity image(sigma_spacial=32, sigma_range = 1.76), which accelerated a lot compared with basic implementation with the same sigma and that takes around 30 minutes on the same image.