## 22.6. 22.6 Rendering with ray tracing

When a virtual world is rendered, we have to identify the surfaces visible in different directions from the virtual eye. The set of possible directions is defined by a rectangle shaped window which is decomposed to a grid corresponding to the pixels of the screen (Figure 22.30). Since a pixel has a unique colour, it is enough to solve the visibility problem in a single point of each pixel, for example, in the points corresponding to pixel centres.

The surface visible at a direction from the eye can be identified by casting a half line, called ray, and identifying its intersection closest to the eye position. This operation is called ray tracing. Ray tracing has many applications. For example, shadow computation tests whether or not a point is occluded from the light source, which requires a ray to be sent from the point at the direction of the light source and the determination whether this ray intersects any surface closer than the light source. Ray tracing is also used by collision detection since a point moving with constant and uniform speed collides that surface which is first intersected by the ray describing the motion of the point.

A ray is defined by the following equation:

where is the place vector of the ray origin, is the direction of the ray, and ray parameter characterizes the distance from the origin. Let us suppose that direction vector has unit length. In this case parameter is the real distance, otherwise it would only be proportional to the distance.

Footnote. In collision detection is not a unit vector, but the velocity of the moving point since this makes ray parameter express the collision time.

If parameter is negative, then the point is behind the eye and is obviously not visible. The identification of the closest intersection with the ray means the determination of the intersection point having the smallest, positive ray parameter. In order to find the closest intersection, the intersection calculation is tried with each surface, and the closest is retained. This algorithm obtaining the first intersection is:

Ray-First-Intersection( )

  1
Initialization to the maximum size in the virtual world.   2
FOR
each object    3
DO

Ray-Surface-Intersection(

)

Negative if no intersection exists.  4
IF

Is the new intersection closer?   5
THEN

Ray parameter of the closest intersection so far.   6
Closest object so far.   7
IF

THEN

Has been intersection at all?   8
THEN

Intersection point using the ray equation.   9
RETURN

10
ELSE

RETURN
“no intersection”       No intersection. 

This algorithm inputs the ray defined by origin and direction , and outputs the ray parameter of the intersection in variable , the intersection point in , and the visible object in . The algorithm calls function Ray-Surface-Intersection for each object, which determines the intersection of the ray and the given object, and indicates with a negative return value if no intersection exists. Function Ray-Surface-Intersection should be implemented separately for each surface type.

### 22.6.1. 22.6.1 Ray surface intersection calculation

The identification of the intersection between a ray and a surface requires the solution of an equation. The intersection point is both on the ray and on the surface, thus it can be obtained by inserting the ray equation into the equation of the surface and solving the resulting equation for the unknown ray parameter.

#### 22.6.1.1.  Intersection calculation for implicit surfaces.

For implicit surfaces of equation , the intersection can be calculated by solving the following scalar equation for : .

Let us take the example of quadrics that include the sphere, the ellipsoid, the cylinder, the cone, the paraboloid, etc. The implicit equation of a general quadric contains a quadratic form:

where is a matrix. Substituting the ray equation into the equation of the surface, we obtain

Rearranging the terms, we get a second order equation for unknown parameter :

where and .

This equation can be solved using the solution formula of second order equations. Now we are interested in only the real and positive roots. If two such roots exist, then the smaller one corresponds to the intersection closer to the origin of the ray.

#### 22.6.1.2.  Intersection calculation for parametric surfaces.

The intersection of parametric surface and the ray is calculated by first solving the following equation for unknown parameters

then checking whether or not is positive and parameters are inside the allowed parameter range of the surface.

Roots of non-linear equations are usually found by numeric methods. On the other hand, the surface can also be approximated by a triangle mesh, which is intersected by the ray. Having obtained the intersection on the coarse mesh, the mesh around this point is refined, and the intersection calculation is repeated with the refined mesh.

#### 22.6.1.3.  Intersection calculation for a triangle.

To compute the ray intersection for a triangle of vertices , , and , first the ray intersection with the plane of the triangle is found. Then it is decided whether or not the intersection point with the plane is inside the triangle. The normal and a place vector of the triangle plane are , and , respectively, thus points of the plane satisfy the following equation:

The intersection of the ray and this plane is obtained by substituting the ray equation (equation (22.24)) into this plane equation, and solving it for unknown parameter . If root is positive, then it is inserted into the ray equation to get the intersection point with the plane. However, if the root is negative, then the intersection is behind the origin of the ray, thus is invalid. Having a valid intersection with the plane of the triangle, we check whether this point is inside the triangle. This is a containment problem, which is discussed in Subsection 22.4.1.

#### 22.6.1.4.  Intersection calculation for an AABB.

The surface of an AABB, that is an axis aligned block, can be subdivided to 6 rectangular faces, or alternatively to 12 triangles, thus its intersection can be solved by the algorithms discussed in the previous subsections. However, realizing that in this special case the three coordinates can be handled separately, we can develop more efficient approaches. In fact, an AABB is the intersection of an -stratum defined by inequality , a -stratum defined by and a -stratum of inequality . For example, the ray parameters of the intersections with the -stratum are:

The smaller of the two parameter values corresponds to the entry at the stratum, while the greater to the exit. Let us denote the ray parameter of the entry by , and the ray parameter of the exit by . The ray is inside the -stratum while the ray parameter is in . Repeating the same calculation for the and -strata as well, three ray parameter intervals are obtained. The intersection of these intervals determine when the ray is inside the AABB. If parameter obtained as the result of intersecting the strata is negative, then the AABB is behind the eye, thus no ray–AABB intersection is possible. If only is negative, then the ray starts at an internal point of the AABB, and the first intersection is at . Finally, if is positive, then the ray enters the AABB from outside at parameter .

The computation of the unnecessary intersection points can be reduced by applying the Cohen – Sutherland line clipping algorithm (subsection 22.4.3). First, the ray is replaced by a line segment where one endpoint is the origin of the ray, and the other endpoint is an arbitrary point on the ray which is farther from the origin than any object of the virtual world. Then this line segment is tried to be clipped by the AABB. If the Cohen – Sutherland algorithm reports that the line segment has no internal part, then the ray has no intersection with the AABB.

### 22.6.2. 22.6.2 Speeding up the intersection calculation

A naive ray tracing algorithm tests each object for a ray to find the closest intersection. If there are objects in the space, the running time of the algorithm is both in the average and in the worst case. The storage requirement is also linear in terms of the number of objects.

The method would be speeded up if we could exclude certain objects from the intersection test without testing them one by one. The reasons of such exclusion include that these objects are “behind” the ray or “not in the direction of the ray”. Additionally, the speed is also expected to improve if we can terminate the search having found an intersection supposing that even if other intersections exist, they are surely farther than the just found intersection point. To make such decisions safely, we need to know the arrangement of objects in the virtual world. This information is gathered during the pre-processing phase. Of course, pre-processing has its own computational cost, which is worth spending if we have to trace a lot of rays.

#### 22.6.2.1.  Bounding volumes.

One of the simplest ray tracing acceleration technique uses bounding volumes, The bounding volume is a shape of simple geometry, typically a sphere or an AABB, which completely contains a complex object. When a ray is traced, first the bounding volume is tried to be intersected. If there is no intersection with the bounding volume, then neither can the contained object be intersected, thus the computation time of the ray intersection with the complex object is saved. The bounding volume should be selected in a way that the ray intersection is computationally cheap, and it is a tight container of the complex object.

The application of bounding volumes does not alter the linear time complexity of the naive ray tracing. However, it can increase the speed by a scalar factor.

On the other hand, bounding volumes can also be organized in a hierarchy putting bounding volumes inside bigger bounding volumes recursively. In this case the ray tracing algorithm traverses this hierarchy, which is possible in sub-linear time.

#### 22.6.2.2.  Space subdivision with uniform grids.

Let us find the AABB of the complete virtual world and subdivide it by an axis aligned uniform grid of cell sizes (Figure 22.31).

In the preprocessing phase, for each cell we identify those objects that are at least partially contained by the cell. The test of an object against a cell can be performed using a clipping algorithm (subsection 22.4.3), or simply checking whether the cell and the AABB of the object overlap.

Uniform-Grid-Construction()

  1  Compute the minimum corner of the AABB ()        and cell sizes   2
FOR
each cell    3
DO
object list of cell empty   4
FOR
each object
Register objects overlapping                     with this cell.  5
DO

IF
cell  and the AABB of object  overlap   6
THEN
add object  to object list of cell


During ray tracing, cells intersected by the ray are visited in the order of their distance from the ray origin. When a cell is processed, only those objects need to be tested for intersection which overlap with this cell, that is, which are registered in this cell. On the other hand, if an intersection is found in the cell, then intersections belonging to other cells cannot be closer to the ray origin than the found intersection. Thus the cell marching can be terminated. Note that when an object registered in a cell is intersected by the ray, we should also check whether the intersection point is also in this cell.

We might meet an object again in other cells. The number of ray–surface intersection can be reduced if the results of ray–surface intersections are stored with the objects and are reused when needed again.

As long as no ray–surface intersection is found, the algorithm traverses those cells which are intersected by the ray. Indices of the first cell are computed from ray origin , minimum corner of the grid, and sizes of the cells:

Uniform-Grid-Enclosing-Cell( )

  1

Integer(

)
2

Integer(

)
3

Integer(

)
4
RETURN



The presented algorithm assumes that the origin of the ray is inside the subspace covered by the grid. Should this condition not be met, then the intersection of the ray and the scene AABB is computed, and the ray origin is moved to this point.

The initial values of ray parameters are computed as the intersection of the ray and the coordinate planes by the Uniform-grid-ray-parameter-initialization algorithm:

Uniform-Grid-Ray-Parameter-Initialization( )

  1
IF

2
THEN

3
ELSE

IF

4
THEN

5
ELSE

The maximum distance.   6
IF

7
THEN

8
ELSE

IF

9
THEN

10
ELSE

11
IF

12
THEN

13
ELSE

IF

14
THEN

15
ELSE

16
RETURN



The next cell of the sequence of the visited cells is determined by the 3D line drawing algorithm (3DDDA algorithm). This algorithm exploits the fact that the ray parameters of the intersection points with planes perpendicular to axis (and similarly to axes and ) are regularly placed at distance (, and , respectively), thus the ray parameter of the next intersection can be obtained with a single addition (Figure 22.31). Ray parameters , , and are stored in global variables, and are incremented by constant values. The smallest from the three ray parameters of the coordinate planes identifies the next intersection with the cell.

The following algorithm computes indices of the next intersected cell, and updates ray parameters :

Uniform-Grid-Next-Cell( )

  1
IF

Next intersection is on the plane perpendicular to axis .   2
THEN

Function  returns the sign.   3          4
ELSE IF

Next intersection is on the plane perpendicular to axis .  5
THEN

6          7
ELSE IF

Next intersection is on the plane perpendicular to axis .  8
THEN

9


To summarize, a complete ray tracing algorithm is presented, which exploits the uniform grid generated during preprocessing and computes the ray-surface intersection closest to the ray origin. The minimum of ray parameters assigned to the coordinate planes, i.e. variable , determines the distance as far as the ray is inside the cell. This parameter is used to decide whether or not a ray-surface intersection is really inside the cell.

Ray-First-Intersection-with-Uniform-Grid( )

  1

Uniform-Grid-Enclosing-Cell(

)
2

Uniform-Grid-Ray-Parameter-Initialization(

)
3
WHILE

are inside the grid   4
DO

Here is the exit from the cell.   5
Initialization: no intersection yet.   6
FOR
each object  registered in cell ()   7
DO

Ray-Surface-Intersection(

)

Negative: no intersection.  8
IF

Is the new intersection closer?   9
THEN

The ray parameter of the closest intersection so far.  10
The first intersected object.  11
IF

Was intersection in the cell?  12
THEN

The position of the intersection.  13
RETURN

Termination.  14
Uniform-Grid-Next-Cell(

)

3DDDA.  15
RETURN
“no intersection” 

#### 22.6.2.3.  Time and storage complexity of the uniform grid algorithm.

The preprocessing phase of the uniform grid algorithm tests each object with each cell, thus runs in time where and are the numbers of objects and cells, respectively. In practice, the resolution of the grid is set to make proportional to since in this case, the average number of objects per cell becomes independent of the total number of objects. Such resolution makes the preprocessing time quadratic, that is . We note that sorting objects before testing them against cells may reduce this complexity, but this optimization is not crucial since not the preprocessing but the ray tracing time is critical. Since in the worst case all objects may overlap with each cell, the storage space is also in .

The ray tracing time can be expressed by the following equation:

where is the time needed to identify the cell containing the origin of the ray, is the number of ray–surface intersection tests until the first intersection is found, is the time required by a single ray–surface intersection test, is the number of visited cells, and is the time needed to step onto the next cell.

To find the first cell, the coordinates of the ray origin should be divided by the cell sizes, and the cell indices are obtained by rounding the results. This step thus runs in constant time. A single ray–surface intersection test also requires constant time. The next cell is determined by the 3DDDA algorithm in constant time as well. Thus the complexity of the algorithm depends only on the number of intersection tests and the number of the visited cells.

Considering a worst case scenario, a cell may contain all objects, requiring intersection test with objects. In the worst case the ray tracing has linear complexity. This means that the uniform grid algorithm needs quadratic preprocessing time and storage, but solves the ray tracing problem still in linear time as the naive algorithm, which is quite disappointing. However, uniform grids are still worth using since worst case scenarios are very unlikely. The fact is that classic complexity measures describing the worst case characteristics are not appropriate to compare the naive algorithm and the uniform grid based ray tracing. For a reasonable comparison, the probabilistic analysis of the algorithms is needed.

#### 22.6.2.4.  Probabilistic model of the virtual world.

To carry out the average case analysis, the scene model, i.e. the probability distribution of the possible virtual world models must be known. In practical situations, this probability distribution is not available, therefore it must be estimated. If the model of the virtual world were too complicated, we would not be able to analytically determine the average, i.e. the expected running time of the ray tracing algorithm. A simple, but also justifiable model is the following: Objects are spheres of the same radius , and sphere centres are uniformly distributed in space.

Since we are interested in the asymptotic behavior when the number of objects is really high, uniform distribution in a finite space would not be feasible. On the other hand, the boundary of the space would pose problems. Thus, instead of dealing with a finite object space, the space should also be expanded as the number of objects grows to sustain constant average spatial object density. This is a classical method in probability theory, and its known result is the Poisson point process.

Definition 22.16 A Poisson point process counts the number of points in subset of space in a way that

• is a Poisson distribution of parameter , where is a positive constant called “intensity” and is the volume of , thus the probability that contains exactly points is

and the expected number of points in volume is ;

• for disjoint sets random variables are independent.

Using the Poisson point process, the probabilistic model of the virtual world is:

1. The object space consists of spheres of the same radius .

2. The sphere centres are the realizations of a Poisson point process of intensity .

Having constructed a probabilistic virtual world model, we can start the analysis of the candidate algorithms assuming that the rays are uniformly distributed in space.

#### 22.6.2.5.  Calculation of the expected number of intersections.

Looking at Figure 22.32 we can see a ray that passes through certain cells of the space partitioning data structure. The collection of those sphere centres where the sphere would have an intersection with a cell is called the candidate space associated with this cell.

Only those spheres of radius can have intersection with the ray whose centres are in a cylinder of radius around the ray. This cylinder is called the intersection space (Figure 22.32). More precisely, the intersection space also includes two half spheres at the bottom and at the top of the cylinder, but these will be ignored.

As the ray tracing algorithm traverses the data structure, it examines each cell that is intersected by the ray. If the cell is empty, then the algorithm does nothing. If the cell is not empty, then it contains, at least partially, a sphere which is tried to be intersected. This intersection succeeds if the centre of the sphere is inside the intersection space and fails if it is outside.

The algorithm should try to intersect objects that are in the candidate space, but this intersection will be successful only if the object is also contained by the intersection space. The probability of the success is the ratio of the projected areas of the intersection space and the candidate space associated with this cell.

From the probability of the successful intersection in a non-empty cell, the probability that the intersection is found in the first, second, etc. cells can also be computed. Assuming statistical independence, the probabilities that the first, second, third, etc. intersection is the first successful intersection are , , , etc., respectively. This is a geometric distribution with expected value . Consequently, the expected number of the ray–object intersection tests is:

If the ray is parallel to one of the sides, then the projected size of the candidate space is where is the edge size of a cell and is the radius of the spheres. The other extreme case happens when the ray is parallel to the diagonal of the cubic cell, where the projection is a rounded hexagon having area . The success probability is then:

According to equation (22.27), the average number of intersection calculations is the reciprocal of this probability:

Note that if the size of the cell is equal to the diameter of the sphere (), then

This result has been obtained assuming that the number of objects converges to infinity. The expected number of intersection tests, however, remains finite and relatively small.

#### 22.6.2.6.  Calculation of the expected number of cell steps.

In the following analysis the conditional expected value theorem will be used. An appropriate condition is the length of the ray segment between its origin and the closest intersection. Using its probability density as a condition, the expected number of visited cells can be written in the following form:

where is the length of the ray and is its probability density.

Since the intersection space is a cylinder if we ignore the half spheres around the beginning and the end, its total volume is . Thus the probability that intersection occurs before is:

Note that this function is the cumulative probability distribution function of . The probability density can be computed as its derivative, thus we obtain:

The expected length of the ray is then:

In order to simplify the analysis, we shall assume that the ray is parallel to one of the coordinate axes. Since all cells have the same edge size , the number of cells intersected by a ray of length can be estimated as . This estimation is quite accurate. If the the ray is parallel to one of the coordinate axes, then the error is at most 1. In other cases the real value can be at most times the given estimation. The estimated expected number of visited cells is then:

For example, if the cell size is similar to the object size (), and the expected number of sphere centres in a cell is , then . Note that the expected number of visited cells is also constant even for infinite number of objects.

#### 22.6.2.7.  Expected running time and storage space.

We concluded that the expected numbers of required intersection tests and visited cells are asymptotically constant, thus the expected time complexity of the uniform grid based ray tracing algorithm is constant after quadratic preprocessing time. The value of the running time can be controlled by cell size according to equations (22.28) and (22.30). Smaller cell sizes reduce the average number of intersection tests, but increase the number of visited cells.

According to the probabilistic model, the average number of objects overlapping with a cell is also constant, thus the storage is proportional to the number of cells. Since the number of cells is set proportional to the number of objects, the expected storage complexity is also linear unlike the quadratic worst-case complexity.

The expected constant running time means that asymptotically the running time is independent of the number of objects, which explains the popularity of the uniform grid based ray tracing algorithm, and also the popularity of the algorithms presented in the next subsections.

#### 22.6.2.8.  Octree.

Uniform grids require many unnecessary cell steps. For example, the empty spaces are not worth partitioning into cells, and two cells are worth separating only if they contain different objects. Adaptive space partitioning schemes are based on these recognitions. The space can be partitioned adaptively following a recursive approach. This results in a hierarchical data structure, which is usually a tree. The type of this tree is the base of the classification of such algorithms.

The adaptive scheme discussed in this subsection uses an octal tree (octree for short), where non-empty nodes have 8 children. An octree is constructed by the following algorithm:

• For each object, an AABB is found, and object AABBs are enclosed by a scene AABB. The scene AABB is the cell corresponding to the root of the octree.

• If the number of objects overlapping with the current cell exceeds a predefined threshold, then the cell is subdivided to 8 cells of the same size by halving the original cell along each coordinate axis. The 8 new cells are the children of the node corresponding to the original cell. The algorithm is recursively repeated for the child cells.

• The recursive tree building procedure terminates if the depth of the tree becomes too big, or when the number of objects overlapping with a cell is smaller than the threshold.

The result of this construction is an octree (Figure 22.33). Overlapping objects are registered in the leaves of this tree.

When a ray is traced, those leaves of the tree should be traversed which are intersected by the ray, and ray–surface intersection test should be executed for objects registered in these leaves:

Ray-First-Intersection-with-Octree( )

  1  intersection of the ray and the scene AABB   2
WHILE

is inside of the scene AABB       Traversal of the tree.   3

Octree-Cell-Search(

)
4    ray parameter of the intersection of the  and the ray   5
Initialization: no ray–surface intersection yet.   6
FOR
each object  registered in    7
DO

Ray-Surface-Intersection(

)

Negative if no intersection exists.   8
IF

Is the new intersection closer?   9
THEN

Ray parameter of the closest intersection so far.  10
First intersected object so far.  11
IF

Has been intersection at all ?  12
THEN

Position of the intersection.  13
RETURN

14
A point in the next cell.  15
RETURN
“no intersection” 

The identification of the next cell intersected by the ray is more complicated for octrees than for uniform grids. The Octree-Cell-Search algorithm determines that leaf cell which contains a given point. At each level of the tree, the coordinates of the point are compared to the coordinates of the centre of the cell. The results of these comparisons determine which child contains the point. Repeating this test recursively, we arrive at a leaf sooner or later.

In order to identify the next cell intersected by the ray, the intersection point of the ray and the current cell is computed. Then, ray parameter of this intersection point is increased “a little” (this little value is denoted by in algorithm Ray-First-Intersection-with-Octree ). The increased ray parameter is substituted into the ray equation, resulting in point that is already in the next cell. The cell containing this point can be identified with Octree-cell-search .

Cells of the octree may be larger than the allowed minimal cell, therefore the octree algorithm requires less number of cell steps than the uniform grid algorithm working on the minimal cells. However, larger cells reduce the probability of the successful intersection tests since in a large cell it is less likely that a random ray intersecting the cell also intersects a contained object. Smaller successful intersection probability, on the other hand, results in greater expected number of intersection tests, which affects the performance negatively. It also means that non-empty octree cells are worth subdividing until the minimum cell size is reached even if the cell contains just a single object. Following this strategy, the size of the non-empty cells are similar, thus the results of the complexity analysis made for the uniform grid remain to be applicable to the octree as well. Since the probability of the successful intersection depends on the size of the non-empty cells, the expected number of needed intersection tests is still given by inequality (22.28). It also means that when the minimal cell size of an octree equals to the cell size of a uniform grid, then the expected number of intersection tests is equal in the two algorithms.

The advantage of the ocree is the ability to skip empty spaces, which reduces the number of cell steps. Its disadvantage is, however, that the time of the next cell identification is not constant. This identification requires the traversal of the tree. If the tree construction is terminated when a cell contains small number of objects, then the number of leaf cells is proportional to the number of objects. The depth of the tree is in , so is the time needed to step onto the next cell.

#### 22.6.2.9.  kd-tree.

An octree adapts to the distribution of the objects. However, the partitioning strategy of octrees always halves the cells without taking into account where the objects are, thus the adaptation is not perfect. Let us consider a partitioning scheme which splits a cell into two cells to make the tree balanced. Such method builds a binary tree which is called binary space partitioning tree, abbreviated as BSP-tree. If the separating plane is always perpendicular to one of the coordinate axes, then the tree is called kd-tree.

The separating plane of a kd-tree node can be placed in many different ways:

• the spatial median method halves the cell into two congruent cells.

• the object median method finds the separating plane to have the same number of objects in the two child cells.

• the cost driven method estimates the average computation time needed when a cell is processed during ray tracing, and minimizes this value by placing the separating plane. An appropriate cost model suggests to separate the cell to make the probabilities of the ray–surface intersection of the two cells similar.

The probability of the ray–surface intersection can be computed using a fundamental theorem of the integral geometry:

Theorem 22.17 If convex solid contains another convex solid , then the probability that a uniformly distributed line intersects solid provided that the line intersected equals to the ratio of the surface areas of objects and .

According to this theorem the cost driven method finds the separating plane to equalize the surface areas in the two children.

Let us now present a general kd-tree construction algorithm. Parameter identifies the current cell, is the current depth of recursion, and stores the orientation of the current separating plane. A is associated with its two children ( and ), and its left-lower-closer and right-upper-farther corners ( and ). Cells also store the list of those objects which overlap with the cell. The orientation of the separation plane is determined by a round-robin scheme implemented by function Round-robin providing a sequence like (). When the following recursive algorithm is called first, it gets the scene AABB in variable and the value of variable is zero:

Kd-Tree-Construction( )

  1
IF
the number of objects overlapping with  is small or  is large   2
THEN

RETURN
3  AABB of  and AABB of AABB of    4
IF

5
THEN

perpendicular separating plane of    6        perpendicular separating plane of    7
ELSE

IF

8
THEN

perpendicular separating plane of    9           perpendicular separating plane of   10
ELSE

IF

11
THEN

perpendicular separating plane of   12           perpendicular separating plane of   13
FOR
each object  of   14
DO

IF
object  is in the AABB of   15
THEN
assign object  to the list of   16
IF
object  is in the AABB of   17
THEN
assign object  to the list of   18
Kd-Tree-Construction(

, Round-Robin(

))
19
Kd-Tree-Construction(

, Round-Robin(

))



Now we discuss an algorithm that traverses the constructed kd-tree and finds the visible object. First we have to test whether the origin of the ray is inside the scene AABB. If it is not, the intersection of the ray and the scene AABB is computed, and the origin of the ray is moved there. The identification of the cell containing the ray origin requires the traversal of the tree. During the traversal the coordinates of the point are compared to the coordinates of the separating plane. This comparison determines which child should be processed recursively until a leaf node is reached. If the leaf cell is not empty, then objects overlapping with the cell are intersected with the ray, and the intersection closest to the origin is retained. The closest intersection is tested to see whether or not it is inside the cell (since an object may overlap more than one cells, it can also happen that the intersection is in another cell). If the intersection is in the current cell, then the needed intersection has been found, and the algorithm can be terminated. If the cell is empty, or no intersection is found in the cell, then the algorithm should proceed with the next cell. To identify the next cell, the ray is intersected with the current cell identifying the ray parameter of the exit point. Then the ray parameter is increased “a little” to make sure that the increased ray parameter corresponds to a point in the next cell. The algorithm keeps repeating these steps as it processes the cells of the tree.

This method has the disadvantage that the cell search always starts at the root, which results in the repetitive traversals of the same nodes of the tree.

This disadvantage can be eliminated by putting the cells to be visited into a stack, and backtracking only to the point where a new branch should be followed. When the ray arrives at a node having two children, the algorithm decides the order of processing the two child nodes. Child nodes are classified as “near” and “far” depending on whether or not the child cell is on the same side of the separating plane as the origin of the ray. If the ray intersects only the near'' child, then the algorithm processes only that subtree which originates at this child. If the ray intersects both children, then the algorithm pushes the “far” node onto the stack and starts processing the “near” node. If no intersection exists in the “near” node, then the stack is popped to obtain the next node to be processed.

The notations of the ray tracing algorithm based on kd-tree traversal are shown by Figure 22.35. The algorithm is the following:

Ray-First-Intersection-with-kd-Tree( )

  1

Ray-AABB-Intersection(

)

Intersection with the scene AABB.  2
IF
no intersection   3
THEN

RETURN
“no intersection”   4
Push(

)
5
WHILE
the stack is not empty       Visit all nodes.   6
DO

Pop(

)
7
WHILE

is not a leaf   8
DO

orientation of the separating plane of the    9               10
Ray parameter of the separating plane.  11
IF

Is  on the left side of the separating plane?  12
THEN

Left.  13
ELSE

Right.  14
IF

or   15
THEN

The ray intersects only the  cell.  16
ELSE

IF

17
THEN

The ray intersects only the  cell.  18
ELSE

Push(

)

The ray intersects both cells.  19
First  is intersected.  20
The ray exists at  from the  cell.                     If the current cell is a leaf. 21
Maximum ray parameter in this cell.  22
FOR
each object  of   23
DO

Ray-surface-intersection(

)

Negative if no intersection exists. 24
IF

Is the new intersection closer to the ray origin?  25
THEN

The ray parameter of the closest intersection so far. 26
The object intersected closest to the ray origin. 27
IF

Has been intersection at all in the cell?  28
THEN

The intersection point.  29
RETURN

Intersection has been found.  30
RETURN
“no intersection”       No intersection. 

Similarly to the octree algorithm, the likelihood of successful intersections can be increased by continuing the tree building process until all empty spaces are cut (Figure 22.36).

Our probabilistic world model contains spheres of same radius , thus the non-empty cells are cubes of edge size . Unlike in uniform grids or octrees, the separating planes of kd-trees are not independent of the objects. Kd-tree splitting planes are rather tangents of the objects. This means that we do not have to be concerned with partially overlapping spheres since a sphere is completely contained by a cell in a kd-tree. The probability of the successful intersection is obtained applying Theorem 22.17. In the current case, the containing convex solid is a cube of edge size , the contained solid is a sphere of radius , thus the intersection probability is:

The expected number of intersection tests is then:

We can conclude that the kd-tree algorithm requires the smallest number of ray-surface intersection tests according to the probabilistic model.

Exercises

22.6-1 Prove that the expected number of intersection tests is constant in all those ray tracing algorithms which process objects in the order of their distance from the ray origin.

22.6-2 Propose a ray intersection algorithm for subdivision surfaces.

22.6-3 Develop a ray intersection method for B-spline surfaces.

22.6-4 Develop a ray intersection algorithm for CSG models assuming that the ray–primitive intersection tests are already available.

22.6-5 Propose a ray intersection algorithm for transformed objects assuming that the algorithm computing the intersection with the non-transformed objects is available (hints: transform the ray).