Transcript 07_mc

Monte Carlo Techniques
Basic Concepts
CIS782
Advanced Computer Graphics
Raghu Machiraju
© Machiraju/Möller
Reading
• Chapter 13, 14, 15 of “Physically Based
Rendering” by Pharr&Humphreys
• Chapter 7 in “Principles of Digital Image
Synthesis,” by A. Glassner
© Machiraju/Möller
Las Vegas Vs Monte Carlo
• Two kinds:
– Las Vegas - always give right
answer (but use elements of
randomness on the way)
– Monte Carlo - give the right
answer on average
• Very difficult to evaluate
• No closed forms
© Machiraju/Möller
Monte Carlo Integration
• Pick a set of evaluation points
• Simplifies evaluation of difficult integrals
• Accuracy grows with O(N-0.5), i.e. in order
to do twice as good we need 4 times as
many samples
• Artifacts manyfest themselves as noise
• Research - minimize error while minimizing
the number of necessary rays
© Machiraju/Möller
Basic Concepts
• X, Y - random variables
– Continuous or discrete
– Apply function f to get Y from X: Y=f(X)
• Example - dice
–
–
–
–
Set of events Xi = {1, 2, 3, 4, 5, 6}
6
F - rolling of dice
pj 1

Probability pi = 1/6 j1
x as a continuous, uniformly distributed
i1
i
random variable:
 pj x   pj

j1
© Machiraju/Möller
j1
Basic Concepts
• Cumulative distribution function (CDF)
P(x) of a random variable X:
P x   PrX  x 
x
 psds

• Dice example
– P(2) = 1/3
– P(4) = 2/3
– P(6)=1
© Machiraju/Möller
Continuous Variable
• Example - lighting
– Probability of sampling illumination based on power
Fi
Fi:
p 
i
F
j
j
• Canonical uniform random variable x
– Takes on all values in [0,1) with equal probability
– Easy 
to create in software (pseudo-random number
generator)
– Can create general random distributions by starting
with x
© Machiraju/Möller
Probability Distribution
Function
• Relative probability of a random variable
taking on a particular value
• Derivative of CDF: px   dP x 
dx
• Non-negative
b
 to one P x  a,b   px dx
• Always integrate
1 xa  0,1
• Uniform random variable: px  
0 otherwise
©
Machiraju/Möller
P x   x
Cond. Probability,
Independence
• We know that the outcome is in A
• What is the probability that it is in B?
B AB
A
Pr(B|A) = Pr(AB)/Pr(A)
Event space
• Independence: knowing A does not help:
Pr(B|A) = Pr(B) ==> Pr(AB)=Pr(A)Pr(B)
© Machiraju/Möller
Expected Value
• Average value of the function f over some
distribution of values p(x) over its domain
D
E p  f x     D f x px dx
• Example - cos over [0,p] where p is uniform
px   1 p


1
p
p
cos x
0
p
sin p  sin 0  0
© Machiraju/Möller
+
-

E p cosx  
Variance
• Expected deviation of f from its expected
value
• Fundamental concept of quantifying the
error in Monte Carlo (MC) methods

V  f x    E  f x   
2

© Machiraju/Möller
2

Properties
E af x   aE  f x 
E  f X i   E  f X i 
i
i
V af x   a 2V  f x 
• Hence we can write:


V  f x  E f x   
2

• For independent random variables:
V
 f X   V f X 
i
i
i

© Machiraju/Möller
i
2
Uniform MC Estimator
• All there is to it, really :)
• Assume we want to compute the integral of
f(x) over [a,b]
• Assuming uniformly distributed random
variables Xi in [a,b] (i.e. p(x) = 1/(b-a))
N
b

a
• Our MC estimator FN: F 
 f X 
N

© Machiraju/Möller
N
i
i1
Simple Integration
N
1
 f xdx   f x x
i
0
i1
N
0
1
1
  f x i 
N i1
1 
Error  O 
N 
© Machiraju/Möller
Trapezoidal Rule
1

0
N 1
x
f x dx    f x i   f x i1
2
i 0
N
1
0

1
  w i f x i 
N i1
0.5 i  0,N
w i  
 1 0  i  N
1 
Error  O 
N 
© Machiraju/Möller
Uniform MC Estimator
• FN is equal to the correct integral:
b  a N

(random variable Y=f(X))
E FN   E 
 f X i 
 N
i1

ba N

E  f X i 

N i1
N
b
ba

f x  px dx


a
N i1
N
b
1
   a f x dx
N i1
© Machiraju/Möller

 f x dx
b
a
General MC Estimator
• Can relax condition for general PDF
• Important for efficient evaluation of integral
N
f X i 
1
(random variable
FN  
Y=f(X)/p(X))
N i1 pX i 
1 N f X 
i
• And hence: E FN   E  

N i1 pX i 
N
b f x 
1
px dx
   a
N i1
px 

 f x dx
b
a
© Machiraju/Möller
Confidence Interval
• We know we should expect the correct
result, but how likely are we going to see it?
• Strong law of large numbers (assuming that
Yi are independent and identically
distributed):
N


1
Prlim Yi  E Y  1
N  N i1

© Machiraju/Möller
Confidence Interval
• Rate of convergence: Chebychev Inequality
V F 
PrF  E F   k 2
k
• Setting

• We have
V F 
 2
k



V F  
PrF  E F  
 

 



© Machiraju/Möller
MC Estimator
• How good is it? What’s our error?
• Our error (root-mean square) is in the
variance, hence
N


1
V FN   V 
 N

1
N
2

i 1
f xi 

pxi  
 f  xi  
V



p
x
i 
i 1 
N

1
 V F 
N
© Machiraju/Möller
MC Estimator
• Hence our overall error:

1

Pr  FN  EFN  

N

V F  


 
• V[F] measures square of RMS error!
• This result is independent of our dimension
© Machiraju/Möller
Distribution of the Average
• Central limit theorem assuming normal
distribution

F 
lim Pr  FN  E F   t

N 
N

t
e

2p
1
 x2 2
dx

• This can be re-arranged as

2
x 2 2
PrFN  I  t FN 
e
dx

p t
• well known Bell curve
N=160
N=40
N=10
© Machiraju/Möller
E[f(x)]
Distribution of the Average
• This can be re-arranged as
PrFN  I  t FN 
2

e
x 2 2
dx
p t
• Hence for t=3 we can conclude
PrFN  I  3 FN  0.997
• I.e. pretty much all results are
within three standard
N=160
deviations
(probabilistic error bound
N=40
- 0.997 confidence)
© Machiraju/Möller
E[f(x)]
N=10
Choosing Samples
• How to sample random variables?
• Assume we can do uniform distribution
• How to do general distributions?
– Inversion
– Rejection
– Transformation
© Machiraju/Möller
Inversion Method
• Idea - we want all the events to be
distributed according to y-axis, not x-axis
1
1
PDF
0
CDF
0
x
x
• Uniform distribution is easy!
1
0
1
PDF
x
© Machiraju/Möller
0
CDF
x
Inversion Method
• Compute CDF (make sure it is normalized!)
P x  
1
x
PDF
1
 psds
CDF

• Compute the inverse
0
P-1(y)
1
x
0
1
CDF
x
P-1

0
x
0
• Obtain a uniformly distributed random number x
• Compute Xi = P-1(x)
© Machiraju/Möller
x
Example - Power Distribution
•
•
•
•
•
px  cx n
0  x 1
Used in BSDF’s
1
n
cx
Make sure it is normalized  dx  1 c  n  1
0 x
n
n 1
P
x

n

1
s
ds

x
Compute the CDF
   

0
1
Invert the CDF
P x   n 1 x
Now we can choosea uniform x distribution
to get a power distribution!


X  n 1 x
© Machiraju/Möller
Example - Exponential Distrib.
•
•
•
•
•
px  ceax 0  x  
E.g. Blinn’s Fresnel Term 
 ax
ce
dx  1 c  a
Make sure it is normalized 
0 x
as
ax
P
x

ae
ds

1
e
Compute the CDF
 

0
1
P x    1 ln 1 x
Invert the CDF
a
Now we can choosea uniform x distribution
to get an exponential
distribution!

1 ln 1 x    1 ln x
X 
a
a
• extend to any funcs by piecewise approx.
© Machiraju/Möller
Rejection Method
• Sometimes
– We cannot integrate p(x)
– We can’t invert a function P(x) (we don’t have
the function description)
• Need to find q(x), such that p(x) < cq(x)
• Dart throwing
– Choose a pair of random variables (X, x)
– test whether x < p(X)/cq(X)
© Machiraju/Möller
Rejection Method
•
•
•
•
Essentially we pick a point (x, xcq(x))
If point lies beneath p(x) then we are ok
Not all points do -> expensive method
1
Example - sampling a
– Circle: p/4=78.5% good samples
p(x)
– Sphere: p/6=52.3% good samples
– Gets worst in higher dimensions
0
© Machiraju/Möller
1
Transforming between Distrib.
• Inversion Method --> transform uniform
random distribution to general distribution
• transform general X (PDF px(x))
to general Y (PDF py(x))
• Case 1: Y=y(X)
• y(x) must be one-to-one, i.e. monotonic
• hence
Py y  PrY  y x  PrX  x  Px x 
© Machiraju/Möller
Transforming between Distrib.
• Hence we have for the PDF’s:
py y dy  Py y   Px x   px x dx
1
dy 
py y     px x 
dx 
• Example: px(x) = 2x; Y = sinX

py y   cos x 
1
2x
2sin 1 y
px x  

cos x
1 y 2
© Machiraju/Möller
Transforming between Distrib.
• y(x) usually not given
• However, if CDF’s are the same, we use
generalization of inversion method:
y x   P
1
y
Px x

© Machiraju/Möller
Multiple Dimensions
• Easily generalized - using the Jacobian of
Y=T(X) p T x  J x 1 p x 
y
T
x
• Example - polar coordinates

x
r
JT x   
y

r
pr,   JT
x  r cos
y  r sin 
x 
cos  r sin  
  

y  
 sin  r cos 


1
px, y   rpx, y 
© Machiraju/Möller
Multiple Dimensions
• Spherical coordinates:
pr,,f  r sin px, y,z
2
• Now looking at spherical directions:
• We want to solid angle to be uniformly
distributed
d  sin ddf
• Hence the density in terms of f and :

p, f ddf  p d
p, f   sin p 
© Machiraju/Möller
Multidimensional Sampling
• Separable case - independently sample X
from px and Y from py: px, y  px xpy y
• Often times this is not possible - compute
the marginal density function p(x) first:
p
x 
 px, ydy
• Then compute conditional density function
px, y 
(p of y given x) py | x 
px 
 sampling with p(x) and p(y|x)
• Use 1D
© Machiraju/Möller
Sampling of Hemisphere
• Uniformly, I.e. p() = c
1

H2
p 
• Sampling  first:
p  

2p
2p
0
0
 p,f df  
1
c
2p
sin 
df  sin 
2p
• Now sampling in f:

p, f  1
pf |   

p  2p
© Machiraju/Möller
Sampling of Hemisphere
• Now we use inversion technique in order to

sample the PDF’s:
P   
 sin d  1 cos
0
P f |   


0
• Inverting these:

1
f
d 
2p
2p
  cos1 x1
f  2px 2
© Machiraju/Möller
Sampling of Hemisphere
• Converting these to Cartesian coords:
  cos1 x1
f  2px 2

x  sin  cos f  cos2px 2  1 x12
y  sin  sin f  sin 2px 2  1 x12
z  cos   x1
• Similar derivation for a full sphere

© Machiraju/Möller
Sampling a Disk
• Uniformly: px, y  
p
pr,   rpx, y  
pr 
• Sampling r first:

1

r
p
2p
 pr, d  2r
0
pr,  1
• Then sampling in : p | r  p r  2p



2
P r  r P  | r 
• Inverting the CDF:
2p
r

x


2
px
1
2

© Machiraju/Möller
Sampling a Disk
• Given method distorts size of compartments
x
• Better method
rx 
y

© Machiraju/Möller
Cosine Weighted Hemisphere
• Our scattering equations are cos-weighted!!
• Hence we would like a sampling
distribution, that reflects that!
• Cos-distributed p() = c.cos
1   p  d
1
c
H2
2p p 2

  c cos sin ddf
0 0
p 2
 2cp
 cos sind
0
© Machiraju/Möller
p
p, f  
1
p
cos sin 
Cosine Weighted Hemisphere
• Could use marginal and conditional
densities, but use Malley’s method instead:
• uniformly generate points on the unit disk
• Generate directions by projecting the points
on the disk up to the hemisphere above it
dw
dA/cos

rejected samples
dw cos
© Machiraju/Möller
dA
Cosine Weighted Hemisphere
•
•
•
•
•
Why does this work?
Unit disk: p(r, f) = r/p
Map to hemisphere: r = sin 
Jacobian of this mapping (r, f) -> (sin , f)
cos  0
Hence:
J x 
 cos
T


 0
p, f   JT pr, f  

© Machiraju/Möller

1
cos sin 
p
Performance Measure
• Key issue of graphics algorithm
time-accuracy tradeoff!
• Efficiency measure of Monte-Carlo:
F  
– V: variance
– T: rendering time
1
V F T F 
• Better algorithm if
– Better variance in same
 time or
– Faster for same variance
• Variance reduction techniques wanted!
© Machiraju/Möller

Russian Roulette
• Don’t evaluate integral if the value is small
(doesn’t add much!)
• Example - lighting integral
Lo p, o  

S
2
f r p, o, i Li p, i cosi d i
• Using N sample direction and a distribution
N
f r p, o, i Li p, i cosi
of p(i) 1

N i1
p i 
• Avoid evaluations where fr is small or 
close to 90 degrees

© Machiraju/Möller
Russian Roulette
• cannot just leave these samples out
– With some probability q we will replace with a
constant c
– With some probability (1-q) we actually do the
normal evaluation, but weigh the result
accordingly
F  qc

x q
F   1 q

otherwise
 c
• The expected value works out fine

E F   qc 
E F   1 q
 qc  E F 
 1 q 
© Machiraju/Möller
Russian Roulette
• Increases variance
• Improves speed dramatically
• Don’t pick q to be high though!!
© Machiraju/Möller

Stratified Sampling - Revisited
•
•
•
•
domain L consists of a bunch of strata Li
Take ni samples in each strata N f X
 i, j 
1
F


i
General MC estimator:
n i j1 pX i, j 
Expected value and variance (assuming vi is
the volume of one strata):
2
1
1
2
i  E f X i, j   L f 
xdx  i   L  f x  i  dx
vi
i
vi
i
• Variance for one strata with ni samples:

© Machiraju/Möller
 i2
ni
Stratified Sampling - Revisited

• Overall estimator / variance:
v i2 i2
2
V F   V  v i Fi  V v i Fi    v i V Fi   
ni
• Assuming number of samples proportional
to volume of strata - ni=viN:
1
V FN   v i i2
N
• Compared to no-strata (Q is the mean of f
over the whole domain L):
2
V F   1
v


 i i  vi i  Q
N
N
© Machiraju/Möller
Stratified Sampling - Revisited
V FN  

1
1
2
v

V FN  

i i
N
N


2
v

 i i  vi i  Q
• Stratified sampling never increases variance
• Right hand side minimized, when strata are

close to the mean of the whole function
• I.e. pick strata so they reflect local
behaviour, not global (I.e. compact)
• Which is better?
© Machiraju/Möller
Stratified Sampling - Revisited
• Improved glossy highlights
Random sampling
stratified sampling
© Machiraju/Möller
Stratified Sampling - Revisited
• Curse of dimensionality
• Alternative - Latin Hypercubes
– Better variance than uniform random
– Worse variance than stratified
© Machiraju/Möller
Quasi Monte Carlo
• Doesn’t use ‘real’ random numbers
• Replaced by low-discrepancy sequences
• Works well for many techniques including
importance sampling
• Doesn’t work as well for Russian Roulette
and rejection sampling
• Better convergence rate than regular MC
© Machiraju/Möller
Bias
b  E F   F
• If b is zero - unbiased, otherwise biased
• Example - pixel filtering

Ix, y 
 f x  s, y  tLs,tdsdt
• Unbiased MC estimator,
with
distribution
p
N

Ix, y  
1
f x  si , y  t i Lsi ,t i 

Np i1
• Biased (regular) filtering:

f x  s , y  t Ls ,t 

Ix, y  
 f x  s , y  t 
i
i
i
© Machiraju/Möller
i
i
i
i
i
Bias
• typically Np  i f x  si, y  ti 
• I.e. the biased estimator is preferred
• Essentially
trading bias for variance

© Machiraju/Möller
Importance Sampling MC
• Can improve our “chances” by sampling
areas, that we expect have a great influence
• called “importance sampling”
• find a (known) function p, that comes close
to the function we want to compute the
integral of,
1
f x 
• then evaluate: I   px  px dx
0
© Machiraju/Möller
Importance Sampling MC
• Crude MC:
n
F   i f x i 
i1
• For importance sampling, actually “probe”
a new function f/p. I.e. we compute our new
 to be:
estimates
f x i 
1
F 
N i1 px i 
N
© Machiraju/Möller
Importance Sampling MC
• For which p does this make any sense? Well
p should be close to f.
N
f x i 
• If p = f, then we would get F  1
1

N
px 
i1
i
• Hence, if we choose p = f/, (I.e. p is the
normalized distribution function of f) then
N
1
f x i 
we’d get:
1
F 
    0 f x dx
N i1 f x i  
© Machiraju/Möller
Optimal Probability Density
• Variance V[f(x)/p(x)] should be small
• Optimal: f(x)/p(x) is constant, variance is 0
p(x)  f (x) and  p(x) dx = 1
• p(x) = f (x) /  f (x) dx
• Optimal selection is impossible since it
needs the integral
• Practice: where f is large p is large
© Machiraju/Möller
Importance Sampling MC
• Since we are finding random samples distributed
by a probability given by p and we are actually
evaluating in our experiments f/p, we find the
variance of these experiments
to be:
2

 f x  
 px dx  I 2
  
px  
0
1
2
imp
1

0
f 2 x 
dx  I 2
px 
• improves error behavior (just plug in p = f/)
© Machiraju/Möller

Multiple Importance Sampling
• Importance strategy for f and g, but how to
sample f*g?, e.g.
Lo ( p, o )   f ( p, i , o )Li ( p, i ) cosi d i

• Should we sample according to f or
according to Li?
• Either one isn’t good
• Use Multiple Importance Sampling (MIS)
© Machiraju/Möller
Multiple Importance Sampling
Importance sampling f
Importance sampling L
© Machiraju/Möller
Multiple
Importance sampling

Multiple Importance Sampling
• In order to evaluate  f (x)g(x)dx
• Pick nf samples according to pf and ng
samples according to pg

• Use new MC estimator:
nf
ng

f (X i )g(X i )w f (X i )
f (Y j )g(Y j )w g (Y j ) 
1





n f  n g i1
p f (X i )
p
(Y
)
g
j
j1

• Balance heuristic vs. power heuristic:
b
n
p
(x)
n s ps (x)


s s
w s (x) 
w s (x) 
b
n
p
(x)
n
p
(x)
i i i © Machiraju/Möller i  i i 
MC for global illumination
• We know the basics of MC
• How to apply for global illumination?
– How to apply to BxDF’s
– How to apply to light source
© Machiraju/Möller
MC for GI - general case
• General problem - evaluate:
Lo ( p, o )   f ( p, i , o )Li ( p, i ) cosi d i

• We don’t know much about f and L, hence
use cos-weighted sampling of hemisphere in
order to find a i
• Use Malley’s method
• Make sure that o and i lie in same
hemisphere
© Machiraju/Möller
MC for GI - microfacet BRDFs
• Typically based on microfacet distribution
(Fresnel and Geometry terms not statistical
measures)
n
• Example - Blinn:D h   n  2 h  N 
• We know how to sample a spherical / power
distribution:
cos h  n 1 x1

f  2px 2
• This sampling is over h, we need a
distribution over i

© Machiraju/Möller
MC for GI - microfacet BRDFs
• This sampling is over h, we need a
distribution over i:
d i  sin i di dfi
d h  sin h dh dfh
• Which yields to
(using that h=2h and fh=fh):
  d df
d h sin
sin  h dh df h
sin  h
h
h
h



d i
sin  i d i df i sin 2 h 2d h df h 4 cos h sin h
1

4 cos h
© Machiraju/Möller
MC for GI - microfacet BRDFs
• Isotropic microfacet model:
ph  
p  
4 o   h 

© Machiraju/Möller
MC for GI - microfacet BRDFs
• Anisotropic model (after Ashikhmin and
Shirley) for a quarter disk:
 e  1 px 
x
1

f  arctan
tan

 e  1  2 


 y

cos h  x


e x cos2 f e y sin 2 f 1
2
1
• If ex= ey, then we get Blinn’s model

© Machiraju/Möller
MC for GI - Specular
• Delta-function - special treatment
1 N

N i1
    i 
hd ( o )
Li ( p, i ) cos  i
N
f ( p, i , o )Li ( p, i ) cos i 1
cosi
 
p i 
N i1
p i 
• Since p is also a delta function
p i       i 
• this simplifies to


hd (o )Li ( p,)
© Machiraju/Möller
MC for GI - Multiple BxDF’s
• Sum up distribution densities
N
1
p    pi  
N i1
• Have three unified samples - the first one
determines according to which BxDF to
distribute the spherical direction

© Machiraju/Möller
Light Sources
• We need to evaluate
– Sp: Cone of directions from point p to light (for
evaluating the rendering equation for direct
illuminations), I.e. i
Lo ( p, o ) 
 f ( p, , )L ( p, ) cos d
i
o
i
i
i
i

– Sr: Generate random rays from the light source
(Bi-directional Path Tracing or Photon
Mapping)
© Machiraju/Möller
Point Lights
•
•
•
•
Source is a point
uniform power in all directions
hard shadows
Sp:
– Delta-light source
– Treat similar to specular BxDF
• Sr: sampling of a uniform sphere
© Machiraju/Möller
Spot Lights
• Like point light, but only emits
light in a cone-like direction
• Sp: like point light, i.e. delta function
• Sr: sampling of a cone
p, f   p pf 
pf   1 2p
1  c 0
 max
sin d  c 1 cos max 
p   1 1 cos  max 
p   c

© Machiraju/Möller
Projection Lights
• Like spot light, but with a
texture in front of it
• Sp: like spot light, i.e. delta function
• Sr: like spot light, i.e. sampling of a cone
p, f   p pf 
pf   1 2p
1  c 0
 max
sin d  c 1 cos max 
p   1 1 cos  max 
p   c

© Machiraju/Möller
Goniophotometric
Lights
• Like point light (hard shadows)
• Non-uniform power in all
directions - given by distribution map
• Sp: like point-light
– Delta-light source
– Treat similar to specular BxDF
• Sr: like point light, i.e. sampling of a
uniform sphere
© Machiraju/Möller
Directional Lights
• Infinite light source,
i.e. only one distinct light direction
• hard shadows
• Sp: like point-light
– Delta function
• Sr:
– create virtual disk of the size of the scene
– sample disk uniformly (e.g. Shirley)
© Machiraju/Möller
Area Lights
• Defined by shape
• Soft shadows
• Sp: distribution over solid angle
– o is the angle between i
and (light) shape normal
– A is the area of the shape
coso
d i  2 dA
r
• Sr:
– Sampling over area of the shape

• Sampling distribution depends on the area of the
shape
1
px  
© Machiraju/Möller
A
Area Lights
• If v(p,p’) determines visibility:
Lo ( p, o ) 

f ( p, i , o )Li ( p, i ) cos i d i

  v p, p f ( p, i , o )Li ( p, i ) cosi
A
cos o
dA
2
r
coso
d i  2 dA
r
1
• Hence: px  
A
1
cos  o
Lo ( p, o ) 
v p, p f ( p, i , o )Li ( p, i ) cos i
2
px 

r
A
 2 v  p, p f ( p, i , o )Li ( p, i ) cos i cos o
r
© Machiraju/Möller
Spherical Lights
• Special area shape
• Not all of the sphere is visible
from outside of the sphere
• Only sample the area, that is
visible from p
• Sp: distribution over solid angle
– Use cone sampling
sin  max
r

pc
• Sr: Simply sample a uniform sphere
© Machiraju/Möller

c
r
p
Infinite Area Lights
• Typically environment light (spherical)
• Encloses the whole scene
• Sp:
– Normal given - cos-weighted sampling
– Otherwise - uniform spherical distribution
• Sr:
– uniformly sample sphere at two points p1 and p2
– The direction p1-p2 is the uniformly distributed
ray
© Machiraju/Möller
Infinite Area Lights
Area light + directional light
Morning skylight
Midday skylight
© Machiraju/Möller
Sunset environment map