High dynamic range imaging
Digital Visual Effectsg 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 the radiance distribution of a scene because it cannot capture the full spectral content and dynamic range
dynamic range.
• Limitations in sensor design prevent cameras
f t i ll i f ti d b l
from capturing all information passed by lens.
Camera pipeline
) , (p L
L i d Ei ( , i)
p
t E dt E X
t
i
i
Assume a static
scene, Thus, L is not a function
lens shutter sensor
t E dt E
X i
t i
i
0 is not a function
of time.
lens shutter
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 by camera makers who consider it part of their proprietary product differentiation. In addition, they are beyond the standard gamma curves they are beyond the standard gamma curves.
The world is high dynamic range
1
1 500 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 Of 1 100 000 i
• Often 1 : 100,000 in a scene
• Typical 1:50, max 1:500 for pictures
10
-610
6R l ld Real world
High dynamic range High dynamic range
spotmeter
Short exposure
10
-610
6Real world dynamic range
10
-610
6radiance
10
610
6Picture
intensity Pixel value 0 to 255 intensity Pixel value 0 to 255
Long exposure
10
-610
6Real world dynamic range
10
-610
6radiance
10
610
6Picture
intensity Pixel value 0 to 255 intensity Pixel value 0 to 255
Camera is not a photometer
• Limited dynamic range
P h lti l ?
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
Sh d
– Shutter speed – Aperture
– Neutral density filters
Shutter speed
• Note: shutter times usually obey a power series each “stop” is a factor of 2
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:
Usually really is:
¼ 1/8 1/16 1/32 1/64 1/128 1/256 1/512
¼, 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
I li ( if i d i i
• 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 method for this task, called Median Threshold Bitmap (MTB) alignment technique.
C id l i t l t l ti It i h
• 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 y g y g the input image using the median of intensities.
Why is MTB better than gradient?
• Edge-detection filters are dependent on image exposures
exposures
• Taking the difference of two edge bitmaps ld t i d i di ti f h th would not give a good indication of where the edges are misaligned.
Search for the optimal offset
• Try all possible offsets
offsets.
• Gradient descent
• Multiscale technique
• log(max_offset) levels
• Try 9 possibilities for
• Try 9 possibilities for the top level
S l b 2 h
• Scale by 2 when
passing down; try its 9 i hb
neighbors
Threshold noise
ignore pixels that are close to the threshold close to the threshold
exclusion bitmap
Efficiency considerations
• XOR for taking difference AND i h l i
• 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 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
curve 255
0
Recovering response curve Image series Image series
•
•
•
• •
•
•
• •
•
•
• • • • • • • • •
•
• 22
•
• 22
•
• 22
•
•
• 22
• 22
•
• 22
•
• 22
•
• 22
•
• 22
•
• 22
t =t =
t =t = t =t =
t =t = t =t =
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 3333 11
3333 3333 3333 3333
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
t =t = 1/2 sec 1/2 sec
255
0
Recovering response curve Image series Image series
•
•
•
• •
•
•
• •
•
•
• • • • • • • • •
•
• 22
•
• 22
•
• 22
•
•
• 22
• 22
•
• 22
•
• 22
•
• 22
•
• 22
•
• 22
t =t =
t =t = t =t =
t =t = t =t =
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 11
•
• 3333 11
3333 3333 3333 3333
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
t =t = 1/2 sec 1/2 sec
X Xij = ln Xij
Idea behind the math
ln2
Idea behind the math
Each line for a scene point.p The offset is essentially determined by the unknown Ei
unknown Ei
Idea behind the math
Note that there is a shift that we can’t recover
Basic idea
• Design an objective function O i i i
• Optimize it
Math for recovering response curve Recovering response curve
• The solution can be only up to a scale, add a constraint
constraint
• Add a hat weighting function
Recovering response curve
• We want
If P 11 N 25 ( i ll 50 i d) 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 SVDq g
How to optimize?
1. Set partial derivatives zero 22.
21 2
1
b b a
a
b ) l t l ti f
(
i 2
N
N 2
N 2 i
i
b x : a b :
x
a ) least-squaresolution of (
min
1
2
i
aN bN
Sparse linear system
256
n
n×p
g(0) :
g(255) lnE1
:
1 lnEn
: :
254
Ax=b
Questions
• Will g(127)=0 always be satisfied? Why or why not?
not?
• How to find the least-square solution for an
d t i d t ?
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 .Axb Axb If there are multiple solutions, we prefer the one with the minimal length .
b Ax x
Least-square solution for a linear system
If we perform SVD on A and rewrite it as
UΣ
TA V
then is the least-square solution.xˆ VΣUTb
pseudo inverse
1/ 1 0 0
1/
r
Σ
0
0 0 0
Proof Proof
Libraries for SVD
• Matlab GSL
• GSL
• Boost
• LAPACK
• ATLAS
• 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));
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 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;
k k+1;
end end
A(k,129) = 1; %% Fix the curve by setting its middle value to 0( , ) ; y g 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 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 raw image = 12-bit CCD snapshot
Easier HDR reconstruction
Exposure (X)
Xij=Ei* ΔtΔtjj
Δt Δt Δt Δt
Portable floatMap (.pfm)
• 12 bytes per pixel, 4 for each channel
sign exponent mantissa
T t h d i il t J ff P k ’
PF 768 512
Text header similar to Jeff Poskanzer’s .ppm image format:
768 512 1
<binary image data>
Floating Point TIFF similar Floating Point TIFF similar
Radiance format (.pic, .hdr, .rad)
32 bits/pixel
Red Green Blue Exponent
(145, 215, 87, 149) = (145 215 87) * 2^(149 128)
(145, 215, 87, 103) = (145 215 87) * 2^(103 128) (145, 215, 87) * 2^(149-128) =
(1190000, 1760000, 713000)
(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 Ward, Greg. Real Pixels, in Graphics Gems IV, edited by James Arvo, Academic Press, 1994
ILM’s OpenEXR (.exr)
• 6 bytes per pixel, 2 for each channel, compressed
sign exponent mantissa
Several lossless compression options 2:1 typical
• Several lossless compression options, 2:1 typical
• Compatible with the “half” datatype in NVidia's Cg
• Supported natively on GeForce FX and Quadro FXSupported natively on GeForce FX and Quadro FX
• Available at http://www.openexr.net/
Radiometric self calibration
• Assume that any response function response function can be modeled as a high order as a high-order polynomial
M XX
M
m
m
mZ
c Z
g X
0
) (
• No need to know exposure time in p advance. Useful
for cheap p Z
cameras
Mitsunaga and Nayar
• To find the coefficients cm to minimize the following
following
N P M M 2
i j m
m j i m j
j m
m ij
m
Z R c Z
c
1 1 0
1 , 1
, 0
j
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
1 1 1
,j i j j
i
Mitsunaga and Nayar
• Again, we can only solve up to a scale. Thus, add a constraint f(1) 1 It reduces to M add a constraint f(1)=1. It reduces to M variables.
H t l it?
• How to solve it?
Mitsunaga and Nayar
• We solve the above iteratively and update the exposure ratio accordingly
exposure ratio accordingly
M (k )
N mMm k ij k m k
j j
Z c
R
0) ( )
( 1
1
i M
m j i k m j
j
Z N c
R
1 0
1 , ) ( 1
,
• How to determine M? Solve up to M=10 and pick th ith th i i l N ti th t
m 0
up the one with the minimal error. Notice that you prefer to have the same order for all
h l U th bi d
channels. Use the combined error.
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 ij
Z tj Eii g(Zijij)
Maximum likelihood
ij
j i ij ij
j ij
i g Z t w Z g Z E t
E ( ) ( ) 2
2 exp 1 )
,
| ,
Pr( ij
ij ij i j
E
i g w Z g Z E t
E
gˆ, ˆ argmin
( ) ( ) 2E ij g, i
Robertson et. al.
ij ij i j
E
i g w Z g Z E t
E g
i
) 2
( ) ( min ˆ arg
ˆ, repeat
assuming is known optimize for g(Z ) E
E ij g, i
assuming is known, optimize for assuming is known, optimize for
) (Zij
g Ei
Ei g(Zij)
until converge
Robertson et. al.
ij ij i j
E
i g w Z g Z E t
E g
i
) 2
( ) ( min ˆ arg
ˆ, repeat
assuming is known optimize for g(Z ) E
E ij g, i
assuming is known, optimize for assuming is known, optimize for
) (Zij
g Ei
Ei g(Zij)
until converge
j
j ij
ij g Z t
Z w E
) ( ) (
j
j ij
i w Z t
E 2
) (
Robertson et. al.
ij ij i j
E
i g w Z g Z E t
E g
i
) 2
( ) ( min ˆ arg
ˆ, repeat
assuming is known optimize for g(Z ) E
E ij g, i
assuming is known, optimize for assuming is known, optimize for
) (Zij
g Ei
Ei g(Zij)
until converge
Em
ij
j i m
t E E
m
g | |
) 1 (
j m
1 ) 128
(
gnormalize so that
1 ) 128
(
g
Space of response curves
Space of response curves HDR Video
• High Dynamic Range Video
Sing Bing Kang Matthew Uyttendaele Simon Sing Bing Kang, Matthew Uyttendaele, Simon Winder, Richard Szeliski
SIGGRAPH 2003 SIGGRAPH 2003
video
Assorted pixel Assorted pixel
Assorted pixel Assignment #1 HDR image assemble
• Work in teams of two T ki i
• Taking pictures
• Assemble HDR images and optionally the response curve.
• Develop your HDR using tone mappingp y g pp g
Taking pictures
• Use a tripod to take multiple photos with different shutter speeds Try to fix anything different shutter speeds. Try to fix anything else. Smaller images are probably good enough.
Th t t f t t i il bl
• 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 – PSRemote
– AHDRIA (obsolete since 2007)
• PSRemote
• PSRemote
– Manual (can use auto-bracket mode) – Not freeNot free
– Supports both jpg and raw
– Support most Canon’s PowerShot camerasSupport most Canon s PowerShot cameras
• AHDRIA
– AutomaticAutomatic – Free
– Only supports jpgy pp jpg – Support less models
2. HDR assembling
• Write a program to convert the captured images into a radiance map and optionally to images into a radiance map and optionally to output the response curve.
W ill id i I/O lib hi h
• We will provide image I/O library which
supports many traditional image formats such
j d d fl t i t i h
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
your photograph.
– Reinhard’s algorithm (HDRShop plugin) Ph t ti
– Photomatix – LogView
F Bil l ( Li l )
– 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
I l i l i h
• 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 executable, a html report, pictures you have taken, HDR image, and an artifact (tone- mapped image)
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
• 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 suggested that you implement at least one bonus (MTB/tone mapping/other HDR construction)
construction)
• You have a total of 10 days of delay without lt f i t ft th t 1 i t penalty for assignments; after that, -1 point per day applies in your final grade until
hi f h j t
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
• 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
• 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.
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.
