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.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= 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 48 | qtAlloc()                      /* allocate a quadtree */
48          if (QT_BLOCK(freet) >= QT_MAX_BLK)
49             return(EMPTY);
50          if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
51 <                   (unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
52 <           return(EMPTY);
51 >                        QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
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+7)>>3));
57 >        if (quad_flag == NULL)
58 >           error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
59      }
60      treetop += 4;
61      return(freet);
62   }
63  
64  
65 + qtClearAllFlags()               /* clear all quadtree branch flags */
66 + {
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 +
76   qtFree(qt)                      /* free a quadtree */
77     register QUADTREE  qt;
78   {
# Line 69 | 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 <            free((char *)quad_block[i],
101 <                  (unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE));
100 >            if (quad_block[i] == NULL)
101 >                break;
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  
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   {
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)
106 – QUADTREE *qtptr;
107 – double bcoord[3];
108 – char *type,*which;
109 – {
117    int i;
111 –  QUADTREE *child;
118  
119 <    if(QT_IS_TREE(*qtptr))
120 <    {  
121 <      i = bary2d_child(bcoord);
122 <      child = QT_NTH_CHILD_PTR(*qtptr,i);
123 <      return(qtLocate_leaf(child,bcoord,type,which));
124 <    }
125 <    else
126 <      return(*qtptr);
119 >  if(QT_IS_TREE(qt))
120 >  {  
121 >    i = baryi_child(bcoord);
122 >
123 >    return(qtLocate_leaf(QT_NTH_CHILD(qt,i),bcoord));
124 >  }
125 >  else
126 >    return(qt);
127   }
128  
129 <
130 <
131 <
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(qtptr,v0,v1,v2,pt,type,which)
135 < QUADTREE *qtptr;
136 < FVECT v0,v1,v2;
134 > qtRoot_point_locate(qt,q0,q1,q2,peq,pt)
135 > QUADTREE qt;
136 > FVECT q0,q1,q2;
137 > FPEQ peq;
138   FVECT pt;
131 – char *type,*which;
139   {
133 –    char d,w;
140      int i,x,y;
141 <    QUADTREE *child;
142 <    QUADTREE qt;
143 <    FVECT n,i_pt;
138 <    double pd,bcoord[3];
141 >    FVECT i_pt;
142 >    double bcoord[3];
143 >    BCOORD bcoordi[3];
144  
145 <    /* Determine if point lies within pyramid (and therefore
141 <       inside a spherical quadtree cell):GT_INTERIOR, on one of the
142 <       pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
143 <       or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
144 <       For each triangle edge: compare the
145 <       point against the plane formed by the edge and the view center
146 <    */
147 <    d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);  
148 <
149 <    /* Not in this triangle */
150 <    if(!d)
151 <    {
152 <      if(which)
153 <        *which = 0;
154 <      return(EMPTY);
155 <    }
156 <
157 <    /* Will return lowest level triangle containing point: It the
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      {  
152        /* Find the intersection point */
153 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
166 <      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;
170 <      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);
164 <      child = QT_NTH_CHILD_PTR(*qtptr,i);
165 <      return(qtLocate_leaf(child,bcoord,type,which));
163 >      i = baryi_child(bcoordi);
164 >
165 >      return(qtLocate_leaf(QT_NTH_CHILD(qt,i),bcoordi));
166      }
167      else
168 <       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);
168 >      return(qt);
169   }
170  
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;
171  
206 –    qt = qtRoot_point_locate(qtptr,v0,v1,v2,pt,type,which);
172  
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--)
214 –    {
215 –        /* Find the triangle that pt falls in (NOTE:FOR now return first 1) */
216 –        id = QT_SET_NEXT_ELEM(optr);
217 –        qtTri_verts_from_id(id,p0,p1,p2);
218 –        d = test_single_point_against_spherical_tri(p0,p1,p2,pt,&w);  
219 –        if(d)
220 –        {
221 –          if(type)
222 –            *type = d;
223 –          if(which)
224 –            *which = w;
225 –          return(id);
226 –        }
227 –    }
228 –    return(EMPTY);
229 – }
173  
231 – QUADTREE
232 – qtSubdivide(qtptr)
233 – QUADTREE *qtptr;
234 – {
235 –  QUADTREE node;
236 –  node = qtAlloc();
237 –  QT_CLEAR_CHILDREN(node);
238 –  *qtptr = node;
239 –  return(node);
240 – }
241 –
242 –
243 – QUADTREE
244 – qtSubdivide_nth_child(qt,n)
245 – QUADTREE qt;
246 – int n;
247 – {
248 –  QUADTREE node;
249 –
250 –  node = qtSubdivide(&(QT_NTH_CHILD(qt,n)));
251 –  
252 –  return(node);
253 – }
254 –
174   /* for triangle v0-v1-v2- returns a,b,c: children are:
175  
176          v2                     0: v0,a,c
# Line 264 | Line 183 | int n;
183          a
184   */
185  
267 – qtSubdivide_tri(v0,v1,v2,a,b,c)
268 – FVECT v0,v1,v2;
269 – FVECT a,b,c;
270 – {
271 –    EDGE_MIDPOINT_VEC3(a,v0,v1);
272 –    EDGE_MIDPOINT_VEC3(b,v1,v2);
273 –    EDGE_MIDPOINT_VEC3(c,v2,v0);
274 – }
186  
187   qtNth_child_tri(v0,v1,v2,a,b,c,i,r0,r1,r2)
188   FVECT v0,v1,v2;
# Line 279 | 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 303 | 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 < qtAdd_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;
252 > int id,n;
253 > {
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 >  return(qt);
258 > }
259 >
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 > int id,n;
266 > {
267 >
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 >
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 > int id;
280   int n;
281   {
314 –  
315 –  char test;
316 –  int i,index;
282    FVECT a,b,c;
283 <  OBJECT os[MAXSET+1],*optr;
319 <  QUADTREE qt;
320 <  int found;
283 >  OBJECT tset[QT_MAXSET+1],*optr,*tptr;
284    FVECT r0,r1,r2;
285 +  int i;
286  
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 –
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);
336 <    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c,n);
337 <    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b,n);
338 <    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2,n);
339 <    found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a,n);
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 <    }
350 < #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 <      {
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 <      }
307 >      if(QT_IS_EMPTY(qt))
308 >        qt = qtaddelem(qt,id);
309        else
310        {
377 –          qtgetset(os,*qtptr);
311            /* If the set is too large: subdivide */
312 <          if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
313 <          {
314 <            *qtptr = qtaddelem(*qtptr,id);
315 < #if 0
316 <            {
317 <              int k;
318 <              qtgetset(os,*qtptr);
319 <              printf("\n%d:\n",os[0]);
320 <              for(k=1; k <= os[0];k++)
321 <                 printf("%d ",os[k]);
322 <              printf("\n");
323 <          }
324 < #endif      
325 <            /*
326 <             insertelem(os,id);
327 <             *qtptr = fullnode(os);
328 <             */
329 <          }
330 <          else
398 <          {
399 <            if (n < QT_MAX_LEVELS)
400 <            {
401 <                 /* If set size exceeds threshold: subdivide cell and
402 <                    reinsert set tris into cell
403 <                    */
404 <                 n++;
405 <                 qtfreeleaf(*qtptr);
406 <                 qtSubdivide(qtptr);
407 <                 found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
408 < #if 0
409 <                 if(!found)
410 <                 {
411 < #ifdef TEST_DRIVER    
412 <      HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
413 < #else
414 <                   eputs("qtAdd_tri():Found in parent but not children\n");
415 < #endif
416 <                 }
417 < #endif
418 <                 for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
419 <                 {
420 <                   id = QT_SET_NEXT_ELEM(optr);
421 <                   qtTri_verts_from_id(id,r0,r1,r2);
422 <                   found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,n);
423 < #ifdef DEBUG
424 <                 if(!found)
425 <                    eputs("qtAdd_tri():Reinsert-in parent but not children\n");
426 < #endif
427 <               }
428 <             }
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);
433 < #if 0
434 <                  {
435 <              int k;
436 <              qtgetset(os,*qtptr);
437 <              printf("\n%d:\n",os[0]);
438 <              for(k=1; k <= os[0];k++)
439 <                 printf("%d ",os[k]);
440 <              printf("\n");
441 <          }
442 < #endif            
443 <                  /*
444 <                    insertelem(os,id);
445 <                    *qtptr = fullnode(os);
446 <                  */
447 <                }
448 <              else
449 <                {
450 < #ifdef DEBUG
451 <                    eputs("qtAdd_tri():two many levels\n");
452 < #endif
453 <                    return(FALSE);
454 <                }
455 <          }
456 <      }
457 <  }
458 <  return(TRUE);
459 < }
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 <
476 <  /* If triangles do not overlap: done */
477 <  if(!test)
478 <    return(FALSE);
479 <
480 <  /* if this is tree: recurse */
481 <  if(QT_IS_TREE(*qtptr))
482 <  {
483 <    qtSubdivide_tri(v0,v1,v2,a,b,c);
484 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t0,t1,t2,v0,a,c,func,arg);
485 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t0,t1,t2,a,v1,b,func,arg);
486 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t0,t1,t2,c,b,v2,func,arg);
487 <    qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t0,t1,t2,b,c,a,func,arg);
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;
499 – FVECT v0,v1,v2;
369   {
501 –  
502 –  char test;
503 –  int i;
370    FVECT a,b,c;
505 –  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  
510 –  /* If triangles do not overlap: done */
511 –  if(!test)
512 –    return(FALSE);
513 –
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 531 | Line 394 | FVECT v0,v1,v2;
394        /* remove id from set */
395        else
396        {
397 <          qtgetset(os,*qtptr);
535 <          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);
406 < #if 0
407 <            {
408 <              int k;
547 <              if(!QT_IS_EMPTY(*qtptr))
548 <              {qtgetset(os,*qtptr);
549 <              printf("\n%d:\n",os[0]);
550 <              for(k=1; k <= os[0];k++)
551 <                 printf("%d ",os[k]);
552 <              printf("\n");
553 <           }
404 >            qt = qtdelelem(qt,id);
405 >      }
406 >  }
407 >  return(qt);
408 > }
409  
410 <          }
411 < #endif
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  
834  

Diff Legend

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