ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm_qtree.c
(Generate patch)

Comparing ray/src/hd/sm_qtree.c (file contents):
Revision 3.4 by gwlarson, Fri Sep 11 11:52:26 1998 UTC vs.
Revision 3.15 by greg, Wed Apr 23 00:52:34 2003 UTC

# Line 1 | Line 1
1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ SGI";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6
4   /*
5   *  sm_qtree.c: adapted from octree.c from radiance code
6   */
# Line 14 | Line 11 | static char SCCSid[] = "$SunId$ SGI";
11   */
12  
13   #include "standard.h"
14 <
14 > #include "sm_flag.h"
15   #include "sm_geom.h"
16 + #include "object.h"
17   #include "sm_qtree.h"
18  
19 < QUADTREE  *quad_block[QT_MAX_BLK];              /* our quadtree */
19 > QUADTREE  *quad_block[QT_MAX_BLK];        /* our quadtree */
20   static QUADTREE  quad_free_list = EMPTY;  /* freed octree nodes */
21 < static QUADTREE  treetop = 0;           /* next free node */
22 < int4 *quad_flag= NULL;
21 > static QUADTREE  treetop = 0;             /* next free node */
22 > int4 *quad_flag= NULL;                    /* quadtree flags */
23  
26 #ifdef TEST_DRIVER
27 extern FVECT Pick_v0[500],Pick_v1[500],Pick_v2[500];
28 extern int Pick_cnt,Pick_tri,Pick_samp;
29 extern FVECT Pick_point[500];
30 #endif
24  
25 + qtremovelast(qt,id)
26 + QUADTREE qt;
27 + OBJECT id;
28 + {
29 + #ifdef DEBUG
30 +  if(qtqueryset(qt)[(*(qtqueryset(qt)))] != id)
31 +    error(CONSISTENCY,"not removed\n");
32 + #endif
33 +  ((*(qtqueryset(qt)))--);
34 + }
35   QUADTREE
36   qtAlloc()                       /* allocate a quadtree */
37   {
# Line 46 | Line 49 | qtAlloc()                      /* allocate a quadtree */
49             return(EMPTY);
50          if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
51                          QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
52 <           return(EMPTY);
52 >           error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
53  
54 <        quad_flag = (int4 *)realloc((char *)quad_flag,
55 <                        (QT_BLOCK(freet)+1)*(QT_BLOCK_SIZE/8));
54 >        /* Realloc the per/node flags */
55 >        quad_flag = (int4 *)realloc((void *)quad_flag,
56 >                        (QT_BLOCK(freet)+1)*((QT_BLOCK_SIZE+7)>>3));
57          if (quad_flag == NULL)
58 <                return(EMPTY);
58 >           error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
59      }
60      treetop += 4;
61      return(freet);
# Line 60 | Line 64 | qtAlloc()                      /* allocate a quadtree */
64  
65   qtClearAllFlags()               /* clear all quadtree branch flags */
66   {
67 <        if (!treetop) return;
68 <        bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*(QT_BLOCK_SIZE/8));
67 >  if (!treetop)
68 >    return;
69 >  
70 >  /* Clear the node flags*/
71 >  bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*((QT_BLOCK_SIZE+7)>>3));
72 >  /* Clear set flags */
73 >  qtclearsetflags();
74   }
75  
67
76   qtFree(qt)                      /* free a quadtree */
77     register QUADTREE  qt;
78   {
# Line 85 | Line 93 | qtFree(qt)                     /* free a quadtree */
93   qtDone()                        /* free EVERYTHING */
94   {
95          register int    i;
96 <
96 >        
97          qtfreeleaves();
98          for (i = 0; i < QT_MAX_BLK; i++)
99          {
100              if (quad_block[i] == NULL)
101                  break;
102 <            free((char *)quad_block[i]);
102 >            free((void *)quad_block[i]);
103              quad_block[i] = NULL;
104          }
105 <        if (i) free((char *)quad_flag);
105 >        /* Free the flags */
106 >        if (i) free((void *)quad_flag);
107          quad_flag = NULL;
108          quad_free_list = EMPTY;
109          treetop = 0;
110   }
111  
112 < QUADTREE
113 < *qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
114 < QUADTREE *qtptr;
115 < double bcoord[3];
116 < FVECT t0,t1,t2;
117 < {
109 <  int i;
110 <  QUADTREE *child;
111 <  FVECT a,b,c;
112 <
113 <    if(QT_IS_TREE(*qtptr))
114 <    {  
115 <      i = bary_child(bcoord);
116 < #ifdef DEBUG_TEST_DRIVER
117 <      qtSubdivide_tri(Pick_v0[Pick_cnt-1],Pick_v1[Pick_cnt-1],
118 <                      Pick_v2[Pick_cnt-1],a,b,c);
119 <      qtNth_child_tri(Pick_v0[Pick_cnt-1],Pick_v1[Pick_cnt-1],
120 <                           Pick_v2[Pick_cnt-1],a,b,c,i,
121 <                           Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
122 <                           Pick_v2[Pick_cnt]);
123 <           Pick_cnt++;
124 < #endif
125 <
126 <      child = QT_NTH_CHILD_PTR(*qtptr,i);
127 <      if(t0)
128 <      {
129 <        qtSubdivide_tri(t0,t1,t2,a,b,c);
130 <        qtNth_child_tri(t0,t1,t2,a,b,c,i,t0,t1,t2);
131 <      }
132 <      return(qtLocate_leaf(child,bcoord,t0,t1,t2));
133 <    }
134 <    else
135 <      return(qtptr);
136 < }
137 <
138 <
139 < QUADTREE
140 < *qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
141 < QUADTREE *qtptr;
142 < FVECT v0,v1,v2;
143 < FVECT pt;
144 < FVECT t0,t1,t2;
145 < {
146 <    int d;
147 <    int i,x,y;
148 <    QUADTREE *child;
149 <    FVECT n,i_pt,a,b,c;
150 <    double pd,bcoord[3];
151 <
152 <    /* Determine if point lies within pyramid (and therefore
153 <       inside a spherical quadtree cell):GT_INTERIOR, on one of the
154 <       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
155 <       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
156 <       For each triangle edge: compare the
157 <       point against the plane formed by the edge and the view center
158 <    */
159 <    d = point_in_stri(v0,v1,v2,pt);  
160 <
161 <    
162 <    /* Not in this triangle */
163 <    if(!d)
164 <      return(NULL);
165 <
166 <    /* Will return lowest level triangle containing point: It the
167 <       point is on an edge or vertex: will return first associated
168 <       triangle encountered in the child traversal- the others can
169 <       be derived using triangle adjacency information
170 <    */
171 <    if(QT_IS_TREE(*qtptr))
172 <    {  
173 <      /* Find the intersection point */
174 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
175 <      intersect_vector_plane(pt,n,pd,NULL,i_pt);
176 <
177 <      i = max_index(n);
178 <      x = (i+1)%3;
179 <      y = (i+2)%3;
180 <      /* Calculate barycentric coordinates of i_pt */
181 <      bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
182 <
183 <      i = bary_child(bcoord);
184 <      child = QT_NTH_CHILD_PTR(*qtptr,i);
185 < #ifdef DEBUG_TEST_DRIVER
186 <      Pick_cnt = 0;
187 <      VCOPY(Pick_v0[0],v0);
188 <      VCOPY(Pick_v1[0],v1);
189 <      VCOPY(Pick_v2[0],v2);
190 <      Pick_cnt++;
191 <      qtSubdivide_tri(Pick_v0[Pick_cnt-1],Pick_v1[Pick_cnt-1],
192 <                      Pick_v2[Pick_cnt-1],a,b,c);
193 <      qtNth_child_tri(Pick_v0[Pick_cnt-1],Pick_v1[Pick_cnt-1],
194 <                      Pick_v2[Pick_cnt-1],a,b,c,i,
195 <                             Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
196 <                             Pick_v2[Pick_cnt]);
197 <           Pick_cnt++;
198 < #endif
199 <      if(t0)
200 <      {
201 <        qtSubdivide_tri(v0,v1,v2,a,b,c);
202 <        qtNth_child_tri(v0,v1,v2,a,b,c,i,t0,t1,t2);
203 <      }
204 <      return(qtLocate_leaf(child,bcoord,t0,t1,t2));
205 <    }
206 <    else
207 <    {
208 <        if(t0)
209 <        {
210 <            VCOPY(t0,v0);
211 <            VCOPY(t1,v1);
212 <            VCOPY(t2,v2);
213 <        }
214 <        return(qtptr);
215 <    }
216 < }
217 <
218 <
219 < QUADTREE
220 < qtSubdivide(qtptr)
221 < QUADTREE *qtptr;
222 < {
223 <  QUADTREE node;
224 <  node = qtAlloc();
225 <  QT_CLEAR_CHILDREN(node);
226 <  *qtptr = node;
227 <  return(node);
228 < }
229 <
230 <
231 < QUADTREE
232 < qtSubdivide_nth_child(qt,n)
233 < QUADTREE qt;
234 < int n;
235 < {
236 <  QUADTREE node;
237 <
238 <  node = qtSubdivide(&(QT_NTH_CHILD(qt,n)));
239 <  
240 <  return(node);
241 < }
242 <
243 < /* for triangle v0-v1-v2- returns a,b,c: children are:
244 <
245 <        v2                     0: v0,a,c
246 <        /\                     1: a,v1,b                  
247 <       /2 \                    2: c,b,v2
248 <     c/____\b                  3: b,c,a
249 <     /\    /\  
250 <    /0 \3 /1 \
251 <  v0____\/____\v1
252 <        a
112 > /*
113 > * bary_parent(coord,i)  : Set coord to be value of pt at one level up in
114 > *                         quadtree, given it was the ith child
115 > * BCOORD coord[3];      :barycentric coordinates of point, also will hold
116 > *                        result as side effect
117 > * int i;                :designates which child coord was
118   */
119 <
120 < qtSubdivide_tri(v0,v1,v2,a,b,c)
256 < FVECT v0,v1,v2;
257 < FVECT a,b,c;
258 < {
259 <    EDGE_MIDPOINT_VEC3(a,v0,v1);
260 <    EDGE_MIDPOINT_VEC3(b,v1,v2);
261 <    EDGE_MIDPOINT_VEC3(c,v2,v0);
262 < }
263 <
264 < qtNth_child_tri(v0,v1,v2,a,b,c,i,r0,r1,r2)
265 < FVECT v0,v1,v2;
266 < FVECT a,b,c;
119 > bary_parent(coord,i)
120 > BCOORD coord[3];
121   int i;
268 FVECT r0,r1,r2;
122   {
123 <
124 <  switch(i){
125 <  case 0:  
126 <  VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
123 >  switch(i) {
124 >  case 0:
125 >    /* update bary for child */
126 >    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
127 >    coord[1] >>= 1;
128 >    coord[2] >>= 1;
129      break;
130 <  case 1:  
131 <  VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
130 >  case 1:
131 >    coord[0] >>= 1;
132 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
133 >    coord[2] >>= 1;
134      break;
135 <  case 2:  
136 <  VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
135 >    
136 >  case 2:
137 >    coord[0] >>= 1;
138 >    coord[1] >>= 1;
139 >    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
140      break;
141 <  case 3:  
142 <  VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
141 >    
142 >  case 3:
143 >    coord[0] = MAXBCOORD4 - (coord[0] >> 1);
144 >    coord[1] = MAXBCOORD4 - (coord[1] >> 1);
145 >    coord[2] = MAXBCOORD4 - (coord[2] >> 1);
146      break;
147 + #ifdef DEBUG
148 +  default:
149 +    eputs("bary_parent():Invalid child\n");
150 +    break;
151 + #endif
152    }
153   }
154  
155 < /* Add triangle "id" to all leaf level cells that are children of
156 < quadtree pointed to by "qtptr" with cell vertices "t1,t2,t3"
157 < that it overlaps (vertex and edge adjacencies do not count
158 < as overlapping). If the addition of the triangle causes the cell to go over
159 < threshold- the cell is split- and the triangle must be recursively inserted
160 < into the new child cells: it is assumed that "v1,v2,v3" are normalized
161 < */
162 <
163 < int
164 < qtRoot_add_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
165 < QUADTREE *qtptr;
166 < FVECT q0,q1,q2;
299 < FVECT t0,t1,t2;
300 < int id;
301 < int n;
155 > /*
156 > * bary_from_child(coord,child,next): Given that coord was the (child) child
157 > *                      of the quadtree node, calculate what the (next) child
158 > *                      is at the same level.
159 > * BCOORD coord[3];   :At current level (child)th child of node,will also hold
160 > *                     result as side effect
161 > * int child,next;    :Which child coord was (child), and which one
162 > *                     is now desired(next)
163 > */
164 > bary_from_child(coord,child,next)
165 > BCOORD coord[3];
166 > int child,next;
167   {
168 <  int test;
169 <  int found;
168 > #ifdef DEBUG
169 >  if(child <0 || child > 3)
170 >  {
171 >    eputs("bary_from_child():Invalid child\n");
172 >    return;
173 >  }
174 >  if(next <0 || next > 3)
175 >  {
176 >    eputs("bary_from_child():Invalid next\n");
177 >    return;
178 >  }
179 > #endif
180 >  if(next == child)
181 >    return;
182  
183 <  test = stri_intersect(q0,q1,q2,t0,t1,t2);
184 <  if(!test)
185 <    return(FALSE);
186 <  
187 <  found = qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n);
188 <  
189 <  return(found);
183 >  switch(child){
184 >  case 0:
185 >      coord[0] = 0;
186 >      coord[1] = MAXBCOORD2 - coord[1];
187 >      coord[2] = MAXBCOORD2 - coord[2];
188 >      break;
189 >  case 1:
190 >      coord[0] = MAXBCOORD2 - coord[0];
191 >      coord[1] = 0;
192 >      coord[2] = MAXBCOORD2 - coord[2];
193 >      break;
194 >  case 2:
195 >      coord[0] = MAXBCOORD2 - coord[0];
196 >      coord[1] = MAXBCOORD2 - coord[1];
197 >      coord[2] = 0;
198 >    break;
199 >  case 3:
200 >    switch(next){
201 >    case 0:
202 >      coord[0] = 0;
203 >      coord[1] = MAXBCOORD2 - coord[1];
204 >      coord[2] = MAXBCOORD2 - coord[2];
205 >      break;
206 >    case 1:
207 >      coord[0] = MAXBCOORD2 - coord[0];
208 >      coord[1] = 0;
209 >      coord[2] = MAXBCOORD2 - coord[2];
210 >      break;
211 >    case 2:
212 >      coord[0] = MAXBCOORD2 - coord[0];
213 >      coord[1] = MAXBCOORD2 - coord[1];
214 >      coord[2] = 0;
215 >      break;
216 >    }
217 >    break;
218 >  }
219   }
220  
221 + /*
222 + * int
223 + * bary_child(coord):  Return which child coord is of node at current level
224 + *                     As side effect recalculate coord at next level
225 + * BCOORD coord[3];   :At current level coordinates of point, will also set
226 + *                     to coords at next level as side effect
227 + */
228   int
229 < qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
230 < QUADTREE *qtptr;
318 < FVECT q0,q1,q2;
319 < FVECT t0,t1,t2;
320 < int id;
321 < int n;
229 > bary_child(coord)
230 > BCOORD coord[3];
231   {
232 <  int i,index,test,found;
324 <  FVECT a,b,c;
325 <  OBJECT os[QT_MAXSET+1],*optr;
326 <  QUADTREE qt;
327 <  FVECT r0,r1,r2;
232 >  int i;
233  
234 <  found = 0;
235 <  /* if this is tree: recurse */
236 <  if(QT_IS_TREE(*qtptr))
237 <  {
238 <    n++;
239 <    qtSubdivide_tri(q0,q1,q2,a,b,c);
240 <    test = stri_intersect(t0,t1,t2,q0,a,c);
336 <    if(test)
337 <      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),q0,a,c,t0,t1,t2,id,n);
338 <    test = stri_intersect(t0,t1,t2,a,q1,b);
339 <    if(test)
340 <      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),a,q1,b,t0,t1,t2,id,n);
341 <    test = stri_intersect(t0,t1,t2,c,b,q2);
342 <    if(test)
343 <      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),c,b,q2,t0,t1,t2,id,n);
344 <    test = stri_intersect(t0,t1,t2,b,c,a);
345 <    if(test)
346 <      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),b,c,a,t0,t1,t2,id,n);
234 >  if(coord[0] > MAXBCOORD4)
235 >  {
236 >      /* update bary for child */
237 >      coord[0] = (coord[0]<< 1) - MAXBCOORD2;
238 >      coord[1] <<= 1;
239 >      coord[2] <<= 1;
240 >      return(0);
241    }
242    else
243 <  {
244 <      /* If this leave node emptry- create a new set */
245 <      if(QT_IS_EMPTY(*qtptr))
246 <        *qtptr = qtaddelem(*qtptr,id);
247 <      else
243 >    if(coord[1] > MAXBCOORD4)
244 >    {
245 >      coord[0] <<= 1;
246 >      coord[1] = (coord[1] << 1) - MAXBCOORD2;
247 >      coord[2] <<= 1;
248 >      return(1);
249 >    }
250 >    else
251 >      if(coord[2] > MAXBCOORD4)
252        {
253 <          /* If the set is too large: subdivide */
254 <        optr = qtqueryset(*qtptr);
255 <
256 <        if(QT_SET_CNT(optr) < QT_SET_THRESHOLD)
359 <          *qtptr = qtaddelem(*qtptr,id);
360 <          else
361 <          {
362 <            if (n < QT_MAX_LEVELS)
363 <            {
364 <                 /* If set size exceeds threshold: subdivide cell and
365 <                    reinsert set tris into cell
366 <                    */
367 <              qtgetset(os,*qtptr);      
368 <
369 <              n++;
370 <              qtfreeleaf(*qtptr);
371 <              qtSubdivide(qtptr);
372 <              found = qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n);
373 <
374 <              for(optr = QT_SET_PTR(os),i = QT_SET_CNT(os); i > 0; i--)
375 <                {
376 <                  id = QT_SET_NEXT_ELEM(optr);
377 <                  qtTri_from_id(id,NULL,NULL,NULL,r0,r1,r2,NULL,NULL,NULL);
378 <                  found=qtAdd_tri(qtptr,q0,q1,q2,r0,r1,r2,id,n);
379 < #ifdef DEBUG
380 <                  if(!found)
381 <                    eputs("qtAdd_tri():Reinsert-in parent but not children\n");
382 < #endif
383 <                }
384 <            }
385 <            else
386 <              if(QT_SET_CNT(optr) < QT_MAXSET)
387 <                  *qtptr = qtaddelem(*qtptr,id);
388 <              else
389 <                {
390 < #ifdef DEBUG
391 <                    eputs("qtAdd_tri():two many levels\n");
392 < #endif
393 <                    return(FALSE);
394 <                }
395 <          }
253 >        coord[0] <<= 1;
254 >        coord[1] <<= 1;
255 >        coord[2] = (coord[2] << 1) - MAXBCOORD2;
256 >        return(2);
257        }
258 <  }
259 <  return(TRUE);
258 >      else
259 >         {
260 >           coord[0] = MAXBCOORD2 - (coord[0] << 1);
261 >           coord[1] = MAXBCOORD2 - (coord[1] << 1);
262 >           coord[2] = MAXBCOORD2 - (coord[2] << 1);
263 >           return(3);
264 >         }
265   }
266  
267  
402 int
403 qtApply_to_tri_cells(qtptr,t0,t1,t2,v0,v1,v2,func,arg)
404 QUADTREE *qtptr;
405 FVECT t0,t1,t2;
406 FVECT v0,v1,v2;
407 int (*func)();
408 int *arg;
409 {
410  int test;
411  FVECT a,b,c;
268  
269 <  /* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
270 <  test = stri_intersect(t0,t1,t2,v0,v1,v2);
271 <
272 <  /* If triangles do not overlap: done */
417 <  if(!test)
418 <    return(FALSE);
419 <
420 <  /* if this is tree: recurse */
421 <  func(qtptr,arg);
422 <
423 <  if(QT_IS_TREE(*qtptr))
424 <  {
425 <     QT_SET_FLAG(*qtptr);
426 <     qtSubdivide_tri(v0,v1,v2,a,b,c);
427 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t0,t1,t2,v0,a,c,func,arg);
428 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t0,t1,t2,a,v1,b,func,arg);
429 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t0,t1,t2,c,b,v2,func,arg);
430 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t0,t1,t2,b,c,a,func,arg);
431 <  }
432 < }
433 <
434 < int
435 < qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
436 < QUADTREE *qtptr;
437 < int id;
438 < FVECT t0,t1,t2;
439 < FVECT v0,v1,v2;
269 > QUADTREE
270 > qtLocate(qt,bcoord)
271 > QUADTREE qt;
272 > BCOORD bcoord[3];
273   {
441  
442  int test;
274    int i;
444  FVECT a,b,c;
445  OBJECT os[QT_MAXSET+1];
446  
447  /* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
448  test = stri_intersect(t0,t1,t2,v0,v1,v2);
275  
276 <  /* If triangles do not overlap: done */
277 <  if(!test)
278 <    return(FALSE);
276 >  if(QT_IS_TREE(qt))
277 >  {  
278 >    i = bary_child(bcoord);
279  
280 <  /* if this is tree: recurse */
455 <  if(QT_IS_TREE(*qtptr))
456 <  {
457 <    qtSubdivide_tri(v0,v1,v2,a,b,c);
458 <    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c);
459 <    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b);
460 <    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2);
461 <    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a);
280 >    return(qtLocate(QT_NTH_CHILD(qt,i),bcoord));
281    }
282    else
283 <  {
465 <      if(QT_IS_EMPTY(*qtptr))
466 <      {
467 < #ifdef DEBUG    
468 <          eputs("qtRemove_tri(): triangle not found\n");
469 < #endif
470 <      }
471 <      /* remove id from set */
472 <      else
473 <      {
474 <          if(!qtinset(*qtptr,id))
475 <          {
476 < #ifdef DEBUG          
477 <              eputs("qtRemove_tri(): tri not in set\n");
478 < #endif
479 <          }
480 <          else
481 <          {
482 <            *qtptr = qtdelelem(*qtptr,id);
483 <        }
484 <      }
485 <  }
486 <  return(TRUE);
283 >    return(qt);
284   }
285  
489
286   int
287   move_to_nbr(b,db0,db1,db2,tptr)
288 < double b[3],db0,db1,db2;
288 > BCOORD b[3];
289 > BDIR db0,db1,db2;
290   double *tptr;
291   {
292    double t,dt;
293 +  BCOORD bc;
294    int nbr;
295 <
295 >  
296    nbr = -1;
297 +  *tptr = 0.0;
298    /* Advance to next node */
299 <  if(!ZERO(db0) && db0 < 0.0)
299 >  if(b[0]==0 && db0 < 0)
300 >    return(0);
301 >  if(b[1]==0 && db1 < 0)
302 >    return(1);
303 >  if(b[2]==0 && db2 < 0)
304 >    return(2);
305 >  if(db0 < 0)
306     {
307 <     t = -b[0]/db0;
307 >     t = ((double)b[0])/-db0;
308       nbr = 0;
309     }
310    else
311 <    t = FHUGE;
312 <  if(!ZERO(db1) && db1 < 0.0 )
311 >    t = MAXFLOAT;
312 >  if(db1 < 0)
313    {
314 <    dt = -b[1]/db1;
314 >     dt = ((double)b[1])/-db1;
315      if( dt < t)
316      {
317        t = dt;
318        nbr = 1;
319      }
320    }
321 <  if(!ZERO(db2) && db2 < 0.0 )
322 <    {
323 <      dt = -b[2]/db2;
324 <      if( dt < t)
321 >  if(db2 < 0)
322 >  {
323 >     dt = ((double)b[2])/-db2;
324 >       if( dt < t)
325        {
326          t = dt;
327          nbr = 2;
# Line 526 | Line 331 | double *tptr;
331    return(nbr);
332   }
333  
334 < int
335 < qtTrace_edge(qtptr,b,db,orig,dir,max_t,func,arg1,arg2)
336 <   QUADTREE *qtptr;
337 <   double b[3],db[3];
338 <   FVECT orig,dir;
339 <   double max_t;
340 <   int (*func)();
341 <   int *arg1,arg2;
334 > qtTrace_ray(qt,b,db0,db1,db2,nextptr,sign,sfactor,func,f)
335 >   QUADTREE qt;
336 >   BCOORD b[3];
337 >   BDIR db0,db1,db2;
338 >   int *nextptr;
339 >   int sign,sfactor;
340 >   FUNC func;
341 >   int *f;
342   {
538
343      int i,found;
344 <    QUADTREE *child;
345 <    int nbr,next;
344 >    QUADTREE child;
345 >    int nbr,next,w;
346      double t;
543 #ifdef DEBUG_TEST_DRIVER
544    
545    FVECT a1,b1,c1;
546    int Pick_parent = Pick_cnt-1;
547    qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
548                           Pick_v2[Pick_parent],a1,b1,c1);
347  
348 < #endif
551 <    if(QT_IS_TREE(*qtptr))
348 >    if(QT_IS_TREE(qt))
349      {
350 <        /* Find the appropriate child and reset the coord */
351 <        i = bary_child(b);
350 >      /* Find the appropriate child and reset the coord */
351 >      i = bary_child(b);
352  
353 <        QT_SET_FLAG(*qtptr);
353 >      for(;;)
354 >      {
355 >        child = QT_NTH_CHILD(qt,i);
356 >        if(i != 3)
357 >          qtTrace_ray(child,b,db0,db1,db2,nextptr,sign,sfactor+1,func,f);
358 >        else
359 >          /* If the center cell- must flip direction signs */
360 >          qtTrace_ray(child,b,-db0,-db1,-db2,nextptr,1-sign,sfactor+1,func,f);
361  
362 <        for(;;)
363 <        {
364 <           child = QT_NTH_CHILD_PTR(*qtptr,i);
365 <
366 <           if(i != 3)
367 <             {
368 <
369 <               db[0] *= 2.0;db[1] *= 2.0; db[2] *= 2.0;
370 <              nbr = qtTrace_edge(child,b,db,orig,dir,max_t,func,arg1,arg2);
371 <               db[0] *= 0.5;db[1] *= 0.5; db[2] *= 0.5;
372 <             }
373 <           else
374 <             {
375 <               db[0] *=-2.0;db[1] *= -2.0; db[2] *= -2.0;
572 <              /* If the center cell- must flip direction signs */
573 <              nbr = qtTrace_edge(child,b,db,orig,dir,max_t,func,arg1,arg2);
574 <              db[0] *=-0.5;db[1] *= -0.5; db[2] *= -0.5;
575 <             }
576 <           if(nbr == QT_DONE)
577 <              return(nbr);
578 <
579 <           /* If in same block: traverse */
580 <           if(i==3)
581 <              next = nbr;
582 <             else
583 <                if(nbr == i)
584 <                   next = 3;
585 <             else
586 <               {
587 <                 /* reset the barycentric coordinates in the parents*/
588 <                 bary_parent(b,i);
589 <                 /* Else pop up to parent and traverse from there */
590 <                 return(nbr);
591 <               }
592 <             bary_from_child(b,i,next);
593 <             i = next;
362 >        if(QT_FLAG_IS_DONE(*f))
363 >          return;
364 >        /* If in same block: traverse */
365 >        if(i==3)
366 >          next = *nextptr;
367 >        else
368 >          if(*nextptr == i)
369 >            next = 3;
370 >          else
371 >         {
372 >           /* reset the barycentric coordinates in the parents*/
373 >           bary_parent(b,i);
374 >           /* Else pop up to parent and traverse from there */
375 >           return(qt);
376           }
377 +        bary_from_child(b,i,next);
378 +        i = next;
379 +      }
380      }
381      else
382     {
383 < #ifdef DEBUG_TEST_DRIVER
384 <           qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
385 <                           Pick_v2[Pick_parent],a1,b1,c1,i,
386 <                           Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
387 <                           Pick_v2[Pick_cnt]);
388 <           Pick_cnt++;
389 < #endif
383 > #ifdef TEST_DRIVER
384 >       if(Pick_cnt < 500)
385 >          Pick_q[Pick_cnt++]=qt;
386 > #endif;
387 >       F_FUNC(func)(qt,F_ARGS(func),f);
388 >     if(QT_FLAG_IS_DONE(*f))
389 >       return;
390 >     /* Advance to next node */
391 >     nbr = move_to_nbr(b,db0,db1,db2,&t);
392  
393 <       if(func(qtptr,arg1,arg2) == QT_DONE)
394 <          return(QT_DONE);
395 <
396 <       /* Advance to next node */
397 <       /* NOTE: Optimize: should only have to check 1/2 */
398 <       nbr = move_to_nbr(b,db[0],db[1],db[2],&t);
399 <
400 <       if(t >= max_t)
401 <          return(QT_DONE);
402 <       if(nbr != -1)
403 <       {
617 <         b[0] += t * db[0];
618 <         b[1] += t * db[1];
619 <         b[2] += t * db[2];
620 <         db[0] *= (1.0 - t);
621 <         db[1] *= (1.0 - t);
622 <         db[2] *= (1.0 - t);
623 <       }
624 <       return(nbr);
393 >     if(nbr==-1)
394 >     {
395 >       QT_FLAG_SET_DONE(*f);
396 >       return;
397 >     }
398 >     b[0] += (BCOORD)(t *db0);
399 >     b[1] += (BCOORD)(t *db1);
400 >     b[2] += (BCOORD)(t *db2);
401 >     *nextptr = nbr;
402 >     return;
403 >    
404     }
405 <    
627 < }
405 > }    
406  
407 + #define TEST_INT(tl,th,d,q0,q1,h,l) \
408 +                  tl=d>q0;th=d<q1; if(tl&&th) goto Lfunc_call; \
409 +                 if(tl) if(l) goto Lfunc_call; else h=TRUE; \
410 +                 if(th) if(h) goto Lfunc_call; else l = TRUE; \
411 +                 if(th) if(h) goto Lfunc_call; else l = TRUE;
412  
413 < int
414 < qtTrace_ray(qtptr,b,db0,db1,db2,orig,dir,func,arg1,arg2)
415 <   QUADTREE *qtptr;
416 <   double b[3],db0,db1,db2;
417 <   FVECT orig,dir;
418 <   int (*func)();
419 <   int *arg1,arg2;
420 < {
413 > /* Leaf node: Do definitive test */
414 > QUADTREE
415 > qtLeaf_insert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
416 >                 scale, s0,s1,s2,sq0,sq1,sq2,f,n)
417 > int root;
418 > QUADTREE qt;
419 > BCOORD q0[3],q1[3],q2[3];
420 > BCOORD t0[3],t1[3],t2[3];
421 > BCOORD dt10[3],dt21[3],dt02[3];
422 > unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
423 > FUNC f;
424 > int n;
425 > {
426 >  double db;
427 >  unsigned int e0,e1,e2;
428 >  char al,ah,bl,bh,cl,ch,testl,testh;
429 >  char test_t0t1,test_t1t2,test_t2t0;
430 >  BCOORD a,b,c;
431  
432 <    int i,found;
433 <    QUADTREE *child;
434 <    int nbr,next;
435 <    double t;
436 < #ifdef DEBUG_TEST_DRIVER
437 <    
438 <    FVECT a1,b1,c1;
439 <    int Pick_parent = Pick_cnt-1;
440 <    qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
441 <                           Pick_v2[Pick_parent],a1,b1,c1);
432 >  /* First check if can trivial accept: if vertex in cell */
433 >  if(s0 & s1 & s2)
434 >  {
435 >    goto Lfunc_call;
436 >  }
437 >  /* Assumption: Could not trivial reject*/
438 >  /* IF any coords lie on edge- must be in if couldnt reject */
439 >  a = q1[0];b= q0[1]; c= q0[2];
440 >  if(t0[0]== a || t0[1] == b || t0[2] == c)
441 >  {
442 >    goto Lfunc_call;
443 >  }
444 >  if(t1[0]== a || t1[1] == b || t1[2] == c)
445 >  {
446 >    goto Lfunc_call;
447 >  }
448 >  if(t2[0]== a || t2[1] == b || t2[2] == c)
449 >  {
450 >    goto Lfunc_call;
451 >  }
452 >  /* Test for edge crossings */
453 >  /* Test t0t1,t1t2,t2t0 against a */
454 >  /* Calculate edge crossings */
455 >  e0  = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
456 >  /* See if edge can be trivially rejected from intersetion testing */
457 >  test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
458 >       ((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
459 >  bl=bh=0;
460 >  if(test_t0t1 && (e0 & 2) )
461 >  {
462 >    /* Must do double calculation to avoid overflow */
463 >    db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
464 >    TEST_INT(testl,testh,db,b,q1[1],bl,bh)
465 >  }
466 >  test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
467 >       ((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
468 >  if(test_t1t2 && (e0 & 1))
469 >  {    /* test t1t2 against a */
470 >      db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
471 >      TEST_INT(testl,testh,db,b,q1[1],bl,bh)
472 >  }
473 >  test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
474 >       ((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
475 >  if(test_t2t0 && (e0 & 4))
476 >  {
477 >      db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
478 >      TEST_INT(testl,testh,db,b,q1[1],bl,bh)
479 >  }
480 >  e1  = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
481 >  cl = ch = 0;
482 >  if(test_t0t1 && (e1 & 2))
483 >  {/* test t0t1 against b */
484 >      db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
485 >      TEST_INT(testl,testh,db,c,q2[2],cl,ch)
486 >  }
487 >  if(test_t1t2 && (e1 & 1))
488 >  {/* test t1t2 against b */
489 >    db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
490 >    TEST_INT(testl,testh,db,c,q2[2],cl,ch)
491 >  }
492 >  if(test_t2t0 && (e1 & 4))
493 >  {/* test t2t0 against b */
494 >    db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
495 >    TEST_INT(testl,testh,db,c,q2[2],cl,ch)
496 >  }
497 >
498 >  /* test edges against c */
499 >  e2  = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
500 >  al = ah = 0;
501 >  if(test_t0t1 && (e2 & 2))
502 >  { /* test t0t1 against c */
503 >    db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
504 >    TEST_INT(testl,testh,db,a,q0[0],al,ah)
505 >   }
506 >  if(test_t1t2 && (e2 & 1))
507 >  {
508 >    db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
509 >    TEST_INT(testl,testh,db,a,q0[0],al,ah)
510 >  }
511 >  if(test_t2t0 && (e2 & 4))
512 >  { /* test t2t0 against c */
513 >    db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
514 >    TEST_INT(testl,testh,db,a,q0[0],al,ah)
515 >  }
516 >  /* Only remaining case is if some edge trivially rejected */
517 >  if(!e0 || !e1 || !e2)
518 >    return(qt);
519  
520 < #endif
521 <    if(QT_IS_TREE(*qtptr))
522 <    {
523 <        /* Find the appropriate child and reset the coord */
524 <        i = bary_child(b);
520 >  /* Only one/none got tested - pick either of other two/any two */
521 >  /* Only need to test one edge */
522 >  if(!test_t0t1 && (e0 & 2))
523 >  {
524 >     db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
525 >     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
526 >  }
527 >  if(!test_t1t2 && (e0 & 1))
528 >    {/* test t1t2 against a */
529 >      db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
530 >      TEST_INT(testl,testh,db,b,q1[1],bl,bh)
531 >   }
532 >  if(!test_t2t0 && (e0 & 4))
533 >  {
534 >    db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
535 >    TEST_INT(testl,testh,db,b,q1[1],bl,bh)
536 >  }
537  
538 <        QT_SET_FLAG(*qtptr);
538 >  return(qt);
539  
540 <        for(;;)
541 <        {
542 <           child = QT_NTH_CHILD_PTR(*qtptr,i);
540 > Lfunc_call:
541 >  qt = f.func(f.argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
542 >              s0,s1,s2,sq0,sq1,sq2,0,f,n);
543 >  return(qt);
544  
545 <           if(i != 3)
663 <              nbr = qtTrace_ray(child,b,db0,db1,db2,orig,dir,func,arg1,arg2);
664 <           else
665 <               /* If the center cell- must flip direction signs */
666 <              nbr =qtTrace_ray(child,b,-db0,-db1,-db2,orig,dir,func,arg1,arg2);
667 <           if(nbr == QT_DONE)
668 <              return(nbr);
545 > }
546  
670           /* If in same block: traverse */
671           if(i==3)
672              next = nbr;
673             else
674                if(nbr == i)
675                   next = 3;
676             else
677               {
678                 /* reset the barycentric coordinates in the parents*/
679                 bary_parent(b,i);
680                 /* Else pop up to parent and traverse from there */
681                 return(nbr);
682               }
683             bary_from_child(b,i,next);
684             i = next;
685         }
686    }
687    else
688   {
689 #ifdef DEBUG_TEST_DRIVER
690           qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
691                           Pick_v2[Pick_parent],a1,b1,c1,i,
692                           Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
693                           Pick_v2[Pick_cnt]);
694           Pick_cnt++;
695 #endif
547  
697       if(func(qtptr,orig,dir,arg1,arg2) == QT_DONE)
698          return(QT_DONE);
548  
549 <       /* Advance to next node */
550 <       /* NOTE: Optimize: should only have to check 1/2 */
551 <       nbr = move_to_nbr(b,db0,db1,db2,&t);
549 > /* Leaf node: Do definitive test */
550 > QUADTREE
551 > qtLeaf_insert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
552 >                 scale, s0,s1,s2,sq0,sq1,sq2,f,n)
553 > int root;
554 > QUADTREE qt;
555 > BCOORD q0[3],q1[3],q2[3];
556 > BCOORD t0[3],t1[3],t2[3];
557 > BCOORD dt10[3],dt21[3],dt02[3];
558 > unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
559 > FUNC f;
560 > int n;
561 > {
562 >  double db;
563 >  unsigned int e0,e1,e2;
564 >  BCOORD a,b,c;
565 >  double p0[2], p1[2],cp;
566 >  char al,ah,bl,bh,cl,ch;
567 >  char testl,testh,test_t0t1,test_t1t2,test_t2t0;
568 >  unsigned int ls0,ls1,ls2;
569 >  
570 >  
571 >  /* May have gotten passed trivial reject if one/two vertices are ON */
572 >  a = q1[0];b= q0[1]; c= q0[2];
573 >  SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
574 >  
575 >  /* First check if can trivial accept: if vertex in cell */
576 >  if(ls0 & ls1 & ls2)
577 >  {
578 >    goto Lfunc_call;
579 >  }
580 >  if(ls0==0 || ls1 == 0 || ls2 ==0)
581 >    return(qt);
582 >  /* Assumption: Could not trivial reject*/
583 >  /* IF any coords lie on edge- must be in if couldnt reject */
584  
585 <       if(nbr != -1)
586 <       {
587 <         b[0] += t * db0;
588 <         b[1] += t * db1;
589 <         b[2] += t * db2;
590 <       }
591 <       return(nbr);
585 >  if(t0[0]== a || t0[1] == b || t0[2] == c)
586 >  {
587 >    goto Lfunc_call;
588 >  }
589 >  if(t1[0]== a || t1[1] == b || t1[2] == c)
590 >  {
591 >    goto Lfunc_call;
592 >  }
593 >  if(t2[0]== a || t2[1] == b || t2[2] == c)
594 >  {
595 >    goto Lfunc_call;
596 >  }
597 >  /* Test for edge crossings */
598 >  /* Test t0t1 against a,b,c */
599 >  /* First test if t0t1 can be trivially rejected */
600 >  /* If both edges are outside bounds- dont need to test */
601 >  
602 >  /* Test t0t1,t1t2,t2t0 against a */
603 >  test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
604 >       ((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
605 >  e0  = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
606 >  bl=bh=0;
607 >  /* Test t0t1,t1t2,t2t0 against a */
608 >  if(test_t0t1 && (e0 & 2) )
609 >  {
610 >      db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
611 >      TEST_INT(testl,testh,db,q1[1],b,bl,bh)
612 >  }
613 >  test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
614 >       ((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
615 >  if(test_t1t2 && (e0 & 1))
616 >  {    /* test t1t2 against a */
617 >      db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
618 >      TEST_INT(testl,testh,db,q1[1],b,bl,bh)
619 >  }
620 > test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
621 >       ((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
622 >  if(test_t2t0 && (e0 & 4))
623 >  {
624 >      db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
625 >      TEST_INT(testl,testh,db,q1[1],b,bl,bh)
626 >  }
627 >  cl = ch = 0;
628 >  e1  = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
629 >  if(test_t0t1 && (e1 & 2))
630 >  {/* test t0t1 against b */
631 >      db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
632 >      TEST_INT(testl,testh,db,q1[2],c,cl,ch)
633 >  }
634 >  if(test_t1t2 && (e1 & 1))
635 >  {/* test t1t2 against b */
636 >    db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
637 >    TEST_INT(testl,testh,db,q1[2],c,cl,ch)
638 >  }
639 >  if(test_t2t0 && (e1 & 4))
640 >  {/* test t2t0 against b */
641 >    db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
642 >    TEST_INT(testl,testh,db,q1[2],c,cl,ch)
643 >  }
644 >  al = ah = 0;
645 >  e2  = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
646 >  if(test_t0t1 && (e2 & 2))
647 >  { /* test t0t1 against c */
648 >    db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
649 >    TEST_INT(testl,testh,db,q0[0],a,al,ah)
650     }
651 <    
652 < }
651 >  if(test_t1t2 && (e2 & 1))
652 >  {
653 >    db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
654 >    TEST_INT(testl,testh,db,q0[0],a,al,ah)
655 >  }
656 >  if(test_t2t0 && (e2 & 4))
657 >  { /* test t2t0 against c */
658 >    db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
659 >    TEST_INT(testl,testh,db,q0[0],a,al,ah)
660 >  }
661 >  /* Only remaining case is if some edge trivially rejected */
662 >  if(!e0 || !e1 || !e2)
663 >    return(qt);
664  
665 < int
666 < qtRoot_trace_ray(qtptr,q0,q1,q2,orig,dir,func,arg1,arg2)
667 <   QUADTREE *qtptr;
668 <   FVECT q0,q1,q2;
669 <   FVECT orig,dir;
670 <   int (*func)();
671 <   int *arg1,arg2;
672 < {
673 <  int i,x,y,nbr;
674 <  QUADTREE *child;
675 <  FVECT n,c,i_pt,d;
676 <  double pd,b[3],db[3],t;
677 <    /* Determine if point lies within pyramid (and therefore
678 <       inside a spherical quadtree cell):GT_INTERIOR, on one of the
679 <       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
680 <       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
681 <       For each triangle edge: compare the
732 <       point against the plane formed by the edge and the view center
733 <    */
734 <  i = point_in_stri(q0,q1,q2,orig);  
735 <    
736 <  /* Not in this triangle */
737 <  if(!i)
738 <     return(INVALID);
739 <  /* Project the origin onto the root node plane */
665 >  /* Only one/none got tested - pick either of other two/any two */
666 >  /* Only need to test one edge */
667 >  if(!test_t0t1 && (e0 & 2))
668 >  {
669 >     db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
670 >     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
671 >  }
672 >  if(!test_t1t2 && (e0 & 1))
673 >    {/* test t1t2 against a */
674 >      db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
675 >      TEST_INT(testl,testh,db,q1[1],b,bl,bh)
676 >   }
677 >  if(!test_t2t0 && (e0 & 4))
678 >  {
679 >    db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
680 >    TEST_INT(testl,testh,db,q1[1],b,bl,bh)
681 >  }
682  
683 <  /* Find the intersection point of the origin */
684 <  tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
685 <  intersect_vector_plane(orig,n,pd,NULL,i_pt);
686 <  /* project the dir as well */
687 <  VADD(c,orig,dir);
688 <  intersect_vector_plane(c,n,pd,&t,c);
683 >  return(qt);
684 > Lfunc_call:
685 >  qt = f.func(f.argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
686 >              s0,s1,s2,sq0,sq1,sq2,1,f,n);
687 >  return(qt);
688 >  }
689  
748    /* map to 2d by dropping maximum magnitude component of normal */
749  i = max_index(n);
750  x = (i+1)%3;
751  y = (i+2)%3;
752  /* Calculate barycentric coordinates of orig */
753  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b);
754  /* Calculate barycentric coordinates of dir */
755  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db);
756  if(t < 0.0)
757     VSUB(db,b,db);
758  else
759     VSUB(db,db,b);
690  
691  
692 < #ifdef DEBUG_TEST_DRIVER
763 <    VCOPY(Pick_v0[Pick_cnt],q0);
764 <    VCOPY(Pick_v1[Pick_cnt],q1);
765 <    VCOPY(Pick_v2[Pick_cnt],q2);
766 <    Pick_cnt++;
767 < #endif
692 > /* ASSUMPTION: that triangle t0,t1,t2 intersects quad cell q0,q1,q2 */
693  
694 <      /* trace the ray starting with this node */
695 <    nbr = qtTrace_ray(qtptr,b,db[0],db[1],db[2],orig,dir,func,arg1,arg2);
696 <    return(nbr);
697 <    
698 < }
694 > /*-------q2--------- sq2
695 >        / \
696 > s1     /sc \ s0
697 >     qb_____qa
698 >     / \   / \
699 > \sq0/sa \ /sb \   /sq1
700 > \ q0____qc____q1/
701 >  \             /
702 >   \     s2    /
703 > */
704  
705 + int Find_id=0;
706  
707 < int
708 < qtRoot_trace_edge(qtptr,q0,q1,q2,orig,dir,max_t,func,arg1,arg2)
709 <   QUADTREE *qtptr;
710 <   FVECT q0,q1,q2;
711 <   FVECT orig,dir;
712 <   double max_t;
713 <   int (*func)();
714 <   int *arg1,arg2;
707 > QUADTREE
708 > qtInsert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
709 >            s0,s1,s2,sq0,sq1,sq2,f,n)
710 > int root;
711 > QUADTREE qt;
712 > BCOORD q0[3],q1[3],q2[3];
713 > BCOORD t0[3],t1[3],t2[3];
714 > BCOORD dt10[3],dt21[3],dt02[3];
715 > BCOORD scale;
716 > unsigned int s0,s1,s2,sq0,sq1,sq2;
717 > FUNC f;
718 > int n;
719   {
720 <  int i,x,y,nbr;
721 <  QUADTREE *child;
722 <  FVECT n,c,i_pt,d;
723 <  double pd,b[3],db[3],t;
724 <    /* Determine if point lies within pyramid (and therefore
725 <       inside a spherical quadtree cell):GT_INTERIOR, on one of the
726 <       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
727 <       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
728 <       For each triangle edge: compare the
794 <       point against the plane formed by the edge and the view center
720 >  BCOORD a,b,c;
721 >  BCOORD va[3],vb[3],vc[3];
722 >  unsigned int sa,sb,sc;
723 >
724 >  /* If a tree: trivial reject and recurse on potential children */
725 >  if(QT_IS_TREE(qt))
726 >  {
727 >    /* Test against new a edge: if entirely to left: only need
728 >       to visit child 0
729      */
730 <  i = point_in_stri(q0,q1,q2,orig);  
731 <    
732 <  /* Not in this triangle */
799 <  if(!i)
800 <     return(-1);
801 <  /* Project the origin onto the root node plane */
730 >    a = q1[0] + scale;
731 >    b = q0[1] + scale;
732 >    c = q0[2] + scale;
733  
734 <  /* Find the intersection point of the origin */
804 <  tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
805 <  intersect_vector_plane(orig,n,pd,NULL,i_pt);
806 <  /* project the dir as well */
807 <  VADD(c,orig,dir);
808 <  intersect_vector_plane(c,n,pd,&t,c);
734 >    SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
735  
736 <    /* map to 2d by dropping maximum magnitude component of normal */
737 <  i = max_index(n);
738 <  x = (i+1)%3;
739 <  y = (i+2)%3;
740 <  /* Calculate barycentric coordinates of orig */
741 <  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b);
742 <  /* Calculate barycentric coordinates of dir */
743 <  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db);
744 <  if(t < 0.0)
745 <     VSUB(db,b,db);
746 <  else
747 <     VSUB(db,db,b);
736 >    if(sa==7)
737 >    {
738 >      vb[1] = q0[1];
739 >      vc[2] = q0[2];
740 >      vc[1] = b;
741 >      vb[2] = c;
742 >      vb[0] = vc[0] = a;
743 >      QT_NTH_CHILD(qt,0) = qtInsert_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
744 >          vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
745 >      return(qt);
746 >    }
747 >   if(sb==7)
748 >   {
749 >     va[0] = q1[0];
750 >     vc[2] = q0[2];
751 >     va[1] = vc[1] = b;
752 >     va[2] = c;
753 >     vc[0] = a;
754 >     QT_NTH_CHILD(qt,1) = qtInsert_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,
755 >             t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
756 >     return(qt);
757 >   }
758 >   if(sc==7)
759 >   {
760 >     va[0] = q1[0];
761 >     vb[1] = q0[1];
762 >     va[1] = b;
763 >     va[2] = vb[2] = c;
764 >     vb[0] = a;
765 >     QT_NTH_CHILD(qt,2) = qtInsert_tri(root,QT_NTH_CHILD(qt,2),vb,va,
766 >       q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
767 >     return(qt);
768 >   }
769  
770 +   va[0] = q1[0];
771 +   vb[1] = q0[1];
772 +   vc[2] = q0[2];
773 +   va[1] = vc[1] = b;
774 +   va[2] = vb[2] = c;
775 +   vb[0] = vc[0] = a;
776 +   /* If outside of all edges: only need to Visit 3  */
777 +   if(sa==0 && sb==0 && sc==0)
778 +   {
779 +      QT_NTH_CHILD(qt,3) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,
780 +       vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
781 +      return(qt);
782 +   }
783  
784 < #ifdef DEBUG_TEST_DRIVER
785 <    VCOPY(Pick_v0[Pick_cnt],q0);
786 <    VCOPY(Pick_v1[Pick_cnt],q1);
787 <    VCOPY(Pick_v2[Pick_cnt],q2);
788 <    Pick_cnt++;
789 < #endif
790 <      /* trace the ray starting with this node */
791 <    nbr = qtTrace_edge(qtptr,b,db,orig,d,max_t,func,arg1,arg2);
792 <    return(nbr);
793 <    
784 >   if(sa)
785 >     QT_NTH_CHILD(qt,0) = qtInsert_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
786 >          t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
787 >   if(sb)
788 >     QT_NTH_CHILD(qt,1) = qtInsert_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
789 >              t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
790 >   if(sc)
791 >        QT_NTH_CHILD(qt,2) = qtInsert_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
792 >              t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
793 >   QT_NTH_CHILD(qt,3) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
794 >             t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
795 >   return(qt);
796 >  }
797 >  /* Leaf node: Do definitive test */
798 >  else
799 >    return(qt = qtLeaf_insert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
800 >                         scale,s0,s1,s2,sq0,sq1,sq2,f,n));
801   }
802  
803  
804 < qtVisit_tri_interior(qtptr,q0,q1,q2,t0,t1,t2,n,func,arg1,arg2)
805 <   QUADTREE *qtptr;
806 <   FVECT q0,q1,q2;
807 <   FVECT t0,t1,t2;
808 <   int n;
809 <   int (*func)();
810 <   int *arg1,arg2;
804 > QUADTREE
805 > qtInsert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
806 >            s0,s1,s2,sq0,sq1,sq2,f,n)
807 > int root;
808 > QUADTREE qt;
809 > BCOORD q0[3],q1[3],q2[3];
810 > BCOORD t0[3],t1[3],t2[3];
811 > BCOORD dt10[3],dt21[3],dt02[3];
812 > BCOORD scale;
813 > unsigned int s0,s1,s2,sq0,sq1,sq2;
814 > FUNC f;
815 > int n;
816   {
817 <    int i,found,test;
818 <    QUADTREE *child;
819 <    FVECT c0,c1,c2,a,b,c;
848 <    OBJECT os[QT_MAXSET+1],*optr;
849 <    int w;
850 <    
851 <    /* If qt Flag set, or qt vertices interior to t0t1t2-descend */
852 < tree_modified:
817 >  BCOORD a,b,c;
818 >  BCOORD va[3],vb[3],vc[3];
819 >  unsigned int sa,sb,sc;
820  
821 <    if(QT_IS_TREE(*qtptr))
822 <    {  
823 <        if(QT_IS_FLAG(*qtptr) ||  point_in_stri(t0,t1,t2,q0))
824 <        {
825 <            QT_SET_FLAG(*qtptr);
826 <            qtSubdivide_tri(q0,q1,q2,a,b,c);
827 <            /* descend to children */
828 <            for(i=0;i < 4; i++)
829 <            {
830 <                child = QT_NTH_CHILD_PTR(*qtptr,i);
831 <                qtNth_child_tri(q0,q1,q2,a,b,c,i,c0,c1,c2);
832 <                qtVisit_tri_interior(child,c0,c1,c2,t0,t1,t2,n+1,
833 <                                     func,arg1,arg2);
834 <            }
835 <        }
821 >  /* If a tree: trivial reject and recurse on potential children */
822 >  if(QT_IS_TREE(qt))
823 >  {
824 >    /* Test against new a edge: if entirely to left: only need
825 >       to visit child 0
826 >    */
827 >    a = q1[0] - scale;
828 >    b = q0[1] - scale;
829 >    c = q0[2] - scale;
830 >
831 >    SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
832 >
833 >    if(sa==0)
834 >    {
835 >      vb[1] = q0[1];
836 >      vc[2] = q0[2];
837 >      vc[1] = b;
838 >      vb[2] = c;
839 >      vb[0] = vc[0] = a;
840 >
841 >      QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
842 >          vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
843 >      return(qt);
844      }
845 <    else
845 >    if(sb==0)
846      {
847 <      /* NOTE THIS IN TRI TEST Could be replaced by a flag */
848 <      if(!QT_IS_EMPTY(*qtptr))
849 <      {
850 <         if(qtinset(*qtptr,arg2))
851 <           if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2)==QT_MODIFIED)
852 <             goto tree_modified;
853 <           else
854 <             return;
880 <      }
881 <      if(point_in_stri(t0,t1,t2,q0) )
882 <          if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2)==QT_MODIFIED)
883 <            goto tree_modified;
847 >      va[0] = q1[0];
848 >      vc[2] = q0[2];
849 >      va[1] = vc[1] = b;
850 >      va[2] = c;
851 >      vc[0] = a;
852 >      QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
853 >         t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
854 >      return(qt);
855      }
856 +    if(sc==0)
857 +    {
858 +      va[0] = q1[0];
859 +      vb[1] = q0[1];
860 +      va[1] = b;
861 +      va[2] = vb[2] = c;
862 +      vb[0] = a;
863 +      QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
864 +         q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
865 +      return(qt);
866 +    }
867 +    va[0] = q1[0];
868 +    vb[1] = q0[1];
869 +    vc[2] = q0[2];
870 +    va[1] = vc[1] = b;
871 +    va[2] = vb[2] = c;
872 +    vb[0] = vc[0] = a;
873 +    /* If outside of all edges: only need to Visit 3  */
874 +    if(sa==7 && sb==7 && sc==7)
875 +    {
876 +      QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,
877 +           vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
878 +      return(qt);
879 +    }
880 +    if(sa!=7)
881 +      QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
882 +           t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
883 +    if(sb!=7)
884 +      QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
885 +           t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
886 +    if(sc!=7)
887 +      QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
888 +           t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
889 +
890 +    QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
891 +              t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
892 +    return(qt);
893 +  }
894 +  /* Leaf node: Do definitive test */
895 +  else
896 +    return(qt = qtLeaf_insert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
897 +                         scale,s0,s1,s2,sq0,sq1,sq2,f,n));
898   }
899  
900  
901  
902  
903 + /*************************************************************************/
904 + /* Visit code for applying functions: NOTE THIS IS THE SAME AS INSERT CODE:
905 +  except sets flags: wanted insert to be as efficient as possible: so
906 +  broke into two sets of routines. Also traverses only to n levels.
907 + */
908  
909 <
910 < int
911 < qtVisit_tri_edges(qtptr,b,db,wptr,sfactor,func,arg1,arg2)
912 <   QUADTREE *qtptr;
913 <   double b[3],db[3][3];
914 <   int *wptr;
915 <   double sfactor;
916 <   int (*func)();
917 <   int *arg1,arg2;
909 > qtVisit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
910 >            s0,s1,s2,sq0,sq1,sq2,f,n)
911 > int root;
912 > QUADTREE qt;
913 > BCOORD q0[3],q1[3],q2[3];
914 > BCOORD t0[3],t1[3],t2[3];
915 > BCOORD dt10[3],dt21[3],dt02[3];
916 > BCOORD scale;
917 > unsigned int s0,s1,s2,sq0,sq1,sq2;
918 > FUNC f;
919 > int n;
920   {
921 <    int i,found;
922 <    QUADTREE *child;
923 <    int nbr,next,w;
904 <    double t;
905 < #ifdef DEBUG_TEST_DRIVER
906 <    FVECT a1,b1,c1;
907 <    int Pick_parent = Pick_cnt-1;
908 <    qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
909 <                           Pick_v2[Pick_parent],a1,b1,c1);
910 < #endif
921 >  BCOORD a,b,c;
922 >  BCOORD va[3],vb[3],vc[3];
923 >  unsigned int sa,sb,sc;
924  
925 <    if(QT_IS_TREE(*qtptr))
926 <    {
927 <        /* Find the appropriate child and reset the coord */
928 <        i = bary_child(b);
925 >  if(n == 0)
926 >    return;
927 >  /* If a tree: trivial reject and recurse on potential children */
928 >  if(QT_IS_TREE(qt))
929 >  {
930 >    QT_SET_FLAG(qt);
931  
932 <        QT_SET_FLAG(*qtptr);
932 >    /* Test against new a edge: if entirely to left: only need
933 >       to visit child 0
934 >    */
935 >    a = q1[0] + scale;
936 >    b = q0[1] + scale;
937 >    c = q0[2] + scale;
938  
939 <        for(;;)
920 <        {
921 <          w = *wptr;
922 <           child = QT_NTH_CHILD_PTR(*qtptr,i);
939 >    SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
940  
941 <           if(i != 3)
942 <             {
943 <
944 <               db[w][0] *= 2.0;db[w][1] *= 2.0; db[w][2] *= 2.0;
945 <               nbr = qtVisit_tri_edges(child,b,db,wptr,sfactor*2.0,
946 <                                      func,arg1,arg2);
947 <               w = *wptr;
948 <               db[w][0] *= 0.5;db[w][1] *= 0.5; db[w][2] *= 0.5;
949 <             }
950 <           else
934 <             {
935 <               db[w][0] *=-2.0;db[w][1] *= -2.0; db[w][2] *= -2.0;
936 <              /* If the center cell- must flip direction signs */
937 <               nbr = qtVisit_tri_edges(child,b,db,wptr,sfactor*(-2.0),
938 <                                 func,arg1,arg2);
939 <               w = *wptr;
940 <               db[w][0] *=-0.5;db[w][1] *= -0.5; db[w][2] *= -0.5;
941 <             }
942 <           if(nbr == QT_DONE)
943 <              return(nbr);
944 <
945 <           /* If in same block: traverse */
946 <           if(i==3)
947 <              next = nbr;
948 <             else
949 <                if(nbr == i)
950 <                   next = 3;
951 <             else
952 <               {
953 <                 /* reset the barycentric coordinates in the parents*/
954 <                 bary_parent(b,i);
955 <                 /* Else pop up to parent and traverse from there */
956 <                 return(nbr);
957 <               }
958 <             bary_from_child(b,i,next);
959 <             i = next;
960 <         }
941 >    if(sa==7)
942 >    {
943 >      vb[1] = q0[1];
944 >      vc[2] = q0[2];
945 >      vc[1] = b;
946 >      vb[2] = c;
947 >      vb[0] = vc[0] = a;
948 >      qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
949 >          vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
950 >      return;
951      }
952 <    else
952 >   if(sb==7)
953     {
954 < #ifdef DEBUG_TEST_DRIVER
955 <           qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
956 <                           Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
957 <                           Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
958 <           Pick_cnt++;
959 < #endif
954 >     va[0] = q1[0];
955 >     vc[2] = q0[2];
956 >     va[1] = vc[1] = b;
957 >     va[2] = c;
958 >     vc[0] = a;
959 >     qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,
960 >             t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
961 >     return;
962 >   }
963 >   if(sc==7)
964 >   {
965 >     va[0] = q1[0];
966 >     vb[1] = q0[1];
967 >     va[1] = b;
968 >     va[2] = vb[2] = c;
969 >     vb[0] = a;
970 >     qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,
971 >       q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
972 >     return;
973 >   }
974  
975 <       if(func(qtptr,arg1,arg2) == QT_DONE)
976 <          return(QT_DONE);
975 >   va[0] = q1[0];
976 >   vb[1] = q0[1];
977 >   vc[2] = q0[2];
978 >   va[1] = vc[1] = b;
979 >   va[2] = vb[2] = c;
980 >   vb[0] = vc[0] = a;
981 >   /* If outside of all edges: only need to Visit 3  */
982 >   if(sa==0 && sb==0 && sc==0)
983 >   {
984 >     qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,
985 >       vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
986 >      return;
987 >   }
988  
989 <       /* Advance to next node */
990 <       /* NOTE: Optimize: should only have to check 1/2 */
991 <       w = *wptr;
992 <       while(1)
993 <       {
994 <         nbr = move_to_nbr(b,db[w][0],db[w][1],db[w][2],&t);
989 >   if(sa)
990 >     qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
991 >          t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
992 >   if(sb)
993 >     qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
994 >              t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
995 >   if(sc)
996 >     qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
997 >              t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
998 >   qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
999 >             t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
1000 >  }
1001 >  /* Leaf node: Do definitive test */
1002 >  else
1003 >    qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1004 >                         scale,s0,s1,s2,sq0,sq1,sq2,f,n);
1005  
981         if(t >= 1.0)
982         {
983           if(w == 2)
984             return(QT_DONE);
985           b[0] += db[w][0];
986           b[1] += db[w][1];
987           b[2] += db[w][2];
988           w++;
989           db[w][0] *= sfactor;
990           db[w][1] *= sfactor;
991           db[w][2] *= sfactor;
992         }
993       else
994         if(nbr != INVALID)
995         {
996           b[0] += t * db[w][0];
997           b[1] += t * db[w][1];
998           b[2] += t * db[w][2];
999           db[w][0] *= (1.0 - t);
1000           db[w][1] *= (1.0 - t);
1001           db[w][2] *= (1.0 - t);
1002           *wptr = w;
1003           return(nbr);
1004         }
1005         else
1006           return(INVALID);
1007       }
1008   }
1009    
1006   }
1007  
1008  
1009 < int
1010 < qtRoot_visit_tri_edges(qtptr,q0,q1,q2,tri,dir,wptr,func,arg1,arg2)
1011 <   QUADTREE *qtptr;
1012 <   FVECT q0,q1,q2;
1013 <   FVECT tri[3],dir[3];
1014 <   int *wptr;
1015 <   int (*func)();
1016 <   int *arg1,arg2;
1009 > qtVisit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
1010 >            s0,s1,s2,sq0,sq1,sq2,f,n)
1011 > int root;
1012 > QUADTREE qt;
1013 > BCOORD q0[3],q1[3],q2[3];
1014 > BCOORD t0[3],t1[3],t2[3];
1015 > BCOORD dt10[3],dt21[3],dt02[3];
1016 > BCOORD scale;
1017 > unsigned int s0,s1,s2,sq0,sq1,sq2;
1018 > FUNC f;
1019 > int n;
1020   {
1021 <  int i,x,y,nbr,w;
1022 <  QUADTREE *child;
1023 <  FVECT n,c,i_pt,d;
1025 <  double pd,b[3][3],db[3][3],t;
1026 <    
1027 <  w = *wptr;
1021 >  BCOORD a,b,c;
1022 >  BCOORD va[3],vb[3],vc[3];
1023 >  unsigned int sa,sb,sc;
1024  
1025 <  /* Project the origin onto the root node plane */
1026 <
1027 <  /* Find the intersection point of the origin */
1028 <  tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
1033 <  /* map to 2d by dropping maximum magnitude component of normal */
1034 <  i = max_index(n);
1035 <  x = (i+1)%3;
1036 <  y = (i+2)%3;
1037 <  /* Calculate barycentric coordinates for current vertex */
1038 <  
1039 <  for(i=0;i < 3; i++)
1025 >  if(n==0)
1026 >    return;
1027 >  /* If a tree: trivial reject and recurse on potential children */
1028 >  if(QT_IS_TREE(qt))
1029    {
1030 <    /* If processing 3rd edge-dont need info for t1 */
1031 <    if(i==1 && w==2)
1032 <      continue;
1033 <    /* project the dir as well */
1034 <    intersect_vector_plane(tri[i],n,pd,NULL,i_pt);
1035 <    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[i]);  
1036 <    VADD(c,tri[i],dir[i]);
1048 <    intersect_vector_plane(c,n,pd,&t,c);
1049 <    /* Calculate barycentric coordinates of dir */
1050 <    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db[i]);
1051 <    if(t < 0.0)
1052 <      VSUB(db[i],b[i],db[i]);
1053 <    else
1054 <      VSUB(db[i],db[i],b[i]);
1055 <  }
1056 < #ifdef DEBUG_TEST_DRIVER
1057 <    VCOPY(Pick_v0[Pick_cnt],q0);
1058 <    VCOPY(Pick_v1[Pick_cnt],q1);
1059 <    VCOPY(Pick_v2[Pick_cnt],q2);
1060 <    Pick_cnt++;
1061 < #endif
1062 <      /* trace the ray starting with this node */
1063 <    nbr = qtVisit_tri_edges(qtptr,b[w],db,wptr,1.0,func,arg1,arg2);
1064 <    return(nbr);
1065 <    
1066 < }
1030 >    QT_SET_FLAG(qt);
1031 >    /* Test against new a edge: if entirely to left: only need
1032 >       to visit child 0
1033 >    */
1034 >    a = q1[0] - scale;
1035 >    b = q0[1] - scale;
1036 >    c = q0[2] - scale;
1037  
1038 +    SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
1039  
1040 <
1070 <
1071 < /* NOTE: SINCE DIR could be unit: then we could use integer math */
1072 < int
1073 < qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
1074 <                   db,wptr,t,sign,sfactor,func,arg1,arg2)
1075 <   QUADTREE *qtptr;
1076 <   double b[3],db0,db1,db2,db[3][3];
1077 <   int *wptr;
1078 <   double t[3];
1079 <   int sign;
1080 <   double sfactor;
1081 <   int (*func)();
1082 <   int *arg1,arg2;
1083 < {
1084 <    int i,found;
1085 <    QUADTREE *child;
1086 <    int nbr,next,w;
1087 <    double t_l,t_g;
1088 < #ifdef DEBUG_TEST_DRIVER
1089 <    FVECT a1,b1,c1;
1090 <    int Pick_parent = Pick_cnt-1;
1091 <    qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
1092 <                           Pick_v2[Pick_parent],a1,b1,c1);
1093 < #endif
1094 <    if(QT_IS_TREE(*qtptr))
1040 >    if(sa==0)
1041      {
1042 <        /* Find the appropriate child and reset the coord */
1043 <        i = bary_child(b);
1042 >      vb[1] = q0[1];
1043 >      vc[2] = q0[2];
1044 >      vc[1] = b;
1045 >      vb[2] = c;
1046 >      vb[0] = vc[0] = a;
1047 >      qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
1048 >          vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
1049 >      return;
1050 >    }
1051 >    if(sb==0)
1052 >    {
1053 >      va[0] = q1[0];
1054 >      vc[2] = q0[2];
1055 >      va[1] = vc[1] = b;
1056 >      va[2] = c;
1057 >      vc[0] = a;
1058 >      qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
1059 >         t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
1060 >      return;
1061 >    }
1062 >    if(sc==0)
1063 >    {
1064 >      va[0] = q1[0];
1065 >      vb[1] = q0[1];
1066 >      va[1] = b;
1067 >      va[2] = vb[2] = c;
1068 >      vb[0] = a;
1069 >      qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
1070 >         q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
1071 >      return;
1072 >    }
1073 >    va[0] = q1[0];
1074 >    vb[1] = q0[1];
1075 >    vc[2] = q0[2];
1076 >    va[1] = vc[1] = b;
1077 >    va[2] = vb[2] = c;
1078 >    vb[0] = vc[0] = a;
1079 >    /* If outside of all edges: only need to Visit 3  */
1080 >    if(sa==7 && sb==7 && sc==7)
1081 >    {
1082 >     qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,
1083 >           vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
1084 >      return;
1085 >    }
1086 >    if(sa!=7)
1087 >      qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
1088 >           t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
1089 >    if(sb!=7)
1090 >      qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
1091 >           t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
1092 >    if(sc!=7)
1093 >      qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
1094 >           t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
1095  
1096 <        QT_SET_FLAG(*qtptr);
1096 >    qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
1097 >              t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
1098 >    return;
1099 >  }
1100 >  /* Leaf node: Do definitive test */
1101 >  else
1102 >    qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1103 >                         scale,s0,s1,s2,sq0,sq1,sq2,f,n);
1104 > }
1105  
1101        for(;;)
1102        {
1103          w = *wptr;
1104           child = QT_NTH_CHILD_PTR(*qtptr,i);
1106  
1106           if(i != 3)
1107               nbr = qtVisit_tri_edges2(child,b,db0,db1,db2,
1108                                        db,wptr,t,sign,
1109                                        sfactor*2.0,func,arg1,arg2);
1110           else
1111             /* If the center cell- must flip direction signs */
1112             nbr = qtVisit_tri_edges2(child,b,-db0,-db1,-db2,
1113                                     db,wptr,t,1-sign,
1114                                     sfactor*2.0,func,arg1,arg2);
1107  
1108 <           if(nbr == QT_DONE)
1109 <             return(nbr);
1110 <           if(*wptr != w)
1111 <           {
1112 <             w = *wptr;
1113 <             db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
1114 <             if(sign)
1115 <               {  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
1116 <           }
1117 <           /* If in same block: traverse */
1118 <           if(i==3)
1119 <              next = nbr;
1120 <             else
1121 <                if(nbr == i)
1122 <                   next = 3;
1123 <             else
1124 <               {
1125 <                 /* reset the barycentric coordinates in the parents*/
1126 <                 bary_parent(b,i);
1127 <                 /* Else pop up to parent and traverse from there */
1136 <                 return(nbr);
1137 <               }
1138 <           bary_from_child(b,i,next);
1139 <           i = next;
1140 <        }
1141 <    }
1142 <    else
1143 <   {
1144 < #ifdef DEBUG_TEST_DRIVER
1145 <           qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
1146 <                           Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
1147 <                           Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
1148 <           Pick_cnt++;
1149 < #endif
1108 > qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1109 >                 scale, s0,s1,s2,sq0,sq1,sq2,f,n)
1110 > int root;
1111 > QUADTREE qt;
1112 > BCOORD q0[3],q1[3],q2[3];
1113 > BCOORD t0[3],t1[3],t2[3];
1114 > BCOORD dt10[3],dt21[3],dt02[3];
1115 > unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
1116 > FUNC f;
1117 > int n;
1118 > {
1119 >  double db;
1120 >  unsigned int e0,e1,e2;
1121 >  char al,ah,bl,bh,cl,ch,testl,testh;
1122 >  char test_t0t1,test_t1t2,test_t2t0;
1123 >  BCOORD a,b,c;
1124 >  
1125 >  /* First check if can trivial accept: if vertex in cell */
1126 >  if(s0 & s1 & s2)
1127 >    goto Lfunc_call;
1128  
1129 <       if(func(qtptr,arg1,arg2) == QT_DONE)
1130 <          return(QT_DONE);
1129 >  /* Assumption: Could not trivial reject*/
1130 >  /* IF any coords lie on edge- must be in if couldnt reject */
1131 >  a = q1[0];b= q0[1]; c= q0[2];
1132 >  if(t0[0]== a || t0[1] == b || t0[2] == c)
1133 >    goto Lfunc_call;
1134 >  if(t1[0]== a || t1[1] == b || t1[2] == c)
1135 >    goto Lfunc_call;
1136 >  if(t2[0]== a || t2[1] == b || t2[2] == c)
1137 >    goto Lfunc_call;
1138 >  
1139 >  /* Test for edge crossings */
1140 >  /* Test t0t1,t1t2,t2t0 against a */
1141 >  /* Calculate edge crossings */
1142 >  e0  = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
1143 >  /* See if edge can be trivially rejected from intersetion testing */
1144 >  test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
1145 >       ((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
1146 >  bl=bh=0;
1147 >  if(test_t0t1 && (e0 & 2) )
1148 >  {
1149 >    /* Must do double calculation to avoid overflow */
1150 >    db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1151 >    TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1152 >  }
1153 >  test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
1154 >       ((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
1155 >  if(test_t1t2 && (e0 & 1))
1156 >  {    /* test t1t2 against a */
1157 >      db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1158 >      TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1159 >  }
1160 >  test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
1161 >       ((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
1162 >  if(test_t2t0 && (e0 & 4))
1163 >  {
1164 >      db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1165 >      TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1166 >  }
1167 >  e1  = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
1168 >  cl = ch = 0;
1169 >  if(test_t0t1 && (e1 & 2))
1170 >  {/* test t0t1 against b */
1171 >      db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
1172 >      TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1173 >  }
1174 >  if(test_t1t2 && (e1 & 1))
1175 >  {/* test t1t2 against b */
1176 >    db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
1177 >    TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1178 >  }
1179 >  if(test_t2t0 && (e1 & 4))
1180 >  {/* test t2t0 against b */
1181 >    db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
1182 >    TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1183 >  }
1184 >
1185 >  /* test edges against c */
1186 >  e2  = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
1187 >  al = ah = 0;
1188 >  if(test_t0t1 && (e2 & 2))
1189 >  { /* test t0t1 against c */
1190 >    db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
1191 >    TEST_INT(testl,testh,db,a,q0[0],al,ah)
1192 >   }
1193 >  if(test_t1t2 && (e2 & 1))
1194 >  {
1195 >    db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1196 >    TEST_INT(testl,testh,db,a,q0[0],al,ah)
1197 >  }
1198 >  if(test_t2t0 && (e2 & 4))
1199 >  { /* test t2t0 against c */
1200 >    db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1201 >    TEST_INT(testl,testh,db,a,q0[0],al,ah)
1202 >  }
1203 >  /* Only remaining case is if some edge trivially rejected */
1204 >  if(!e0 || !e1 || !e2)
1205 >    return;
1206  
1207 <       /* Advance to next node */
1208 <       w = *wptr;
1209 <       while(1)
1210 <       {
1211 <         nbr = move_to_nbr(b,db0,db1,db2,&t_l);
1207 >  /* Only one/none got tested - pick either of other two/any two */
1208 >  /* Only need to test one edge */
1209 >  if(!test_t0t1 && (e0 & 2))
1210 >  {
1211 >     db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1212 >     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1213 >  }
1214 >  if(!test_t1t2 && (e0 & 1))
1215 >    {/* test t1t2 against a */
1216 >      db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1217 >      TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1218 >   }
1219 >  if(!test_t2t0 && (e0 & 4))
1220 >  {
1221 >    db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1222 >    TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1223 >  }
1224  
1225 <         t_g = t_l/sfactor;
1161 < #ifdef DEBUG
1162 <         if(t[w] <= 0.0)
1163 <           eputs("qtVisit_tri_edges2():negative t\n");
1164 < #endif
1165 <         if(t_g >= t[w])
1166 <         {
1167 <           if(w == 2)
1168 <             return(QT_DONE);
1225 >  return;
1226  
1227 <           b[0] += (t[w])*sfactor*db0;
1171 <           b[1] += (t[w])*sfactor*db1;
1172 <           b[2] += (t[w])*sfactor*db2;
1173 <           w++;
1174 <           db0 = db[w][0];
1175 <           db1 = db[w][1];
1176 <           db2 = db[w][2];
1177 <           if(sign)
1178 <          {  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
1179 <         }
1180 <       else
1181 <         if(nbr != INVALID)
1182 <         {
1183 <           b[0] += t_l * db0;
1184 <           b[1] += t_l * db1;
1185 <           b[2] += t_l * db2;
1227 > Lfunc_call:
1228  
1229 <           t[w] -= t_g;
1230 <           *wptr = w;
1231 <           return(nbr);
1232 <         }
1191 <         else
1192 <           return(INVALID);
1193 <       }
1194 <   }
1195 <    
1229 >  f.func(f.argptr,root,qt,n);
1230 >  if(!QT_IS_EMPTY(qt))
1231 >    QT_LEAF_SET_FLAG(qt);
1232 >
1233   }
1234  
1235  
1199 int
1200 qtRoot_visit_tri_edges2(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2)
1201   QUADTREE *qtptr;
1202   FVECT q0,q1,q2;
1203   FVECT tri[3],i_pt;
1204   int *wptr;
1205   int (*func)();
1206   int *arg1,arg2;
1207 {
1208  int x,y,z,nbr,w,i,j;
1209  QUADTREE *child;
1210  FVECT n,c,d,v[3];
1211  double pd,b[4][3],db[3][3],et[3],t[3],exit_pt;
1212    
1213  w = *wptr;
1236  
1237 <  /* Project the origin onto the root node plane */
1237 > /* Leaf node: Do definitive test */
1238 > QUADTREE
1239 > qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1240 >                 scale, s0,s1,s2,sq0,sq1,sq2,f,n)
1241 > int root;
1242 > QUADTREE qt;
1243 > BCOORD q0[3],q1[3],q2[3];
1244 > BCOORD t0[3],t1[3],t2[3];
1245 > BCOORD dt10[3],dt21[3],dt02[3];
1246 > unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
1247 > FUNC f;
1248 > int n;
1249 > {
1250 >  double db;
1251 >  unsigned int e0,e1,e2;
1252 >  BCOORD a,b,c;
1253 >  double p0[2], p1[2],cp;
1254 >  char al,ah,bl,bh,cl,ch;
1255 >  char testl,testh,test_t0t1,test_t1t2,test_t2t0;
1256 >  unsigned int ls0,ls1,ls2;
1257 >  
1258 >  /* May have gotten passed trivial reject if one/two vertices are ON */
1259 >  a = q1[0];b= q0[1]; c= q0[2];
1260 >  SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
1261 >  
1262 >  /* First check if can trivial accept: if vertex in cell */
1263 >  if(ls0 & ls1 & ls2)
1264 >    goto Lfunc_call;
1265  
1266 <  /* Find the intersection point of the origin */
1267 <  tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
1268 <  /* map to 2d by dropping maximum magnitude component of normal */
1269 <  z = max_index(n);
1270 <  x = (z+1)%3;
1271 <  y = (z+2)%3;
1272 <  /* Calculate barycentric coordinates for current vertex */
1273 <  if(w != -1)
1266 >  if(ls0==0 || ls1 == 0 || ls2 ==0)
1267 >    return;
1268 >  /* Assumption: Could not trivial reject*/
1269 >  /* IF any coords lie on edge- must be in if couldnt reject */
1270 >
1271 >  if(t0[0]== a || t0[1] == b || t0[2] == c)
1272 >    goto Lfunc_call;
1273 >  if(t1[0]== a || t1[1] == b || t1[2] == c)
1274 >    goto Lfunc_call;
1275 >  if(t2[0]== a || t2[1] == b || t2[2] == c)
1276 >    goto Lfunc_call;
1277 >
1278 >  /* Test for edge crossings */
1279 >  /* Test t0t1 against a,b,c */
1280 >  /* First test if t0t1 can be trivially rejected */
1281 >  /* If both edges are outside bounds- dont need to test */
1282 >  
1283 >  /* Test t0t1,t1t2,t2t0 against a */
1284 >  test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
1285 >       ((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
1286 >  e0  = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
1287 >  bl=bh=0;
1288 >  /* Test t0t1,t1t2,t2t0 against a */
1289 >  if(test_t0t1 && (e0 & 2) )
1290    {
1291 <    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);  
1292 <    intersect_vector_plane(tri[w],n,pd,&(et[w]),v[w]);
1228 <    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[w][x],v[w][y],b[w]);  
1291 >      db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
1292 >      TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1293    }
1294 <  else
1295 <  /* Just starting: b[0] is the origin point: guaranteed to be valid b*/
1294 >  test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
1295 >       ((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
1296 >  if(test_t1t2 && (e0 & 1))
1297 >  {    /* test t1t2 against a */
1298 >      db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1299 >      TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1300 >  }
1301 >  test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
1302 >       ((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
1303 >  if(test_t2t0 && (e0 & 4))
1304    {
1305 <    w = 0;
1306 <    intersect_vector_plane(tri[0],n,pd,&(et[0]),v[0]);
1235 <    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);  
1236 <    VCOPY(b[3],b[0]);
1305 >      db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1306 >      TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1307    }
1308 +  cl = ch = 0;
1309 +  e1  = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
1310 +  if(test_t0t1 && (e1 & 2))
1311 +  {/* test t0t1 against b */
1312 +      db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
1313 +      TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1314 +  }
1315 +  if(test_t1t2 && (e1 & 1))
1316 +  {/* test t1t2 against b */
1317 +    db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
1318 +    TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1319 +  }
1320 +  if(test_t2t0 && (e1 & 4))
1321 +  {/* test t2t0 against b */
1322 +    db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
1323 +    TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1324 +  }
1325 +  al = ah = 0;
1326 +  e2  = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
1327 +  if(test_t0t1 && (e2 & 2))
1328 +  { /* test t0t1 against c */
1329 +    db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
1330 +    TEST_INT(testl,testh,db,q0[0],a,al,ah)
1331 +   }
1332 +  if(test_t1t2 && (e2 & 1))
1333 +  {
1334 +    db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1335 +    TEST_INT(testl,testh,db,q0[0],a,al,ah)
1336 +  }
1337 +  if(test_t2t0 && (e2 & 4))
1338 +  { /* test t2t0 against c */
1339 +    db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1340 +    TEST_INT(testl,testh,db,q0[0],a,al,ah)
1341 +  }
1342 +  /* Only remaining case is if some edge trivially rejected */
1343 +  if(!e0 || !e1 || !e2)
1344 +    return;
1345  
1346 <
1347 <  j = (w+1)%3;
1348 <  intersect_vector_plane(tri[j],n,pd,&(et[j]),v[j]);
1242 <  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);  
1243 <  if(et[j] < 0.0)
1346 >  /* Only one/none got tested - pick either of other two/any two */
1347 >  /* Only need to test one edge */
1348 >  if(!test_t0t1 && (e0 & 2))
1349    {
1350 <      VSUB(db[w],b[3],b[j]);
1351 <      t[w] = FHUGE;
1350 >     db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1351 >     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1352    }
1353 <  else
1354 < {
1355 <   VSUB(db[w],b[j],b[3]);
1356 <   t[w] = 1.0;
1252 <   move_to_nbr(b[3],db[w][0],db[w][1],db[w][2],&exit_pt);
1253 <   if(exit_pt >= 1.0)
1254 <   {
1255 <     for(;j < 3;j++)
1256 <      {
1257 <          i = (j+1)%3;
1258 <          if(i!= w)
1259 <            {
1260 <              intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
1261 <              bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
1262 <                     v[i][y],b[i]);
1263 <            }
1264 <          if(et[i] < 0.0)
1265 <             {
1266 <                 VSUB(db[j],b[j],b[i]);
1267 <                 t[j] = FHUGE;
1268 <                 break;
1269 <             }
1270 <          else
1271 <             {
1272 <                 VSUB(db[j],b[i],b[j]);
1273 <                 t[j] = 1.0;
1274 <             }
1275 <          move_to_nbr(b[j],db[j][0],db[j][1],db[j][2],&exit_pt);
1276 <          if(exit_pt < 1.0)
1277 <             break;
1278 <      }
1353 >  if(!test_t1t2 && (e0 & 1))
1354 >    {/* test t1t2 against a */
1355 >      db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1356 >      TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1357     }
1358 < }
1359 <  *wptr = w;
1360 <  /* trace the ray starting with this node */
1361 <    nbr = qtVisit_tri_edges2(qtptr,b[3],db[w][0],db[w][1],db[w][2],
1362 <                             db,wptr,t,0,1.0,func,arg1,arg2);
1285 <    if(nbr != INVALID && nbr != QT_DONE)
1286 <    {
1287 <      i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
1288 <      i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
1289 <      i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
1290 <    }
1291 <    return(nbr);
1292 <    
1293 < }
1358 >  if(!test_t2t0 && (e0 & 4))
1359 >  {
1360 >    db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1361 >    TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1362 >  }
1363  
1364 +  return;
1365  
1366 + Lfunc_call:
1367 +  f.func(f.argptr,root,qt,n);
1368 +  if(!QT_IS_EMPTY(qt))
1369 +    QT_LEAF_SET_FLAG(qt);
1370 + }
1371  
1372  
1373 + QUADTREE
1374 + qtInsert_point(root,qt,parent,q0,q1,q2,bc,scale,f,n,doneptr)
1375 + int root;
1376 + QUADTREE qt,parent;
1377 + BCOORD q0[3],q1[3],q2[3],bc[3],scale;
1378 + FUNC f;
1379 + int n,*doneptr;
1380 + {
1381 +  BCOORD a,b,c;
1382 +  BCOORD va[3],vb[3],vc[3];
1383  
1384 +  if(QT_IS_TREE(qt))
1385 +  {  
1386 +    a = q1[0] + scale;
1387 +    b = q0[1] + scale;
1388 +    c = q0[2] + scale;
1389 +    if(bc[0] > a)
1390 +    {
1391 +      vc[0] = a; vc[1] = b; vc[2] = q0[2];
1392 +      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1393 +      QT_NTH_CHILD(qt,0) = qtInsert_point(root,QT_NTH_CHILD(qt,0),
1394 +                                  qt,q0,vc,vb,bc,scale>>1,f,n+1,doneptr);
1395 +      if(!*doneptr)
1396 +      {
1397 +        f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1398 +        if(!*doneptr)
1399 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1400 +        if(!*doneptr)
1401 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1402 +      }      
1403 +      return(qt);
1404 +    }
1405 +    if(bc[1] > b)
1406 +    {
1407 +      vc[0] = a; vc[1] = b; vc[2] = q0[2];
1408 +      va[0] = q1[0]; va[1] = b; va[2] = c;
1409 +      QT_NTH_CHILD(qt,1) =qtInsert_point(root,QT_NTH_CHILD(qt,1),
1410 +                                 qt,vc,q1,va,bc,scale >>1,f,n+1,doneptr);
1411 +      if(!*doneptr)
1412 +      {
1413 +        f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1414 +        if(!*doneptr)
1415 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1416 +        if(!*doneptr)
1417 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1418 +      }      
1419 +      return(qt);
1420 +    }
1421 +    if(bc[2] > c)
1422 +    {
1423 +      va[0] = q1[0]; va[1] = b; va[2] = c;
1424 +      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1425 +      QT_NTH_CHILD(qt,2) = qtInsert_point(root,QT_NTH_CHILD(qt,2),
1426 +                                  qt,vb,va,q2,bc,scale>>1,f,n+1,doneptr);
1427 +      if(!*doneptr)
1428 +      {
1429 +        f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1430 +        if(!*doneptr)
1431 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1432 +        if(!*doneptr)
1433 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1434 +      }
1435 +      return(qt);
1436 +    }
1437 +    else
1438 +    {
1439 +      va[0] = q1[0]; va[1] = b; va[2] = c;
1440 +      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1441 +      vc[0] = a; vc[1] = b; vc[2] = q0[2];
1442 +      QT_NTH_CHILD(qt,3) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,3),
1443 +                                      qt,va,vb,vc,bc,scale>>1,f,n+1,doneptr);
1444 +      if(!*doneptr)
1445 +      {
1446 +        f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1447 +        if(!*doneptr)
1448 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1449 +        if(!*doneptr)
1450 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1451 +      }
1452 +      return(qt);
1453 +    }
1454 +  }
1455 +  else
1456 +  {
1457 +    qt = f.func(f.argptr,root,qt,parent,q0,q1,q2,bc,scale,0,f,n,doneptr);
1458 +    return(qt);
1459 +  }
1460 + }
1461  
1462  
1463 + QUADTREE
1464 + qtInsert_point_rev(root,qt,parent,q0,q1,q2,bc,scale,f,n,doneptr)
1465 + int root;
1466 + QUADTREE qt,parent;
1467 + BCOORD q0[3],q1[3],q2[3],bc[3];
1468 + BCOORD scale;
1469 + FUNC f;
1470 + int n,*doneptr;
1471 + {
1472 +  BCOORD a,b,c;
1473 +  BCOORD va[3],vb[3],vc[3];
1474 +
1475 +  if(QT_IS_TREE(qt))
1476 +  {  
1477 +    a = q1[0] - scale;
1478 +    b = q0[1] - scale;
1479 +    c = q0[2] - scale;
1480 +    if(bc[0] < a)
1481 +    {
1482 +      vc[0] = a; vc[1] = b; vc[2] = q0[2];
1483 +      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1484 +      QT_NTH_CHILD(qt,0) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,0),
1485 +                                  qt,q0,vc,vb,bc,scale>>1,f,n+1,doneptr);
1486 +      if(!*doneptr)
1487 +      {
1488 +        f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1489 +        if(!*doneptr)
1490 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1491 +        if(!*doneptr)
1492 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1493 +      }
1494 +      return(qt);
1495 +    }
1496 +    if(bc[1] < b)
1497 +    {
1498 +      vc[0] = a; vc[1] = b; vc[2] = q0[2];
1499 +      va[0] = q1[0]; va[1] = b; va[2] = c;
1500 +      QT_NTH_CHILD(qt,1) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,1),
1501 +                                   qt,vc,q1,va,bc,scale>>1,f,n+1,doneptr);
1502 +      if(!*doneptr)
1503 +      {
1504 +        f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1505 +        if(!*doneptr)
1506 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1507 +        if(!*doneptr)
1508 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1509 +      }
1510 +      return(qt);
1511 +    }
1512 +    if(bc[2] < c)
1513 +    {
1514 +      va[0] = q1[0]; va[1] = b; va[2] = c;
1515 +      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1516 +      QT_NTH_CHILD(qt,2) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,2),
1517 +                                qt,vb,va,q2,bc,scale>>1,f,n+1,doneptr);
1518 +      if(!*doneptr)
1519 +      {
1520 +        f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1521 +        if(!*doneptr)
1522 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1523 +        if(!*doneptr)
1524 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1525 +      }
1526 +      return(qt);
1527 +    }
1528 +    else
1529 +    {
1530 +      va[0] = q1[0]; va[1] = b; va[2] = c;
1531 +      vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1532 +      vc[0] = a; vc[1] = b; vc[2] = q0[2];
1533 +      QT_NTH_CHILD(qt,3) = qtInsert_point(root,QT_NTH_CHILD(qt,3),
1534 +                                   qt,va,vb,vc,bc,scale>>1,f,n+1,doneptr);
1535 +      if(!*doneptr)
1536 +      {
1537 +        f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1538 +        if(!*doneptr)
1539 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1540 +        if(!*doneptr)
1541 +          f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1542 +      }
1543 +      return(qt);
1544 +    }
1545 +  }
1546 +  else
1547 +  {
1548 +    qt = f.func(f.argptr,root,qt,parent,q0,q1,q2,bc,scale,1,f,n,doneptr);
1549 +    return(qt);
1550 +  }
1551 + }
1552  
1553  
1554  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines