Monte Carlo Integration

Why: we want to solve an integral, but it can be too difficult to solve analytically.

What&How: estimate the integral of a function by averaging random samples of the function's value.

Define the Monte Carlo estimator for the definite integral of given function $f(x)$

• Definite integral $\int^b_af(x)dx$
• Random variable $X_i\sim p(x)$
• Monte Carlo estimator $\displaystyle F_N=\frac{1}{N}\sum^N_{i=1}\frac{f(X_i)}{p(X_i)}$

Path Tracing

Whitted-Style Ray Tracing's Problem:

• In Glossy reflection
• Color Bleeding

Whitted-Style Ray Tracing is Wrong. But the rendering equation is correct. $L_o(p,\omega_o)=L_e(p,\omega_o)+\int_{\Omega^+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)d\omega_i$

But it involves: Solving an integral over the hemisphere, and Recursive execution.

A Simple Monte Carlo Solution

Suppose we want to render one pixel (point) in the following scene for direct illumination only

Abuse the concept of Rendering Equation a little bit $L_o(p,\omega_o)=\int_{\Omega^+}L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)d\omega_i$ (assume all directions are pointing outwards)

It's still just an integration over directions. So, we can solve it using Monte Carlo integration.

We want to compute the radiance at p towards the camera, Monte Carlo integration: $\displaystyle\large\int_a^bf(x)dx\approx\frac{1}{N}\sum^N_{k=1}\frac{f(X_k)}{p(X_k)}$ $X_k\sim p(x)$

$f(x)$: $\large L_i(p,\omega_i)f_r(p,\omega_i,\omega_o)(n\cdot\omega_i)$

$p(\omega_i)=1/2\pi$ (assume uniformly sampling the hemisphere)

So, in general

A correct shading algorithm for direct illumination.

shade(p, wo)
    Randomly choose N directions wi~pdf
    Lo = 0.0
    For each wi
        Trace a ray r(p, wi)
        If ray r hit the light
            Lo += (1/N) * L_i * f_r * cosine / pdf(wi)
    Return Lo

Introducing Global Illumination

shade(p, wo)
    Randomly choose N directions wi~pdf
    Lo = 0.0
    For each wi
        Trace a ray r(p, wi)
        If ray r hit the light
            Lo += (1/N) * L_i * f_r * cosine / pdf(wi)
        Else If ray r hit an object at q
            Lo += (1/N) * shade(q, -wi) * f_r * cosine / pdf(wi)
    Return Lo

Problem:

• Explosion of #rays as #bounces go up: #rays = $N^{\#bounces}$: From now on, we always assume that only 1 ray is traced at each shading point:

  shade(p, wo)
      Randomly choose ONE direction wi~pdf(w)
      Trace a ray r(p, wi)
      If ray r hit the light
          Return L_i * f_r * cosine / pdf(wi)
      Else If ray r hit an object at q
          Return shade(q, -wi) * f_r * cosine / pdf(wi)
  

  This is path tracing. (Distributed Ray Tracing if N != 1)

  But this will be noisy! No problem, just trace more paths through each pixel and average their radiance.

  Ray Generation:

  ray_generation(camPos, pixel)
      Uniformly choose N sample positions within the pixel
      pixel_radiance = 0.0
      For each sample in the pixel
          Shoot a ray r(camPos, cam_to_sample)
          If ray r hit the scene at p
              pixel_radiance += 1/N * shade(p, sample_to_cam)
      Return pixel_radiance
  

• The recursive algorithm will never stop: Cutting #bounces == cutting energy.

  Solution: Russian Roulette (RR).

  Previously, we always shoot a ray at a shading point and get the shading result Lo

  Suppose we manually set a probability P (0<P<1)

  With probability P, shoot a ray and return the shading result divided by P: Lo / P

  With probability 1-P, don't shoot a ray and you'll get 0

  In this way, you can still expect to get Lo:

  $E=P*(Lo/P)+(1-P)*0=Lo$

  shade(p, wo)
      Manually specify a probability P_RR
      Randomly select ksi in a uniform dist. in [0,1]
      If (ksi > P_RR) return 0.0;
      
      Randomly choose ONE direction wi~pdf(w)
      Trace a ray r(p, wi)
      If ray r hit the light
          Return L_i * f_r * cosine / pdf(wi) / P_RR
      Else If ray r hit an object at q
          Return shade(q, -wi) * f_r * cosine / pdf(wi) / P_RR
  

Now we already have a correct version of path tracing. But it's not really efficient.

Sampling the Light

Understanding the reason of being inefficient

there will be 1 ray hitting the light. So a lot of rays are "wasted" if we uniformly sample the hemisphere at the shading point.

Monte Carlo methods allows any sampling methods, so we can sample the light (therefore no rays are "wasted")

Assume uniformly sampling on the light: pdf = 1 / A (because $\int$ pdf dA = 1)

But the rendering equation integrates on the solid angle: Lo = $\int$ Li fr cos d$\omega$

Recall Monte Carlo Integration: Sample on x & integrate on x

• Need to make the rendering equation as an integral of dA
• Need the relationship between d$\omega$ and dA

Projected area on the unit sphere $\Large d\omega=\frac{dA\cos\theta^\prime}{||x^\prime-x||^2}$

Then we can rewrite the rendering equation as

Now we consider the radiance coming from two parts:

• light source (direct, no need to have RR)
• other reflectors (indirect, RR)

shade(p, wo)
    # Contribution from the light source.
    Uniformly sample the light at x' (pdf_light = 1 / A)
    L_dir = L_i * f_r * cosΘ * cosΘ' / |x' - p|^2 / pdf_light
    
    # Contribution from other reflectors.
    L_indir = 0.0
    Test Russian Roulette with probability P_RR
    Uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi)
    Trace a ray r(p, wi)
    If ray r hit a non-emitting object at q
        L_indir = shade(q, -wi) * f_r * cosΘ / pdf_hemi / P_RR
    
    Return L_dir + L_indir

How do we know if the sample on the light is not blocked or not

# Contribution from the light source.
L_dir = 0.0
Uniformly sample the light at x' (pdf_light = 1 / A)
Shoot a ray from p to x'
If the ray is not blocked in the middle
    L_dir = L_i * f_r * cosΘ * cosΘ' / |x' - p|^2 / pdf_light