Monday, June 29, 2015

Ray tracing: refraction

As my new ray tracer is spectral I am rewriting some of the basics from scratch including dielectric code.   I figured I could cut and paste some stuff out of ancient "ray tracing news" but ran into some problems with normal vectors.   So I am re-deriving this for my own benefit.
Sketch created with limnu iphone app from limnu.com
All renderers have some convention about which way surface normals point and whether that depends on incoming rays.   I am a fan of surfaces being boundaries between materials, and a dielectric has a surface normal N pointing to the side with refractive index n.   One thing that means is that in the diagram N might be pointing "down" along the dotted line.

First what about reflection?   It is usually the classic formula:

S = D - 2*N*dot(D, N)

Does that formula still work if N points along the dotted line?   Happily, yes-- try substituting -N in there and you get a "-" on the N and the dot product and they cancel.   Also, D need not be a unit vector.  I rarely require that in my ray tracers so I can use instancing if I want.   N however does need to be a unit vector.

Now what about refraction?   It obeys Snell's law:

n sin(theta) = n' sin(theta')

We can use the same trick as in reflection where we use N and D as a 2D basis.   First let's make life easy on ourselves and assume D is a unit vector and N points "up" like in the picture, and make an orthogonal basis (let's change the sin() and cos() with c, c', s, s')

D = -cN + sQ

Here Q is to the right and perpendicular to N.   We don't know what is is yet, but it we know what it is by rearranging the terms above to be get

Q = (D +  cN) / s


We also know that

T = -c'N + s'*Q

We know c = -dot(D,N), and s = sqrt(1 - c^2)

Further we know that s' = (n/n') s

Right there is a place we can get into trouble-- n can be high (over 2 for diamond) so if the more dense medium is where D is (the top medium), we could have s > 1.   In that case we get all reflection and no refraction (total internal reflection).   Here's a nice image of that:
The water acts as a perfect mirror due to total internal reflection (from wikipedia)
OK plug and chug the algebra:

T = -c'N + s'Q

T =  -c'N + s'(D +  cN) / s

We know s'/s = n'/n.

We also know c' = sqrt(1-s'^2) = sqrt(1 - (n/n')^2 (1-c^2) ).


T = (n'/n) D +(c(n'/n) - sqrt(1 - (n/n')^2 (1-c^2) ))

putting back in c = -dot(D, N) we get

T = (n'/n) D +(-dot(D, N)(n'/n) - sqrt(1 - (n/n')^2 (1-dot(D, N)^2) )) N

Ok!  That agrees with the ray tracing news formula.   Now what happens when N points down?  Only the last term is different:

T = (n'/n) D +(-dot(D, N)(n'/n) + sqrt(1 - (n/n')^2 (1-dot(D, N)^2) )) N

Well isn't that annoying?   That sign difference is just sgn(dot(D, N)).   I never know how to best code that sort of thing, but when in doubt, be very readable!

if (dot(D,N) > 0) temp_normal = -N    // now proceed with original formula


2 comments:

Eric Haines said...

Heh, I of course like seeing that the Ray Tracing News is a source for this bit of info. That's both great and sad - I would have hoped by now there would have been some grad student who got possessed by the ray tracing bug and scrawled all these sorts of bits (how to form a camera matrix, what the center of a pixel is labeled, reflection and refraction formulae, etc.) and put them on a nice web page somewhere.

I was surprised Phil Dutre's old document http://people.cs.kuleuven.be/~philip.dutre/GI/TotalCompendium.pdf doesn't have this vector formula. Morgan McGuire's Graphics Codex might well have it, http://graphicscodex.com/, but drat he's making dozens of dollars off it and flying around in his sky-yacht with the earnings.

Of course, now we need a new equation for metamaterials, which can have a negative index of refraction for some wavelengths: https://en.wikipedia.org/wiki/Negative_refraction

Haralambi Todorov said...

I was deriving the equations on my own and I think, I saw a small mistake in part of your equations, namely:
T = (n'/n)D +(c(n'/n) - sqrt(1 - (n/n')^2 (1-c^2) )) N
T = (n'/n)D +(-dot(D, N)(n'/n) - sqrt(1 - (n/n')^2 (1-dot(D, N)^2) )) N
and
T = (n'/n)D +(-dot(D, N)(n'/n) + sqrt(1 - (n/n')^2 (1-dot(D, N)^2) )) N
Shouldn't it be: n/n', because of n*sin(theta) = n'*sin(theta'), sin(theta') = (n/n')sin(theta)?

That aside, I want to really thank you for bringing this precious blog as well as the numerous books you've published on different topics of computer graphics. As a CS student, who is fascinated by the subject, I found them invaluable!

Wish you all the best,
Harry!