Wednesday, November 20, 2013

Scattering in a constant medium.

I had a sign error in this one.  Here's the derivation of the right answer.

Suppose we have a participating media like smoke or fog that has constant density, and we send a photon-like particle into it and it interacts with the volume randomly so we want to randomly decide where, if anywhere, it gets scattered/absorbed.

Suppose we parameterize a ray along which the particle travels with "t" going from 0 to infinity.  For any given small part of the ray we have the following:

if at t the particle has not scattered, the probability it scatters n the small dt length region on the interval [t,t+dt] is C*dt.  C is proportional to the density of the medium: double the density of particles and the probability of scattering is doubled (because dt is small so the particles can't "line up" and each might block the particle independently).

For the cumulative probability density function P(t) which is the probability the particle has scattered at or before t, we can see:

P(t + dt) = P(t) + (1-P(t))*C*dt

That last term has the (1-P(t)) as the probability the particle is "live" at t, so it has the potential to scatter.


C*(1-P(t)) = (P(t + dt) - P(t)) / dt

This right hand side is the definition of derivative so we have the differential equation:

P'(t) = C*(1-P(t))

This yields

P(t) = 1 - exp(-C*t)

If we have a random number r = random(0,1), we find the t of random scattering by solving:

r = 1 - exp(-C*t)

-C*t = log(1-r)

t = - (1/C)*log(1-r)

Note that r and 1-r have the same distribution (replacing "random(0,1)" with "1-random(0,1)" has the same random behavior) so this is often written

t = -(1/C)*log(r)

Random directions with cosine distribution

I messed up this formula in my course this week.  I promised to fix it here.  Here it is!

Suppose you want to emit ray from a Lambertian (ideal diffuse) surface.  This has the probability of a direction proportional to cosine of the angle with the normal.  If we use a spherical coordinate system with (theta,phi) where phi is the angle "around" the z-axis and theta is the angle "down from" the z axis, we can see from symmetry that all phi are equally likely so

phi = 2*PI*random(0,1)

Now theta you can't just treat as a 1D random variable with density proportional to cosine of theta.  The "differential measure" on the sphere is sin(theta)*dtheta*dphi.  The sine of theta goes along with the theta term so the pdf of theta as a 1d density function is proportional to sin(theta)cos(theta)

p(t) = k*sin(t)*cos(t)    // t is theta

We know INT p(t) dt = 1, with t ranging 0 to PI/2 (hemisphere).  We will have to compute the cumulative probability density function and can normalize then.  The key bit is this: if we call r = random(0,1), we can get a Theta by solving for:

r = INT_0^Theta p(t) dt

The basic idea behind that is that if r = 0.7 for example, we want 70% of the "mass" of the density function to be below the Theta we choose.

r = INT_0^Theta k*sin(t)*cos(t)

r = 0.5*k*sin^2(Theta)

given Sin goes 0 to 1 the normalization falls out easily:

sin^2(Theta) = r

In spherical coordinates we have z = cos(Theta), so a unit vector with cosine density has the Cartesian coordunates:

x = cos(2*PI*r1)*sqrt(r)
y = sin(2*PI*r2)*sqrt(r)
z = sqrt(1-r)