src/hd/sm_stree.c

/* Copyright (c) 1998 Silicon Graphics, Inc. */

#ifndef lint
static char SCCSid[] = "$SunId$ SGI";
#endif

/*
 * sm_stree.c
 *  An stree (spherical quadtree) is defined by an octahedron in 
 *  canonical form,and a world center point. Each face of the 
 *  octahedron is adaptively subdivided as a planar triangular quadtree.
 *  World space geometry is projected onto the quadtree faces from the
 *  sphere center.
 */
#include "standard.h"
#include "sm_list.h"
#include "sm_flag.h"
#include "sm_geom.h"
#include "sm_qtree.h"
#include "sm_stree.h"

#ifdef TEST_DRIVER
extern FVECT Pick_point[500],Pick_v0[500],Pick_v1[500],Pick_v2[500];
extern int Pick_cnt;
#endif
/* octahedron coordinates */
FVECT stDefault_base[6] = {  {1.,0.,0.},{0.,1.,0.}, {0.,0.,1.},
                            {-1.,0.,0.},{0.,-1.,0.},{0.,0.,-1.}};
/* octahedron triangle vertices */
int stBase_verts[8][3] = { {0,1,2},{3,1,2},{0,4,2},{3,4,2},
                           {0,1,5},{3,1,5},{0,4,5},{3,4,5}};
/* octahedron triangle nbrs ; nbr i is the face opposite vertex i*/
int stBase_nbrs[8][3] =  { {1,2,4},{0,3,5},{3,0,6},{2,1,7},
                           {5,6,0},{4,7,1},{7,4,2},{6,5,3}};
int stRoot_indices[8][3] = {{1,1,1},{-1,1,1},{1,-1,1},{-1,-1,1},
                            {1,1,-1},{-1,1,-1},{1,-1,-1},{-1,-1,-1}};
/*
 +z   y                -z   y
      |                     |
 1    |   0             5   |   4
______|______ x      _______|______ x
 3    |   2             7   |   6   
      |                     |

Nbrs
  +z   y                -z   y
     /0|1\                 /1|0\
 5  /  |  \ 4             /  |  \
   /(1)|(0)\           1 /(5)|(4)\ 0
  /    |    \           /    |    \
 /2   1|0   2\         /2   0|1   2\
/------|------\x      /------|------\x
\0    1|2    0/       \0    2|2    1/ 
 \     |     /         \     |     /
 7\ (3)|(2) / 6       3 \ (7)|(6) / 2
   \   |   /             \   |   /
    \ 2|1 /               \ 1|0 /
*/


stInit(st)
STREE *st;
{
  int i,j;

    ST_TOP_QT(st) = qtAlloc();
    ST_BOTTOM_QT(st) = qtAlloc();
    /* Clear the children */

   QT_CLEAR_CHILDREN(ST_TOP_QT(st));
   QT_CLEAR_CHILDREN(ST_BOTTOM_QT(st));
}

/* Frees the children of the 2 quadtrees rooted at st,
   Does not free root nodes: just clears
*/
stClear(st)
   STREE *st;
{
    qtDone();
    stInit(st);
}

/* Allocates a stree structure  and creates octahedron base */
STREE
*stAlloc(st)
STREE *st;
{
  int i,m;
  FVECT v0,v1,v2;
  FVECT n;
  
  if(!st)
    if(!(st = (STREE *)malloc(sizeof(STREE))))
       error(SYSTEM,"stAlloc(): Unable to allocate memory\n");

  /* Allocate the top and bottom quadtree root nodes */
  stInit(st);
  
  
  /* will go ********************************************/
  /* Set the octahedron base */
  ST_SET_BASE(st,stDefault_base);

  /* Calculate octahedron face and edge normals */
  for(i=0; i < ST_NUM_ROOT_NODES; i++)
  {
      VCOPY(v0,ST_NTH_V(st,i,0));
      VCOPY(v1,ST_NTH_V(st,i,1));
      VCOPY(v2,ST_NTH_V(st,i,2));
      tri_plane_equation(v0,v1,v2, &ST_NTH_PLANE(st,i),FALSE);
      m = max_index(FP_N(ST_NTH_PLANE(st,i)),NULL);
      FP_X(ST_NTH_PLANE(st,i)) = (m+1)%3; 
      FP_Y(ST_NTH_PLANE(st,i)) = (m+2)%3;
      FP_Z(ST_NTH_PLANE(st,i)) = m;
      VCROSS(ST_EDGE_NORM(st,i,0),v0,v1);
      VCROSS(ST_EDGE_NORM(st,i,1),v1,v2);
      VCROSS(ST_EDGE_NORM(st,i,2),v2,v0);
  }

  /*****************************************************************/
  return(st);
}

#define BARY_INT(v,b)  if((v)>2.0) (b) = MAXBCOORD;else \
  if((v)<-2.0) (b)=-MAXBCOORD;else (b)=(BCOORD)((v)*MAXBCOORD2);

vert_to_qt_frame(root,v,b)
int root;
FVECT v;
BCOORD b[3];
{
  int i;
  double scale;
  double d0,d1,d2;

  if(STR_NTH_INDEX(root,0)==-1)
    d0 = -v[0];
  else
    d0 = v[0];
  if(STR_NTH_INDEX(root,1)==-1)
    d1 = -v[1];
  else
    d1 = v[1];
  if(STR_NTH_INDEX(root,2)==-1)
    d2 = -v[2];
  else
    d2 = v[2];

  /* Plane is now x+y+z = 1 - intersection of pt ray is qtv/den */
  scale = 1.0/ (d0 + d1 + d2);
  d0 *= scale;
  d1 *= scale;
  d2 *= scale;

  BARY_INT(d0,b[0])
  BARY_INT(d1,b[1])
  BARY_INT(d2,b[2])
}


ray_to_qt_frame(root,v,dir,b,db)
int root;
FVECT v,dir;
BCOORD b[3],db[3];
{
  int i;
  double scale;
  double d0,d1,d2;
  double dir0,dir1,dir2;

  if(STR_NTH_INDEX(root,0)==-1)
  {
    d0 = -v[0];
    dir0 = -dir[0];
  }
  else
  {
    d0 = v[0];
    dir0 = dir[0];
  }
  if(STR_NTH_INDEX(root,1)==-1)
  {
    d1 = -v[1];
    dir1 = -dir[1];
  }
  else
  {
    d1 = v[1];
    dir1 = dir[1];
  }
  if(STR_NTH_INDEX(root,2)==-1)
  {
    d2 = -v[2];
    dir2 = -dir[2];
  }
  else
  {
    d2 = v[2];
    dir2 = dir[2];
  }
  /* Plane is now x+y+z = 1 - intersection of pt ray is qtv/den */
  scale = 1.0/ (d0 + d1 + d2);
  d0 *= scale;
  d1 *= scale;
  d2 *= scale;

  /* Calculate intersection point of orig+dir: This calculation is done
     after the origin is projected into the plane in order to constrain
     the projection( i.e. the size of the projection of the unit direction
     vector translated to the origin depends on how close
     the origin is to the view center
     */
  /* Must divide by at least root2 to insure that projection will fit 
     int [-2,2] bounds: assumed length is 1: therefore greatest projection
     from endpoint of triangle is at 45 degrees or projected length of root2
  */
  dir0 = d0 + dir0*0.5;
  dir1 = d1 + dir1*0.5;
  dir2 = d2 + dir2*0.5;

  scale = 1.0/ (dir0 + dir1 + dir2);
  dir0 *= scale;
  dir1 *= scale;
  dir2 *= scale;

  BARY_INT(d0,b[0])
  BARY_INT(d1,b[1])
  BARY_INT(d2,b[2])
  BARY_INT(dir0,db[0])
  BARY_INT(dir1,db[1])
  BARY_INT(dir2,db[2])

  db[0]  -= b[0];
  db[1]  -= b[1];
  db[2]  -= b[2];
}

qt_frame_to_vert(root,b,v)
int root;
BCOORD b[3];
FVECT v;
{
  int i;
  double d0,d1,d2;

  d0 = b[0]/(double)MAXBCOORD2;
  d1 = b[1]/(double)MAXBCOORD2;
  d2 = b[2]/(double)MAXBCOORD2;
  
  if(STR_NTH_INDEX(root,0)==-1)
    v[0] = -d0;
  else
    v[0] = d0;
  if(STR_NTH_INDEX(root,1)==-1)
    v[1] = -d1;
  else
    v[1] = d1;
  if(STR_NTH_INDEX(root,2)==-1)
    v[2] = -d2;
  else
    v[2] = d2;
}


/* Return quadtree leaf node containing point 'p'*/
QUADTREE
stPoint_locate(st,p)
    STREE *st;
    FVECT p;
{
    QUADTREE qt;
    BCOORD bcoordi[3];
    int i;

    /* Find root quadtree that contains p */
    i = stLocate_root(p);
    qt = ST_ROOT_QT(st,i);
    
     /* Will return lowest level triangle containing point: It the
       point is on an edge or vertex: will return first associated
       triangle encountered in the child traversal- the others can
       be derived using triangle adjacency information
    */
    if(QT_IS_TREE(qt))
    {  
      vert_to_qt_frame(i,p,bcoordi);
      i = bary_child(bcoordi);
      return(qtLocate(QT_NTH_CHILD(qt,i),bcoordi));
    }
    else
      return(qt);
}

static unsigned int nbr_b[8][3] ={{2,4,16},{1,8,32},{8,1,64},{4,2,128},
                           {32,64,1},{16,128,2},{128,16,4},{64,32,8}};
unsigned int
stTri_cells(st,v)
     STREE *st;
     FVECT v[3];
{
  unsigned int cells,cross;
  unsigned int vcell[3];
  double t0,t1;
  int i,inext;

  /* First find base cells that tri vertices are in (0-7)*/
  vcell[0] = stLocate_root(v[0]);
  vcell[1] = stLocate_root(v[1]);
  vcell[2] = stLocate_root(v[2]);

  /* If all in same cell- return that bit only */
  if(vcell[0] == vcell[1] && vcell[1] == vcell[2])
    return( 1 << vcell[0]);

  cells = 0;
  for(i=0;i<3; i++)
  {
    if(i==2)
      inext = 0;
    else
      inext = i+1;
    /* Mark cell containing initial vertex */
    cells |= 1 << vcell[i];

    /* Take the exclusive or: will have bits set where edge crosses axis=0*/
    cross = vcell[i] ^ vcell[inext];
    /* If crosses 2 planes: then have 2 options for edge crossing-pick closest
     otherwise just hits two*/
    /* Neighbors are zyx */
    switch(cross){
    case 3: /* crosses x=0 and y=0 */
      t0 = -v[i][0]/(v[inext][0]-v[i][0]);
      t1 = -v[i][1]/(v[inext][1]-v[i][1]);
      if(t0==t1)
        break;
      else if(t0 < t1)
        cells |= nbr_b[vcell[i]][0];
          else
            cells |= nbr_b[vcell[i]][1];
      break;
    case 5: /* crosses x=0 and z=0 */
      t0 = -v[i][0]/(v[inext][0]-v[i][0]);
      t1 = -v[i][2]/(v[inext][2]-v[i][2]);
      if(t0==t1)
        break;
      else if(t0 < t1)
        cells |= nbr_b[vcell[i]][0];
          else
            cells |=nbr_b[vcell[i]][2];

      break;
    case 6:/* crosses  z=0 and y=0 */
      t0 = -v[i][2]/(v[inext][2]-v[i][2]);
      t1 = -v[i][1]/(v[inext][1]-v[i][1]);
      if(t0==t1)
        break;
      else if(t0 < t1)
      {
        cells |= nbr_b[vcell[i]][2];
      }
      else
      {
        cells |=nbr_b[vcell[i]][1];
      }
      break;
    case 7:
      error(CONSISTENCY," Insert:Edge shouldnt be able to span 3 cells");
      break;
    }
  }
  return(cells);
}


stRoot_trace_ray(qt,root,orig,dir,nextptr,func,f)
   QUADTREE qt;
   int root;
   FVECT orig,dir;
   int *nextptr;
   FUNC func;
   int *f;
{
  double br[3];
  BCOORD bi[3],dbi[3];
  
  /* Project the origin onto the root node plane */
  /* Find the intersection point of the origin */
  ray_to_qt_frame(root,orig,dir,bi,dbi);

  /* trace the ray starting with this node */
  qtTrace_ray(qt,bi,dbi[0],dbi[1],dbi[2],nextptr,0,0,func,f);
  if(!QT_FLAG_IS_DONE(*f))
    qt_frame_to_vert(root,bi,orig);

}

/* Trace ray 'orig-dir' through stree and apply 'func(qtptr,f,argptr)' at each
   node that it intersects
*/
int
stTrace_ray(st,orig,dir,func)
   STREE *st;
   FVECT orig,dir;
   FUNC func;
{
    int next,last,i,f=0;
    QUADTREE qt;
    FVECT o,n,v;
    double pd,t,d;

    VCOPY(o,orig);
#ifdef TEST_DRIVER
       Pick_cnt=0;
#endif;
    /* Find the qt node that o falls in */
    i = stLocate_root(o);
    qt = ST_ROOT_QT(st,i);
    
    stRoot_trace_ray(qt,i,o,dir,&next,func,&f);

    if(QT_FLAG_IS_DONE(f))
      return(TRUE);
    /*    
    d = DOT(orig,dir)/sqrt(DOT(orig,orig));
    VSUM(v,orig,dir,-d);
    */
    /* Crossed over to next cell: id = nbr */
    while(1) 
      {
        /* test if ray crosses plane between this quadtree triangle and
           its neighbor- if it does then find intersection point with 
           ray and plane- this is the new origin
           */
        if(next == INVALID)
          return(FALSE);
        /*
        if(DOT(o,v) < 0.0)
          return(FALSE);
          */
        i = stBase_nbrs[i][next];
        qt = ST_ROOT_QT(st,i);
        stRoot_trace_ray(qt,i,o,dir,&next,func,&f);
        if(QT_FLAG_IS_DONE(f))
          return(TRUE);
      }
}


stVisit_poly(st,verts,l,root,func)
STREE *st;
FVECT *verts;
LIST *l;
unsigned int root;
FUNC func;
{
  int id0,id1,id2;
  FVECT tri[3];

  id0 = pop_list(&l);
  id1 = pop_list(&l);
  while(l)
  {
    id2 = pop_list(&l);
    VCOPY(tri[0],verts[id0]);
    VCOPY(tri[1],verts[id1]);
    VCOPY(tri[2],verts[id2]);
    stRoot_visit_tri(st,root,tri,func);
    id1 = id2;
  }
}

stVisit_clip(st,i,verts,vcnt,l,cell,func)
     STREE *st;
     int i;
     FVECT *verts;
     int *vcnt;
     LIST *l;
     unsigned int cell;
     FUNC func;
{

  LIST *labove,*lbelow,*endb,*enda;
  int last = -1;
  int id,last_id;
  int first,first_id;
  unsigned int cellb;

  labove = lbelow = NULL;
  enda = endb = NULL;
  while(l)
  {
    id = pop_list(&l);
    if(ZERO(verts[id][i]))
    {
      if(last==-1)
      {/* add below and above */
        first = 2;
        first_id= id;
      }
      lbelow=add_data(lbelow,id,&endb);
      labove=add_data(labove,id,&enda);
      last_id = id;
      last = 2;
      continue;
    }
    if(verts[id][i] < 0)
    {
      if(last != 1)
      {
        lbelow=add_data(lbelow,id,&endb);
        if(last==-1)
        {
          first = 0;
          first_id = id;
        }
        last_id = id;
        last = 0;
        continue;
      }
      /* intersect_edges */
      intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
      /*newpoint goes to above and below*/
      lbelow=add_data(lbelow,*vcnt,&endb);
      lbelow=add_data(lbelow,id,&endb);
      labove=add_data(labove,*vcnt,&enda);
      last = 0;
      last_id = id;
      (*vcnt)++;
    }
    else
    {
      if(last != 0)
      {
        labove=add_data(labove,id,&enda);
        if(last==-1)
        {
          first = 1;
          first_id = id;
        }
        last_id = id;
        last = 1;
        continue;
      }
      /* intersect_edges */
      /*newpoint goes to above and below*/
      intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
      lbelow=add_data(lbelow,*vcnt,&endb);
      labove=add_data(labove,*vcnt,&enda);
      labove=add_data(labove,id,&enda);
      last_id = id;
      (*vcnt)++;
      last = 1;
    }
  }
  if(first != 2 && first != last)
  {
    intersect_edge_coord_plane(verts[id],verts[first_id],i,verts[*vcnt]);
    /*newpoint goes to above and below*/
    lbelow=add_data(lbelow,*vcnt,&endb);
    labove=add_data(labove,*vcnt,&enda);
    (*vcnt)++;

  }
  if(i==2)
  {
    if(lbelow)
    {
      if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
      {
        cellb = cell | (1 << i);
        stVisit_poly(st,verts,lbelow,cellb,func);
      }
      else
        free_list(lbelow);
    }
    if(labove)
     {
      if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
        stVisit_poly(st,verts,labove,cell,func);
      else
        free_list(labove);
     }
  }
  else
  {
    if(lbelow)
    {
      if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
        {
          cellb = cell | (1 << i);
          stVisit_clip(st,i+1,verts,vcnt,lbelow,cellb,func);
        }
      else
        free_list(lbelow);
    }
    if(labove)
     {
       if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
         stVisit_clip(st,i+1,verts,vcnt,labove,cell,func);
       else
         free_list(labove);
     }
  }

}

stVisit(st,tri,func)
   STREE *st;
   FVECT tri[3];
   FUNC func;
{
    int r0,r1,r2;
    LIST *l;

    r0 = stLocate_root(tri[0]);
    r1 = stLocate_root(tri[1]);
    r2 = stLocate_root(tri[2]);
    if(r0 == r1 && r1==r2)
      stRoot_visit_tri(st,r0,tri,func);
    else
      {
        FVECT verts[ST_CLIP_VERTS];
        int cnt;

        VCOPY(verts[0],tri[0]);
        VCOPY(verts[1],tri[1]);
        VCOPY(verts[2],tri[2]); 
        
        l = add_data(NULL,0,NULL);
        l = add_data(l,1,NULL);
        l = add_data(l,2,NULL);
        cnt = 3;
        stVisit_clip(st,0,verts,&cnt,l,0,func);
      }
}


/* New Insertion code!!! */


BCOORD qtRoot[3][3] = { {MAXBCOORD2,0,0},{0,MAXBCOORD2,0},{0,0,MAXBCOORD2}};


convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02)
int root;
FVECT tri[3];
BCOORD b0[3],b1[3],b2[3];
BCOORD db10[3],db21[3],db02[3];
{
  /* Project the vertex into the qtree plane */
  vert_to_qt_frame(root,tri[0],b0);
  vert_to_qt_frame(root,tri[1],b1);
  vert_to_qt_frame(root,tri[2],b2);

  /* calculate triangle edge differences in new frame */
  db10[0] = b1[0] - b0[0]; db10[1] = b1[1] - b0[1]; db10[2] = b1[2] - b0[2];
  db21[0] = b2[0] - b1[0]; db21[1] = b2[1] - b1[1]; db21[2] = b2[2] - b1[2];
  db02[0] = b0[0] - b2[0]; db02[1] = b0[1] - b2[1]; db02[2] = b0[2] - b2[2];
}


QUADTREE
stRoot_insert_tri(st,root,tri,f)
   STREE *st;
   int root;
   FVECT tri[3];
   FUNC f;
{
  BCOORD b0[3],b1[3],b2[3];
  BCOORD db10[3],db21[3],db02[3];
  unsigned int s0,s1,s2,sq0,sq1,sq2;
  QUADTREE qt;

  /* Map the triangle vertices into the canonical barycentric frame */
  convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);

  /* Calculate initial sidedness info */
  SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
  SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);

  qt = ST_ROOT_QT(st,root);
  /* Visit cells that triangle intersects */
  qt = qtInsert_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
       b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f,0);

  return(qt);
}

stRoot_visit_tri(st,root,tri,f)
   STREE *st;
   int root;
   FVECT tri[3];
   FUNC f;
{
  BCOORD b0[3],b1[3],b2[3];
  BCOORD db10[3],db21[3],db02[3];
  unsigned int s0,s1,s2,sq0,sq1,sq2;
  QUADTREE qt;

  /* Map the triangle vertices into the canonical barycentric frame */
  convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);

  /* Calculate initial sidedness info */
  SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
  SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);

  qt = ST_ROOT_QT(st,root);
  QT_SET_FLAG(ST_QT(st,root));
  /* Visit cells that triangle intersects */
  qtVisit_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
       b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f);

}

stInsert_tri(st,tri,f)
   STREE *st;
   FVECT tri[3];
   FUNC f;
{
  unsigned int cells,which;
  int root;
  

  /* calculate entry/exit points of edges through the cells */
  cells = stTri_cells(st,tri);

  /* For each cell that quadtree intersects: Map the triangle vertices into
     the canonical barycentric frame of (1,0,0), (0,1,0),(0,0,1). Insert 
     by first doing a trivial reject on the interior nodes, and then a 
     tri/tri intersection at the leaf nodes.
  */
  for(root=0,which=1; root < ST_NUM_ROOT_NODES; root++,which <<= 1)
  {
    /* For each of the quadtree roots: check if marked as intersecting tri*/
    if(cells & which)
      /* Visit tri cells */
      ST_ROOT_QT(st,root) = stRoot_insert_tri(st,root,tri,f);
  }
}


Revision:	3.9
Committed:	Tue Jan 5 16:52:39 1999 UTC (25 years, 3 months ago) by gwlarson
Content type:	text/plain
Branch:	MAIN
Changes since 3.8:	+0 -0 lines
Log Message:	fixed logic in smTest to handle nearby and coincident points, added base tris to rendering, made list of new triangles to speed up rendering
#	User	Rev	Content
1	gwlarson	3.1	/* Copyright (c) 1998 Silicon Graphics, Inc. */
2
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ SGI";
5			#endif
6
7			/*
8			* sm_stree.c
9	gwlarson	3.6	* An stree (spherical quadtree) is defined by an octahedron in
10			* canonical form,and a world center point. Each face of the
11			* octahedron is adaptively subdivided as a planar triangular quadtree.
12			* World space geometry is projected onto the quadtree faces from the
13			* sphere center.
14	gwlarson	3.1	*/
15			#include "standard.h"
16	gwlarson	3.8	#include "sm_list.h"
17	gwlarson	3.6	#include "sm_flag.h"
18	gwlarson	3.1	#include "sm_geom.h"
19	gwlarson	3.6	#include "sm_qtree.h"
20	gwlarson	3.1	#include "sm_stree.h"
21
22	gwlarson	3.4	#ifdef TEST_DRIVER
23			extern FVECT Pick_point[500],Pick_v0[500],Pick_v1[500],Pick_v2[500];
24			extern int Pick_cnt;
25			#endif
26	gwlarson	3.6	/* octahedron coordinates */
27			FVECT stDefault_base[6] = { {1.,0.,0.},{0.,1.,0.}, {0.,0.,1.},
28			{-1.,0.,0.},{0.,-1.,0.},{0.,0.,-1.}};
29			/* octahedron triangle vertices */
30	gwlarson	3.8	int stBase_verts[8][3] = { {0,1,2},{3,1,2},{0,4,2},{3,4,2},
31			{0,1,5},{3,1,5},{0,4,5},{3,4,5}};
32	gwlarson	3.6	/* octahedron triangle nbrs ; nbr i is the face opposite vertex i*/
33	gwlarson	3.8	int stBase_nbrs[8][3] = { {1,2,4},{0,3,5},{3,0,6},{2,1,7},
34			{5,6,0},{4,7,1},{7,4,2},{6,5,3}};
35			int stRoot_indices[8][3] = {{1,1,1},{-1,1,1},{1,-1,1},{-1,-1,1},
36			{1,1,-1},{-1,1,-1},{1,-1,-1},{-1,-1,-1}};
37			/*
38			+z y -z y
39			\| \|
40			1 \| 0 5 \| 4
41			______\|______ x _______\|______ x
42			3 \| 2 7 \| 6
43			\| \|
44	gwlarson	3.1
45	gwlarson	3.8	Nbrs
46			+z y -z y
47			/0\|1\ /1\|0\
48			5 / \| \ 4 / \| \
49			/(1)\|(0)\ 1 /(5)\|(4)\ 0
50			/ \| \ / \| \
51			/2 1\|0 2\ /2 0\|1 2\
52			/------\|------\x /------\|------\x
53			\0 1\|2 0/ \0 2\|2 1/
54			\ \| / \ \| /
55			7\ (3)\|(2) / 6 3 \ (7)\|(6) / 2
56			\ \| / \ \| /
57			\ 2\|1 / \ 1\|0 /
58			*/
59	gwlarson	3.6
60	gwlarson	3.8
61	gwlarson	3.6	stInit(st)
62	gwlarson	3.1	STREE *st;
63			{
64	gwlarson	3.8	int i,j;
65
66			ST_TOP_QT(st) = qtAlloc();
67			ST_BOTTOM_QT(st) = qtAlloc();
68			/* Clear the children */
69
70			QT_CLEAR_CHILDREN(ST_TOP_QT(st));
71			QT_CLEAR_CHILDREN(ST_BOTTOM_QT(st));
72	gwlarson	3.1	}
73
74	gwlarson	3.6	/* Frees the children of the 2 quadtrees rooted at st,
75			Does not free root nodes: just clears
76			*/
77	gwlarson	3.1	stClear(st)
78	gwlarson	3.6	STREE *st;
79	gwlarson	3.1	{
80	gwlarson	3.6	qtDone();
81			stInit(st);
82	gwlarson	3.1	}
83
84	gwlarson	3.6	/* Allocates a stree structure and creates octahedron base */
85	gwlarson	3.1	STREE
86			*stAlloc(st)
87			STREE *st;
88			{
89	gwlarson	3.6	int i,m;
90			FVECT v0,v1,v2;
91			FVECT n;
92
93	gwlarson	3.1	if(!st)
94	gwlarson	3.6	if(!(st = (STREE *)malloc(sizeof(STREE))))
95			error(SYSTEM,"stAlloc(): Unable to allocate memory\n");
96	gwlarson	3.1
97	gwlarson	3.6	/* Allocate the top and bottom quadtree root nodes */
98			stInit(st);
99
100	gwlarson	3.8
101			/* will go ********************************************/
102	gwlarson	3.6	/* Set the octahedron base */
103			ST_SET_BASE(st,stDefault_base);
104	gwlarson	3.1
105	gwlarson	3.6	/* Calculate octahedron face and edge normals */
106			for(i=0; i < ST_NUM_ROOT_NODES; i++)
107			{
108			VCOPY(v0,ST_NTH_V(st,i,0));
109			VCOPY(v1,ST_NTH_V(st,i,1));
110			VCOPY(v2,ST_NTH_V(st,i,2));
111			tri_plane_equation(v0,v1,v2, &ST_NTH_PLANE(st,i),FALSE);
112			m = max_index(FP_N(ST_NTH_PLANE(st,i)),NULL);
113			FP_X(ST_NTH_PLANE(st,i)) = (m+1)%3;
114			FP_Y(ST_NTH_PLANE(st,i)) = (m+2)%3;
115			FP_Z(ST_NTH_PLANE(st,i)) = m;
116	gwlarson	3.7	VCROSS(ST_EDGE_NORM(st,i,0),v0,v1);
117			VCROSS(ST_EDGE_NORM(st,i,1),v1,v2);
118			VCROSS(ST_EDGE_NORM(st,i,2),v2,v0);
119	gwlarson	3.6	}
120	gwlarson	3.8
121			/*****************************************************************/
122	gwlarson	3.1	return(st);
123			}
124
125	gwlarson	3.8	#define BARY_INT(v,b) if((v)>2.0) (b) = MAXBCOORD;else \
126			if((v)<-2.0) (b)=-MAXBCOORD;else (b)=(BCOORD)((v)*MAXBCOORD2);
127	gwlarson	3.1
128	gwlarson	3.8	vert_to_qt_frame(root,v,b)
129			int root;
130			FVECT v;
131			BCOORD b[3];
132			{
133			int i;
134			double scale;
135			double d0,d1,d2;
136
137			if(STR_NTH_INDEX(root,0)==-1)
138			d0 = -v[0];
139			else
140			d0 = v[0];
141			if(STR_NTH_INDEX(root,1)==-1)
142			d1 = -v[1];
143			else
144			d1 = v[1];
145			if(STR_NTH_INDEX(root,2)==-1)
146			d2 = -v[2];
147			else
148			d2 = v[2];
149
150			/* Plane is now x+y+z = 1 - intersection of pt ray is qtv/den */
151			scale = 1.0/ (d0 + d1 + d2);
152			d0 *= scale;
153			d1 *= scale;
154			d2 *= scale;
155
156			BARY_INT(d0,b[0])
157			BARY_INT(d1,b[1])
158			BARY_INT(d2,b[2])
159			}
160
161
162
163
164			ray_to_qt_frame(root,v,dir,b,db)
165			int root;
166			FVECT v,dir;
167			BCOORD b[3],db[3];
168			{
169			int i;
170			double scale;
171			double d0,d1,d2;
172			double dir0,dir1,dir2;
173
174			if(STR_NTH_INDEX(root,0)==-1)
175			{
176			d0 = -v[0];
177			dir0 = -dir[0];
178			}
179			else
180			{
181			d0 = v[0];
182			dir0 = dir[0];
183			}
184			if(STR_NTH_INDEX(root,1)==-1)
185			{
186			d1 = -v[1];
187			dir1 = -dir[1];
188			}
189			else
190			{
191			d1 = v[1];
192			dir1 = dir[1];
193			}
194			if(STR_NTH_INDEX(root,2)==-1)
195			{
196			d2 = -v[2];
197			dir2 = -dir[2];
198			}
199			else
200			{
201			d2 = v[2];
202			dir2 = dir[2];
203			}
204			/* Plane is now x+y+z = 1 - intersection of pt ray is qtv/den */
205			scale = 1.0/ (d0 + d1 + d2);
206			d0 *= scale;
207			d1 *= scale;
208			d2 *= scale;
209
210			/* Calculate intersection point of orig+dir: This calculation is done
211			after the origin is projected into the plane in order to constrain
212			the projection( i.e. the size of the projection of the unit direction
213			vector translated to the origin depends on how close
214			the origin is to the view center
215			*/
216			/* Must divide by at least root2 to insure that projection will fit
217			int [-2,2] bounds: assumed length is 1: therefore greatest projection
218			from endpoint of triangle is at 45 degrees or projected length of root2
219			*/
220			dir0 = d0 + dir0*0.5;
221			dir1 = d1 + dir1*0.5;
222			dir2 = d2 + dir2*0.5;
223
224			scale = 1.0/ (dir0 + dir1 + dir2);
225			dir0 *= scale;
226			dir1 *= scale;
227			dir2 *= scale;
228
229			BARY_INT(d0,b[0])
230			BARY_INT(d1,b[1])
231			BARY_INT(d2,b[2])
232			BARY_INT(dir0,db[0])
233			BARY_INT(dir1,db[1])
234			BARY_INT(dir2,db[2])
235
236			db[0] -= b[0];
237			db[1] -= b[1];
238			db[2] -= b[2];
239			}
240
241			qt_frame_to_vert(root,b,v)
242			int root;
243			BCOORD b[3];
244			FVECT v;
245			{
246			int i;
247			double d0,d1,d2;
248
249			d0 = b[0]/(double)MAXBCOORD2;
250			d1 = b[1]/(double)MAXBCOORD2;
251			d2 = b[2]/(double)MAXBCOORD2;
252
253			if(STR_NTH_INDEX(root,0)==-1)
254			v[0] = -d0;
255			else
256			v[0] = d0;
257			if(STR_NTH_INDEX(root,1)==-1)
258			v[1] = -d1;
259			else
260			v[1] = d1;
261			if(STR_NTH_INDEX(root,2)==-1)
262			v[2] = -d2;
263			else
264			v[2] = d2;
265			}
266
267
268	gwlarson	3.6	/* Return quadtree leaf node containing point 'p'*/
269	gwlarson	3.3	QUADTREE
270	gwlarson	3.6	stPoint_locate(st,p)
271	gwlarson	3.1	STREE *st;
272	gwlarson	3.2	FVECT p;
273	gwlarson	3.1	{
274	gwlarson	3.8	QUADTREE qt;
275			BCOORD bcoordi[3];
276	gwlarson	3.6	int i;
277	gwlarson	3.1
278	gwlarson	3.6	/* Find root quadtree that contains p */
279	gwlarson	3.8	i = stLocate_root(p);
280			qt = ST_ROOT_QT(st,i);
281	gwlarson	3.1
282	gwlarson	3.8	/* Will return lowest level triangle containing point: It the
283			point is on an edge or vertex: will return first associated
284			triangle encountered in the child traversal- the others can
285			be derived using triangle adjacency information
286			*/
287			if(QT_IS_TREE(qt))
288			{
289			vert_to_qt_frame(i,p,bcoordi);
290			i = bary_child(bcoordi);
291			return(qtLocate(QT_NTH_CHILD(qt,i),bcoordi));
292			}
293			else
294			return(qt);
295	gwlarson	3.1	}
296
297	gwlarson	3.8	static unsigned int nbr_b[8][3] ={{2,4,16},{1,8,32},{8,1,64},{4,2,128},
298			{32,64,1},{16,128,2},{128,16,4},{64,32,8}};
299			unsigned int
300			stTri_cells(st,v)
301			STREE *st;
302			FVECT v[3];
303	gwlarson	3.1	{
304	gwlarson	3.8	unsigned int cells,cross;
305			unsigned int vcell[3];
306			double t0,t1;
307			int i,inext;
308	gwlarson	3.2
309	gwlarson	3.8	/* First find base cells that tri vertices are in (0-7)*/
310			vcell[0] = stLocate_root(v[0]);
311			vcell[1] = stLocate_root(v[1]);
312			vcell[2] = stLocate_root(v[2]);
313
314			/* If all in same cell- return that bit only */
315			if(vcell[0] == vcell[1] && vcell[1] == vcell[2])
316			return( 1 << vcell[0]);
317
318			cells = 0;
319			for(i=0;i<3; i++)
320	gwlarson	3.1	{
321	gwlarson	3.8	if(i==2)
322			inext = 0;
323			else
324			inext = i+1;
325			/* Mark cell containing initial vertex */
326			cells \|= 1 << vcell[i];
327	gwlarson	3.1
328	gwlarson	3.8	/* Take the exclusive or: will have bits set where edge crosses axis=0*/
329			cross = vcell[i] ^ vcell[inext];
330			/* If crosses 2 planes: then have 2 options for edge crossing-pick closest
331			otherwise just hits two*/
332			/* Neighbors are zyx */
333			switch(cross){
334			case 3: /* crosses x=0 and y=0 */
335			t0 = -v[i][0]/(v[inext][0]-v[i][0]);
336			t1 = -v[i][1]/(v[inext][1]-v[i][1]);
337			if(t0==t1)
338			break;
339			else if(t0 < t1)
340			cells \|= nbr_b[vcell[i]][0];
341			else
342			cells \|= nbr_b[vcell[i]][1];
343			break;
344			case 5: /* crosses x=0 and z=0 */
345			t0 = -v[i][0]/(v[inext][0]-v[i][0]);
346			t1 = -v[i][2]/(v[inext][2]-v[i][2]);
347			if(t0==t1)
348			break;
349			else if(t0 < t1)
350			cells \|= nbr_b[vcell[i]][0];
351			else
352			cells \|=nbr_b[vcell[i]][2];
353	gwlarson	3.1
354	gwlarson	3.8	break;
355			case 6:/* crosses z=0 and y=0 */
356			t0 = -v[i][2]/(v[inext][2]-v[i][2]);
357			t1 = -v[i][1]/(v[inext][1]-v[i][1]);
358			if(t0==t1)
359			break;
360			else if(t0 < t1)
361			{
362			cells \|= nbr_b[vcell[i]][2];
363			}
364			else
365			{
366			cells \|=nbr_b[vcell[i]][1];
367			}
368			break;
369			case 7:
370			error(CONSISTENCY," Insert:Edge shouldnt be able to span 3 cells");
371			break;
372			}
373	gwlarson	3.1	}
374	gwlarson	3.8	return(cells);
375	gwlarson	3.1	}
376
377	gwlarson	3.6
378	gwlarson	3.8	stRoot_trace_ray(qt,root,orig,dir,nextptr,func,f)
379			QUADTREE qt;
380			int root;
381			FVECT orig,dir;
382			int *nextptr;
383			FUNC func;
384			int *f;
385	gwlarson	3.4	{
386	gwlarson	3.8	double br[3];
387			BCOORD bi[3],dbi[3];
388
389			/* Project the origin onto the root node plane */
390			/* Find the intersection point of the origin */
391			ray_to_qt_frame(root,orig,dir,bi,dbi);
392	gwlarson	3.3
393	gwlarson	3.8	/* trace the ray starting with this node */
394			qtTrace_ray(qt,bi,dbi[0],dbi[1],dbi[2],nextptr,0,0,func,f);
395			if(!QT_FLAG_IS_DONE(*f))
396			qt_frame_to_vert(root,bi,orig);
397	gwlarson	3.6
398	gwlarson	3.4	}
399
400	gwlarson	3.6	/* Trace ray 'orig-dir' through stree and apply 'func(qtptr,f,argptr)' at each
401			node that it intersects
402			*/
403	gwlarson	3.4	int
404	gwlarson	3.8	stTrace_ray(st,orig,dir,func)
405	gwlarson	3.4	STREE *st;
406			FVECT orig,dir;
407	gwlarson	3.8	FUNC func;
408	gwlarson	3.4	{
409	gwlarson	3.6	int next,last,i,f=0;
410	gwlarson	3.8	QUADTREE qt;
411	gwlarson	3.7	FVECT o,n,v;
412			double pd,t,d;
413	gwlarson	3.4
414			VCOPY(o,orig);
415	gwlarson	3.7	#ifdef TEST_DRIVER
416			Pick_cnt=0;
417			#endif;
418	gwlarson	3.8	/* Find the qt node that o falls in */
419			i = stLocate_root(o);
420			qt = ST_ROOT_QT(st,i);
421	gwlarson	3.6
422	gwlarson	3.8	stRoot_trace_ray(qt,i,o,dir,&next,func,&f);
423	gwlarson	3.6
424			if(QT_FLAG_IS_DONE(f))
425			return(TRUE);
426	gwlarson	3.8	/*
427	gwlarson	3.7	d = DOT(orig,dir)/sqrt(DOT(orig,orig));
428			VSUM(v,orig,dir,-d);
429	gwlarson	3.8	*/
430	gwlarson	3.6	/* Crossed over to next cell: id = nbr */
431			while(1)
432			{
433			/* test if ray crosses plane between this quadtree triangle and
434			its neighbor- if it does then find intersection point with
435			ray and plane- this is the new origin
436			*/
437			if(next == INVALID)
438			return(FALSE);
439	gwlarson	3.8	/*
440			if(DOT(o,v) < 0.0)
441			return(FALSE);
442			*/
443	gwlarson	3.6	i = stBase_nbrs[i][next];
444	gwlarson	3.8	qt = ST_ROOT_QT(st,i);
445			stRoot_trace_ray(qt,i,o,dir,&next,func,&f);
446	gwlarson	3.6	if(QT_FLAG_IS_DONE(f))
447			return(TRUE);
448			}
449	gwlarson	3.4	}
450
451
452	gwlarson	3.8	stVisit_poly(st,verts,l,root,func)
453			STREE *st;
454			FVECT *verts;
455			LIST *l;
456			unsigned int root;
457			FUNC func;
458	gwlarson	3.4	{
459	gwlarson	3.8	int id0,id1,id2;
460			FVECT tri[3];
461	gwlarson	3.4
462	gwlarson	3.8	id0 = pop_list(&l);
463			id1 = pop_list(&l);
464			while(l)
465			{
466			id2 = pop_list(&l);
467			VCOPY(tri[0],verts[id0]);
468			VCOPY(tri[1],verts[id1]);
469			VCOPY(tri[2],verts[id2]);
470			stRoot_visit_tri(st,root,tri,func);
471			id1 = id2;
472			}
473			}
474	gwlarson	3.6
475	gwlarson	3.8	stVisit_clip(st,i,verts,vcnt,l,cell,func)
476			STREE *st;
477			int i;
478			FVECT *verts;
479			int *vcnt;
480			LIST *l;
481			unsigned int cell;
482			FUNC func;
483			{
484
485			LIST labove,lbelow,endb,enda;
486			int last = -1;
487			int id,last_id;
488			int first,first_id;
489			unsigned int cellb;
490
491			labove = lbelow = NULL;
492			enda = endb = NULL;
493			while(l)
494	gwlarson	3.4	{
495	gwlarson	3.8	id = pop_list(&l);
496			if(ZERO(verts[id][i]))
497			{
498			if(last==-1)
499			{/* add below and above */
500			first = 2;
501			first_id= id;
502			}
503			lbelow=add_data(lbelow,id,&endb);
504			labove=add_data(labove,id,&enda);
505			last_id = id;
506			last = 2;
507			continue;
508			}
509			if(verts[id][i] < 0)
510			{
511			if(last != 1)
512			{
513			lbelow=add_data(lbelow,id,&endb);
514			if(last==-1)
515			{
516			first = 0;
517			first_id = id;
518			}
519			last_id = id;
520			last = 0;
521			continue;
522			}
523			/* intersect_edges */
524			intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
525			/newpoint goes to above and below/
526			lbelow=add_data(lbelow,*vcnt,&endb);
527			lbelow=add_data(lbelow,id,&endb);
528			labove=add_data(labove,*vcnt,&enda);
529			last = 0;
530			last_id = id;
531			(*vcnt)++;
532			}
533			else
534			{
535			if(last != 0)
536			{
537			labove=add_data(labove,id,&enda);
538			if(last==-1)
539			{
540			first = 1;
541			first_id = id;
542			}
543			last_id = id;
544			last = 1;
545			continue;
546			}
547			/* intersect_edges */
548			/newpoint goes to above and below/
549			intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
550			lbelow=add_data(lbelow,*vcnt,&endb);
551			labove=add_data(labove,*vcnt,&enda);
552			labove=add_data(labove,id,&enda);
553			last_id = id;
554			(*vcnt)++;
555			last = 1;
556			}
557			}
558			if(first != 2 && first != last)
559			{
560			intersect_edge_coord_plane(verts[id],verts[first_id],i,verts[*vcnt]);
561			/newpoint goes to above and below/
562			lbelow=add_data(lbelow,*vcnt,&endb);
563			labove=add_data(labove,*vcnt,&enda);
564			(*vcnt)++;
565	gwlarson	3.6
566	gwlarson	3.4	}
567	gwlarson	3.8	if(i==2)
568			{
569			if(lbelow)
570			{
571			if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
572			{
573			cellb = cell \| (1 << i);
574			stVisit_poly(st,verts,lbelow,cellb,func);
575			}
576			else
577			free_list(lbelow);
578			}
579			if(labove)
580			{
581			if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
582			stVisit_poly(st,verts,labove,cell,func);
583			else
584			free_list(labove);
585			}
586			}
587			else
588			{
589			if(lbelow)
590			{
591			if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
592			{
593			cellb = cell \| (1 << i);
594			stVisit_clip(st,i+1,verts,vcnt,lbelow,cellb,func);
595			}
596			else
597			free_list(lbelow);
598			}
599			if(labove)
600			{
601			if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
602			stVisit_clip(st,i+1,verts,vcnt,labove,cell,func);
603			else
604			free_list(labove);
605			}
606			}
607
608	gwlarson	3.4	}
609
610	gwlarson	3.8	stVisit(st,tri,func)
611	gwlarson	3.4	STREE *st;
612	gwlarson	3.8	FVECT tri[3];
613			FUNC func;
614	gwlarson	3.4	{
615	gwlarson	3.8	int r0,r1,r2;
616			LIST *l;
617	gwlarson	3.7
618	gwlarson	3.8	r0 = stLocate_root(tri[0]);
619			r1 = stLocate_root(tri[1]);
620			r2 = stLocate_root(tri[2]);
621			if(r0 == r1 && r1==r2)
622			stRoot_visit_tri(st,r0,tri,func);
623			else
624			{
625			FVECT verts[ST_CLIP_VERTS];
626			int cnt;
627	gwlarson	3.3
628	gwlarson	3.8	VCOPY(verts[0],tri[0]);
629			VCOPY(verts[1],tri[1]);
630			VCOPY(verts[2],tri[2]);
631
632			l = add_data(NULL,0,NULL);
633			l = add_data(l,1,NULL);
634			l = add_data(l,2,NULL);
635			cnt = 3;
636			stVisit_clip(st,0,verts,&cnt,l,0,func);
637			}
638	gwlarson	3.3	}
639	gwlarson	3.7
640
641	gwlarson	3.8	/* New Insertion code!!! */
642
643
644			BCOORD qtRoot[3][3] = { {MAXBCOORD2,0,0},{0,MAXBCOORD2,0},{0,0,MAXBCOORD2}};
645
646
647			convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02)
648			int root;
649			FVECT tri[3];
650			BCOORD b0[3],b1[3],b2[3];
651			BCOORD db10[3],db21[3],db02[3];
652			{
653			/* Project the vertex into the qtree plane */
654			vert_to_qt_frame(root,tri[0],b0);
655			vert_to_qt_frame(root,tri[1],b1);
656			vert_to_qt_frame(root,tri[2],b2);
657
658			/* calculate triangle edge differences in new frame */
659			db10[0] = b1[0] - b0[0]; db10[1] = b1[1] - b0[1]; db10[2] = b1[2] - b0[2];
660			db21[0] = b2[0] - b1[0]; db21[1] = b2[1] - b1[1]; db21[2] = b2[2] - b1[2];
661			db02[0] = b0[0] - b2[0]; db02[1] = b0[1] - b2[1]; db02[2] = b0[2] - b2[2];
662			}
663
664
665			QUADTREE
666			stRoot_insert_tri(st,root,tri,f)
667			STREE *st;
668			int root;
669			FVECT tri[3];
670			FUNC f;
671			{
672			BCOORD b0[3],b1[3],b2[3];
673			BCOORD db10[3],db21[3],db02[3];
674			unsigned int s0,s1,s2,sq0,sq1,sq2;
675			QUADTREE qt;
676
677			/* Map the triangle vertices into the canonical barycentric frame */
678			convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);
679
680			/* Calculate initial sidedness info */
681			SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
682			SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);
683
684			qt = ST_ROOT_QT(st,root);
685			/* Visit cells that triangle intersects */
686			qt = qtInsert_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
687			b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f,0);
688
689			return(qt);
690			}
691
692			stRoot_visit_tri(st,root,tri,f)
693			STREE *st;
694			int root;
695			FVECT tri[3];
696			FUNC f;
697			{
698			BCOORD b0[3],b1[3],b2[3];
699			BCOORD db10[3],db21[3],db02[3];
700			unsigned int s0,s1,s2,sq0,sq1,sq2;
701			QUADTREE qt;
702
703			/* Map the triangle vertices into the canonical barycentric frame */
704			convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);
705
706			/* Calculate initial sidedness info */
707			SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
708			SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);
709
710			qt = ST_ROOT_QT(st,root);
711			QT_SET_FLAG(ST_QT(st,root));
712			/* Visit cells that triangle intersects */
713			qtVisit_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
714			b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f);
715
716			}
717
718			stInsert_tri(st,tri,f)
719			STREE *st;
720			FVECT tri[3];
721			FUNC f;
722			{
723			unsigned int cells,which;
724			int root;
725
726
727			/* calculate entry/exit points of edges through the cells */
728			cells = stTri_cells(st,tri);
729
730			/* For each cell that quadtree intersects: Map the triangle vertices into
731			the canonical barycentric frame of (1,0,0), (0,1,0),(0,0,1). Insert
732			by first doing a trivial reject on the interior nodes, and then a
733			tri/tri intersection at the leaf nodes.
734			*/
735			for(root=0,which=1; root < ST_NUM_ROOT_NODES; root++,which <<= 1)
736			{
737			/* For each of the quadtree roots: check if marked as intersecting tri*/
738			if(cells & which)
739			/* Visit tri cells */
740			ST_ROOT_QT(st,root) = stRoot_insert_tri(st,root,tri,f);
741			}
742			}
743	gwlarson	3.7
744
745	gwlarson	3.6
746
747
748
749	gwlarson	3.5
750	gwlarson	3.3
751	gwlarson	3.4
752
753