ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/hd/sm_qtree.c
Revision: 3.15
Committed: Wed Apr 23 00:52:34 2003 UTC (22 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.14: +1 -1 lines
Log Message: 
Added (void *) cast to realloc calls

File Contents

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