Monte Carlo Integration

Monte Carlo Integration

Digital Image Synthesis Yung-Yu Chuang

11/29/2005

with slides by Pat Hanrahan and Torsten Moller

Introduction

• The integral equations generally don’t have

analytic solutions, so we must turn to numerical methods.

• Standard methods like Trapezoidal integration or Gaussian quadrature are not effective for high-dimensional and discontinuous integrals.

= ) ω p,

(

L

L

( p , ω

)

i i

i i

o

, ω ) ( p, ω ) cos θ ω

ω p,

2

f ( L

d

∫

+

Simple integration

Trapezoidal rule

Randomized algorithms

• Las Vegas v.s. Monte Carlo

• Las Vegas: gives the right answer by using randomness.

• Monte Carlo: gives the right answer on the average. Results depend on random numbers used, but statistically likely to be close to the right answer.

Monte Carlo integration

• Monte Carlo integration: uses sampling to estimate the values of integrals. It only

requires to be able to evaluate the integrand at arbitrary points, making it easy to implement

and applicable to many problems.

• If n samples are used, its converges at the rate of O(n^-1/2). That is, to cut the error in half, it is necessary to evaluate four times as many

samples.

• Images by Monte Carlo methods are often noisy.

Monte Carlo methods

Basic concepts

• X is a random variable

• Applying a function to a random variable gives another random variable, Y=f(X).

• CDF (cumulative distribution function)

• PDF (probability density function): nonnegative, sum to 1

• canonical uniform random variable ξ (provided by standard library and easy to transform to

other distributions)

dx x x dP

p ( )

) ( ≡

} Pr{

)

(x X x

P ≡ ≤

Discrete probability distributions

Continuous probability distributions

Expected values

• Average value of a function f(x) over some distribution of values p(x) over its domain D

• Example: cos function over [0, π], p is uniform

[ ]

∫

p f x f x p x dx

E ( ) ( ) ( )

[ ]

∫

cos )

cos(x x dx

E _p

Variance

• Expected deviation from the expected value

• Fundamental concept of quantifying the error in Monte Carlo methods

[ ]

( )

[

]

)]

(

[ f x E f x E f x

V = −

Properties

[

]

[

]

E =

[ ]

∑

∑

⎦

⎢ ⎤

⎣

⎡

i

i i

i E f X

X f

E ( ) ( )

[

]

[

]

V =

[

]

[ (

)

]

^[

^]

V = −

Monte Carlo estimator

• Assume that we want to evaluate the integral of f(x) over [a,b]

• Given a uniform random variable Xi over [a,b], Monte Carlo estimator says that the expected value E[F_N] of the

estimator F_N equals the integral

∑

= − ^N

i

i

N f X

N a F b

1

) (

General Monte Carlo estimator

• Given a random variable X drawn from an arbitrary PDF p(x), then the estimator is

• Although the converge rate of MC estimator is O(N^1/2), slower than other integral methods, its converge rate is independent of the dimension, making it the only practical method for high

dimensional integral

∑

= ^N

i i

i

N p X

X f

F N

1 ( )

) (

1

Choosing samples

• How to sample an arbitrary distribution from a variable of uniform distribution?

– Inversion – Rejection – Transform

Inversion method

Inversion method

• Compute CDF P(x)

• Compute P^-1(x)

• Obtain ξ

• Compute X_i=P^-1(ξ)

• Compute CDF P(x)

• Compute P^-1(x)

• Obtain ξ

• Compute X_i=P^-1(ξ)

Example: exponential distribution

, for example, Blinn’s Fresnel term

Example: power function

Sampling a circle

Sampling a circle

Rejection method

Rejection method

• Sometimes, we can’t integrate into CDF or invert CDF

• Rejection method is a fart-throwing method without performing the above steps

1. Find q(x) so that p(x)<cq(x) 2. Dart throwing

a. Choose a pair (X, ξ), where X is sampled from q(x) b. If (ξ<p(X)/cq(X)) return X

• Essentially, we pick a point (X, ξcq(X)). If it lies beneath p(X) then we are fine.

Example: sampling a unit sphere

void RejectionSampleDisk(float *x, float *y) { float sx, sy;

do {

sx = 1.f -2.f * RandomFloat();

sy = 1.f -2.f * RandomFloat();

} while (sx*sx + sy*sy > 1.f)

*x = sx; *y = sy;

}

π/4～78.5% good samples, gets worse in higher dimensions, for example, for sphere, π/6～52.3%

Transforming between distributions

• Transform a random variable X from distribution p_x(x) to a random variable Y with distribution

p_y(x)

• Y=y(X), y is one-to-one, i.e. monotonic

• Hence,

• PDF:

) ( }

Pr{

)}

( Pr{

)

(

P_y

= ≤ = ≤ =

dx x dP dx

y

dP_y

( )

( )

=

)

(x

dx dy dy

y dP

dx y dy

p_y ^y

( )

)

( =

) ( )

(

1

x dx p

y dy

p_{y x}

−

⎟ ⎠

⎜ ⎞

⎝

= ⎛

Example

X Y

x x

p_x

sin 2 )

(

=

=

2 1 1

1 sin 2

cos ) 2

( )

(cos )

(

y y x

x x p

x y

p_{y x}

= −

=

=

Transform

• Given X with p_x(x) and p_y(y), try to use X to generate Y.

Multiple dimensions

Multiple dimensions

Multidimensional sampling

Sampling a hemisphere

Sampling a hemisphere

Sampling a hemisphere

Sampling a disk

Sampling a disk

Shirley’s mapping

Sampling a triangle

Sampling a triangle

Multiple importance sampling

Multiple importance sampling

Multiple importance sampling

Monte Carlo for rendering equation

• Importance sampling: sample ω_i according to BxDF f and L (especially for light sources)

• If we don’t need anything about f and L, use

cosine-weighted sampling of hemisphere to find a sampled ω_i

= ) ω p,

(

L

L

( p , ω

)

i i

i i

o

, ω ) ( p, ω ) cos θ ω

ω p,

2

f ( L

d

∫

+

Cosine weighted hemisphere

θ ω ) cos

( ∝

p

∫

=

( )

1

p ω d ω

∫ ∫

=

θ θ θ φ

0 0

2

cos sin

1 c d d

φ θ

θ

ω d d

d = sin

∫

=

0

cos sin

2

1 c π

θ θ d θ

π¹

= c

θ π θ

φ

θ 1 cos sin )

,

( =

p

Cosine weighted hemisphere

θ π θ

φ

θ 1 cos sin )

,

( =

p

θ θ

θ φ

θ π θ

θ

1 cos sin 2 cos sin sin 2 )

(

0

= =

= ∫ ^d

p

π θ

φ θ θ

φ 2

1 )

(

) , ) (

|

( = =

p p p

2

2 1 2 cos

) 1

( θ = − θ + = ξ

P

2

)

|

( ξ

π θ φ

φ = =

P

) 2

1 ( 2 cos

1

1

1

ξ

θ =

− 2 πξ

φ =

Cosine weighted hemisphere

• Malley’s method: uniformly generates points on the unit disk and then generates directions by projecting them up to the hemisphere above it.

Vector CosineSampleHemisphere(float u1,float u2){

Vector ret;

ConcentricSampleDisk(u1, u2, &ret.x, &ret.y);

ret.z = sqrtf(max(0.f,1.f - ret.x*ret.x - ret.y*ret.y));

return ret;

}

Cosine weighted hemisphere

• Why Malley’s method works

• Unit disk sampling

• Map to hemisphere here

•

) ( )

(

1

x dx p

y dy

p_{y x}

−

⎟ ⎠

⎜ ⎞

⎝

= ⎛

φ

r

p ( , ) =

) , (sin )

,

( r φ → θ φ

) , ( ),

,

( φ = θ φ

= r X

Y

Sampling reflection functions

Spectrum BxDF::Sample_f(const Vector &wo,

Vector *wi, float u1, float u2, float *pdf){

*wi = CosineSampleHemisphere(u1, u2);

if (wo.z < 0.) wi->z *= -1.f;

*pdf = Pdf(wo, *wi);

return f(wo, *wi);

}

For those who modified Sample_f, Pdf must be changed accordingly

float BxDF::Pdf(Vector &wo, Vector &wi) { return SameHemisphere(wo, wi) ?

fabsf(wi.z) * INV_PI : 0.f;

} This function is useful for multiple importance sampling.

Sampling microfacet model

Too complicated to sample according to f, sample D instead.

Sampling Blinn microfacet model

• Blinn distribution:

• Generate ω_h according to the above function

• Convert ω_h to ω_i

e h

h

e n

D ( )

2 ) 2

( + ⋅

= ω

ω π

1 1

cos θ

=

ξ

2 πξ

φ

=

i

h h

o o

i

ω ω ω ω

ω = − + 2 ( ⋅ )

Sampling Blinn microfacet model

• Convert half-angle Pdf to incoming-angle Pdf, that is, change from a density in term of ω_h to one in terms of ω_i

Sampling anisotropic microfacet model

• Anisotropic model (after Ashikhmin and Shirley) for the first quadrant of the unit hemisphere

φ

ω

ω ) ( 1 )( 1 ) ( )

(

e

e

n

D = + + ⋅

Sampling BSDF (mixture of BxDFs)

• Sample from the density that is the sum of individual densities

• Three uniform random numbers are used, the first one determines which BxDF to be sampled and then sample that BxDF using the other two random numbers

Sampling light sources

• Direct illumination from light sources makes an important contribution, so it is crucial to be

able to

– Sp: samples directions from a point p to the light – Sr: Generates random rays from the light source

(for bidirectional light transport algorithms such as bidirectional path tracing and photon

mapping)

Point lights

• Sp: delta distribution, treat similar to specular BxDF

• Sr: sampling of a uniform sphere

Spotlights

• Sp: the same as a point light

• Sr: sampling of a cone (do not consider the falloff)

Directional lights

• Sp: like point lights

• Sr: create a virtual disk of the same radius as scene’s bounding sphere and then sample the disk uniformly.

Area lights

• Defined by shape

• Add shape sampling functions for Shape

• Sp: uses a density with respect to solid angle from the point p

Point Sample(Point &P, float u1, float u2, Normal *Ns)

• Sr: generates points on the shape according to a density with respect to surface area

Point Sample(float u1, float u2, Normal *Ns)

Infinite area lights

• Sp:

– normal given: cosine weighted sampling – Otherwise: uniform spherical sampling

– Does not take directional radiance distribution into account

• Sr:

– Uniformly sample two points p₁ and p₂ on the sphere – Use p1 as the origin and p₂-p₁ as the direction

– It can be shown that p₂-p₁ is uniformly distributed (Li et. al. 2003)