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.3 by gwlarson, Tue Aug 25 11:03:27 1998 UTC vs.
Revision 3.7 by gwlarson, Tue Oct 6 18:18:54 1998 UTC

# Line 14 | Line 14 | static char SCCSid[] = "$SunId$ SGI";
14   */
15  
16   #include "standard.h"
17 <
17 > #include "sm_flag.h"
18   #include "sm_geom.h"
19   #include "sm_qtree.h"
20 – #include "object.h"
20  
21   QUADTREE  *quad_block[QT_MAX_BLK];              /* our quadtree */
22   static QUADTREE  quad_free_list = EMPTY;  /* freed octree nodes */
23   static QUADTREE  treetop = 0;           /* next free node */
24 < int4 *quad_flag;
24 > int4 *quad_flag= NULL;
25  
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 +
31 +
32 + #endif
33 + int Incnt=0;
34 +
35   QUADTREE
36   qtAlloc()                       /* allocate a quadtree */
37   {
# Line 41 | 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 >        /* Realloc the per/node flags */
55          quad_flag = (int4 *)realloc((char *)quad_flag,
56 <                        (QT_BLOCK(freet)+1)*QT_BLOCK_SIZE/(8*sizeof(int4)));
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 54 | Line 64 | qtAlloc()                      /* allocate a quadtree */
64  
65   qtClearAllFlags()               /* clear all quadtree branch flags */
66   {
67 <        if (!treetop) return;
68 <        bzero((char *)quad_flag,
69 <                (QT_BLOCK(treetop-1)+1)*QT_BLOCK_SIZE/(8*sizeof(int4)));
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  
62 –
76   qtFree(qt)                      /* free a quadtree */
77     register QUADTREE  qt;
78   {
# Line 80 | 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)
# Line 88 | Line 102 | qtDone()                       /* free EVERYTHING */
102              free((char *)quad_block[i]);
103              quad_block[i] = NULL;
104          }
105 +        /* Free the flags */
106          if (i) free((char *)quad_flag);
107          quad_flag = NULL;
108          quad_free_list = EMPTY;
109          treetop = 0;
110   }
111  
97 –
112   QUADTREE
113 < qtCompress(qt)                  /* recursively combine nodes */
114 < register QUADTREE  qt;
113 > qtLocate_leaf(qt,bcoord)
114 > QUADTREE qt;
115 > BCOORD bcoord[3];
116   {
102 –    register int  i;
103 –    register QUADTREE  qres;
104 –
105 –    if (!QT_IS_TREE(qt))        /* not a tree */
106 –       return(qt);
107 –    qres = QT_NTH_CHILD(qt,0) = qtCompress(QT_NTH_CHILD(qt,0));
108 –    for (i = 1; i < 4; i++)
109 –       if((QT_NTH_CHILD(qt,i) = qtCompress(QT_NTH_CHILD(qt,i))) != qres)
110 –          qres = qt;
111 –    if(!QT_IS_TREE(qres))
112 –    {   /* all were identical leaves */
113 –        QT_NTH_CHILD(qt,0) = quad_free_list;
114 –        quad_free_list = qt;
115 –    }
116 –    return(qres);
117 – }
118 –
119 –
120 – QUADTREE
121 – *qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
122 – QUADTREE *qtptr;
123 – double bcoord[3];
124 – FVECT t0,t1,t2;
125 – {
117    int i;
127 –  QUADTREE *child;
128 –  FVECT a,b,c;
118  
119 <    if(QT_IS_TREE(*qtptr))
131 <    {  
132 <      i = bary2d_child(bcoord);
133 <      child = QT_NTH_CHILD_PTR(*qtptr,i);
134 <      if(t0)
135 <      {
136 <        qtSubdivide_tri(t0,t1,t2,a,b,c);
137 <        qtNth_child_tri(t0,t1,t2,a,b,c,i,t0,t1,t2);
138 <      }
139 <      return(qtLocate_leaf(child,bcoord,t0,t1,t2));
140 <    }
141 <    else
142 <      return(qtptr);
143 < }
144 <
145 <
146 <
147 < int
148 < qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n)
149 < QUADTREE *qtptr;
150 < double bcoord[3];
151 < int id;
152 < FVECT v0,v1,v2;
153 < int n;
154 < {
155 <  int i;
156 <  QUADTREE *child;
157 <  OBJECT os[MAXSET+1],*optr;
158 <  int found;
159 <  FVECT r0,r1,r2;
160 <  
161 <  if(QT_IS_TREE(*qtptr))
119 >  if(QT_IS_TREE(qt))
120    {  
121 <    QT_SET_FLAG(*qtptr);
122 <    i = bary2d_child(bcoord);
123 <    child = QT_NTH_CHILD_PTR(*qtptr,i);
166 <    return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,++n));
121 >    i = baryi_child(bcoord);
122 >
123 >    return(qtLocate_leaf(QT_NTH_CHILD(qt,i),bcoord));
124    }
125    else
126 <    {
170 <      /* If this leave node emptry- create a new set */
171 <      if(QT_IS_EMPTY(*qtptr))
172 <        *qtptr = qtaddelem(*qtptr,id);
173 <      else
174 <      {
175 <          qtgetset(os,*qtptr);
176 <          /* If the set is too large: subdivide */
177 <          if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
178 <            *qtptr = qtaddelem(*qtptr,id);
179 <          else
180 <          {
181 <            if (n < QT_MAX_LEVELS)
182 <            {
183 <                 /* If set size exceeds threshold: subdivide cell and
184 <                    reinsert set tris into cell
185 <                    */
186 <                 n++;
187 <                 qtfreeleaf(*qtptr);
188 <                 qtSubdivide(qtptr);
189 <                 QT_SET_FLAG(*qtptr);
190 <                 found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n);
191 < #if 0
192 <                 if(!found)
193 <                 {
194 < #ifdef TEST_DRIVER    
195 <      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
196 < #else
197 <                   eputs("qtAdd_tri():Found in parent but not children\n");
198 < #endif
199 <                 }
200 < #endif
201 <                 for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
202 <                 {
203 <                   id = QT_SET_NEXT_ELEM(optr);
204 <                   qtTri_verts_from_id(id,r0,r1,r2);
205 <                   found= qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,++n);
206 < #ifdef DEBUG
207 <                 if(!found)
208 <                    eputs("qtAdd_tri():Reinsert-in parent but not children\n");
209 < #endif
210 <               }
211 <             }
212 <            else
213 <              if(QT_SET_CNT(os) < QT_MAX_SET)
214 <                {
215 <                  *qtptr = qtaddelem(*qtptr,id);
216 <                }
217 <              else
218 <                {
219 < #ifdef DEBUG
220 <                    eputs("qtAdd_tri():two many levels\n");
221 < #endif
222 <                    return(FALSE);
223 <                }
224 <          }
225 <      }
226 <  }
227 <  return(TRUE);
126 >    return(qt);
127   }
128  
129 <
130 < int
131 < qtAdd_tri_from_point(qtptr,v0,v1,v2,pt,id)
132 < QUADTREE *qtptr;
133 < FVECT v0,v1,v2;
129 > /*
130 >   Return the quadtree node containing pt. It is assumed that pt is in
131 >   the root node qt with ws vertices q0,q1,q2 and plane equation peq.
132 > */
133 > QUADTREE
134 > qtRoot_point_locate(qt,q0,q1,q2,peq,pt)
135 > QUADTREE qt;
136 > FVECT q0,q1,q2;
137 > FPEQ peq;
138   FVECT pt;
236 – int id;
139   {
238 –    char d,w;
140      int i,x,y;
141 <    QUADTREE *child;
142 <    QUADTREE qt;
143 <    FVECT i_pt,n,a,b,c,r0,r1,r2;
243 <    double pd,bcoord[3];
244 <    OBJECT os[MAXSET+1],*optr;
245 <    int found;
246 <        
247 <    /* Determine if point lies within pyramid (and therefore
248 <       inside a spherical quadtree cell):GT_INTERIOR, on one of the
249 <       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
250 <       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
251 <       For each triangle edge: compare the
252 <       point against the plane formed by the edge and the view center
253 <    */
254 <    d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);  
141 >    FVECT i_pt;
142 >    double bcoord[3];
143 >    BCOORD bcoordi[3];
144  
145 <    /* Not in this triangle */
257 <    if(!d)
258 <      return(FALSE);
259 <
260 <    /* Will return lowest level triangle containing point: It the
145 >     /* Will return lowest level triangle containing point: It the
146         point is on an edge or vertex: will return first associated
147         triangle encountered in the child traversal- the others can
148         be derived using triangle adjacency information
149      */
150 <    if(QT_IS_TREE(*qtptr))
150 >    if(QT_IS_TREE(qt))
151      {  
267 –      QT_SET_FLAG(*qtptr);
152        /* Find the intersection point */
153 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
270 <      intersect_vector_plane(pt,n,pd,NULL,i_pt);
153 >      intersect_vector_plane(pt,peq,NULL,i_pt);
154  
155 <      i = max_index(n);
156 <      x = (i+1)%3;
274 <      y = (i+2)%3;
155 >      x = FP_X(peq);
156 >      y = FP_Y(peq);
157        /* Calculate barycentric coordinates of i_pt */
158 <      bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
158 >      bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],bcoord);
159 >      
160 >      /* convert to integer coordinate */
161 >      convert_dtol(bcoord,bcoordi);
162  
163 <      i = bary2d_child(bcoord);
279 <      child = QT_NTH_CHILD_PTR(*qtptr,i);
280 <      /* NOTE: Optimize: only subdivide for specified child */
281 <      qtSubdivide_tri(v0,v1,v2,a,b,c);
282 <      qtNth_child_tri(v0,v1,v2,a,b,c,i,v0,v1,v2);
283 <      return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,1));
284 <    }
285 <    else
286 <   {
287 <      /* If this leave node emptry- create a new set */
288 <      if(QT_IS_EMPTY(*qtptr))
289 <        *qtptr = qtaddelem(*qtptr,id);
290 <      else
291 <      {
292 <          qtgetset(os,*qtptr);
293 <          /* If the set is too large: subdivide */
294 <          if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
295 <            *qtptr = qtaddelem(*qtptr,id);
296 <          else
297 <          {
298 <                 /* If set size exceeds threshold: subdivide cell and
299 <                    reinsert set tris into cell
300 <                    */
301 <                 qtfreeleaf(*qtptr);
302 <                 qtSubdivide(qtptr);
303 <                 found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,1);
304 < #if 0
305 <                 if(!found)
306 <                 {
307 < #ifdef TEST_DRIVER    
308 <      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
309 < #else
310 <                   eputs("qtAdd_tri():Found in parent but not children\n");
311 < #endif
312 <                 }
313 < #endif
314 <                 for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
315 <                 {
316 <                   id = QT_SET_NEXT_ELEM(optr);
317 <                   qtTri_verts_from_id(id,r0,r1,r2);
318 <                   found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,1);
319 < #ifdef DEBUG
320 <                 if(!found)
321 <                    eputs("qtAdd_tri():Reinsert-in parent but not children\n");
322 < #endif
323 <               }
324 <             }
325 <      }
326 <  }
327 <  return(TRUE);
328 < }
163 >      i = baryi_child(bcoordi);
164  
165 <
331 < QUADTREE
332 < *qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
333 < QUADTREE *qtptr;
334 < FVECT v0,v1,v2;
335 < FVECT pt;
336 < FVECT t0,t1,t2;
337 < {
338 <    char d,w;
339 <    int i,x,y;
340 <    QUADTREE *child;
341 <    QUADTREE qt;
342 <    FVECT n,i_pt,a,b,c;
343 <    double pd,bcoord[3];
344 <
345 <    /* Determine if point lies within pyramid (and therefore
346 <       inside a spherical quadtree cell):GT_INTERIOR, on one of the
347 <       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
348 <       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
349 <       For each triangle edge: compare the
350 <       point against the plane formed by the edge and the view center
351 <    */
352 <    d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);  
353 <
354 <    /* Not in this triangle */
355 <    if(!d)
356 <    {
357 <      return(NULL);
165 >      return(qtLocate_leaf(QT_NTH_CHILD(qt,i),bcoordi));
166      }
359 –
360 –    /* Will return lowest level triangle containing point: It the
361 –       point is on an edge or vertex: will return first associated
362 –       triangle encountered in the child traversal- the others can
363 –       be derived using triangle adjacency information
364 –    */
365 –    if(QT_IS_TREE(*qtptr))
366 –    {  
367 –      /* Find the intersection point */
368 –      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
369 –      intersect_vector_plane(pt,n,pd,NULL,i_pt);
370 –
371 –      i = max_index(n);
372 –      x = (i+1)%3;
373 –      y = (i+2)%3;
374 –      /* Calculate barycentric coordinates of i_pt */
375 –      bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
376 –
377 –      i = bary2d_child(bcoord);
378 –      child = QT_NTH_CHILD_PTR(*qtptr,i);
379 –      if(t0)
380 –      {
381 –        qtSubdivide_tri(v0,v1,v2,a,b,c);
382 –        qtNth_child_tri(v0,v1,v2,a,b,c,i,t0,t1,t2);
383 –      }
384 –      return(qtLocate_leaf(child,bcoord,t0,t1,t2));
385 –    }
167      else
168 <      return(qtptr);
168 >      return(qt);
169   }
170  
390 – int
391 – qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
392 – QUADTREE *qtptr;
393 – FVECT v0,v1,v2;
394 – FVECT pt;
395 – char *type,*which;
396 – {
397 –    OBJECT os[MAXSET+1],*optr;
398 –    char d,w;
399 –    int i,id;
400 –    FVECT p0,p1,p2;
171  
402 –    qtptr = qtRoot_point_locate(qtptr,v0,v1,v2,pt,NULL,NULL,NULL);
172  
404 –    if(!qtptr || QT_IS_EMPTY(*qtptr))
405 –       return(EMPTY);
406 –    
407 –    /* Get the set */
408 –    qtgetset(os,*qtptr);
409 –    for (i = QT_SET_CNT(os),optr = QT_SET_PTR(os); i > 0; i--)
410 –    {
411 –        /* Find the triangle that pt falls in (NOTE:FOR now return first 1) */
412 –        id = QT_SET_NEXT_ELEM(optr);
413 –        qtTri_verts_from_id(id,p0,p1,p2);
414 –        d = test_single_point_against_spherical_tri(p0,p1,p2,pt,&w);  
415 –        if(d)
416 –        {
417 –          if(type)
418 –            *type = d;
419 –          if(which)
420 –            *which = w;
421 –          return(id);
422 –        }
423 –    }
424 –    return(EMPTY);
425 – }
173  
427 – QUADTREE
428 – qtSubdivide(qtptr)
429 – QUADTREE *qtptr;
430 – {
431 –  QUADTREE node;
432 –  node = qtAlloc();
433 –  QT_CLEAR_CHILDREN(node);
434 –  *qtptr = node;
435 –  return(node);
436 – }
437 –
438 –
439 – QUADTREE
440 – qtSubdivide_nth_child(qt,n)
441 – QUADTREE qt;
442 – int n;
443 – {
444 –  QUADTREE node;
445 –
446 –  node = qtSubdivide(&(QT_NTH_CHILD(qt,n)));
447 –  
448 –  return(node);
449 – }
450 –
174   /* for triangle v0-v1-v2- returns a,b,c: children are:
175  
176          v2                     0: v0,a,c
# Line 460 | Line 183 | int n;
183          a
184   */
185  
463 – qtSubdivide_tri(v0,v1,v2,a,b,c)
464 – FVECT v0,v1,v2;
465 – FVECT a,b,c;
466 – {
467 –    EDGE_MIDPOINT_VEC3(a,v0,v1);
468 –    EDGE_MIDPOINT_VEC3(b,v1,v2);
469 –    EDGE_MIDPOINT_VEC3(c,v2,v0);
470 – }
186  
187   qtNth_child_tri(v0,v1,v2,a,b,c,i,r0,r1,r2)
188   FVECT v0,v1,v2;
# Line 475 | Line 190 | FVECT a,b,c;
190   int i;
191   FVECT r0,r1,r2;
192   {
193 <  switch(i){
194 <  case 0:  
195 <  VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
196 <    break;
197 <  case 1:  
198 <  VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
199 <    break;
200 <  case 2:  
201 <  VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
202 <    break;
203 <  case 3:  
204 <  VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
205 <    break;
193 >
194 >  if(!a)
195 >  {
196 >    /* Caution: r's must not be equal to v's:will be incorrect */
197 >    switch(i){
198 >    case 0:  
199 >      VCOPY(r0,v0);
200 >      EDGE_MIDPOINT_VEC3(r1,v0,v1);
201 >      EDGE_MIDPOINT_VEC3(r2,v2,v0);
202 >      break;
203 >    case 1:  
204 >      EDGE_MIDPOINT_VEC3(r0,v0,v1);
205 >      VCOPY(r1,v1);    
206 >      EDGE_MIDPOINT_VEC3(r2,v1,v2);
207 >      break;
208 >    case 2:  
209 >      EDGE_MIDPOINT_VEC3(r0,v2,v0);
210 >      EDGE_MIDPOINT_VEC3(r1,v1,v2);
211 >      VCOPY(r2,v2);
212 >      break;
213 >    case 3:  
214 >      EDGE_MIDPOINT_VEC3(r0,v1,v2);
215 >      EDGE_MIDPOINT_VEC3(r1,v2,v0);
216 >      EDGE_MIDPOINT_VEC3(r2,v0,v1);
217 >      break;
218 >    }
219    }
220 +  else
221 +  {
222 +    switch(i){
223 +    case 0:  
224 +      VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
225 +      break;
226 +    case 1:  
227 +      VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
228 +      break;
229 +    case 2:  
230 +      VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
231 +      break;
232 +    case 3:  
233 +      VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
234 +      break;
235 +    }
236 +  }
237   }
238  
239   /* Add triangle "id" to all leaf level cells that are children of
# Line 499 | Line 244 | threshold- the cell is split- and the triangle must be
244   into the new child cells: it is assumed that "v1,v2,v3" are normalized
245   */
246  
247 < int
248 < qtRoot_add_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
249 < QUADTREE *qtptr;
250 < int id;
247 > QUADTREE
248 > qtRoot_add_tri(qt,q0,q1,q2,t0,t1,t2,id,n)
249 > QUADTREE qt;
250 > FVECT q0,q1,q2;
251   FVECT t0,t1,t2;
252 < FVECT v0,v1,v2;
508 < int n;
252 > int id,n;
253   {
254 <  char test;
255 <  int found;
254 >  if(stri_intersect(q0,q1,q2,t0,t1,t2))
255 >    qt = qtAdd_tri(qt,q0,q1,q2,t0,t1,t2,id,n);
256  
257 <  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
514 <  if(!test)
515 <    return(FALSE);
516 <  
517 <  found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
518 <  return(found);
257 >  return(qt);
258   }
259  
260 < int
261 < qtRoot_add_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n)
262 < QUADTREE *qtptr;
263 < int id;
260 > QUADTREE
261 > qtRoot_remove_tri(qt,q0,q1,q2,t0,t1,t2,id,n)
262 > QUADTREE qt;
263 > FVECT q0,q1,q2;
264   FVECT t0,t1,t2;
265 < FVECT v0,v1,v2;
527 < int n;
265 > int id,n;
266   {
529 –  char test;
530 –  int found;
267  
268 <  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
269 <  if(!test)
270 <    return(FALSE);
535 <  
536 <  found = qtAdd_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n);
537 <  return(found);
268 >  if(stri_intersect(q0,q1,q2,t0,t1,t2))
269 >    qt = qtRemove_tri(qt,q0,q1,q2,t0,t1,t2,id,n);
270 >  return(qt);
271   }
272  
273 < int
274 < qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
275 < QUADTREE *qtptr;
276 < int id;
273 >
274 > QUADTREE
275 > qtAdd_tri(qt,q0,q1,q2,t0,t1,t2,id,n)
276 > QUADTREE qt;
277 > FVECT q0,q1,q2;
278   FVECT t0,t1,t2;
279 < FVECT v0,v1,v2;
279 > int id;
280   int n;
281   {
548 –  
549 –  char test;
550 –  int i,index;
282    FVECT a,b,c;
283 <  OBJECT os[MAXSET+1],*optr;
553 <  QUADTREE qt;
554 <  int found;
283 >  OBJECT tset[QT_MAXSET+1],*optr,*tptr;
284    FVECT r0,r1,r2;
285 +  int i;
286  
557 –  found = 0;
287    /* if this is tree: recurse */
288 <  if(QT_IS_TREE(*qtptr))
288 >  if(QT_IS_TREE(qt))
289    {
290 +    QT_SET_FLAG(qt);
291      n++;
292 <    qtSubdivide_tri(v0,v1,v2,a,b,c);
563 <    test = spherical_tri_intersect(t0,t1,t2,v0,a,c);
564 <    if(test)
565 <      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c,n);
566 <    test = spherical_tri_intersect(t0,t1,t2,a,v1,b);
567 <    if(test)
568 <      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b,n);
569 <    test = spherical_tri_intersect(t0,t1,t2,c,b,v2);
570 <    if(test)
571 <      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2,n);
572 <    test = spherical_tri_intersect(t0,t1,t2,b,c,a);
573 <    if(test)
574 <      found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a,n);
292 >    qtSubdivide_tri(q0,q1,q2,a,b,c);
293  
294 < #if 0
295 <    if(!found)
296 <    {
297 < #ifdef TEST_DRIVER    
298 <      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
299 < #else
300 <      eputs("qtAdd_tri():Found in parent but not children\n");
301 < #endif
302 <    }
585 < #endif
294 >    if(stri_intersect(t0,t1,t2,q0,a,c))
295 >      QT_NTH_CHILD(qt,0) = qtAdd_tri(QT_NTH_CHILD(qt,0),q0,a,c,t0,t1,t2,id,n);
296 >    if(stri_intersect(t0,t1,t2,a,q1,b))
297 >       QT_NTH_CHILD(qt,1) = qtAdd_tri(QT_NTH_CHILD(qt,1),a,q1,b,t0,t1,t2,id,n);
298 >    if(stri_intersect(t0,t1,t2,c,b,q2))
299 >      QT_NTH_CHILD(qt,2) = qtAdd_tri(QT_NTH_CHILD(qt,2),c,b,q2,t0,t1,t2,id,n);
300 >    if(stri_intersect(t0,t1,t2,b,c,a))
301 >      QT_NTH_CHILD(qt,3) = qtAdd_tri(QT_NTH_CHILD(qt,3),b,c,a,t0,t1,t2,id,n);
302 >    return(qt);
303    }
304    else
305    {
306        /* If this leave node emptry- create a new set */
307 <      if(QT_IS_EMPTY(*qtptr))
308 <        *qtptr = qtaddelem(*qtptr,id);
307 >      if(QT_IS_EMPTY(qt))
308 >        qt = qtaddelem(qt,id);
309        else
310        {
594 –          qtgetset(os,*qtptr);
311            /* If the set is too large: subdivide */
312 <          if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
313 <            *qtptr = qtaddelem(*qtptr,id);
314 <          else
315 <          {
316 <            if (n < QT_MAX_LEVELS)
317 <            {
318 <                 /* If set size exceeds threshold: subdivide cell and
319 <                    reinsert set tris into cell
320 <                    */
321 <                 n++;
322 <                 qtfreeleaf(*qtptr);
323 <                 qtSubdivide(qtptr);
324 <                 found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
325 < #if 0
326 <                 if(!found)
327 <                 {
328 < #ifdef TEST_DRIVER    
329 <      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
330 < #else
615 <                   eputs("qtAdd_tri():Found in parent but not children\n");
616 < #endif
617 <                 }
618 < #endif
619 <                 for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
620 <                 {
621 <                   id = QT_SET_NEXT_ELEM(optr);
622 <                   qtTri_verts_from_id(id,r0,r1,r2);
623 <                   found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,n);
624 < #ifdef DEBUG
625 <                 if(!found)
626 <                    eputs("qtAdd_tri():Reinsert-in parent but not children\n");
627 < #endif
628 <               }
629 <             }
312 >        optr = qtqueryset(qt);
313 >
314 >        if(QT_SET_CNT(optr) < QT_SET_THRESHOLD)
315 >          qt = qtaddelem(qt,id);
316 >        else
317 >        {
318 >          if (n < QT_MAX_LEVELS)
319 >          {
320 >            /* If set size exceeds threshold: subdivide cell and
321 >               reinsert set tris into cell
322 >               */
323 >            /* CAUTION:If QT_SET_THRESHOLD << QT_MAXSET, and dont add
324 >               more than a few triangles before expanding: then are safe here
325 >               otherwise must check to make sure set size is < MAXSET,
326 >               or  qtgetset could overflow os.
327 >               */
328 >            tptr = qtqueryset(qt);
329 >            if(QT_SET_CNT(tptr) > QT_MAXSET)
330 >              tptr = (OBJECT *)malloc((QT_SET_CNT(tptr)+1)*sizeof(OBJECT));
331              else
332 <              if(QT_SET_CNT(os) < QT_MAX_SET)
333 <                {
334 <                  *qtptr = qtaddelem(*qtptr,id);
634 < #if 0
635 <                  {
636 <              int k;
637 <              qtgetset(os,*qtptr);
638 <              printf("\n%d:\n",os[0]);
639 <              for(k=1; k <= os[0];k++)
640 <                 printf("%d ",os[k]);
641 <              printf("\n");
642 <          }
643 < #endif            
644 <                  /*
645 <                    insertelem(os,id);
646 <                    *qtptr = fullnode(os);
647 <                  */
648 <                }
649 <              else
650 <                {
651 < #ifdef DEBUG
652 <                    eputs("qtAdd_tri():two many levels\n");
653 < #endif
654 <                    return(FALSE);
655 <                }
656 <          }
657 <      }
658 <  }
659 <  return(TRUE);
660 < }
332 >              tptr = tset;
333 >            if(!tptr)
334 >              goto memerr;
335  
336 +            qtgetset(tptr,qt);  
337 +            n++;
338 +            qtfreeleaf(qt);
339 +            qtSubdivide(qt);
340 +            qt = qtAdd_tri(qt,q0,q1,q2,t0,t1,t2,id,n);
341  
342 < int
343 < qtApply_to_tri_cells(qtptr,t0,t1,t2,v0,v1,v2,func,arg)
344 < QUADTREE *qtptr;
345 < FVECT t0,t1,t2;
346 < FVECT v0,v1,v2;
347 < int (*func)();
348 < char *arg;
349 < {
350 <  char test;
351 <  FVECT a,b,c;
352 <
353 <  /* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
354 <  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
355 <
677 <  /* If triangles do not overlap: done */
678 <  if(!test)
679 <    return(FALSE);
680 <
681 <  /* if this is tree: recurse */
682 <  if(QT_IS_TREE(*qtptr))
683 <  {
684 <    qtSubdivide_tri(v0,v1,v2,a,b,c);
685 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t0,t1,t2,v0,a,c,func,arg);
686 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t0,t1,t2,a,v1,b,func,arg);
687 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t0,t1,t2,c,b,v2,func,arg);
688 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t0,t1,t2,b,c,a,func,arg);
342 >            for(optr = QT_SET_PTR(tptr),i = QT_SET_CNT(tptr); i > 0; i--)
343 >              {
344 >                id = QT_SET_NEXT_ELEM(optr);
345 >                if(!qtTri_from_id(id,r0,r1,r2))
346 >                  continue;
347 >                qt  = qtAdd_tri(qt,q0,q1,q2,r0,r1,r2,id,n);
348 >              }
349 >            if(tptr != tset)
350 >              free(tptr);
351 >            }
352 >            else
353 >                qt = qtaddelem(qt,id);
354 >        }
355 >      }
356    }
357 <  else
358 <    return(func(qtptr,arg));
357 >  return(qt);
358 > memerr:
359 >    error(SYSTEM,"qtAdd_tri():Unable to allocate memory");
360   }
361  
362  
363 < int
364 < qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
365 < QUADTREE *qtptr;
363 > QUADTREE
364 > qtRemove_tri(qt,id,q0,q1,q2,t0,t1,t2)
365 > QUADTREE qt;
366   int id;
367 + FVECT q0,q1,q2;
368   FVECT t0,t1,t2;
700 – FVECT v0,v1,v2;
369   {
702 –  
703 –  char test;
704 –  int i;
370    FVECT a,b,c;
706 –  OBJECT os[MAXSET+1];
371    
372    /* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
373 <  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
373 >  if(!stri_intersect(t0,t1,t2,q0,q1,q2))
374 >    return(qt);
375  
711 –  /* If triangles do not overlap: done */
712 –  if(!test)
713 –    return(FALSE);
714 –
376    /* if this is tree: recurse */
377 <  if(QT_IS_TREE(*qtptr))
377 >  if(QT_IS_TREE(qt))
378    {
379 <    qtSubdivide_tri(v0,v1,v2,a,b,c);
380 <    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c);
381 <    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b);
382 <    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2);
383 <    qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a);
379 >    qtSubdivide_tri(q0,q1,q2,a,b,c);
380 >    QT_NTH_CHILD(qt,0) = qtRemove_tri(QT_NTH_CHILD(qt,0),id,t0,t1,t2,q0,a,c);
381 >    QT_NTH_CHILD(qt,1) = qtRemove_tri(QT_NTH_CHILD(qt,1),id,t0,t1,t2,a,q1,b);
382 >    QT_NTH_CHILD(qt,2) = qtRemove_tri(QT_NTH_CHILD(qt,2),id,t0,t1,t2,c,b,q2);
383 >    QT_NTH_CHILD(qt,3) = qtRemove_tri(QT_NTH_CHILD(qt,3),id,t0,t1,t2,b,c,a);
384 >    return(qt);
385    }
386    else
387    {
388 <      if(QT_IS_EMPTY(*qtptr))
388 >      if(QT_IS_EMPTY(qt))
389        {
390   #ifdef DEBUG    
391            eputs("qtRemove_tri(): triangle not found\n");
# Line 732 | Line 394 | FVECT v0,v1,v2;
394        /* remove id from set */
395        else
396        {
397 <          qtgetset(os,*qtptr);
736 <          if(!inset(os,id))
397 >          if(!qtinset(qt,id))
398            {
399   #ifdef DEBUG          
400                eputs("qtRemove_tri(): tri not in set\n");
401   #endif
402            }
403            else
404 <          {
405 <            *qtptr = qtdelelem(*qtptr,id);
404 >            qt = qtdelelem(qt,id);
405 >      }
406 >  }
407 >  return(qt);
408 > }
409 >
410 >
411 > QUADTREE
412 > qtVisit_tri_interior(qt,q0,q1,q2,t0,t1,t2,n0,n1,n2,n,func,f,argptr)
413 >   QUADTREE qt;
414 >   FVECT q0,q1,q2;
415 >   FVECT t0,t1,t2;
416 >   FVECT n0,n1,n2;
417 >   int n;
418 >   int (*func)(),*f;
419 >   int *argptr;
420 > {
421 >    FVECT a,b,c;
422 >    
423 >    /* If qt Flag set, or qt vertices interior to t0t1t2-descend */
424 > tree_modified:
425 >
426 >    if(QT_IS_TREE(qt))
427 >    {  
428 >        if(QT_IS_FLAG(qt) ||  point_in_stri_n(n0,n1,n2,q0))
429 >        {
430 >            QT_SET_FLAG(qt);
431 >            qtSubdivide_tri(q0,q1,q2,a,b,c);
432 >            /* descend to children */
433 >            QT_NTH_CHILD(qt,0) = qtVisit_tri_interior(QT_NTH_CHILD(qt,0),
434 >                                  q0,a,c,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
435 >            QT_NTH_CHILD(qt,1) = qtVisit_tri_interior(QT_NTH_CHILD(qt,1),
436 >                                  a,q1,b,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
437 >            QT_NTH_CHILD(qt,2) = qtVisit_tri_interior(QT_NTH_CHILD(qt,2),
438 >                                  c,b,q2,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
439 >            QT_NTH_CHILD(qt,3) = qtVisit_tri_interior(QT_NTH_CHILD(qt,3),
440 >                                  b,c,a,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
441          }
442 +    }
443 +    else
444 +      if((!QT_IS_EMPTY(qt) && QT_LEAF_IS_FLAG(qt)) ||
445 +         point_in_stri_n(n0,n1,n2,q0))
446 +      {
447 +        func(&qt,f,argptr,q0,q1,q2,t0,t1,t2,n);
448 +        if(QT_FLAG_IS_MODIFIED(*f))
449 +       {
450 +         QT_SET_FLAG(qt);
451 +         goto tree_modified;
452 +       }
453 +        if(QT_IS_LEAF(qt))
454 +           QT_LEAF_SET_FLAG(qt);
455 +        else
456 +          if(QT_IS_TREE(qt))
457 +            QT_SET_FLAG(qt);
458        }
459 +    return(qt);
460 + }
461 +
462 +
463 +
464 + int
465 + move_to_nbri(b,db0,db1,db2,tptr)
466 + BCOORD b[3];
467 + BDIR db0,db1,db2;
468 + TINT *tptr;
469 + {
470 +  TINT t,dt;
471 +  BCOORD bc;
472 +  int nbr;
473 +  
474 +  nbr = -1;
475 +  *tptr = 0;
476 +  /* Advance to next node */
477 +  if(b[0]==0 && db0 < 0)
478 +    return(0);
479 +  if(b[1]==0 && db1 < 0)
480 +    return(1);
481 +  if(b[2]==0 && db2 < 0)
482 +    return(2);
483 +
484 +  if(db0 < 0)
485 +   {
486 +     bc = b[0]<<SHIFT_MAXBCOORD;
487 +     t = bc/-db0;
488 +     nbr = 0;
489 +   }
490 +  else
491 +    t = HUGET;
492 +  if(db1 < 0)
493 +  {
494 +     bc = b[1] <<SHIFT_MAXBCOORD;
495 +     dt = bc/-db1;
496 +    if( dt < t)
497 +    {
498 +      t = dt;
499 +      nbr = 1;
500 +    }
501    }
502 <  return(TRUE);
502 >  if(db2 < 0)
503 >  {
504 >     bc = b[2] << SHIFT_MAXBCOORD;
505 >     dt = bc/-db2;
506 >       if( dt < t)
507 >      {
508 >        t = dt;
509 >        nbr = 2;
510 >      }
511 >    }
512 >  *tptr = t;
513 >  return(nbr);
514   }
515 +
516 + QUADTREE
517 + qtVisit_tri_edges(qt,b,db0,db1,db2,db,wptr,nextptr,t,sign,sfactor,func,f,argptr)
518 +   QUADTREE qt;
519 +   BCOORD b[3];
520 +   BDIR db0,db1,db2,db[3][3];
521 +   int *wptr,*nextptr;
522 +   TINT t[3];
523 +   int sign,sfactor;
524 +   int (*func)();
525 +   int *f,*argptr;
526 + {
527 +    int i,found;
528 +    QUADTREE child;
529 +    int nbr,next,w;
530 +    TINT t_g,t_l,t_i,l;
531 +
532 +    if(QT_IS_TREE(qt))
533 +    {
534 +      /* Find the appropriate child and reset the coord */
535 +      i = baryi_child(b);
536 +
537 +      QT_SET_FLAG(qt);
538 +
539 +      for(;;)
540 +      {
541 +        w = *wptr;
542 +        child = QT_NTH_CHILD(qt,i);
543 +        if(i != 3)
544 +          QT_NTH_CHILD(qt,i) =
545 +            qtVisit_tri_edges(child,b,db0,db1,db2,db,wptr,nextptr,t,sign,
546 +                                sfactor+1,func,f,argptr);
547 +        else
548 +          /* If the center cell- must flip direction signs */
549 +          QT_NTH_CHILD(qt,i) =
550 +            qtVisit_tri_edges(child,b,-db0,-db1,-db2,db,wptr,nextptr,t,1-sign,
551 +                                  sfactor+1,func,f,argptr);
552 +
553 +        if(QT_FLAG_IS_DONE(*f))
554 +          return(qt);
555 +        if(*wptr != w)
556 +        {
557 +          w = *wptr;
558 +          db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
559 +          if(sign)
560 +            {  db0 *= -1;db1 *= -1; db2 *= -1;}
561 +        }
562 +        /* If in same block: traverse */
563 +        if(i==3)
564 +          next = *nextptr;
565 +        else
566 +          if(*nextptr == i)
567 +            next = 3;
568 +          else
569 +         {
570 +           /* reset the barycentric coordinates in the parents*/
571 +           baryi_parent(b,i);
572 +           /* Else pop up to parent and traverse from there */
573 +           return(qt);
574 +         }
575 +        baryi_from_child(b,i,next);
576 +        i = next;
577 +      }
578 +    }
579 +    else
580 +   {
581 +     func(&qt,f,argptr);
582 +     if(QT_FLAG_IS_DONE(*f))
583 +     {
584 +       if(!QT_IS_EMPTY(qt))
585 +         QT_LEAF_SET_FLAG(qt);
586 +       return(qt);
587 +     }
588 +    
589 +     if(!QT_IS_EMPTY(qt))
590 +       QT_LEAF_SET_FLAG(qt);
591 +     /* Advance to next node */
592 +     w = *wptr;
593 +     while(1)
594 +       {
595 +         nbr = move_to_nbri(b,db0,db1,db2,&t_i);
596 +        
597 +         t_g = t_i >> sfactor;
598 +                
599 +         if(t_g >= t[w])
600 +         {
601 +           if(w == 2)
602 +           {
603 +             QT_FLAG_SET_DONE(*f);
604 +             return(qt);
605 +           }
606 +           /* The edge fits in the cell- therefore the result of shifting
607 +              db by sfactor is guaranteed to be less than MAXBCOORD
608 +              */
609 +           /* Caution: (t[w]*db) must occur before divide by MAXBCOORD
610 +              since t[w] will always be < MAXBCOORD
611 +              */
612 +           l = t[w] << sfactor;
613 +           /* NOTE: Change divides to Shift and multiply by sign*/
614 +           b[0] += (l*db0)/MAXBCOORD;
615 +           b[1] += (l*db1)/MAXBCOORD;
616 +           b[2] += (l*db2)/MAXBCOORD;
617 +           w++;
618 +           db0 = db[w][0]; db1 = db[w][1]; db2 = db[w][2];
619 +           if(sign)
620 +             {  db0 *= -1;db1 *= -1; db2 *= -1;}
621 +         }
622 +         else
623 +         {
624 +           /* Caution: (t_i*db) must occur before divide by MAXBCOORD
625 +              since t_i will always be < MAXBCOORD*/
626 +           /* NOTE: Change divides to Shift and  by sign*/
627 +           b[0] += (t_i *db0) / MAXBCOORD;
628 +           b[1] += (t_i *db1) / MAXBCOORD;
629 +           b[2] += (t_i *db2) / MAXBCOORD;
630 +          
631 +           t[w] -= t_g;
632 +           *wptr = w;
633 +           *nextptr = nbr;
634 +           return(qt);
635 +         }
636 +       }
637 +   }    
638 + }
639 +
640 +
641 + QUADTREE
642 + qtRoot_visit_tri_edges(qt,q0,q1,q2,peq,tri,i_pt,wptr,nextptr,func,f,argptr)
643 +   QUADTREE qt;
644 +   FVECT q0,q1,q2;
645 +   FPEQ  peq;
646 +   FVECT tri[3],i_pt;
647 +   int *wptr,*nextptr;
648 +   int (*func)();
649 +   int *f,*argptr;
650 + {
651 +  int x,y,z,w,i,j,first;
652 +  QUADTREE child;
653 +  FVECT c,d,v[3];
654 +  double b[4][3],db[3][3],et[3],exit_pt;
655 +  BCOORD bi[3];
656 +  TINT t[3];
657 +  BDIR dbi[3][3];
658 +  
659 + first =0;
660 +  w = *wptr;
661 +  if(w==-1)
662 +  {
663 +    first = 1;
664 +    *wptr = w = 0;
665 +  }
666 +  /* Project the origin onto the root node plane */
667 +
668 +  t[0] = t[1] = t[2] = 0;
669 +  /* Find the intersection point of the origin */
670 +  /* map to 2d by dropping maximum magnitude component of normal */
671 +
672 +  x = FP_X(peq);
673 +  y = FP_Y(peq);
674 +  z = FP_Z(peq);
675 +  /* Calculate barycentric coordinates for current vertex */
676 +  if(!first)
677 +    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);  
678 +  else
679 +  /* Just starting: b[0] is the origin point: guaranteed to be valid b*/
680 +  {
681 +    intersect_vector_plane(tri[0],peq,&(et[0]),v[0]);
682 +    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);  
683 +    VCOPY(b[3],b[0]);
684 +  }
685 +
686 +  j = (w+1)%3;
687 +  intersect_vector_plane(tri[j],peq,&(et[j]),v[j]);
688 +  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);  
689 +  if(et[j] < 0.0)
690 +  {
691 +      VSUB(db[w],b[3],b[j]);
692 +      t[w] = HUGET;
693 +  }
694 +  else
695 + {
696 +   /* NOTE: for stability: do not increment with ipt- use full dir and
697 +      calculate t: but for wrap around case: could get same problem?
698 +      */
699 +   VSUB(db[w],b[j],b[3]);
700 +   /* Check if the point is out side of the triangle: if so t[w] =HUGET */
701 +   if((fabs(b[j][0])>(FTINY+1.0)) ||(fabs(b[j][1])>(FTINY+1.0)) ||
702 +      (fabs(b[j][2])>(FTINY+1.0))||(b[j][0] <-FTINY) ||
703 +      (b[j][1]<-FTINY) ||(b[j][2]<-FTINY))
704 +     t[w] = HUGET;
705 +   else
706 +   {
707 +       /* The next point is in the triangle- so db values will be in valid
708 +          range and t[w]= MAXT
709 +        */  
710 +       t[w] = MAXT;
711 +       if(j != 0)
712 +         for(;j < 3;j++)
713 +         {
714 +           i = (j+1)%3;
715 +           if(!first || i != 0)
716 +           {
717 +             intersect_vector_plane(tri[i],peq,&(et[i]),v[i]);
718 +             bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
719 +                    v[i][y],b[i]);
720 +           }
721 +           if(et[i] < 0.0)
722 +           {
723 +             VSUB(db[j],b[j],b[i]);
724 +             t[j] = HUGET;
725 +             break;
726 +           }
727 +           else
728 +           {
729 +             VSUB(db[j],b[i],b[j]);
730 +             if((fabs(b[j][0])>(FTINY+1.0))||(fabs(b[j][1])>(FTINY+1.0)) ||
731 +                (fabs(b[j][2])>(FTINY+1.0))||(b[i][0] <-FTINY) ||
732 +                (b[i][1]<-FTINY) ||(b[i][2]<-FTINY))
733 +               {
734 +                 t[j] = HUGET;
735 +                 break;
736 +               }
737 +             else
738 +               t[j] = MAXT;
739 +           }
740 +         }
741 +   }
742 + }
743 +  bary_dtol(b[3],db,bi,dbi,t,w);
744 +  
745 +  /* trace the ray starting with this node */
746 +  qt = qtVisit_tri_edges(qt,bi,dbi[w][0],dbi[w][1],dbi[w][2],
747 +                     dbi,wptr,nextptr,t,0,0,func,f,argptr);
748 +  if(!QT_FLAG_IS_DONE(*f))
749 +  {
750 +    b[3][0] = (double)bi[0]/(double)MAXBCOORD;
751 +    b[3][1] = (double)bi[1]/(double)MAXBCOORD;
752 +    b[3][2] = (double)bi[2]/(double)MAXBCOORD;
753 +    i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
754 +    i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
755 +    i_pt[z] = (-FP_N(peq)[x]*i_pt[x] - FP_N(peq)[y]*i_pt[y]-FP_D(peq))/FP_N(peq)[z];
756 +  }
757 +  return(qt);
758 +    
759 + }
760 +
761 +
762 + QUADTREE
763 + qtRoot_trace_ray(qt,q0,q1,q2,peq,orig,dir,nextptr,func,f,argptr)
764 +   QUADTREE qt;
765 +   FVECT q0,q1,q2;
766 +   FPEQ  peq;
767 +   FVECT orig,dir;
768 +   int *nextptr;
769 +   int (*func)();
770 +   int *f,*argptr;
771 + {
772 +  int x,y,z,nbr,w,i;
773 +  QUADTREE child;
774 +  FVECT v,i_pt;
775 +  double b[2][3],db[3],et[2],d,br[3];
776 +  BCOORD bi[3];
777 +  TINT t[3];
778 +  BDIR dbi[3][3];
779 +  
780 +  /* Project the origin onto the root node plane */
781 +  t[0] = t[1] = t[2] = 0;
782 +
783 +  VADD(v,orig,dir);
784 +  /* Find the intersection point of the origin */
785 +  /* map to 2d by dropping maximum magnitude component of normal */
786 +  x = FP_X(peq);
787 +  y = FP_Y(peq);
788 +  z = FP_Z(peq);
789 +
790 +  /* Calculate barycentric coordinates for current vertex */
791 +  intersect_vector_plane(orig,peq,&(et[0]),i_pt);
792 +  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[0]);  
793 +
794 +  intersect_vector_plane(v,peq,&(et[1]),i_pt);
795 +  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[1]);  
796 +  if(et[1] < 0.0)
797 +    VSUB(db,b[0],b[1]);
798 +  else
799 +   VSUB(db,b[1],b[0]);
800 +  t[0] = HUGET;
801 +  convert_dtol(b[0],bi);
802 +   if(et[1]<0.0 || (fabs(b[1][0])>(FTINY+1.0)) ||(fabs(b[1][1])>(FTINY+1.0)) ||
803 +      (fabs(b[1][2])>(FTINY+1.0))||(b[1][0] <-FTINY) ||
804 +      (b[1][1]<-FTINY) ||(b[1][2]<-FTINY))
805 +     {
806 +       max_index(db,&d);
807 +       for(i=0; i< 3; i++)
808 +         dbi[0][i] = (BDIR)(db[i]/d*MAXBDIR);
809 +     }
810 +   else
811 +       for(i=0; i< 3; i++)
812 +         dbi[0][i] = (BDIR)(db[i]*MAXBDIR);
813 +  w=0;
814 +  /* trace the ray starting with this node */
815 +  qt = qtVisit_tri_edges(qt,bi,dbi[0][0],dbi[0][1],dbi[0][2], dbi,&w,
816 +                         nextptr,t,0,0,func,f,argptr);
817 +  if(!QT_FLAG_IS_DONE(*f))
818 +  {
819 +    br[0] = (double)bi[0]/(double)MAXBCOORD;
820 +    br[1] = (double)bi[1]/(double)MAXBCOORD;
821 +    br[2] = (double)bi[2]/(double)MAXBCOORD;
822 +    orig[x] = br[0]*q0[x] + br[1]*q1[x] + br[2]*q2[x];
823 +    orig[y] = br[0]*q0[y] + br[1]*q1[y] + br[2]*q2[y];
824 +    orig[z]=(-FP_N(peq)[x]*orig[x] -
825 +             FP_N(peq)[y]*orig[y]-FP_D(peq))/FP_N(peq)[z];
826 +  }
827 +  return(qt);
828 +    
829 + }
830 +
831  
832  
833  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines