High dynamic range imaging
Digital Visual Effects Yung-Yu Chuang
with slides by Fredo Durand, Brian Curless, Steve Seitz, Paul Debevec and Alexei Efros
Camera is an imperfect device
• Camera is an imperfect device for measuring the radiance distribution of a scene because it cannot capture the full spectral content and dynamic range.
• Limitations in sensor design prevent cameras from capturing all information passed by lens.
Camera pipeline
lens shutter sensor
t E dt
E
X
it
t
i
i
0
L i
d
Ei ( , i )i
i
) ,
( p
L p
Assume a static scene, Thus, L is not a function of time.
Camera pipeline
12 bits 8 bits
Real-world response functions
In general, the response function is not provided by camera makers who consider it part of their proprietary product differentiation. In addition, they are beyond the standard gamma curves.
The world is high dynamic range
1
1,500
25,000
400,000
2,000,000,000
The world is high dynamic range
Real world dynamic range
• Eye can adapt from ~ 10-6 to 106 cd/m2
• Often 1 : 100,000 in a scene
• Typical 1:50, max 1:500 for pictures
10
-610
6Real world
High dynamic range
spotmeter
Short exposure
10
-610
610
-610
6Real world radiance
Picture intensity
dynamic range
Pixel value 0 to 255
Long exposure
10
-610
610
-610
6Real world radiance
Picture intensity
dynamic range
Pixel value 0 to 255
Camera is not a photometer
• Limited dynamic range
Perhaps use multiple exposures?
• Unknown, nonlinear response
Not possible to convert pixel values to radiance
• Solution:
– Recover response curve from multiple exposures, then reconstruct the radiance map
Varying exposure
• Ways to change exposure – Shutter speed
– Aperture
– Neutral density filters
Shutter speed
• Note: shutter times usually obey a power series – each “stop” is a factor of 2
• ¼, 1/8, 1/15, 1/30, 1/60, 1/125, 1/250, 1/500, 1/1000 sec
Usually really is:
¼, 1/8, 1/16, 1/32, 1/64, 1/128, 1/256, 1/512, 1/1024 sec
Varying shutter speeds
HDRI capturing from multiple exposures
• Capture images with multiple exposures
• Image alignment (even if you use tripod, it is suggested to run alignment)
• Response curve recovery
• Ghost/flare removal
Image alignment
• We will introduce a fast and easy-to-implement method for this task, called Median Threshold Bitmap (MTB) alignment technique.
• Consider only integral translations. It is enough empirically.
• The inputs are N grayscale images. (You can either use the green channel or convert into grayscale by Y=(54R+183G+19B)/256)
• MTB is a binary image formed by thresholding the input image using the median of intensities.
Why is MTB better than gradient?
• Edge-detection filters are dependent on image exposures
• Taking the difference of two edge bitmaps
would not give a good indication of where the edges are misaligned.
Search for the optimal offset
• Try all possible offsets.
• Gradient descent
• Multiscale technique
• log(max_offset) levels
• Try 9 possibilities for the top level
• Scale by 2 when
passing down; try its 9 neighbors
Threshold noise
ignore pixels that are close to the threshold
exclusion bitmap
Efficiency considerations
• XOR for taking difference
• AND with exclusion maps
• Bit counting by table lookup
Results
Success rate = 84%. 10% failure due to rotation.
3% for excessive motion and 3% for too much high-frequency content.
Recovering response curve
12 bits 8 bits
Recovering response curve
• We want to obtain the inverse of the response curve
0 255
t =t = 1/4 sec 1/4 sec
t =t = 1 sec 1 sec
t =t = 1/8 sec 1/8 sec
t =t = 2 sec 2 sec
Image series Image series
t =t = 1/2 sec 1/2 sec
Recovering response curve
• • 11
• •
11 • •
11
• •
11 • •
11
• •
11 • •
11
• •
11 • •
11
• • 11
• • 33
• • 33
• • 33
• • 33
• • 33
• • 33
• • 33
• •
33 • •
33
• • 33
• • 22
• • 22
• • 22
• • 22
• • 22
• • 22
• • 22
• •
22 • •
22
• • 22
0 255
t =t = 1/4 sec 1/4 sec
t =t = 1 sec 1 sec
t =t = 1/8 sec 1/8 sec
t =t = 2 sec 2 sec
Image series Image series
t =t = 1/2 sec 1/2 sec
Recovering response curve
• • 11
• •
11 • •
11
• •
11 • •
11
• •
11 • •
11
• •
11 • •
11
• • 11
• • 33
• • 33
• • 33
• • 33
• • 33
• • 33
• • 33
• •
33 • •
33
• • 33
• • 22
• • 22
• • 22
• • 22
• • 22
• • 22
• • 22
• •
22 • •
22
• • 22
Xij = ln Xij
Idea behind the math
ln2
Idea behind the math
Each line for a scene point.
The offset is essentially determined by the
unknown Ei
Idea behind the math
Note that there is a shift that we can’t recover
Basic idea
• Design an objective function
• Optimize it
Math for recovering response curve
Recovering response curve
• The solution can be only up to a scale, add a constraint
• Add a hat weighting function
Recovering response curve
• We want
If P=11, N~25 (typically 50 is used)
• We prefer that selected pixels are well
distributed and sampled from constant regions.
They picked points by hand.
• It is an overdetermined system of linear equations and can be solved using SVD
How to optimize?
1. Set partial derivatives zero 2.
N 2 1
N 2 1
i i
b : b b x
a : a a b
x
a ) least -squaresolution of (
min
1 N 2 i
Sparse linear system
Ax=b
256 n
n×p
1 254
g(0)
g(255) lnE1
lnEn :
::
Questions
• Will g(127)=0 always be satisfied? Why or why not?
• How to find the least-square solution for an over-determined system?
Least-square solution for a linear system
b Ax
n
m n m n
m
They are often mutually incompatible. We instead find x to minimize the norm of the residual vector . If there are multiple solutions, we prefer the one with the minimal length .
b Ax b
Ax x
Least-square solution for a linear system
If we perform SVD on A and rewrite it as
then is the least-square solution.
UΣ T
A V
b U VΣ
x ˆ
Tpseudo inverse
0 0
0 0
/ 1
0 0
/ 1 1
r
Σ
Proof
Proof
Libraries for SVD
• Matlab
• GSL
• Boost
• LAPACK
• ATLAS
Matlab code
Matlab code
function [g,lE]=gsolve(Z,B,l,w) n = 256;
A = zeros(size(Z,1)*size(Z,2)+n+1,n+size(Z,1));
b = zeros(size(A,1),1);
k = 1; %% Include the data-fitting equations for i=1:size(Z,1)
for j=1:size(Z,2) wij = w(Z(i,j)+1);
A(k,Z(i,j)+1) = wij; A(k,n+i) = -wij; b(k,1) = wij * B(i,j);
k=k+1;
end end
A(k,129) = 1; %% Fix the curve by setting its middle value to 0 k=k+1;
for i=1:n-2 %% Include the smoothness equations
A(k,i)=l*w(i+1); A(k,i+1)=-2*l*w(i+1); A(k,i+2)=l*w(i+1);
k=k+1;
end
x = A\b; %% Solve the system using SVD g = x(1:n);
lE = x(n+1:size(x,1));
Recovered response function
Constructing HDR radiance map
combine pixels to reduce noise and obtain a more reliable estimation
Reconstructed radiance map
What is this for?
• Human perception
• Vision/graphics applications
Automatic ghost removal
before after
Weighted variance
Moving objects and high-contrast edges render high variance.
Region masking
Thresholding; dilation; identify regions;
Best exposure in each region
Lens flare removal
before after
Easier HDR reconstruction
raw image =
12-bit CCD snapshot
Easier HDR reconstruction
Xij=Ei* ΔtΔtjj
Exposure (X)
ΔtΔt
• 12 bytes per pixel, 4 for each channel
sign exponent mantissa
PF
768 512 1
<binary image data>
Floating Point TIFF similar
Text header similar to Jeff Poskanzer’s .ppm image format:
Portable floatMap (.pfm)
(145, 215, 87, 149) =
(145, 215, 87) * 2^(149-128) = (1190000, 1760000, 713000)
(145, 215, 87, 103) =
(145, 215, 87) * 2^(103-128) =
(0.00000432, 0.00000641, 0.00000259)
Ward, Greg. "Real Pixels," in Graphics Gems IV, edited by James Arvo, Academic Press, 1994
Radiance format (.pic, .hdr, .rad)
Red Green Blue Exponent
32 bits/pixel
ILM’s OpenEXR (.exr)
• 6 bytes per pixel, 2 for each channel, compressed
sign exponent mantissa
• Several lossless compression options, 2:1 typical
• Compatible with the “half” datatype in NVidia's Cg
• Supported natively on GeForce FX and Quadro FX
• Available at http://www.openexr.net/
Radiometric self calibration
• Assume that any response function can be modeled as a high-order polynomial
• No need to know exposure time in advance. Useful for cheap
cameras Z
X
M
m
m mZ c Z
g X
0
) (
Mitsunaga and Nayar
• To find the coefficients cm to minimize the following
Ni
P j
M m
m j i m j
j M
m
m ij
m
Z R c Z
c
1 1
2
0
1 , 1
, 0
A guess for the ratio of
1 1
1
,
j j j
i
j i j
i ij
t t t
E t E X
X
Mitsunaga and Nayar
• Again, we can only solve up to a scale. Thus, add a constraint f(1)=1. It reduces to M
variables.
• How to solve it?
Mitsunaga and Nayar
• We solve the above iteratively and update the exposure ratio accordingly
• How to determine M? Solve up to M=10 and pick up the one with the minimal error. Notice that you prefer to have the same order for all
channels. Use the combined error.
N i M
m
m j i k
m M
m
m k ij
k k m
j j
Z c
Z c
R N
1
0
1 , )
( 0
) ( )
(
1 ,
1
Robertson et. al.
) ( i j
ij f E t
Z
j i
ij
ij f Z E t
Z
g( ) 1( )
Given and , the goal is to find both and
Zij t j Ei g(Zij )
ij
j i
ij ij
j ij
i g Z t w Z g Z E t
E ( ) ( ) 2
2 exp 1
) ,
| , Pr(
ij
j i
ij E ij
i g w Z g Z E t
E g
i
2
, ( ) ( )
min ˆ arg
ˆ,
Maximum likelihood
Robertson et. al.
repeat
assuming is known, optimize for assuming is known, optimize for
until converge
) (Zij
g Ei
Ei g(Zij )
ij
j i
ij E ij
i g w Z g Z E t
E g
i
2
, ( ) ( )
min ˆ arg
ˆ,
Robertson et. al.
repeat
assuming is known, optimize for assuming is known, optimize for
until converge
) (Zij
g Ei
Ei g(Zij )
ij
j i
ij E ij
i g w Z g Z E t
E g
i
2
, ( ) ( )
min ˆ arg
ˆ,
j
j ij
j
j ij
ij
i w Z t
t Z
g Z
w
E 2
) (
) (
) (
Robertson et. al.
repeat
assuming is known, optimize for assuming is known, optimize for
until converge
) (Zij
g Ei
Ei g(Zij )
ij
j i
ij E ij
i g w Z g Z E t
E g
i
2
, ( ) ( )
min ˆ arg
ˆ,
Em
ij
j i
m
t E E
m
g | |
) 1 (
1 )
128
(
gnormalize so that
Space of response curves
Space of response curves
HDR Video
• High Dynamic Range Video
Sing Bing Kang, Matthew Uyttendaele, Simon Winder, Richard Szeliski
SIGGRAPH 2003
video
Assorted pixel
Assorted pixel
Assorted pixel
Assignment #1 HDR image assemble
• Work in teams of two
• Taking pictures
• Assemble HDR images and optionally the response curve.
• Develop your HDR using tone mapping
Taking pictures
• Use a tripod to take multiple photos with
different shutter speeds. Try to fix anything
else. Smaller images are probably good enough.
• There are two sets of test images available on the web.
• We have tripods and a Canon PowerShot G7 for you to borrow starting from the next week.
• Try not touching the camera during capturing.
But, how?
1. Taking pictures
• Use a laptop and a remote capturing program.
– PSRemote
– AHDRIA (obsolete since 2007)
• PSRemote
– Manual (can use auto-bracket mode) – Not free
– Supports both jpg and raw
– Support most Canon’s PowerShot cameras
• AHDRIA
– Automatic – Free
– Only supports jpg – Support less models
2. HDR assembling
• Write a program to convert the captured
images into a radiance map and optionally to output the response curve.
• We will provide image I/O library which
supports many traditional image formats such as .jpg and .png, and float-point images such as .hdr and .exr.
• Paul Debevec’s method. You will need SVD for this method.
• Recover from CCD snapshots. You will need dcraw.c.
3. Tone mapping
• Apply some tone mapping operation to develop your photograph.
– Reinhard’s algorithm (HDRShop plugin) – Photomatix
– LogView
– Fast Bilateral (.exr Linux only) – PFStmo (Linux only)
pfsin a.hdr | pfs_fattal02 | pfsout o.hdr
Bells and Whistles
• Other methods for HDR assembling algorithms
• Implement tone mapping algorithms
• Implement MTB alignment algorithm
• Others
Submission
• You have to turn in your complete source, the executable, a html report, pictures you have taken, HDR image, and an artifact (tone-
mapped image).
• Report page contains:
description of the project, what do you learn,
algorithm, implementation details, results, bells and whistles…
• The class will have vote on artifacts.
• Submission mechanism will be announced later.
Reference software
• Photomatix
• AHDRIA/AHDRIC
• HDRShop
• RASCAL
Assignment #1 HDR image assemble
• It is for warming up and considered easier; it is suggested that you implement at least one
bonus (MTB/tone mapping/other HDR construction)
• You have a total of 10 days of delay without
penalty for assignments; after that, -1 point per day applies in your final grade until reaching
zero for each project.
References
References
• Paul E. Debevec, Jitendra Malik, Recovering High Dynamic Range Radiance Maps from Photographs, SIGGRAPH 1997.
• Tomoo Mitsunaga, Shree Nayar, Radiometric Self Calibration, CVPR 1999.
• Mark Robertson, Sean Borman, Robert Stevenson, Estimation-Theoretic Approach to Dynamic Range Enhancement using Multiple Exposures, Journal of Electronic Imaging 2003.
• Michael Grossberg, Shree Nayar, Determining the Camera Response from Images: What Is Knowable, PAMI 2003.
• Michael Grossberg, Shree Nayar, Modeling the Space of Camera Response Functions, PAMI 2004.
• Srinivasa Narasimhan, Shree Nayar, Enhancing Resolution Along Multiple Imaging Dimensions Using Assorted Pixels, PAMI 2005.
• G. Krawczyk, M. Goesele, H.-P. Seidel, Photometric Calibration of High Dynamic Range Cameras, MPI Research Report 2005.
• G. Ward, Fast Robust Image Registration for Compositing High Dynamic Range Photographs from Hand-held Exposures, jgt 2003.
