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 */


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(
                   (unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
           return(EMPTY);
    }
    treetop += 4;
    return(freet);
}


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++)
        {
            free((char *)quad_block[i],
                  (unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE));
            quad_block[i] = 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
qtPoint_locate(qtptr,v1,v2,v3,pt,type,which,p0,p1,p2)
QUADTREE *qtptr;
FVECT v1,v2,v3;
FVECT pt;
char *type,*which;
FVECT p0,p1,p2;
{
    char d,w;
    int i;
    QUADTREE *child;
    QUADTREE qt;
    FVECT a,b,c;
    FVECT t1,t2,t3;

    /* 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(v1,v2,v3,pt,&w);  

    /* Not in this triangle */
    if(!d)
    {
      if(which)
        *which = 0;
      return(EMPTY);
    }

    /* 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))
    {  
      qtSubdivide_tri(v1,v2,v3,a,b,c);
      child = QT_NTH_CHILD_PTR(*qtptr,0);
      if(!QT_IS_EMPTY(*child))
      {
        qt = qtPoint_locate(child,v1,a,c,pt,type,which,p0,p1,p2);
        if(!QT_IS_EMPTY(qt))
          return(qt);
      }
      child = QT_NTH_CHILD_PTR(*qtptr,1);
      if(!QT_IS_EMPTY(*child))
      {
        qt = qtPoint_locate(child,a,b,c,pt,type,which,p0,p1,p2);
        if(!QT_IS_EMPTY(qt))
          return(qt);
      }
      child = QT_NTH_CHILD_PTR(*qtptr,2);
      if(!QT_IS_EMPTY(*child))
      {
        qt = qtPoint_locate(child,a,v2,b,pt,type,which,p0,p1,p2);
        if(!QT_IS_EMPTY(qt))
          return(qt);
      }
      child = QT_NTH_CHILD_PTR(*qtptr,3);
      if(!QT_IS_EMPTY(*child))
      {
        qt = qtPoint_locate(child,c,b,v3,pt,type,which,p0,p1,p2);
        if(!QT_IS_EMPTY(qt))
          return(qt);
      }
    }
    else
       if(!QT_IS_EMPTY(*qtptr))
       {  
           /* map GT_VERTEX,GT_EDGE,GT_FACE GT_INTERIOR of pyramid to
              spherical triangle primitives
              */
         if(type)
           *type = d;
         if(which)
           *which = w;
         VCOPY(p0,v1);
         VCOPY(p1,v2);
         VCOPY(p2,v3);
         return(*qtptr);
       }
    return(EMPTY);
}

int
qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
QUADTREE *qtptr;
FVECT v0,v1,v2;
FVECT pt;
char *type,*which;
{
    QUADTREE qt;
    OBJECT os[MAXSET+1],*optr;
    char d,w;
    int i,id;
    FVECT p0,p1,p2;
  
    qt = qtPoint_locate(qtptr,v0,v1,v2,pt,type,which,p0,p1,p2);
    if(QT_IS_EMPTY(qt))
       return(EMPTY);
    
    /* Get the set */
    qtgetset(os,qt);
    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 v1-v2-v3- returns a,b,c: children are:

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

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

qtNth_child_tri(v1,v2,v3,a,b,c,i,r1,r2,r3)
FVECT v1,v2,v3;
FVECT a,b,c;
int i;
FVECT r1,r2,r3;
{
  VCOPY(r1,a); VCOPY(r2,b);    VCOPY(r3,c);
  switch(i){
  case 0:  
    VCOPY(r2,r1);
    VCOPY(r1,v1);
    break;
  case 1:  
    break;
  case 2:  
    VCOPY(r3,r2);
    VCOPY(r2,v2);
    break;
  case 3:  
    VCOPY(r1,r3);
    VCOPY(r3,v3);
    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
qtAdd_tri(qtptr,id,t1,t2,t3,v1,v2,v3,n)
QUADTREE *qtptr;
int id;
FVECT t1,t2,t3;
FVECT v1,v2,v3;
int n;
{
  
  char test;
  int i,index;
  FVECT a,b,c;
  OBJECT os[MAXSET+1],*optr;
  QUADTREE qt;
  int found;
  FVECT r1,r2,r3;

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

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

  /* if this is tree: recurse */
  if(QT_IS_TREE(*qtptr))
  {
    n++;
    qtSubdivide_tri(v1,v2,v3,a,b,c);
    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t1,t2,t3,v1,a,c,n);
    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t1,t2,t3,a,b,c,n);
    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t1,t2,t3,a,v2,b,n);
    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t1,t2,t3,c,b,v3,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);
#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
          /*
          os[0] = 0;
          insertelem(os,id);
          qt = fullnode(os);
          *qtptr = qt;
          */
      }
      else
      {
          qtgetset(os,*qtptr);
          /* If the set is too large: subdivide */
          if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
          {
            *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
          {
            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,t1,t2,t3,v1,v2,v3,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,r1,r2,r3);
                   found=qtAdd_tri(qtptr,id,r1,r2,r3,v1,v2,v3,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,t1,t2,t3,v1,v2,v3,func,arg)
QUADTREE *qtptr;
FVECT t1,t2,t3;
FVECT v1,v2,v3;
int (*func)();
char *arg;
{
  char test;
  FVECT a,b,c;

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

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

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


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

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

  /* if this is tree: recurse */
  if(QT_IS_TREE(*qtptr))
  {
    qtSubdivide_tri(v1,v2,v3,a,b,c);
    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t1,t2,t3,v1,a,c);
    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t1,t2,t3,a,b,c);
    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t1,t2,t3,a,v2,b);
    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t1,t2,t3,c,b,v3);
  }
  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);
#if 0
            {
              int k;
              if(!QT_IS_EMPTY(*qtptr))
              {qtgetset(os,*qtptr);
              printf("\n%d:\n",os[0]);
              for(k=1; k <= os[0];k++)
                 printf("%d ",os[k]);
              printf("\n");
           }

          }
#endif
        }
      }
  }
  return(TRUE);
}
Revision:	3.1
Committed:	Wed Aug 19 17:45:24 1998 UTC (25 years, 8 months ago) by gwlarson
Content type:	text/plain
Branch:	MAIN
Log Message:	Initial revision
#	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
26
27
28			QUADTREE
29			qtAlloc() /* allocate a quadtree */
30			{
31			register QUADTREE freet;
32
33			if ((freet = quad_free_list) != EMPTY)
34			{
35			quad_free_list = QT_NTH_CHILD(freet, 0);
36			return(freet);
37			}
38			freet = treetop;
39			if (QT_BLOCK_INDEX(freet) == 0)
40			{
41			if (QT_BLOCK(freet) >= QT_MAX_BLK)
42			return(EMPTY);
43			if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
44			(unsigned)QT_BLOCK_SIZE4sizeof(QUADTREE))) == NULL)
45			return(EMPTY);
46			}
47			treetop += 4;
48			return(freet);
49			}
50
51
52			qtFree(qt) /* free a quadtree */
53			register QUADTREE qt;
54			{
55			register int i;
56
57			if (!QT_IS_TREE(qt))
58			{
59			qtfreeleaf(qt);
60			return;
61			}
62			for (i = 0; i < 4; i++)
63			qtFree(QT_NTH_CHILD(qt, i));
64			QT_NTH_CHILD(qt, 0) = quad_free_list;
65			quad_free_list = qt;
66			}
67
68
69			qtDone() /* free EVERYTHING */
70			{
71			register int i;
72
73			for (i = 0; i < QT_MAX_BLK; i++)
74			{
75			free((char *)quad_block[i],
76			(unsigned)QT_BLOCK_SIZE4sizeof(QUADTREE));
77			quad_block[i] = NULL;
78			}
79			quad_free_list = EMPTY;
80			treetop = 0;
81			}
82
83			QUADTREE
84			qtCompress(qt) /* recursively combine nodes */
85			register QUADTREE qt;
86			{
87			register int i;
88			register QUADTREE qres;
89
90			if (!QT_IS_TREE(qt)) /* not a tree */
91			return(qt);
92			qres = QT_NTH_CHILD(qt,0) = qtCompress(QT_NTH_CHILD(qt,0));
93			for (i = 1; i < 4; i++)
94			if((QT_NTH_CHILD(qt,i) = qtCompress(QT_NTH_CHILD(qt,i))) != qres)
95			qres = qt;
96			if(!QT_IS_TREE(qres))
97			{ /* all were identical leaves */
98			QT_NTH_CHILD(qt,0) = quad_free_list;
99			quad_free_list = qt;
100			}
101			return(qres);
102			}
103
104			QUADTREE
105			qtPoint_locate(qtptr,v1,v2,v3,pt,type,which,p0,p1,p2)
106			QUADTREE *qtptr;
107			FVECT v1,v2,v3;
108			FVECT pt;
109			char type,which;
110			FVECT p0,p1,p2;
111			{
112			char d,w;
113			int i;
114			QUADTREE *child;
115			QUADTREE qt;
116			FVECT a,b,c;
117			FVECT t1,t2,t3;
118
119			/* Determine if point lies within pyramid (and therefore
120			inside a spherical quadtree cell):GT_INTERIOR, on one of the
121			pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
122			or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
123			For each triangle edge: compare the
124			point against the plane formed by the edge and the view center
125			*/
126			d = test_single_point_against_spherical_tri(v1,v2,v3,pt,&w);
127
128			/* Not in this triangle */
129			if(!d)
130			{
131			if(which)
132			*which = 0;
133			return(EMPTY);
134			}
135
136			/* Will return lowest level triangle containing point: It the
137			point is on an edge or vertex: will return first associated
138			triangle encountered in the child traversal- the others can
139			be derived using triangle adjacency information
140			*/
141
142			if(QT_IS_TREE(*qtptr))
143			{
144			qtSubdivide_tri(v1,v2,v3,a,b,c);
145			child = QT_NTH_CHILD_PTR(*qtptr,0);
146			if(!QT_IS_EMPTY(*child))
147			{
148			qt = qtPoint_locate(child,v1,a,c,pt,type,which,p0,p1,p2);
149			if(!QT_IS_EMPTY(qt))
150			return(qt);
151			}
152			child = QT_NTH_CHILD_PTR(*qtptr,1);
153			if(!QT_IS_EMPTY(*child))
154			{
155			qt = qtPoint_locate(child,a,b,c,pt,type,which,p0,p1,p2);
156			if(!QT_IS_EMPTY(qt))
157			return(qt);
158			}
159			child = QT_NTH_CHILD_PTR(*qtptr,2);
160			if(!QT_IS_EMPTY(*child))
161			{
162			qt = qtPoint_locate(child,a,v2,b,pt,type,which,p0,p1,p2);
163			if(!QT_IS_EMPTY(qt))
164			return(qt);
165			}
166			child = QT_NTH_CHILD_PTR(*qtptr,3);
167			if(!QT_IS_EMPTY(*child))
168			{
169			qt = qtPoint_locate(child,c,b,v3,pt,type,which,p0,p1,p2);
170			if(!QT_IS_EMPTY(qt))
171			return(qt);
172			}
173			}
174			else
175			if(!QT_IS_EMPTY(*qtptr))
176			{
177			/* map GT_VERTEX,GT_EDGE,GT_FACE GT_INTERIOR of pyramid to
178			spherical triangle primitives
179			*/
180			if(type)
181			*type = d;
182			if(which)
183			*which = w;
184			VCOPY(p0,v1);
185			VCOPY(p1,v2);
186			VCOPY(p2,v3);
187			return(*qtptr);
188			}
189			return(EMPTY);
190			}
191
192			int
193			qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
194			QUADTREE *qtptr;
195			FVECT v0,v1,v2;
196			FVECT pt;
197			char type,which;
198			{
199			QUADTREE qt;
200			OBJECT os[MAXSET+1],*optr;
201			char d,w;
202			int i,id;
203			FVECT p0,p1,p2;
204
205			qt = qtPoint_locate(qtptr,v0,v1,v2,pt,type,which,p0,p1,p2);
206			if(QT_IS_EMPTY(qt))
207			return(EMPTY);
208
209			/* Get the set */
210			qtgetset(os,qt);
211			for (i = QT_SET_CNT(os),optr = QT_SET_PTR(os); i > 0; i--)
212			{
213			/* Find the triangle that pt falls in (NOTE:FOR now return first 1) */
214			id = QT_SET_NEXT_ELEM(optr);
215			qtTri_verts_from_id(id,p0,p1,p2);
216			d = test_single_point_against_spherical_tri(p0,p1,p2,pt,&w);
217			if(d)
218			{
219			if(type)
220			*type = d;
221			if(which)
222			*which = w;
223			return(id);
224			}
225			}
226			return(EMPTY);
227			}
228
229			QUADTREE
230			qtSubdivide(qtptr)
231			QUADTREE *qtptr;
232			{
233			QUADTREE node;
234			node = qtAlloc();
235			QT_CLEAR_CHILDREN(node);
236			*qtptr = node;
237			return(node);
238			}
239
240
241			QUADTREE
242			qtSubdivide_nth_child(qt,n)
243			QUADTREE qt;
244			int n;
245			{
246			QUADTREE node;
247
248			node = qtSubdivide(&(QT_NTH_CHILD(qt,n)));
249
250			return(node);
251			}
252
253			/* for triangle v1-v2-v3- returns a,b,c: children are:
254
255			v3 0: v1,a,c
256			/\ 1: a,b,c
257			/3 \ 2: a,v2,b
258			c/____\b 3: c,b,v3
259			/\ /\
260			/0 \1 /2 \
261			v1/____\/____\v2
262			a
263			*/
264
265			qtSubdivide_tri(v1,v2,v3,a,b,c)
266			FVECT v1,v2,v3;
267			FVECT a,b,c;
268			{
269			EDGE_MIDPOINT_VEC3(a,v1,v2);
270			normalize(a);
271			EDGE_MIDPOINT_VEC3(b,v2,v3);
272			normalize(b);
273			EDGE_MIDPOINT_VEC3(c,v3,v1);
274			normalize(c);
275			}
276
277			qtNth_child_tri(v1,v2,v3,a,b,c,i,r1,r2,r3)
278			FVECT v1,v2,v3;
279			FVECT a,b,c;
280			int i;
281			FVECT r1,r2,r3;
282			{
283			VCOPY(r1,a); VCOPY(r2,b); VCOPY(r3,c);
284			switch(i){
285			case 0:
286			VCOPY(r2,r1);
287			VCOPY(r1,v1);
288			break;
289			case 1:
290			break;
291			case 2:
292			VCOPY(r3,r2);
293			VCOPY(r2,v2);
294			break;
295			case 3:
296			VCOPY(r1,r3);
297			VCOPY(r3,v3);
298			break;
299			}
300			}
301
302			/* Add triangle "id" to all leaf level cells that are children of
303			quadtree pointed to by "qtptr" with cell vertices "t1,t2,t3"
304			that it overlaps (vertex and edge adjacencies do not count
305			as overlapping). If the addition of the triangle causes the cell to go over
306			threshold- the cell is split- and the triangle must be recursively inserted
307			into the new child cells: it is assumed that "v1,v2,v3" are normalized
308			*/
309
310			int
311			qtAdd_tri(qtptr,id,t1,t2,t3,v1,v2,v3,n)
312			QUADTREE *qtptr;
313			int id;
314			FVECT t1,t2,t3;
315			FVECT v1,v2,v3;
316			int n;
317			{
318
319			char test;
320			int i,index;
321			FVECT a,b,c;
322			OBJECT os[MAXSET+1],*optr;
323			QUADTREE qt;
324			int found;
325			FVECT r1,r2,r3;
326
327			/* test if triangle (t1,t2,t3) overlaps cell triangle (v1,v2,v3) */
328			test = spherical_tri_intersect(t1,t2,t3,v1,v2,v3);
329
330			/* If triangles do not overlap: done */
331			if(!test)
332			return(FALSE);
333			found = 0;
334
335			/* if this is tree: recurse */
336			if(QT_IS_TREE(*qtptr))
337			{
338			n++;
339			qtSubdivide_tri(v1,v2,v3,a,b,c);
340			found \|= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t1,t2,t3,v1,a,c,n);
341			found \|= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t1,t2,t3,a,b,c,n);
342			found \|= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t1,t2,t3,a,v2,b,n);
343			found \|= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t1,t2,t3,c,b,v3,n);
344
345			#if 0
346			if(!found)
347			{
348			#ifdef TEST_DRIVER
349			HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
350			#else
351			eputs("qtAdd_tri():Found in parent but not children\n");
352			#endif
353			}
354			#endif
355			}
356			else
357			{
358			/* If this leave node emptry- create a new set */
359			if(QT_IS_EMPTY(*qtptr))
360			{
361			qtptr = qtaddelem(qtptr,id);
362			#if 0
363			{
364			int k;
365			qtgetset(os,*qtptr);
366			printf("\n%d:\n",os[0]);
367			for(k=1; k <= os[0];k++)
368			printf("%d ",os[k]);
369			printf("\n");
370			}
371			#endif
372			/*
373			os[0] = 0;
374			insertelem(os,id);
375			qt = fullnode(os);
376			*qtptr = qt;
377			*/
378			}
379			else
380			{
381			qtgetset(os,*qtptr);
382			/* If the set is too large: subdivide */
383			if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
384			{
385			qtptr = qtaddelem(qtptr,id);
386			#if 0
387			{
388			int k;
389			qtgetset(os,*qtptr);
390			printf("\n%d:\n",os[0]);
391			for(k=1; k <= os[0];k++)
392			printf("%d ",os[k]);
393			printf("\n");
394			}
395			#endif
396			/*
397			insertelem(os,id);
398			*qtptr = fullnode(os);
399			*/
400			}
401			else
402			{
403			if (n < QT_MAX_LEVELS)
404			{
405			/* If set size exceeds threshold: subdivide cell and
406			reinsert set tris into cell
407			*/
408			n++;
409			qtfreeleaf(*qtptr);
410			qtSubdivide(qtptr);
411			found = qtAdd_tri(qtptr,id,t1,t2,t3,v1,v2,v3,n);
412			#if 0
413			if(!found)
414			{
415			#ifdef TEST_DRIVER
416			HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
417			#else
418			eputs("qtAdd_tri():Found in parent but not children\n");
419			#endif
420			}
421			#endif
422			for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
423			{
424			id = QT_SET_NEXT_ELEM(optr);
425			qtTri_verts_from_id(id,r1,r2,r3);
426			found=qtAdd_tri(qtptr,id,r1,r2,r3,v1,v2,v3,n);
427			#ifdef DEBUG
428			if(!found)
429			eputs("qtAdd_tri():Reinsert-in parent but not children\n");
430			#endif
431			}
432			}
433			else
434			if(QT_SET_CNT(os) < QT_MAX_SET)
435			{
436			qtptr = qtaddelem(qtptr,id);
437			#if 0
438			{
439			int k;
440			qtgetset(os,*qtptr);
441			printf("\n%d:\n",os[0]);
442			for(k=1; k <= os[0];k++)
443			printf("%d ",os[k]);
444			printf("\n");
445			}
446			#endif
447			/*
448			insertelem(os,id);
449			*qtptr = fullnode(os);
450			*/
451			}
452			else
453			{
454			#ifdef DEBUG
455			eputs("qtAdd_tri():two many levels\n");
456			#endif
457			return(FALSE);
458			}
459			}
460			}
461			}
462			return(TRUE);
463			}
464
465
466			int
467			qtApply_to_tri_cells(qtptr,t1,t2,t3,v1,v2,v3,func,arg)
468			QUADTREE *qtptr;
469			FVECT t1,t2,t3;
470			FVECT v1,v2,v3;
471			int (*func)();
472			char *arg;
473			{
474			char test;
475			FVECT a,b,c;
476
477			/* test if triangle (t1,t2,t3) overlaps cell triangle (v1,v2,v3) */
478			test = spherical_tri_intersect(t1,t2,t3,v1,v2,v3);
479
480			/* If triangles do not overlap: done */
481			if(!test)
482			return(FALSE);
483
484			/* if this is tree: recurse */
485			if(QT_IS_TREE(*qtptr))
486			{
487			qtSubdivide_tri(v1,v2,v3,a,b,c);
488			qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t1,t2,t3,v1,a,c,func,arg);
489			qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t1,t2,t3,a,b,c,func,arg);
490			qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t1,t2,t3,a,v2,b,func,arg);
491			qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t1,t2,t3,c,b,v3,func,arg);
492			}
493			else
494			return(func(qtptr,arg));
495			}
496
497
498			int
499			qtRemove_tri(qtptr,id,t1,t2,t3,v1,v2,v3)
500			QUADTREE *qtptr;
501			int id;
502			FVECT t1,t2,t3;
503			FVECT v1,v2,v3;
504			{
505
506			char test;
507			int i;
508			FVECT a,b,c;
509			OBJECT os[MAXSET+1];
510
511			/* test if triangle (t1,t2,t3) overlaps cell triangle (v1,v2,v3) */
512			test = spherical_tri_intersect(t1,t2,t3,v1,v2,v3);
513
514			/* If triangles do not overlap: done */
515			if(!test)
516			return(FALSE);
517
518			/* if this is tree: recurse */
519			if(QT_IS_TREE(*qtptr))
520			{
521			qtSubdivide_tri(v1,v2,v3,a,b,c);
522			qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t1,t2,t3,v1,a,c);
523			qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t1,t2,t3,a,b,c);
524			qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t1,t2,t3,a,v2,b);
525			qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t1,t2,t3,c,b,v3);
526			}
527			else
528			{
529			if(QT_IS_EMPTY(*qtptr))
530			{
531			#ifdef DEBUG
532			eputs("qtRemove_tri(): triangle not found\n");
533			#endif
534			}
535			/* remove id from set */
536			else
537			{
538			qtgetset(os,*qtptr);
539			if(!inset(os,id))
540			{
541			#ifdef DEBUG
542			eputs("qtRemove_tri(): tri not in set\n");
543			#endif
544			}
545			else
546			{
547			qtptr = qtdelelem(qtptr,id);
548			#if 0
549			{
550			int k;
551			if(!QT_IS_EMPTY(*qtptr))
552			{qtgetset(os,*qtptr);
553			printf("\n%d:\n",os[0]);
554			for(k=1; k <= os[0];k++)
555			printf("%d ",os[k]);
556			printf("\n");
557			}
558
559			}
560			#endif
561			}
562			}
563			}
564			return(TRUE);
565			}