ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm_qtree.c
Revision: 3.17
Committed: Mon Jun 30 14:59:12 2003 UTC (20 years, 10 months ago) by schorsch
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R6P1, rad3R6
Changes since 3.16: +5 -2 lines
Log Message: 
Replaced most outdated BSD function calls with their posix equivalents, and cleaned up a few other platform dependencies.

File Contents

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