1 files changed, 901 insertions, 0 deletions

diff --git a/mp/src/raytrace/raytrace.cpp b/mp/src/raytrace/raytrace.cpp
new file mode 100644
index 00000000..7b951dfb
--- /dev/null
+++ b/mp/src/raytrace/raytrace.cpp
@@ -0,0 +1,901 @@
+//========= Copyright Valve Corporation, All rights reserved. ============//
+// $Id$
+
+#include "raytrace.h"
+#include <filesystem_tools.h>
+#include <cmdlib.h>
+#include <stdio.h>
+
+static bool SameSign(float a, float b)
+{
+	int32 aa=*((int *) &a);
+	int32 bb=*((int *) &b);
+	return ((aa^bb)&0x80000000)==0;
+}
+
+int FourRays::CalculateDirectionSignMask(void) const
+{
+	// this code treats the floats as integers since all it cares about is the sign bit and
+	// floating point compares suck.
+
+	int ret;
+	int ormask;
+	int andmask;
+	int32 const *treat_as_int=((int32 const *) (&direction));
+
+	ormask=andmask=*(treat_as_int++);
+	ormask|=*treat_as_int;
+	andmask&=*(treat_as_int++);
+	ormask|=*(treat_as_int);
+	andmask&=*(treat_as_int++);
+	ormask|=*(treat_as_int);
+	andmask&=*(treat_as_int++);
+	if (ormask>=0)
+		ret=0;
+	else
+	{
+		if (andmask<0)
+			ret=1;
+		else return -1;
+	}
+	ormask=andmask=*(treat_as_int++);
+	ormask|=*treat_as_int;
+	andmask&=*(treat_as_int++);
+	ormask|=*(treat_as_int);
+	andmask&=*(treat_as_int++);
+	ormask|=*(treat_as_int);
+	andmask&=*(treat_as_int++);
+	if (ormask<0)
+	{
+		if (andmask<0)
+			ret|=2;
+		else return -1;
+	}
+	ormask=andmask=*(treat_as_int++);
+	ormask|=*treat_as_int;
+	andmask&=*(treat_as_int++);
+	ormask|=*(treat_as_int);
+	andmask&=*(treat_as_int++);
+	ormask|=*(treat_as_int);
+	andmask&=*(treat_as_int++);
+	if (ormask<0)
+	{
+		if (andmask<0)
+			ret|=4;
+		else return -1;
+	}
+	return ret;
+}
+
+
+
+
+void RayTracingEnvironment::MakeRoomForTriangles( int ntris )
+{
+	//OptimizedTriangleList.EnsureCapacity( ntris );
+	if (! (Flags & RTE_FLAGS_DONT_STORE_TRIANGLE_COLORS))
+		TriangleColors.EnsureCapacity( ntris );
+}
+
+
+void RayTracingEnvironment::AddTriangle(int32 id, const Vector &v1,
+										const Vector &v2, const Vector &v3,
+										const Vector &color)
+{
+	AddTriangle( id, v1, v2, v3, color, 0, 0 );
+}
+
+void RayTracingEnvironment::AddTriangle(int32 id, const Vector &v1,
+										const Vector &v2, const Vector &v3,
+										const Vector &color, uint16 flags, int32 materialIndex)
+{
+	CacheOptimizedTriangle tmptri;
+	tmptri.m_Data.m_GeometryData.m_nTriangleID = id;
+	tmptri.Vertex( 0 ) = v1;
+	tmptri.Vertex( 1 ) = v2;
+	tmptri.Vertex( 2 ) = v3;
+	tmptri.m_Data.m_GeometryData.m_nFlags = flags;
+	OptimizedTriangleList.AddToTail(tmptri);
+	if (! ( Flags & RTE_FLAGS_DONT_STORE_TRIANGLE_COLORS) )
+		TriangleColors.AddToTail(color);
+	if ( !( Flags & RTE_FLAGS_DONT_STORE_TRIANGLE_MATERIALS) )
+		TriangleMaterials.AddToTail(materialIndex);
+// 	printf("add triange from (%f %f %f),(%f %f %f),(%f %f %f) id %d\n",
+// 		   XYZ(v1),XYZ(v2),XYZ(v3),id);
+}
+
+void RayTracingEnvironment::AddQuad(
+	int32 id, const Vector &v1, const Vector &v2, const Vector &v3,
+	const Vector &v4,							// specify vertices in cw or ccw order
+	const Vector &color)
+{
+	AddTriangle(id,v1,v2,v3,color);
+	AddTriangle(id+1,v1,v3,v4,color);
+}
+
+
+void RayTracingEnvironment::AddAxisAlignedRectangularSolid(int id,Vector minc, Vector maxc,
+														   const Vector &color)
+{
+
+	// "far" face
+	AddQuad(id,
+			Vector(minc.x,maxc.y,maxc.z),
+			Vector(maxc.x,maxc.y,maxc.z),Vector(maxc.x,minc.y,maxc.z),
+			Vector(minc.x,minc.y,maxc.z),color);
+	// "near" face
+	AddQuad(id,
+			Vector(minc.x,maxc.y,minc.z),
+			Vector(maxc.x,maxc.y,minc.z),Vector(maxc.x,minc.y,minc.z),
+			Vector(minc.x,minc.y,minc.z),color);
+
+	// "left" face
+	AddQuad(id,
+			Vector(minc.x,maxc.y,maxc.z),
+			Vector(minc.x,maxc.y,minc.z),
+			Vector(minc.x,minc.y,minc.z),
+			Vector(minc.x,minc.y,maxc.z),color);
+	// "right" face
+	AddQuad(id,
+			Vector(maxc.x,maxc.y,maxc.z),
+			Vector(maxc.x,maxc.y,minc.z),
+			Vector(maxc.x,minc.y,minc.z),
+			Vector(maxc.x,minc.y,maxc.z),color);
+	
+	// "top" face
+	AddQuad(id,
+			Vector(minc.x,maxc.y,maxc.z),
+			Vector(maxc.x,maxc.y,maxc.z),
+			Vector(maxc.x,maxc.y,minc.z),
+			Vector(minc.x,maxc.y,minc.z),color);
+	// "bot" face
+	AddQuad(id,
+			Vector(minc.x,minc.y,maxc.z),
+			Vector(maxc.x,minc.y,maxc.z),
+			Vector(maxc.x,minc.y,minc.z),
+			Vector(minc.x,minc.y,minc.z),color);
+}
+
+
+
+static Vector GetEdgeEquation(Vector p1, Vector p2, int c1, int c2, Vector InsidePoint)
+{
+	float nx=p1[c2]-p2[c2];
+	float ny=p2[c1]-p1[c1];
+	float d=-(nx*p1[c1]+ny*p1[c2]);
+// 	assert(fabs(nx*p1[c1]+ny*p1[c2]+d)<0.01);
+// 	assert(fabs(nx*p2[c1]+ny*p2[c2]+d)<0.01);
+
+	// use the convention that negative is "outside"
+	float trial_dist=InsidePoint[c1]*nx+InsidePoint[c2]*ny+d;
+	if (trial_dist<0)
+	{
+		nx = -nx;
+		ny = -ny;
+		d = -d;
+		trial_dist = -trial_dist;
+	}
+	nx /= trial_dist;										// scale so that it will be =1.0 at the oppositve vertex
+	ny /= trial_dist;
+	d /= trial_dist;
+
+	return Vector(nx,ny,d);
+}
+
+void CacheOptimizedTriangle::ChangeIntoIntersectionFormat(void)
+{
+	// lose the vertices and use edge equations instead
+
+	// grab the whole original triangle to we don't overwrite it
+	TriGeometryData_t srcTri = m_Data.m_GeometryData;
+
+	m_Data.m_IntersectData.m_nFlags = srcTri.m_nFlags;
+	m_Data.m_IntersectData.m_nTriangleID = srcTri.m_nTriangleID;
+
+	Vector p1 = srcTri.Vertex( 0 );
+	Vector p2 = srcTri.Vertex( 1 );
+	Vector p3 = srcTri.Vertex( 2 );
+
+	Vector e1 = p2 - p1;
+	Vector e2 = p3 - p1;
+
+	Vector N = e1.Cross( e2 );
+	N.NormalizeInPlace();
+	// now, determine which axis to drop
+	int drop_axis = 0;
+	for(int c=1 ; c<3 ; c++)
+		if ( fabs(N[c]) > fabs( N[drop_axis] ) )
+			drop_axis = c;
+
+	m_Data.m_IntersectData.m_flD = N.Dot( p1 );
+	m_Data.m_IntersectData.m_flNx = N.x;
+	m_Data.m_IntersectData.m_flNy = N.y;
+	m_Data.m_IntersectData.m_flNz = N.z;
+
+	// decide which axes to keep
+	int nCoordSelect0 = ( drop_axis + 1 ) % 3;
+	int nCoordSelect1 = ( drop_axis + 2 ) % 3;
+
+	m_Data.m_IntersectData.m_nCoordSelect0 = nCoordSelect0;
+	m_Data.m_IntersectData.m_nCoordSelect1 = nCoordSelect1;
+
+
+	Vector edge1 = GetEdgeEquation( p1, p2, nCoordSelect0, nCoordSelect1, p3 );
+	m_Data.m_IntersectData.m_ProjectedEdgeEquations[0] = edge1.x;
+	m_Data.m_IntersectData.m_ProjectedEdgeEquations[1] = edge1.y;
+	m_Data.m_IntersectData.m_ProjectedEdgeEquations[2] = edge1.z;
+
+	Vector edge2 = GetEdgeEquation( p2, p3, nCoordSelect0, nCoordSelect1, p1 );
+	m_Data.m_IntersectData.m_ProjectedEdgeEquations[3] = edge2.x;
+	m_Data.m_IntersectData.m_ProjectedEdgeEquations[4] = edge2.y;
+	m_Data.m_IntersectData.m_ProjectedEdgeEquations[5] = edge2.z;
+
+
+}
+
+int n_intersection_calculations=0;
+
+int CacheOptimizedTriangle::ClassifyAgainstAxisSplit(int split_plane, float split_value)
+{
+	// classify a triangle against an axis-aligned plane
+	float minc=Vertex(0)[split_plane];
+	float maxc=minc;
+	for(int v=1;v<3;v++)
+	{
+		minc=min(minc,Vertex(v)[split_plane]);
+		maxc=max(maxc,Vertex(v)[split_plane]);
+	}
+
+	if (minc>=split_value)
+		return PLANECHECK_POSITIVE;
+	if (maxc<=split_value)
+		return PLANECHECK_NEGATIVE;
+	if (minc==maxc)
+		return PLANECHECK_POSITIVE;
+	return PLANECHECK_STRADDLING;
+}
+
+#define MAILBOX_HASH_SIZE 256
+#define MAX_TREE_DEPTH 21
+#define MAX_NODE_STACK_LEN (40*MAX_TREE_DEPTH)
+
+struct NodeToVisit {
+	CacheOptimizedKDNode const *node;
+	fltx4 TMin;
+	fltx4 TMax;
+};
+
+
+static fltx4 FourEpsilons={1.0e-10,1.0e-10,1.0e-10,1.0e-10};
+static fltx4 FourZeros={1.0e-10,1.0e-10,1.0e-10,1.0e-10};
+static fltx4 FourNegativeEpsilons={-1.0e-10,-1.0e-10,-1.0e-10,-1.0e-10};
+
+static float BoxSurfaceArea(Vector const &boxmin, Vector const &boxmax)
+{
+	Vector boxdim=boxmax-boxmin;
+	return 2.0*((boxdim[0]*boxdim[2])+(boxdim[0]*boxdim[1])+(boxdim[1]*boxdim[2]));
+}
+
+void RayTracingEnvironment::Trace4Rays(const FourRays &rays, fltx4 TMin, fltx4 TMax,
+									   RayTracingResult *rslt_out,
+									   int32 skip_id, ITransparentTriangleCallback *pCallback)
+{
+	int msk=rays.CalculateDirectionSignMask();
+	if (msk!=-1)
+		Trace4Rays(rays,TMin,TMax,msk,rslt_out,skip_id, pCallback);
+	else
+	{
+		// sucky case - can't trace 4 rays at once. in the worst case, need to trace all 4
+		// separately, but usually we will still get 2x, Since our tracer only does 4 at a
+		// time, we will have to cover up the undesired rays with the desired ray
+
+		//!! speed!! there is room for some sse-ization here
+		FourRays tmprays;
+		tmprays.origin=rays.origin;
+
+		uint8 need_trace[4]={1,1,1,1};
+		for(int try_trace=0;try_trace<4;try_trace++)
+		{
+			if (need_trace[try_trace])
+			{
+				need_trace[try_trace]=2;			// going to trace it
+				// replicate the ray being traced into all 4 rays
+				tmprays.direction.x=ReplicateX4(rays.direction.X(try_trace));
+				tmprays.direction.y=ReplicateX4(rays.direction.Y(try_trace));
+				tmprays.direction.z=ReplicateX4(rays.direction.Z(try_trace));
+				// now, see if any of the other remaining rays can be handled at the same time.
+				for(int try2=try_trace+1;try2<4;try2++)
+					if (need_trace[try2])
+					{
+						if (
+							SameSign(rays.direction.X(try2),
+									 rays.direction.X(try_trace)) &&
+							SameSign(rays.direction.Y(try2),
+									 rays.direction.Y(try_trace)) &&
+							SameSign(rays.direction.Z(try2),
+									 rays.direction.Z(try_trace)))
+						{
+							need_trace[try2]=2;
+							tmprays.direction.X(try2) = rays.direction.X(try2);
+							tmprays.direction.Y(try2) = rays.direction.Y(try2);
+							tmprays.direction.Z(try2) = rays.direction.Z(try2);
+						}
+					}
+				// ok, now trace between 1 and 3 rays, and output the results
+				RayTracingResult tmpresults;
+				msk=tmprays.CalculateDirectionSignMask();
+				assert(msk!=-1);
+				Trace4Rays(tmprays,TMin,TMax,msk,&tmpresults,skip_id, pCallback);
+				// now, move results to proper place
+				for(int i=0;i<4;i++)
+					if (need_trace[i]==2)
+					{
+						need_trace[i]=0;
+						rslt_out->HitIds[i]=tmpresults.HitIds[i];
+						SubFloat(rslt_out->HitDistance, i) = SubFloat(tmpresults.HitDistance, i);
+						rslt_out->surface_normal.X(i) = tmpresults.surface_normal.X(i);
+						rslt_out->surface_normal.Y(i) = tmpresults.surface_normal.Y(i);
+						rslt_out->surface_normal.Z(i) = tmpresults.surface_normal.Z(i);
+					}
+				
+			}
+		}
+	}
+}
+
+
+void RayTracingEnvironment::Trace4Rays(const FourRays &rays, fltx4 TMin, fltx4 TMax,
+									   int DirectionSignMask, RayTracingResult *rslt_out,
+									   int32 skip_id, ITransparentTriangleCallback *pCallback)
+{
+	rays.Check();
+
+	memset(rslt_out->HitIds,0xff,sizeof(rslt_out->HitIds));
+
+	rslt_out->HitDistance=ReplicateX4(1.0e23);
+
+	rslt_out->surface_normal.DuplicateVector(Vector(0.,0.,0.));
+	FourVectors OneOverRayDir=rays.direction;
+	OneOverRayDir.MakeReciprocalSaturate();
+	
+	// now, clip rays against bounding box
+	for(int c=0;c<3;c++)
+	{
+		fltx4 isect_min_t=
+			MulSIMD(SubSIMD(ReplicateX4(m_MinBound[c]),rays.origin[c]),OneOverRayDir[c]);
+		fltx4 isect_max_t=
+			MulSIMD(SubSIMD(ReplicateX4(m_MaxBound[c]),rays.origin[c]),OneOverRayDir[c]);
+		TMin=MaxSIMD(TMin,MinSIMD(isect_min_t,isect_max_t));
+		TMax=MinSIMD(TMax,MaxSIMD(isect_min_t,isect_max_t));
+	}
+	fltx4 active=CmpLeSIMD(TMin,TMax);					// mask of which rays are active
+	if (! IsAnyNegative(active) )
+		return;												// missed bounding box
+
+	int32 mailboxids[MAILBOX_HASH_SIZE];					// used to avoid redundant triangle tests
+	memset(mailboxids,0xff,sizeof(mailboxids));				// !!speed!! keep around?
+
+	int front_idx[3],back_idx[3];							// based on ray direction, whether to
+															// visit left or right node first
+
+	if (DirectionSignMask & 1)
+	{
+		back_idx[0]=0;
+		front_idx[0]=1;
+	}
+		else
+	{
+		back_idx[0]=1;
+		front_idx[0]=0;
+	}
+	if (DirectionSignMask & 2)
+	{
+		back_idx[1]=0;
+		front_idx[1]=1;
+	}
+	else
+	{
+		back_idx[1]=1;
+		front_idx[1]=0;
+	}
+	if (DirectionSignMask & 4)
+	{
+		back_idx[2]=0;
+		front_idx[2]=1;
+	}
+	else
+	{
+		back_idx[2]=1;
+		front_idx[2]=0;
+	}
+		
+	NodeToVisit NodeQueue[MAX_NODE_STACK_LEN];
+	CacheOptimizedKDNode const *CurNode=&(OptimizedKDTree[0]);
+	NodeToVisit *stack_ptr=&NodeQueue[MAX_NODE_STACK_LEN];
+	while(1)
+	{
+		while (CurNode->NodeType() != KDNODE_STATE_LEAF)		// traverse until next leaf
+		{	   
+			int split_plane_number=CurNode->NodeType();
+			CacheOptimizedKDNode const *FrontChild=&(OptimizedKDTree[CurNode->LeftChild()]);
+			
+			fltx4 dist_to_sep_plane=						// dist=(split-org)/dir
+				MulSIMD(
+					SubSIMD(ReplicateX4(CurNode->SplittingPlaneValue),
+							   rays.origin[split_plane_number]),OneOverRayDir[split_plane_number]);
+			fltx4 active=CmpLeSIMD(TMin,TMax);			// mask of which rays are active
+
+			// now, decide how to traverse children. can either do front,back, or do front and push
+			// back.
+			fltx4 hits_front=AndSIMD(active,CmpGeSIMD(dist_to_sep_plane,TMin));
+			if (! IsAnyNegative(hits_front))
+			{
+				// missed the front. only traverse back
+				//printf("only visit back %d\n",CurNode->LeftChild()+back_idx[split_plane_number]);
+				CurNode=FrontChild+back_idx[split_plane_number];
+				TMin=MaxSIMD(TMin, dist_to_sep_plane);
+
+			}
+			else
+			{
+				fltx4 hits_back=AndSIMD(active,CmpLeSIMD(dist_to_sep_plane,TMax));
+				if (! IsAnyNegative(hits_back) )
+				{
+					// missed the back - only need to traverse front node
+					//printf("only visit front %d\n",CurNode->LeftChild()+front_idx[split_plane_number]);
+					CurNode=FrontChild+front_idx[split_plane_number];
+					TMax=MinSIMD(TMax, dist_to_sep_plane);
+				}
+				else
+				{
+					// at least some rays hit both nodes.
+					// must push far, traverse near
+ 					//printf("visit %d,%d\n",CurNode->LeftChild()+front_idx[split_plane_number],
+ 					//	   CurNode->LeftChild()+back_idx[split_plane_number]);
+					assert(stack_ptr>NodeQueue);
+					--stack_ptr;
+					stack_ptr->node=FrontChild+back_idx[split_plane_number];
+					stack_ptr->TMin=MaxSIMD(TMin,dist_to_sep_plane);
+					stack_ptr->TMax=TMax;
+					CurNode=FrontChild+front_idx[split_plane_number];
+					TMax=MinSIMD(TMax,dist_to_sep_plane);
+				}
+			}
+		}
+		// hit a leaf! must do intersection check
+		int ntris=CurNode->NumberOfTrianglesInLeaf();
+		if (ntris)
+		{
+			int32 const *tlist=&(TriangleIndexList[CurNode->TriangleIndexStart()]);
+			do
+			{
+				int tnum=*(tlist++);
+				//printf("try tri %d\n",tnum);
+				// check mailbox
+				int mbox_slot=tnum & (MAILBOX_HASH_SIZE-1);
+				TriIntersectData_t const *tri = &( OptimizedTriangleList[tnum].m_Data.m_IntersectData );
+				if ( ( mailboxids[mbox_slot] != tnum ) && ( tri->m_nTriangleID != skip_id ) )
+				{
+					n_intersection_calculations++;
+					mailboxids[mbox_slot] = tnum;
+					// compute plane intersection
+
+
+					FourVectors N;
+					N.x = ReplicateX4( tri->m_flNx );
+					N.y = ReplicateX4( tri->m_flNy );
+					N.z = ReplicateX4( tri->m_flNz );
+
+					fltx4 DDotN = rays.direction * N;
+					// mask off zero or near zero (ray parallel to surface)
+					fltx4 did_hit = OrSIMD( CmpGtSIMD( DDotN,FourEpsilons ),
+											CmpLtSIMD( DDotN, FourNegativeEpsilons ) );
+
+					fltx4 numerator=SubSIMD( ReplicateX4( tri->m_flD ), rays.origin * N );
+
+					fltx4 isect_t=DivSIMD( numerator,DDotN );
+					// now, we have the distance to the plane. lets update our mask
+					did_hit = AndSIMD( did_hit, CmpGtSIMD( isect_t, FourZeros ) );
+					//did_hit=AndSIMD(did_hit,CmpLtSIMD(isect_t,TMax));
+					did_hit = AndSIMD( did_hit, CmpLtSIMD( isect_t, rslt_out->HitDistance ) );
+
+					if ( ! IsAnyNegative( did_hit ) )
+						continue;
+
+					// now, check 3 edges
+					fltx4 hitc1 = AddSIMD( rays.origin[tri->m_nCoordSelect0],
+										MulSIMD( isect_t, rays.direction[ tri->m_nCoordSelect0] ) );
+					fltx4 hitc2 = AddSIMD( rays.origin[tri->m_nCoordSelect1],
+										   MulSIMD( isect_t, rays.direction[tri->m_nCoordSelect1] ) );
+					
+					// do barycentric coordinate check
+					fltx4 B0 = MulSIMD( ReplicateX4( tri->m_ProjectedEdgeEquations[0] ), hitc1 );
+
+					B0 = AddSIMD(
+						B0,
+						MulSIMD( ReplicateX4( tri->m_ProjectedEdgeEquations[1] ), hitc2 ) );
+					B0 = AddSIMD(
+						B0, ReplicateX4( tri->m_ProjectedEdgeEquations[2] ) );
+
+					did_hit = AndSIMD( did_hit, CmpGeSIMD( B0, FourZeros ) );
+
+					fltx4 B1 = MulSIMD( ReplicateX4( tri->m_ProjectedEdgeEquations[3] ), hitc1 );
+					B1 = AddSIMD(
+						B1,
+						MulSIMD( ReplicateX4( tri->m_ProjectedEdgeEquations[4]), hitc2 ) );
+
+					B1 = AddSIMD(
+						B1, ReplicateX4( tri->m_ProjectedEdgeEquations[5] ) );
+					
+					did_hit = AndSIMD( did_hit, CmpGeSIMD( B1, FourZeros ) );
+
+					fltx4 B2 = AddSIMD( B1, B0 );
+					did_hit = AndSIMD( did_hit, CmpLeSIMD( B2, Four_Ones ) );
+
+					if ( ! IsAnyNegative( did_hit ) )
+						continue;
+
+					// if the triangle is transparent
+					if ( tri->m_nFlags & FCACHETRI_TRANSPARENT )
+					{
+						if ( pCallback )
+						{
+							// assuming a triangle indexed as v0, v1, v2
+							// the projected edge equations are set up such that the vert opposite the first
+							// equation is v2, and the vert opposite the second equation is v0
+							// Therefore we pass them back in 1, 2, 0 order
+							// Also B2 is currently B1 + B0 and needs to be 1 - (B1+B0) in order to be a real
+							// barycentric coordinate.  Compute that now and pass it to the callback
+							fltx4 b2 = SubSIMD( Four_Ones, B2 );
+							if ( pCallback->VisitTriangle_ShouldContinue( *tri, rays, &did_hit, &B1, &b2, &B0, tnum ) )
+							{
+								did_hit = Four_Zeros;
+							}
+						}
+					}
+					// now, set the hit_id and closest_hit fields for any enabled rays
+					fltx4 replicated_n = ReplicateIX4(tnum);
+					StoreAlignedSIMD((float *) rslt_out->HitIds,
+								 OrSIMD(AndSIMD(replicated_n,did_hit),
+										   AndNotSIMD(did_hit,LoadAlignedSIMD(
+															 (float *) rslt_out->HitIds))));
+					rslt_out->HitDistance=OrSIMD(AndSIMD(isect_t,did_hit),
+									 AndNotSIMD(did_hit,rslt_out->HitDistance));
+
+					rslt_out->surface_normal.x=OrSIMD(
+						AndSIMD(N.x,did_hit),
+						AndNotSIMD(did_hit,rslt_out->surface_normal.x));
+					rslt_out->surface_normal.y=OrSIMD(
+						AndSIMD(N.y,did_hit),
+						AndNotSIMD(did_hit,rslt_out->surface_normal.y));
+					rslt_out->surface_normal.z=OrSIMD(
+						AndSIMD(N.z,did_hit),
+						AndNotSIMD(did_hit,rslt_out->surface_normal.z));
+					
+				}
+			} while (--ntris);
+			// now, check if all rays have terminated
+			fltx4 raydone=CmpLeSIMD(TMax,rslt_out->HitDistance);
+			if (! IsAnyNegative(raydone))
+			{
+				return;
+			}
+		}
+		
+ 		if (stack_ptr==&NodeQueue[MAX_NODE_STACK_LEN])
+		{
+			return;
+		}
+		// pop stack!
+		CurNode=stack_ptr->node;
+		TMin=stack_ptr->TMin;
+		TMax=stack_ptr->TMax;
+		stack_ptr++;
+	}
+}
+
+
+int RayTracingEnvironment::MakeLeafNode(int first_tri, int last_tri)
+{
+	CacheOptimizedKDNode ret;
+	ret.Children=KDNODE_STATE_LEAF+(TriangleIndexList.Count()<<2);
+	ret.SetNumberOfTrianglesInLeafNode(1+(last_tri-first_tri));
+	for(int tnum=first_tri;tnum<=last_tri;tnum++)
+		TriangleIndexList.AddToTail(tnum);
+	OptimizedKDTree.AddToTail(ret);
+	return OptimizedKDTree.Count()-1;
+}
+
+
+void RayTracingEnvironment::CalculateTriangleListBounds(int32 const *tris,int ntris,
+														Vector &minout, Vector &maxout)
+{
+	minout = Vector( 1.0e23, 1.0e23, 1.0e23);
+	maxout = Vector( -1.0e23, -1.0e23, -1.0e23);
+	for(int i=0; i<ntris; i++)
+	{
+		CacheOptimizedTriangle const &tri=OptimizedTriangleList[tris[i]];
+		for(int v=0; v<3; v++)
+			for(int c=0; c<3; c++)
+			{
+				minout[c]=min(minout[c],tri.Vertex(v)[c]);
+							  maxout[c]=max(maxout[c],tri.Vertex(v)[c]);
+			}
+	}
+}
+
+
+// Both the "quick" and regular kd tree building algorithms here use the "surface area heuristic":
+// the relative probability of hitting the "left" subvolume (Vl) from a split is equal to that
+// subvolume's surface area divided by its parent's surface area (Vp) : P(Vl | V)=SA(Vl)/SA(Vp).
+// The same holds for the right subvolume, Vp. Nl is the number of triangles in the left volume,
+// and Nr in the right volume. if Ct is the cost of traversing one tree node, and Ci is the cost of
+// intersection with the primitive, than the cost of splitting is estimated as:
+//
+//    Ct+Ci*((SA(Vl)/SA(V))*Nl+(SA(Vr)/SA(V)*Nr)).
+// and the cost of not splitting is
+//    Ci*N
+//
+//  This both provides a metric to minimize when computing how and where to split, and also a
+//  termination criterion.
+//
+// the "quick" method just splits down the middle, while the slow method splits at the best
+// discontinuity of the cost formula. The quick method splits along the longest axis ; the 
+// regular algorithm tries all 3 to find which one results in the minimum cost
+// 
+// both methods use the additional optimization of "growing" empty nodes - if the split results in
+// one side being devoid of triangles, the empty side is "grown" as much as possible.
+//
+
+#define COST_OF_TRAVERSAL 75								// approximate #operations
+#define COST_OF_INTERSECTION 167							// approximate #operations
+
+
+float RayTracingEnvironment::CalculateCostsOfSplit(
+	int split_plane,int32 const *tri_list,int ntris,
+	Vector MinBound,Vector MaxBound, float &split_value,
+	int &nleft, int &nright, int &nboth)
+{
+	// determine the costs of splitting on a given axis, and label triangles with respect to
+	// that axis by storing the value in coordselect0. It will also return the number of
+	// tris in the left, right, and nboth groups, in order to facilitate memory
+	nleft=nboth=nright=0;
+	
+	// now, label each triangle. Since we have not converted the triangles into
+	// intersection fromat yet, we can use the CoordSelect0 field of each as a temp.
+	nleft=0;
+	nright=0;
+	nboth=0;
+	float min_coord=1.0e23,max_coord=-1.0e23;
+
+	for(int t=0;t<ntris;t++)
+	{
+		CacheOptimizedTriangle &tri=OptimizedTriangleList[tri_list[t]];
+		// determine max and min coordinate values for later optimization
+		for(int v=0;v<3;v++)
+		{
+			min_coord = min( min_coord, tri.Vertex(v)[split_plane] );
+			max_coord = max( max_coord, tri.Vertex(v)[split_plane] );
+		}
+		switch(tri.ClassifyAgainstAxisSplit(split_plane,split_value))
+		{
+			case PLANECHECK_NEGATIVE:
+				nleft++;
+				tri.m_Data.m_GeometryData.m_nTmpData0 = PLANECHECK_NEGATIVE;
+				break;
+
+			case PLANECHECK_POSITIVE:
+				nright++;
+				tri.m_Data.m_GeometryData.m_nTmpData0 = PLANECHECK_POSITIVE;
+				break;
+
+			case PLANECHECK_STRADDLING:
+				nboth++;
+				tri.m_Data.m_GeometryData.m_nTmpData0 = PLANECHECK_STRADDLING;
+				break;
+		}
+	}
+	// now, if the split resulted in one half being empty, "grow" the empty half
+	if (nleft && (nboth==0) && (nright==0))
+		split_value=max_coord;
+	if (nright && (nboth==0) && (nleft==0))
+		split_value=min_coord;
+
+	// now, perform surface area/cost check to determine whether this split was worth it
+	Vector LeftMins=MinBound;
+	Vector LeftMaxes=MaxBound;
+	Vector RightMins=MinBound;
+	Vector RightMaxes=MaxBound;
+	LeftMaxes[split_plane]=split_value;
+	RightMins[split_plane]=split_value;
+	float SA_L=BoxSurfaceArea(LeftMins,LeftMaxes);
+	float SA_R=BoxSurfaceArea(RightMins,RightMaxes);
+	float ISA=1.0/BoxSurfaceArea(MinBound,MaxBound);
+	float cost_of_split=COST_OF_TRAVERSAL+COST_OF_INTERSECTION*(nboth+
+		(SA_L*ISA*(nleft))+(SA_R*ISA*(nright)));
+	return cost_of_split;
+}
+
+
+#define NEVER_SPLIT 0
+
+void RayTracingEnvironment::RefineNode(int node_number,int32 const *tri_list,int ntris,
+									   Vector MinBound,Vector MaxBound, int depth)
+{
+	if (ntris<3)											// never split empty lists
+	{
+		// no point in continuing
+		OptimizedKDTree[node_number].Children=KDNODE_STATE_LEAF+(TriangleIndexList.Count()<<2);
+		OptimizedKDTree[node_number].SetNumberOfTrianglesInLeafNode(ntris);
+
+#ifdef DEBUG_RAYTRACE
+		OptimizedKDTree[node_number].vecMins = MinBound;
+		OptimizedKDTree[node_number].vecMaxs = MaxBound;
+#endif
+
+		for(int t=0;t<ntris;t++)
+			TriangleIndexList.AddToTail(tri_list[t]);
+		return;
+	}
+
+	float best_cost=1.0e23;
+	int best_nleft=0,best_nright=0,best_nboth=0;
+	float best_splitvalue=0;
+	int split_plane=0;
+
+	int tri_skip=1+(ntris/10);								// don't try all trinagles as split
+															// points when there are a lot of them
+	for(int axis=0;axis<3;axis++)
+	{
+		for(int ts=-1;ts<ntris;ts+=tri_skip)
+		{
+			for(int tv=0;tv<3;tv++)
+			{
+				int trial_nleft,trial_nright,trial_nboth;
+				float trial_splitvalue;
+				if (ts==-1)
+					trial_splitvalue=0.5*(MinBound[axis]+MaxBound[axis]);
+				else
+				{
+					// else, split at the triangle vertex if possible
+					CacheOptimizedTriangle &tri=OptimizedTriangleList[tri_list[ts]];
+					trial_splitvalue = tri.Vertex(tv)[axis];
+					if ((trial_splitvalue>MaxBound[axis]) || (trial_splitvalue<MinBound[axis]))
+						continue;							// don't try this vertex - not inside
+					
+				}
+//				printf("ts=%d tv=%d tp=%f\n",ts,tv,trial_splitvalue);
+				float trial_cost=
+					CalculateCostsOfSplit(axis,tri_list,ntris,MinBound,MaxBound,trial_splitvalue,
+										  trial_nleft,trial_nright, trial_nboth);
+// 				printf("try %d cost=%f nl=%d nr=%d nb=%d sp=%f\n",axis,trial_cost,trial_nleft,trial_nright, trial_nboth,
+// 					   trial_splitvalue);
+				if (trial_cost<best_cost)
+				{
+					split_plane=axis;
+					best_cost=trial_cost;
+					best_nleft=trial_nleft;
+					best_nright=trial_nright;
+					best_nboth=trial_nboth;
+					best_splitvalue=trial_splitvalue;
+					// save away the axis classification of each triangle
+					for(int t=0 ; t < ntris; t++)
+					{
+						CacheOptimizedTriangle &tri=OptimizedTriangleList[tri_list[t]];
+						tri.m_Data.m_GeometryData.m_nTmpData1 = tri.m_Data.m_GeometryData.m_nTmpData0;
+					}
+				}
+				if (ts==-1)
+					break;
+			}
+		}
+
+	}
+	float cost_of_no_split=COST_OF_INTERSECTION*ntris;
+	if ( (cost_of_no_split<=best_cost) || NEVER_SPLIT || (depth>MAX_TREE_DEPTH))
+	{
+		// no benefit to splitting. just make this a leaf node
+		OptimizedKDTree[node_number].Children=KDNODE_STATE_LEAF+(TriangleIndexList.Count()<<2);
+		OptimizedKDTree[node_number].SetNumberOfTrianglesInLeafNode(ntris);
+#ifdef DEBUG_RAYTRACE
+		OptimizedKDTree[node_number].vecMins = MinBound;
+		OptimizedKDTree[node_number].vecMaxs = MaxBound;
+#endif
+		for(int t=0;t<ntris;t++)
+			TriangleIndexList.AddToTail(tri_list[t]);
+	}
+	else
+	{
+// 		printf("best split was %d at %f (mid=%f,n=%d, sk=%d)\n",split_plane,best_splitvalue,
+// 			   0.5*(MinBound[split_plane]+MaxBound[split_plane]),ntris,tri_skip);
+		// its worth splitting!
+		// we will achieve the splitting without sorting by using a selection algorithm.
+		int32 *new_triangle_list;
+		new_triangle_list=new int32[ntris];
+
+		// now, perform surface area/cost check to determine whether this split was worth it
+		Vector LeftMins=MinBound;
+		Vector LeftMaxes=MaxBound;
+		Vector RightMins=MinBound;
+		Vector RightMaxes=MaxBound;
+		LeftMaxes[split_plane]=best_splitvalue;
+		RightMins[split_plane]=best_splitvalue;
+		
+		int n_left_output=0;
+		int n_both_output=0;
+		int n_right_output=0;
+		for(int t=0;t<ntris;t++)
+		{
+			CacheOptimizedTriangle &tri=OptimizedTriangleList[tri_list[t]];
+			switch( tri.m_Data.m_GeometryData.m_nTmpData1 )
+			{
+				case PLANECHECK_NEGATIVE:
+//					printf("%d goes left\n",t);
+					new_triangle_list[n_left_output++]=tri_list[t];
+					break;
+				case PLANECHECK_POSITIVE:
+					n_right_output++;
+//					printf("%d goes right\n",t);
+					new_triangle_list[ntris-n_right_output]=tri_list[t];
+					break;
+				case PLANECHECK_STRADDLING:
+//					printf("%d goes both\n",t);
+					new_triangle_list[best_nleft+n_both_output]=tri_list[t];
+					n_both_output++;
+					break;
+
+					
+			}
+		}
+		int left_child=OptimizedKDTree.Count();
+		int right_child=left_child+1;
+// 		printf("node %d split on axis %d at %f, nl=%d nr=%d nb=%d lc=%d rc=%d\n",node_number,
+// 			   split_plane,best_splitvalue,best_nleft,best_nright,best_nboth,
+// 			   left_child,right_child);
+		OptimizedKDTree[node_number].Children=split_plane+(left_child<<2);
+		OptimizedKDTree[node_number].SplittingPlaneValue=best_splitvalue;
+#ifdef DEBUG_RAYTRACE
+		OptimizedKDTree[node_number].vecMins = MinBound;
+		OptimizedKDTree[node_number].vecMaxs = MaxBound;
+#endif
+		CacheOptimizedKDNode newnode;
+		OptimizedKDTree.AddToTail(newnode);
+		OptimizedKDTree.AddToTail(newnode);
+		// now, recurse!
+		if ( (ntris<20) && ((best_nleft==0) || (best_nright==0)) )
+			depth+=100;
+		RefineNode(left_child,new_triangle_list,best_nleft+best_nboth,LeftMins,LeftMaxes,depth+1);
+		RefineNode(right_child,new_triangle_list+best_nleft,best_nright+best_nboth,
+				   RightMins,RightMaxes,depth+1);
+		delete[] new_triangle_list;
+	}	
+}
+
+
+void RayTracingEnvironment::SetupAccelerationStructure(void)
+{
+	CacheOptimizedKDNode root;
+	OptimizedKDTree.AddToTail(root);
+	int32 *root_triangle_list=new int32[OptimizedTriangleList.Count()];
+	for(int t=0;t<OptimizedTriangleList.Count();t++)
+		root_triangle_list[t]=t;
+	CalculateTriangleListBounds(root_triangle_list,OptimizedTriangleList.Count(),m_MinBound,
+								m_MaxBound);
+	RefineNode(0,root_triangle_list,OptimizedTriangleList.Count(),m_MinBound,m_MaxBound,0);
+	delete[] root_triangle_list;
+
+	// now, convert all triangles to "intersection format"
+	for(int i=0;i<OptimizedTriangleList.Count();i++)
+		OptimizedTriangleList[i].ChangeIntoIntersectionFormat();
+}
+
+
+
+void RayTracingEnvironment::AddInfinitePointLight(Vector position, Vector intensity)
+{
+	LightDesc_t mylight(position,intensity);
+	LightList.AddToTail(mylight);
+	
+}
+
+