• 沒有找到結果。

# Monte Carlo Integration

N/A
N/A
Protected

Academic year: 2022

Share "Monte Carlo Integration"

Copied!
59
0
0

(1)

### Monte Carlo Integration

Digital Image Synthesis Yung-Yu Chuang

11/29/2005

with slides by Pat Hanrahan and Torsten Moller

(2)

### 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.

o

o

e

o

i i

i i

o

2

i

s

(3)

(4)

(5)

### 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.

(6)

### 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.

(7)

(8)

### 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 ≡ ≤

(9)

(10)

(11)

### 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

=

D

p f x f x p x dx

E ( ) ( ) ( )

=

0π π1 = 0

cos )

cos(x x dx

E p

(12)

### Variance

• Expected deviation from the expected value

• Fundamental concept of quantifying the error in Monte Carlo methods

( ) ( ) 2

## ]

)]

(

[ f x E f x E f x

V = −

(13)

af (x)

aE

f (x)

E =

=

⎢ ⎤

i

i i

i E f X

X f

E ( ) ( )

af (x)

a2V

f (x)

V =

f (x)

E

f (x)

2

## ]

E

f (x)

2

V = −

(14)

### 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[FN] of the

estimator FN equals the integral

=

= − N

i

i

N f X

N a F b

1

) (

(15)

### 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(N1/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

(16)

### Choosing samples

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

– Inversion – Rejection – Transform

(17)

(18)

### Inversion method

• Compute CDF P(x)

• Compute P-1(x)

• Obtain ξ

• Compute Xi=P-1(ξ)

(19)

• Compute CDF P(x)

• Compute P-1(x)

• Obtain ξ

• Compute Xi=P-1(ξ)

### Example: exponential distribution

, for example, Blinn’s Fresnel term

(20)

(21)

(22)

(23)

(24)

### 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.

(25)

### 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%

(26)

### Transforming between distributions

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

py(x)

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

• Hence,

• PDF:

y Y y x X x P x

Py

x

dx x dP dx

y

dPy

x

px

dx dy dy

y dP

dx y dy

py y

1

x dx p

y dy

py x

(27)

X Y

x x

px

2 1 1

y y x

x x p

x y

py x

(28)

### Transform

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

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

(37)

(38)

(39)

(40)

(41)

(42)

(43)

### 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

o

o

e

o

i i

i i

o

2

i

s

(44)

2

H

2π π

0 0

2

2

0

π

π1

(45)

π

2

0

1

2

1

1

2

(46)

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

}

(47)

### Cosine weighted hemisphere

• Why Malley’s method works

• Unit disk sampling

• Map to hemisphere here

1

x dx p

y dy

py x

πr

(48)

### 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.

(49)

### Sampling microfacet model

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

(50)

### Sampling Blinn microfacet model

• Blinn distribution:

• Generate ωh according to the above function

• Convert ωh to ωi

e h

h

1 1

h

e+

2

h

ωo ωh ω

i

h h

o o

i

(51)

### 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

(52)

### Sampling anisotropic microfacet model

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

φ

φ

cos2 sin2

h

x

y h

ex ey

+

(53)

### 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

(54)

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

(55)

### Point lights

• Sp: delta distribution, treat similar to specular BxDF

• Sr: sampling of a uniform sphere

(56)

### Spotlights

• Sp: the same as a point light

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

(57)

### 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.

(58)

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

(59)

### 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 p1 and p2 on the sphere – Use p1 as the origin and p2-p1 as the direction

– It can be shown that p2-p1 is uniformly distributed (Li et. al. 2003)

Biases in Pricing Continuously Monitored Options with Monte Carlo (continued).. • If all of the sampled prices are below the barrier, this sample path pays max(S(t n ) −

Again, the Vinaya instructs that one who using bronze bowl (other bronze utensils are not committed `Duskrta`, a most venial sin that the transgressor would

Directed numbers 2.1 understand the concept of directed numbers 9 Students are required to represent the directed numbers on the number line.. Students are required to

Population: the form of the distribution is assumed known, but the parameter(s) which determines the distribution is unknown.. Sample: Draw a set of random sample from the

Then, it is easy to see that there are 9 problems for which the iterative numbers of the algorithm using ψ α,θ,p in the case of θ = 1 and p = 3 are less than the one of the

y A stochastic process is a collection of &#34;similar&#34; random variables ordered over time.. variables ordered

Biases in Pricing Continuously Monitored Options with Monte Carlo (continued).. • If all of the sampled prices are below the barrier, this sample path pays max(S(t n ) −

Biases in Pricing Continuously Monitored Options with Monte Carlo (continued). • If all of the sampled prices are below the barrier, this sample path pays max(S(t n ) −