CSC418 - Notes
Topic 9: RayTracing

Navigation: back up next notes exercises

Key concepts & Readings

OpenGL

to be written

Notes

Introduction

In ray tracing, a ray of light is traced in a backwards direction. That is, we start from the eye or camera and trace the ray through a pixel in the image plane into the scene and determine what it hits. The pixel is then set to the color values returned by the ray. Each ray that hits an object can spawn other rays (reflected ray & refracted ray).


Ray tracing tree


raytrace( ray ) {
     find closest intersection
     cast shadow ray, calculate colour_local
     colour_reflect = raytrace( reflected_ray )
     colour_refract = raytrace( refracted_ray )
     colour = k1*colour_local + k2*colour_reflect
                                + k3*colour_refract
     return( colour )
}

Comments

Spawned Rays

Reflected Ray

Look at geometry of reflection

Note that for our rays: Rin = -Rd
The projection of Rin onto N is N cos (q) so
Rout = N cos(q) + S
Rin + S = N cos(q)
S = N cos(q) - Rin
Rout = 2Ncos(q) - Rin = 2N(N · Rin) - Rin

Refracted Ray

na sin(qin) = nb sin(qtransmitted)
qtransmitted = arc sin(na sin(qin) / nb)
to be written

Shadow Ray

Object Intersection

Reference: A. Glassner, etal., An Introduction to Ray Tracing, Academic Press, 1989.

Ray-Sphere Intersection

A ray is defined by: R(t) = R0 + t * Rd , t > 0 with R0 = [X0, Y0, Z0] and Rd = [Xd, Yd, Zd]

A sphere can be defined by its center and radius with Sc = [xc, yc, zc] So a sphere of radius Sr is:

S = the set of points[xs, ys, zs], where (xs - xc)2 + (ys - yc)2 + (zs - zc)2 = Sr2

To solve algebraically, substitute the ray equation into sphere equation and solve for t.

For ray:
X = X0 + Xd * t
Y = Y0 + Yd * t
Z = Z0 + Zd * t

putting X, Y, Z into the sphere equation for Xs, Ys, Zs

(X0 + Xd * t - Xc)2 + (Y0 + Yd * t - Yc)2 + (Z0 + Zd * t - Zc)2 = Sr^2

or A*t^2 + B*t + C = 0
with: A = Xd^2 + Yd^2 + Zd^2
B = 2 * (Xd * (X0 - Xc) + Yd * (Y0 - Yc) + Zd * (Z0 - Zc))
C = (X0 - Xc)^2 + (Y0 - Yc)^2 + (Z0 - Zc)^2 - Sr^2

Note: If |Rd| = 1 (normalized), then A = 1. So we can compute Sr^2 once.
So with A = 1, the solution of the quadratic equation is

t0, t1 = (- B + (B^2 - 4*C)^1/2) / 2 where t0 is for (-) and t1 is for (+)

If the discriminant is < 0.0 then there is no real root and no intersection. If there is a real root (Disc. > = 0.0) then the smaller positive root is the closest intersection point. So we can just compute t0 and if it is positive, then we are done, else compute t1. The intersection point is: 

Ri = [xi, yi, zi] = [x0 + xd * ti ,  y0 + yd * ti,  z0 + zd * ti]

Unit N at surface
SN = [(xi - xc)/Sr,   (yi - yc)/Sr,   (zi - zc)/Sr]

Note: For code optimization, "*" is faster than "/" so may compute 1/Sr and multiply. Look at the algorithm analysis:

Summary of steps:
compute A, B, C: 7 (*), 8 (+/-)
compute discriminant: 2 (*), 1 compare
compute t0 and determine if positive: 1 (-), 1 (*) /(/), 1 sqrt, 1 compare
possibly compute t1 and determine if positive: 1 (-), 1 (*), 1 compare
compute intersection point: 3 (*), 3 (+)
compute SN: 3 (-), 3 (*)
Total worst case = 17 (+/-), 17 (*), 1 sqrt, 3 compares.
Best case = steps 1, 2 = 9 (+/-), 9 (*), 1 compare.
Best case hit: 1, 2, 3, 5, 6 --> 16 (*), 16 (+/-), 1 sqrt, 3 compare.

Ray-Plane Intersection

A plane is defined by the equation: Ax + By + Cz + D = 0, or the vector [A B C D]. A, B, and C, define the normal to the plane. If A2 + B2 + C2 = 1 then the unit normal Pn = [A B C]. If A, B, and C define a unit normal, then the distance from the origin [0 0 0] to the plane is D.

A ray is defined by: R0 = [X0, Y0, Z0]

Rd = [Xd, Yd, Zd]

so R(t) = R0 + t * Rd , t > 0

To determine if there is an intersection with the plane, substitute for R(t) into the plane equation and get:

A(X0 + Xd * t) + B(Y0 + Yd * t) + (Z0 + Zd * t) + D = 0

which yields:

t = -(AX0 + BY0 + CZ0 + D) / (AXd + BYd + CZd)

= -(Pn· R0 + D) / (Pn · Rd)

First compute Pn· Rd = Vd. If Vd = 0 (incident angle, q = 900) then the ray is parallel to the plane and there is no intersection (if ray in in the plane then we ignore it). If Vd > 0 then the normal of the plane is pointing away from the ray. If we use one-sided planar objects then could stop if Vd > 0, else continue.

Now compute second dot product V0 = -(Pn· R0 + D) and compute t = V0 / Vd . If t < 0 then the ray intersects plane behind origin, i.e. no intersection of interest, else compute intersection point:

Pi = [Xi Yi Zi] = [X0 + Xd * t Y0 + Yd * t Z0 + Zd * t]

Now we usually want surface normal for the surface facing the ray, so if Vd > 0 (normal facing away) then reverse sign of ray.

So algorithm analysis:

  1. Compute Vd and compare to 0: 3 "*"s, 2 "+"s, 1 comparison.
  2. Compute V0 and t and compare to 0: 3 or 4 "*"s, 3 "+"s, 1 comparison.
  3. Compute intersection joint: 3 "*"s, 3 "+"s.
  4. Compare Vprd to 0 and reverse normal: 1 comparison.

Total = 10 "*"s, 8 "+"s, 3 comparisons.

Example calculation:

Given a plane [1 0 0 -7] ( plane with x = -7)

Ray with R0 = [ 2 3 4 ], Rd = [ 0.577 0.577 0.577]

Compute Vd = Pn · Rd = 0.577 > 0 , therefore the plane points away from the ray

Next, compute V0 = -(Pn · R0 + D) = 5

t = 5/0.577 = 8.66 > 0 , therefore the intersection point is not behind the ray, so compute the coordinates of the intersection point:

Xi = 2 + 0.577 * 8.66 = 7

Yi = 3 + 0.577 * 8.66 = 8

Ri = [ 7 8 9 ] Zi = 4 + 0.577 * 8.66 = 9

Must reverse normal so it is N = [-1 0 0]

Ray-Triangle Intersection

to be written

Monte-Carlo Sampling

Hit or Miss Method

  1. Assume we have a function f(x) and we want to find the integral of f(x) from x = 0.0 to x = 1.0.
  2. Draw a unit square.
  3. Choose N random points within the square.
  4. For each point we determine if it lies below the f(x) curve, if it does then we increment a counter Hit.
  5. The area under the curve is given by Hit/N. N = Hit + Miss.

Quadrature Method

  1. Let G = òg(x) dx be the integral to be evaluated.
  2. Rewrite g(x) as f(x)g'(x) where f(x) is a probability distribtion function and g'(x) = g(x)/f(x). Note that for a uniform random distribution in the interval [a,b] the pdf = 1/(b-a).
  3. Sample f (x) for a xi, i.e. obtain some xi from f(x).
  4. For each such sample xi evaluate g'(x).
  5. Perform steps 3 and 4 N times.
  6. The average = 1/N SNi=1 g'(x). The estimate has an associated error, i.e. = G + error where error is approximately equal to (variance(g')/N)^1/2.

Exercises

to be written