Raytrace 2: Assignment due Wednesday April 29 before class

This is the second part of ray tracing, and we are continuing to add cool features. For next week I'm going to ask you to do as many of the following as you feel up to, depending on your level of ambition:

I realize that there is a lot of stuff here. I decided to err on the side of giving you lots of options, so that you have lots of possibilities to explore of things to learn about ray tracing.

Shadows:

Implementing shadows is really quite easy to do with ray tracing. As we discussed in class, when you're at a single surface point and looping over light sources, fire a ray from the surface point into the direction of each light source. If the ray toward a light source hits any objects, then don't illuminate that surface point by that light source (ie: in your Phong shader, don't add in any Diffuse or Specular contributions for that light source).

There is one little trick though. You don't want the shadow ray to accidentally hit the object itself, due to numerical inaccuracy. So to play it safe, you want to start the ray just a little ways along the direction to the light source. Formally, here is how you form the shadow ray to a point light source:

w' = normalize(L - S)
v' = S + ε w'
where ε is an appropriately small number such as 0.001.

If your light source is just a directional light at infinity, rather than an actual point in the scene, then it's even easier, since then you can just set w' = L

Reflecting a ray off a surface:

To do reflections you will want to be able to bounce a ray (v + t w) off a surface to create a reflecting ray (v' + t w'). As we discussed in class, given surface normal vector n you can compute w' by:

w' = w - 2 (wn) n

and then the reflecting ray's origin can be given by moving just a little off of the surface point S:

v' = S + ε w'

Refracting a ray into a surface:

We didn't go into detail about refraction in class, other to point out that it is implemented by Snell's Law:

n1 sin(θ1) = n2 sin(θ2)

where n1 and n2 are the respective indices of refraction of the medium that the ray comes from and the medium that the ray passes into.

It turns out not to be all that hard to implement refraction. To compute the direction of the exit ray w2, take the dot product of entering ray direction w1 with the surface normal n. This will let you decompose w1 into two components c1 and s1 which are, respectively, perpendicular to and parallel to the surface.

The vector c1 = (w1n) n is the component of w1 perpendicular to the surface. This component has length cos(θ1).

The vector s1 = w1 - c1. is the component of w1 parallel to the surface. This component has length sin(θ1).

You can then just adjust the lengths of these two component vectors to cos(θ2) and sin(θ2), respectively. When you sum the length-adjusted component vectors, you will have w2.

As we discussed in class, the index of refraction in air is essentially 1.0, and the index of refraction in most optically clear solids and fluids such as water or glass is around 1.5.

Ray tracing to convex polyhedra:

As we said in class, to implement ray tracing to convex polyhedra:

Using pre-exisitng polyhedral models:

Although we didn't cover this in class, you may already have a convex polyhedron described as vertices and faces (eg: a pyramid or rectangular prism or icosohedron), and wish to convert this shape into an intersection of half spaces so you can ray trace to it. You would do this by taking each of your shape's polygonal faces in turn, and computing the half space for that face. You can figure out the linear inequality (a,b,c,d) that goes through the face of any polygon, by considering any three of the polygon's vertices P,Q,R. Remember that we have been following the convention that the points P,Q,R go around counterclockwise, as seen from the outside (ie: from the positive-valued half-space) of the shape. The method is as follows:

  • Take the difference vectors u = Q - P and v = R - Q.

  • Take their cross product g = u × v, where cross product is defined by ( uyvz - uzvy , uzvx - uxvz , uxvy - uyvx ).

  • Vector g gives the direction into which the gradient function increases the fastest. So you can use the components of g as your a,b,c coefficients. Also, if you normalize g to unit length, you get the surface normal n.

    You then just need to choose the fourth coefficient d such that one of your polygon vertices, say Q, lies on plane defined by ax+by+cz+d=0. In other words: gxQx + gyQy + gzQz + d = 0, so you can just set d = -(gQ).

Transforming polyhedral shapes:

To transform, by some matrix M, a polyhedral shape to which we are ray tracing, we need to transform each of the polyhedron's component half-spaces (a b c d).

To do this, you are going to use exactly the same technique that you used for transforming surface normals. But I'm going to go over the math here again anyway, just to make sure it's clear.

In the discussion that follows, I'm going use transpose notation (x y z 1)T when referring to the homogeneous decription of a point in space, to make it clear that I am describing the point as a column vector.

When transforming a half-space (a b c d) by matrix M, the key thing to observe is that if you transform each point

(x y z 1)TM(x y z 1)T,

then your corresponding transformation (a b c d) → (a' b' c' d') has to ensure that the transformed point still lies on the same side of the transformed plane:

(a' b' c' d') (M(x y z 1)T) = (a b c d) (x y z 1)T.

In other words, if we ask the question "is this point inside the original half-space?", we need to get the same answer after we transform both the point and the half-space.

As we have discussed in class, the way to do this is to post-multiply by the matrix inverse, so that:

(a b c d) → ( (a b c d) M-1 ).

This works because linear transformations are associative:

( (a b c d) M-1 )   ( M(x y z 1)T ) =
(a b c d)   ( M-1 M )   (x y z 1)T =
(a b c d) (x y z 1)T

Ray tracing to shapes bounded by second order polynomial surfaces:

As we discussed in class this last week, a second order polynomial shape is defined by a second order polynomial. That's a polynomial where the sum of the powers of any term is not greater than two. Examples are x2 + y2 - 1 and xy + 3.

In contrast, a polynomial like x2z + y2 - 1 is not second order, because the sum of the powers in the first term x2z is too high: 2 + 1 = 3.

The most general possible second order polynomial has ten terms:

ax2 + by2 + cz2 + dxy + eyz + fzx + gx + hy + iz + j
We can define a solid shape as all the places where such a polynomial evaluates to less than or equal to zero. Points on the surface of the shape give a zero value, points inside of the shape give a negative value, and points on the outside of the shape give a positive value.

Given any point (x,y,z), you can compute the partial derivatives of this polynomial to get the gradient vector (a vector which points perpendicularly out of the surface). If you normalize this gradient vector to unit length, you get the surface normal, which you can use to do shading, reflection and refraction.

As you probably know from your freshman calculus class, you can get the x,y,z components of this gradient vector by taking the derivative of the polynomial with respect to x,y and z, respectively:

( 2ax+dy+fz+g , 2by+dx+ez+h , 2cz+ey+fx+i )

If you want to transform a second order polynomial by a linear transformation matrix, it is useful to describe the evaluation of the polynomial at any point (x y z 1)T by using the notation of linear transformations. As we discussed in class, the following notation (originally devised by Jim Blinn) does the trick, by forming the polynomial coefficients into a 4×4 matrix:

x   y   z   1
a   d   f   g
0   b   e   h
0   0   c   i
0   0   0   j
x
y
z
1

The above expression multiplies out to exactly ax2 + by2 + cz2 + dxy + eyz + fzx + gx + hy + iz + j.

Now that we have the equation in this form, we can transform the original point and see what happens. First we transform (x y z 1)T to (M (x y z 1)T). We also need to transform the point in the other place it appears - on the left side of the coefficients matrix. When you transpose vectors and matrices, you end up reversing the order that they multiply together, so (x y z 1) is transformed to ((x y z 1) MT).

Following the same logic that we used for transforming half-spaces, we need to insert a matrix both to the left and to the right of the coefficients matrix, so that the result of the transformed equation doesn't change:

As we discussed in class, we do this as follows:

x   y   z   1
MT
M-1T
a   d   f   g
0   b   e   h
0   0   c   i
0   0   0   j
M-1
M
x
y
z
1

Of the three terms above, the middle term gives the transformed coefficients matrix. After you do the above matrix multiplication, you will end up with a 4×4 matrix that has non-zero values in all 16 places. In order to turn this into the ten coefficients of the general second order polynomial ax2+ by2+...+j, you need to add the six values on the lower left of the matrix into their corresponding places in the upper right:

a d f g
D b e h
F E c i
G H I j
a d+D f+F g+G
0  b  e+E h+H
0  0   c  i+I
0  0   0   j

Once you have this method of transformation, then you can start with just a few primitive polynomials, and use matrix transformations to get al; the others:

Shape Equation of bounding surface Polynomial coefficients
Sphere x2 + y2 + z2 = 1 1 1 1 0 0 0 0 0 0 -1
Infinite cylinder x2 + y2 = 1 1 1 0 0 0 0 0 0 0 -1
Infinite cone x2 + y2 = z2 1 1 -1 0 0 0 0 0 0 0
Infinite paraboloid x2 + y2 = z 1 1 0 0 0 0 0 0 -1 0

Given these shapes, and matrix-transformed versions of them, you can create other shapes by boolean intersection. For example, as we discussed in class, you can create a unit cylinder by intersecting the infinite volume bounded by the cylindrical surface x2+y2=1 with the two half spaces z≥-1 and z≤1

Boolean intersections:

As we have discussed in class, given any convex shapes (such as spheres, cylinders, or boxes), it is very easy to do boolean intersection via raytracing. When you trace a ray to any convex shape, either the ray misses the shape, or else the ray intersects the shape along some interval [tin , tout] along the ray.

Given a collection of convex shapes, if the ray misses any one of them, then the ray misses their intersection. If the ray hits all of them, then the intersection is given by the interval [max(tin) , min(tout)], where max(tin) is the maximum over all the tin, and where min(tin) is the minimum over all the tout.

For example, you can create a unit length cylinder by intersecting the infinite cylinder x2 + y2 - 1 = 0 with the cube that goes from -1 to +1 in x, y and z. Then you can apply a matrix to the component shapes in order to transform this cylinder to any other finite length cylinder.

Boolean unions and differences (making non-convex shapes):

Non-convex shapes are more difficult to deal with than are convex shapes, because the same ray can enter and leave a shape multiple times, to it is not sufficient just to speak of a single interval [tin , tout]. Instead, we need to describe the places where the ray is inside the shape by a sequence of t values: [ t0 t1 t2 ... tn-1 ], where the even numbered elements represent places where the ray enters the shape, and the odd numbered elements represent places where the ray exits from the shape. For example, [0.3   1.4   2.7   3.5] represents two successive segments where the ray is inside the shape: a first segment between t=0.3 and t=1.4, followed by a second segment between t=2.7 and t=3.5.

This notation makes it possible to describe non-convex things like the result of boolean difference. For example, suppose shape A contains a hollow cavity with shape B, and suppose that a particular ray intersects shape A at [t(A)in , t(A)out], and that the same ray intersects shape B at [t(B)in , t(B)out].

If t(A)in <   t(B)in <   t(B)out <   t(A)in, then the boolean difference A - B is described along this ray by:

[ t(A)in   t(B)in   t(B)out   t(A)out ].

In general you can implement boolean difference by replacing:

[ t0   t2   t3 ... tn-1 ]   →   [ 0   t0   t1   t2 ... tn-1   ]

Notice that if t0 was already zero, then the first interval becomes null, and you can just drop the first two elements of the resulting sequence. Similarly, if if tn-1 was already ∞ then you can just drop the last two elements of the resulting sequence.

In order to implement boolean union, you only need to implement boolean intersection and boolean difference, since de Morgan's law gives us:

A ∪ B = ¬( ¬A ∩ ¬B )

The hardest part is implementing intersection for non-convex shapes. To intersect two non-convex shapes along a ray you need to merge together their two sequences. This is done by progressively moving along the two sequences, looking for values of t where the ray is inside both shapes, and using those values of t to construct a new sequence. Here is a sketch of the algorithm:

   Create an empty destination sequence.

   Starting from the beginning of both sequences,
   Repeat until at the end of both sequences:

      Choose the smallest remaining t value, from whichever of the
      two input sequences you find it.

      If this was an entering (even numbered) t:
         If this t value is inside the other shape:
            then append it to the destination sequence.

      If this was an exiting (odd numbered) t:
         If this t value is outside the other shape:
            then append it to the destination sequence.


HOMEWORK:

Your job for this assignment is to use as many different ray tracing techniques as you feel comfortable implementing, and come up with an interesting collection of reflecting, refracting and shaded shapes, rendered by ray tracing, to produce cool and compelling pictures or animations. Try to use a variety of shapes, and make sure to create something fun!