Mattias Rabb showed me a really neat trick a few years ago that this post will share. He found it in a paper I had never seen by Coleman in 1967. Kudos to Mattias for not only finding this paper, but managing to find this trick in it (I would not have understood it had he not told me what was there). It gives a way to compute an unbiased random hitpoint for a light particle and a nonuniform density.

Recall that for a uniform density volume this can be done analytically. Given a canonical random number X (e.g., a call to drand48()), and a medium with interaction coefficient s (the probability of a collision in a small distance d is sd), we have distance = -log(1-X)/s (see Dutre et al., p 241 for a derivation of this classic result).

In a photon tracing program we would use this as follows:

X = drand48()

d = -log(1-X)/s

newposition = rayorigin + d*raydirection // assumes unit length direction

scatter or absorb

No suppose s(x,y,z) is a function of position like a call to noise or whatever. Usually we would resort to ray maching with small steps that assume s is locally constant. This is a pain in the same family as ray epsilons. But the Coleman paper has a lovely trick when s(x,y,z) is always less than a known smax. Pretend the medium is all constant smax, take a step as above, and sometimes reject the result:

d = 0

while (1) {

X = drand48()

d += -log(1-X)/smax

Y = drand48()

newposition = rayorigin + d*raydirection

if (s(newposition)/smax > Y) break

}

This takes a bunch of random steps and is unbaised (I would love to see a clean proof of this-- I don't follow Coleman's). I think it saved us something like 5X in our adjoint photon tracer.