Implementation

1. Image alignment method: median threshold bitmap algorithm

2. HDR creation method: Paul Debevec's method

3. Tone mapping method: gradient domain high dynamic range ompression

I. Images Alignment Method: Median Threshold Bitmap

Here are the procedures of MTB algorithm that we had implemented.

1. Convert the input image source1 and source2 into grayscale image img1 and img2 using the following formula:

grey = (54*red + 183*green + 19*blue) / 256

2. Convert img1 and img2 into binary bitmaps b1 and b2 using the median threshold of each image. If the greyscale value of a pixel is bigger than or equal to the threshold, set the bit value in the bitmap to 1, otherwise to 0.

3. Repeat procedure 2 to the image pyramid (we down sampling the images by the ratio of 2, and build a image pyramid with depth equal to 6).

4. Shift b2 to nine directions ((1,1), (1,0), (1,-1), (0,1), (0,0), (0,-1), (-1,1), (-1,0), (-1,-1)) in each level, and computing the differences with b1. Find the direction that brings the smallest error. We will not count the differences of a pixel if the greyscale value in img1 or img2 is close to our threshold (if the difference is less than 4).

5. We start from the smallest level, and counting the shift direction that brings the smallest error. Using this direction, and enlarging the direction twice in the next larger level, and doing procedure 4 again. Keeping this recusion until reachs the biggest level, and we can get the final shift direction in x and y direction.

Here is the test of two different exposured images:

 input image 1 intput image 2 result shifted image (click the images to see original size image)

Our MTB program shifted the image2 with (28,-17) to fit image1. We filled the area that don't be filled to black.

II. High Dynamic Range Image Generation

We implemented the HDR algorithm of Paul E. Debevec's paper. We uses multiple jpeg image as input, and recovered the response function from them. Matlab is used in this part, and we applied the Matlab function that the author provided in the paper.

The following are two example:

1. Stair

 8s 5s 2.5s 1s 1/2s 1/4s 1/8s 1/15s 1/30s 1/60s 1/320s 1/400s 1/640s

response curve
R
G
B

We used 13 jpeg images and 50 sample pixels in this example. The result curves are quite smooth except for the top of the R channel.

2. Bridge

 30s 25s 20s 15s 13s 10s 8s 6s 5s 4s 3s 2.5s 2s

response curve
R
G
B

We used 13 jpeg images and 500 sample pixels in this example. The result response function curves still not monotonic increasing, I think this is due to the scene. The flowing river, and the passing cars on the bridge, and there was also landing airplane appearing in my photo sequence. All of those moving objects bring noise to the optimization, and leading to the zigzag response curves.

III. Tone Mapping

After producing the HDR image, because of the range of luminance is too high, we need to do tone mapping, to lower the range.

The tone mapping algorithm we use operates on gradient domain. The algorithm proceeds as follows:

1.      Take the logarithm of the image luminance.

For color image, we set the luminance Y = 0.299*R+0.587*G+0.114*B. Then we store log(Y) into a matrix, H.

2.      Compute the attenuation map, PHI.

Construct a Gaussian pyramid H0, H1,…, Hd, where H0 is the full resolution HDR image and Hd is the coarsest level in the pyramid. At each level k we compute a scaling factor:

, where α=0.1 and β is between 0.8 and 0.9.

The process of computing PHI is given by following equations:

L is an upsampling operator with linear interpolation.

3.      Compute the attenuated gradient field, G.

4.      Find the image I which gradient is closet to G.

Since G is not necessarily integrable, there might not exist an image I such that G=I. So we can only find an approximate solution I which minimize

5.      Compute the output image luminance, Lout.

At the beginning we convert Y into logarithm domain, so we need to convert our result back. That is, Lout=exp(I).

6.      For color image, recover the color of the image.

, where s controls the color saturation, and we set s between 0.4 and 0.6.

We use library gil for image IO, and solving image I by library LSQR.

Here is some of our result:

 Input image Result 1 Result 2 Belgium.hdr β = 0.8, s = 0.4 β = 0.9, s = 0.6 stairs.hdr β = 0.85, s = 0.5 β = 0.9, s = 0.6 bridge.hdr β = 0.8, s = 0.6 β = 0.9, s = 0.6

Artifacts

References

2. Paul E. Debevec, Jitendra Malik, Recovering High Dynamic Range Radiance Maps from Photographs, SIGGRAPH 1997.

3. Raanan Fattal, Dani Lischinski, Michael Werman, Gradient Domain High Dynamic Range Compression, SIGGRAPH 2002.