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.2 by gwlarson, Thu Aug 20 16:47:21 1998 UTC vs.
Revision 3.6 by gwlarson, Wed Sep 16 18:16:29 1998 UTC

# Line 17 | Line 17 | static char SCCSid[] = "$SunId$ SGI";
17  
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= 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 + #endif
31  
27 –
32   QUADTREE
33   qtAlloc()                       /* allocate a quadtree */
34   {
# Line 41 | Line 45 | qtAlloc()                      /* allocate a quadtree */
45          if (QT_BLOCK(freet) >= QT_MAX_BLK)
46             return(EMPTY);
47          if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
48 <                   (unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
48 >                        QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
49             return(EMPTY);
50 +
51 +        quad_flag = (int4 *)realloc((char *)quad_flag,
52 +                        (QT_BLOCK(freet)+1)*(QT_BLOCK_SIZE/8));
53 +        if (quad_flag == NULL)
54 +                return(EMPTY);
55      }
56      treetop += 4;
57      return(freet);
58   }
59  
60  
61 + qtClearAllFlags()               /* clear all quadtree branch flags */
62 + {
63 +        if (!treetop) return;
64 +        bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*(QT_BLOCK_SIZE/8));
65 + }
66 +
67 +
68   qtFree(qt)                      /* free a quadtree */
69     register QUADTREE  qt;
70   {
# Line 70 | Line 86 | qtDone()                       /* free EVERYTHING */
86   {
87          register int    i;
88  
89 +        qtfreeleaves();
90          for (i = 0; i < QT_MAX_BLK; i++)
91          {
92 <            free((char *)quad_block[i],
93 <                  (unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE));
92 >            if (quad_block[i] == NULL)
93 >                break;
94 >            free((char *)quad_block[i]);
95              quad_block[i] = NULL;
96          }
97 +        if (i) free((char *)quad_flag);
98 +        quad_flag = NULL;
99          quad_free_list = EMPTY;
100          treetop = 0;
101   }
102  
103   QUADTREE
104 < qtCompress(qt)                  /* recursively combine nodes */
85 < register QUADTREE  qt;
86 < {
87 <    register int  i;
88 <    register QUADTREE  qres;
89 <
90 <    if (!QT_IS_TREE(qt))        /* not a tree */
91 <       return(qt);
92 <    qres = QT_NTH_CHILD(qt,0) = qtCompress(QT_NTH_CHILD(qt,0));
93 <    for (i = 1; i < 4; i++)
94 <       if((QT_NTH_CHILD(qt,i) = qtCompress(QT_NTH_CHILD(qt,i))) != qres)
95 <          qres = qt;
96 <    if(!QT_IS_TREE(qres))
97 <    {   /* all were identical leaves */
98 <        QT_NTH_CHILD(qt,0) = quad_free_list;
99 <        quad_free_list = qt;
100 <    }
101 <    return(qres);
102 < }
103 <
104 < QUADTREE
105 < qtLocate_leaf(qtptr,bcoord,type,which)
104 > *qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
105   QUADTREE *qtptr;
106   double bcoord[3];
107 < char *type,*which;
107 > FVECT t0,t1,t2;
108   {
109    int i;
110    QUADTREE *child;
111 +  FVECT a,b,c;
112  
113      if(QT_IS_TREE(*qtptr))
114      {  
115 <      i = bary2d_child(bcoord);
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 <      return(qtLocate_leaf(child,bcoord,type,which));
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);
135 >      return(qtptr);
136   }
137  
138  
124 –
125 –
139   QUADTREE
140 < qtRoot_point_locate(qtptr,v0,v1,v2,pt,type,which)
140 > *qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
141   QUADTREE *qtptr;
142   FVECT v0,v1,v2;
143   FVECT pt;
144 < char *type,*which;
144 > FVECT t0,t1,t2;
145   {
146 <    char d,w;
146 >    int d;
147      int i,x,y;
148      QUADTREE *child;
149 <    QUADTREE qt;
137 <    FVECT n,i_pt;
149 >    FVECT n,i_pt,a,b,c;
150      double pd,bcoord[3];
151  
152      /* Determine if point lies within pyramid (and therefore
# Line 144 | Line 156 | char *type,*which;
156         For each triangle edge: compare the
157         point against the plane formed by the edge and the view center
158      */
159 <    d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);  
159 >    d = point_in_stri(v0,v1,v2,pt);  
160  
161 +    
162      /* Not in this triangle */
163      if(!d)
164 <    {
152 <      if(which)
153 <        *which = 0;
154 <      return(EMPTY);
155 <    }
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
# Line 165 | Line 174 | char *type,*which;
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);
177 >      i = max_index(n,NULL);
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 = bary2d_child(bcoord);
183 >      i = bary_child(bcoord);
184        child = QT_NTH_CHILD_PTR(*qtptr,i);
185 <      return(qtLocate_leaf(child,bcoord,type,which));
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
179 –       if(!QT_IS_EMPTY(*qtptr))
180 –       {  
181 –           /* map GT_VERTEX,GT_EDGE,GT_FACE GT_INTERIOR of pyramid to
182 –              spherical triangle primitives
183 –              */
184 –         if(type)
185 –           *type = d;
186 –         if(which)
187 –           *which = w;
188 –         return(*qtptr);
189 –       }
190 –    return(EMPTY);
191 – }
192 –
193 – int
194 – qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
195 – QUADTREE *qtptr;
196 – FVECT v0,v1,v2;
197 – FVECT pt;
198 – char *type,*which;
199 – {
200 –    QUADTREE qt;
201 –    OBJECT os[MAXSET+1],*optr;
202 –    char d,w;
203 –    int i,id;
204 –    FVECT p0,p1,p2;
205 –
206 –    qt = qtRoot_point_locate(qtptr,v0,v1,v2,pt,type,which);
207 –
208 –    if(QT_IS_EMPTY(qt))
209 –       return(EMPTY);
210 –    
211 –    /* Get the set */
212 –    qtgetset(os,qt);
213 –    for (i = QT_SET_CNT(os),optr = QT_SET_PTR(os); i > 0; i--)
207      {
208 <        /* Find the triangle that pt falls in (NOTE:FOR now return first 1) */
209 <        id = QT_SET_NEXT_ELEM(optr);
210 <        qtTri_verts_from_id(id,p0,p1,p2);
211 <        d = test_single_point_against_spherical_tri(p0,p1,p2,pt,&w);  
212 <        if(d)
220 <        {
221 <          if(type)
222 <            *type = d;
223 <          if(which)
224 <            *which = w;
225 <          return(id);
208 >        if(t0)
209 >        {
210 >            VCOPY(t0,v0);
211 >            VCOPY(t1,v1);
212 >            VCOPY(t2,v2);
213          }
214 +        return(qtptr);
215      }
228 –    return(EMPTY);
216   }
217  
218 +
219   QUADTREE
220   qtSubdivide(qtptr)
221   QUADTREE *qtptr;
# Line 279 | Line 267 | FVECT a,b,c;
267   int i;
268   FVECT r0,r1,r2;
269   {
270 +
271    switch(i){
272    case 0:  
273    VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
# Line 304 | Line 293 | into the new child cells: it is assumed that "v1,v2,v3
293   */
294  
295   int
296 < qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
296 > qtRoot_add_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
297   QUADTREE *qtptr;
298 < int id;
298 > FVECT q0,q1,q2;
299   FVECT t0,t1,t2;
300 < FVECT v0,v1,v2;
300 > int id;
301   int n;
302   {
303 +  int test;
304 +  int found;
305 +
306 +  test = stri_intersect(q0,q1,q2,t0,t1,t2);
307 +  if(!test)
308 +    return(FALSE);
309    
310 <  char test;
311 <  int i,index;
310 >  found = qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n);
311 >  
312 >  return(found);
313 > }
314 >
315 > int
316 > qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
317 > QUADTREE *qtptr;
318 > FVECT q0,q1,q2;
319 > FVECT t0,t1,t2;
320 > int id;
321 > int n;
322 > {
323 >  int i,index,test,found;
324    FVECT a,b,c;
325 <  OBJECT os[MAXSET+1],*optr;
325 >  OBJECT os[QT_MAXSET+1],*optr;
326    QUADTREE qt;
320 –  int found;
327    FVECT r0,r1,r2;
328  
323 –  /* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
324 –  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
325 –
326 –  /* If triangles do not overlap: done */
327 –  if(!test)
328 –    return(FALSE);
329    found = 0;
330 –
330    /* if this is tree: recurse */
331    if(QT_IS_TREE(*qtptr))
332    {
333      n++;
334 <    qtSubdivide_tri(v0,v1,v2,a,b,c);
335 <    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c,n);
336 <    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b,n);
337 <    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2,n);
338 <    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a,n);
339 <
340 < #if 0
341 <    if(!found)
342 <    {
343 < #ifdef TEST_DRIVER    
344 <      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
345 < #else
346 <      eputs("qtAdd_tri():Found in parent but not children\n");
348 < #endif
349 <    }
350 < #endif
334 >    qtSubdivide_tri(q0,q1,q2,a,b,c);
335 >    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);
347    }
348    else
349    {
350        /* If this leave node emptry- create a new set */
351        if(QT_IS_EMPTY(*qtptr))
352 <      {
357 <          *qtptr = qtaddelem(*qtptr,id);
358 < #if 0
359 <          {
360 <              int k;
361 <              qtgetset(os,*qtptr);
362 <              printf("\n%d:\n",os[0]);
363 <              for(k=1; k <= os[0];k++)
364 <                 printf("%d ",os[k]);
365 <              printf("\n");
366 <          }
367 < #endif
368 <          /*
369 <          os[0] = 0;
370 <          insertelem(os,id);
371 <          qt = fullnode(os);
372 <          *qtptr = qt;
373 <          */
374 <      }
352 >        *qtptr = qtaddelem(*qtptr,id);
353        else
354        {
377 –          qtgetset(os,*qtptr);
355            /* If the set is too large: subdivide */
356 <          if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
357 <          {
358 <            *qtptr = qtaddelem(*qtptr,id);
359 < #if 0
383 <            {
384 <              int k;
385 <              qtgetset(os,*qtptr);
386 <              printf("\n%d:\n",os[0]);
387 <              for(k=1; k <= os[0];k++)
388 <                 printf("%d ",os[k]);
389 <              printf("\n");
390 <          }
391 < #endif      
392 <            /*
393 <             insertelem(os,id);
394 <             *qtptr = fullnode(os);
395 <             */
396 <          }
356 >        optr = qtqueryset(*qtptr);
357 >
358 >        if(QT_SET_CNT(optr) < QT_SET_THRESHOLD)
359 >          *qtptr = qtaddelem(*qtptr,id);
360            else
361            {
362              if (n < QT_MAX_LEVELS)
# Line 401 | Line 364 | int n;
364                   /* If set size exceeds threshold: subdivide cell and
365                      reinsert set tris into cell
366                      */
367 <                 n++;
368 <                 qtfreeleaf(*qtptr);
369 <                 qtSubdivide(qtptr);
370 <                 found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
371 < #if 0
372 <                 if(!found)
373 <                 {
374 < #ifdef TEST_DRIVER    
375 <      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
376 < #else
377 <                   eputs("qtAdd_tri():Found in parent but not children\n");
378 < #endif
416 <                 }
417 < #endif
418 <                 for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
419 <                 {
420 <                   id = QT_SET_NEXT_ELEM(optr);
421 <                   qtTri_verts_from_id(id,r0,r1,r2);
422 <                   found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,n);
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,r0,r1,r2,NULL,NULL,NULL,NULL,NULL,NULL);
378 >                  found=qtAdd_tri(qtptr,q0,q1,q2,r0,r1,r2,id,n);
379   #ifdef DEBUG
380 <                 if(!found)
380 >                  if(!found)
381                      eputs("qtAdd_tri():Reinsert-in parent but not children\n");
382   #endif
383 <               }
384 <             }
383 >                }
384 >            }
385              else
386 <              if(QT_SET_CNT(os) < QT_MAX_SET)
431 <                {
386 >              if(QT_SET_CNT(optr) < QT_MAXSET)
387                    *qtptr = qtaddelem(*qtptr,id);
433 – #if 0
434 –                  {
435 –              int k;
436 –              qtgetset(os,*qtptr);
437 –              printf("\n%d:\n",os[0]);
438 –              for(k=1; k <= os[0];k++)
439 –                 printf("%d ",os[k]);
440 –              printf("\n");
441 –          }
442 – #endif            
443 –                  /*
444 –                    insertelem(os,id);
445 –                    *qtptr = fullnode(os);
446 –                  */
447 –                }
388                else
389                  {
390   #ifdef DEBUG
# Line 465 | Line 405 | QUADTREE *qtptr;
405   FVECT t0,t1,t2;
406   FVECT v0,v1,v2;
407   int (*func)();
408 < char *arg;
408 > int *arg;
409   {
410 <  char test;
410 >  int test;
411    FVECT a,b,c;
412  
413    /* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
414 <  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
414 >  test = stri_intersect(t0,t1,t2,v0,v1,v2);
415  
416    /* 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 <    qtSubdivide_tri(v0,v1,v2,a,b,c);
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    }
489 –  else
490 –    return(func(qtptr,arg));
432   }
433  
493 –
434   int
435   qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
436   QUADTREE *qtptr;
# Line 499 | Line 439 | FVECT t0,t1,t2;
439   FVECT v0,v1,v2;
440   {
441    
442 <  char test;
442 >  int test;
443    int i;
444    FVECT a,b,c;
445 <  OBJECT os[MAXSET+1];
445 >  OBJECT os[QT_MAXSET+1];
446    
447    /* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
448 <  test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
448 >  test = stri_intersect(t0,t1,t2,v0,v1,v2);
449  
450    /* If triangles do not overlap: done */
451    if(!test)
# Line 531 | Line 471 | FVECT v0,v1,v2;
471        /* remove id from set */
472        else
473        {
474 <          qtgetset(os,*qtptr);
535 <          if(!inset(os,id))
474 >          if(!qtinset(*qtptr,id))
475            {
476   #ifdef DEBUG          
477                eputs("qtRemove_tri(): tri not in set\n");
# Line 541 | Line 480 | FVECT v0,v1,v2;
480            else
481            {
482              *qtptr = qtdelelem(*qtptr,id);
483 < #if 0
483 >        }
484 >      }
485 >  }
486 >  return(TRUE);
487 > }
488 >
489 >
490 > int
491 > move_to_nbr(b,db0,db1,db2,tptr)
492 > double b[3],db0,db1,db2;
493 > double *tptr;
494 > {
495 >  double t,dt;
496 >  int nbr;
497 >
498 >  nbr = -1;
499 >  /* Advance to next node */
500 >  if(!ZERO(db0) && db0 < 0.0)
501 >   {
502 >     t = -b[0]/db0;
503 >     nbr = 0;
504 >   }
505 >  else
506 >    t = FHUGE;
507 >  if(!ZERO(db1) && db1 < 0.0 )
508 >  {
509 >    dt = -b[1]/db1;
510 >    if( dt < t)
511 >    {
512 >      t = dt;
513 >      nbr = 1;
514 >    }
515 >  }
516 >  if(!ZERO(db2) && db2 < 0.0 )
517 >    {
518 >      dt = -b[2]/db2;
519 >      if( dt < t)
520 >      {
521 >        t = dt;
522 >        nbr = 2;
523 >      }
524 >    }
525 >  *tptr = t;
526 >  return(nbr);
527 > }
528 >
529 > int
530 > qtTrace_ray(qtptr,b,db0,db1,db2,orig,dir,func,arg1,arg2)
531 >   QUADTREE *qtptr;
532 >   double b[3],db0,db1,db2;
533 >   FVECT orig,dir;
534 >   int (*func)();
535 >   int *arg1,arg2;
536 > {
537 >
538 >    int i,found;
539 >    QUADTREE *child;
540 >    int nbr,next;
541 >    double t;
542 > #ifdef DEBUG_TEST_DRIVER
543 >    
544 >    FVECT a1,b1,c1;
545 >    int Pick_parent = Pick_cnt-1;
546 >    qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
547 >                           Pick_v2[Pick_parent],a1,b1,c1);
548 >
549 > #endif
550 >    if(QT_IS_TREE(*qtptr))
551 >    {
552 >        /* Find the appropriate child and reset the coord */
553 >        i = bary_child(b);
554 >
555 >        QT_SET_FLAG(*qtptr);
556 >
557 >        for(;;)
558 >        {
559 >           child = QT_NTH_CHILD_PTR(*qtptr,i);
560 >
561 >           if(i != 3)
562 >              nbr = qtTrace_ray(child,b,db0,db1,db2,orig,dir,func,arg1,arg2);
563 >           else
564 >               /* If the center cell- must flip direction signs */
565 >              nbr =qtTrace_ray(child,b,-db0,-db1,-db2,orig,dir,func,arg1,arg2);
566 >           if(nbr == QT_DONE)
567 >              return(nbr);
568 >
569 >           /* If in same block: traverse */
570 >           if(i==3)
571 >              next = nbr;
572 >             else
573 >                if(nbr == i)
574 >                   next = 3;
575 >             else
576 >               {
577 >                 /* reset the barycentric coordinates in the parents*/
578 >                 bary_parent(b,i);
579 >                 /* Else pop up to parent and traverse from there */
580 >                 return(nbr);
581 >               }
582 >             bary_from_child(b,i,next);
583 >             i = next;
584 >         }
585 >    }
586 >    else
587 >   {
588 > #ifdef DEBUG_TEST_DRIVER
589 >           qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
590 >                           Pick_v2[Pick_parent],a1,b1,c1,i,
591 >                           Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
592 >                           Pick_v2[Pick_cnt]);
593 >           Pick_cnt++;
594 > #endif
595 >
596 >       if(func(qtptr,orig,dir,arg1,arg2) == QT_DONE)
597 >          return(QT_DONE);
598 >
599 >       /* Advance to next node */
600 >       /* NOTE: Optimize: should only have to check 1/2 */
601 >       nbr = move_to_nbr(b,db0,db1,db2,&t);
602 >
603 >       if(nbr != -1)
604 >       {
605 >         b[0] += t * db0;
606 >         b[1] += t * db1;
607 >         b[2] += t * db2;
608 >       }
609 >       return(nbr);
610 >   }
611 >    
612 > }
613 >
614 > int
615 > qtRoot_trace_ray(qtptr,q0,q1,q2,orig,dir,func,arg1,arg2)
616 >   QUADTREE *qtptr;
617 >   FVECT q0,q1,q2;
618 >   FVECT orig,dir;
619 >   int (*func)();
620 >   int *arg1,arg2;
621 > {
622 >  int i,x,y,nbr;
623 >  QUADTREE *child;
624 >  FVECT n,c,i_pt,d;
625 >  double pd,b[3],db[3],t;
626 >    /* Determine if point lies within pyramid (and therefore
627 >       inside a spherical quadtree cell):GT_INTERIOR, on one of the
628 >       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
629 >       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
630 >       For each triangle edge: compare the
631 >       point against the plane formed by the edge and the view center
632 >    */
633 >  i = point_in_stri(q0,q1,q2,orig);  
634 >    
635 >  /* Not in this triangle */
636 >  if(!i)
637 >     return(INVALID);
638 >  /* Project the origin onto the root node plane */
639 >
640 >  /* Find the intersection point of the origin */
641 >  tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
642 >  intersect_vector_plane(orig,n,pd,NULL,i_pt);
643 >  /* project the dir as well */
644 >  VADD(c,orig,dir);
645 >  intersect_vector_plane(c,n,pd,&t,c);
646 >
647 >    /* map to 2d by dropping maximum magnitude component of normal */
648 >  i = max_index(n,NULL);
649 >  x = (i+1)%3;
650 >  y = (i+2)%3;
651 >  /* Calculate barycentric coordinates of orig */
652 >  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b);
653 >  /* Calculate barycentric coordinates of dir */
654 >  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db);
655 >  if(t < 0.0)
656 >     VSUB(db,b,db);
657 >  else
658 >     VSUB(db,db,b);
659 >
660 >
661 > #ifdef DEBUG_TEST_DRIVER
662 >    VCOPY(Pick_v0[Pick_cnt],q0);
663 >    VCOPY(Pick_v1[Pick_cnt],q1);
664 >    VCOPY(Pick_v2[Pick_cnt],q2);
665 >    Pick_cnt++;
666 > #endif
667 >
668 >      /* trace the ray starting with this node */
669 >    nbr = qtTrace_ray(qtptr,b,db[0],db[1],db[2],orig,dir,func,arg1,arg2);
670 >    return(nbr);
671 >    
672 > }
673 >
674 > qtVisit_tri_interior(qtptr,q0,q1,q2,t0,t1,t2,n,func,arg1,arg2,arg3)
675 >   QUADTREE *qtptr;
676 >   FVECT q0,q1,q2;
677 >   FVECT t0,t1,t2;
678 >   int n;
679 >   int (*func)();
680 >   int *arg1,arg2,*arg3;
681 > {
682 >    int i,found,test;
683 >    QUADTREE *child;
684 >    FVECT c0,c1,c2,a,b,c;
685 >    OBJECT os[QT_MAXSET+1],*optr;
686 >    int w;
687 >    
688 >    /* If qt Flag set, or qt vertices interior to t0t1t2-descend */
689 > tree_modified:
690 >
691 >    if(QT_IS_TREE(*qtptr))
692 >    {  
693 >        if(QT_IS_FLAG(*qtptr) ||  point_in_stri(t0,t1,t2,q0))
694 >        {
695 >            QT_SET_FLAG(*qtptr);
696 >            qtSubdivide_tri(q0,q1,q2,a,b,c);
697 >            /* descend to children */
698 >            for(i=0;i < 4; i++)
699              {
700 <              int k;
701 <              if(!QT_IS_EMPTY(*qtptr))
702 <              {qtgetset(os,*qtptr);
703 <              printf("\n%d:\n",os[0]);
704 <              for(k=1; k <= os[0];k++)
705 <                 printf("%d ",os[k]);
706 <              printf("\n");
707 <           }
700 >                child = QT_NTH_CHILD_PTR(*qtptr,i);
701 >                qtNth_child_tri(q0,q1,q2,a,b,c,i,c0,c1,c2);
702 >                qtVisit_tri_interior(child,c0,c1,c2,t0,t1,t2,n+1,
703 >                                     func,arg1,arg2,arg3);
704 >            }
705 >        }
706 >    }
707 >    else
708 >    {
709 >      /* NOTE THIS IN TRI TEST Could be replaced by a flag */
710 >      if(!QT_IS_EMPTY(*qtptr))
711 >      {
712 >         if(qtinset(*qtptr,arg2))
713 >           if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2,arg3)==QT_MODIFIED)
714 >             goto tree_modified;
715 >           else
716 >             return;
717 >      }
718 >      if(point_in_stri(t0,t1,t2,q0) )
719 >          if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2,arg3)==QT_MODIFIED)
720 >            goto tree_modified;
721 >    }
722 > }
723  
724 <          }
724 >
725 >
726 >
727 >
728 >
729 > /* NOTE: SINCE DIR could be unit: then we could use integer math */
730 > int
731 > qtVisit_tri_edges(qtptr,b,db0,db1,db2,
732 >                   db,wptr,t,sign,sfactor,func,arg1,arg2,arg3)
733 >   QUADTREE *qtptr;
734 >   double b[3],db0,db1,db2,db[3][3];
735 >   int *wptr;
736 >   double t[3];
737 >   int sign;
738 >   double sfactor;
739 >   int (*func)();
740 >   int *arg1,arg2,*arg3;
741 > {
742 >    int i,found;
743 >    QUADTREE *child;
744 >    int nbr,next,w;
745 >    double t_l,t_g;
746 > #ifdef DEBUG_TEST_DRIVER
747 >    FVECT a1,b1,c1;
748 >    int Pick_parent = Pick_cnt-1;
749 >    qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
750 >                           Pick_v2[Pick_parent],a1,b1,c1);
751   #endif
752 +    if(QT_IS_TREE(*qtptr))
753 +    {
754 +        /* Find the appropriate child and reset the coord */
755 +        i = bary_child(b);
756 +
757 +        QT_SET_FLAG(*qtptr);
758 +
759 +        for(;;)
760 +        {
761 +          w = *wptr;
762 +           child = QT_NTH_CHILD_PTR(*qtptr,i);
763 +
764 +           if(i != 3)
765 +               nbr = qtVisit_tri_edges(child,b,db0,db1,db2,
766 +                                        db,wptr,t,sign,
767 +                                        sfactor*2.0,func,arg1,arg2,arg3);
768 +           else
769 +             /* If the center cell- must flip direction signs */
770 +             nbr = qtVisit_tri_edges(child,b,-db0,-db1,-db2,
771 +                                     db,wptr,t,1-sign,
772 +                                     sfactor*2.0,func,arg1,arg2,arg3);
773 +
774 +           if(nbr == QT_DONE)
775 +             return(nbr);
776 +           if(*wptr != w)
777 +           {
778 +             w = *wptr;
779 +             db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
780 +             if(sign)
781 +               {  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
782 +           }
783 +           /* If in same block: traverse */
784 +           if(i==3)
785 +              next = nbr;
786 +             else
787 +                if(nbr == i)
788 +                   next = 3;
789 +             else
790 +               {
791 +                 /* reset the barycentric coordinates in the parents*/
792 +                 bary_parent(b,i);
793 +                 /* Else pop up to parent and traverse from there */
794 +                 return(nbr);
795 +               }
796 +           bary_from_child(b,i,next);
797 +           i = next;
798          }
799 +    }
800 +    else
801 +   {
802 + #ifdef DEBUG_TEST_DRIVER
803 +           qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
804 +                           Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
805 +                           Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
806 +           Pick_cnt++;
807 + #endif
808 +
809 +       if(func(qtptr,arg1,arg2,arg3) == QT_DONE)
810 +          return(QT_DONE);
811 +
812 +       /* Advance to next node */
813 +       w = *wptr;
814 +       while(1)
815 +       {
816 +         nbr = move_to_nbr(b,db0,db1,db2,&t_l);
817 +
818 +         t_g = t_l/sfactor;
819 + #ifdef DEBUG
820 +         if(t[w] <= 0.0)
821 +           eputs("qtVisit_tri_edges():negative t\n");
822 + #endif
823 +         if(t_g >= t[w])
824 +         {
825 +           if(w == 2)
826 +             return(QT_DONE);
827 +
828 +           b[0] += (t[w])*sfactor*db0;
829 +           b[1] += (t[w])*sfactor*db1;
830 +           b[2] += (t[w])*sfactor*db2;
831 +           w++;
832 +           db0 = db[w][0];
833 +           db1 = db[w][1];
834 +           db2 = db[w][2];
835 +           if(sign)
836 +          {  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
837 +         }
838 +       else
839 +         if(nbr != INVALID)
840 +         {
841 +           b[0] += t_l * db0;
842 +           b[1] += t_l * db1;
843 +           b[2] += t_l * db2;
844 +
845 +           t[w] -= t_g;
846 +           *wptr = w;
847 +           return(nbr);
848 +         }
849 +         else
850 +           return(INVALID);
851 +       }
852 +   }
853 +    
854 + }
855 +
856 +
857 + int
858 + qtRoot_visit_tri_edges(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2,arg3)
859 +   QUADTREE *qtptr;
860 +   FVECT q0,q1,q2;
861 +   FVECT tri[3],i_pt;
862 +   int *wptr;
863 +   int (*func)();
864 +   int *arg1,arg2,*arg3;
865 + {
866 +  int x,y,z,nbr,w,i,j;
867 +  QUADTREE *child;
868 +  FVECT n,c,d,v[3];
869 +  double pd,b[4][3],db[3][3],et[3],t[3],exit_pt;
870 +    
871 +  w = *wptr;
872 +
873 +  /* Project the origin onto the root node plane */
874 +
875 +  /* Find the intersection point of the origin */
876 +  tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
877 +  /* map to 2d by dropping maximum magnitude component of normal */
878 +  z = max_index(n,NULL);
879 +  x = (z+1)%3;
880 +  y = (z+2)%3;
881 +  /* Calculate barycentric coordinates for current vertex */
882 +  if(w != -1)
883 +  {
884 +    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);  
885 +    intersect_vector_plane(tri[w],n,pd,&(et[w]),v[w]);
886 +    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[w][x],v[w][y],b[w]);  
887 +  }
888 +  else
889 +  /* Just starting: b[0] is the origin point: guaranteed to be valid b*/
890 +  {
891 +    w = 0;
892 +    intersect_vector_plane(tri[0],n,pd,&(et[0]),v[0]);
893 +    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);  
894 +    VCOPY(b[3],b[0]);
895 +  }
896 +
897 +
898 +  j = (w+1)%3;
899 +  intersect_vector_plane(tri[j],n,pd,&(et[j]),v[j]);
900 +  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);  
901 +  if(et[j] < 0.0)
902 +  {
903 +      VSUB(db[w],b[3],b[j]);
904 +      t[w] = FHUGE;
905 +  }
906 +  else
907 + {
908 +   /* NOTE: for stability: do not increment with ipt- use full dir and
909 +      calculate t: but for wrap around case: could get same problem?
910 +      */
911 +   VSUB(db[w],b[j],b[3]);
912 +   t[w] = 1.0;
913 +   move_to_nbr(b[3],db[w][0],db[w][1],db[w][2],&exit_pt);
914 +   if(exit_pt >= 1.0)
915 +   {
916 +     for(;j < 3;j++)
917 +      {
918 +          i = (j+1)%3;
919 +          if(i!= w)
920 +            {
921 +              intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
922 +              bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
923 +                     v[i][y],b[i]);
924 +            }
925 +          if(et[i] < 0.0)
926 +             {
927 +                 VSUB(db[j],b[j],b[i]);
928 +                 t[j] = FHUGE;
929 +                 break;
930 +             }
931 +          else
932 +             {
933 +                 VSUB(db[j],b[i],b[j]);
934 +                 t[j] = 1.0;
935 +             }
936 +          move_to_nbr(b[j],db[j][0],db[j][1],db[j][2],&exit_pt);
937 +          if(exit_pt < 1.0)
938 +             break;
939        }
940 +   }
941 + }
942 +  *wptr = w;
943 +  /* trace the ray starting with this node */
944 +    nbr = qtVisit_tri_edges(qtptr,b[3],db[w][0],db[w][1],db[w][2],
945 +                             db,wptr,t,0,1.0,func,arg1,arg2,arg3);
946 +    if(nbr != INVALID && nbr != QT_DONE)
947 +    {
948 +      i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
949 +      i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
950 +      i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
951 +    }
952 +    return(nbr);
953 +    
954 + }
955 +
956 + int
957 + move_to_nbri(b,db0,db1,db2,tptr)
958 + BCOORD b[3];
959 + BDIR db0,db1,db2;
960 + TINT *tptr;
961 + {
962 +  TINT t,dt;
963 +  BCOORD bc;
964 +  int nbr;
965 +  
966 +  nbr = -1;
967 +  /* Advance to next node */
968 +  if(db0 < 0)
969 +   {
970 +     bc = MAXBCOORD*b[0];
971 +     t = bc/-db0;
972 +     nbr = 0;
973 +   }
974 +  else
975 +    t = HUGET;
976 +  if(db1 < 0)
977 +  {
978 +     bc = MAXBCOORD*b[1];
979 +     dt = bc/-db1;
980 +    if( dt < t)
981 +    {
982 +      t = dt;
983 +      nbr = 1;
984 +    }
985    }
986 <  return(TRUE);
986 >  if(db2 < 0)
987 >  {
988 >     bc = MAXBCOORD*b[2];
989 >     dt = bc/-db2;
990 >       if( dt < t)
991 >      {
992 >        t = dt;
993 >        nbr = 2;
994 >      }
995 >    }
996 >  *tptr = t;
997 >  return(nbr);
998   }
999 + /* NOTE: SINCE DIR could be unit: then we could use integer math */
1000 + int
1001 + qtVisit_tri_edgesi(qtptr,b,db0,db1,db2,
1002 +                   db,wptr,t,sign,sfactor,func,arg1,arg2,arg3)
1003 +   QUADTREE *qtptr;
1004 +   BCOORD b[3];
1005 +   BDIR db0,db1,db2,db[3][3];
1006 +   int *wptr;
1007 +   TINT t[3];
1008 +   int sign;
1009 +   int sfactor;
1010 +   int (*func)();
1011 +   int *arg1,arg2,*arg3;
1012 + {
1013 +    int i,found;
1014 +    QUADTREE *child;
1015 +    int nbr,next,w;
1016 +    TINT t_g,t_l,t_i;
1017 + #ifdef DEBUG_TEST_DRIVER
1018 +    FVECT a1,b1,c1;
1019 +    int Pick_parent = Pick_cnt-1;
1020 +    qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
1021 +                           Pick_v2[Pick_parent],a1,b1,c1);
1022 + #endif
1023 +    if(QT_IS_TREE(*qtptr))
1024 +    {
1025 +        /* Find the appropriate child and reset the coord */
1026 +        i = baryi_child(b);
1027 +
1028 +        QT_SET_FLAG(*qtptr);
1029 +
1030 +        for(;;)
1031 +        {
1032 +          w = *wptr;
1033 +           child = QT_NTH_CHILD_PTR(*qtptr,i);
1034 +
1035 +           if(i != 3)
1036 +               nbr = qtVisit_tri_edgesi(child,b,db0,db1,db2,
1037 +                                        db,wptr,t,sign,
1038 +                                        sfactor+1,func,arg1,arg2,arg3);
1039 +           else
1040 +             /* If the center cell- must flip direction signs */
1041 +             nbr = qtVisit_tri_edgesi(child,b,-db0,-db1,-db2,
1042 +                                     db,wptr,t,1-sign,
1043 +                                     sfactor+1,func,arg1,arg2,arg3);
1044 +
1045 +           if(nbr == QT_DONE)
1046 +             return(nbr);
1047 +           if(*wptr != w)
1048 +           {
1049 +             w = *wptr;
1050 +             db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
1051 +             if(sign)
1052 +               {  db0 *= -1;db1 *= -1; db2 *= -1;}
1053 +           }
1054 +           /* If in same block: traverse */
1055 +           if(i==3)
1056 +              next = nbr;
1057 +             else
1058 +                if(nbr == i)
1059 +                   next = 3;
1060 +             else
1061 +               {
1062 +                 /* reset the barycentric coordinates in the parents*/
1063 +                 baryi_parent(b,i);
1064 +                 /* Else pop up to parent and traverse from there */
1065 +                 return(nbr);
1066 +               }
1067 +           baryi_from_child(b,i,next);
1068 +           i = next;
1069 +        }
1070 +    }
1071 +    else
1072 +   {
1073 + #ifdef DEBUG_TEST_DRIVER
1074 +           qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
1075 +                           Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
1076 +                           Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
1077 +           Pick_cnt++;
1078 + #endif
1079 +
1080 +       if(func(qtptr,arg1,arg2,arg3) == QT_DONE)
1081 +          return(QT_DONE);
1082 +
1083 +       /* Advance to next node */
1084 +       w = *wptr;
1085 +       while(1)
1086 +       {
1087 +           nbr = move_to_nbri(b,db0,db1,db2,&t_i);
1088 +
1089 +         t_g = t_i >> sfactor;
1090 +                
1091 +         if(t_g >= t[w])
1092 +         {
1093 +           if(w == 2)
1094 +             return(QT_DONE);
1095 +
1096 +           /* The edge fits in the cell- therefore the result of shifting db by
1097 +              sfactor is guaranteed to be less than MAXBCOORD
1098 +            */
1099 +           /* Caution: (t[w]*db) must occur before divide by MAXBCOORD
1100 +              since t[w] will always be < MAXBCOORD
1101 +           */
1102 +           b[0] = b[0] + (t[w]*db0*(1 << sfactor))/MAXBCOORD;
1103 +           b[1] = b[1] + (t[w]*db1*(1 << sfactor))/MAXBCOORD;
1104 +           b[2] = b[2] + (t[w]*db2*(1 << sfactor))/MAXBCOORD;
1105 +           w++;
1106 +           db0 = db[w][0]; db1 = db[w][1]; db2 = db[w][2];
1107 +           if(sign)
1108 +           {  db0 *= -1;db1 *= -1; db2 *= -1;}
1109 +         }
1110 +       else
1111 +         if(nbr != INVALID)
1112 +         {
1113 +           /* Caution: (t_i*db) must occur before divide by MAXBCOORD
1114 +              since t_i will always be < MAXBCOORD
1115 +           */
1116 +           b[0] = b[0] + (t_i *db0) / MAXBCOORD;
1117 +           b[1] = b[1] + (t_i *db1) / MAXBCOORD;
1118 +           b[2] = b[2] + (t_i *db2) / MAXBCOORD;
1119 +
1120 +           t[w] -= t_g;
1121 +           *wptr = w;
1122 +           return(nbr);
1123 +         }
1124 +         else
1125 +           return(INVALID);
1126 +       }
1127 +   }    
1128 + }
1129 +
1130 +
1131 + int
1132 + qtRoot_visit_tri_edgesi(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2,arg3)
1133 +   QUADTREE *qtptr;
1134 +   FVECT q0,q1,q2;
1135 +   FVECT tri[3],i_pt;
1136 +   int *wptr;
1137 +   int (*func)();
1138 +   int *arg1,arg2,*arg3;
1139 + {
1140 +  int x,y,z,nbr,w,i,j;
1141 +  QUADTREE *child;
1142 +  FVECT n,c,d,v[3];
1143 +  double pd,b[4][3],db[3][3],et[3],exit_pt;
1144 +  BCOORD bi[3];
1145 +  TINT t[3];
1146 +  BDIR dbi[3][3];
1147 +  w = *wptr;
1148 +
1149 +  /* Project the origin onto the root node plane */
1150 +
1151 +  t[0] = t[1] = t[2] = 0;
1152 +  /* Find the intersection point of the origin */
1153 +  tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
1154 +  /* map to 2d by dropping maximum magnitude component of normal */
1155 +  z = max_index(n,NULL);
1156 +  x = (z+1)%3;
1157 +  y = (z+2)%3;
1158 +  /* Calculate barycentric coordinates for current vertex */
1159 +  if(w != -1)
1160 +  {
1161 +    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);  
1162 +    intersect_vector_plane(tri[w],n,pd,&(et[w]),v[w]);
1163 +    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[w][x],v[w][y],b[w]);  
1164 +  }
1165 +  else
1166 +  /* Just starting: b[0] is the origin point: guaranteed to be valid b*/
1167 +  {
1168 +    w = 0;
1169 +    intersect_vector_plane(tri[0],n,pd,&(et[0]),v[0]);
1170 +    bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);  
1171 +    VCOPY(b[3],b[0]);
1172 +  }
1173 +
1174 +
1175 +  j = (w+1)%3;
1176 +  intersect_vector_plane(tri[j],n,pd,&(et[j]),v[j]);
1177 +  bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);  
1178 +  if(et[j] < 0.0)
1179 +  {
1180 +      VSUB(db[w],b[3],b[j]);
1181 +      t[w] = HUGET;
1182 +  }
1183 +  else
1184 + {
1185 +   /* NOTE: for stability: do not increment with ipt- use full dir and
1186 +      calculate t: but for wrap around case: could get same problem?
1187 +      */
1188 +   VSUB(db[w],b[j],b[3]);
1189 +   /* Check if the point is out side of the triangle: if so t[w] =HUGET */
1190 +   if((fabs(b[j][0])>(FTINY+1.0)) ||(fabs(b[j][1])>(FTINY+1.0)) ||
1191 +      (fabs(b[j][2])>(FTINY+1.0)))
1192 +      t[w] = HUGET;
1193 +   else
1194 +   {
1195 +       /* The next point is in the triangle- so db values will be in valid
1196 +          range and t[w]= MAXT
1197 +        */  
1198 +       t[w] = MAXT;
1199 +       if(j != 0)
1200 +          for(;j < 3;j++)
1201 +             {
1202 +                 i = (j+1)%3;
1203 +                 if(i!= w)
1204 +                 {
1205 +                     intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
1206 +                     bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
1207 +                            v[i][y],b[i]);
1208 +                 }
1209 +                 if(et[i] < 0.0)
1210 +                 {
1211 +                     VSUB(db[j],b[j],b[i]);
1212 +                     t[j] = HUGET;
1213 +                     break;
1214 +                 }
1215 +                 else
1216 +                 {
1217 +                     VSUB(db[j],b[i],b[j]);
1218 +                     if((fabs(b[j][0])>(FTINY+1.0)) ||
1219 +                        (fabs(b[j][1])>(FTINY+1.0)) ||
1220 +                        (fabs(b[j][2])>(FTINY+1.0)))
1221 +                        {
1222 +                            t[j] = HUGET;
1223 +                            break;
1224 +                        }
1225 +                     else
1226 +                        t[j] = MAXT;
1227 +                 }
1228 +             }
1229 +   }
1230 + }
1231 +  *wptr = w;
1232 +  bary_dtol(b[3],db,bi,dbi,t);
1233 +
1234 +  /* trace the ray starting with this node */
1235 +    nbr = qtVisit_tri_edgesi(qtptr,bi,dbi[w][0],dbi[w][1],dbi[w][2],
1236 +                             dbi,wptr,t,0,0,func,arg1,arg2,arg3);
1237 +    if(nbr != INVALID && nbr != QT_DONE)
1238 +    {
1239 +        b[3][0] = (double)bi[0]/(double)MAXBCOORD;
1240 +        b[3][1] = (double)bi[1]/(double)MAXBCOORD;
1241 +        b[3][2] = (double)bi[2]/(double)MAXBCOORD;
1242 +        i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
1243 +        i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
1244 +        i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
1245 +    }
1246 +    return(nbr);
1247 +    
1248 + }
1249 +
1250 +
1251 +
1252 +
1253 +
1254  
1255  
1256  

Diff Legend

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