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
qtLocate_leaf(qtptr,bcoord,type,which)
QUADTREE *qtptr;
double bcoord[3];
char *type,*which;
{
  int i;
  QUADTREE *child;

    if(QT_IS_TREE(*qtptr))
    {  
      i = bary2d_child(bcoord);
      child = QT_NTH_CHILD_PTR(*qtptr,i);
      return(qtLocate_leaf(child,bcoord,type,which));
    }
    else
      return(*qtptr);
}


QUADTREE
qtRoot_point_locate(qtptr,v0,v1,v2,pt,type,which)
QUADTREE *qtptr;
FVECT v0,v1,v2;
FVECT pt;
char *type,*which;
{
    char d,w;
    int i,x,y;
    QUADTREE *child;
    QUADTREE qt;
    FVECT n,i_pt;
    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)
    {
      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))
    {  
      /* 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);
      return(qtLocate_leaf(child,bcoord,type,which));
    }
    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;
         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 = qtRoot_point_locate(qtptr,v0,v1,v2,pt,type,which);

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

  /* 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);
  found = 0;

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