# Update

Submitted answers here are tour de force, and very educational comp-sci-wise, -- thanks folks! Tech is developing though, and here is the current implementation in the Wolfram Language. Please study the relevant docs and options of these functions (ImageExposureCombine and ColorToneMapping) and also related functions, for example ImageFocusCombine. So, we got 16 images under different exposures:

imgs=Import["http://www.pauldebevec.com/Research/HDR/SourceImages/Memorial_SourceImages.zip","*.png"];

Thumbnail /@ imgs


This one-liner does it (with some arbitrary choice of parameters, feel free to experiment):

ColorToneMapping[ImageExposureCombine[imgs,"HDR"],.7]


# Original post

The High dynamic range imaging (HDR or HDRI) direction in photography and image processing became very popular recently. Besides obvious photo art applications (see examples), there are many great opportunities in computer vision and graphics.The German Max Planck Institute has a page dedicated to the pipeline of related technologies, publications and discussions.

I wonder if our community has members who are knowledgeable enough to reproduce a correct algorithm in Mathematica's top level image processing functions. I hope the answers may also serve as a good intro to methodology from the programmatic point of view. I am sure code will be quite short and consider this to be a very practical down to earth question. There are many Photoshop-like tutorials online. I won't link to them because they are quite easy to find and I cannot judge which one reflects the right approach. Read a concept article about HDR here.

I found this interesting research site with Mathematica notebook and manual provided (read down the page). Unfortunately it gives only part of the process – tone mapping, and is written in Mathematica version 6, while there have been many image processing upgrades in the current version 8.

Precise question formulation:

• Consider a set of images taken at different exposures – like these public images
• Write Mathematica code that produces an HDR image – like the one here labeled “Local tone mapping”
• Implement any HDR algorithm (there are a few) by any Mathematica means possible

• I wanted to ask this question myself, but I couldn't figure out how to ask it and lazily never did. +1 Commented Aug 12, 2012 at 8:28
• Hey! This is wonderful. You can get very nice effects just performing basic image operations with those four images in Wikipedia :) i.sstatic.net/4fhd1.png Commented Aug 12, 2012 at 8:46
• @Mr.Wizard Thanks for endorsing the question - I was hoping to set creative powers in motion ;-) Commented Aug 12, 2012 at 9:11
• @Verde Thanks! You're absolutely right - there are quite a few other techniques that avoid HDR and extract useful information form different exposures in some other ways. Simplest way is to use masking to "highlight" the best areas in different images and image-add them. Cool image btw ;-) Commented Aug 12, 2012 at 9:15
• Version 8 has the intriguing symbols CreateHDRI and ToneMapHDRI in the Image context. (I don't know what they do, if anything.) Commented Aug 12, 2012 at 17:23

## 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.

• "Not brilliant"? I think the HDR effect here is actually better than all three of those shown on Lischinski's website! Commented Aug 13, 2012 at 16:52
• +1 Very interesting, nice image, thank you! I could not figure out how to execute your code - some variables are not defined - like sgd for example. Also your "local contrast enhancement" - is it a known HDR method or where you just nicely improvising with mma? Commented Aug 13, 2012 at 23:31
• @VitaliyKaurov, sorry - I changed variable names before posting and missed a couple. It should now work. Regarding the local contrast enhancement, I think it's something that is often used with HDR images but the implementation is improvised. Commented Aug 14, 2012 at 8:26
• @Oleksandr while I think this is a fine beginning I cannot agree with your assessment. The final image here is very "lumpy" with heavy variations in luminance across areas that should be relatively uniform; avoiding this unnatural result, at least perceptually, is a primary goal of (photographic) dynamic range compression. (Please don't be offended Simon; this is a complex subject.) Commented Aug 16, 2012 at 13:50
• @Mr.Wizard yes, i got the tools ;) Commented Aug 17, 2012 at 16:43

For starters i tried an easy intuitive approach, namely, combining the best parts from each image adjusted for the different exposure times all into one HDR image.

Let's start by importing all the images

imageurls = "http://upload.wikimedia.org/wikipedia/commons/thumb/" <> # & /@
{"0/09/StLouisArchMultExpEV-4.72.JPG/320px-StLouisArchMultExpEV-4.72.JPG",
"c/c3/StLouisArchMultExpEV-1.82.JPG/320px-StLouisArchMultExpEV-1.82.JPG",
"8/89/StLouisArchMultExpEV%2B1.51.JPG/320px-StLouisArchMultExpEV-1.51.JPG",
"8/8f/StLouisArchMultExpEV%2B4.09.JPG/320px-StLouisArchMultExpEV-4.09.JPG"};

FPImage = Image[#, "Real"]& (* change image rep. to floating point *)
inputimages = SortBy[
Composition[FPImage, Import] /@ imageurls,
(* sort images by increasing exposure *)
Composition[Mean, Flatten, ImageData, First, ColorSeparate[#, "Intensity"] &]
]


where we changed the internal representation to floating point to avoid rounding errors later and already sorted in ascending order by exposure times.

Next we compare neighbour pairs of images and try to guess their exposure ratios. (Probably there is a more robust and effective method, but for first try it seems to work and give reasonably exact values).

ScaledImageSquaredError[img1_Image, img2_Image, s_?NumericQ] :=
Composition[Total, Flatten, ImageData, ImageApply[#^2 &, #] &,
First, ColorSeparate[#, "Intensity"] &][
ImageDifference[img2, ImageClip@ImageApply[s # &, img1]]
];

(* Guess exposure ratio. img2 must be the brighter than img1 *)
GuessExposureRatio[img1_Image, img2_Image] := NArgMin[
{ScaledImageSquaredError[img1, img2, \[Alpha]],
1 < \[Alpha] < 10},
\[Alpha], Method -> {"NelderMead", "RandomSeed" -> 24}
]

exposureratios = GuessExposureRatio @@@ Partition[inputimages, 2, 1]
(* {4.66667, 5.89289, 3.74074} *)

(* brightness scaling factors to get all images to the same exposure level *)
unifyingexposures = #/Max[#] &@Reverse@FoldList[Times, 1, Reverse@exposureratios]
(* {1., 0.214286, 0.0363634, 0.00972091} *)


Next we want to generate masks for each image, marking the usable portion of each image. Our criterium will be the brightness of each pixel. We don't want the low-brightness range because of possibly bad signal-to-noise ratios and we don't want the brightnesses near 1.0 because of possible clipped values.

Therefore we define a (kind of arbitrary and adjustable) brightness weight curve

brightnessweights = Function[b,
Piecewise[
{{Interpolation[{{{0}, 0, 0}, {{0.2}, 1, 0}, {{0.95}, 1, 0}, {{1}, 0, 0}},
InterpolationOrder -> 3][b],
0 <= b <= 1}}, 0
]
]
Plot[brightnessweights[b], {b, 0, 1}, PlotRange -> All, ImageSize -> Small]


and see how the different exposure regions overlap and combine

LogLinearPlot[
Evaluate[brightnessweights[b/#] & /@ unifyingexposures],
{b, 0.0002, 1},
PlotRange -> All, Filling -> Axis
]


We have good overlap, but not too much, and no gaps in brightness levels, so that looks nice.

Next, we'll construct the weight maps by applying the brightness curve to the brightnesses of each image pixel and then normalizing by the sum of all weight maps at each pixel.

GenerateWeightMask = Composition[
ImageApply[brightnessweights, #] &, First, ColorSeparate[#, "Intensity"] &
]
ImageTotal = Fold[ImageAdd, First@#, Rest@#] &

nrimages = Length[inputimages];


Now we multiply our input images with the weight masks and adjust for the different exposure levels and then combine all those into the final HDR composite image.

sameexposureimages = Function[{s, i}, ImageApply[s # &, i]] @@@
Transpose[{unifyingexposures, inputimages}]
hdrimage = ImageTotal[
]


which we can now watch at different artificial exposures

Manipulate[ImageMultiply[hdrimage, \[Alpha]], {{\[Alpha], 1}, 1, 100}]


Also we could now start compressing the dynamic range again by applying tone mapping algorithms and the like.

Update An example for an easy global tone mapping algorithm can be achieved by applying #/(#+C) & to every pixel intensity which looks like this:

Manipulate[
ImageApply[#/(# + \[Alpha]) &, hdrimage], {{\[Alpha], 0.2}, 0.001, 0.2}
]


• +1 This is great, thank you! Your hdrimage, though cleverly constructed, looses some information from the original image set - am I correct? Even so I like the idea. I wonder if it is possible to design similar thing with true HDR loss-less representation and proceed with some cool local tone mapping. Possible or not - thank you anyway! Commented Aug 13, 2012 at 23:20
• Thanks for your feedback! Most of the information of the original images is still there in the hdrimage spanning about 10 f-stops, by combining the best parts of each image. So you could export it as .exr and would actually get a useful HDR representation out of it. Commented Aug 14, 2012 at 3:47
• The information that is lost are basically the very low and very high brightness parts in each image making it easier to get rid of the artifacts that come with them. To get even more information out of the images the way to go probably would be maximum entropy methods combined with a noise-model for the sensor. A lot of work but surely an interesting project on its own! Commented Aug 14, 2012 at 3:49

## Version 9 answer - use built-in functionality

The symbols ImageCreateHDRI and ImageToneMapHDRI were present in version 8 but didn't seem to do anything. In version 9 there is functioning code behind them. This is all undocumented, and therefore liable to change before it is officially released, but here is what I've managed to dig up.

ImageCreateHDRI

This function takes a list of images at different exposures and composes them into a single high dynamic range image. The basic usage is simply

ImageCreateHDRI[{image1, image2, ...}, options]


The following default options are defined:

Options[ImageCreateHDRI]
(* {"ExposureTimes" -> {}, "EV" -> {},
"GenerateCameraResponseFunction" -> True, "AlignTranslation" -> False} *)


"ExposureTimes" and "EV" allow you to provide a list of known exposure times or exposure values. If both are left at the default empty list setting, the exposure times will be computed from the images.

"AlignTranslation" can be set to True to align the images prior to combining them. As the name suggests, the alignment is translation only, scale and rotation errors are not handled.

"GenerateCameraResponseFunction" presumably generates a camera response function :-) This takes suboptions, which default to  {"Variance" -> 16., "PredefinedResponse" -> "Linear"}. Other settings for the predefined response curve are "Log10" and "Gamma". I'm not sure what the "Variance" setting does, or what happens if "GenerateCameraResponseFunction" is set to False.

ImageToneMapHDRI

This function takes a single HDR image and applies one of a variety of tone-mapping algorithms, which are specified in a Method option. The default is:

Options[ImageToneMapHDRI]
(* {Method -> "PhotographicToneReproduction"} *)


The available algorithms are:

algos = ImageHDRImageProcessingDump\$ToneMapAlgos
"PhotographicToneReproduction", "ContrastDomain"} *)


Each algorithm has a set of sub-options, the default settings are:

{#, ImageHDRImageProcessingDumpdefaultsuboptions[#]}& /@ algos

(* {{"AdaptiveLog", {"Bias" -> 0.8, "Gamma" -> 1., "Saturation" -> 1.}},
{"Photoreceptor", {"LuminousIntensity" -> 0., "ChromaticAdaptation" -> 0., "LightAdaptation" -> 1., "OverallContrast" -> 0.}},
{"FattalGradientDomain", {"Threshold" -> 0.1, "Order" -> 0.8, "ParamX" -> Automatic, "Saturation" -> 0.8}},
{"PhotographicToneReproduction", {"Key" -> 0.18, "Sharpening" -> 8., "UseLocalOperator" -> False, "NumberOfScales" -> 8, "SmallestScale" -> 1, "LargestScale" -> 43, "Saturation" -> 0.8}},
{"ContrastDomain", {"ContrastScaleFactor" -> 0.3, "UseContrastEqualization" -> True, "Saturation" -> 0.8}}} *)


Not being an expert in HDR tone-mapping algorithms, I can't really comment on what these parameters control, but it is easy enough to play around with the values and see the results.

Example

Here is an example of processing the Stanford Memorial Church images with the new functionality:

(* load the source images *)
files = Reverse@FileNames["memorial*.png"];
images = ImagePad[Import[#], {{-2, -12}, {-35, -30}}] & /@ files;

(* create the HDR image *)
i = ImageCreateHDRI[images];

(* apply a tone-mapping algorithm *)
ImageToneMapHDRI[i, Method -> {"FattalGradientDomain", {"Threshold" -> 0.3}}]


• +1 Simon, persistently great. "FattalGradientDomain" would be my choice too ;-) Commented Nov 3, 2013 at 19:40

Exposure fusion is often a powerful alternative to the reconstruction of HDR images if we are only interested in the tone mapped result. HDR reconstructions are often not possible because exposure times are not available or hard to estimate from the given set of images. Many exposure fusion algorithms are also very fast and therefore suitable for real time applications. So I thought it can't hurt to post another such algorithm.

The idea behind the algorithm is to blend the images in a way so that the level of detail in the final image is maximal.

First each image is weighted by three quality measures, well-exposedness, saturation and contrast. The three weight maps are combined and normalised over all images. Three example weights are shown.

images = Image[Import["..\\memorial00" <> ToString@# <> ".png"], "Real"]
& /@ Range[61, 76];

WellExposedness[img_] := Module[
{ch = ColorSeparate@img},
ImageMultiply @@ (ImageApply[Exp[-(# - 0.5)^2/0.08] &, #] & /@ ch)]
SaturationMeasure[img_] :=
ImageApply[StandardDeviation@# + 10^-10 &, img]
ContrastMeasure[img_] :=
ImageApply[Abs@# + 10^-10 &,
LaplacianFilter[ColorConvert[img, "Grayscale"], 10]]

ImageWeightMap[img_, {a_, b_, c_}] := ImageMultiply @@ {
ImageApply[#^a &, WellExposedness[img]],
ImageApply[#^a &, SaturationMeasure[img]],
ImageApply[#^a &, ContrastMeasure[img]]}

NormalizeWeightmaps[weightmaps_] := Module[
{weightmapSum = ImageApply[1/# &, ImageAdd @@ weightmaps]},
ImageMultiply[weightmapSum, #] & /@ weightmaps]

ImageAssemble[{{WellExposedness@images[[1]],
SaturationMeasure@images[[1]],
ContrastMeasure@images[[1]]}}]


Simple blending of the images using the generated weight maps can introduce ugly boundaries between the image parts as shown in the following example where only 2 images of the set are blended.

SimpleExposureBlend[images_, {exp_, sat_, con_}] := Module[
{imgN = Length@images,
wm = ImageWeightMap[#, {exp, sat, con}] & /@ images,
wmNorm},
wmNorm = NormalizeWeightmaps[wm];
]

SimpleExposureBlend[images, {1, 1, 1}]


This issue can be avoided by using multi-scale pyramid blending to combine the images.

ImagePyramidDepth = Length[NestWhileList[Ceiling[#/2] &, #, Ceiling[#] != 1 &]] - 1 &;
GaussPyramid[img_, depth_, kernel_: GaussianMatrix[2], re_: "Gaussian"] :=
Module[{},NestList[ImageResize[ImageConvolve[#, kernel], Scaled[1/2],
Resampling -> re] &, img, depth]];
LaplacePyramid[img_, depth_, kernel_: GaussianMatrix[2], re_: "Gaussian"] :=
Module[{Gpy = GaussPyramid[img, depth, kernel, re]},
Join[ImageSubtract[#1,
ImageResize[#2, Scaled[2], Resampling -> re]] & @@@
Partition[Gpy, 2, 1], {Last@Gpy}]]
CollapseLaplacePyramid[lapPyr_, re_: "Gaussian"] :=
ImageResize[#2, ImageDimensions@#1, Resampling -> re]] &,
First@lapPyr, Rest@lapPyr]

MultiLevelExposureBlend[images_, {exp_, sat_, con_}] := Module[
{imgN = Length@images,
wm = ImageWeightMap[#, {exp, sat, con}] & /@ images,
pyrDepth = ImagePyramidDepth@Min@ImageDimensions@images[[1]],
wmNorm, gp, lp, blend},
wmNorm = NormalizeWeightmaps[wm];

gp = GaussPyramid[#, pyrDepth] & /@ wmNorm;
lp = LaplacePyramid[#, pyrDepth] & /@ images;
blend = ImageAdd @@ # & /@ Transpose[
ImageMultiply @@@ # & /@ (Transpose[{#1, #2}] & @@@ Transpose[{gp, lp}])];
CollapseLaplacePyramid@blend
]

i1 = MultiLevelExposureBlend[images, {1, 1, 1}];
`