Implementing and Using BSP Trees
Fast 3-D sorting
Nathan Dwyer
Nathan is a programmer at Starware in Bellevue, WA and can be reached at [email protected].
With the arrival of graphic-intensive applications like DOOM, high-speed 3-D display engines are becoming a necessity for PC systems. These 3-D engines owe their speed, in part, to binary space partitioning (BSP) trees, data structures that provide an extremely fast method for one component of the rendering process.
Rendering 3-D shapes usually consists of several steps. First, the shape is oriented by rotation about its center. Next, its location is determined relative to a viewpoint and view direction. Finally, the shape's component polygons are sorted, and each is projected and rendered. Writing a fast 3-D engine consists of speeding up each step as much as possible. BSP trees address the issue of sorting.
Sorting the polygons may seem trivial, but is usually fairly complicated. The difficulty is in determining which polygons lie behind which. If all the polygons are small and far apart, the task is simple. In that case, it's usually sufficient to compute the distance from the viewpoint to the center of each polygon, rendering the farther polygons first. This technique is called the "painter's algorithm" because polygons are rendered in the same order a painter would paint them on a flat wall. Consider Figure 1, however. In the first case, the center of the large polygon is closest to the viewpoint, yet the large polygon should clearly be drawn before the small one. In the second case, where the three polygons overlap in a circular manner, the painter's algorithm will be incorrect, no matter what order the polygons are drawn in.
A number of algorithms have been developed to attack this problem, but the most popular rely on two techniques: z-buffering and subdivision.
Z-buffering algorithms represent each display pixel by one element in an array. The array is used to keep track of the distance from the viewpoint to the polygon to be drawn at that point. As a polygon is drawn to the screen, the distance at each pixel is computed and compared with its corresponding array element. A new polygon pixel is rendered only if it is closer to the viewpoint than the last polygon drawn at that pixel. While this is an extremely elegant solution, it requires a lot of memory. Each pixel may require 16 or 32 bits to accurately represent the distance, and the extra distance computation can be time intensive.
Subdivision algorithms analyze a group of polygons, breaking them into smaller pieces that don't overlap. This may entail drawing the polygons concurrently, splitting each scanline into nonoverlapping, horizontal sections. Alternatively, each polygon may be checked against every other polygon for overlap. Both approaches require significant extra processing.
BSP trees were developed to determine the visibility of surfaces. They were later adapted to represent arbitrary 3-D shapes. BSP trees trade some (possibly expensive) preprocessing for very fast polygon sorting, making them best suited to static shapes: floor plans, walls, and other objects that don't change.
Generating a BSP tree is straightforward. From a list of polygons, choose one to be the root. Then separate the remaining polygons into two groups: those in front and those behind the plane of the root polygon, with respect to the root polygon's normal vector. Next, recursively build BSP trees out of each group. The root of each subgroup becomes a child of the root containing it. If any polygon isn't completely in front of or behind the plane of the root, split it into two polygons that are. This defines the criteria for choosing a root polygon: The best choice causes the fewest number of splits in the remaining polygons.
Almost all the difficulty in creating the BSP tree lies in splitting the polygons intersected by a root plane. Three pieces of code bear examination: the plane-plane intersection code, the code that tests for inaccurate intersections, and the code that actually splits a polygon.
Plane-plane intersection can be a bear. If brute-force linear algebra is applied, much complexity results. Luckily, plane-plane intersection can be simplified to line-plane intersection, which is much easier. The function BSPNode::_SplitPoly() in node.cpp (Listing Two) tests each line segment in the boundary of a polygon to see if it intersects the plane of a second (Listing One is the bsp.h include file.) The math involved is described in the accompanying text box entitled, "Intersecting a Line and a Plane." The complete source code, including sample data files, that implements this approach to BSP is provided electronically; see "Availability," page 3.
If a single vertex or edge from the polygon is contained in the plane, the polygon need not be split. If the variable t is very close to 1.0 or 0.0, a vertex is contained in the plane. While it would be nice to simply skip these occurrences, a plane may cut across the polygon and just happen to pass through a vertex. To guard against this possibility, the Boolean variable LastSideParallel indicates whether the last polygon edge was parallel to the plane. If so (and t equals 0.0), the intersection can be safely ignored.
The function BSPNode::Split() in node.cpp contains the code that creates two polygons from one. It calls BSPNode::_SplitPoly() to fill in a sparse array of CPoint objects with intersection points. BSPNode::_SplitPoly() returns a set of flags that indicates which boundary segments were split. BSPNode::Split() builds two polygons from this information by following the original list of points around the polygon and adding them to the first new polygon. When an intersection point is encountered, it's added to both new polygons. Subsequent points are added to the second new polygon until another intersection point is encountered. Once again, the intersection point is added to both new polygons, and the first polygon becomes the destination for subsequent points. Once the point list is exhausted, the new polygons have complete point lists of their own. The two new polygons have the same original normal, but their centers must be recalculated.
Sorting the polygons contained in a BSP tree is simply an in-order, depth-first traversal of the tree. At each node, a dot product is computed to determine if the viewpoint is in front of or behind the node's polygon. If in front, the polygons behind the node are drawn first, then the node's polygon, then the polygons in front of the node. If behind, the node's children are visited in reverse order.
To see how this works, look at Figure 2. If the viewpoint is at position A, it's behind polygon 1. Therefore, polygon 2 is drawn first, then polygon 1, and finally polygon 3, which is correct. If the viewpoint is at position B, the polygon order becomes 3, 1, 2, which is also correct.
All of the code to traverse the tree is contained in the function BSPNode::Traverse() in node.cpp. The vector from each polygon's center to the viewpoint is computed first. The dot product between this vector and the polygon's normal vector indicates the side on which the viewpoint is located. Here, the center point is used for clarity, but any vertex of the polygon may be used with the same result.
Let's take a moment to discuss file formats. The code reads and writes text files so the input and output may be understood easily. An input file of unprocessed polygons has the format shown in Example 1(a). Once the BSP tree has been created, it's saved in the format shown in Example 1(b), which can also be read back in.
There are many extensions to this basic data structure that can further increase rendering speed. Polygons can be back-face culled very easily. If the viewpoint is behind a node, that node's polygon should not be drawn. Bounding boxes can be generated for each level of the tree, allowing branches which don't intersect the viewcone to be skipped entirely. It should also be noted that limited polygon movement is possible, as long as the polygon never leaves its original plane and doesn't cross the plane of a root polygon.
While BSP trees are primarily limited to static shapes and require extra preprocessing, they provide a substantial speed increase without a lot of memory overhead. Preprocessing is fine if the goal is speed in the final presentation. Finally, since the implementation of BSP-tree algorithms is not overly difficult, the result is easy maintenance, as well as a large speed increase for a small amount of work.
Intersecting a Line and a Plane
A plane is defined by a point A and a normal N. Given a point X, X is in the plane if it satisfies the requirement in Example 2(a). A line is defined as Example 2(b), which can be rewritten as the set of equations in Example 2(c). Substituting these equations into the plane equation and solving for t results in Example 2(d). If t is in the interval [0, 1], then the line intersects the plane. If the denominator is 0, it is because the dot product of the plane's normal and the direction of the line is 0. This means the line and the plane are parallel. Either the line is contained in the plane or it doesn't intersect it at all.
--N.D.
Figure 1: Sorting polygons.
Figure 2: In-order, depth-first traversal of the tree.
Example 1: (a) Format of input file of unprocessed polygons; (b) format after BSP tree has been created.
(a) ushort NumPolygons
for each polygon:
ushort NumPoints
for each point:
long X Y Z
(b) for each node:
ushort NodeIndex
ushort FrontIndex, BackIndex
ushort NumPoints
for each point:
long X Y Z
Example 2: Math required to intersect a line and a plane.
(a) N·(X-A)=0
N·X-N·A=0
n1*x+n2*y+n3*z-(N·A)=0
(b) X=X1+t(X2-X1)
(c) x=x1+t*i
y=y1+t*j
z=z1+t*k
(d) t=-(n1*x1+n2*y1+n2*z1+(-(N·A)))
____________________________
(n1*i+n2*j+n3*k )
Listing One
/* BSP.H -- Type definitions for BSP code */
#include <fstream.h>
#define TRUE 1
#define FALSE 0
typedef unsigned char bool;
typedef unsigned short ushort;
typedef unsigned long ulong;
class CPoint
{
public:
double x, y, z;
public:
CPoint( void );
void Read( ifstream& Input );
void Write( ofstream& Output ) const;
double Magnitude( void ) const;
void Normalize( void );
double DotProduct( const CPoint& Pnt ) const;
CPoint CrossProduct( const CPoint& Pnt ) const;
short operator==( const CPoint& Pnt ) const;
CPoint operator-( const CPoint& Pnt ) const;
};
class BSPNodeList;
class BSPNode
{
private:
ushort Index;
BSPNode *FrontNode, *BackNode;
short PntCnt;
CPoint *PntList;
CPoint Center;
CPoint Normal;
double D;
ulong _SplitPoly( BSPNode *Plane, CPoint *SplitPnts );
void _ComputeCenter( void );
void _ComputeNormal( void );
void _ComputeD( void );
public:
BSPNode( void );
~BSPNode( void );
void ReadPoly( ifstream& Input );
void ReadNode( ifstream& Input, const BSPNodeList& NodePool );
void WriteNode( ofstream& Output, short& CurIndex );
CPoint GetCenter( void ) { return Center; }
CPoint GetNormal( void ) { return Normal; }
bool Intersects( BSPNode *Plane );
BSPNode *Split( BSPNode *Plane );
BSPNode *GetFront( void ) { return FrontNode; }
BSPNode *GetBack( void ) { return BackNode; }
void SetFront( BSPNode *Node ) { FrontNode = Node; }
void SetBack( BSPNode *Node) { BackNode = Node; }
void Traverse( const CPoint& CameraLoc );
};
class BSPNodeHeader
{
friend class BSPListIterator;
friend class BSPNodeList;
private:
BSPNodeHeader *Next;
BSPNode *Data;
BSPNodeHeader( BSPNode *DataNode ) { Data=DataNode; Next=NULL;}
void SetNext( BSPNodeHeader *NextNode ) { Next = NextNode; }
BSPNodeHeader *GetNext( void ) { return Next; }
BSPNode *GetData( void ) { return Data; }
};
class BSPNodeList
{
friend class BSPListIterator;
private:
BSPNodeHeader *FirstNode, *LastNode;
public:
BSPNodeList( void );
~BSPNodeList( void );
void ReadPolys( ifstream& Input );
bool Empty( void ) { return FirstNode == NULL; }
void Insert( BSPNode *Node );
void Remove( BSPNode *Node );
};
class BSPListIterator
{
private:
BSPNodeHeader *Current;
public:
BSPListIterator( void );
BSPNode *First( const BSPNodeList *List );
BSPNode *Next( void );
};
class BSPTree
{
private:
BSPNodeList Nodes;
BSPNode *Root;
BSPNode *_FindRoot( BSPNodeList *List );
BSPNode *_BuildBSPTree( BSPNodeList *List );
public:
void ReadPolyList( ifstream& Input );
void ReadTree( ifstream& Input );
void WriteTree( ofstream& Output );
void BuildTree( void );
void Traverse( const CPoint& CameraLoc );
};
Listing Two
/* BSPNode.cpp -- Implementation of BSPNode class */
#include <assert.h>
#include "bsp.h"
//-------------- Private Functions
ulong BSPNode::_SplitPoly( BSPNode *Plane, CPoint *SplitPnts )
// This is limited to convex polygons with no more than 32 sides
{
ulong Sides = 0;
bool LastSideParallel = FALSE;
if( !( Normal == Plane->Normal ) )
{
CPoint EdgeDelta;
double numer, denom, t;
for( ushort vertex=0; vertex<PntCnt; vertex++ )
{
ushort prevVertex = vertex ? vertex - 1 : PntCnt - 1;
EdgeDelta = PntList[ vertex ] - PntList[ prevVertex ];
denom = EdgeDelta.DotProduct( Plane->Normal );
if( denom )
{
numer = PntList[ prevVertex ].DotProduct
( Plane->Normal ) + Plane->D;
t = - numer / denom;
if( !( LastSideParallel && t == 0.0 ) )
{
if( t >= 0.0 && t < 0.999999 )
{
Sides |= 1 << vertex;
if( SplitPnts )
{
SplitPnts[ vertex ].x = PntList[
prevVertex ].x + SplitPnts[
vertex ].y = PntList[
prevVertex ].y + t * EdgeDelta.y;
SplitPnts[ vertex ].z = PntList[
prevVertex ].z + t * EdgeDelta.z;
}
}
}
}
LastSideParallel = ( denom == 0 );
}
}
return Sides;
}
void BSPNode::_ComputeCenter( void )
{
Center.x = Center.y = Center.z = 0.0;
for( ushort i=0; i<PntCnt; i++ )
{
Center.x += PntList[ i ].x;
Center.y += PntList[ i ].y;
Center.z += PntList[ i ].z;
}
Center.x /= PntCnt;
Center.y /= PntCnt;
Center.z /= PntCnt;
}
void BSPNode::_ComputeNormal( void )
{
CPoint a, b;
assert( PntCnt >= 3 );
a = PntList[ 0 ] - PntList[ 1 ];
b = PntList[ 2 ] - PntList[ 1 ];
Normal = a.CrossProduct( b );
Normal.Normalize();
}
void BSPNode::_ComputeD( void )
{
D = -Normal.DotProduct( Center );
}
//-------------- Public Functions
BSPNode::BSPNode( void ) :
FrontNode( NULL ),
BackNode( NULL ),
Index( 0 ),
PntCnt( 0 ),
PntList( NULL )
{
}
BSPNode::~BSPNode( void )
{
if( PntList )
delete[] PntList;
};
void BSPNode::ReadPoly( ifstream& Input )
{
Input >> PntCnt;
assert( PntCnt >= 3 );
PntList = new CPoint[ PntCnt ];
for( ushort i=0; i<PntCnt; i++ )
PntList[ i ].Read( Input );
_ComputeCenter();
_ComputeNormal();
_ComputeD();
}
void BSPNode::ReadNode( ifstream& Input, const BSPNodeList& NodePool )
{
short FrontIndex, BackIndex;
Input >> Index >> FrontIndex >> BackIndex;
if( !Input )
return;
Input >> PntCnt;
PntList = new CPoint[ PntCnt ];
for( short i=0; i<PntCnt; i++ )
PntList[ i ].Read( Input );
_ComputeCenter();
_ComputeNormal();
_ComputeD();
// Find Children
BSPListIterator Iter;
BSPNode *TestNode = Iter.First( &NodePool );
FrontNode = BackNode = NULL;
while( TestNode )
{
if( TestNode->Index == FrontIndex )
FrontNode = TestNode;
else if( TestNode->Index == BackIndex )
BackNode = TestNode;
if( FrontNode && BackNode )
return;
TestNode = Iter.Next();
}
}
void BSPNode::WriteNode( ofstream& Output, short& CurIndex )
{
if( FrontNode )
FrontNode->WriteNode( Output, CurIndex );
if( BackNode )
BackNode->WriteNode( Output, CurIndex );
// write index and child indices
Index = CurIndex++;
Output << Index << '\n';
Output << ( FrontNode ? FrontNode->Index : -1 ) << ' ';
Output << ( BackNode ? BackNode->Index : -1 ) << '\n';
// write point list
Output << PntCnt << '\n';
for( short i=0; i<PntCnt; i++ )
{
Output << '\t';
PntList[ i ].Write( Output );
}
Output << '\n';
}
bool BSPNode::Intersects( BSPNode *Plane )
{
return ( _SplitPoly( Plane, NULL ) != 0 );
}
BSPNode *BSPNode::Split( BSPNode *Plane )
{
BSPNode *NewNode;
CPoint *SplitPnts;
ulong Splits;
SplitPnts = new CPoint[ PntCnt ];
Splits = _SplitPoly( Plane, SplitPnts );
if( Splits )
{
CPoint *NewPoly1, *NewPoly2;
ushort Poly1Index = 0, Poly2Index = 0;
ushort Destination = 0;
NewPoly1 = new CPoint[ PntCnt ];
NewPoly2 = new CPoint[ PntCnt ];
for( ushort i=0; i<PntCnt; i++ )
{
// Handle split points
if( Splits & ( 1 << i ) )
{
NewPoly1[ Poly1Index++ ] = SplitPnts[ i ];
NewPoly2[ Poly2Index++ ] = SplitPnts[ i ];
Destination ^= 1;
}
if( Destination )
NewPoly1[ Poly1Index++ ] = PntList[ i ];
else
NewPoly2[ Poly2Index++ ] = PntList[ i ];
}
// Make New node
NewNode = new BSPNode;
NewNode->PntCnt = Poly1Index;
NewNode->PntList = NewPoly1;
NewNode->Normal = Normal;
NewNode->_ComputeCenter();
delete[] PntList;
PntCnt = Poly2Index;
PntList = NewPoly2;
_ComputeCenter();
}
delete SplitPnts;
return NewNode;
}
void BSPNode::Traverse( const CPoint& CameraLoc )
{
CPoint VecToCam = CameraLoc - Center;
if( VecToCam.DotProduct( Normal ) < 0 )
{
if( FrontNode )
FrontNode->Traverse( CameraLoc );
// Process 'this' i.e. render it to screen
if( BackNode )
BackNode->Traverse( CameraLoc );
}
else
{
if( BackNode )
BackNode->Traverse( CameraLoc );
// Process 'this' i.e. render it to screen
if( FrontNode )
FrontNode->Traverse( CameraLoc );
}
}
Copyright © 1995, Dr. Dobb's Journal
