Recursive Ray Tracing

by William Shoaff with lots of help

Contents

• Contents

Recursive Ray Tracing

• Ray tracing handles global reflections and refractions (transmissions), hidden surface removal, shadows, and other special effects
• We want to consider
  • How to build a ray tracer
  • How to make it efficient
  • How to improve the image quality with special effects
• Whitted presented the first general recursive ray tracing model in 1980 (others had used some of the ideas before)
• Ray tracing follows specularly reflected and refracted rays through a scene
• The rays are infinitely thin, so aliasing is a problem
• Ray tracing can be used as a basic technique for volume rendering
• Ray tracing is recursive; secondary rays are followed recursively from primary rays

• The basic algorithm fires an eye ray $\vec{E}$ through the center of a pixel to find the intersection with the closest object
• A (specular) reflected ray is followed from the intersection
• A (specular) refracted ray is followed from the intersection
• To calculate shadows, fire a shadow ray (or shadow feeler) $\vec{S}$ from the intersection point to each light source

• Reflection and refraction are handled by secondary rays from the point of intersection
• The (specular) reflection ray is given by
  
  $$\vec{R} = \vec{I}-2\vec{N}(\vec{N}\cdot \vec{I})$$
  
  
  where $\vec{N}$ and $\vec{I}$ are the unit surface normal and unit incident vector hitting the surface
• Note that $\vec{I}$ is first an eye ray, then for secondary rays it can be a reflected or refracted ray; its role changes as rays are cast from ray-object intersection points
• The (specular) refraction (transmission) ray is given by
  
  $$\vec{T} = \frac{\eta_1}{\eta_2}\vec{I}- \left[\cos\theta_t+\frac{\eta_1}{\eta_2}(\vec{N}\cdot \vec{I})\right]\vec{N}$$
  
  
  where $\eta_1$ and $\eta_2$ are the indices of refraction in the incident and transmitted media (e.g. air to water), and $\theta_t$ is the angle between the transmitted ray and the negative surface normal
• We can calculate $\cos\theta_t$ as
  
  $$\cos\theta_t=\sqrt{1-\left(\frac{\eta_1}{\eta_2}\right)^2(1-(\vec{N}\cdot\vec{I})^2)}$$
  
  

• Each reflection and refraction ray may spawn their own shadow, reflection and refraction rays, recursively generating a ray tree
• The ray tree can be built to any depth depending on
  • The ray does not hit an object
  • The amount of work (time) we are willing to spend
  • The amount of contributed illumination from the rays
• Shadow rays will be ignored for the moment

The Intensity Equation for Ray Tracing

• The ray tree is evaluated bottom-up with each node's intensity a function of its children's intensity
• At each intersection point
  
  $$I=I_{\mbox{local}}+k_{rg}I_{\mbox{reflected}}+k_{tg}I_{\mbox{transmitted}}$$
  
  
  where
  • $I_{\mbox{local}}$ is the intensity at the intersection
  • k_rg is the global reflectivity (between 0 and 1)
  • $I_{\mbox{reflected}}$ is the reflected ray's returned intensity
  • k_tg is the global transmissivity (between 0 and 1)
  • $I_{\mbox{transmitted}}$ is the transmitted ray's returned intensity

• The local intensity can be computed from any illumination model, e.g.
  
  $$I_{\mbox{local}} = I_ak_a + I_{p}\left[k_d(\vec{N}\cdot\vec{L}) + k_s(\vec{N}\cdot\vec{H})^{\alpha}\right]$$
  
  
  (assuming the light reflects light to the viewer)
• I, $I_{\mbox{local}}$, $I_{\mbox{reflected}}$, $I_{\mbox{transmitted}}$ are $\langle R,\,G,\,B\rangle$ vectors
• Many of the terms in depend on wavelength: $\eta_1$, $\eta_2$, k_a, k_d, k_s, k_rg, k_tg; they can have different values for each component of the intensity vector
• Distance attenuation factors may be necessary (especially with the transmitted term)
• Only specular reflected and transmitted rays are considered
• Diffuse reflection and spread of specular reflection is considered only in the local intensity term

Pseudocode for a Simple Ray Tracer

`;=3000 raytrace() $\{$    Select the eye position $(e_x,\, e_y,\,e_z)$;     Select window $(x_{\min},\,y_{\min}),\,(x_{\max},\,y_{\max})$;     for $(y = y_{\min};\, y \leq y_{\max};\, y\mbox{\rm ++})\: \{$         for $(x = x_{\min};\, x \leq x_{\max};\, x\mbox{\rm ++})\: \{$             eyeray = $\langle e_x - x - 1/2,\, e_y - y - 1/2,\,e_z\rangle$;             pixel = trace $(0,\, 1,\, \mbox{\it eyeray})$;         $\}$    $\}$$\}$

`;=3000 color trace $(\mbox{\it level},\, \mbox{\it weight},\, \mbox{\it ray})$ int level; float weight; Ray ray; $\{$    if $(\mbox{\it intersect}(\mbox{\it ray},\,\mbox{\it model},\,\mbox{\it hit}))$ $\:\{$        obj = hit.object;         $\vec{N} =$ obj.normal;         P= raypoint(ray, hit.t);         return $(\mbox{\it shade}(\mbox{\it level},\mbox{\it weight},\mbox{\it ray},\mbox{\it obj}, \vec{N}, P))$;     $\}$    else return $(\mbox{\it background})$; $\}$

`=3000 color shade(level, weight, ray, obj, $\vec{N}$, P) begin     color = obj.k_a*I_a;     for each light do         $\vec{L} =$ unit vector from P to light;         $\vec{V} =$ unit vector from P to viewer; $\vec{H}=\frac{\vec{L}+\vec{V}}{\Vert\vec{L}+\vec{V}\Vert}$;         if $((\vec{N}\cdot\vec{L})>0)$ and (ray is not in shadow)             color += $I_{p}(\vec{N}\cdot\vec{L})*\mbox{\it obj}.k_{d}I_{p}(\vec{N}\cdot\vec{H})^{\mbox{\it obj}.\alpha}*\mbox{\it obj}.k_{s}$;     endfor;     if (level+1 < maxlevel$\,)$ then         if (obj.k_rg*weight > minweight) then             $\vec{R} =$ ReflectedDirection $(\mbox{ray}, \vec{N})$;             rcolor = trace(level+1, weigth*obj.k_rg, $\vec{R})$;             color += obj.k_rg*rcolor; endif;         if (obj.k_tg*weigth > minweight) then             $\vec{T} =$ RefractedDirection $(\mbox{\it ray},\,\vec{N})$;             tcolor = trace(level+1, weigth*obj.k_tg, $\vec{T})$;             color += obj.k_tg*tcolor; endif;     endif;     return color; end color

Ray Casting For Intersections

• A basic ray tracer intersects the eye ray with every object in the scene
• The nearest object hit spawns secondary rays that are also intersected with every object in the scene (except the current one)
• It has been estimated that 95% of the time in ray tracing is spend on intersection tests
• Quadric surfaces are often used in ray tracing because ray-quadric surface intersections are easy to compute (via quadratic equation)
• Spheres (a quadric surface) is often used as a bounding volume
• Ray-polygon intersections are easy (but expensive for a fine polygon mesh)
• Ray-box intersections are also fairly easy, and (generalized) boxes make good bounding volumes
• Intersection tests for bicubic patches, fractals, surfaces of revolutions, and other objects have been developed
• Intersection tests are independent of the rest of the ray tracing program so new objects can be added once the intersection code has been written

Ray/Sphere Intersection - Algebraic Solution

• Ray-sphere intersection calculations are easy
• Let the sphere be given by
  
  (x-a)² + (y-b)² + (z-c)² = r²
  
  
  where $C=(a,\,b,\,c)$ is the center of the sphere and r is its radius
• Let the ray be given by
  
  $$\begin{eqnarray*}x & = & x_o + x_d t \\ y & = & y_o + y_d t \\ z & = & z_o + z_d t \end{eqnarray*}$$
  
  where $R_o=(x_o,\,y_o,\,z_o)$ is the ray's origin and $\vec{R}_d=\langle x_d,\,y_d,\,z_d\rangle$ is the unit direction vector for the ray

• The plugging parametric line into sphere equation yield the quadratic equation
  
  At² + Bt + C = 0
  
  
  where
  
  A = x_d² + y_d² + z_d²,
  
  
  
  
  B = 2[x_d (x_o - a) + y_d (y_o - b) + z_d (z_o - c)],
  
  
  and
  
  C = (x_o - a)² + (y_o - b)² + (z_o - c)² - r²
  
  

• If B²-4AC<0, the ray and sphere do not intersect
• If B²-4AC=0, ray grazes sphere
• If B²-4AC>0, the smallest positive t corresponds to the intersection
  
  $$t_{\pm} = \frac{-B\pm\sqrt{B^2-4AC}}{2A}$$
  
  

• If the intersection is $(x_i,\,y_i,\,z_i)$, then the surface normal at the intersection is
  
  $$\vec{N}=\left\langle \frac{x_i-a}{r},\,\frac{y_i-b}{r},\,\frac{z_i-c}{r}\right\rangle$$
  
  

Ray/Sphere Intersection - Geometric Solution

• A more efficient solution can be determined geometrically
• Find if the ray's intersection is outside the sphere, that is if
  
  $$v^2\equiv (\vec{R_oC}\cdot\vec{R_oC})=(a-x_o)^2+(b-y_o)^2+(c-z_o)^2 \geq r^2$$
  
  

• Find the closest approach of the ray to the sphere's center
  
  $$\begin{eqnarray*}t_{ca} & = & v\cos\theta \\ & = & v\frac{(\vec{R_oC}\cdot\vec... ...a-x_o,\,b-y_o,\,c-z_o\rangle\cdot\langle x_d,\,y_d,\,z_d\rangle \end{eqnarray*}$$
  
  (this assumes $\langle x_d,\,y_d,\,z_d\rangle$ is a unit vector)
• If the ray is outside and points away from the sphere it misses the sphere, that is, if
  
  $$v^2\geq r^2\qquad \mbox{and}\qquad t_{ca} < 0$$
  
  
  the sphere is missed

• If the origin is inside the sphere or $t_{ca}\geq 0$, find the ``squared distance'' from t_ca to the sphere's surface
  • This is the half chord distance squared
  • If the ray misses the sphere, it can be negative
  
  
  t_hc²=r²-(v²-t_ca²)
  
  

• If t_hc < 0, the ray misses the sphere
• If $t_{hc}\geq 0$, then
  
  $$t=t_{ca}-\sqrt{t_{hc}^2} \quad \mbox {if $v^2\geq r^2$}$$
  
  
  and
  
  $$t=t_{ca}+\sqrt{t_{hc}^2} \quad \mbox {if $v^2< r^2$}$$
  
  

• The code can be written as

`=3000          $t_{ca}=\vec{R_oC}\cdot \vec{R}_d$;          $t^2_{hc} = r^2-((\vec{R_oC}\cdot\vec{R_oC})-t^2_{ca})$;         if (t²_hc < 0) then no intersection;         else begin              $t=\sqrt{t^2_{hc}}$;              $I=R_o+(t_{ca}-t)\vec{R}_d$;         end

Computing Ray Intersection With A Polygon

• Assume the polygon is described by vertices V₀, ..., V_n-1, $n\geq 3$
• Let $(x_i,\,y_i,\,z_i)$ be the coordinates of vertex V_i
• The normal to the polygon is computed as
  
  $$\vec{N}=\vec{V_0V_1}\times \vec{V_0V_2}$$
  
  

• If P is in the plane of the polygon, then
  
  $$\vec{N}\cdot P + D = 0$$
  
  

• Let the ray be given by
  
  $$R(t) = R_o + \vec{R}_d t$$
  
  

• The intersection of the plane and the ray is at parameter value
  
  $$t=-\frac{D+\vec{N}\cdot R_o}{\vec{N}\cdot \vec{R}_d}$$
  
  

• This calculation requires 12 floating point operations and 3 tests
  • If $\vec{N}\cdot \vec{R}_d =0$ the intersection is rejected
  • If $t\leq 0$, the intersection is rejected
  • If a closer intersection has been found the intersection is rejected

Polygon Inside/Outside Testing

• If the intersection is not rejected, determine if it is inside the polygon
• The solution presented here is based on triangles, and works for convex polygons
• A polygon with n vertices can be broken up into n-2 triangles
• Given point P we can write
  
  $$\vec{V_0P} = \alpha\vec{V_0V_1} + \beta\vec{V_0V_2}$$
  
  

• Point P is inside triangle $(V_0,\,V_1,\,V_2)$ if
  
  $$\alpha \geq 0,\,\beta\geq 0,\,\mbox{and}\,\alpha+\beta \leq 1$$
  
  

• The above equation for P can be written as
  
  $$\begin{eqnarray*}x_P-x_0 & = & \alpha(x_1-x_0) + \beta(x_2-x_0) \\ y_P-y_0 & =... ...a(y_2-y_0) \\ z_P-z_0 & = & \alpha(z_1-z_0) + \beta(z_2-z_0) \end{eqnarray*}$$
  
  

• To reduce the system, project it onto either the xy, xz or yz planes
• If the polygon is perpendicular to one of these planes, it projection is a line (or if it is nearly perpendicular, its projection has little area)
• To make the projection as large as possible, find the dominant axis of the normal vector and project onto the plane perpendicular to this axis
  
  $$i_0 = \left\{\begin{array}{rr} 0 & \mbox{if $\vert N_x\vert ... ...\vert,\,\vert N_y\vert,\,\vert N_z\vert\}$} \end{array}\right.$$
  
  

• Let i₁ and i₂ be indices different from i₀

• Let $(u,\,v)$ denote the plane defined by i₁ and i₂
• The coordinates of $\vec{V_0P}$, $\vec{V_0V_1}$, $\vec{V_0V_2}$ are
  
  $$u_0=P_{i_1}-V_{0_{i_1}}\qquad u_1=V_{1_{i_1}}-V_{0_{i_1}}\qquad u_2=V_{2_{i_1}}-V_{0_{i_1}}$$
  
  
  
  
  $$v_0=P_{i_2}-V_{0_{i_2}}\qquad v_1=V_{1_{i_2}}-V_{0_{i_2}}\qquad v_2=V_{2_{i_2}}-V_{0_{i_2}}$$
  
  

• And the equations for $\alpha$ and $\beta$ reduce to
  
  $$\begin{eqnarray*}u_0 & = & \alpha u_1 + \beta u_2 \\ v_0 & = & \alpha v_1 + \beta v_2 \end{eqnarray*}$$
  
  
• Which has solution
  
  $$\alpha = \frac{\det\left(\begin{array}{cc} u_0 & u_2 \\ v_0... ...(\begin{array}{cc} u_1 & u_2 \\ v_1 & v_2 \end{array}\right)}$$
  
  

• Here's a more general algorithm
• Reduce dimensionality - project polygon and intersection point orthographically onto coordinate plane which yields largest projection
• Translate polygon so intersection point is at origin
• Draw ray from origin along positive u axis
• Count how many times it intersects polygon
• Even number of intersections $\Rightarrow$ outside; odd $\Rightarrow$ inside

`;=3000 Polygon-Inside-Outside-Test $\{$    for i=0 to N         project $(X_i,\,Y_i,\,Z_i)$ onto dominant plane             creating a list of vertices $(U_i,\,V_i)$;     translate intersection point to the origin;         call the new points $(U'_i,\, V'_i)$;     Number-Cross =0;     if V'₀ < 0 then Sign-Hold = -1;     else Sign-Hold =1;     for a = 0 to N-1          $b = (a+1) \bmod N$;         if V'_b < 0 then Next-Sign-Hold =-1;         else Next-Sign-Hold =1;         if Sign-Hold $\neq$ Next-Sign-Hold then             if U'_a >0 and U'_b >0 then                 Number-Cross += 1;             if U'_a >0 or U'_b > 0 then                 compute edge intersection with U' axis;                 if U'_a - V'_a(U'_b-U'_a)/(V'_b-V'_a) > 0 then                     Number-Cross += 1;         Sign-Hold = Next-Sign-Hold;     if Number-Cross is odd then         intersection is inside the polygon; $\}$

Ray-Box Intersection

• Boxes are useful for bounding volumes (they usually fit the objects tighter)
• Hierarchies of boxes (boxes inside of boxes) are often used to partition space
• Generalized boxes are formed from three pairs of parallel planes
• Here we consider upright rectanguluar boxes
• The ray-box intersection is just the Liang-Barsky clipping algorithm in three dimensions
• For each of x, y, and z in turn, find the near and far parameter value where intersections occur
• The largest near intersection and smallest far parameter values are retained

• Set $t_{near}=-\infty$ and $t_{far}=\infty$
• Let $\langle x_d,\,y_d,\,z_d\rangle$ the ray direction vector and $(x_{b1},\,y_{b1},\,z_{b1})$, $(x_{b2},\,y_{b2},\,z_{b2})$ opposite corners defining the box
• If x_d=0, the ray is parallel to the x planes, so if x_o < x_b1 or x_o>x_b2 no intersection
• Otherwise, the two intersections with the x planes of the box are $t_1=\frac{x_{b1}-x_o}{x_d}$ and $t_2=\frac{x_{b2}-x_o}{x_d}$
• If t₁ > t₂ swap $(t_1,\,t_2)$
• If $t_1 > t_{\mbox{near}}$, then $t_{\mbox{near}}=t_1$
• If $t_2 < t_{\mbox{far}}$,then $t_{\mbox{far}}=t_2$
• If $t_{\mbox{near}} > t_{\mbox{far}}$, then the box is missed
• If $t_{\mbox{far}} < 0$, the box is behind the ray, so it is missed
• Repeat for y and z pairs of parallel planes

Shadows in Ray Tracing

• Whether an intersection point is in shadow or not is determined by casting a ray to the light source
• If the shadow feeler intersects an object before the light, then the point is in shadow and $I_{\mbox{local}}$ is decreased (or set to 0)
• If there are n light sources then a total of n+2 rays must be spawned from the nearest intersection (a reflected, a refracted, and n shadow feelers)
• A light buffer can be used to accelerate the shadow testing
• Light buffers exploit the fact that point can be determined to be in shadow without finding the closest occluing object
• For each light source, six faces of pointers, surrounding the light point to potentially blocing objects

• The light buffers are constructed as a pre-processing step, before ray tracing begins
• The candidate lists are created by projecting each object onto the six faces (of a unit cube) surrounding the light
• If a direction cell is partially or totally covered by the projected object it is added to the candidate list
• For polygon objects this can be performed efficiently by a modified scan-line algorithm to the projected edges
• Once all lists are created, they are sorted into ascending order according to depth from the light

• Observations that simplify the candidate lists:
  • All polygons which face away from the light and are part of opaque solids can be culled
  • Any list containing exactly one polygon can be deleted (a polygon cannot occlude itself)
  • If a projection completely covers a direction cell, the candidate list can be marked as fully-occluded at the polygon depth and all polygons of greater depth can be eliminated
• To determine if a point on a given surface is in shadow:
  • Check the surface orientation with respect to the light source
  • If it faces away, the polygon is in shadow
  • Otherwise, retrieve the list of potentially occluding objects from the light buffer using the direction of the shadow feeler
  • The objects on the list are tested for intersection in order of increasing depth until
    • An occlusion is found (point in shadow)
    • An object is found whose depth exceeds the point being tested (point not in shadow)

Efficient Considerations in Ray Tracing

• Ray tracing can be very inefficient
  • 12 years passed between first published paper that used ray tracing (1968) and ray tracing started to be heavily used
  • This delay was primarily due to the cost (in time) required to ray trace scenes
• Three basic approaches to making ray tracing more efficient are:
  • To compute fewer intersection -- adaptive depth control and adaptive sampling
  • To compute intersections faster -- bounding volume, light buffers, space subdivision
  • To use generalized rays -- beams, cones, pencils
• Adaptive depth control allows the depth of the ray tree to change, getting deeper when the returned intensity is large and more shallow then the returned intensity is small
• Bounding volumes can be use quickly discard objects that can not be intersected by a ray
• Beam, cone, and pencil tracing use generalized rays

Adaptive Depth Control

• The percentage of a scene that consists of highly reflective and transparent surfaces is, in general, small
• It is inefficient to trace all rays to a maximum depth
• The ray tree can be pruned when little intensity will result from the pruned subtree
• Rays are attenuated as they pass through a scene
• Adaptive depth control uses properties of the materials the ray hits to determine this attenuation
  • Reflected rays are attenuated by the global reflection coefficient k_rg
  • Refracted rays are attenuated by the global refraction coefficient k_tg
• Hall and Greenberg report that adaptive depth control resulted in an average depth of 1.71, in their example scenes, where the maximum depth was set at 15

Bounding Volumes in Ray Tracing

• Bounding volumes, or extents, or enclosures are the most basic tool for ray tracing acceleration
• A bounding volume encloses a given object (or objects)
• It should be easier to compute the ray-bounding volume intersection than the ray-object intersection
• Whitted initially used spheres as bounding volumes since they are the simplest shapes to test for intersection
• Using simple bounding volumes substitutes simple intersection checks for more costly ones, but does not decrease the number of intersections
• The simplicity of the intersection test is not the only criterion; the void area of bounding volume should also be considered

• The cost associated with an object and its bounding volume can be defined as
  
  $$\mbox{cost}=n*B+m*I$$
  
  
  where
  • n is the number of rays tested against the bounding volume,
  • B is the cost of each bounding volume test
  • m is the number of rays that hit the bounding volume ($m\leq n$)
  • I is the cost of intersecting the objects within the bounding volume
• Assuming n and I are fixed, that is, the number of rays fired n depends on the number of pixels and the cost of intersecting objects I depends on the scene, neither of which can change
• We want to select bounding volumes making B and m small

• Reducing B generally increases m
• For example, B can be reduced by using sphere, but this loosens the fit of the bounding volume causing more intersections (increasing m)
• On the other hand, a tigher fit, decreases m but increases B
• Spheres, cylinders, axis aligned (upright) boxes, and oriented boxes are typically used as bounding volumes
• Each bounding volume can be given a complexity rating (spheres = low, oriented boxes = high)
• The tightness of fit can be measured by the projected void area

Hierarchical Bounding Volumes

• Hierarchical bounding volumes can attain a theoretical time complexity that is logarithmic in the number of objects
• Rubin and Whitted used bounding boxes to create a hierarchy of bounding volumes
• A bounding volume hierarchy is assumed to be a tree or directed acyclic graph (DAG), with an arbitrary branching factor
• Thus bounding volumes may enclose any number of other bounding volumes
• Leaf nodes are primitive objects (spheres, polygons, boxes, etc)

• Setting up the hierarchy may require human involvment
• The code below can be used to test bounding volume hierarchies for intersections with rays

`;=3000 bounding-volume-hierarchy-intersect(ray, node) Ray ray; Node node; $\{$    if (node is a leaf) then         intersect(ray, node.object);     else if $(\mbox{\it intersect(ray, boundingvolume(node))})$         for (each child of node)             bounding-volume-hierarchy-intersect(ray, child); $\}$

• Kay and Kajiya give desirable properties of hierarchical schemes
  • Any subtree should contain objects that are near each other; the deeper the subtree, the nearer the objects
  • The volume of each node should be minimal
  • The sum of the volume of all bounding volumes should be minimal
  • The tree construction should concentrate on the nodes nearer the root; this makes pruning more efficient
  • The time spend constructing the hierarchy should more than pay for itself in the time spend rendering the image

• Kay and Kajiya also introduce of the concept of parallel slabs to approximate the convex hull of an object
• The slabs better approximate the shape of the object; they have less void area
• At least three independent sets of parallel slabs are used to enclose a 3 dimensional object
• Bounding boxes use three sets of slabs; each set is parallel to one of the coordinate axes
• Slabs with arbitrary orientation are characterized by a normal vector $\vec{N}=\langle A,\,B,\,C\rangle$ and a displacement D (from the origin)
  
  Ax+By+Cz+D = 0
  
  

• Kay and Kajiya use a predetermined set of normals to save computation costs and storage
• As more normals are added to the set, the void area is reduced, but more storage and computation are required
• For an object, the slab bounding volume is determined by a (sub)set of the normals and corresponding pairs of $D_{\mbox{near}},\,D_{\mbox{far}})$ of displacements
• For polygons, the displacements are easy to find; for each normal $\vec{N}_i$, compute
  
  $$d_{ij} = (\vec{N}_i\cdot\langle x_j,\,y_j,\,z_j\rangle)$$
  
  
  where $(x_j,\,y_j,\,z_j)$ is the jth vertex of the polygon
• Then
  
  $$D_{\mbox{near}}=\min\{d_{i1},\ldots,d_{im}\}\qquad D_{\mbox{far}}=\max\{d_{i1},\ldots,d_{im}\}$$
  
  

• Spheres and other surfaces are more difficult to handle

Spatial Subdivision for Ray Tracing

• Bounding volume techniques select volumes based on a given set of objects
• Spatial subdivision selects sets of objects based on given volumes
• Spatial subdivision divides space into nonoverlapping regions and considers the objects in those regions
• Rather than check a ray against all objects or set of bounding volumes, we check whether the region in which the ray is traveling contains objects
• If the region is empty, the ray is followed into the next region
• If the region contains objects, the ray is tested for intersection with those objects

• The preprocessing step of subdividing space eliminate two major overheads associated with ray tracing
  • Rays do not need to be checked for intersection with all objects; only the closest object in the first region pierced by the ray needs to be found
  • Ray-object intersections are checked only for objects near the path of the ray
• Rendering time is constant (independent of the scene complexity) for spatial subdivision techniques

Auxiliary Data Structures for Spatial Subdivision

• All spatial subdivision approaches preprocess space setting up auxiliary data structures that contain information about objects occupying regions of space
• Rays are traced through this data structure
• Octrees establish a hierarchical data structure of nonoverlapping cubes in space
• Binary space partitioning trees use an indexing scheme for branches that differs from octrees
• Space can be partitioned into uniform sized regions or the regions can differ is size
• Nonuniform spatial subdivision discretize space into regions of varying size to conform to features of the environment
  • More subdivisions can be performed in densely populated regions
  • Fewer (or larger) regions cover space that have few or no objects

Octrees

• Glassner first applied octrees to ray tracing
• Octrees are hierarchical data structures that can be efficiently indexed into regions of space
• Rectangular volumes are recursively subdivided into 8 sub-octants until at the leaf (or voxel) some simplicity criterion is met
• The parent node, labeled 1, encloses the ``world''
• The voxels (leaf nodes) are marked as empty, full, or mixed

• Each (non-empty) voxel contains
  • A name (numeric identifier),
  • A subdivion flag,
  • A list of objects whose surfaces penetrate that volume
• Each ray that pierces the voxel is tested for intersection against the objects in this list

• Given a point $(x,\,y,\,z)$, we can find the node it belongs to using the code below

`;=3000 findnode $(x,\,y,\,z)$$\{$    node =1;     while ( node-subdivided = TRUE ) do         if $(x > \mbox{\it node-center-x})$             if $(y > \mbox{\it node-center-y})$                 if $(z > \mbox{\it node-center-z})$                     node $= 10*\mbox{\it node}+6$;                 else node $= 10*\mbox{\it node}+2$;             else                 if $(z > \mbox{\it node-center-z})$                     node $= 10*\mbox{\it node}+8$;                 else node $= 10*\mbox{\it node}+4$;         else             if $(y > \mbox{\it node-center-y})$                 if $(z > \mbox{\it node-center-z})$                     node $= 10*\mbox{\it node}+5$;                 else node $= 10*\mbox{\it node}+1$;             else                 if $(z > \mbox{\it node-center-z})$                     node $= 10*\mbox{\it node}+7$;                 else node $= 10*\mbox{\it node}+3$;     end while; return node $\}$

• Given the list of objects in the ``world'' here's how Glassner decides to subdivide a given node
  • To create a child of a parent, for each object that passes through the child's volume, add the object to the child's list
  • We can determine if an object passes through a volume by intersecting the object with each of the 6 planes bounding the voxel
  • If any intersection point lies within the square of the face, the object is placed on the child's list
  • An additional test for proper containment, of the object within the voxel, is made by testing if a single point on the object is inside the voxel
  • Once the list of objects is constructed, if the list has too many objects, the child node is recursively subdivided

• The intersection of rays with the environment using octrees can be performed by the code below
• We must also find a point in the next voxel if the ray hits no objects in the current voxel
• The exit point of the ray and the current voxel can be calculated using the standard ray-box intersection test
• The next point can be found by moving a small distance from this exit point
• The length of the shortest voxel edge in the entire octree can be used for this purpose, call this length minlen
• If the exit point is interior to a face, move the exit point directly away from this face by minlen/2
• If the exit point is along an edge or vertex, move the exit point by minlen/2 for each face containing the exit point

`;=3000 octree-intersect (ray $\vec{R})$ $\{$     $Q=\vec{R}.\mbox{\it origin}$;     do begin /* walk through the voxels         voxel = findnode(Q);         for each object in the voxel's list do             intersect(ray, object);         if no intersection has been found then             Q= point in next voxel pierced by ray;     until intersection found or Q outside world $\}$

• There are potential pitfalls that must be avoided when implementing spatial subdivision
• Two of them arise from the fact that an object may intersect several voxels
• The first is repeated ray-object intersection tests between the same ray and object
• As the ray passes from voxel 1 to 2, it finds object A in object list of each
• One way to avoid repeated intersection tests is to use a mailbox
  • Each object has a mailbox; each ray is tagged by a unique id
  • When an object is tested for intersection, the results of the test and the ray tag are stored in the object's mailbox
  • Before testing every object, the tag stored in the mailbox is compared against that of the current ray, if they match the object has been previously tested and the intersection result can be retrieved without recalculation
  • After testing object A in voxel 1, it can be tested in voxel 2 using the mailbox
• The second pitfall is more dangerous -- it can lead to wrong results
• Notice object B is on the object list for voxel 3 and B is intersected by the ray
• If the affirmative intersection test with object B causes the ray walking to stop at voxel 3, the closer object will not be found

Uniform Spatial Subdivision

• Fujimoto et al, introduced the idea of subdividing space into uniformly sized voxels
• They call this organization of space SEADS (Spatially Enumerated Auxiliary Data Structure)
• The subdivision is independent of the structure of the environment
• The uniform subdivision is a disadvantage in that it take more space and uses more time subdividing areas that may be simple
• However, rays can be traversed through the voxels very efficiently by a 3D-DDA calculation
• This speed of voxel space traversal can offset the disadvantage of the uniform space subdivision
• A point is indexed into a node directly, e.g., assume a $1024\times1024\times1024$ grid, then point $(x,\,y,\,z)$ belongs to node $(\lfloor x \rfloor,\,\lfloor y \rfloor,\,\lfloor z \rfloor)$
• Here all voxels pierced by the ray must be identified (not as in the DDA where only the ``nearest'' pixel is found)
• First, we consider the extension of the 2D-DDA that identifies all pixels pierced by a line, not just those closest to the ideal line

Extending the DDA Algorithm

• Let $R_o+t\vec{R}_d = \langle x_o+tx_d,\,y_o+ty_d\rangle$ denote the ray
• Let s₁=y_d/x_d denote the slope of the ray, and for simplicity assume $0\leq s_1 \leq 1$
• The pixels are identified by their lower left hand corners, and the ray starts in pixel $(\lfloor x_o,\,y_o\rfloor)$
• Depending on the position of the ray's origin and the slope s₁ any of the right, diagonal, or up pixel can be pierced next by the ray
• Let e denote the ``error'' in y at the left hand edge of the current pixel, that is when $x_o+tx_d=\lfloor x_o \rfloor$, or $t=(\lfloor x_o \rfloor - x_o)/x_d$, thus
  
  $$\begin{eqnarray*}e & = & y_o+ ty_d \\ & = & y_o + \frac{y_d}{x_d}(\lfloor x_o \rfloor - x_o) = y_o + s_1(\lfloor x_o \rfloor - x_o) \end{eqnarray*}$$
  
  

• If $e+s_1 < \lfloor y_o\rfloor + 1$, the right pixel is pierced next
• If $e+s_1 = \lfloor y_o\rfloor + 1$, the diagonal pixel is pierced next
• If $e+s_1 > \lfloor y_o\rfloor + 1$, the up pixel is pierced next
• Every time through the loop of the 2D-DDA algorithm, the right or diagonal pixel will be identified, but a special test must be made for the up pixel

`;=3000 2d-dda-extended $(\mbox{\it point}\:R_o,\,\mbox{\it ray}\:R_d)$$\{$     $x = \lfloor x_o\rfloor$;      $y = \lfloor y_o\rfloor$;      s₁ = y_d/x_d;      e = y_o+s₁*(x-x_o);     do $\:\{$        pierced$(x,\,y)$;         e= e+s₁;         if $(e \geq y+1.0)$ then             y=y+1;             if (e > y) pierced$(x,\,y)$;         x = x+1.0;     $\}$ until ray leaves scene; $\}$

3D-DDA Algorithm for Uniform Octree Traversal

• The 3D-DDA used two synchronized 2D-DDA's working in mutually perpendicular planes
• Let s₁=y_d/x_d and s₂=z_d/x_d and suppose both these slopes are between -1 and 1
• Consider the 2D-DDA extended so every pixel hit by a ray is identified
• Assume $R_o=(x_o,\,y_o,\,z_o)$ is the ray's origin and $\vec{R}_d=\langle x_d,\,y_d,\,z_d\rangle$ is its direction
• Let s₁ = y_d/x_d denote the xy slope and s₂=z_d/x_d the xz slope, assume both are between 0 and 1
• A 2D-DDA is run in the xy plane and xz plane simultaneously; together they identify each voxel pierced

Binary Space Partitioning

• We've seen BSP trees as a data structure that can be correctly walked from any view to determine visible surfaces
• Here they are used to partition space into regions where ray will be traced
• BSP trees use a different indexing scheme into spatial regions than octrees
• BSP trees nodes are constructed with explicit pointers to their two children; this increases storage but leads to faster traversal than nonuniform octrees
• In the standard implementation for ray tracing, space is partitioned by planes parallel to the x, y and z axes

Summary of Space Partitioning Methods

• Bounding volumes have exactly one path to each object (they form a tree)
• Nonuniform spatial subdivision (octrees) results in a DAG; there is one path to any leaf volume (circles in the graphs), but leaf volumes can contain many objects and some objects (squares in the graphs) can belong to many leaf volumes
  • Octrees are good for scenes whose occupancy density varies
  • Traversal of octrees tends to be slower than BSP trees and SEADS because the trees tend to be unbalanced
  • BSP trees tend to have smaller depth than octrees because the tree is balanced
• Uniform spatial subdivision (SEADS) results in a bipartite graph
  • Volumes are indexed by direct index calculation, not a path through other volumes
  • SEADS has fast traversal but enormous memory costs

Generalized Rays

• A generalized ray considers an entire family of rays bundled as beams, cones, or pencils
• Some sacrifice is required to use each of these types of generalized ray
  • The types of primitive objects may need to be restricted
  • The computation of ``exact'' intersections may need to be abandoned
• The use of generalized rays leads to advantages such as
  • Increased efficiency by exploiting coherence
  • Effective antialiasing
  • Additional optical effects

Cone Tracing

• A right circular cone of rays can be defined by an apex, center line, and spread angle
• The primitive objects are restricted to spheres, planes, and polygons
• A cone intersects a sphere if
  
  $$D\tan\theta +\frac{R}{\cos\theta} > d$$
  
  
  where
  • $\theta$ is the spread angle of the cone
  • R is the radius of the sphere
  • d is the distance between the center of the sphere and the closest approach of the cone's axis
  • D is the distance from the cone origin to the closest approach on the cone axis
• For antialiasing, the intersection calculation needs to detect how much of the cone is blocked by the object
• For reflection and refraction, the new center line is computed as the reflected or refracted center line of the incident cone

Beam Tracing

• Heckbert and Hanrahan exploit coherence by observing that neighbors of a particular ray tend to follow the same path
• They do this by recursively applying a version of the Weiler-Atherton hidden surface removal algorithm
• The basic steps in the Weiler-Atherton algorithm are
  • A preliminary rough depth sort of polygons based on the nearest z value
  • The first polygon in the sorted list is used to clip the remainder of the list into new lists of inside and outside polygons
  • Removal of any polygons on the inside list behind the current (front) polygon
  • If any polygons remain on the inside list, a recursive call to subdivide the region again with the front polygon and the inside list
• Only polygons can be handled by this algorithm

• The initial beam is the viewing frustrum
• An beam tree is built from intersections of this beam with the environment
• A beam may intersect many surfaces, so each node in the list contains a list of surfaces intersected by the beam
• New beams are generated from beam-object intersections by clipping the incident beam and defining a new virtual eye point
• ``Pencils'' of rays have also be suggested
• A pencil of rays is a collection parallel to and near an axial ray

Ray Coherence

• A similar approach to cone, beam and pencil tracing exploits ray coherence without introducing new geometrical entities (cones, beams, pencils)
• The idea is to use the path (and ray tree) of the previous ray to construct the ray tree for the current ray
• As the current tree is constructed, information from the corresponding branch of the previous tree is use to predict then next object hit
• New intervening objects must be detected

• A ``safety zone'' is constructed around each ray
• If the current ray does not exit the safety zone of the previous ray and intersects the same object as the previous ray, then it can not hit any intervening object
• Thus, two basic tests are made: (1) does the ray hit the same object as the previous ray; (2) does the ray leave the previous ray's safety zone
• Alternatively, we can identify the potential blockers (in a ``cache'')
• A cache miss occurs when the current ray misses the previously hit object or hits a potential blocker

Distributed Ray Tracing

• Cook, Porter, and Carpenter introduced this concept in 1984
• Distributed ray tracing is a Monte Carlo technique that stochastically distributes the direction of rays
• It supersamples the image, firing more than one type
  • Sampling the reflected ray according to a specular distribution function produces gloss (blurred reflection)
  • Sampling the transmitted ray produces translucency (blurred transparency)
  • Sampling the solid angle of the light source produces penumbras (soft shadows)
  • Sampling the camera lens area produces depth of field (focusing)
  • Sampling in time produces motion blur

• Sharp reflections, refraction, shadows ... can be considered aliasing artifacts and distributed ray tracing is an antialiasing approach
• Suppose the 16 eye rays are fired through a pixel; consider the pixel is composed of 16 subpixels
• A sample point is placed in the middle of each subpixel and then noise is added to the x and y locations independently
• Color intensities returned from each sample can be simply averaged (box filtered), or weighted averages can be used

• The goal in shading (not just in ray tracing) is to evaluate the intensity I of the reflected light at a point on the surface
• The intensity is an integral over the sphere around the surface of an illumination function L and a reflection function R
  
  $$I(\phi_r,\,\theta_r) = \int_{\phi_i}\int_{\theta_i}L(\phi_i,... ...) R(\phi_i,\,\theta_i,\,\phi_r,\,\theta_r)\,d\phi_i\,d\theta_i$$
  
  

• The reflection function R describes the surface and the illumination function L describes the light in the environment
• Distributed ray tracing uses Monte Carlo integration techniques to approximate this integral

• Spatial dimension are divided up into $4 \times 4$ grids, the center of each subpixel is jittered
• To jitter in a nonspatial dimension, randomly created ``prototype'' patterns in screen space are associated with sample points within a certain range
• The exact location is then determined by jittering
• For example, sampling in time to produce motion blur, the frame time is divided into slices and a randomly assigned slice of time is assigned to each sample point
• To assign times in a pixel with a $4 \times 4$ grid of sample points, the integers 1 to 16 can be randomly distributed

• The sample in the ith row and j column would yield a prototype time
  
  $$t_{ij}= \frac{P_{ij}-0.5}{16.0}$$
  
  
  where P_ij is the value in the prototype pattern
• A random jitter $\pm 1/32$ is added to t_ij to obtain the actual time, e.g. the sample in the upper left subpixel would have a time $6/16 \leq t \leq 7/16$

• Sometimes the samples need to be weighted; for example, we may want to weight the reflected samples according to a specular distribution function
• Importance sampling is used to distribute the samples; the sample points are distributed so that chance of a location being sampled is proportional to the filter at that location
• The filter is divided into regions of equal area, each region corresponds to one sample point at the center of the region

• Usually the filter is known ahead of time, so the centers and jitter magnitudes can be precomputed and stored in a lookup table
• For example, for the reflection ray each surface has an associated lookup table for a presampled reflection function $R(\phi_i,\,\theta_i,\,\phi_r,\,\theta_r)$
• Given the angle between the surface normal and the incident ray, the lookup table gives the range of reflection angles plus a jitter magnitude
• For depth of field, the focal point is determined from the center of the pixel
• A lens is ``placed at the pixel'' and a point on the lens found from a jittering lens location lookup table
• The ray is then traced to the first hit
• The effect of depth of field is controlled by the diameter of the lens, F/n where F is the focal length and n is the f-stop

Summary of Distributed Ray Tracing

• Determine the spatial screen location of the ray by jittering
• Determine the time for the ray and move the objects and ``camera'' accordingly
• Determine the focal point by following the ray from the center of the lens through the screen location to the distance of the focal distance
• Determine the visible point for this ray using ray-object intersection techniques
• Trace a reflection ray whose direction is determined by jittering a set of directions distributed according to a specular reflection function (use a lookup table directions)
• Trace a transmitted ray if the object is transparent; the direction is determined by jittering a set of directions distributed according to a specular transmission function (use a lookup table directions)
• Trace the shadow rays; the location of the target on the light is determined by its distribution function

Adaptive Supersampling

• Whitted use an adaptive supersampling algorithm to antialiasing ray traced images
• The technique fires rays at pixel corners as well as their center
• If all five rays return about the same color, they can be averaged to produce the pixel's color
• If they return different colors, we'll subdivide the pixel into smaller regions and treat the subpixel just as we did the entire pixel
• The number of subdivisions can be stopped at some maximum level
• This technique is easy, not to slow, and often works well to antialias a ray traced image
• Object too small may still be missed by all rays
• If objects are arranged in a regular pattern and some are missing this causes visible effects
• In animation, objects may appear and disappear as they are hit and missed by rays

Backward Ray Tracing

• Backward ray tracing follows rays from the light source instead of following rays from the eye
• Since most rays from the light will never reach the viewer, this can be wasteful, but it can also lead to more realistic images by including diffuse effects