High dynamic range imaging

(1)

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 _i

t

t i

i  ^



 

0





 L i  d E_i ( , _i)

i



) , (p  L

 p

Assume a static scene, Thus, L is not a function of time.

Camera pipeline

12 bits 8 bits

(2)

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 10⁶ cd/m²

• Often 1 : 100,000 in a scene

• Typical 1:50, max 1:500 for pictures

10

^-6

10

⁶

Real world

High dynamic range

spotmeter

(3)

Short exposure

10

^-6

10

⁶

10

^-6

10

⁶

Real world radiance

Picture intensity

dynamic range

Pixel value 0 to 255

Long exposure

10

^-6

10

⁶

10

^-6

10

⁶

Real 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

(4)

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.

(5)

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

(6)

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

(7)

t = 1/4 sec

t = 1 sec

t = 1/8 sec

t = 2 sec

Image series

t = 1/2 sec

Recovering response curve

• 1

• 3

•

• 3 3

• 3

• 2

•

• 2 2

• 2

0 255

t = 1/4 sec

t = 1 sec

t = 1/8 sec

t = 2 sec

Image series

t = 1/2 sec

Recovering response curve

• 1

• 3

•

• 3 3

• 3

• 2

•

• 2 2

• 2

X_ij= ln X_ij

Idea behind the math

ln2

Idea behind the math

Each line for a scene point.

The offset is essentially determined by the unknown E_i

(8)

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

(9)

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 to zero 2.



































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)

lnE₁

lnE_n :

::

(10)

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ˆ  ^ ^T

pseudo inverse

















0 0

/ 1

0 0

/ 1 ₁











r

 Σ

Proof

(11)

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));

(12)

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

(13)

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

(14)

Lens flare removal

before after

Easier HDR reconstruction

raw image = 12-bit CCD snapshot

Easier HDR reconstruction

X^ij=E_i* Δt_j

Exposure (X)

Δ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)

(15)

(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

mZ

c Z

g X

0

) (

Mitsunaga and Nayar

• To find the coefficients c_m to minimize the following

  

    



 

 



 



^N

i 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

(16)

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 m k

j j

Z c

Z c R N

1 0

1 , ) ( 0

) ( )

( 1 ,

1 Space of response curves Space of response curves

(17)

Robertson et. al.

) ( _i _j

ij f E t

Z  

j i ij

ij f Z E t

Z

g( ) ^¹( ) 

Given and , the goal is to find both and ^ij

Z t_j Ei g(Z_ij)

 

_







  







ij

j i ij ij

j ij

i g Z t w Z g Z E t

E ( ) ( ) ²

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

) (Z_ij

g E_i

Ei g(Z_ij)

 



^ ^



ij

j i ij E ij

i g w Z g Z E t

E g

i

2

, ( ) ( )

min ˆ arg

ˆ,

Robertson et. al.

repeat

) (Z_ij

g E_i

Ei g(Z_ij)

 



^ ^



ij

j i ij E ij

i g w Z g Z E t

E g

i

2

, ( ) ( )

min ˆ arg

ˆ,

Robertson et. al.

repeat

) (Z_ij

g E_i

Ei g(Z_ij)

 



^ ^



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 ₂

) (

) ( ) (

(18)

Robertson et. al.

repeat

) (Z_ij

g E_i

Ei g(Z_ij)

 



^ ^



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

Patch-Based HDR

Deep learning HDR assembly Deep learning HDR assembly

(19)

Deep single-image HDR reconstruction

DrTMO

HDRCNN

ExpandNet

Learning to reverse the pipeline

Camera pipeline

Learning to reverse the camera pipeline

Comparison

input DrTMO HDRCNN

ExpandNet Ours

Truth

Input

76

(20)

Result

77

Input

78

Result

79

HDR Video

• High Dynamic Range Video

Sing Bing Kang, Matthew Uyttendaele, Simon Winder, Richard Szeliski

SIGGRAPH 2003

video

(21)

Assorted pixel Assorted pixel

Assorted pixel A Versatile HDR Video System

video

(22)

A Versatile HDR Video System HDR becomes common practice

• Many cameras has bracket exposure modes

• For example, since iPhone 4, iPhone has HDR option. But, it could be more exposure blending rather than true HDR.

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.