Listing 1

/** 
 * Gets the volume containing all objects that can cast a shadow 
 * into the view frustum.
 * @param frustum view frustum 
 * @return shadow caster volume 
 */ 
Polyhedron<double> PointLightSource::shadowCasterVolume( Polyhedron<double>const& frustum ) const 
{ 
	// If lightsource lies inside frustum, then simply return frustum.
	if ( frustum.pointInside( _position ) ) { 
		return frustum; 
	} 
	
	// Create a new polyhedron, using all frustum planes for which the light 
	// source lies on the positive side of (or lies exactly on the plane), 
	// and creating new planes for each edge that is shared by one plane that 
	// points towards the light and one that points away. 
	
	// First copy all planes pointing towards the light source 
	Polyhedron<double> shadowVol; 
	Polyhedron<double>::PlanesItors planesItors = frustum.getPlanes(); 
	for ( Polyhedron<double>::PlanesItor p=planesItors.first; p!=planesItors.second; ++p ) 
	{ 
		if ( p->dist(_position) >= 0 ) { 
			shadowVol.addPlane( *p ); 
		} 
	} 
	
	// Calculate the average of the frustum vertices. 
	// This is used for determining the direction of the 
	// normals for the new planes. 
	Polyhedron<double>::VerticesItors const verticesItors = frustum.getVertices(); 
	Polyhedron<double>::Vertex verticesAverage( 0.f ); 
	double numVertices = 0; 
	for ( Polyhedron<double>::VerticesItor v=verticesItors.first; v!=verticesItors.second; ++v ) 
	{ 
		verticesAverage += *v; 
		++numVertices; 
	} 
	verticesAverage /= numVertices; 
	
	// Now add new planes for all silloute edges 
	Polyhedron<double>::EdgesItors const edgesItors = frustum.getEdges(); 
	for ( Polyhedron<double>::EdgesItor e=edgesItors.first; e!=edgesItors.second; ++e ) 
	{ 
		// One plane faces light and other faces away (or is parallel) ? 
		bool const p0pos = e->plane[0]->dist(_position) >= 0; 
		bool const p1pos = e->plane[1]->dist(_position) >= 0; 
		if ( p0pos != p1pos ) { 
		
			// Create new plane 
			Vector<3, double> n( (e->end[0]-_position).cross(e->end[0]-e->end[1]).norm() ); 
			Plane<double> p( n, e->end[0] ); 
			// Ensure the normal is pointing in the correct direction. 
			if ( p.dist( verticesAverage ) < 0 ) { 
				p = -p; 
			} 
			
			// Save the plane 
			shadowVol.addPlane( p ); 
		} 
	} 
	
	return shadowVol; 
} 

Listing 2

/** 
 * Convex polyhedron. 
 * 
 * Note that there is no checking to enforce convexity. 
 * 
 * Polyhedron is defined as the intersection of the 
 * positive half-spaces for a list of planes. 
 * That is, plane normals must point inwards. 
 * 
 * Degenerate infinite polyhedra are supported, the 
 * simplest case of this being a half-space. 
 */ 
template < class T= double > 
class Polyhedron { 

	typedef std::vector< Plane<T> > Planes; 
	
public : 

	/** 
	 * Constructor 
	 */ 
	Polyhedron() 
	: _recalculateVertices( true ) 
	, _recalculateEdges( true ) 
	{ 
	} 
	
	/** 
	* Adds a plane to polyhedron. 
	* @param p plane 
	*/ 
	void addPlane( Plane<T> const& p ) { 
		_planes.push_back( p ); 
		_recalculateVertices = true ; 
	} 
	
	/// Const iterator over a collection of planes 
	typedef typename Planes::const_iterator PlanesItor; 
	
	/// Pair of PlanesItor's, where first is begin, and second is end. 
	typedef std::pair<PlanesItor,PlanesItor> PlanesItors; 
	
	/** 
	 * Get the begin and end const iterators over the planes. 
	 * @return plane begin and end iterators 
	 */ 
	PlanesItors getPlanes() const { 
		return PlanesItors( _planes.begin(), _planes.end() ); 
	} 
	
	/** 
	 * Test if a point lies in (or on) the polyhedron. 
	 * @param point point 
	 * @param error each "point" is given a radius of 
	 * max(error, error*plane's distance to origin) 
	 * @return true iff point lies in or on the polyhedron 
	 */ 
	bool pointInside( Vector<3,T> const& point, T error=1e-3 ) const { 
		for ( PlanesItor i=_planes.begin(); i!=_planes.end(); ++i ) { 
			if ( i->dist(point) < -std::max( error, std::fabs(error*i->d()) ) ) 
			{ 
				return false ; 
			} 
		} 
		return true ; 
	} 
	
private : 
	
	/// Vertex of the polyhedron 
	struct PlanesVertex { 
		PlanesItor 	plane[3]; 	///< planes that meet at the vertex 
		Vector<3,T> 	pos; 		///< vertex position 
	}; 
	
	/// Collection of PlanesVertex's 
	typedef std::vector<PlanesVertex> PlanesVertices; 
	
public : 

	/// Vertex of the polyhedron 
	typedef Vector<3,T> Vertex; 

private : 

	struct VertexCmp { 
		/** 
		 * Vertex comparison operator for use in std::set. 
		 * @param v0 first vertex 
		 * @param v1 second vertex 
		 * @return true iff v0 "is less than" v1 
		 */ 
		bool operator ()( Vertex const & v0, Vertex const & v1 ) const { 
			if ( v0[0] < v1[0] ) return true ; 
			if ( v0[0] > v1[0] ) return false ; 
			if ( v0[1] < v1[1] ) return true ; 
			if ( v0[1] > v1[1] ) return false ; 
			return v0[2] < v1[2]; 
		} 
	}; 
	
	/// Collection of Vertex's 
	typedef std::set<Vertex,VertexCmp> Vertices; 
	
public : 
	
	/// Const iterator over a collection of vertices 
	typedef typename Vertices::const_iterator VerticesItor; 
	
	/// Pair of VerticesItor's, where first is begin, and second is end. 
	typedef std::pair<VerticesItor,VerticesItor> VerticesItors; 
	
	/** 
	 * Gets the begin and end const iterators overt the vertices 
	 * @return vertex begin and end iterators 
	 */ 
	VerticesItors getVertices() const { 
		if ( _recalculateVertices ) { 
			_recalculateVertices = false ; 
			_recalculateEdges = true ; 
			
			// First calculate _planesVertices 
			_planesVertices.clear(); 
			for ( PlanesItor p0=_planes.begin(); p0!=_planes.end(); ++p0 ) { 
				PlanesItor p1 = p0; 
				for ( ++p1; p1!=_planes.end(); ++p1 ) { 
					PlanesItor p2 = p1; 
					for ( ++p2; p2!=_planes.end(); ++p2 ) { 
						// Test if the planes intersect at a point 
						Matrix<3,3,T> m; 
						Vector<3,T> const& n0( p0->n() ); 
						Vector<3,T> const& n1( p1->n() ); 
						Vector<3,T> const& n2( p2->n() ); 
						m[0][0] = n0[0]; m[0][1] = n0[1]; m[0][2] = n0[2]; 
						m[1][0] = n1[0]; m[1][1] = n1[1]; m[1][2] = n1[2]; 
						m[2][0] = n2[0]; m[2][1] = n2[1]; m[2][2] = n2[2]; 
						
						// Find the point of intersection 
						typename Matrix<3,3,T>::Lup const lup( m, false ); 
						if ( !lup.failed() ) { 
							Vector<3,T> const d(-p0->d(),-p1->d(),-p2->d()); 
							Vector<3,T> const pos( lup.solve(d) ); 
							
							// Now check that the point is in or on 
							// the polyhedron 
							if ( pointInside( pos ) ) { 
								// Save the vertex 
								PlanesVertex const pv = { {p0,p1,p2}, pos }; 
								_planesVertices.push_back( pv ); 
							} 
						} 
					} 
				} 
			} 
			
			// Now calculate unique set of vertex positions 
			_vertices.clear(); 
			for ( typename PlanesVertices::const_iterator pv=_planesVertices.begin(); pv!=_planesVertices.end(); ++pv ) 
			{ 
				_vertices.insert( Vertex(pv->pos) ); 
			} 
		} 
		return VerticesItors( _vertices.begin(), _vertices.end() ); 
	} 
	
	/// Edge of polyhedron 
	struct Edge { 
		PlanesItor 	plane[2]; 	///< planes that meet to form the edge 
		Vector<3,T> 	end[2]; 	///< vertex positions 
	}; 
	
private : 

	struct EdgeCmp { 
		/** 
		 * Edge comparison operator for use in std::set. 
		 * @param e0 first edge 
		 * @param e1 second edge 
		 * @return true iff e0 "is less than" e1 
		 */ 
		bool operator ()( Edge const & e0, Edge const & e1 ) const { 
			for ( unsigned i=0; i<2; ++i ) { 
				for ( unsigned j=0; j<3; ++j ) { 
					if ( e0.plane[i]->n()[j] < e1.plane[i]->n()[j] ) 
						return true ; 
					if ( e0.plane[i]->n()[j] > e1.plane[i]->n()[j] ) 
						return false ; 
				} 
				if ( e0.plane[i]->d() < e1.plane[i]->d() ) return true ; 
				if ( e0.plane[i]->d() > e1.plane[i]->d() ) return false ; 
			} 
			for ( unsigned i=0; i<2; ++i ) { 
				for ( unsigned j=0; j<3; ++j ) { 
					if ( e0.end[i][j] < e1.end[i][j] ) return true ; 
					if ( e0.end[i][j] > e1.end[i][j] ) return false ; 
				} 
			} 
			return false ; 
		} 
	}; 
	
	/// Collection of Edge's 
	typedef std::set<Edge,EdgeCmp> Edges; 
	
public : 

	/// Const iterator over collection of edges 
	typedef typename Edges::const_iterator EdgesItor; 
	
	/// Pair of EdgesItor's, where first is begin, and second is end. 
	typedef std::pair<EdgesItor,EdgesItor> EdgesItors; 
	
	/** 
	 * Gets the begin and end const iterators overt the edges 
	 * @return edge begin and end iterators 
	 */ 
	EdgesItors getEdges() const { 
		getVertices(); 
		if ( _recalculateEdges ) { 
			_recalculateEdges = false ; 
			
			_edges.clear(); 
			
			for ( typename PlanesVertices::const_iterator v0=_planesVertices.begin(); v0!=_planesVertices.end(); ++v0 ) 
			{ 
				typename PlanesVertices::const_iterator v1 = v0; 
				for ( ++v1; v1!=_planesVertices.end(); ++v1 ) { 

					// Check that the two vertices are not the same point 
					// (which will happen if more than three planes meet 
					// at the same point) 
					if ( v0->pos != v1->pos ) { 
						
						// For every common pair of planes, add an edge 
						#define _POLYHEDRON_CMP_PLANE(n) \ 
							( v0->plane[n] == v1->plane[0] || \ 
							v0->plane[n] == v1->plane[1] || \ 
							v0->plane[n] == v1->plane[2] ) 
						#define _POLYHEDRON_ADD_EDGE(i,j) \ 
							do { \ 
								Edge const e = { { v0->plane[i], \ 
										   v0->plane[j] }, \ 
										 { v0->pos, v1->pos } }; \ 
								_edges.insert( e ); \ 
							} while (0) 
						if ( _POLYHEDRON_CMP_PLANE(0) ) { 
							if ( _POLYHEDRON_CMP_PLANE(1) ) { 
								_POLYHEDRON_ADD_EDGE(0,1); 
							} 
							else if ( _POLYHEDRON_CMP_PLANE(2) ) { 
								_POLYHEDRON_ADD_EDGE(0,2); 
							} 
						} 
						else if ( _POLYHEDRON_CMP_PLANE(1) ) { 
							if ( _POLYHEDRON_CMP_PLANE(2) ) { 
								_POLYHEDRON_ADD_EDGE(1,2); 
							} 
						} 
						#undef _POLYHEDRON_ADD_EDGE 
						#undef _POLYHEDRON_CMP_PLANE 
					} 
				} 
			} 
		} 
		return EdgesItors( _edges.begin(), _edges.end() ); 
	} 
private : 

	Planes 			_planes; 
	mutable bool 		_recalculateVertices; 
	mutable PlanesVertices 	_planesVertices; 
	mutable Vertices 	_vertices; 
	mutable bool 		_recalculateEdges; 
	mutable Edges 		_edges; 
}; 

Listing 3

/** 
 * LUP decomposition of a matrix. 
 * For a matrix non-singular square matrix M, it's decomposition satisfies 
 * PM = LU 
 */ 
template < int ROWS, int COLS, class T > 
class Matrix<ROWS,COLS,T>::Lup { 
BOOST_STATIC_ASSERT( ROWS == COLS ); 

public : 

	/** 
	 * Constructor 
	 * @param m matrix to decompose 
	 * @param canThrow if set, MatrixSingularException will be thrown if 
	 * LUP decomposition fails (the default), otherwise, 
	 * use the failed() method to test whether or not 
	 * the decomposition was successful 
	 */ 
	explicit Lup( Matrix const & m, bool canThrow= true ) 
	: _failed( true ) 
	, _lu( m ) 
	{ 
		// This is an implementation of the LUP-DECOMPOSITION algorithm 
		// presented in Cohen, Leiserson, Rivest and Stein. "Introduction To 
		// Algorithms", 2nd Ed. MIT Press. 2001. 
		
		// Initialise _p to identity transform 
		for ( unsigned i=0; i<ROWS; ++i ) { 
			_p[i] = i; 
		} 

		// Loop through each column 
		// While the last iteration does not modify the matrix _lu in anyway, 
		// it is required to check for singularity. 
		for ( int k=0; k<COLS; ++k ) { 
		
			// Find the row with the largest absolute value element 
			T maxElem = 0; 
			int maxRow; 
			for ( int i=k; i<ROWS; ++i ) { 
				T const absmik = std::fabs( _lu[i][k] ); 
				if ( absmik > maxElem ) { 
					maxElem = absmik; 
					maxRow = i; 
				} 
			} 
			if ( equals(maxElem,0) ) { 	// maxElem approx equal to zero 
				if ( canThrow ) { 
					throw MatrixSingularException(); 
				} 
				else { 
					return ; 
				} 
			} 
			
			// Swap rows k and maxRow 
			std::swap( _p[k], _p[maxRow] ); 
			_lu.swapRows( k, maxRow ); 
			
			// Perform LU decomposition 
			T const mkk = _lu[k][k]; 
			for ( int i=k+1; i<ROWS; ++i ) { 
				_lu[i][k] /= mkk; 
				T const mik = _lu[i][k]; 
				for ( int j=k+1; j<COLS; ++j ) { 
					_lu[i][j] -= mik * _lu[k][j]; 
				} 
			} 
		} 
		_failed = false ; 
	} 
	
	/** 
	 * Test whether or not the LUP decomposition was successful. 
	 * Note that the constructor was not called with canThrow=false, 
	 * then failed() will always return false. 
	 * @return false iff LUP decomposition was successful 
	 */ 
	bool failed() const { 
		return _failed; 
	} 
	
	/** 
	 * Solve decomposed set of linear equations. 
	 * @param c vector containing constants 
	 * @return solution 
	 */ 
	Vector<COLS,T> solve( Vector<COLS,T> const & c=(Vector<COLS,T>((T)0)) ) const { 
	
		// This is an implementation of the LUP-SOLVE algorithm presented in 
		// Cohen, Leiserson, Rivest and Stein. "Introduction To Algorithms", 
		// 2nd Ed. MIT Press. 2001. 
		
		Vector<COLS,T> x; 
		
		// Forward substitution to solve Ly = Pc, where y = Ux 
		// Bit confusing, but we are actually storing y in x, to save 
		// creating a second vector! 
		for ( int i=0; i<ROWS; ++i ) { 
			T sum = 0; 
			for ( int j=0; j<i; ++j ) { 
				sum += _lu[i][j] * x[j]; 
			} 
			x[i] = c[_p[i]] - sum; 
		} 
		
		// Backwards substitution to solve Ux = y 
		for ( int i=ROWS-1; i>=0; --i ) { 
			T sum = 0; 
			for ( int j=i+1; j<ROWS; ++j ) { 
				sum += _lu[i][j] * x[j]; 
			} 
			x[i] = (x[i] - sum) / _lu[i][i]; 
		} 
		
		return x; 
	} 
	
private : 

	bool 	_failed; 	///< set iff LUP decomposition failed 
				///< (input matrix was singular) 
	Matrix 	_lu; 		///< LU matrices combined (since L is 
				///< unit lower-triangular, the 1's on it's diagonal 
				///< are implied) 
	int 	_p[ROWS]; 	///< P matrix converted into array form, 
				///< where P[i][j] = (j==_p[i]) 
};