src/hd/sm_stree.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#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 "object.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;
  }
}
/* Assumption: know crosses plane:dont need to check for 'on' case */
intersect_edge_coord_plane(v0,v1,w,r)
FVECT v0,v1;
int w;
FVECT r;
{
  FVECT dv;
  int wnext;
  double t;

  VSUB(dv,v1,v0);
  t = -v0[w]/dv[w];
  r[w] = 0.0;
  wnext = (w+1)%3;
  r[wnext] = v0[wnext] + dv[wnext]*t;
  wnext = (w+2)%3;
  r[wnext] = v0[wnext] + dv[wnext]*t;
}


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


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