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

File Contents

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