.NH 1
Interrogating the Model:  Ray-tracing
.NH 2
What Applications Want and Need
.PP
The objective of a given application will to a large degree determine the
most ``natural'' form in which the model might be presented.
For example, extracting just the \fIedges\fR of the objects in a model
would be suitable for a program attempting to construct a wire-frame display
of the model.
Another family of applications exists which needs to be able to find
the intersection of small object paths (eg, photons) with the
model.
Generally, these alternatives are motivated by the representation of
a physical process being simulated, and each alternative
is useful for a whole family of applications.
By choosing the ray sampling density
within the Nyquist limit, these applications are well satisfied by
extracting ray/geometry intersection information,
the well known ``ray-tracing'' algorithm.
However, a mathematical ray has as its cross section a point,
while physical objects have significant cross-sectional area.
This lack of cross-sectional area can result in sampling inaccuracies, and
has lead some applications developers to try to
minimize these difficulties.  Applications which simulate
particles or small rocks approaching the model might benefit
from cylinder/geometry intersection
capability, and applications which shine beams of light on the model
such as spotlights or
even highly collimated light such as laser light might benefit from
cone/geometry intersection capabilities.
.[
Amanatides ray tracing cones
.]
.[
Kirk simulation natural features using cone tracing
.]
Applications which are attempting to simulate wave effects might
be well expressed in terms of plane/geometry intersection curves,
and structural analysis routines would probably prefer to obtain
the geometry as a collection of connected hyperpatches.
.PP
While very recent research has begun to explore techniques for
intersecting cylinders, cones, and planes with geometry,
.[
Kajiya ray tracing procedurally defined objects
.]
ray-tracing is by far the most well developed approach.
Fortunately, most applications can function well with approximate,
sampled data about the model, rather than needing a precise
representation.
Data with statistical validity
can be obtained by sampling the model with
an adequate number of rays and computing the ray/geometry intersections.
This approach is also by far the easiest one to implement,
as the one-dimensional nature of a mathematical ray makes the
intersection equations relatively straightforward.
.PP
Given that one has decided to force all applications to interrogate
the model strictly by computing ray/geometry intersections,
there are still several implementation strategies that could be
adopted.  The most traditional approach has been batch-oriented,
with the user defining a set of ``viewing angles'',
turning loose a big batch job to compute all the ray intersections,
and then post-processing all the ray data into some meaningful form.
However, the major drawback of this approach is that the application
has no immediate control over ray paths, making another batch run
necessary for each level of reflection, etc.
.PP
In order to be successful, applications need: (1) interactive control of ray
paths, to naturally implement reflection, refraction, and fragmenting into
multiple subsidiary rays, and (2) the ability to fire rays in arbitrary
directions from arbitrary points.
Nearly all non-batch implementations have closely
coupled a specific application (typically a model of illumination) with the
ray-tracing code, allowing efficient and effective control of the ray paths.
The most flexible approach of all is to implement the ray-tracing
capability as a general-purpose library, and make the functionality
available as needed to any application.
.NH 2
References
.PP
The origins of modern ray-tracing come from work at MAGI
under contract to BRL, initiated in the early 1960s.
The initial results were reported by MAGI
.[
geometric description technique MAGI
.]
in 1967.
Extensions to the early developments were undertaken by a
DoD Joint Technical Coordinating Group effort, resulting in
publications in 1970
.[
MAGIC User Manual
.]
and 1971.
.[
MAGIC Analyst Manual
.]
The detailed presentation of the
fundamental analysis and implementations of the ray-tracing algorithm
found in these two documents forms an excellent and thorough
review of the principles of solid modeling.
.PP
More recently, interest in ray-tracing developed in the academic
community, with Kay's
.[
Kay Ray Tracing 1979
.]
thesis in 1979
being a notable early work.
One of the central papers in the ray-tracing literature is
the work of Whitted.
.[
Whitted Improved Illumination Model
.]
Model sampling techniques can be improved to provide substantially
more realistic images by using the ``Distributed Ray Tracing'' strategy.
.[
Cook Carpenter Distributed Ray Tracing
.]
For an excellent, concise discussion of ray-tracing,
consult pages 363-381 of Rogers.
.[
Rogers Procedural Elements
.]
.NH 2
Ray -vs- Geometry
.PP
The process of intersecting one ray with the model geometry
can be decomposed into two distinct tasks.
First, the ray must be intersected with each and every primitive
solid in the model.  Each successful intersection will result in
a set of line segments being returned.
Secondly, the boolean combinations expressed in the
model's directed acyclic graph must be evaluated over this
collection of line segments.
.PP
We will examine the details of intersecting a
ray with a halfspace, a generalized ARB, and an ellipsoid.
These three analyses provide all the information necessary to
enable one to intersect a ray with any solid possessing
a quadratic surface.  For example, a truncated general cone
can be dealt with as an infinite cylinder (a quadratic)
bounded by two planes.  The extension to surfaces of higher
order is straightforward, and depends only on having a
reliable root-finding algorithm.
Information about ray-tracing of spline objects can be found in a recent
paper by Sweeny.
.[
Sweeny Ray Tracing
.]
.PP
Rays begin at a point $bold P vec$,
and proceed infinitely in a given direction $bold D vec$.
Any point $bold A vec$ on a ray may be expressed as a linear
combination of $bold P vec$ and $bold D vec$:
.EQ C
bold A vec ~=~ bold P vec ~+~ k ~*~ bold D vec
.EN
where valid solutions for $k$ are in the range $[ 0, inf )$.
$( bold D sub x , bold D sub y , bold D sub z )$ are the
\fIdirection cosines\fR for the ray, being the cosine of
the angle between the ray and the appropriate axis.
While not necessary, in this analysis it is assumed that
$bold D vec$ is of unit length, i.e. $| bold D vec | = 1$.
.NH 2
Ray/Halfspace Intersection
.PP
Consider a plane with arbitrary orientation.  This plane
partitions space into two half-spaces.  Define 
$bold N vec$ to be an outward-pointing unit-length vector
which is normal to the plane.
Let point $bold A vec$ lie anywhere on the plane.
Let the ray in question be defined such that any point $bold X vec$
on the ray may be expressed as
$bold X vec = bold P vec + k bold D vec$.
.PP
Determine if the ray is parallel to the surface of the plane, and
if not, determine if the ray is entering or exiting the half-space
by comparing the direction of the ray with the direction of the
normal $bold N vec$:
.EQ
bold N vec cdot bold D vec ~=~ left {
lpile { 0 above >0 above <0 }
~~ roman lpile {
parallel,~no~solution above
exiting~half-space,~segment~t=[ 0 , k ] above
entering~half-space,~segment~t=[ k , inf ) }
.EN
.PP
If a solution exists, i.e. $bold N vec cdot bold D vec != 0$,
find parameter $k$ for intersection of ray and plane.
To accomplish this, a useful form for the equation of the plane is
$bold N vec cdot ( bold X vec - bold A vec ) = 0$,
where the vector from any point on the plane ($bold A vec$)
is tested for perpendicularity with the normal $bold N vec$.
To find $k$, we substitute the parametric version of $bold X vec$, giving:
.EQ
bold N vec cdot ( bold P vec + k * bold D vec - bold A vec ) = 0
.EN
.EQ
bold N vec cdot bold P vec + k ( bold N vec cdot bold D vec ) -
bold N vec cdot bold A vec = 0
.EN
.br
.EQ
k = { bold N vec cdot bold A vec - bold N vec cdot bold P vec }
over { bold N vec cdot bold D vec }
.EN
.KF
.PS
scale=300
box invis ht 201 wid 526 with .sw at 0,0
circle rad 8 at 496,41
circle rad 8 at 88,65
line -> from 88,65 to 24,129 
line  from 16,41 to 24,25 
line  from 40,49 to 48,33 
line  from 64,57 to 72,41 
line  from 112,73 to 120,57 
line  from 136,81 to 144,65 
line  from 160,89 to 168,73 
line  from 184,97 to 192,81 
line  from 208,105 to 216,89 
line  from 232,113 to 240,97 
line  from 280,129 to 288,113 
line  from 304,137 to 312,121 
line  from 328,145 to 336,129 
line  from 352,153 to 360,137 
line  from 16,41 to 352,153 
line -> from 496,41 to 16,201 
spline  from 256,121\
to 248,82\
to 313,82\
to 341,79\
to 352,65\
to 352,41
"\f1\s16\&N\f1\s0" at 8,133
"\f1\s16\&A\f1\s0" at 112,29
"\f1\s16\&k\f1\s0" at 344,13
"\f1\s16\&P\f1\s0" at 520,61
"\f1\s16\&D\f1\s0" at 400,93
spline  from 496,41\
to 474,16\
to 450,16\
to 400,57\
to 376,57\
to 352,41
circle rad 8 at 256,121
.PE
.ce 1
\fBFigure 3.1 \(em Ray/Halfspace Intersection\fR
.KE
.NH 2
Ray/ARB Intersection
.PP
An ARB is a convex volume bounded by
4 (pyramid), 5 (wedge), or 6 (box) planes.
This analysis depends on the properties of objects with convex hulls.
Let the ray in question be defined such that any point $bold X vec$
on the ray may be expressed as
$bold X vec = bold P vec + k bold D vec$.
Intersect the ray with each of the
planes bounding the ARB as discussed above,
and record the values of the
parametric distance $k$ along the ray.
With outward pointing normal vectors, note that the
ray \fIenters\fR the half-space defined by a plane when
$bold D vec cdot bold N vec < 0$,
is \fIparallel\fR to the plane when $bold D vec cdot bold N vec = 0$,
and \fIexits\fR otherwise.
Find the \fIentry\fR point farthest away from the starting point $bold P vec$,
i.e. it has the largest value of $k$ among the entry points.
The ray enters the solid at this point.
Similarly, find the \fIexit\fR point closest to point
$bold P vec$, i.e. it has the smallest value of $k$ among the exit points.
The ray exits the solid here.
.KF
.PS
scale=300
box invis ht 319 wid 574 with .sw at 0,0
line  from 261,18 to 436,311 
line  from 384,226 to 376,232 
line  from 375,210 to 367,216 
line  from 367,195 to 359,200 
line  from 356,179 to 348,185 
line  from 346,163 to 338,169 
line  from 337,148 to 329,154 
line  from 327,131 to 319,137 
line  from 318,115 to 310,121 
line -> from 355,17 to 29,302 
"\f1\s16\&P\f1\s0" at 383,13
"\f1\s16\&D\f1\s0" at 337,66
spline  from 262,99\
to 268,121\
to 257,134\
to 232,138\
to 223,144\
to 219,161\
to 224,168
line  from 192,92 to 189,102 
line  from 208,94 to 205,104 
line  from 224,96 to 221,106 
line  from 240,97 to 237,107 
line  from 303,103 to 300,113 
line  from 287,101 to 284,111 
line  from 383,249 to 383,241 
line  from 367,249 to 367,241 
line  from 351,249 to 351,241 
line  from 335,249 to 335,241 
line  from 319,249 to 319,241 
line  from 303,249 to 303,241 
line  from 287,249 to 287,241 
line  from 272,249 to 272,241 
line  from 256,249 to 256,241 
line  from 240,249 to 240,241 
line  from 224,249 to 224,241 
line  from 208,249 to 208,241 
line  from 192,249 to 192,241 
circle rad 8 at 291,71
circle rad 8 at 262,99
circle rad 8 at 175,174
circle rad 8 at 89,249
circle rad 8 at 355,17
line  from 0,249 to 534,249 
line  from 8,74 to 574,130 
spline  from 175,174\
to 189,185\
to 202,176\
to 202,167\
to 211,154\
to 218,157\
to 226,160
line  from 177,234 to 188,234 
line  from 177,219 to 188,219 
line  from 176,201 to 187,201 
line  from 175,187 to 186,187 
line  from 174,155 to 185,155 
line  from 173,139 to 184,139 
line  from 173,123 to 184,123 
line  from 172,107 to 183,107 
line  from 178,319 to 168,11 
.PE
.ce 1
\fBFigure 3.2 \(em Ray/ARB Intersection\fR
.KE
.NH 2
Ray/Ellipse Intersection
.PP
Let the point $W$ belong to the set of points
on the surface of an ellipsoid which is defined by
$bold V vec , bold A vec , bold B vec , bold C vec$, where $bold V vec$
is the center of the ellipsoid, and $bold A vec$, $bold B vec$,
and $bold C vec$ are mutually perpendicular and define the orientation
and shape of the ellipsoid.
By affine transformations, map the points on the
original ellipsoid into the set of points $W hat$ in the ``hat'' space
which lie on the surface of a unit sphere, located at the origin.
For any point $X$, the corresponding point in the ``hat'' space is
$X hat = S(R( X - V ))$,
where
.EQ
R( X ) = left [ pile {
{ A vec / ( | A vec |) }
 above { B vec / ( | B vec | ) }
 above { C vec / ( | C vec | ) }
} right ] ~*~ X
.EN
and
.EQ
S(X) =	 left [ matrix {
col { 1 / | A vec | above 0 above 0 }
col { 0 above 1 / | B vec | above  0 }
col { 0 above 0 above  1 / | C vec | }
} right ] ~*~ X
.EN
Intersect this ellipse with a ray;
points on the ray $L$ are of the form $P + k * D$.
This can be be expressed as the intersection of
the ray in ``hat'' space with the unit sphere, where
points on the ray $L hat$ are of the form $P hat + k hat * D hat$.
The mapping between regular and ``hat'' space for vectors is:
.EQ
D hat = S( R( D ) )
.EN
and for points is:
.EQ
P hat = S( R( P - V ) )
.EN
With $W$ defined to be the point where ray $L$ intersects the ellipsoid,
and $W hat$ is the point where $L hat$ intersects the unit sphere, then
.EQ
W hat = P hat + k hat * D hat
.EN
and
.EQ
W hat = S ( R ( W - V ) )
.EN
.LP
In order to find the point $W hat$, we refer to the definition for
a sphere with unit radius located at the origin
.EQ
x sup 2 + y sup 2 + z sup 2 ~=~ 1
.EN
and substitute in the parametric definition of $W hat$, giving
.EQ
( P hat sub x + k hat * D hat sub x ) sup 2 ~+~
( P hat sub y + k hat * D hat sub y ) sup 2 ~+~
( P hat sub z + k hat * D hat sub z ) sup 2 ~=~ 1
.EN
Regrouping terms as a polynomial in $k hat$ gives
.EQ
( D hat sub x sup 2 + D hat sub y sup 2 + D hat sub z sup 2 )
 k hat sup 2 ~+~
2 * ( D hat sub x P hat sub x +
     D hat sub y P hat sub y +
     D hat sub z P hat sub z ) k hat
~+~ P hat sub x sup 2 + P hat sub y sup 2 + P hat sub z sup 2 - 1 ~=~ 0
.EN
.br
.EQ
( D hat cdot D hat ) k hat sup 2 ~+~ 2 ( D hat cdot P hat ) k hat ~+~
[ ( P hat cdot P hat ) - 1 ] ~=~ 0
.EN
Solving the quadratic equation for $k hat$ yields 0, 1, or 2 roots
.EQ
k hat = {
  -2 ( D hat cdot P hat ) +- sqrt {
    [ 2 ( D hat cdot P hat ) ] sup 2
      - 4 ( D hat cdot D hat ) [ ( P hat cdot P hat ) - 1 ]
    }
} over {
    2 ( D hat cdot D hat )
}
.EN
.br
.EQ
k hat = {
     - ( D hat cdot P hat ) +- sqrt {
        ( D hat cdot P hat ) sup 2
           - ( D hat cdot D hat ) [ ( P hat cdot P hat ) - 1 ]
    }
} over {
    ( D hat cdot D hat )
}
.EN
.LP
If a solution for $k hat$ exists, these values of $k hat$
define the intersection point(s) $W hat$.  In order to determine
$W$, the relationship between $k hat$ and $k$ needs to be found.
.EQ
W hat mark = P hat + k hat * D hat
.EN
.br
.EQ
W hat lineup = S ( R ( W - V ) )
.EN
can be rearranged in terms of $W$
.EQ
W lineup = V + R sup -1 ( S sup -1 ( W hat ) ) 
.EN
Then, substituting in the parametric form of $W hat$
.EQ
lineup = V + R sup -1 ( S sup -1 [ P hat + k hat * D hat ] ) 
.EN
.br
.EQ
lineup = V + R sup -1 ( S sup -1 [ S( R( P - V ) ) + k hat * S( R( D ) ) ] )
.EN
.br
.EQ
lineup = V + R sup -1 ( R( P - V ) + k hat * R( D ) )
.EN
.br
.EQ
lineup = V + ( P - V ) + k hat * D
.EN
.br
.EQ
W lineup = P + k hat * D
.EN
but,
.EQ
W lineup = P + k * D
.EN
so $k = k hat$.
This useful result is not surprising because
affine transformations
applied to straight lines yield straight lines, and
$R$ and $S$ are both affine transformation matrices.
The fact that $k$ is constant for both $W$ and $W hat$ means that
the solution for $k hat$ can be determined from $D hat$ and $P hat$
without needing to transform the results out of the ``hat'' space.
.LP
Given that we can find the point $W hat$
where ray $L hat$ intersects the unit sphere and
the point $W$ where ray $L$ intersects the ellipsoid, what is the vector
normal to the tangent plane at that point?
The tangent plane on the unit sphere at $W hat$ has
normal vector $N hat ~=~ W hat$.
(The vector from the origin to any point on the surface of the sphere
is pointing in the direction of the tangent at that point).
.EQ
N hat mark = W hat = S ( R ( W - V ) )
.EN
The tangent plane at $W hat$ transforms back to the tangent plane at $W$,
with normal vector $N vec$:
.EQ
N vec lineup = { [ { ( R sup -1  *  S sup -1 ) } sup T ] } sup -1 ( N hat )
.EN
because if $H vec$ is perpendicular to plane $Q$, and matrix $M$ maps
from $Q$ to $Q hat$,
then ${ [ { M } sup T ] } sup -1 ( H vec )$ is perpendicular to $Q hat$.
.EQ
N vec lineup = { [ { ( S sup -1 ) } sup T * { ( R sup -1 ) } sup T ] } sup -1
  ( N hat )
.EN
Because $( R sup -1 ) sup T = R$ and $S sup T = S$,
.EQ
N vec lineup = { [ S sup -1 * R ] } sup -1 ( N hat )
.EN
.EQ
lineup = R sup -1  * S ( N hat )
.EN
.EQ
lineup =  R sup -1 ( S( S( R( W - V ) ) ) )
.EN
Note that $| N vec | != 1$.
.NH 2
Boolean Evaluation
.PP
In the context of the three dimensional model, each ray passes
through zero or more objects, but the objects are
defined in terms of a boolean combination of primitive solids.
The first part of the task is to intersect the ray with
each and every primitive solid in the model.
For each solid,
either the ray will miss, in which case nothing further need be done,
or the ray will intersect the solid in one or more line segments.
When all of the primitive solids have been intersected, the
segments have to be sorted, and ``woven'' together into
intervals.  The original ray is broken into a new interval
at each point where a primitive solid was entered or exited;
the list of primitive solids intersected is constant within each interval.
Note that the correct entry and exit normals must be selected
each time a segment is broken into an additional interval.
.PP
After all the segments are woven together into intervals, each
interval needs to be applied to the boolean expression tree.
For all primitive solids in the expression tree,
if a segment is present within an interval,
the expression tree sees a TRUE input, otherwise a FALSE input is presented.
For all intervals in which the expression tree returns a TRUE output,
the ray has passed through an actual solid object, and that interval
is added to the solution list for this ray/model intersection.
Note that by keeping the interval list sorted front-to-back
in increasing $k$, and by
incrementally weaving new segments into the interval list as they
are computed, applications such as lighting models which are
only interested in the first intersection with a solid object
can be processed much more efficiently by stopping the intersection
process after the first hit.
.KF
.PS
scale=300
box invis ht 700 wid 719 with .sw at 0,0
box ht 8 wid 248 with .nw at 371,82 
box ht 8 wid 304 with .nw at 315,114 
line  from 467,536 to 467,0 
line  from 619,538 to 619,2 
box ht 89 wid 200 with .nw at 267,407 
box ht 2 wid 697 with .nw at 1,370 
line -> from 697,369 to 719,369 
ellipse ht 162 wid 300 at 467,363
circle rad 184 at 187,376
line  from 3,536 to 3,0 
line  from 267,536 to 267,0 
line  from 315,536 to 315,0 
line  from 371,536 to 371,0 
box ht 8 wid 368 with .nw at 3,176 
box ht 8 wid 200 with .nw at 267,144 
.PE
.ce 1
\fBFigure 3.3 \(em Boolean Evaluation\fR
.KE
.NH 2
Bounding Volumes & Space Partitioning
.PP
Part of the computational expense of ray-tracing derives from the
statement that ``the first part of the task is to intersect the ray with
each and every primitive solid in the model'' results in a
an enormous amount of computation.
Much of the recent research in ray-tracing techniques has focused
on strategies for reducing the the amount of computation required.
Overall, the goal is to develop methods to avoid computing the
actual ray/solid intersections when the ray does not pass near the solid.
.PP
The first method for eliminating unnecessary ray/solid intersections
involves enclosing the solid within a \fIbounding solid\fR.
Before computing the intersection of a ray with the actual primitive
solid, the ray is first intersected with the bounding solid
to reject obvious misses.  This is only effective if the cost of
computing the intersection with the bounding solid is significantly
less than the cost of intersecting the ray with the primitive solid.
Spheres, axis-aligned right-parallelpipeds (RPPs), and boxes
are the three most popularly used bounding volumes.
An important issue is the measure of the quality of
a bounding volume.  The desire is to minimize
$VOL sub bound ~-~ VOL sub geom$, but
each bounding volume is a poor fit for some objects.
The sphere is a poor bound for any object which is long and thin.
Axis-aligned RPPs are poor bounds for any object where the
greatest width is not axis aligned.  The torus represents
a solid that is very hard to bound well at all.
.PP
The best bounding sphere has it's center located at the centroid
of the primitive solid, with the radius $r$ set to the smallest
possible value that encloses the entire primitive solid.
Then, to determine if a ray needs to be intersected with the
primitive solid, it is merely necessary to determine if
the distance from the sphere center $bold V vec$
to line $bold P vec ~+~ k * bold D vec ~<= r$, i.e.
.EQ
left | {
( bold V vec - bold P vec ) times
{ bold D vec } over { | bold D vec | }
} right |
~<=~ r
.EN
When $| bold D vec | = 1$, this can be quite fast to compute.
Another way to speed the calculation of this would be to
compare the magnitude squared to the radius squared,
eliminating the necessity for taking a square root in finding
the distance.
.KF
.PS
scale=300
box invis ht 458 wid 614 with .sw at 0,0
spline  from 303,255\
to 303,287\
to 271,311\
to 223,287\
to 196,321
"\f1\s16\&D\f1\s0" at 174,20
"\f1\s16\&P\f1\s0" at 6,20
"\f1\s16\&dist\f1\s0" at 366,204
"\f1\s16\&V\f1\s0" at 334,292
"\f1\s16\&r\f1\s0" at 198,340
circle rad 8 at 38,8
line  from 302,256 to 100,256 
circle rad 8 at 302,256
line -> from 38,8 to 614,200 
line  from 366,144 to 374,120 
line  from 342,136 to 366,144 
line  from 302,256 to 350,112 
circle rad 202 at 302,256
spline  from 100,256\
to 100,288\
to 132,312\
to 180,288\
to 204,321
.PE
.ce 1
\fBFigure 3.4 \(em Bounding Sphere\fR
.KE
.PP
A bounding RPP
is defined by six numbers giving the extent of the RPP along each of
the three axes, i.e. $MIN sub { x,y,z }$ and $MAX sub { x,y,z }$.
Cyrus and Beck
.[
Cyrus Beck clipping
.]
provide an efficient algorithm for intersecting
the ray with the six planes of the RPP.  By taking advantage of
the fact that all six planes are axis-aligned, the plane intersection
equations are greatly simplified, and the following algorithm
(written in pseudo-\fBC\fR) provides the fastest possible
ray/RPP intersection, with a computational cost very close
to the cost of computing the distance from a bounding sphere.
(Assuming that the division by $D sub i$ is replaced with a multiplication
by the inverse, which is calculated once per ray).
.br
.nf
for( i = x, y, z )  {
    if( $D sub i$ == 0 )  {
        if( $MIN sub i$ > $P sub i$ )  return(MISS);
        if( $MAX sub i$ < $P sub i$ )  return(MISS);
    } else if( $D sub i$ > 0 )  {
        /* Heading toward larger numbers */
        if( $MAX sub i$ < $P sub i$ )  return(MISS);
        k = $(MAX sub i - P sub i ) / D sub i$;
        if( maxK > k )  maxK = k;
        k = $(MIN sub i - P sub i ) / D sub i$;
        if( minK < k )  minK = k;
    } else {
        /* $D sub i$ < 0, heading to smaller numbers */
        if( $MIN sub i > P sub i$ )  return(MISS);
        k = $(MIN sub i - P sub i ) / D sub i$;
        if( maxK > k )  maxK = k;
        k = $(MAX sub i - P sub i ) / D sub i$;
        if( minK < k )  minK = k;
    }
}
if( minK > maxK )  return(MISS);
.fi
If the ray intersects the RPP at all, it does so in the range (minK, maxK).
.PP
Space partitioning is an efficiency strategy which
only intersects rays with bounding volumes that are "near" the path of
the ray, and avoids finding intersections in irrelevant areas of the model.
.[
Glassner space subdivision
.]
.[
Fujimoto accelerated ray tracing
.]
.[
Fujimoto fast elaboration geometry
.]
.PP
Bin-trees
enclose all of model space with 6 axis-aligned planes (an RPP), and
split the RPP with a single plane carefully chosen to reduce complexity
of at least one side of the tree by one solid.
Lists of the contents of both new RPPs are made, and then
both new RPPs are recursed into until the list of contents is ``small''.
When a ray is to be fired,
find the first RPP the ray enters by walking the tree,
making binary decisions at each node.
When a leaf node is found in the space partitioning tree,
fire rays at solids listed (or their bounding volumes).
.PP
Oct-Trees operate similarly, except that at each node space is
uniformly split by planes in X, Y, and Z.
.NH 2
RT Library Interface
.PP
In order to give all applications interactive control over
the ray paths, and to allow the rays to be fired in arbitrary directions
from arbitrary points, BRL has implemented its second generation
ray-tracing capability as a set of library routines.
The RT library exists to allow application programs to
intersect rays with model geometry.  There are two parts to the
interface: preparation routines and the actual ray-tracing routine.
Three "preparation" routines exist;  the first routine which must be called is
dir_build(), which opens the database file, and builds the
in-core database table of contents.
The second routine to be called is get_tree(), which
adds a database sub-tree to the active model space.
get_tree() can be called multiple times
to load different parts of the database
into the active model space.
The third routine is rt_prep(), which
computes the space partitioning data structures, and does other
initialization chores, prior to actual ray-tracing.
Calling this routine is optional,
as it will be called by shootray() if needed.
rt_prep() is provided as a separate routine to
facilitate more accurate timing of the preparation and ray-tracing phases of
applications.
.PP
To compute the intersection of a ray with the geometry in the
active model space, the application must call shootray() once for each
ray.
Ray-path selection for perspective, reflection, refraction, etc,
is entirely determined by the applications program, and passed as a parameter
to shootray() in the RT ``application'' structure, which contains five
major elements:
the vector a_ray.r_pt ($bold P vec$) which is
the starting point of the ray to be fired,
the vector a_ray.r_dir ($bold D vec$) which is
the unit-length direction vector of the ray,
the pointer *a_hit() which is the address of an application-provided
routine to call on those rays where some geometry is hit by the ray,
the pointer *a_miss() which is the address of an application-provided
routine to call on those rays where the ray does not hit any geometry,
the flag a_onehit which is set non-zero to stop ray-tracing
as soon as the ray has intersected at least one piece of geometry
(useful for lighting models),
plus various locations for applications to store state
(recursion level, colors, etc).
Note that
the return from the application provided a_hit()/a_miss() routine
is the formal return of the function shootray().
The shootray() function is prepared for full recursion so that
the application provided a_hit()/a_miss() routines can themselves
fire additional rays by calling shootray() recursively before deciding their
own return value.
In addition, the function shootray() is fully prepared to be operating
in parallel with other instances of itself in the same address space,
allowing the application to take advantage of parallel hardware capabilities
where such exist.  This has actually been successfully demonstrated on
the Denelcor HEP H-1000 machine and on the Alliant FX/8.
.NH 2
Sample RT Application
.PP
A simple application program that fires one ray at a model and prints
the result is included below, to demonstrate the simplicity of the
interface to the RT library.
.LP
.nf
struct application ap;
main() {
      dir_build("model.g");
      get_tree("car");
      rt_prep();
      ap.a_point = [ 100, 0, 0 ];
      ap.a_dir = [ -1, 0, 0 ];
      ap.a_hit = &hit_geom;
      ap.a_miss = &miss_geom;
      ap.a_onehit = 1;
      shootray( &ap );
}
hit_geom(app, part)
struct application *app;
struct partition *part;
{
      printf("Hit %s", part->pt_forw->pt_regionp->reg_name);
}
miss_geom(){
      printf("Missed");
}
.fi
.NH 2
RT Shootray
.PP
To show the simplicity of the overall ray-tracing algorithm, and
to provide the reader with enough information to be able to implement
a full ray-tracing package, the RT shootray() function is
presented.
.LP
.nf
shootray(app)
struct application *app;
{
      if( ray outside model RPP )  return( a_miss() );
      do "push" from box to box  {
            for all solids in box {
                  check bounding RPP
                  invoke SHOOT on solid, collect segment list
            }
            bool_weave():  add segs to interval list
            if( app->a_onehit > 0 )  {
                  bool_final()
                  if( any geometry hit )  break;
            }
      }
      bool_final();
      if( any geometry hit )
            call a_hit();
      else
            call a_miss();
}
.fi
.PP
The overall strategy here is simple, and only the addition of the space
partitioning ``box'' (voxel) notion adds any complexity.  If the ray
does not even intersect the model RPP, then call a_miss().
Otherwise, examine every box (voxel) that the ray pierces, in increasing
distance from the starting point $bold P vec$.  For every solid inside
each voxel, call the SHOOT routine appropriate for that type of solid,
and collect the resultant lists of intersection segments.  After
processing all the solids in the voxel, call bool_weave() to sort
the segments into the interval list, creating additional intervals
where segments overlap.  If the a_onehit flag is set, call bool_final()
to evaluate the boolean expressions on the sorted interval list;  if
any interval evaluates as TRUE, then break the loop and call a_hit().
When all voxels have been examined, call bool_final(), and call a_hit()
or a_miss() as appropriate.
.PP
It is important to note that the BRL MGED and RT software is available
from BRL at no cost
by writing to the author.  In this way, the author hopes to prevent
the implementation of more ray-tracers, and focus attention onto
more fruitful areas of research, such as developing more sophisticated
algorithms for the analysis of solid models.
