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;

  qtDone();

  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));
}

stFree(st)
STREE *st;
{
  qtDone();
  free(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);
  
  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);
  }
}

stInsert_samp(st,p,f)
   STREE *st;
   FVECT p;
   FUNC f;
{

    QUADTREE qt;
    BCOORD bcoordi[3];
    int i,done;

    /* 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
    */
    vert_to_qt_frame(i,p,bcoordi);
    ST_ROOT_QT(st,i) =  qtInsert_point(i,qt,EMPTY,qtRoot[0],qtRoot[1],
                          qtRoot[2],bcoordi,MAXBCOORD2>>1,f,0,&done);

}


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