Paranoia2/utils/p2rad/raytracer.cpp

575 lines
16 KiB
C++

/***
*
* Copyright (c) 1996-2002, Valve LLC. All rights reserved.
*
* This product contains software technology licensed from Id
* Software, Inc. ("Id Technology"). Id Technology (c) 1996 Id Software, Inc.
* All Rights Reserved.
*
****/
// raytracer.cpp
#include "qrad.h"
#include "raytracer.h"
#include "..\..\engine\studio.h"
#include "model_trace.h"
#ifdef HLRAD_RAYTRACE
void CWorldRayTrace :: AddTriangle( tface_t *tf )
{
if( !tf->shadow ) return;
m_TriangleList.AddToTail( tf );
}
float CWorldRayTrace :: BoxSurfaceArea( const vec3_t boxmin, const vec3_t boxmax )
{
vec3_t boxdim;
vec_t area;
VectorSubtract( boxmax, boxmin, boxdim );
area = 2.0 * ((boxdim[0] * boxdim[2]) + (boxdim[0] * boxdim[1]) + (boxdim[1] * boxdim[2]));
if( area == 0.0f ) area = FLT_EPSILON; // to prevent INF
assert( area != 0.0 );
return area;
}
int CWorldRayTrace :: MakeLeafNode( int first_tri, int last_tri )
{
KDNode ret;
ret.m_iChild = KDNODE_STATE_LEAF + (m_TriangleIndexList.Count() << 2);
ret.SetNumberOfTrianglesInLeafNode( 1 + (last_tri - first_tri));
for( int tnum = first_tri; tnum <= last_tri; tnum++ )
m_TriangleIndexList.AddToTail( tnum );
m_KDTree.AddToTail( ret );
return m_KDTree.Count() - 1;
}
void CWorldRayTrace :: CalculateTriangleListBounds( const int *tris, int ntris, vec3_t minout, vec3_t maxout )
{
VectorSet( minout, 1.0e23, 1.0e23, 1.0e23 );
VectorSet( maxout, -1.0e23, -1.0e23, -1.0e23 );
for( int i = 0; i < ntris; i++ )
{
const tface_t *face = m_TriangleList[tris[i]];
for( int c = 0; c < 3; c++ )
{
minout[c] = Q_min( minout[c], mesh->verts[face->a].point[c] );
maxout[c] = Q_max( maxout[c], mesh->verts[face->a].point[c] );
minout[c] = Q_min( minout[c], mesh->verts[face->b].point[c] );
maxout[c] = Q_max( maxout[c], mesh->verts[face->b].point[c] );
minout[c] = Q_min( minout[c], mesh->verts[face->c].point[c] );
maxout[c] = Q_max( maxout[c], mesh->verts[face->c].point[c] );
}
}
}
float CWorldRayTrace :: CalculateCostsOfSplit( int plane, int *list, int ntris, vec3_t absmin, vec3_t absmax, float &value, int &nleft, int &nright, int &nboth )
{
float min_coord = 1.0e23;
float max_coord = -1.0e23;
// 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.
for( int t = 0; t < ntris; t++ )
{
tface_t *face = m_TriangleList[list[t]];
// determine max and min coordinate values for later optimization
min_coord = Q_min( min_coord, mesh->verts[face->a].point[plane] );
max_coord = Q_max( max_coord, mesh->verts[face->a].point[plane] );
min_coord = Q_min( min_coord, mesh->verts[face->b].point[plane] );
max_coord = Q_max( max_coord, mesh->verts[face->b].point[plane] );
min_coord = Q_min( min_coord, mesh->verts[face->c].point[plane] );
max_coord = Q_max( max_coord, mesh->verts[face->c].point[plane] );
switch( face->ClassifyAgainstAxisSplit( mesh->verts, plane, value ))
{
case PLANECHECK_NEGATIVE:
face->pcheck0 = PLANECHECK_NEGATIVE;
nleft++;
break;
case PLANECHECK_POSITIVE:
face->pcheck0 = PLANECHECK_POSITIVE;
nright++;
break;
case PLANECHECK_STRADDLING:
face->pcheck0 = PLANECHECK_STRADDLING;
nboth++;
break;
default:
ASSERT( 0 );
break;
}
}
// now, if the split resulted in one half being empty, "grow" the empty half
if( nleft && (nboth == 0) && (nright == 0))
value = max_coord;
if( nright && (nboth == 0) && (nleft == 0))
value = min_coord;
vec3_t leftMins, leftMaxs;
vec3_t rightMins, rightMaxs;
// now, perform surface area/cost check to determine whether this split was worth it
VectorCopy( absmin, leftMins );
VectorCopy( absmax, leftMaxs );
VectorCopy( absmin, rightMins );
VectorCopy( absmax, rightMaxs );
leftMaxs[plane] = value;
rightMins[plane] = value;
float SA_L = BoxSurfaceArea( leftMins, leftMaxs );
float SA_R = BoxSurfaceArea( rightMins, rightMaxs );
float ISA = 1.0 / BoxSurfaceArea( absmin, absmax );
return COST_OF_TRAVERSAL + COST_OF_INTERSECTION * (nboth + (SA_L * ISA * nleft) + (SA_R * ISA * nright));
}
void CWorldRayTrace :: RefineNode( int node_number, int *list, int ntris, vec3_t absmin, vec3_t absmax, int depth )
{
// never split empty lists
if( ntris < 3 )
{
// no point in continuing
m_KDTree[node_number].m_iChild = KDNODE_STATE_LEAF + (m_TriangleIndexList.Count() << 2);
m_KDTree[node_number].SetNumberOfTrianglesInLeafNode( ntris );
for( int t = 0; t < ntris; t++ )
m_TriangleIndexList.AddToTail( list[t] );
return;
}
float best_splitvalue = 0;
float best_cost = 1.0e23;
int best_nleft = 0;
int best_nright = 0;
int best_nboth = 0;
int split_plane = 0;
// don't try all trinagles as split
int tri_skip = 1 + (ntris / 10);
// 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.5f * (absmin[axis] + absmax[axis]);
}
else
{
// else, split at the triangle vertex if possible
tface_t *face = m_TriangleList[list[ts]];
int vert = -1;
switch( tv )
{
case 0: vert = face->a; break;
case 1: vert = face->b; break;
case 2: vert = face->c; break;
}
assert( vert != -1 );
trial_splitvalue = mesh->verts[vert].point[axis];
if(( trial_splitvalue > absmax[axis] ) || ( trial_splitvalue < absmin[axis] ))
continue; // don't try this vertex - not inside
}
float trial_cost = CalculateCostsOfSplit( axis, list, ntris, absmin, absmax,
trial_splitvalue, trial_nleft, trial_nright, trial_nboth );
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++ )
{
tface_t *face = m_TriangleList[list[t]];
face->pcheck1 = face->pcheck0;
}
}
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
m_KDTree[node_number].m_iChild = KDNODE_STATE_LEAF + (m_TriangleIndexList.Count() << 2);
m_KDTree[node_number].SetNumberOfTrianglesInLeafNode( ntris );
for( int t = 0; t < ntris; t++ )
m_TriangleIndexList.AddToTail( list[t] );
}
else
{
// its worth splitting!
// we will achieve the splitting without sorting by using a selection algorithm.
int *new_triangle_list = new int[ntris];
vec3_t leftMins, leftMaxs;
vec3_t rightMins, rightMaxs;
int n_right_output = 0;
int n_left_output = 0;
int n_both_output = 0;
// now, perform surface area/cost check to determine whether this split was worth it
VectorCopy( absmin, leftMins );
VectorCopy( absmax, leftMaxs );
VectorCopy( absmin, rightMins );
VectorCopy( absmax, rightMaxs );
leftMaxs[split_plane] = best_splitvalue;
rightMins[split_plane] = best_splitvalue;
#ifdef _DEBUG
memset( new_triangle_list, 0, sizeof( int ) * ntris );
#endif
for( int t = 0; t < ntris; t++ )
{
tface_t *face = m_TriangleList[list[t]];
switch( face->pcheck1 )
{
case PLANECHECK_NEGATIVE:
new_triangle_list[n_left_output++] = list[t];
break;
case PLANECHECK_POSITIVE:
n_right_output++;
new_triangle_list[ntris - n_right_output] = list[t];
break;
case PLANECHECK_STRADDLING:
new_triangle_list[best_nleft + n_both_output] = list[t];
n_both_output++;
break;
default:
ASSERT( 0 );
break;
}
}
int left_child = m_KDTree.Count();
int right_child = left_child + 1;
m_KDTree[node_number].m_iChild = split_plane + (left_child<<2);
m_KDTree[node_number].m_flSplitValue = best_splitvalue;
KDNode newnode;
m_KDTree.AddToTail( newnode );
m_KDTree.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, leftMaxs, depth + 1 );
RefineNode( right_child, new_triangle_list + best_nleft, best_nright + best_nboth, rightMins, rightMaxs, depth + 1 );
delete[] new_triangle_list;
}
}
void CWorldRayTrace :: BuildTree( tmesh_t *src )
{
KDNode root;
m_KDTree.AddToTail( root );
mesh = src;
int *root_triangle_list = new int[m_TriangleList.Count()];
for( int t = 0; t < m_TriangleList.Count(); t++ )
root_triangle_list[t] = t;
CalculateTriangleListBounds( root_triangle_list, m_TriangleList.Count(), m_AbsMins, m_AbsMaxs );
RefineNode( 0, root_triangle_list, m_TriangleList.Count(), m_AbsMins, m_AbsMaxs, 0 );
delete[] root_triangle_list;
}
void CWorldRayTrace :: TraceRay( const vec3_t start, const vec3_t stop, trace_t *trace )
{
vec3_t direction;
vec_t dist;
VectorSubtract( stop, start, direction );
dist = VectorNormalize( direction );
trace->fraction = dist;
TraceRay( 0.0f, dist, start, direction, trace );
trace->fraction = trace->fraction / dist;
}
_forceinline float Reciprocal( const float a )
{
float ret = 1.0 / a;
// newton iteration is: Y(n+1) = 2*Y(n)-a*Y(n)^2
ret = (( ret + ret ) - ( a * ( ret * ret )));
return ret;
}
bool CWorldRayTrace :: TraceTexture( const tface_t *face, float u, float v, float w )
{
// calculate st from uvw (barycentric) coordinates
float s = w * mesh->verts[face->a].st[0] + u * mesh->verts[face->b].st[0] + v * mesh->verts[face->c].st[0];
float t = w * mesh->verts[face->a].st[1] + u * mesh->verts[face->b].st[1] + v * mesh->verts[face->c].st[1];
// convert ST to real pixels position
int x = fix_coord( s * face->texture->width, face->texture->width - 1 );
int y = fix_coord( t * face->texture->height, face->texture->height - 1 );
// test pixel
if( face->texture->data[(face->texture->width * y) + x] != 255 )
{
// at this point we hit the opaque pixel
return true;
}
return false;
}
void CWorldRayTrace :: TraceRay( vec_t p1f, vec_t p2f, const vec3_t start, const vec3_t direction, trace_t *trace )
{
vec3_t OneOverRayDir;
int c;
VectorCopy( direction, OneOverRayDir );
// add epsilon to avoid division by zero
for( c = 0; c < 3; c++ )
{
if( OneOverRayDir[c] == 0.0f )
OneOverRayDir[c] = FLT_EPSILON;
OneOverRayDir[c] = Reciprocal( OneOverRayDir[c] );
}
// now, clip rays against bounding box
for( c = 0; c < 3; c++ )
{
float isect_min_t = ( m_AbsMins[c] - start[c] ) * OneOverRayDir[c];
float isect_max_t = ( m_AbsMaxs[c] - start[c] ) * OneOverRayDir[c];
p1f = Q_min( p1f, Q_min( isect_min_t, isect_max_t ));
p2f = Q_max( p2f, Q_max( isect_min_t, isect_max_t ));
}
if( p1f > p2f ) return;
// based on ray direction, whether to
// visit left or right node first
int front_idx[3], back_idx[3];
int signbits = SignbitsForPlane( direction );
if( FBitSet( signbits, 1 ))
{
back_idx[0] = 0;
front_idx[0] = 1;
}
else
{
back_idx[0] = 1;
front_idx[0] = 0;
}
if( FBitSet( signbits, 2 ))
{
back_idx[1] = 0;
front_idx[1] = 1;
}
else
{
back_idx[1] = 1;
front_idx[1] = 0;
}
if( FBitSet( signbits, 4 ))
{
back_idx[2] = 0;
front_idx[2] = 1;
}
else
{
back_idx[2] = 1;
front_idx[2] = 0;
}
NodeToVisit NodeQueue[MAX_NODE_STACK_LEN];
NodeToVisit *stack_ptr = &NodeQueue[MAX_NODE_STACK_LEN];
int mailboxids[MAILBOX_HASH_SIZE]; // used to avoid redundant triangle tests
KDNode const *CurNode = &(m_KDTree[0]);
memset( mailboxids, 0xff, sizeof( mailboxids )); // !!speed!! keep around?
while( 1 )
{
// traverse until next leaf
while( CurNode->NodeType() != KDNODE_STATE_LEAF )
{
KDNode const *FrontChild = &(m_KDTree[CurNode->LeftChild()]);
int split_plane_number = CurNode->NodeType();
// dist = (split - org) / dir
float dist_to_sep_plane = ( CurNode->m_flSplitValue - start[split_plane_number] ) * OneOverRayDir[split_plane_number];
bool active = ( p1f <= p2f ); // mask of which rays are active
bool hits_front = active && ( dist_to_sep_plane >= p1f );
// now, decide how to traverse children. can either do front, back, or do front and push back.
if( !hits_front )
{
// missed the front. only traverse back
CurNode = FrontChild + back_idx[split_plane_number];
p1f = Q_max( p1f, dist_to_sep_plane );
}
else
{
bool hits_back = active && ( dist_to_sep_plane <= p2f );
if( !hits_back )
{
// missed the back - only need to traverse front node
CurNode = FrontChild + front_idx[split_plane_number];
p2f = Q_min( p2f, dist_to_sep_plane );
}
else
{
// at least some rays hit both nodes.
// must push far, traverse near
assert( stack_ptr > NodeQueue );
stack_ptr--;
stack_ptr->node = FrontChild + back_idx[split_plane_number];
stack_ptr->p1f = Q_max( p1f, dist_to_sep_plane );
stack_ptr->p2f = p2f;
CurNode = FrontChild + front_idx[split_plane_number];
p2f = Q_min( p2f, dist_to_sep_plane );
}
}
}
// hit a leaf! must do intersection check
int ntris = CurNode->NumberOfTrianglesInLeaf();
if( ntris )
{
int const *tlist = &(m_TriangleIndexList[CurNode->TriangleIndexStart()]);
do
{
int tnum = *(tlist++);
// check mailbox
int mbox_slot = tnum & (MAILBOX_HASH_SIZE - 1);
const tface_t *face = m_TriangleList[tnum];
if( mailboxids[mbox_slot] != tnum )
{
mailboxids[mbox_slot] = tnum;
// compute plane intersection
float DDotN = DotProduct( direction, face->normal );
// mask off zero or near zero (ray parallel to surface)
bool did_hit = (( DDotN > FLT_EPSILON ) || ( DDotN < -FLT_EPSILON ));
if( !FBitSet( mesh->flags, FMESH_VERTEX_LIGHTING|FMESH_MODEL_LIGHTMAPS ))
{
if( !FBitSet( face->texture->flags, STUDIO_NF_TWOSIDE ))
{
// NOTE: probably we should normalize the n but for speed reasons i don't do it
if( DDotN < -FLT_EPSILON ) continue;
}
}
if( !did_hit ) continue; // to prevent division by zero
float numerator = face->NdotP1 - DotProduct( start, face->normal );
float isect_t = numerator / DDotN; // fraction
// now, we have the distance to the plane. lets update our mask
did_hit = did_hit && ( isect_t > m_flBackFrac );
did_hit = did_hit && ( isect_t < trace->fraction );
if( !did_hit ) continue;
// now, check 3 edges
float hitc1 = start[face->pcoord0] + ( isect_t * direction[face->pcoord0] );
float hitc2 = start[face->pcoord1] + ( isect_t * direction[face->pcoord1] );
// do barycentric coordinate check
float B0 = face->edge1[0] * hitc1 + face->edge1[1] * hitc2 + face->edge1[2];
did_hit = did_hit && ( B0 >= 0.0f );
float B1 = face->edge2[0] * hitc1 + face->edge2[1] * hitc2 + face->edge2[2];
did_hit = did_hit && ( B1 >= 0.0f );
float B2 = B0 + B1;
did_hit = did_hit && ( B2 <= 1.0f );
if( !did_hit ) continue;
// if the triangle is transparent
if( face->texture->data )
{
// 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
if( !TraceTexture( face, 1.0 - B2, B0, B1 ))
continue; // passed through alpha-pixel
}
if( trace->fraction > isect_t )
{
trace->contents = face->contents;
trace->fraction = isect_t;
}
}
} while( --ntris );
// now, check if all rays have terminated
if( p2f > trace->fraction ) return;
}
if( stack_ptr == &NodeQueue[MAX_NODE_STACK_LEN] )
return;
// pop stack!
CurNode = stack_ptr->node;
p1f = stack_ptr->p1f;
p2f = stack_ptr->p2f;
stack_ptr++;
}
}
#endif