Czech Technical University in Prague

Faculty of Electrical Engineering

DIPLOMA THESIS
Luk as Cerman
High Dynamic Range Images from
Multiple Exposures
Department of Computer Science and Engineering
Supervisor: Prof. Ing. V aclav Hlav ac, CSc.
Study programme: Electrical Engineering & Informatics
Study eld: Computer Science and Engineering
Acknowledgment
I wish to thank my supervisor V aclav Hlav ac for careful leading and valuable
help with writing this thesis. Thanks to Jan

Cech for advice with geometric
image registration. Thanks to my parents for support during my studies.
Thanks to my brother Martin for tolerance during late night computer ses-
sions.
Prohlasen
Prohlasuji, ze jsem svou diplomovou pr aci vypracoval samostatne a pouzil
pouze podklady uvedene v prilozenem seznamu.
Nemam z avazn y duvod proti uzit tohoto skolnho dla ve smyslu 60
Zakona c.121/2000 Sb., o pr avu autorskem, o pr avech souvisejcch s pravem
autorsk ym a o zmene nekter ych zakonu (autorsk y z akon).
Statement
This is to certify that I wrote this thesis on my own and that the references
include all the sources of information I have utilized. I agree with lending
this thesis.
Prague January 26, 2006
Luk as Cerman
Abstract
The limited dynamic range of a camera may be extended by composing dif-
ferently exposed images of the same scene. The nonlinear camera has to be
calibrated radiometrically rst. We present several state of the art calibration
methods and show that the calibration process can be dicult for real cam-
eras. The problems with calibration can be overcomed if linear 12-bit RAW
images are used as oered by modern mid-class and professional cameras. An-
other aspect aecting the radiometric quality of the composed high dynamic
range (HDR) image is the knowledge of exposure times. This knowledge also
helps in registration of dierently exposed hand held shots. We present a
new method for estimating exposure times from the histograms of images
captured using a linear response camera or generated from a RAW image.
The actual process of HDR image composition needs perfectly aligned images
on its input. We present a method for registering dierently exposed images
captured with a hand held camera. The presented registration method en-
ables capture of HDR images without the need of mounting a camera on the
tripod. The methods are implemented in Matlab and tested on synthetic and
real data. The individual methods are incorporated to HDR Library which
is contributed to public domain.
Contents
1 Motivation 10
2 Diploma Thesis Goals 11
2.1 Original Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Revised Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 The Diploma Thesis Structure . . . . . . . . . . . . . . . . . . 12
3 Problem Specication 12
3.1 Combining Dierently Exposed Images . . . . . . . . . . . . . 12
3.2 Image Set Acquisition . . . . . . . . . . . . . . . . . . . . . . 13
3.3 Estimate of the Camera Response Function . . . . . . . . . . . 13
3.4 Image Registration . . . . . . . . . . . . . . . . . . . . . . . . 15
3.5 High Dynamic Range Image Composition . . . . . . . . . . . . 16
3.6 Tonemapping . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
4 State of the Art 17
4.1 Image Set Acquisition . . . . . . . . . . . . . . . . . . . . . . 17
4.2 Estimate of the Camera Response Function . . . . . . . . . . . 18
4.3 Image Registration . . . . . . . . . . . . . . . . . . . . . . . . 20
4.4 High Dynamic Range Image Composition . . . . . . . . . . . . 21
4.5 Tonemapping . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5 The HDR Library 22
5.1 Image Set Acquisition . . . . . . . . . . . . . . . . . . . . . . 23
5.2 Estimate of the Camera Response Function . . . . . . . . . . . 23
5.2.1 Brightness Transfer Function Recovery . . . . . . . . . 23
5.2.2 Equation System for Camera Response Recovery . . . 23
5.2.3 Nonparametric Model . . . . . . . . . . . . . . . . . . 24
5.2.4 Polynomial Model . . . . . . . . . . . . . . . . . . . . . 24
5.2.5 Linear Model . . . . . . . . . . . . . . . . . . . . . . . 25
5.2.6 Split Linear Model . . . . . . . . . . . . . . . . . . . . 25
5.2.7 Color Handling . . . . . . . . . . . . . . . . . . . . . . 27
5.3 Estimate of the Exposure Times . . . . . . . . . . . . . . . . . 27
5.3.1 Exposure Ratio Estimate for Nonlinear Camera . . . . 28
5.3.2 Exposure Ratio Estimate for Linear Camera . . . . . . 28
5.3.3 Limitations of Exposure Times Recovery . . . . . . . . 29
5.4 Image Registration . . . . . . . . . . . . . . . . . . . . . . . . 30
5.5 High Dynamic Range Image Composition . . . . . . . . . . . . 31
5.6 Tonemapping . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
8
6 Results 31
6.1 Use Luminance or Green Channel? . . . . . . . . . . . . . . . 31
6.2 Response Recovery from Synthetic Data . . . . . . . . . . . . 33
6.2.1 Synthetic Data Without Noise . . . . . . . . . . . . . . 33
6.2.2 Noisy Synthetic Data . . . . . . . . . . . . . . . . . . . 33
6.2.3 Summary of Response Recovery Test . . . . . . . . . . 35
6.3 Verication of Exposure Ratio Estimate . . . . . . . . . . . . 36
6.4 Verication of RAW Readout Linearity . . . . . . . . . . . . . 36
6.5 HDR Image with Hand Held Camera . . . . . . . . . . . . . . 38
6.6 Sample Scene . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
7 Future Work 42
8 Conclusions 42
A DVD Contents 44
B HDR Library Documentation 46
References 50
List of Figures
1 The same scene exposed for shadows and for lights. . . . . . . 10
2 Scheme of a camera. . . . . . . . . . . . . . . . . . . . . . . . 13
3 Spatially varying lter. . . . . . . . . . . . . . . . . . . . . . . 18
4 Block scheme of HDR Library design. . . . . . . . . . . . . . . 22
5 Approximation by the split linear model. . . . . . . . . . . . . 26
6 Crosshistograms from green channel and luminosity. . . . . . . 32
7 Response recovery from synthetic data. . . . . . . . . . . . . . 34
8 Exposure ratio estimate. . . . . . . . . . . . . . . . . . . . . . 37
9 Plot proong Nikon D70 RAW image linearity. . . . . . . . . . 37
10 HDR image with hand held camera. . . . . . . . . . . . . . . . 39
11 HDR composition example. . . . . . . . . . . . . . . . . . . . 41
9
Figure 1: The same scene exposed for shadows (a) and for lights (b).
1 Motivation
Since the early years of photography, the photographic process has been
constrained by the dynamic range of a photosensitive material which is able
to record only a limited range of light intensity.
Imagine a high contrast scene containing dark shadows in some parts and
very bright zones in the other parts of the scene. If someone sets the exposure
with the aim to capture details in the dark zones correctly, like in Figure 1a,
then details in the bright zones are lost by overexposure. The material is
saturated and is not able to record dierences between light intensities above
the saturation level. Likewise, if someone exposes to capture bright zones
correctly, like in Figure 1b, then details in the underexposed dark zones are
lost in sensor noise or a lm grain.
Many techniques for setting proper exposure were developed during years,
probably the most popular is the legendary Zone System for black and white
photography by Ansel Adams [1]. This is a sophisticated method of under-
standing and controlling the exposure, development of the negative and the
print process. He divided a spectrum of print values into eleven gray zones
ranging from the dark black with no details and no texture to the pure white
and provided a way to get subjects in the scene to the desired gray zones
in the print. If we leave the world of classic black and white photography
and move towards the modern automatic cameras lm or digital we see
they all have a very sophisticated metering and exposure programs devel-
oped by the analysis of thousands professional photographs.
1
However all
such systems, no matter how sophisticated, do not increase the limited dy-
namic range of the sensing material. They are only trying to exploit it best
by the correct exposure and image processing.
1
For example the latest Nikon metering systems take in account brightness and color
information from 1 005 segments as well as focus informations from the CPU equipped lens
and compares these data to data gathered from more than 30 000 analyzed photographs.
10
Obviously, the limited dynamic range is not a problem of artistic photog-
raphy only, but aects industrial and scientic processes having images as
their inputs. Such processes become more robust if they can utilize the whole
dynamic range of a high contrast scene. If we had wide dynamic range cam-
eras or methods for capturing wide dynamic range scenes with usual cameras
were developed then we would provide more accurate results. That is why
we are interested in the dynamic range widening.
2 Diploma Thesis Goals
2.1 Original Goals
The ocial Czech Technical University document with diploma thesis as-
signment from September 6, 2004 states: Find known approaches by others
for capturing and processing of HDR image in literature, select a suitable
one, potentially improve it, implement it and test it experimentally on real
data. Consider two application scenarios, e.g., (1) use of HDR in digital
2
still photography; (2) processing of a videosequence captured by a camera
positioned in a car having the view similar to a driver.
2.2 Revised Goals
After reading available literature and getting more into the problem, we de-
cided to leave topic of videosequence processing and concern more on digital
still photography. Improving dynamic range of a video camera would require
hardware solution we are not equipped for. Though software methods for
still images may be used for videosequence like Kang et al. demonstrated
in [9] we found them rather impractical. These methods require a lot of
processing time which prevents them from being used in real-time video.
The problem of improving dynamic range of a camera by varying expo-
sures is well described in the literature and many useful approaches have
been developed (see Sections 3 and 4). We decided to implement some of
these methods and make a set of Matlab functions for HDR images process-
ing an HDR Library. This library should cover the whole HDR processing
pipeline image acquisition, response function recovery, registration, HDR
image composition and tonemapping.
3
2
By word digital we mean images captured directly by using electronic light sensors
(CCD, CMOS . . . ) or digitized (scanned) images captured by classical photochemical
methods.
3
Tonemapping is a process of a dynamic range compression which preserves local con-
11
2.3 The Diploma Thesis Structure
Section 1 is a brief introduction to the problem associated with the limited
dynamic range of real word cameras. Section 2 describes goals of this work
and the structure of this document. Detailed description of the problem
will be given in Section 3, where dierent aspects of the problem will be
pointed out. We will mention the need of radiometric camera calibration
and image registration, the process of HDR image composition and the need
of tonemapping. Existing approaches will be briey summarized in Section 4.
Methods we have chosen for our HDR library as well as our implementation
of these methods will be described in Section 5. Our results with sample
images and graphs will be presented in Section 6. Possible improvements of
our approach and known problems will be discussed in Section 7. Conclusions
will be given in Section 8.
3 Problem Specication
Although the high dynamic range image can be captured directly using the
special high dynamic range sensors like Kodak KAC-9628 [18], IMS Chips
HDRC sensors [17], Silicon Vision Products [20], SMaL Camera [21] and
Pixim [19], these solution are quite expensive. However, it can be shown that
the high dynamic range image can be obtained with widely available cameras
by fusing multiple dierently exposed images of the same scene taken with
these low dynamic range cameras.
3.1 Combining Dierently Exposed Images
The idea of fusing images with dierent exposures is not new. In the early
sixties, Charles Wycko developed a wide dynamic range lm composed of
three layers of the same spectral sensitivity but diering in their overall
sensitivity to light. The three layers were colored exactly like in the color
negative lm cyan, magenta and yellow so after printing on the color
paper Wycko created a pseudocolor high dynamic range image where the
individual gray
4
tones were represented by dierent colors.
Wyckos most well-knows pictures are perhaps his motion pictures of
nuclear explosions one could clearly see the faint glow of a bomb just before
it exploded (which would appear as blue, since it only exposed the fast top
trast while reducing global contrast of the image.
4
Note that all three layers have the same spectral sensitivity, so the recorded informa-
tion is achromatic.
12
M
In Camera
I E L
Development Film
Sensor Processing ADC
Electronic
Shutter Lens Scene
Film Camera
Digital Camera
Pixel
Raw Pixel
Values
Values
Figure 2: Scheme of a camera illustrating how is scene radiance transformed
to pixel values.
layer), as well as the details in the highlights of the explosion (which appeared
white since they exposed all 3 layers the details discernible primarily on
account of the slow bottom layer).
5
The same principle may be applied using a standard lm or an image
sensor. Multiple images of the same scene diering only in exposure can be
captured. Each captured image records dierent range of light intensities. If
there is enough overlapping exposures ranging from dark images, recording
information in the bright zones, to light images, recording information in the
dark zones, a high dynamic range image can be reconstructed by properly
combining the information from captured images. This combined image will
represent the whole scene radiance in one oating point map.
3.2 Image Set Acquisition
As discussed in [6], the process of choosing proper exposures and number of
images in the set aects the properties and quality of composed HDR image.
Though it is possible to choose exposures and images number manually, it
may by useful to provide some automatic control over exposure values and
the number of used images in combination with computer driven camera.
3.3 Estimate of the Camera Response Function
If we want to compose a high dynamic range image from dierently exposed
images of the same scene, we must know what the camera measures. Gen-
5
Sentence cited from [11]
13
erally, pixel values are proportional to the scene radiance, transformed by a
nonlinear mapping called the response function. This nonlinearity is induced
by dierent parts of the imaging process. Let us explore the camera scheme
in Figure 2.
As the light characterized by scene radiance L is passing the lens, it
is attenuated by the lens aperture and yields sensor irradiance
E = L

4
_
d
h
_
2
V () , (1)
where h is focal length of the lens, d is diameter of its aperture and is
the angle between the principal ray and the optical axis. V () stands for
vignetting light fallo in the image periphery and can be approximated
by
6
V () = cos
4
. (2)
The total quantity of light accumulated on the image sensor during integra-
tion time t which is controlled by the shutter produce sensor exposure
I = Et . (3)
Typical sensors like CCD or CMOS are designed to produce electrical
signals which are linearly proportional to the sensor exposure up to some
saturation level. Above this level, the sensor is not able to distinguish be-
tween dierent exposure values.
In the case of a lm photography, the lm itself is designed to have a
nonlinear response and the development process causes further nonlinearities.
After exposure, the accumulated charge is converted to integer values
using analog-digital converters connected to each photosensitive element on
the image sensor. The process of digitization brings quantization noise to
recorded data. In addition, the digital values are in the case of the typ-
ical camera further transformed to mimic lm characteristics or scaled by
a gamma function to properly display on CRT monitors. At this point, we
must note that many modern mid-class and professional cameras oer an op-
tion to obtain so called RAW image a pure output from the image sensor.
The RAW image should not suer from any in-camera software postprocess-
ing. The RAW image is aected only by sensor saturation and quantization
noise.
The individual sources of nonlinearity in the imaging process are not
important and the whole process can be represented by one nonlinear function
6
In fact, cos
4
law stands for optical vignetting, typical imaging system exhibits more
types of vignetting, see [27] for details.
14
a camera response function f. We see that measured pixel values M are
proportional to the sensor exposure I scaled by the response function f,
M = f(I) . (4)
We can see now that if we want to fuse dierently exposed images to the
one high dynamic range image then we have to transform all of them to the
same domain irradiance E. There is an assumption that the scene radiance
did not change between the individual exposures. To map pixel values M
to irradiance values E we need to nd inverse function g = f
1
so that the
sensor exposure
I = g(M) . (5)
Knowing the integration time t we obtain irradiance E using expression (3).
In fact, we do not get irradiance in the terms of radiometric quantity
because sensor does not have a uniform spectral response. The sensor mea-
sures quantity of light I() integrated over the full range of wavelengths
weighted by its spectral response R(),
_

I()R()d . (6)
However, a typical sensor will not have the same spectral response as the
human eye, so the measured quantity cannot be luminance either. In the
rest of this text, we will use terms radiance and irradiance, while reminding
the reader to remember that the quantities we are dealing with are weighted
by the spectral response of the imaging system. The color channels can be
treated separately in case of color photography.
3.4 Image Registration
To compose a high dynamic range image from a set of dierently exposed
images of the same scene, we have to relate pixels in the individual exposures
to the corresponding pixels in the other exposures.
When shooting a static scene from a tripod there is no problem all cor-
responding pixels are located at the same coordinates across all exposures.
On the other hand, if we want to take advantage of high dynamic range
imaging with hand-held camera then we need to perform a dense registra-
tion. It means nding pixel-to-pixel correspondences between the individual
exposures.
Unlike a usual registration problem, there is a need to register images
with large exposure dierences and small spacial dierence among captured
images. The spatial dierence is only caused by a little camera shift and
15
shakes between the individual shots. The knowledge of the camera response
and exposure times may help this process by transforming all dierently
exposed images to the same domain irradiance E.
Hand held capturing is also limited by the need of suciently long expo-
sure times preventing the individual images from vibration blur. This lim-
itation may by removed by applying multichannel deconvolution to blurred
images but it is out of the scope of this work.
3.5 High Dynamic Range Image Composition
Denote the measured value of the pixel i observed in the image j by M
ij
,
the sensor irradiance at point i computed from the image j by E
ij
, the cor-
responding pixel i of the high dynamic range map by E
i
and the integration
time of image j by t
j
.
Once the inverse response g is recovered, we can quickly convert pixel
values M
ij
to irradiance values E
ij
, assuming the integration time t
j
is
known. From (3) and (5), we obtain
E
ij
=
g(M
ij
)
t
j
. (7)
Having images registered and properly aligned, we can fuse the recon-
structed irradiance values E
ij
to a high dynamic range map using
E
i
=

j
E
ij
(M
ij
)

j
(M
ij
)
, (8)
where (M
ij
) is some weighting criterion. Reconstructing values E
i
using
irradiance computed from all exposures in the set helps reducing noise in the
nal data. The choice of weighting criterion signicantly aects results and
will be discussed later.
3.6 Tonemapping
It is not possible to show a high dynamic range image directly on the usual
display device. The display devices like CRT/LCD monitors or paper prints
have a very limited dynamic range and cannot present the high contrast
image without a signicant lost of details. The high dynamic range image
needs to be tonemapped rst. More formally, the global contrast of the
high dynamic range image has to be reduced while the local contrast should
be retained as it carries the most valuable information. The example of a
16
high dynamic range image linearly mapped to the device together with a
tonemapped version of the same image can be seen in Section 6.6 where a
sample high contrast scene will be presented.
Image size and processing speed are another factors advancing the use of
tonemapped images over the high dynamic range images. The high dynamic
range image is usually represented by a oating point map where each pixel
consumes 12 bytes when represented by a single precision oating point
numbers or even 24 bytes in the case of double precision representation.
The six megapixel image would consume as much as 432MB of memory in
the double precision format. Additionally, the numeric calculations with
oating point numbers are often more time consuming than the xed point
calculations. Unlike the high dynamic range images, the tonemapped images
may be represented by integers without signicant loss of precision.
Another reason for the use of tonemapped images is the behavior of some
image processing lters which might not work correctly when applied to the
very contrast high dynamic range image.
4 State of the Art
A lot of theoretical work in the problem of high dynamic range imaging has
been published during the past fteen years. In the following sections, we
will overview the most interesting ones which are related to our work.
4.1 Image Set Acquisition
Debevec and Malik discuss [3] the number of images required. They consider
two aspects of the process a response function recovery and a high dynamic
range image composition.
Recovering the response function requires a minimum of two photographs,
provided they contain suciently many dierent radiance values.
The composition of high dynamic range image requires so many images
that every pixel falls at least in one image in the working range of the imaging
system.
7
More precisely, if R is the dynamic range of the scene and F is
the width of dynamic range corresponding to the working range of imaging
system then the minimum number of images is
_
R
F
_
Stumpfel et al. described [22] a simple program which automatically
captures the whole set with the computer driven camera. The rst image
7
The working range of the imaging system corresponds to the middle section of the
response curve, where small changes in exposure cause signicantly large changes in a
pixel value.
17
e e
e
3
e
2
4 1
Figure 3: Left side illustrates an example of spatially varying lter attached
to the image sensor. Every little square is the neutral density lter geomet-
rically aligned to underlaying sensor pixel. The arrangement of the lter
groups the sensor pixels to the exposure quaternions which are illustrated
on the right side. Every pixel e
i
receives dierent exposure according to the
density of lters in the quaternion.
is captured with the automatic exposure set. The image is analyzed for
presence of underexposed or overexposed pixels and it is decided whether
to take additional exposures. Further images are taken in manual mode
with exposures diering by a factor of two. In manual mode, aperture and
integration time are controlled by the connected computer.
Nayar and Mitsunaga describe [13] a quite dierent approach. They at-
tached neutral density lter mosaic (see Figure 3) to the image sensor allow-
ing to capture multiple exposures simultaneously at the expense of slightly
reduced spatial resolution. The spatial resolution is not reduced so dramat-
ically as it might look at the rst glance. Mid-tone radiances fall in the
working range at all pixels in the exposure quaternion so they could be re-
constructed well. The only data missing is of the saturated fast pixels in
the highly exposed quaternions and underexposed slow pixels of the low ex-
posed quaternions. This method does not require image registration, which
tends to be computationally very expensive. Thus the approach might be
interesting for generating high dynamic range video.
4.2 Estimate of the Camera Response Function
The camera response function can be estimated by taking a picture of the
uniformly illuminated chart containing patches of known reectance, such
as the Gretag Macbeth chart [2]. However this process is quite complicated
and can be performed only in laboratory conditions which are not always
accessible. In addition, the camera response may be aected by temperature
changes requiring frequent recalibration.
18
Fortunately, it has been shown that a set of dierently exposed images
contains usually enough information to recover a response using the images
themselves [3], [12], [23], [10] and [6]. If there is not enough information
there is not suciently many dierent radiance values in the images and we
cannot make more additional exposures then we may overcome this problem
by multiply shooting
8
a chart containing patches of diverse reectance.
9
The fundamental equation for all chartless recovery methods comes from
terms (3) and (5). Suppose we take two images A and B of the same scene
with dierent exposure times t
A
and t
B
.
10
If the sensor irradiance had
not changed between the two shots then we can write using (3)
I
A
t
A
=
I
B
t
B
, (9)
where I
A
is sensor exposure at a point in image A and I
B
is sensor exposure at
the corresponding point in image B. Taking exposure ratio k = t
B
/t
A
we
can express the relation between the images using I
B
= kI
A
. Now if we come
to measured pixel values M
A
and M
B
we get the basic equation using (5)
g(M
B
) = kg(M
A
) . (10)
Recovering g from (10) is a dicult task. However we can place certain
restrictions on g. The minimum irradiance is obviously 0 and will produce
no response of imaging system
11
, hence
g(0) = 0 . (11)
The maximum irradiance is an unrecoverable parameter. We get rid of it by
normalizing domain of g to go from 0 to 1 and setting
g(1) = 1 , (12)
while having image intensities normalized to the interval 0, 1. The mono-
tonicity of monotonicity of g may be supposed to.
8
Each time with dierent exposure.
9
Note that unlike chart recovery method described in previous paragraph, this method
does not require the use of uniformly lit chart nor the knowledge of each patch reectance.
The only assumption is that it contains enough patches with dierent reectances.
10
We do not consider changes in aperture because of problems with vignetting, though
our equations may be easily modied to enable changes in aperture.
11
A real imaging system exhibits some little but nonzero response called dark current.
This eect may by reduced by taking an image a dark frame with a lens covered and
subtracting it from an ordinary image.
19
Grossberg and Nayar showed shown in [7] that recovering g from (10)
exhibits a fractal ambiguity. To break this ambiguity, further a priori as-
sumptions on g are required. Mann and Pickard [11] assumed g is a gamma
function. Debevec and Malik [3] used a smoothness constraint to break this
ambiguity and recover response using nonparametric model of g, sampled
at certain values, and represent it by a vector. Mitsunaga and Nayar [12]
assumed that the response can be modeled by a low degree polynomial. Fur-
ther, they are able to start from rough estimates of exposure ratios k and by
iterative approach recover both g and k.
The previous methods used pixel values (M
A
, M
B
) selected from images
at corresponding points to constrain the solution of g. Selection of these
points had a signicant impact on the quality of the solution. Mann et al.
showed in [10] that all information for response recovery available in the
image pair can be represented by the brightness transfer function T.
12
The
function describes how the brightness transfers from one image to another,
M
B
= T(M
A
). Having brightness transfer function, equation (10) can be
rewritten to
g(T(M)) = kg(M) . (13)
Additionally we can construct brightness transfer function T
1
, which de-
scribes the transfer of pixel values from image B to image A, T
1
(I
B
) = I
A
.
Having T
1
we get another constraint using (10)
g(M) = kg(T
1
(M)) . (14)
Grossberg and Nayar showed in [7] that brightness transfer function can
be recovered using the images histograms only, requiring no registration.
They modied the solution of Mitsunaga and Nayar [12] to constraint polyno-
mial by equations (13) and (14) utilizing brightness transfer function instead
the values of some selected pixels.
4.3 Image Registration
Techniques used to register dierently exposed images before a high dynamic
range image composition are similar to those used in mosaic and panorama
stitching applications. Schechner and Nayar [16] mentioned traditional reg-
istration methods based on optimizing the sum of square or absolute dier-
ences, correlation or mutual information between frames, over the general
motion parameters.
12
He calls it comparagraph, but we will keep to the terminology of Grossberg and Nayar.
20
Kang et al. [9] implemented a sophisticated registration algorithm which
was used in their high dynamic range video application. Before registra-
tion, they boosted the brightness of darker frames to bring all frames in the
video stream to the same brightness level. To register exposure compensated
frames, they used motion modeling together with hierarchical homography.
The motion model is used to compensate camera movement while the hier-
archical homography is used to compensate for camera shake.
Ward described [26] a fast image registration method for high dynamic
range image composition. The median thresholded images are registered
using his approach. The images thresholded by their median are invariant to
the changes in exposure. The registration process is performed using the fast
bit manipulation routines applied on an image pyramid. However, to keep
this method fast it is possible to compensate for image shifts only and not
for rotation or general homography which would require image warping. For
the same reason, it is not possible to register shifts at subpixel precision.
4.4 High Dynamic Range Image Composition
All solutions we have seen base the high dynamic range image composition on
the equation (8). However they dier in the choice of the weighting criterion
(M) used in (8). Mann and Pickard [11], Mann et al. [10] use derivative of
the camera response g. Debevec and Malik [3] use a hat function
(M) =
_
M for M
1
2
1 M for M >
1
2
. (15)
Finally Mitsunaga and Nayar [12] base their weighting criterion on the signal
to noise ratio
SNR = I
dM
dI
1

N
(M)
=
g(M)

N
(M)g

(M)
, (16)
assuming noise
N
is independent of the measurement pixel value, they dene
weighting function as
(M) =
g(M)
g

(M)
. (17)
4.5 Tonemapping
Tonemapping algorithms belong to computer graphics rather than computer
vision. Many useful algorithms with dierent properties were described dur-
ing the past decade. Very comprehensive collection of the tonemapping al-
gorithms is available at www pages of Martin

Cadk [25]. The algorithms
21
Exposures
Registration
Tonemapper
HDR
Stitch
Response
Recovery
Acquisition
Image
Camera
Output
Data Flow
Module
Exposures
Images
Images
Registered Images
HDR Image
HDR Image
LDR Image
Response
Optimized Exposures
Response
Images
Exposures
Images
Radiance
Figure 4: Block scheme of HDR Library design.
are presented there with samples and references to their detailed description
and in case of some algorithms the implementation source code is provided.
We found the most interesting works of Thumblin and Turk LCIS [24],
Durand and Dorsey Fast Bilateral Filtering [4], Reinhard et al. Pho-
tographic Tone Reproduction for Digital Images [14] and Fattal et al.
Gradient Domain High Dynamic Range Compression [5].
5 The HDR Library
We implemented the experimental library to test dierent methods. The
basic structure of HDR Library is illustrated in Figure 4. The individual
blocks are designed to be more or less independent of the others so one
could easily replace a specic block with his own implementation. The only
requirement is to keep the data formats on the block interfaces.
Additionally one might use only some bocks if he does not need the func-
tionality of the others, for example he needs chartless camera calibration
only. Eventually our blocks may be incorporated into another image pro-
cessing system.
22
All blocks were implemented as a set of Matlab functions, which may in
case of need, use some external programs called from Matlab.
5.1 Image Set Acquisition
Currently, we acquire source images manually without any sort of automation
support. We vary the exposure using the manual settings of the camera or
using a bracketing function the ability of the camera to take three images
with a xed, adjustable exposure dierence.
5.2 Estimate of the Camera Response Function
For a comparison and a performance testing, we decided to implement more
than one response recovery method. We implemented both versions of poly-
nomial approximation described in Section 4.2 the original one by Mit-
sunaga and Nayar [12] and the histogram based modication by Grossberg
and Nayar [7]. We also implemented the original nonparametric method of
Debevec and Malik [3] according to their sample implementation, which was
published in the appendix of their article.
Additionally, we modied the original method of Debevec and Malik [3]
to recover a nonparametric model of g from image histograms. We followed
the idea of Grossberg and Nayar [7] originally used for the polynomial model
of Mitsunaga and Nayar.
5.2.1 Brightness Transfer Function Recovery
The brightness transfer function represents all information available in the
image pair required for a response recovery. With Grossberg and Nayar [7],
we recover brightness transfer function T from image histograms using
T(M) = H
1
B
(H
A
(M)) , (18)
where H
A
is cumulative histogram of image A, H
1
B
is the inversion of cu-
mulative histogram of image B and M ranges across all levels corresponding
to all histogram bins. Number of histogram bins which is simultaneously a
number of T samples is the only adjustable parameter of this algorithm.
H
1
B
(n) is inverted in the sense: What will be the brightness value of the
most bright pixel in a set containing n darkest pixels of the image B?
5.2.2 Equation System for Camera Response Recovery
We order images in a set by increasing exposure. Consequently we get N
pairs of adjacent images having exposure ratios k
j
, where j goes from 1 to N.
23
Having reconstructed a brightness transfer functions T
j
and T
1
j
for each
image pair j at Z discrete brightness levels, we use (13) and (14) to form a
system of 2NZ equations
g(T
j
(M
z
))
_
h
j
(M
z
) = k
j
g(M
z
)
_
h
j
(M
z
)
g(M
z
)
_
h
j+1
(M
z
) = k
j
g(T
1
j
(M
z
))
_
h
j+1
(M
z
) , (19)
where M
z
goes across all brightness levels used for T reconstruction, h
i
is
the histogram of the image i. We use square roots to weigh the least square
sense solution of the equations system (19).
5.2.3 Nonparametric Model
Using the idea of Debevec and Malik [3], we may recover camera response g
at a xed number of points (e.g. Z) and treat is as a vector g
s
. The value of
g at z-th brightness level may be calculated by
g(M
z
) = g
s
s(z) , (20)
where s(z) is Z-dimensional vector having nulled all components but the z-th
which is set to 1. Setting (20) to (19), we get a linear of linear equations. To
break fractal ambiguity, we force smoothness of g by extending the equations
system with equations
(g(M
z1
) 2g(M
z
) +g(M
z+1
)) = 0 , (21)
where z goes from 2 to Z1 and weights the smoothness terms (21) relative
to the data tting terms (19). Additionally we include the normalizing term
g(1)

P =

P (22)
to the nal system of equations to x value of g(1). P is the number of pixels
in any image from the set. We rewrite system to the matrix form and solve
it in Matlab using well known x = A\b.
5.2.4 Polynomial Model
Following Mitsunaga and Nayar [12] with Grossberg and Nayar [7] we imple-
mented an approximation of camera response g by order Q polynomial
g(M) =
Q

q=0
c
q
M
q
. (23)
24
From constraint (11) we get
c
0
= 0 , (24)
from constraint (12) and terms (23), (24) we get
c
1
= 1
Q

q=2
c
q
. (25)
The polynomial model (23) can be rewritten to the following form using the
terms (24) and (25)
g(M) = M +
Q

q=2
c
q
(M
q
M) . (26)
Setting (26) to (19), we get an overdetermined system of 2NZ equations
linear in polynomial coecients c
q
. We solve it using Matlab after rewriting
it in the matrix form.
5.2.5 Linear Model
We also implemented the linear model approximation, which is suitable to
be used together with RAW images produced by modern cameras. Due to
constraints (11) and (12), there is no need for parameter estimate and g may
be represented by a linear function going from (0, 0) to (1, 1)
g(M) = M . (27)
5.2.6 Split Linear Model
Sometimes it is useful to model g with a split linear function
g(M) =
_
a
1
M +b
1
for M < d
a
1
M +b
1
for M d
, (28)
where d (0, 1). From (11) and (12) we get
b
1
= 0 (29)
b
2
= (1 a
2
) . (30)
Fixing d at some value, we may include term (28) to equation (19) and
solve the linear system for parameters a
1
and a
2
. We have implemented an
iterative search for parameters a
1
, a
2
and d with complexity proportional to
desired precision of d. It can be described by the following pseudocode. See
Figure 5 for illustration.
25
I
1
y
2
1
y
y
2
d d M 1
I
0
1
d 1 M d
1
0
y
Figure 5: Illustration to the iterative search for split line model parameters.
On both graphs, the two parts of the split line do not meet at the same
point. The graph on the left side illustrates the situation where d

estimate
is smaller than the real d, we see that the left part of the split line has a
greater value in point d

than the right part. The graph on the right side

illustrates the situation where d

estimate is greater than the real d. We see

that the left part of the split line has a smaller value in point d

than the
right part.
Split Line Solver
Input: Precision
Output: a
1
, a
2
and d
begin
l = 0 Lower bound on d.
u = 1 Upper bound on d.
repeat Precision times
d = (l +u)/2
(a
1
, a
2
) = SolveForLineSlopes(d)
y
1
= a
1
d
y
2
= a
2
(d 1) + 1
if (y
1
= y
2
) then
return Solution found prematurely.
elif (y
1
> y
2
)
l = d
else
u = d
end
end
end
26
5.2.7 Color Handling
The chartless recovery methods recover camera response function up to the
unknown scale. To remove the unknown scale parameter we x the values by
normalizing the camera response function to go from 0 to 1. In case of color
images, we recover the response function of each color channel separately.
To produce a color balanced high dynamic range image we need to relate
the individual channels each to the other. We manage this by requiring that
mid-grey pixels in the image space are achromatic in the recovered radiance
space
g
r
(0.5) = g
b
(0.5) = g
g
(0.5) . (31)
5.3 Estimate of the Exposure Times
Exposure times reported by the camera on the status display or in the EXIF
data are not the times actually used as noted by Debevec and Malik [3] or
Mitsunaga and Nayar [12]. The typical camera provides the exposure times
usually rounded to one half stop or one third stop to conform with photo-
graphic tradition in marking exposure times. For example, if the camera
claims 1/30s or 1/60s was used then the real value could be 1/32s or 1/64s.
Even worse, some cheap cameras may have a signicant variance in repro-
ducing the same exposure setting. The knowledge of the real exposure times
positively aects the radiometric quality of the composed high dynamic range
image.
While the produced high dynamic range image is normalized to have
brightness values ranging from 0 to 1 there is no need to know the absolute
value of the exposure times. It is sucient to estimate the relative exposures
against the darkest image. We set the exposure time of the darkest image
t
1
to be unit,
t
1
= 1 . (32)
The relative exposures of the other images can be computed from the expo-
sure ratios k
j
of the consecutive image pairs using the formula
t
j
= k
j1
t
j1
. (33)
Recovering exposure times from images themselves can also compensate
for uniform changes in illumination between the individual shots provided
the change is uniform across the whole scene. The dierence in the scene
illumination is contained in the recovered exposure times.
27
5.3.1 Exposure Ratio Estimate for Nonlinear Camera
Mitsunaga and Nayar showed [12] that the exposure ratios can be estimated
together with the camera response function using an iterative scheme. We in-
corporated their solution to our HDR library to improve quality of generated
high dynamic range images. In the iterative process, the current exposure ra-
tios k
(i1)
j
are used to compute next camera response model g
(i)
. This model
is used to update exposure ratio estimates using (13) and (14)
k
(i)
j
=
Z

z=1
_

1
(M
z
)
g
(i)
(T(Mz))
g
(i)
(Mz)
+
2
(M
z
)
g
(i)
(Mz)
g
(i)
(T
1
(Mz))
_
Z

z=1

1
(M
z
) +
2
(M
z
)
, (34)
where

1
(M
z
) =
3
_
(M
z
)(T(M
z
))h
j
(M
z
) , (35)

2
(M
z
) =
3
_
(M
z
)(T
1
(M
z
))h
j+1
(M
z
) . (36)
The initial exposure ratios k
(0)
j
are computed from the exposure times re-
ported by the camera.
The combined criterion puts together our condence in the components of
the weighted expression. We use a geometric mean to combine the individual
criterions. Such combination suppresses the dierence in absolute values
of the criterions. Criterion (M
z
) expresses our condence in g in point
M
z
whereas histograms h
i
express our condence in recovered brightness
transfer functions at respective points. The actual choice of criterion (M
z
)
is discussed later in Section 5.5.
5.3.2 Exposure Ratio Estimate for Linear Camera
When capturing with a linear response camera there is no need for response
recovery. Hence, we use a dierent approach to estimate the exposure ratio
from the captured images. Unlike the iterative scheme of Mitsunaga and
Nayar [12] our method does not require the knowledge of rough exposure
estimates. It only needs the brightness transfer function which can be recov-
ered using the image histograms only. The recovery of brightness transfer
function is described in Section 5.2.1.
In our approach, we t the brightness transfer function using a linear
function. The derivative of this linear function is the unknown exposureratio.
28
If the camera has a linear response then the brightness values of the image
pair are related by the equations
k
j
M
A
= T
j
(M
A
)
M
B
= k
j
T
1
j
(M
B
) , (37)
where k
j
is the exposure ratio of j-th image pair. T
j
is the brightness transfer
function relating brightness values M
A
of the image j to the brightness values
of the image j + 1. T
1
j
is the inverse brightness transfer function relating
brightness values M
B
of the image j + 1 to the brightness values of the
image j.
From (37), M
A
{M: T
j
(M) < U} and M
B
{M: M < U}, where U is
some saturation level, we can construct system S
0
of equations linear in k
j
k
j
M
A
_
h
j
(M
A
) = T
j
(M
A
)
_
h
j
(M
A
)
M
B
_
h
j+1
(M
B
) = k
j
T
1
j
(M
B
)
_
h
j+1
(M
B
) . (38)
We use square roots of image histograms h
j
, h
j+1
to weigh the solution of
the equations system in the least square sense. We solve the linear system
using Matlab. Because of noise in the data, there are outliers aecting the
least square t. An iterative scheme is applied to remove the outliers.
First, the initial ratio k
(0)
j
is estimated using the system S
0
. Next, the
following steps are repeated until the change of k
j
is smaller than the required
precision.
1. The average square error V
i
of the solution S
i
is computed.
2. Using current estimate k
(i)
j
, the square errors of all equations in the
system S
0
are computed.
3. Equations having errors smaller than 3V
i
the inliers are incorporated
in the next linear system S
i+1
.
4. The system S
i+1
is used to estimate the next exposure ratio k
(i+1)
j
.
5.3.3 Limitations of Exposure Times Recovery
The process of exposure times recovery requires that the images suciently
overlap in the brightness domain. If there is no information to relate bright-
ness values in one image to the brightness values in the other image then the
exposure ratio of the two images cannot be estimated.
29
5.4 Image Registration
The HDR library provides the image registration module. The image reg-
istration algorithm is necessary in the case of image capturing with a hand
held camera is considered.
Before the registration itself, the camera response function and exposure
times are estimated from the histograms of unregistered images. This is done
using methods described in Section 5.2 and Section 5.3.
Only consecutive image pairs are registered. To warp disjunct images a
composed transformation is created from the transformations coupling the
registered pairs. At the start, each image pair is linearized using the esti-
mated camera response function. Next, the darker image is multiplied by
the exposure ratio of the two images. Values above some saturation level
are discarded in both images. Similarly values below some noise level are
discarded in both images.
The rough estimate of the image shift (x
0
, y
0
) at a pixel precision is
found using correlation in Fourier space
(x
0
, y
0
) = max(F
1
(F(I
1
) F(I
2
))) , (39)
where I
1
and I
2
are the two images in the pair and F is the complex Fourier
spectrum.
The estimates (x
0
, y
0
) and angle
0
= 0 are used to initialize the local
search for the image shift (x, y) and rotation with subpixel precision.
The search is performed by the Matlab function fminsearch which uses the
iterative simplex search method. The optimized criterion is the square of
normalized dierence summed across all corresponding pixel candidates

i
=

x,y
_
I
1
(x, y) W
2
(x, y)
n(x, y)
_
2
, (40)
where W
2
is image I
2
warped using current (x
i
, y
i
) and angle
i
= 0 in
the iteration i. n(x, y) is the normalizing function balancing the importance
of bright and dark pixels. We used the average value of the potentially
corresponding pixels as the normalization term
n(x, y) =
_
_
_
I
1
(x,y)+W
2
(x,y)
2
for
I
1
(x,y)+W
2
(x,y)
2
t
t for
I
1
(x,y)+W
2
(x,y)
2
< t
, (41)
where t is normalization threshold preventing the boost of importance of the
very dark noisy pixels.
30
5.5 High Dynamic Range Image Composition
We create high dynamic range images according to expression (8). For non-
parametric and polynomial model we use SNR (17) as the weighting criterion
. The linear and split linear responses are weighted by a hat function.
Composed HDR image may be saved to disk as a 32-bit oating point
TIFF. Because of Matlab is unable to handle 32-bit oating point TIFFs
we had to implement it by ourself using raw2tiff utility which is a part of
popular libtiff library. The extern raw2tiff program is called from our
Matlab function responsible for high dynamic range image saving. First the
oating point high dynamic range map is saved to the temporary le as a
raw data pixel by pixel, channel by channel. The raw le is then converted
to TIFF using the raw2tiff utility.
5.6 Tonemapping
We used the iCAM algorithm [8] in the tonemapping module. We have
chosen this algorithm because its implementation was freely available in the
form of ready to use collection of Matlab functions. If someone will not
be satised with results of this method he/she might implement one of the
methods mentioned in Section 4.5 and replace iCAM module with his/her
own implementation.
6 Results
The individual parts of the HDR library were tested to verify their func-
tionality. In this section, the realized experiments will be described and the
attained results will be presented.
6.1 Use Luminance or Green Channel?
For simplicity, we used grayscale images in our very rst experiments. By
properly combining all three color channels we converted all the input color
images to luminance images. We tried to estimate the camera response func-
tion from these luminance images but we did not get satisfying results. Even
for RAW image data, which should be linear, we got nonlinear response. Af-
ter spending some time with experiments and tries to improve implemented
estimation algorithms, the estimated response was still nonlinear.
Finally we realized that the problem is not in the estimation algorithm nor
in nonlinearity of real RAW data, but in the use of the luminance combined
31
M
A
M
B
0 0.25 0.5 0.75 1
0
0.25
0.5
0.75
1
M
A
M
B
0 0.25 0.5 0.75 1
0
0.25
0.5
0.75
1
Figure 6: Crosshistograms of image A against B computed from the green
channel (left) and luminosity (right). Red line shows the recovered brightness
transfer function.
from the three dierent color channels. In Figure 6, one can see crosshis-
tograms of two RAW images A and B of the same scene photographed with
dierent exposures. It shows how the brightness values in less exposed im-
age M
A
are related to the brightness values in the more exposed image M
B
.
The crosshistogram on the left side was computed from the green channel
while the crosshistogram on the right side was computed from the luminance
image. The red lines show the brightness transfer function recovered from
the image histograms according to description given in Section 5.2.1. One
can see that the transfer of the brightness values in the green channel is
strictly linear. On the other side the transfer of brightness values between
the luminance images exhibits some sort of piecewise linear characteristic.
There might be applications requiring grayscale high dynamic range im-
ages but the source photographs were captured using a color camera. In
such cases, we would recommend to use data from the green channel only,
according to our experience described above. If the information from the all
color channels is required then the color high dynamic range image might
be computed from the source images and then converted to luminance using
usual algorithm.
32
6.2 Response Recovery from Synthetic Data
To verify the usability of response recovery methods described in articles [3],
[12] and [7] we tested our implementation of these methods on synthetic
images created with the known response function.
6.2.1 Synthetic Data Without Noise
We created six virtual grayscale images generated from a green channel of
a high dynamic range image. For this purpose, we used the high dynamic
range image of the sample scene which will be presented later in Section 6.6.
The individual images were generated using formula
M
ij
= g
1
(s(E
i
t
j
)) , (42)
where E
i
is normalized irradiance value of pixel i. t
j
is exposure time
of image j, we used values {1, 4, 16, 64, 256, 1024} which provide two stops
dierence between the individual exposures. s(I) is saturation function
s(I) =
_
I for I 1
1 for I > 1
. (43)
g
1
is the known camera response function. We used a 3rd order polynomial.
M
ij
is the generated brightness value of pixel i in the virtual image j. Note
that the rst image is the exact copy of the original high dynamic range
image.
The results of dierent response recovery algorithms are plotted in Figure
7b. The true response function g actually used is painted by the green line.
The original method by Debevec and Malik [3] is shown blue, the modied
version of their algorithm utilizing only image histograms is displayed in red.
We used 10th order polynomial to estimate camera response using the
polynomial model of Mitsunaga and Nayar [12]. The original method of
Mitsunaga and Nayar is painted in cyan. The histogram based version of
this method described by Grossberg and Nayar [7] is shown in black.
It is obvious that the results of all methods are very similar and close to
the true response g when applied to data without noise and quantization.
6.2.2 Noisy Synthetic Data
The previous experiment veried the ability of all tested algorithms to recover
camera response, but only for images without noise and quantization errors.
However, the real situation is rather dierent. To make our experiment more
natural, the images in the testing series were quantized to 256 discrete levels.
33
(a)
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
N
o
r
m
a
l
i
z
e
d

I
r
r
a
d
i
a
n
c
e
Image Brightness
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
N
o
r
m
a
l
i
z
e
d

I
r
r
a
d
i
a
n
c
e
Image Brightness
(b) (c)
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
N
o
r
m
a
l
i
z
e
d

I
r
r
a
d
i
a
n
c
e
Image Brightness
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
N
o
r
m
a
l
i
z
e
d

I
r
r
a
d
i
a
n
c
e
Image Brightness
(d) (e)
Figure 7: (a) Synthetic nonlinear images. Plots: Green line represents inverse
of camera response function actually used. The other lines represent curves
estimated by using dierent algorithms non parametric model of Debevec
and Malik (blue), histogram based version of the same algorithm (red), poly-
nomial model by Mitsunaga and Nayar (cyan) and histogram based version
of this algorithm (black). (b) response estimated from data without noise
and quantization. (c)(e) response estimated from noisy and quantized data.
34
A normally distributed noise with zero mean and unit variance was added to
each pixel value. Noisy values smaller than 0 were clamped to 0. Similarly,
values over 1 were clamped to 1.
We used the same algorithms to estimate the response function from
the noisy data like in the previous experiment. The resulting estimates are
plotted in Figure 7c. It can be seen that all methods used failed to recover
the true camera response shown by the green line.
Duiker, Hawkins, and Debevec provided the program mkhdr based on the
algorithm by Debevec and Malik [3]. The program works quite well even if
quantization errors and noise are present. When we closely examined the
response function generated by mkhdr, we learned that it is cut above some
unknown level. The level was dierent for each color channel and varied
between dierent image series. According to this new knowledge, we decided
to recover response function only from values in the range 0, 0.98. The value
of response function at the interval (0.98, 1 was set to be constant value of
g(0.98). The results of such an improved algorithm were very similar to those
of mkhdr utility. The results of this improved algorithm with noisy synthetic
data from previous experiment are plotted in Figure 7d. We see that the two
methods based on the nonparametric model of Debevec and Malik provided
results that were very close to the real response function even for noisy and
quantized data. However, the polynomial approximation still failed to recover
the true camera response.
In the last experiment, we tried to estimate the response function using
a polynomial of the same order (third) as the real response function. It can
be seen in Figure 7e that the third order polynomial models tted well to
the real response function.
6.2.3 Summary of Response Recovery Test
It is obvious from the experiments presented above that the naive implemen-
tation of the algorithms described in articles [3], [12] and [7] does not give
reasonable results. There are many unknown aspects aecting the quality of
results provided by these algorithms.
We showed that the method of Debevec and Malik [3] can be strengthened
by ignoring the brightness values greater than some value near the saturation
point. Additionally, we showed that the polynomial based methods are, in
practice, useful to model polynomial response of a known order only.
Regarding the problems associated with the camera response recovery,
we would recommend to use a linear camera or a modern camera with the
RAW output available to capture the high dynamic range images. The usage
of RAW data has one additional advantage the 12-bit precision.
35
6.3 Verication of Exposure Ratio Estimate
The test using both the synthetic and the real images was conducted to verify
our algorithm estimating exposure ratio of the two dierently exposed images
captured with a linear response camera. The estimation algorithm itself was
described in section Section 5.3.2.
For the synthetic test, we created six virtual grayscale photographs from
a high dynamic range image using the similar approach to the one described
in Section 6.2. We used a linear response function instead of a third order
polynomial. Similarly we quantized brightness values to 256 discrete levels
and added normally distributed noise with mean zero and unit variance.
The left plot in Figure 8 illustrates the estimate of the exposure ratio for
the two darkest synthetic images. The darkest images were chosen because
they suer from noise more than brighter images. The initial guess was
3.1378. After seven iterations removing outliers the estimate converged to
the value 4.0178. The real exposure ratio used to generate the images was 4.
The right plot in Figure 8 illustrates the estimate of the exposure ratio
for the second and third image from the sample set presented in Section
6.6. The initial guess was 2.1406. It converged to the value 2.1868 after
four iterations. The ratio computed from the exposure times reported by the
camera was 1.8750. In Figure 8, it can be seen that the estimated value ts
the brightness transfer function much better than the value computed from
camera data.
6.4 Verication of RAW Readout Linearity
To overcome problems associated with the estimation of the camera response
function described in Section 6.2, one might use the RAW output available
in many modern digital cameras. The RAW image represents data read
directly from the image sensor which is, in most cases, designed to have
a strictly linear characteristics. Having the linear data, the problematical
calibration step may be skipped. In addition, using RAW images has one
extra advantage over JPEG shots the 12-bit precision of recorded brightness
values. The four extra bits in every pixel value result in a lesser quantization
error and a wider dynamic range in the input images.
We used images captured by Nikon D70 single-lens-reex camera in our
experiments. To prove that the RAW images provided by the Nikon D70
camera are truly linear we performed a simple experiment. We uniformly
lit a sheet of clear paper by a lamp having natural daylight spectrum. We
mounted the camera on the tripod and adjusted its position in front of the
36
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
0
0.2
0.4
0.6
0.8
1
M
A
M
B
0 0.1 0.2 0.3 0.4 0.5 0.6
0
0.2
0.4
0.6
0.8
1
M
A
M
B
Figure 8: Plots show a linear t to the blue painted brightness transfer
function (BTF). The red line shows initial estimate, the green line (it is
covered by the black line in the left plot) shows the result after a nal iteration
with outliers removed. The left plot is the BTF of the two darkest synthetic
images. The real BTF used to generate these images is shown in black. The
right plot is a real world example. The black line expresses BTF as it would
correspond to exposure times reported by the camera.
0 5 10 15 20 25 30
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Relative Exposure
I
m
a
g
e

B
r
i
g
t
h
n
e
s
s
Figure 9: Plot of average brightness intensity of each color channel computed
from 15 dierently exposed images of a uniformly lit white paper plotted
against the relative exposure of each respective image.
37
paper sheet so that the lens axes was almost
13
perpendicular to sheet plane
and the photographed paper covered the whole image frame.
We took a fteen images diering in exposure by one third of stop the
desired relative exposure dierence of two consecutive images was
3

2. The
images taken were ranging from the deep dark one to the almost saturated
one. The real exposure times may dier from the values reported by the cam-
era. We used algorithm described in Section 5.3.2 to estimate the relative
exposure times which were actually used. We averaged the measured bright-
ness values in each color channel of the every captured image. We set the
averaged brightness values to the relation with the relative exposure times.
The resulting graph expressing the dependence of the image brightness on
the exposure time is shown in Figure 9. It can be seen that characteristics
of the dependence is perfect linear.
6.5 HDR Image with Hand Held Camera
Figure 10 demonstrates the performance of the implemented registration al-
gorithm. There is a tonemapped version of the high dynamic range image on
the left side of Figure 10. The image was stitched from three photographs
captured with a hand held camera. Exposure time was set so that the dif-
ference between the individual images was two stops. The photographs were
captured at a wide end of the zoom lens a 28mm focal length in 35mm
camera equivalent. The use of a short focal length makes the registration
process dicult due to the wide eld of view and the lens distortion.
The right side of Figure 10 shows enlarged details of dierent images of
the photographed scene. The top three images are details from the middle
section of the scene while the bottom three images are details from the lower
left corner. First image in each series was taken from the second source
photograph and provides a sharpness comparison with the high dynamic
range stitches. Middle image in each series was taken from the high dynamic
range image stitched from photographs registered by modeling shift only.
The registration process was performed using the correlation in the Fourier
domain at a pixel resolution. Bottom image in each series was taken from
the high dynamic range image stitched from photographs registered by shift
and rotation modeling at a subpixel resolution. This was done by a nonlinear
optimization using a simplex method.
From the details shown, it can be seen that the registration by modeling
shift at a pixel resolution is not sucient. On the other side, the images
13
We could not place it directly in front of the paper to prevent shadows casted by the
camera from obscuring the paper sheet.
38
Figure 10: The HDR image taken with a hand held camera (tonemapped).
There are enlarged parts on the right side. The top half shows the middle
section of (1) the original mid-exposed photograph, (2) HDR image registered
for shift only, (3) HDR image registered for shift and rotation at a subpixel
resolution. The bottom half shows the lower left corner of (1) the original
mid-exposed photograph, (2) HDR image registered for shift only, (3) HDR
image registered for shift and rotation at a subpixel resolution.
39
registered at a subpixel resolution by modeling shift and rotation provide
a high dynamic range image which is only slightly softer than the source
images. The stitched high dynamic range image provides almost the same
amount of detail compared to the information contained in the all source
images. However, the total amount of detail in the high dynamic range
image is much higher when comparing to one source image only. This is
caused by the missing information in the overexposed or underexposed areas
of the source image.
6.6 Sample Scene
We prepared Figure 11 to demonstrate the composition of a high dynamic
range image from a set of dierently exposed images of the same scene. It
also shows the unique features of the high dynamic range image.
Figure 11a represents eleven dierently exposed images of the high con-
trast scene. The individual images were captured in RAW by a Nikon D70
camera mounted on the tripod. The scene and illumination had not changed
between the individual images. The exposure dierence between all images is
one stop. The usage of linear RAW data enabled us to skip the complicated
camera response recovery and provided us the 12-bit precision of captured
brightness values.
The image in Figure 11b is the composed high dynamic range image
linearly mapped to the dynamic range of the display device. The most pixels
in the scaled high dynamic range map are near the zero brightness level. The
only bright pixels are those of the wire in the bulb reected by the mirror.
Exploring the pixel values of the high dynamic range image, one can nd
that the brightest values are 3.3292 10
5
times greater than the darkest ones.
The image in Figure 11c is the high dynamic range image tonemapped
using iCAM algorithm [8]. The tonemapping algorithm reduces a global
contrast of the high dynamic range image while preserving a local contrast
of the image. The shadow details can be observed together with highlight
details in the tonemapped image.
The image in Figure 11d is the high dynamic range image displayed by
mapping the lower 0.1% of its dynamic range to the dynamic range of the dis-
play device. The image demonstrates nonuniform utilization of the dynamic
range. The most data are contained in the low 0.1% of the dynamic range.
The biggest part 99.9% of the dynamic range is used for the brightest
areas in image near the lamp and the mirror reecting the lamp.
The image in Figure 11e is a reconstruction of the second source image
used to create the high dynamic range image. It demonstrates that the high
40
(a)
(b) (c)
(d) (e)
Figure 11: (a) Source images. (b) HDR image linearly mapped to device. (c)
HDR image tonemapped using iCAM. (d) Low 0.1% of HDR image dynamic
range linearly mapped to device. (e) Reconstruction of the second source
image from HDR image.
41
dynamic range image contains information from the all images so that the
individual images can be reconstructed from it.
7 Future Work
At the end of writing this text we found two articles presenting more so-
phisticated methods for camera response recovery. The methods promise to
provide more accurate results than the calibration methods implemented in
this work. The articles are works by Tsin et al. [23] and Robertson et al.
[15]. It might be interesting to implement these methods and compare them
to the methods implemented in this work.
It might be also useful to implement a module for automatic image cap-
turing. In combination with computer driven camera the whole sequence of
dierently exposed images might be captured automatically. This task might
be accomplished using SDKs provided by some camera manufactures.
A ghosts appear in the image if some parts of the scene changes be-
tween exposures. If the ghosted areas could be detected the ghosts might
be removed by using only data from the best exposed image instead of the
weighted sum (8) which combines data from all images. The basic idea of
automatic ghost detection is to mark pixel as ghosted when the pixels in its
neighborhood has a large disagreement of irradiance values computed from
individual exposures using (7).
8 Conclusions
We showed that the recovery of camera response function using the dier-
ently exposed images is a dicult task. We implemented state-of-the-art
algorithms according to descriptions given by their authors. We learned that
they are very sensitive to the noise in data. In agreement with our results,
we see the way around the problem in the use of linear 12-bit RAW images
which are now oered by many modern mid-class and professional cameras.
We proposed a robust method recovering the exposure times from photo-
graphs captured by a linear response camera. The novel method improves
quality of stitched high dynamic range images by providing more accurate
exposure values than the values provided by the camera. While the method
estimates the exposure ratios from the image histograms it could be used with
unregistered photographs captured with a hand held camera. The knowledge
of proper exposure times may help later with image registration. The method
42
is suitable to be used with linear RAW images. It was implemented in Matlab
and tested on both synthetic and real data.
Next, we showed that capturing HDR images is possible without mount-
ing the camera on a tripod and described a method to register captured
images diering in exposure. The presented method expands the usability
of HDR widening to hand held photography. Method was implemented in
Matlab and tested on both synthetic and real data.
Finally, we implemented an HDR Library allowing anyone to exploit the
advantages of high dynamic range imaging. The HDR Library covers all
parts of high dynamic range capturing process. The HDR Library was de-
signed to be easily extensible by anyone wanting to improve it in the future.
The individual parts of the HDR Library were implemented as Matlab func-
tions. The HDR Library was contributed to public domain and is available
to anyone interested in high dynamic range imaging.
43
A DVD Contents
The attached DVD contains, among others, the implemented HDR Library,
Matlab demo, sample high dynamic range images and the electronic version of
the diploma thesis. The CVWW conference paper Exposure time estimation
for high dynamic range imaging with hand held camera based on results of
this work is also included. The detailed description of DVD contents follows.
All content can be viewed a executed directly from the DVD with the
exception of the image registration demo. The registration demo is noninter-
active and writes stitched high dynamic range images to directory out. To
run registration demo, copy DVD to your disk and ensure that the directory
out is writable.
On Linux/Unix the libti/raw2ti must be properly installed to enable
writing of high dynamic range images. For windows the libti is included on
the DVD and does not require installation.
To view high dynamic range images in 32-bit oating point TIFF use
included HDRView.
3rd-party Third party software.
dcraw Tool for converting RAW images to linear TIFF.
hdrview Viewer for high dynamic range images.
icam-hdr-matlab Matlab implementation of iCAM tonemapping.
libti Support for the Tag Image File Format (TIFF), only for windows
platform.
code HDR Library source code.
cvww CVWW conference paper based on this work.
data Dierently exposed photographs used to create sample high dynamic
range images.
demo Demos illustrating use of HDR Library.
calibration.m Demonstrates function of implemented response recovery
methods. The response is recovered from six synthetic images.
registration full.m Demonstrates stitching of three dierently exposed
hand held shots at full resolution. Images are registered for shift and
rotation at a subpixel resolution. Warning: This demo takes a lot of
time to proceed.
registration small.m Demonstrates of stitching of three dierently ex-
posed hand held shots. The smaller versions of original photographs
are used to speedup this demo. For the same reason, image shift only
is determined at a pixel resolution using correlation in the Fourier
domain.
44
samplescene.m Illustrates stitching of eleven dierently exposed pho-
tographs captured using a camera mounted on the tripod.
out Output directory for image registration demo.
report Electronic version of the diploma thesis.
samples Sample high dynamic range images.
45
B HDR Library Documentation
hdrBTFFromHist
Recovers the brightness transfer function of two dierently exposed images
from their histograms.
Synopsis
btf = hdrBTFFromHist( h1, h2 )
Arguments
btf Samples of recovered brightness transfer function. Number of samples
is equal to the number of bins in supplied the histograms.
h1 Histogram of the rst image.
h2 Histogram of the second image, must have the same number of bins
as h1.
hdrExposures
Estimates exposure times of images captured by the linear camera.
Synopsis
exposures = hdrExposures( imageSet, saturationLevel )
Arguments
exposures Vector of estimated exposure times. Time values are normal-
ized so that the darkest image has unit exposure time.
imageSet Set of dierently exposed images, see hdrReadSet for details.
saturationLevel Ignore brightness values above saturationLevel, acceptable
values are from interval 0, 1, default is 0.98.
hdrInverseResponse
Estimates inversion of camera response function from a set of dierently
exposed photographs.
Synopsis
g = hdrInverseResponse( imageSet, model, saturationLevel, args )
Arguments
g.response(c).model Identify model used to recover the response of chan-
nel c.
g.response(c).func Calling feval(g.response(c).func, g.response(c), x) eval-
uates response of channel c at each element of x.
46
g.response(c).cert Calling feval(g.response(c).cert, g.response(c), x) evalu-
ates response certainty of channel c at each element of x.
g.wb Color balance data.
imageSet Set of dierently exposed images, see hdrReadSet for details.
model Identify which model should be used for response recovery. Can
be one of Debevec, DebevecHistogram, Nayar, NayarHistogram,
SplitLine and Linear. Default is Debevec.
saturationLevel Ignore the brightness values above saturationLevel, accept-
able values are from interval 0, 1, default is 0.98.
args Optional arguments for some models. It can be order of the polyno-
mial for Nayar and NayarHistogram models or it can specify desired
precision for SplitLine model.
hdrReadSet
Reads set of dierently exposed images of the same scene.
Synopsis
imageSet = hdrReadSet( imagenames, exposures )
Arguments
imageSet.images Matrix (i,j,c,k) of k dierently exposed images of size
ij each with c channels.
imageSet.exposures Vector of exposure times.
imagenames Cell array of image names. All images must be of the same
size.
exposures Vector of exposure times.
hdrRegister
Register images in the set.
Synopsis
H = hdrRegister( imageSet, g, subpixel, noiseLevel, normTresh )
Arguments
H(:,:,j) 33 matrix dening transformation of j-th image.
imageSet Set of dierently exposed images, see hdrReadSet for details.
g.model Identify model used for response recovery.
g.func Calling feval(g.func, g, x) evaluates response at each element of x.
g.cert Calling feval(g.cert, g, x) evaluate response certainty at each ele-
ment of x.
47
noiseLevel Discard pixel values below noiseLevel, acceptable values are
from interval 0, 1, default is 0.008.
normTresh Normalization threshold, aect optimizing criterion, accept-
able values are from interval 0, 1, default is 0.04.
hdrSpy
Displays tonemapped portion of high dynamic range image. Portion is se-
lected interactively by user.
Synopsis
im = hdrSpy( hdri )
Arguments
im Tonemapped portion of high dynamic range image.
hdri High dynamic range image.
hdrStitch
Stitches dierently exposed images into one high dynamic range image.
Synopsis
hdri = hdrStitch( imageSet, inverseResponse )
Arguments
hdri Stitched high dynamic range image.
imageSet Set of dierently exposed images, see hdrReadSet for details.
inverseResponse Inverse of camera response function, see hdrInverseRe-
sponse for details.
H(:,:,j) 3x3 matrix dening transformation of j-th image.
hdrToneMap
Performs tonemapping of linear high dynamic range image.
Synopsis
mapped = hdrToneMap( hdri )
Arguments
hdri High dynamic range image.
scale Aects brightness/saturation of tonemapped image, default is 50.
48
hdrToRadiance
Convert images in the set to radiance maps using the inverse camera response.
Synopsis
radianceMap = hdrToRadiance( imageSet, inverseResponse )
Arguments
radianceMap Computed radiance maps. The dimension are same as the
dimensions of imageSet.images.
imageSet Set of dierently exposed images, see hdrReadSet for details.
inverseResponse Inverse of the camera response function, see hdrInverse-
Response for details.
hdrWrite
Saves high dynamic range image to disk in a 32-bit oating point TIFF. On
Linux/Unix the libti/raw2ti must be properly installed to enable writing
of high dynamic range images. For windows the libti is included on the
DVD and does not require installation.
Synopsis
hdrWrite( hdri, le )
Arguments
hdri High dynamic range image.
le Filename of le to save high dynamic range image to.
49
References
[1] A. Adams. The Negative. Bulnch, 1995.
[2] Young-Chang Chang and J. F. Reid. RGB calibration for color image
analysis in machine vision. IEEE Transactions on image processing,
5(10):14141422, 1996.
[3] P. E. Debevec and J. Malik. Recovering high dynamic range radiance
maps from photographs. In SIGGRAPH, pages 369378, 1997.
[4] F. Durand and J. Dorsey. Fast bilateral ltering for the display of high-
dynamic-range images. In SIGGRAPH, pages 257266, 2002.
[5] R. Fattal, D. Lischinski, and M. Werman. Gradient domain high dy-
namic range compression. In SIGGRAPH, pages 249256, 2002.
[6] M. D. Grossberg and Nayar S. K. High dynamic range from multiple
images: Which exposures to combine? In CPMCV, 2004.
[7] M. D. Grossberg and S. K. Nayar. What can be known about the
radiometric response function from images? In European Conference on
Computer Vision, Part IV, pages 189205, 2002.
[8] G. M. Johnson and M. D. Fairchild. Rendering HDR images. In Color
Imaging Conference, pages 3641, 2003.
[9] S. B. Kang, M. Uyttendaele, S. Winder, and R. Szeliski. High dynamic
range video. ACM Trans. Graph., 22(3):319325, 2003.
[10] S. Mann, C. Manders, and J. Fung. Painting with looks: Photographic
images from video using quantimetric processing. In ACM international
conference on Multimedia, pages 117126, 2002.
[11] S. Mann and R. Picard. Being undigital with digital cameras: Extend-
ing dynamic range by combining dierently exposed pictures. In IS&T
46th annual conference, pages 422428, 1995.
[12] T. Mitsunaga and S. K. Nayar. Radiometric self calibration. In Com-
puter Vision and Pattern Recognition, volume 2, pages 374380, 1997.
[13] S. K. Nayar and T. Mitsunaga. High dynamic range imaging: Spatially
varying pixel exposures. In Computer Vision and Pattern Recognition,
volume 1, pages 472479, 2000.
50
[14] E. Reinhard, M. Stark, P. Shirley, and J. Ferwerda. Photographic tone
reproduction for digital images. In SIGGRAPH, pages 267276, 2002.
[15] M. A. Robertson, S. Borman, and R. L. Stevenson. Estimation-theoretic
approach to dynamic range enhancement using multiple exposures. Elec-
tronic Imaging, 12(2):219228, 2003.
[16] Y. Y. Schechner and S. K. Nayar. Generalized mosaicing: High dynamic
range in a wide eld of view. Computer Vision, 53(3), 2003.
[17] IMS Chips HDRC. http://www.ims-chips.com/home.php3?id=d0841.
[18] Kodak KAC-9628. http://www.kodak.com /global/en/digital/ccd/
products / cmos / KAC- 9628 / indexKAC- 9628 . jhtml ? id=0.1.10.4.11 /
&lc=en.
[19] Pixim. http://www.pixim.com/html/tech about.htm.
[20] Silicon Vision. http://www.siliconvision.de.
[21] SMaL Camera. http://www.smalcamera.com/technology.html.
[22] J. Stumpfel, C. Tchou, A. Jones, T. Hawkins, A. Wenger, and P. E.
Debevec. Direct HDR capture of the sun and sky. In AFRIGRAPH,
pages 145149, 2004.
[23] Y. Tsin, V. Ramesh, and T. Kanade. Statistical calibration of CCD
imaging process. In International Conference on Computer Vision, vol-
ume 1, July 2001.
[24] Jack Tumblin and Greg Turk. LCIS: A boundary hierarchy for detail-
preserving contrast reduction. In SIGGRAPH, pages 8390, 1999.
[25] M.

Cadk. Tone mapping operators .
http://www.cgg.cvut.cz/cadikm/tmo.
[26] G. Ward. Fast, robust image registration for compositing high dynamic
range photographcs from hand-held exposures. Graphics Tools, 8(2):17
30, 2003.
[27] Wikipedia. Vignetting. http://en.wikipedia.org/wiki/Vignetting, 2005.
51