H3cJP/source-engine
1
0
Fork 0
mirror of https://github.com/nillerusr/source-engine.git synced 2024-12-22 22:27:05 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity
efeda62add
source-engine/raytrace/raytrace.cpp

902 lines
28 KiB
C++
Raw Normal View History

1
2020-04-22 16:56:21 +00:00
//========= 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]);
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);
}
Reference in New Issue Copy Permalink