Edit.
I have produced an image which is "cleaner" looking than my original attempt, and the processing is faster too.
As before we start by loading the images in order from darkest to brightest, and cropping away the artifacts from alignment.
files = Reverse@FileNames["memorial*.png"];
images = ImagePad[Import[#], {{-2, -12}, {-35, -30}}] & /@ files;
HDR image construction:
Like Thies's approach, this uses averaging over multiple images to obtain pixel values. Intensity data is extracted from the images and low or high values are set to zero to flag potentially noisy or saturated pixels. The exposure ratio between two images is estimated by considering only those pixels which are non-zero in both images. After compensating for exposure differences, the pixel values from all images are combined into a mean image, again using only the non-zero pixels. Finally the mean image is split into HSB components.
data = ImageData[First@ColorSeparate[#, "Intensity"]] & /@ images;
data = Map[Clip[#, {0.1, 0.97}, {0, 0}] &, data, {3}];
exposureratios = Module[{x, A, g},
First@Fit[Cases[Flatten[#, {{2, 3}, {1}}], {Except[0], Except[0]}, 1],
x, x] & /@ Partition[data, 2, 1]];
exposurecompensation = 1/FoldList[Times, 1, exposureratios];
data = MapThread[Times, {exposurecompensation, Unitize[data] (ImageData /@ images)}];
data = Transpose[data, {3, 1, 2, 4}];
meanimage = Map[Mean[Cases[#, Except[{0., 0., 0.}]]] &, data, {2}];
{h, s, b} = ColorSeparate[ColorConvert[ImageAdjust@Image[meanimage], "RGB"], "HSB"];
Tone mapping:
We now have a brightness channel containing a range of values from 0 to 1, but with a very non-uniform distribution.
ImageHistogram[b]
First I do a histogram equalisation on the brightness data:
cdf = Rescale@Accumulate@BinCounts[Flatten@ImageData@b, {0, 1, 0.00025}];
cdffunc = ListInterpolation[cdf, {{0, 1}}];
histeq = Map[cdffunc, ImageData[b], {2}];
ImageHistogram[Image@histeq]
Next I apply a sort of double-sided gamma adjustment to reduce the number of very low and very high values (we don't want too many deep shadow or bright highlights).
b2 = Image[1 - (1 - (histeq^0.25))^0.5];
ImageHistogram[b2]
Final image:
Finally I apply a built-in Sharpen filter to the new brightness channel, to boost local contrast a little bit, and apply a gamma adjustment to the saturation channel to make it a little more colourful. The HSB channels are then recombined into the final colour image.
ColorCombine[{h, ImageAdjust[s, {0, 0, 0.75}], Sharpen[b2]}, "HSB"]
Original version
Here's an attempt at the Stanford Memorial Church image using a local contrast filter to do the tone mapping.
First load the images and crop to remove the artifacts around the edges of some of them:
files = Reverse @ FileNames["memorial*.png"];
images = ImagePad[Import[#], -40] & /@ files;
Next create small grayscale versions and use these to estimate the brightness scaling between the images
small = ImageData[ImageResize[ColorConvert[#, "Grayscale"], 50]] & /@ images;
imageratios = FoldList[Times, 1,
Table[a /. Last@
FindMinimum[Total[(small[[i]] - a small[[i + 1]])^2, -1], {a, 1}],
{i, Length@small - 1}]]
Now select the "best" image from which to take each pixel value, and scale that value accordingly. I've defined the "best" image for a given pixel as the one for which the median of the {R,G,B} numbers is closest to 0.5.
data = Transpose[ImageData /@ images, {3, 1, 2, 4}];
bestimage = Map[Module[{best},
best = Ordering[(Median /@ # - 0.5)^2, 1][[1]];
#[[best]]*imageratios[[best]]] &, data, {2}];
Next apply a local contrast enhancement to the brightness channel of the image. This is quite simple and slow. For each pixel the filter sorts the unique values in the pixel's neighbourhood and finds the pixel's position in that list. The pixel value is set to the fractional list position. For example if a pixel is the brightest one in its neighbourhood, it gets a value of 1. The size value in the localcontrast function must match the range parameter in the ImageFilter.
localcontrast = With[{size = 20}, Compile[{{x, _Real, 2}},
Block[{a, b, val}, val = x[[size + 1, size + 1]];
a = Union[Flatten[x]];
b = Position[a, val][[1, 1]];
b/Length[a]]]];
{h, s, b} = ColorSeparate[ColorConvert[Image[bestimage], "RGB"], "HSB"];
newb = ImageFilter[localcontrast, b, 20];
Finally combine the contrast-enhanced brightness channel with the original saturation and hue to get the final image:
ColorCombine[{h, s, newb}, "HSB"]
It's not brilliant, but I think the general HDRI effect is there. The contrast enhancement could probably be toned down a bit by increasing the size parameter, though it'll be slower.
CreateHDRIandToneMapHDRIin theImagecontext. (I don't know what they do, if anything.) – Simon Woods Aug 12 '12 at 17:23