We are presenting our results from an implementation project done at Marc Alexa's Computational Photography seminar at TU Berlin in January 2007. The aim of this project was to implement a recent computational photography technique, we have chosen HDR imaging and tonemapping.
Creating HDR images usually requires a series of exposures of the same scene, then the following three steps are applied to compute a tonemapped version of the resulting hdr image:
We have implemented both papers in Matlab, see Downloads for the code.
Both papers are very well written and explain their subject in great detail. During implementation however, we noticed some important subtleties not explicity mentioned that we would like to share with the reader:
The hdr radiance map is computed as a weighted sum of the recovered radiance values from the pixels of each exposure. Saturated pixels, i.e. pixels with Zij = 255 need to be ignored in that computation since the camera thresholds all radiance values bigger than a certain value to Z = 255. We therefore cannot recover the exact radiance value of Zij in case Zij = 255. Note that the corresponding pixels at the other two color channels also have to be exluded from the weighted sum.
Given a set of pictures with an exposure range of e.g. 1s, 1/4s, 1/15s, 1/60s and 1/250s we expect that a pixel that is saturated in the image exposed with 1/250s is saturated in all the other exposures as well. We found that due to sensor noise it is quite possible that this is not the case and we may find a non-saturated pixel with e.g. Z=254 in one of the other exposures. This unfortunately has a strong influence on the result of the weighted sum and results in artifacts in the resulting hdr radiance map. To solve this problem we are processing the images in descending order. When encountering a saturated pixel, we assume that it should have been saturated in all previous pictures (longer exposures) as well. Therefore, we cancel the influence of all the previous images on the weighting function. If a pixel is saturated in all exposures we approximate the resulting radiance value only from the image with the shortest exposure time.
The Reinhard local tonemapping operator finds the biggest possible neighbourhood for each pixel that still has a fairly even lumiance variation. Given a range of eight discrete scales sm with s1 = 1.6, s2 = 1.6^2, ..., s8 = 1.6^8, the corresponding V(x,y,sm) = 0.01, 0.015, 0.017, 0.02, 0.04, 0.05, 0.07, 0.03 and a thresholding value eps = 0.05. According to Reinhard's paper we need to "seek the first scale sm where: |V(x,y,sm)| < eps". This however leads to selecting the first scale s1 in most cases, especially in smooth areas of the image and basically degrades the operator to its global counterpart.
Instead, in our implementation, we seek the biggest scale s_biggest with |V(x,y,s_biggest)| < eps starting from the lowest scale while all smaller scales s_smaller with s_smaller < s_bigger also fulfill |V(x,y,s_smaller)| < eps. In our example this would result in selecting s5.
We have computed the camera response curve of a Canon EOS 350D from a set of nine pictures with exposure settings ranging from 1/4000 second to 15 seconds each four times bigger than the last one.
Since it is impossible to display a hdr radiance map on a computer screen we are giving three simple compressed representations:
![Linear mapping of a hdr picture into range [0,1]](windowLinear.jpg)
![Logarithmic mapping of a hdr picture into range [0,1]](windowLog.jpg)
The first (left) image is mapped with the Reinhard global operator, the second (right) image with the local operator.
You may contact me at m.eitz at tu-berlin dot de