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, 9 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

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: sm_qtree.c,v 3.16 2003/06/20 00:25:49 greg Exp $";
3 #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 <string.h>
14
15 #include "standard.h"
16 #include "sm_flag.h"
17 #include "sm_geom.h"
18 #include "object.h"
19 #include "sm_qtree.h"
20
21 QUADTREE *quad_block[QT_MAX_BLK]; /* our quadtree */
22 static QUADTREE quad_free_list = EMPTY; /* freed octree nodes */
23 static QUADTREE treetop = 0; /* next free node */
24 int32 *quad_flag= NULL; /* quadtree flags */
25
26
27 qtremovelast(qt,id)
28 QUADTREE qt;
29 OBJECT id;
30 {
31 #ifdef DEBUG
32 if(qtqueryset(qt)[(*(qtqueryset(qt)))] != id)
33 error(CONSISTENCY,"not removed\n");
34 #endif
35 ((*(qtqueryset(qt)))--);
36 }
37 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 QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
54 error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
55
56 /* Realloc the per/node flags */
57 quad_flag = (int32 *)realloc((void *)quad_flag,
58 (QT_BLOCK(freet)+1)*((QT_BLOCK_SIZE+7)>>3));
59 if (quad_flag == NULL)
60 error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
61 }
62 treetop += 4;
63 return(freet);
64 }
65
66
67 qtClearAllFlags() /* clear all quadtree branch flags */
68 {
69 if (!treetop)
70 return;
71
72 /* Clear the node flags*/
73 memset((char *)quad_flag, '\0',
74 (QT_BLOCK(treetop-4)+1)*((QT_BLOCK_SIZE+7)>>3));
75 /* Clear set flags */
76 qtclearsetflags();
77 }
78
79 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
100 qtfreeleaves();
101 for (i = 0; i < QT_MAX_BLK; i++)
102 {
103 if (quad_block[i] == NULL)
104 break;
105 free((void *)quad_block[i]);
106 quad_block[i] = NULL;
107 }
108 /* Free the flags */
109 if (i) free((void *)quad_flag);
110 quad_flag = NULL;
111 quad_free_list = EMPTY;
112 treetop = 0;
113 }
114
115 /*
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 QUADTREE
273 qtLocate(qt,bcoord)
274 QUADTREE qt;
275 BCOORD bcoord[3];
276 {
277 int i;
278
279 if(QT_IS_TREE(qt))
280 {
281 i = bary_child(bcoord);
282
283 return(qtLocate(QT_NTH_CHILD(qt,i),bcoord));
284 }
285 else
286 return(qt);
287 }
288
289 int
290 move_to_nbr(b,db0,db1,db2,tptr)
291 BCOORD b[3];
292 BDIR db0,db1,db2;
293 double *tptr;
294 {
295 double t,dt;
296 BCOORD bc;
297 int nbr;
298
299 nbr = -1;
300 *tptr = 0.0;
301 /* Advance to next node */
302 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 if(db0 < 0)
309 {
310 t = ((double)b[0])/-db0;
311 nbr = 0;
312 }
313 else
314 t = MAXFLOAT;
315 if(db1 < 0)
316 {
317 dt = ((double)b[1])/-db1;
318 if( dt < t)
319 {
320 t = dt;
321 nbr = 1;
322 }
323 }
324 if(db2 < 0)
325 {
326 dt = ((double)b[2])/-db2;
327 if( dt < t)
328 {
329 t = dt;
330 nbr = 2;
331 }
332 }
333 *tptr = t;
334 return(nbr);
335 }
336
337 qtTrace_ray(qt,b,db0,db1,db2,nextptr,sign,sfactor,func,f)
338 QUADTREE qt;
339 BCOORD b[3];
340 BDIR db0,db1,db2;
341 int *nextptr;
342 int sign,sfactor;
343 FUNC func;
344 int *f;
345 {
346 int i,found;
347 QUADTREE child;
348 int nbr,next,w;
349 double t;
350
351 if(QT_IS_TREE(qt))
352 {
353 /* Find the appropriate child and reset the coord */
354 i = bary_child(b);
355
356 for(;;)
357 {
358 child = QT_NTH_CHILD(qt,i);
359 if(i != 3)
360 qtTrace_ray(child,b,db0,db1,db2,nextptr,sign,sfactor+1,func,f);
361 else
362 /* If the center cell- must flip direction signs */
363 qtTrace_ray(child,b,-db0,-db1,-db2,nextptr,1-sign,sfactor+1,func,f);
364
365 if(QT_FLAG_IS_DONE(*f))
366 return;
367 /* 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 bary_parent(b,i);
377 /* Else pop up to parent and traverse from there */
378 return(qt);
379 }
380 bary_from_child(b,i,next);
381 i = next;
382 }
383 }
384 else
385 {
386 #ifdef TEST_DRIVER
387 if(Pick_cnt < 500)
388 Pick_q[Pick_cnt++]=qt;
389 #endif;
390 F_FUNC(func)(qt,F_ARGS(func),f);
391 if(QT_FLAG_IS_DONE(*f))
392 return;
393 /* Advance to next node */
394 nbr = move_to_nbr(b,db0,db1,db2,&t);
395
396 if(nbr==-1)
397 {
398 QT_FLAG_SET_DONE(*f);
399 return;
400 }
401 b[0] += (BCOORD)(t *db0);
402 b[1] += (BCOORD)(t *db1);
403 b[2] += (BCOORD)(t *db2);
404 *nextptr = nbr;
405 return;
406
407 }
408 }
409
410 #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 QUADTREE
418 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 {
438 goto Lfunc_call;
439 }
440 /* 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 {
445 goto Lfunc_call;
446 }
447 if(t1[0]== a || t1[1] == b || t1[2] == c)
448 {
449 goto Lfunc_call;
450 }
451 if(t2[0]== a || t2[1] == b || t2[2] == c)
452 {
453 goto Lfunc_call;
454 }
455 /* 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 {
465 /* 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 }
469 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
523 /* 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
541 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
574 /* 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 {
581 goto Lfunc_call;
582 }
583 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 {
590 goto Lfunc_call;
591 }
592 if(t1[0]== a || t1[1] == b || t1[2] == c)
593 {
594 goto Lfunc_call;
595 }
596 if(t2[0]== a || t2[1] == b || t2[2] == c)
597 {
598 goto Lfunc_call;
599 }
600 /* 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 {
613 db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
614 TEST_INT(testl,testh,db,q1[1],b,bl,bh)
615 }
616 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 test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
624 ((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
668 /* 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 {
672 db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
673 TEST_INT(testl,testh,db,q1[1],b,bl,bh)
674 }
675 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 }
680 if(!test_t2t0 && (e0 & 4))
681 {
682 db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
683 TEST_INT(testl,testh,db,q1[1],b,bl,bh)
684 }
685
686 return(qt);
687 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 }
805
806
807 QUADTREE
808 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 int n;
819 {
820 BCOORD a,b,c;
821 BCOORD va[3],vb[3],vc[3];
822 unsigned int sa,sb,sc;
823
824 /* 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
834 SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
835
836 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 else
899 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 broke into two sets of routines. Also traverses only to n levels.
910 */
911
912 qtVisit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
913 s0,s1,s2,sq0,sq1,sq2,f,n)
914 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 int n;
923 {
924 BCOORD a,b,c;
925 BCOORD va[3],vb[3],vc[3];
926 unsigned int sa,sb,sc;
927
928 if(n == 0)
929 return;
930 /* If a tree: trivial reject and recurse on potential children */
931 if(QT_IS_TREE(qt))
932 {
933 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 vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
953 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 t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
964 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 q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
975 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 vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
989 return;
990 }
991
992 if(sa)
993 qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
994 t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
995 if(sb)
996 qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
997 t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
998 if(sc)
999 qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
1000 t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
1001 qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
1002 t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
1003 }
1004 /* Leaf node: Do definitive test */
1005 else
1006 qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1007 scale,s0,s1,s2,sq0,sq1,sq2,f,n);
1008
1009 }
1010
1011
1012 qtVisit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
1013 s0,s1,s2,sq0,sq1,sq2,f,n)
1014 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 int n;
1023 {
1024 BCOORD a,b,c;
1025 BCOORD va[3],vb[3],vc[3];
1026 unsigned int sa,sb,sc;
1027
1028 if(n==0)
1029 return;
1030 /* 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
1041 SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
1042
1043 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 vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
1052 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 t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
1063 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 q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
1074 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 vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
1087 return;
1088 }
1089 if(sa!=7)
1090 qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
1091 t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
1092 if(sb!=7)
1093 qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
1094 t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
1095 if(sc!=7)
1096 qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
1097 t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
1098
1099 qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
1100 t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
1101 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 scale,s0,s1,s2,sq0,sq1,sq2,f,n);
1107 }
1108
1109
1110
1111 qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1112 scale, s0,s1,s2,sq0,sq1,sq2,f,n)
1113 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 int n;
1121 {
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
1128 /* 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
1232 f.func(f.argptr,root,qt,n);
1233 if(!QT_IS_EMPTY(qt))
1234 QT_LEAF_SET_FLAG(qt);
1235
1236 }
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 scale, s0,s1,s2,sq0,sq1,sq2,f,n)
1244 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 int n;
1252 {
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 f.func(f.argptr,root,qt,n);
1371 if(!QT_IS_EMPTY(qt))
1372 QT_LEAF_SET_FLAG(qt);
1373 }
1374
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
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565