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,n)
STREE *st;
FVECT *verts;
LIST *l;
unsigned int root;
FUNC func;
int n;
{
  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,n);
    id1 = id2;
  }
}

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

  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,n);
      }
      else
        free_list(lbelow);
    }
    if(labove)
     {
      if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
        stVisit_poly(st,verts,labove,cell,func,n);
      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,n);
        }
      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,n);
       else
         free_list(labove);
     }
  }

}

stVisit(st,tri,func,n)
   STREE *st;
   FVECT tri[3];
   FUNC func;
   int n;
{
    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,n);
    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,n);
      }
}


/* 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,n)
   STREE *st;
   int root;
   FVECT tri[3];
   FUNC f;
   int n;
{
  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,n);

}

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.11
Committed:	Fri Mar 5 16:32:50 1999 UTC (25 years, 1 month ago) by gwlarson
Content type:	text/plain
Branch:	MAIN
Changes since 3.10:	+17 -13 lines
Log Message:	added variable level traversal to node visiting routines to support display list 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.11	/*
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.11	stVisit_poly(st,verts,l,root,func,n)
453	gwlarson	3.8	STREE *st;
454			FVECT *verts;
455			LIST *l;
456			unsigned int root;
457			FUNC func;
458	gwlarson	3.11	int n;
459	gwlarson	3.4	{
460	gwlarson	3.8	int id0,id1,id2;
461			FVECT tri[3];
462	gwlarson	3.4
463	gwlarson	3.8	id0 = pop_list(&l);
464			id1 = pop_list(&l);
465			while(l)
466			{
467			id2 = pop_list(&l);
468			VCOPY(tri[0],verts[id0]);
469			VCOPY(tri[1],verts[id1]);
470			VCOPY(tri[2],verts[id2]);
471	gwlarson	3.11	stRoot_visit_tri(st,root,tri,func,n);
472	gwlarson	3.8	id1 = id2;
473			}
474			}
475	gwlarson	3.6
476	gwlarson	3.11	stVisit_clip(st,i,verts,vcnt,l,cell,func,n)
477	gwlarson	3.8	STREE *st;
478			int i;
479			FVECT *verts;
480			int *vcnt;
481			LIST *l;
482			unsigned int cell;
483			FUNC func;
484	gwlarson	3.11	int n;
485	gwlarson	3.8	{
486
487			LIST labove,lbelow,endb,enda;
488			int last = -1;
489			int id,last_id;
490			int first,first_id;
491			unsigned int cellb;
492
493			labove = lbelow = NULL;
494			enda = endb = NULL;
495			while(l)
496	gwlarson	3.4	{
497	gwlarson	3.8	id = pop_list(&l);
498			if(ZERO(verts[id][i]))
499			{
500			if(last==-1)
501			{/* add below and above */
502			first = 2;
503			first_id= id;
504			}
505			lbelow=add_data(lbelow,id,&endb);
506			labove=add_data(labove,id,&enda);
507			last_id = id;
508			last = 2;
509			continue;
510			}
511			if(verts[id][i] < 0)
512			{
513			if(last != 1)
514			{
515			lbelow=add_data(lbelow,id,&endb);
516			if(last==-1)
517			{
518			first = 0;
519			first_id = id;
520			}
521			last_id = id;
522			last = 0;
523			continue;
524			}
525			/* intersect_edges */
526			intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
527			/newpoint goes to above and below/
528			lbelow=add_data(lbelow,*vcnt,&endb);
529			lbelow=add_data(lbelow,id,&endb);
530			labove=add_data(labove,*vcnt,&enda);
531			last = 0;
532			last_id = id;
533			(*vcnt)++;
534			}
535			else
536			{
537			if(last != 0)
538			{
539			labove=add_data(labove,id,&enda);
540			if(last==-1)
541			{
542			first = 1;
543			first_id = id;
544			}
545			last_id = id;
546			last = 1;
547			continue;
548			}
549			/* intersect_edges */
550			/newpoint goes to above and below/
551			intersect_edge_coord_plane(verts[last_id],verts[id],i,verts[*vcnt]);
552			lbelow=add_data(lbelow,*vcnt,&endb);
553			labove=add_data(labove,*vcnt,&enda);
554			labove=add_data(labove,id,&enda);
555			last_id = id;
556			(*vcnt)++;
557			last = 1;
558			}
559			}
560			if(first != 2 && first != last)
561			{
562			intersect_edge_coord_plane(verts[id],verts[first_id],i,verts[*vcnt]);
563			/newpoint goes to above and below/
564			lbelow=add_data(lbelow,*vcnt,&endb);
565			labove=add_data(labove,*vcnt,&enda);
566			(*vcnt)++;
567	gwlarson	3.6
568	gwlarson	3.4	}
569	gwlarson	3.8	if(i==2)
570			{
571			if(lbelow)
572			{
573			if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
574			{
575			cellb = cell \| (1 << i);
576	gwlarson	3.11	stVisit_poly(st,verts,lbelow,cellb,func,n);
577	gwlarson	3.8	}
578			else
579			free_list(lbelow);
580			}
581			if(labove)
582			{
583			if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
584	gwlarson	3.11	stVisit_poly(st,verts,labove,cell,func,n);
585	gwlarson	3.8	else
586			free_list(labove);
587			}
588			}
589			else
590			{
591			if(lbelow)
592			{
593			if(LIST_NEXT(lbelow) && LIST_NEXT(LIST_NEXT(lbelow)))
594			{
595			cellb = cell \| (1 << i);
596	gwlarson	3.11	stVisit_clip(st,i+1,verts,vcnt,lbelow,cellb,func,n);
597	gwlarson	3.8	}
598			else
599			free_list(lbelow);
600			}
601			if(labove)
602			{
603			if(LIST_NEXT(labove) && LIST_NEXT(LIST_NEXT(labove)))
604	gwlarson	3.11	stVisit_clip(st,i+1,verts,vcnt,labove,cell,func,n);
605	gwlarson	3.8	else
606			free_list(labove);
607			}
608			}
609
610	gwlarson	3.4	}
611
612	gwlarson	3.11	stVisit(st,tri,func,n)
613	gwlarson	3.4	STREE *st;
614	gwlarson	3.8	FVECT tri[3];
615			FUNC func;
616	gwlarson	3.11	int n;
617	gwlarson	3.4	{
618	gwlarson	3.8	int r0,r1,r2;
619			LIST *l;
620	gwlarson	3.7
621	gwlarson	3.8	r0 = stLocate_root(tri[0]);
622			r1 = stLocate_root(tri[1]);
623			r2 = stLocate_root(tri[2]);
624			if(r0 == r1 && r1==r2)
625	gwlarson	3.11	stRoot_visit_tri(st,r0,tri,func,n);
626	gwlarson	3.8	else
627			{
628			FVECT verts[ST_CLIP_VERTS];
629			int cnt;
630	gwlarson	3.3
631	gwlarson	3.8	VCOPY(verts[0],tri[0]);
632			VCOPY(verts[1],tri[1]);
633			VCOPY(verts[2],tri[2]);
634
635			l = add_data(NULL,0,NULL);
636			l = add_data(l,1,NULL);
637			l = add_data(l,2,NULL);
638			cnt = 3;
639	gwlarson	3.11	stVisit_clip(st,0,verts,&cnt,l,0,func,n);
640	gwlarson	3.8	}
641	gwlarson	3.3	}
642	gwlarson	3.7
643
644	gwlarson	3.8	/* New Insertion code!!! */
645
646
647			BCOORD qtRoot[3][3] = { {MAXBCOORD2,0,0},{0,MAXBCOORD2,0},{0,0,MAXBCOORD2}};
648
649
650			convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02)
651			int root;
652			FVECT tri[3];
653			BCOORD b0[3],b1[3],b2[3];
654			BCOORD db10[3],db21[3],db02[3];
655			{
656			/* Project the vertex into the qtree plane */
657			vert_to_qt_frame(root,tri[0],b0);
658			vert_to_qt_frame(root,tri[1],b1);
659			vert_to_qt_frame(root,tri[2],b2);
660
661			/* calculate triangle edge differences in new frame */
662			db10[0] = b1[0] - b0[0]; db10[1] = b1[1] - b0[1]; db10[2] = b1[2] - b0[2];
663			db21[0] = b2[0] - b1[0]; db21[1] = b2[1] - b1[1]; db21[2] = b2[2] - b1[2];
664			db02[0] = b0[0] - b2[0]; db02[1] = b0[1] - b2[1]; db02[2] = b0[2] - b2[2];
665			}
666
667
668			QUADTREE
669			stRoot_insert_tri(st,root,tri,f)
670			STREE *st;
671			int root;
672			FVECT tri[3];
673			FUNC f;
674			{
675			BCOORD b0[3],b1[3],b2[3];
676			BCOORD db10[3],db21[3],db02[3];
677			unsigned int s0,s1,s2,sq0,sq1,sq2;
678			QUADTREE qt;
679
680			/* Map the triangle vertices into the canonical barycentric frame */
681			convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);
682
683			/* Calculate initial sidedness info */
684			SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
685			SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);
686
687			qt = ST_ROOT_QT(st,root);
688			/* Visit cells that triangle intersects */
689			qt = qtInsert_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
690			b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f,0);
691
692			return(qt);
693			}
694
695	gwlarson	3.11	stRoot_visit_tri(st,root,tri,f,n)
696	gwlarson	3.8	STREE *st;
697			int root;
698			FVECT tri[3];
699			FUNC f;
700	gwlarson	3.11	int n;
701	gwlarson	3.8	{
702			BCOORD b0[3],b1[3],b2[3];
703			BCOORD db10[3],db21[3],db02[3];
704			unsigned int s0,s1,s2,sq0,sq1,sq2;
705			QUADTREE qt;
706
707			/* Map the triangle vertices into the canonical barycentric frame */
708			convert_tri_to_frame(root,tri,b0,b1,b2,db10,db21,db02);
709
710			/* Calculate initial sidedness info */
711			SIDES_GTR(b0,b1,b2,s0,s1,s2,qtRoot[1][0],qtRoot[0][1],qtRoot[0][2]);
712			SIDES_GTR(b0,b1,b2,sq0,sq1,sq2,qtRoot[0][0],qtRoot[1][1],qtRoot[2][2]);
713
714			qt = ST_ROOT_QT(st,root);
715			QT_SET_FLAG(ST_QT(st,root));
716			/* Visit cells that triangle intersects */
717			qtVisit_tri(root,qt,qtRoot[0],qtRoot[1],qtRoot[2],
718	gwlarson	3.11	b0,b1,b2,db10,db21,db02,MAXBCOORD2 >> 1,s0,s1,s2, sq0,sq1,sq2,f,n);
719	gwlarson	3.8
720			}
721
722			stInsert_tri(st,tri,f)
723			STREE *st;
724			FVECT tri[3];
725			FUNC f;
726			{
727			unsigned int cells,which;
728			int root;
729
730
731			/* calculate entry/exit points of edges through the cells */
732			cells = stTri_cells(st,tri);
733
734			/* For each cell that quadtree intersects: Map the triangle vertices into
735			the canonical barycentric frame of (1,0,0), (0,1,0),(0,0,1). Insert
736			by first doing a trivial reject on the interior nodes, and then a
737			tri/tri intersection at the leaf nodes.
738			*/
739			for(root=0,which=1; root < ST_NUM_ROOT_NODES; root++,which <<= 1)
740			{
741			/* For each of the quadtree roots: check if marked as intersecting tri*/
742			if(cells & which)
743			/* Visit tri cells */
744			ST_ROOT_QT(st,root) = stRoot_insert_tri(st,root,tri,f);
745			}
746			}
747	gwlarson	3.7
748
749	gwlarson	3.6
750
751
752
753	gwlarson	3.5
754	gwlarson	3.3
755	gwlarson	3.4
756
757