src/hd/sm_qtree.c

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

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

/*
*  sm_qtree.c: adapted from octree.c from radiance code
*/
/*
 *  octree.c - routines dealing with octrees and cubes.
 *
 *     7/28/85
 */

#include "standard.h"

#include "sm_geom.h"
#include "sm_qtree.h"
#include "object.h"

QUADTREE  *quad_block[QT_MAX_BLK];              /* our quadtree */
static QUADTREE  quad_free_list = EMPTY;  /* freed octree nodes */
static QUADTREE  treetop = 0;           /* next free node */
int4 *quad_flag;

QUADTREE
qtAlloc()                       /* allocate a quadtree */
{
    register QUADTREE  freet;

    if ((freet = quad_free_list) != EMPTY)
    {
        quad_free_list = QT_NTH_CHILD(freet, 0);
        return(freet);
    }
    freet = treetop;
    if (QT_BLOCK_INDEX(freet) == 0)
    {
        if (QT_BLOCK(freet) >= QT_MAX_BLK)
           return(EMPTY);
        if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
                        QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
           return(EMPTY);
        quad_flag = (int4 *)realloc((char *)quad_flag,
                        (QT_BLOCK(freet)+1)*QT_BLOCK_SIZE/(8*sizeof(int4)));
        if (quad_flag == NULL)
                return(EMPTY);
    }
    treetop += 4;
    return(freet);
}


qtClearAllFlags()               /* clear all quadtree branch flags */
{
        if (!treetop) return;
        bzero((char *)quad_flag,
                (QT_BLOCK(treetop-1)+1)*QT_BLOCK_SIZE/(8*sizeof(int4)));
}


qtFree(qt)                      /* free a quadtree */
   register QUADTREE  qt;
{
        register int  i;

        if (!QT_IS_TREE(qt))
        {
          qtfreeleaf(qt);
          return;
        }
        for (i = 0; i < 4; i++)
           qtFree(QT_NTH_CHILD(qt, i));
        QT_NTH_CHILD(qt, 0) = quad_free_list;
        quad_free_list = qt;
}


qtDone()                        /* free EVERYTHING */
{
        register int    i;

        for (i = 0; i < QT_MAX_BLK; i++)
        {
            if (quad_block[i] == NULL)
                break;
            free((char *)quad_block[i]);
            quad_block[i] = NULL;
        }
        if (i) free((char *)quad_flag);
        quad_flag = NULL;
        quad_free_list = EMPTY;
        treetop = 0;
}


QUADTREE
qtCompress(qt)                  /* recursively combine nodes */
register QUADTREE  qt;
{
    register int  i;
    register QUADTREE  qres;

    if (!QT_IS_TREE(qt))        /* not a tree */
       return(qt);
    qres = QT_NTH_CHILD(qt,0) = qtCompress(QT_NTH_CHILD(qt,0));
    for (i = 1; i < 4; i++)
       if((QT_NTH_CHILD(qt,i) = qtCompress(QT_NTH_CHILD(qt,i))) != qres)
          qres = qt;
    if(!QT_IS_TREE(qres))
    {   /* all were identical leaves */
        QT_NTH_CHILD(qt,0) = quad_free_list;
        quad_free_list = qt;
    }
    return(qres);
}


QUADTREE
*qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
QUADTREE *qtptr;
double bcoord[3];
FVECT t0,t1,t2;
{
  int i;
  QUADTREE *child;
  FVECT a,b,c;

    if(QT_IS_TREE(*qtptr))
    {  
      i = bary2d_child(bcoord);
      child = QT_NTH_CHILD_PTR(*qtptr,i);
      if(t0)
      {
        qtSubdivide_tri(t0,t1,t2,a,b,c);
        qtNth_child_tri(t0,t1,t2,a,b,c,i,t0,t1,t2);
      }
      return(qtLocate_leaf(child,bcoord,t0,t1,t2));
    }
    else
      return(qtptr);
}


int
qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n)
QUADTREE *qtptr;
double bcoord[3];
int id;
FVECT v0,v1,v2;
int n;
{
  int i;
  QUADTREE *child;
  OBJECT os[MAXSET+1],*optr;
  int found;
  FVECT r0,r1,r2;
  
  if(QT_IS_TREE(*qtptr))
  {  
    QT_SET_FLAG(*qtptr);
    i = bary2d_child(bcoord);
    child = QT_NTH_CHILD_PTR(*qtptr,i);
    return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,++n));
  }
  else
    {
      /* If this leave node emptry- create a new set */
      if(QT_IS_EMPTY(*qtptr))
        *qtptr = qtaddelem(*qtptr,id);
      else
      {
          qtgetset(os,*qtptr);
          /* If the set is too large: subdivide */
          if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
            *qtptr = qtaddelem(*qtptr,id);
          else
          {
            if (n < QT_MAX_LEVELS)
            {
                 /* If set size exceeds threshold: subdivide cell and
                    reinsert set tris into cell
                    */
                 n++;
                 qtfreeleaf(*qtptr);
                 qtSubdivide(qtptr);
                 QT_SET_FLAG(*qtptr);
                 found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n);
#if 0
                 if(!found)
                 {
#ifdef TEST_DRIVER     
      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
#else 
                   eputs("qtAdd_tri():Found in parent but not children\n");
#endif
                 }
#endif
                 for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
                 {
                   id = QT_SET_NEXT_ELEM(optr);
                   qtTri_verts_from_id(id,r0,r1,r2);
                   found= qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,++n);
#ifdef DEBUG
                 if(!found)
                    eputs("qtAdd_tri():Reinsert-in parent but not children\n");
#endif
               }
             }
            else
              if(QT_SET_CNT(os) < QT_MAX_SET)
                {
                  *qtptr = qtaddelem(*qtptr,id);
                }
              else
                {
#ifdef DEBUG
                    eputs("qtAdd_tri():two many levels\n");
#endif
                    return(FALSE);
                }
          }
      }
  }
  return(TRUE);
}


int
qtAdd_tri_from_point(qtptr,v0,v1,v2,pt,id)
QUADTREE *qtptr;
FVECT v0,v1,v2;
FVECT pt;
int id;
{
    char d,w;
    int i,x,y;
    QUADTREE *child;
    QUADTREE qt;
    FVECT i_pt,n,a,b,c,r0,r1,r2;
    double pd,bcoord[3];
    OBJECT os[MAXSET+1],*optr;
    int found;
        
    /* Determine if point lies within pyramid (and therefore
       inside a spherical quadtree cell):GT_INTERIOR, on one of the 
       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3), 
       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
       For each triangle edge: compare the
       point against the plane formed by the edge and the view center
    */
    d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);  

    /* Not in this triangle */
    if(!d)
      return(FALSE);

    /* 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(*qtptr))
    {  
      QT_SET_FLAG(*qtptr);
      /* Find the intersection point */
      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
      intersect_vector_plane(pt,n,pd,NULL,i_pt);

      i = max_index(n);
      x = (i+1)%3;
      y = (i+2)%3;
      /* Calculate barycentric coordinates of i_pt */
      bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);

      i = bary2d_child(bcoord);
      child = QT_NTH_CHILD_PTR(*qtptr,i);
      /* NOTE: Optimize: only subdivide for specified child */
      qtSubdivide_tri(v0,v1,v2,a,b,c);
      qtNth_child_tri(v0,v1,v2,a,b,c,i,v0,v1,v2);
      return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,1));
    }
    else
   {
      /* If this leave node emptry- create a new set */
      if(QT_IS_EMPTY(*qtptr))
        *qtptr = qtaddelem(*qtptr,id);
      else
      {
          qtgetset(os,*qtptr);
          /* If the set is too large: subdivide */
          if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
            *qtptr = qtaddelem(*qtptr,id);
          else
          {
                 /* If set size exceeds threshold: subdivide cell and
                    reinsert set tris into cell
                    */
                 qtfreeleaf(*qtptr);
                 qtSubdivide(qtptr);
                 found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,1);
#if 0
                 if(!found)
                 {
#ifdef TEST_DRIVER     
      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
#else 
                   eputs("qtAdd_tri():Found in parent but not children\n");
#endif
                 }
#endif
                 for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
                 {
                   id = QT_SET_NEXT_ELEM(optr);
                   qtTri_verts_from_id(id,r0,r1,r2);
                   found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,1);
#ifdef DEBUG
                 if(!found)
                    eputs("qtAdd_tri():Reinsert-in parent but not children\n");
#endif
               }
             }
      }
  }
  return(TRUE);
}


QUADTREE
*qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
QUADTREE *qtptr;
FVECT v0,v1,v2;
FVECT pt;
FVECT t0,t1,t2;
{
    char d,w;
    int i,x,y;
    QUADTREE *child;
    QUADTREE qt;
    FVECT n,i_pt,a,b,c;
    double pd,bcoord[3];

    /* Determine if point lies within pyramid (and therefore
       inside a spherical quadtree cell):GT_INTERIOR, on one of the 
       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3), 
       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
       For each triangle edge: compare the
       point against the plane formed by the edge and the view center
    */
    d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);  

    /* Not in this triangle */
    if(!d)
    {
      return(NULL);
    }

    /* 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(*qtptr))
    {  
      /* Find the intersection point */
      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
      intersect_vector_plane(pt,n,pd,NULL,i_pt);

      i = max_index(n);
      x = (i+1)%3;
      y = (i+2)%3;
      /* Calculate barycentric coordinates of i_pt */
      bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);

      i = bary2d_child(bcoord);
      child = QT_NTH_CHILD_PTR(*qtptr,i);
      if(t0)
      {
        qtSubdivide_tri(v0,v1,v2,a,b,c);
        qtNth_child_tri(v0,v1,v2,a,b,c,i,t0,t1,t2);
      }
      return(qtLocate_leaf(child,bcoord,t0,t1,t2));
    }
    else
      return(qtptr);
}

int
qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
QUADTREE *qtptr;
FVECT v0,v1,v2;
FVECT pt;
char *type,*which;
{
    OBJECT os[MAXSET+1],*optr;
    char d,w;
    int i,id;
    FVECT p0,p1,p2;

    qtptr = qtRoot_point_locate(qtptr,v0,v1,v2,pt,NULL,NULL,NULL);

    if(!qtptr || QT_IS_EMPTY(*qtptr))
       return(EMPTY);
    
    /* Get the set */
    qtgetset(os,*qtptr);
    for (i = QT_SET_CNT(os),optr = QT_SET_PTR(os); i > 0; i--)
    {
        /* Find the triangle that pt falls in (NOTE:FOR now return first 1) */
        id = QT_SET_NEXT_ELEM(optr);
        qtTri_verts_from_id(id,p0,p1,p2);
        d = test_single_point_against_spherical_tri(p0,p1,p2,pt,&w);  
        if(d)
        {
          if(type)
            *type = d;
          if(which)
            *which = w;
          return(id);
        }
    }
    return(EMPTY);
}

QUADTREE
qtSubdivide(qtptr)
QUADTREE *qtptr;
{
  QUADTREE node;
  node = qtAlloc();
  QT_CLEAR_CHILDREN(node);
  *qtptr = node;
  return(node);
}


QUADTREE
qtSubdivide_nth_child(qt,n)
QUADTREE qt;
int n;
{
  QUADTREE node;

  node = qtSubdivide(&(QT_NTH_CHILD(qt,n)));
  
  return(node);
}

/* for triangle v0-v1-v2- returns a,b,c: children are:

        v2                     0: v0,a,c
        /\                     1: a,v1,b                  
       /2 \                    2: c,b,v2
     c/____\b                  3: b,c,a
     /\    /\  
    /0 \3 /1 \
  v0____\/____\v1
        a
 */

qtSubdivide_tri(v0,v1,v2,a,b,c)
FVECT v0,v1,v2;
FVECT a,b,c;
{
    EDGE_MIDPOINT_VEC3(a,v0,v1);
    EDGE_MIDPOINT_VEC3(b,v1,v2);
    EDGE_MIDPOINT_VEC3(c,v2,v0);
}

qtNth_child_tri(v0,v1,v2,a,b,c,i,r0,r1,r2)
FVECT v0,v1,v2;
FVECT a,b,c;
int i;
FVECT r0,r1,r2;
{
  switch(i){
  case 0:  
  VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
    break;
  case 1:  
  VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
    break;
  case 2:  
  VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
    break;
  case 3:  
  VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
    break;
  }
}

/* Add triangle "id" to all leaf level cells that are children of 
quadtree pointed to by "qtptr" with cell vertices "t1,t2,t3"
that it overlaps (vertex and edge adjacencies do not count
as overlapping). If the addition of the triangle causes the cell to go over
threshold- the cell is split- and the triangle must be recursively inserted
into the new child cells: it is assumed that "v1,v2,v3" are normalized
*/

int
qtRoot_add_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
QUADTREE *qtptr;
int id;
FVECT t0,t1,t2;
FVECT v0,v1,v2;
int n;
{
  char test;
  int found;

  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
  if(!test)
    return(FALSE);
  
  found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
  return(found);
}

int
qtRoot_add_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n)
QUADTREE *qtptr;
int id;
FVECT t0,t1,t2;
FVECT v0,v1,v2;
int n;
{
  char test;
  int found;

  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
  if(!test)
    return(FALSE);
  
  found = qtAdd_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n);
  return(found);
}

int
qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
QUADTREE *qtptr;
int id;
FVECT t0,t1,t2;
FVECT v0,v1,v2;
int n;
{
  
  char test;
  int i,index;
  FVECT a,b,c;
  OBJECT os[MAXSET+1],*optr;
  QUADTREE qt;
  int found;
  FVECT r0,r1,r2;

  found = 0;
  /* if this is tree: recurse */
  if(QT_IS_TREE(*qtptr))
  {
    n++;
    qtSubdivide_tri(v0,v1,v2,a,b,c);
    test = spherical_tri_intersect(t0,t1,t2,v0,a,c);
    if(test)
      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c,n);
    test = spherical_tri_intersect(t0,t1,t2,a,v1,b);
    if(test)
      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b,n);
    test = spherical_tri_intersect(t0,t1,t2,c,b,v2);
    if(test)
      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2,n);
    test = spherical_tri_intersect(t0,t1,t2,b,c,a);
    if(test)
      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a,n);

#if 0
    if(!found)
    {
#ifdef TEST_DRIVER     
      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
#else 
      eputs("qtAdd_tri():Found in parent but not children\n");
#endif
    }
#endif
  }
  else
  {
      /* If this leave node emptry- create a new set */
      if(QT_IS_EMPTY(*qtptr))
        *qtptr = qtaddelem(*qtptr,id);
      else
      {
          qtgetset(os,*qtptr);
          /* If the set is too large: subdivide */
          if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
            *qtptr = qtaddelem(*qtptr,id);
          else
          {
            if (n < QT_MAX_LEVELS)
            {
                 /* If set size exceeds threshold: subdivide cell and
                    reinsert set tris into cell
                    */
                 n++;
                 qtfreeleaf(*qtptr);
                 qtSubdivide(qtptr);
                 found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
#if 0
                 if(!found)
                 {
#ifdef TEST_DRIVER     
      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
#else 
                   eputs("qtAdd_tri():Found in parent but not children\n");
#endif
                 }
#endif
                 for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
                 {
                   id = QT_SET_NEXT_ELEM(optr);
                   qtTri_verts_from_id(id,r0,r1,r2);
                   found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,n);
#ifdef DEBUG
                 if(!found)
                    eputs("qtAdd_tri():Reinsert-in parent but not children\n");
#endif
               }
             }
            else
              if(QT_SET_CNT(os) < QT_MAX_SET)
                {
                  *qtptr = qtaddelem(*qtptr,id);
#if 0
                  {
              int k;
              qtgetset(os,*qtptr);
              printf("\n%d:\n",os[0]);
              for(k=1; k <= os[0];k++)
                 printf("%d ",os[k]);
              printf("\n");
          }
#endif            
                  /*
                    insertelem(os,id);
                    *qtptr = fullnode(os);
                  */
                }
              else
                {
#ifdef DEBUG
                    eputs("qtAdd_tri():two many levels\n");
#endif
                    return(FALSE);
                }
          }
      }
  }
  return(TRUE);
}


int
qtApply_to_tri_cells(qtptr,t0,t1,t2,v0,v1,v2,func,arg)
QUADTREE *qtptr;
FVECT t0,t1,t2;
FVECT v0,v1,v2;
int (*func)();
char *arg;
{
  char test;
  FVECT a,b,c;

  /* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);

  /* If triangles do not overlap: done */
  if(!test)
    return(FALSE);

  /* if this is tree: recurse */
  if(QT_IS_TREE(*qtptr))
  {
    qtSubdivide_tri(v0,v1,v2,a,b,c);
    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t0,t1,t2,v0,a,c,func,arg);
    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t0,t1,t2,a,v1,b,func,arg);
    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t0,t1,t2,c,b,v2,func,arg);
    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t0,t1,t2,b,c,a,func,arg);
  }
  else
    return(func(qtptr,arg));
}


int
qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
QUADTREE *qtptr;
int id;
FVECT t0,t1,t2;
FVECT v0,v1,v2;
{
  
  char test;
  int i;
  FVECT a,b,c;
  OBJECT os[MAXSET+1];
  
  /* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);

  /* If triangles do not overlap: done */
  if(!test)
    return(FALSE);

  /* if this is tree: recurse */
  if(QT_IS_TREE(*qtptr))
  {
    qtSubdivide_tri(v0,v1,v2,a,b,c);
    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c);
    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b);
    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2);
    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a);
  }
  else
  {
      if(QT_IS_EMPTY(*qtptr))
      {
#ifdef DEBUG     
          eputs("qtRemove_tri(): triangle not found\n");
#endif
      }
      /* remove id from set */
      else
      {
          qtgetset(os,*qtptr);
          if(!inset(os,id))
          {
#ifdef DEBUG          
              eputs("qtRemove_tri(): tri not in set\n");
#endif
          }
          else
          {
            *qtptr = qtdelelem(*qtptr,id);
        }
      }
  }
  return(TRUE);
}


Revision:	3.3
Committed:	Tue Aug 25 11:03:27 1998 UTC (25 years, 8 months ago) by gwlarson
Content type:	text/plain
Branch:	MAIN
Changes since 3.2:	+280 -91 lines
Log Message:	fixed problem with picking (ray tracking) of tetrahedron
#	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_qtree.c: adapted from octree.c from radiance code
9			*/
10			/*
11			* octree.c - routines dealing with octrees and cubes.
12			*
13			* 7/28/85
14			*/
15
16			#include "standard.h"
17
18			#include "sm_geom.h"
19			#include "sm_qtree.h"
20			#include "object.h"
21
22			QUADTREE quad_block[QT_MAX_BLK]; / our quadtree */
23			static QUADTREE quad_free_list = EMPTY; /* freed octree nodes */
24			static QUADTREE treetop = 0; /* next free node */
25	gwlarson	3.3	int4 *quad_flag;
26	gwlarson	3.1
27			QUADTREE
28			qtAlloc() /* allocate a quadtree */
29			{
30			register QUADTREE freet;
31
32			if ((freet = quad_free_list) != EMPTY)
33			{
34			quad_free_list = QT_NTH_CHILD(freet, 0);
35			return(freet);
36			}
37			freet = treetop;
38			if (QT_BLOCK_INDEX(freet) == 0)
39			{
40			if (QT_BLOCK(freet) >= QT_MAX_BLK)
41			return(EMPTY);
42			if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
43	gwlarson	3.3	QT_BLOCK_SIZE4sizeof(QUADTREE))) == NULL)
44	gwlarson	3.1	return(EMPTY);
45	gwlarson	3.3	quad_flag = (int4 )realloc((char )quad_flag,
46			(QT_BLOCK(freet)+1)QT_BLOCK_SIZE/(8sizeof(int4)));
47			if (quad_flag == NULL)
48			return(EMPTY);
49	gwlarson	3.1	}
50			treetop += 4;
51			return(freet);
52			}
53
54
55	gwlarson	3.3	qtClearAllFlags() /* clear all quadtree branch flags */
56			{
57			if (!treetop) return;
58			bzero((char *)quad_flag,
59			(QT_BLOCK(treetop-1)+1)QT_BLOCK_SIZE/(8sizeof(int4)));
60			}
61
62
63	gwlarson	3.1	qtFree(qt) /* free a quadtree */
64			register QUADTREE qt;
65			{
66			register int i;
67
68			if (!QT_IS_TREE(qt))
69			{
70			qtfreeleaf(qt);
71			return;
72			}
73			for (i = 0; i < 4; i++)
74			qtFree(QT_NTH_CHILD(qt, i));
75			QT_NTH_CHILD(qt, 0) = quad_free_list;
76			quad_free_list = qt;
77			}
78
79
80			qtDone() /* free EVERYTHING */
81			{
82			register int i;
83
84			for (i = 0; i < QT_MAX_BLK; i++)
85			{
86	gwlarson	3.3	if (quad_block[i] == NULL)
87			break;
88			free((char *)quad_block[i]);
89	gwlarson	3.1	quad_block[i] = NULL;
90			}
91	gwlarson	3.3	if (i) free((char *)quad_flag);
92			quad_flag = NULL;
93	gwlarson	3.1	quad_free_list = EMPTY;
94			treetop = 0;
95			}
96
97	gwlarson	3.3
98	gwlarson	3.1	QUADTREE
99			qtCompress(qt) /* recursively combine nodes */
100			register QUADTREE qt;
101			{
102			register int i;
103			register QUADTREE qres;
104
105			if (!QT_IS_TREE(qt)) /* not a tree */
106			return(qt);
107			qres = QT_NTH_CHILD(qt,0) = qtCompress(QT_NTH_CHILD(qt,0));
108			for (i = 1; i < 4; i++)
109			if((QT_NTH_CHILD(qt,i) = qtCompress(QT_NTH_CHILD(qt,i))) != qres)
110			qres = qt;
111			if(!QT_IS_TREE(qres))
112			{ /* all were identical leaves */
113			QT_NTH_CHILD(qt,0) = quad_free_list;
114			quad_free_list = qt;
115			}
116			return(qres);
117			}
118
119	gwlarson	3.3
120	gwlarson	3.1	QUADTREE
121	gwlarson	3.3	*qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
122	gwlarson	3.1	QUADTREE *qtptr;
123	gwlarson	3.2	double bcoord[3];
124	gwlarson	3.3	FVECT t0,t1,t2;
125	gwlarson	3.2	{
126			int i;
127			QUADTREE *child;
128	gwlarson	3.3	FVECT a,b,c;
129	gwlarson	3.2
130			if(QT_IS_TREE(*qtptr))
131			{
132			i = bary2d_child(bcoord);
133			child = QT_NTH_CHILD_PTR(*qtptr,i);
134	gwlarson	3.3	if(t0)
135			{
136			qtSubdivide_tri(t0,t1,t2,a,b,c);
137			qtNth_child_tri(t0,t1,t2,a,b,c,i,t0,t1,t2);
138			}
139			return(qtLocate_leaf(child,bcoord,t0,t1,t2));
140	gwlarson	3.2	}
141			else
142	gwlarson	3.3	return(qtptr);
143	gwlarson	3.2	}
144
145
146
147	gwlarson	3.3	int
148			qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n)
149			QUADTREE *qtptr;
150			double bcoord[3];
151			int id;
152			FVECT v0,v1,v2;
153			int n;
154			{
155			int i;
156			QUADTREE *child;
157			OBJECT os[MAXSET+1],*optr;
158			int found;
159			FVECT r0,r1,r2;
160
161			if(QT_IS_TREE(*qtptr))
162			{
163			QT_SET_FLAG(*qtptr);
164			i = bary2d_child(bcoord);
165			child = QT_NTH_CHILD_PTR(*qtptr,i);
166			return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,++n));
167			}
168			else
169			{
170			/* If this leave node emptry- create a new set */
171			if(QT_IS_EMPTY(*qtptr))
172			qtptr = qtaddelem(qtptr,id);
173			else
174			{
175			qtgetset(os,*qtptr);
176			/* If the set is too large: subdivide */
177			if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
178			qtptr = qtaddelem(qtptr,id);
179			else
180			{
181			if (n < QT_MAX_LEVELS)
182			{
183			/* If set size exceeds threshold: subdivide cell and
184			reinsert set tris into cell
185			*/
186			n++;
187			qtfreeleaf(*qtptr);
188			qtSubdivide(qtptr);
189			QT_SET_FLAG(*qtptr);
190			found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n);
191			#if 0
192			if(!found)
193			{
194			#ifdef TEST_DRIVER
195			HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
196			#else
197			eputs("qtAdd_tri():Found in parent but not children\n");
198			#endif
199			}
200			#endif
201			for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
202			{
203			id = QT_SET_NEXT_ELEM(optr);
204			qtTri_verts_from_id(id,r0,r1,r2);
205			found= qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,++n);
206			#ifdef DEBUG
207			if(!found)
208			eputs("qtAdd_tri():Reinsert-in parent but not children\n");
209			#endif
210			}
211			}
212			else
213			if(QT_SET_CNT(os) < QT_MAX_SET)
214			{
215			qtptr = qtaddelem(qtptr,id);
216			}
217			else
218			{
219			#ifdef DEBUG
220			eputs("qtAdd_tri():two many levels\n");
221			#endif
222			return(FALSE);
223			}
224			}
225			}
226			}
227			return(TRUE);
228			}
229	gwlarson	3.2
230	gwlarson	3.3
231			int
232			qtAdd_tri_from_point(qtptr,v0,v1,v2,pt,id)
233			QUADTREE *qtptr;
234			FVECT v0,v1,v2;
235			FVECT pt;
236			int id;
237			{
238			char d,w;
239			int i,x,y;
240			QUADTREE *child;
241			QUADTREE qt;
242			FVECT i_pt,n,a,b,c,r0,r1,r2;
243			double pd,bcoord[3];
244			OBJECT os[MAXSET+1],*optr;
245			int found;
246
247			/* Determine if point lies within pyramid (and therefore
248			inside a spherical quadtree cell):GT_INTERIOR, on one of the
249			pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
250			or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
251			For each triangle edge: compare the
252			point against the plane formed by the edge and the view center
253			*/
254			d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);
255
256			/* Not in this triangle */
257			if(!d)
258			return(FALSE);
259
260			/* 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(*qtptr))
266			{
267			QT_SET_FLAG(*qtptr);
268			/* Find the intersection point */
269			tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
270			intersect_vector_plane(pt,n,pd,NULL,i_pt);
271
272			i = max_index(n);
273			x = (i+1)%3;
274			y = (i+2)%3;
275			/* Calculate barycentric coordinates of i_pt */
276			bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
277
278			i = bary2d_child(bcoord);
279			child = QT_NTH_CHILD_PTR(*qtptr,i);
280			/* NOTE: Optimize: only subdivide for specified child */
281			qtSubdivide_tri(v0,v1,v2,a,b,c);
282			qtNth_child_tri(v0,v1,v2,a,b,c,i,v0,v1,v2);
283			return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,1));
284			}
285			else
286			{
287			/* If this leave node emptry- create a new set */
288			if(QT_IS_EMPTY(*qtptr))
289			qtptr = qtaddelem(qtptr,id);
290			else
291			{
292			qtgetset(os,*qtptr);
293			/* If the set is too large: subdivide */
294			if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
295			qtptr = qtaddelem(qtptr,id);
296			else
297			{
298			/* If set size exceeds threshold: subdivide cell and
299			reinsert set tris into cell
300			*/
301			qtfreeleaf(*qtptr);
302			qtSubdivide(qtptr);
303			found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,1);
304			#if 0
305			if(!found)
306			{
307			#ifdef TEST_DRIVER
308			HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
309			#else
310			eputs("qtAdd_tri():Found in parent but not children\n");
311			#endif
312			}
313			#endif
314			for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
315			{
316			id = QT_SET_NEXT_ELEM(optr);
317			qtTri_verts_from_id(id,r0,r1,r2);
318			found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,1);
319			#ifdef DEBUG
320			if(!found)
321			eputs("qtAdd_tri():Reinsert-in parent but not children\n");
322			#endif
323			}
324			}
325			}
326			}
327			return(TRUE);
328			}
329
330
331	gwlarson	3.2	QUADTREE
332	gwlarson	3.3	*qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
333	gwlarson	3.2	QUADTREE *qtptr;
334			FVECT v0,v1,v2;
335	gwlarson	3.1	FVECT pt;
336	gwlarson	3.3	FVECT t0,t1,t2;
337	gwlarson	3.1	{
338			char d,w;
339	gwlarson	3.2	int i,x,y;
340	gwlarson	3.1	QUADTREE *child;
341			QUADTREE qt;
342	gwlarson	3.3	FVECT n,i_pt,a,b,c;
343	gwlarson	3.2	double pd,bcoord[3];
344	gwlarson	3.1
345			/* Determine if point lies within pyramid (and therefore
346			inside a spherical quadtree cell):GT_INTERIOR, on one of the
347			pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
348			or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
349			For each triangle edge: compare the
350			point against the plane formed by the edge and the view center
351			*/
352	gwlarson	3.2	d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);
353	gwlarson	3.1
354			/* Not in this triangle */
355			if(!d)
356			{
357	gwlarson	3.3	return(NULL);
358	gwlarson	3.1	}
359
360			/* Will return lowest level triangle containing point: It the
361			point is on an edge or vertex: will return first associated
362			triangle encountered in the child traversal- the others can
363			be derived using triangle adjacency information
364			*/
365			if(QT_IS_TREE(*qtptr))
366			{
367	gwlarson	3.2	/* Find the intersection point */
368			tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
369			intersect_vector_plane(pt,n,pd,NULL,i_pt);
370
371			i = max_index(n);
372			x = (i+1)%3;
373			y = (i+2)%3;
374			/* Calculate barycentric coordinates of i_pt */
375			bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
376
377			i = bary2d_child(bcoord);
378			child = QT_NTH_CHILD_PTR(*qtptr,i);
379	gwlarson	3.3	if(t0)
380			{
381			qtSubdivide_tri(v0,v1,v2,a,b,c);
382			qtNth_child_tri(v0,v1,v2,a,b,c,i,t0,t1,t2);
383			}
384			return(qtLocate_leaf(child,bcoord,t0,t1,t2));
385	gwlarson	3.1	}
386			else
387	gwlarson	3.3	return(qtptr);
388	gwlarson	3.1	}
389
390			int
391			qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
392			QUADTREE *qtptr;
393			FVECT v0,v1,v2;
394			FVECT pt;
395			char type,which;
396			{
397			OBJECT os[MAXSET+1],*optr;
398			char d,w;
399			int i,id;
400			FVECT p0,p1,p2;
401	gwlarson	3.2
402	gwlarson	3.3	qtptr = qtRoot_point_locate(qtptr,v0,v1,v2,pt,NULL,NULL,NULL);
403	gwlarson	3.2
404	gwlarson	3.3	if(!qtptr \|\| QT_IS_EMPTY(*qtptr))
405	gwlarson	3.1	return(EMPTY);
406
407			/* Get the set */
408	gwlarson	3.3	qtgetset(os,*qtptr);
409	gwlarson	3.1	for (i = QT_SET_CNT(os),optr = QT_SET_PTR(os); i > 0; i--)
410			{
411			/* Find the triangle that pt falls in (NOTE:FOR now return first 1) */
412			id = QT_SET_NEXT_ELEM(optr);
413			qtTri_verts_from_id(id,p0,p1,p2);
414			d = test_single_point_against_spherical_tri(p0,p1,p2,pt,&w);
415			if(d)
416			{
417			if(type)
418			*type = d;
419			if(which)
420			*which = w;
421			return(id);
422			}
423			}
424			return(EMPTY);
425			}
426
427			QUADTREE
428			qtSubdivide(qtptr)
429			QUADTREE *qtptr;
430			{
431			QUADTREE node;
432			node = qtAlloc();
433			QT_CLEAR_CHILDREN(node);
434			*qtptr = node;
435			return(node);
436			}
437
438
439			QUADTREE
440			qtSubdivide_nth_child(qt,n)
441			QUADTREE qt;
442			int n;
443			{
444			QUADTREE node;
445
446			node = qtSubdivide(&(QT_NTH_CHILD(qt,n)));
447
448			return(node);
449			}
450
451	gwlarson	3.2	/* for triangle v0-v1-v2- returns a,b,c: children are:
452	gwlarson	3.1
453	gwlarson	3.2	v2 0: v0,a,c
454			/\ 1: a,v1,b
455			/2 \ 2: c,b,v2
456			c/____\b 3: b,c,a
457	gwlarson	3.1	/\ /\
458	gwlarson	3.2	/0 \3 /1 \
459			v0____\/____\v1
460	gwlarson	3.1	a
461			*/
462
463	gwlarson	3.2	qtSubdivide_tri(v0,v1,v2,a,b,c)
464			FVECT v0,v1,v2;
465	gwlarson	3.1	FVECT a,b,c;
466			{
467	gwlarson	3.2	EDGE_MIDPOINT_VEC3(a,v0,v1);
468			EDGE_MIDPOINT_VEC3(b,v1,v2);
469			EDGE_MIDPOINT_VEC3(c,v2,v0);
470	gwlarson	3.1	}
471
472	gwlarson	3.2	qtNth_child_tri(v0,v1,v2,a,b,c,i,r0,r1,r2)
473			FVECT v0,v1,v2;
474	gwlarson	3.1	FVECT a,b,c;
475			int i;
476	gwlarson	3.2	FVECT r0,r1,r2;
477	gwlarson	3.1	{
478			switch(i){
479			case 0:
480	gwlarson	3.2	VCOPY(r0,v0); VCOPY(r1,a); VCOPY(r2,c);
481	gwlarson	3.1	break;
482			case 1:
483	gwlarson	3.2	VCOPY(r0,a); VCOPY(r1,v1); VCOPY(r2,b);
484	gwlarson	3.1	break;
485			case 2:
486	gwlarson	3.2	VCOPY(r0,c); VCOPY(r1,b); VCOPY(r2,v2);
487	gwlarson	3.1	break;
488			case 3:
489	gwlarson	3.2	VCOPY(r0,b); VCOPY(r1,c); VCOPY(r2,a);
490	gwlarson	3.1	break;
491			}
492			}
493
494			/* Add triangle "id" to all leaf level cells that are children of
495			quadtree pointed to by "qtptr" with cell vertices "t1,t2,t3"
496			that it overlaps (vertex and edge adjacencies do not count
497			as overlapping). If the addition of the triangle causes the cell to go over
498			threshold- the cell is split- and the triangle must be recursively inserted
499			into the new child cells: it is assumed that "v1,v2,v3" are normalized
500			*/
501
502			int
503	gwlarson	3.3	qtRoot_add_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
504			QUADTREE *qtptr;
505			int id;
506			FVECT t0,t1,t2;
507			FVECT v0,v1,v2;
508			int n;
509			{
510			char test;
511			int found;
512
513			test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
514			if(!test)
515			return(FALSE);
516
517			found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
518			return(found);
519			}
520
521			int
522			qtRoot_add_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n)
523			QUADTREE *qtptr;
524			int id;
525			FVECT t0,t1,t2;
526			FVECT v0,v1,v2;
527			int n;
528			{
529			char test;
530			int found;
531
532			test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
533			if(!test)
534			return(FALSE);
535
536			found = qtAdd_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n);
537			return(found);
538			}
539
540			int
541	gwlarson	3.2	qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
542	gwlarson	3.1	QUADTREE *qtptr;
543			int id;
544	gwlarson	3.2	FVECT t0,t1,t2;
545			FVECT v0,v1,v2;
546	gwlarson	3.1	int n;
547			{
548
549			char test;
550			int i,index;
551			FVECT a,b,c;
552			OBJECT os[MAXSET+1],*optr;
553			QUADTREE qt;
554			int found;
555	gwlarson	3.2	FVECT r0,r1,r2;
556	gwlarson	3.1
557			found = 0;
558			/* if this is tree: recurse */
559			if(QT_IS_TREE(*qtptr))
560			{
561			n++;
562	gwlarson	3.2	qtSubdivide_tri(v0,v1,v2,a,b,c);
563	gwlarson	3.3	test = spherical_tri_intersect(t0,t1,t2,v0,a,c);
564			if(test)
565			found \|= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c,n);
566			test = spherical_tri_intersect(t0,t1,t2,a,v1,b);
567			if(test)
568			found \|= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b,n);
569			test = spherical_tri_intersect(t0,t1,t2,c,b,v2);
570			if(test)
571			found \|= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2,n);
572			test = spherical_tri_intersect(t0,t1,t2,b,c,a);
573			if(test)
574			found \|= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a,n);
575	gwlarson	3.1
576			#if 0
577			if(!found)
578			{
579			#ifdef TEST_DRIVER
580			HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
581			#else
582			eputs("qtAdd_tri():Found in parent but not children\n");
583			#endif
584			}
585			#endif
586			}
587			else
588			{
589			/* If this leave node emptry- create a new set */
590			if(QT_IS_EMPTY(*qtptr))
591	gwlarson	3.3	qtptr = qtaddelem(qtptr,id);
592	gwlarson	3.1	else
593			{
594			qtgetset(os,*qtptr);
595			/* If the set is too large: subdivide */
596			if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
597			qtptr = qtaddelem(qtptr,id);
598			else
599			{
600			if (n < QT_MAX_LEVELS)
601			{
602			/* If set size exceeds threshold: subdivide cell and
603			reinsert set tris into cell
604			*/
605			n++;
606			qtfreeleaf(*qtptr);
607			qtSubdivide(qtptr);
608	gwlarson	3.2	found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
609	gwlarson	3.1	#if 0
610			if(!found)
611			{
612			#ifdef TEST_DRIVER
613			HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
614			#else
615			eputs("qtAdd_tri():Found in parent but not children\n");
616			#endif
617			}
618			#endif
619			for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
620			{
621			id = QT_SET_NEXT_ELEM(optr);
622	gwlarson	3.2	qtTri_verts_from_id(id,r0,r1,r2);
623			found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,n);
624	gwlarson	3.1	#ifdef DEBUG
625			if(!found)
626			eputs("qtAdd_tri():Reinsert-in parent but not children\n");
627			#endif
628			}
629			}
630			else
631			if(QT_SET_CNT(os) < QT_MAX_SET)
632			{
633			qtptr = qtaddelem(qtptr,id);
634			#if 0
635			{
636			int k;
637			qtgetset(os,*qtptr);
638			printf("\n%d:\n",os[0]);
639			for(k=1; k <= os[0];k++)
640			printf("%d ",os[k]);
641			printf("\n");
642			}
643			#endif
644			/*
645			insertelem(os,id);
646			*qtptr = fullnode(os);
647			*/
648			}
649			else
650			{
651			#ifdef DEBUG
652			eputs("qtAdd_tri():two many levels\n");
653			#endif
654			return(FALSE);
655			}
656			}
657			}
658			}
659			return(TRUE);
660			}
661
662
663			int
664	gwlarson	3.2	qtApply_to_tri_cells(qtptr,t0,t1,t2,v0,v1,v2,func,arg)
665	gwlarson	3.1	QUADTREE *qtptr;
666	gwlarson	3.2	FVECT t0,t1,t2;
667			FVECT v0,v1,v2;
668	gwlarson	3.1	int (*func)();
669			char *arg;
670			{
671			char test;
672			FVECT a,b,c;
673
674	gwlarson	3.2	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
675			test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
676	gwlarson	3.1
677			/* If triangles do not overlap: done */
678			if(!test)
679			return(FALSE);
680
681			/* if this is tree: recurse */
682			if(QT_IS_TREE(*qtptr))
683			{
684	gwlarson	3.2	qtSubdivide_tri(v0,v1,v2,a,b,c);
685			qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t0,t1,t2,v0,a,c,func,arg);
686			qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t0,t1,t2,a,v1,b,func,arg);
687			qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t0,t1,t2,c,b,v2,func,arg);
688			qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t0,t1,t2,b,c,a,func,arg);
689	gwlarson	3.1	}
690			else
691			return(func(qtptr,arg));
692			}
693
694
695			int
696	gwlarson	3.2	qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
697	gwlarson	3.1	QUADTREE *qtptr;
698			int id;
699	gwlarson	3.2	FVECT t0,t1,t2;
700			FVECT v0,v1,v2;
701	gwlarson	3.1	{
702
703			char test;
704			int i;
705			FVECT a,b,c;
706			OBJECT os[MAXSET+1];
707
708	gwlarson	3.2	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
709			test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
710	gwlarson	3.1
711			/* If triangles do not overlap: done */
712			if(!test)
713			return(FALSE);
714
715			/* if this is tree: recurse */
716			if(QT_IS_TREE(*qtptr))
717			{
718	gwlarson	3.2	qtSubdivide_tri(v0,v1,v2,a,b,c);
719			qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c);
720			qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b);
721			qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2);
722			qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a);
723	gwlarson	3.1	}
724			else
725			{
726			if(QT_IS_EMPTY(*qtptr))
727			{
728			#ifdef DEBUG
729			eputs("qtRemove_tri(): triangle not found\n");
730			#endif
731			}
732			/* remove id from set */
733			else
734			{
735			qtgetset(os,*qtptr);
736			if(!inset(os,id))
737			{
738			#ifdef DEBUG
739			eputs("qtRemove_tri(): tri not in set\n");
740			#endif
741			}
742			else
743			{
744			qtptr = qtdelelem(qtptr,id);
745			}
746			}
747			}
748			return(TRUE);
749			}
750	gwlarson	3.3
751	gwlarson	3.2
752
753
754
755
756
757
758
759
760
761