ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm_qtree.c
Revision: 3.10
Committed: Tue Jan 5 16:52:39 1999 UTC (25 years, 3 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.9: +0 -0 lines
Log Message: 
fixed logic in smTest to handle nearby and coincident points, added base tris to rendering, made list of new triangles to speed up rendering

File Contents

# User Rev Content
1 gwlarson 3.1 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2    
3     #ifndef lint
4     static char SCCSid[] = "$SunId$ SGI";
5     #endif
6    
7     /*
8     * sm_qtree.c: adapted from octree.c from radiance code
9     */
10     /*
11     * octree.c - routines dealing with octrees and cubes.
12     *
13     * 7/28/85
14     */
15    
16     #include "standard.h"
17 gwlarson 3.7 #include "sm_flag.h"
18 gwlarson 3.1 #include "sm_geom.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 gwlarson 3.4 int4 *quad_flag= NULL;
25 gwlarson 3.1
26 gwlarson 3.4 #ifdef TEST_DRIVER
27     extern FVECT Pick_v0[500],Pick_v1[500],Pick_v2[500];
28     extern int Pick_cnt,Pick_tri,Pick_samp;
29     extern FVECT Pick_point[500];
30 gwlarson 3.8 extern int Pick_q[500];
31 gwlarson 3.7
32 gwlarson 3.4 #endif
33 gwlarson 3.7 int Incnt=0;
34 gwlarson 3.4
35 gwlarson 3.1 QUADTREE
36     qtAlloc() /* allocate a quadtree */
37     {
38     register QUADTREE freet;
39    
40     if ((freet = quad_free_list) != EMPTY)
41     {
42     quad_free_list = QT_NTH_CHILD(freet, 0);
43     return(freet);
44     }
45     freet = treetop;
46     if (QT_BLOCK_INDEX(freet) == 0)
47     {
48     if (QT_BLOCK(freet) >= QT_MAX_BLK)
49     return(EMPTY);
50     if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
51 gwlarson 3.3 QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
52 gwlarson 3.7 error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
53 gwlarson 3.4
54 gwlarson 3.7 /* Realloc the per/node flags */
55 gwlarson 3.3 quad_flag = (int4 *)realloc((char *)quad_flag,
56 gwlarson 3.7 (QT_BLOCK(freet)+1)*((QT_BLOCK_SIZE+7)>>3));
57 gwlarson 3.3 if (quad_flag == NULL)
58 gwlarson 3.7 error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
59 gwlarson 3.1 }
60     treetop += 4;
61     return(freet);
62     }
63    
64    
65 gwlarson 3.3 qtClearAllFlags() /* clear all quadtree branch flags */
66     {
67 gwlarson 3.7 if (!treetop)
68     return;
69    
70     /* Clear the node flags*/
71     bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*((QT_BLOCK_SIZE+7)>>3));
72     /* Clear set flags */
73     qtclearsetflags();
74 gwlarson 3.3 }
75    
76 gwlarson 3.1 qtFree(qt) /* free a quadtree */
77     register QUADTREE qt;
78     {
79     register int i;
80    
81     if (!QT_IS_TREE(qt))
82     {
83     qtfreeleaf(qt);
84     return;
85     }
86     for (i = 0; i < 4; i++)
87     qtFree(QT_NTH_CHILD(qt, i));
88     QT_NTH_CHILD(qt, 0) = quad_free_list;
89     quad_free_list = qt;
90     }
91    
92    
93     qtDone() /* free EVERYTHING */
94     {
95     register int i;
96 gwlarson 3.7
97 gwlarson 3.4 qtfreeleaves();
98 gwlarson 3.1 for (i = 0; i < QT_MAX_BLK; i++)
99     {
100 gwlarson 3.3 if (quad_block[i] == NULL)
101     break;
102     free((char *)quad_block[i]);
103 gwlarson 3.1 quad_block[i] = NULL;
104     }
105 gwlarson 3.7 /* Free the flags */
106 gwlarson 3.3 if (i) free((char *)quad_flag);
107     quad_flag = NULL;
108 gwlarson 3.1 quad_free_list = EMPTY;
109     treetop = 0;
110     }
111    
112     QUADTREE
113 gwlarson 3.9 qtLocate(qt,bcoord)
114 gwlarson 3.7 QUADTREE qt;
115     BCOORD bcoord[3];
116 gwlarson 3.2 {
117     int i;
118    
119 gwlarson 3.7 if(QT_IS_TREE(qt))
120     {
121 gwlarson 3.9 i = bary_child(bcoord);
122 gwlarson 3.4
123 gwlarson 3.9 return(qtLocate(QT_NTH_CHILD(qt,i),bcoord));
124 gwlarson 3.7 }
125     else
126     return(qt);
127 gwlarson 3.2 }
128    
129 gwlarson 3.6 int
130 gwlarson 3.9 move_to_nbr(b,db0,db1,db2,tptr)
131 gwlarson 3.6 BCOORD b[3];
132     BDIR db0,db1,db2;
133 gwlarson 3.9 double *tptr;
134 gwlarson 3.6 {
135 gwlarson 3.9 double t,dt;
136 gwlarson 3.6 BCOORD bc;
137     int nbr;
138    
139     nbr = -1;
140 gwlarson 3.9 *tptr = 0.0;
141 gwlarson 3.6 /* Advance to next node */
142 gwlarson 3.7 if(b[0]==0 && db0 < 0)
143     return(0);
144     if(b[1]==0 && db1 < 0)
145     return(1);
146     if(b[2]==0 && db2 < 0)
147     return(2);
148 gwlarson 3.6 if(db0 < 0)
149     {
150 gwlarson 3.9 t = ((double)b[0])/-db0;
151 gwlarson 3.6 nbr = 0;
152     }
153     else
154 gwlarson 3.9 t = MAXFLOAT;
155 gwlarson 3.6 if(db1 < 0)
156     {
157 gwlarson 3.9 dt = ((double)b[1])/-db1;
158 gwlarson 3.6 if( dt < t)
159     {
160     t = dt;
161     nbr = 1;
162     }
163     }
164     if(db2 < 0)
165     {
166 gwlarson 3.9 dt = ((double)b[2])/-db2;
167 gwlarson 3.6 if( dt < t)
168     {
169     t = dt;
170     nbr = 2;
171     }
172     }
173     *tptr = t;
174     return(nbr);
175     }
176 gwlarson 3.7
177 gwlarson 3.9 qtTrace_ray(qt,b,db0,db1,db2,nextptr,sign,sfactor,func,f)
178 gwlarson 3.7 QUADTREE qt;
179 gwlarson 3.6 BCOORD b[3];
180 gwlarson 3.9 BDIR db0,db1,db2;
181     int *nextptr;
182 gwlarson 3.7 int sign,sfactor;
183 gwlarson 3.9 FUNC func;
184     int *f;
185 gwlarson 3.4 {
186     int i,found;
187 gwlarson 3.7 QUADTREE child;
188 gwlarson 3.4 int nbr,next,w;
189 gwlarson 3.9 double t;
190 gwlarson 3.7
191     if(QT_IS_TREE(qt))
192 gwlarson 3.4 {
193 gwlarson 3.7 /* Find the appropriate child and reset the coord */
194 gwlarson 3.9 i = bary_child(b);
195 gwlarson 3.4
196 gwlarson 3.7 for(;;)
197     {
198     child = QT_NTH_CHILD(qt,i);
199     if(i != 3)
200 gwlarson 3.9 qtTrace_ray(child,b,db0,db1,db2,nextptr,sign,sfactor+1,func,f);
201 gwlarson 3.7 else
202     /* If the center cell- must flip direction signs */
203 gwlarson 3.9 qtTrace_ray(child,b,-db0,-db1,-db2,nextptr,1-sign,sfactor+1,func,f);
204 gwlarson 3.7
205     if(QT_FLAG_IS_DONE(*f))
206 gwlarson 3.9 return;
207 gwlarson 3.7 /* If in same block: traverse */
208     if(i==3)
209     next = *nextptr;
210     else
211     if(*nextptr == i)
212     next = 3;
213     else
214     {
215     /* reset the barycentric coordinates in the parents*/
216 gwlarson 3.9 bary_parent(b,i);
217 gwlarson 3.7 /* Else pop up to parent and traverse from there */
218     return(qt);
219     }
220 gwlarson 3.9 bary_from_child(b,i,next);
221 gwlarson 3.7 i = next;
222     }
223 gwlarson 3.4 }
224     else
225     {
226 gwlarson 3.8 #ifdef TEST_DRIVER
227     if(Pick_cnt < 500)
228     Pick_q[Pick_cnt++]=qt;
229     #endif;
230 gwlarson 3.9 F_FUNC(func)(qt,F_ARGS(func),f);
231 gwlarson 3.7 if(QT_FLAG_IS_DONE(*f))
232 gwlarson 3.9 return;
233     /* Advance to next node */
234     nbr = move_to_nbr(b,db0,db1,db2,&t);
235    
236     if(nbr==-1)
237 gwlarson 3.7 {
238 gwlarson 3.9 QT_FLAG_SET_DONE(*f);
239     return;
240 gwlarson 3.7 }
241 gwlarson 3.9 b[0] += (BCOORD)(t *db0);
242     b[1] += (BCOORD)(t *db1);
243     b[2] += (BCOORD)(t *db2);
244     *nextptr = nbr;
245     return;
246 gwlarson 3.7
247 gwlarson 3.9 }
248     }
249 gwlarson 3.4
250    
251 gwlarson 3.9
252    
253    
254    
255    
256    
257     #define TEST_INT(tl,th,d,q0,q1,h,l) \
258     tl=d>q0;th=d<q1; if(tl&&th) goto Lfunc_call; \
259     if(tl) if(l) goto Lfunc_call; else h=TRUE; \
260     if(th) if(h) goto Lfunc_call; else l = TRUE; \
261     if(th) if(h) goto Lfunc_call; else l = TRUE;
262    
263     /* Leaf node: Do definitive test */
264 gwlarson 3.7 QUADTREE
265 gwlarson 3.9 qtLeaf_insert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
266     scale, s0,s1,s2,sq0,sq1,sq2,f,n)
267     int root;
268     QUADTREE qt;
269     BCOORD q0[3],q1[3],q2[3];
270     BCOORD t0[3],t1[3],t2[3];
271     BCOORD dt10[3],dt21[3],dt02[3];
272     unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
273     FUNC f;
274     int n;
275     {
276     double db;
277     unsigned int e0,e1,e2;
278     char al,ah,bl,bh,cl,ch,testl,testh;
279     char test_t0t1,test_t1t2,test_t2t0;
280     BCOORD a,b,c;
281    
282     /* First check if can trivial accept: if vertex in cell */
283     if(s0 & s1 & s2)
284     goto Lfunc_call;
285    
286     /* Assumption: Could not trivial reject*/
287     /* IF any coords lie on edge- must be in if couldnt reject */
288     a = q1[0];b= q0[1]; c= q0[2];
289     if(t0[0]== a || t0[1] == b || t0[2] == c)
290     goto Lfunc_call;
291     if(t1[0]== a || t1[1] == b || t1[2] == c)
292     goto Lfunc_call;
293     if(t2[0]== a || t2[1] == b || t2[2] == c)
294     goto Lfunc_call;
295 gwlarson 3.7
296 gwlarson 3.9 /* Test for edge crossings */
297     /* Test t0t1,t1t2,t2t0 against a */
298     /* Calculate edge crossings */
299     e0 = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
300     /* See if edge can be trivially rejected from intersetion testing */
301     test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
302     ((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
303     bl=bh=0;
304     if(test_t0t1 && (e0 & 2) )
305 gwlarson 3.7 {
306 gwlarson 3.9 /* Must do double calculation to avoid overflow */
307     db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
308     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
309 gwlarson 3.7 }
310 gwlarson 3.9 test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
311     ((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
312     if(test_t1t2 && (e0 & 1))
313     { /* test t1t2 against a */
314     db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
315     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
316     }
317     test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
318     ((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
319     if(test_t2t0 && (e0 & 4))
320     {
321     db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
322     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
323     }
324     e1 = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
325     cl = ch = 0;
326     if(test_t0t1 && (e1 & 2))
327     {/* test t0t1 against b */
328     db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
329     TEST_INT(testl,testh,db,c,q2[2],cl,ch)
330     }
331     if(test_t1t2 && (e1 & 1))
332     {/* test t1t2 against b */
333     db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
334     TEST_INT(testl,testh,db,c,q2[2],cl,ch)
335     }
336     if(test_t2t0 && (e1 & 4))
337     {/* test t2t0 against b */
338     db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
339     TEST_INT(testl,testh,db,c,q2[2],cl,ch)
340     }
341    
342     /* test edges against c */
343     e2 = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
344     al = ah = 0;
345     if(test_t0t1 && (e2 & 2))
346     { /* test t0t1 against c */
347     db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
348     TEST_INT(testl,testh,db,a,q0[0],al,ah)
349     }
350     if(test_t1t2 && (e2 & 1))
351     {
352     db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
353     TEST_INT(testl,testh,db,a,q0[0],al,ah)
354     }
355     if(test_t2t0 && (e2 & 4))
356     { /* test t2t0 against c */
357     db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
358     TEST_INT(testl,testh,db,a,q0[0],al,ah)
359     }
360     /* Only remaining case is if some edge trivially rejected */
361     if(!e0 || !e1 || !e2)
362     return(qt);
363 gwlarson 3.4
364 gwlarson 3.9 /* Only one/none got tested - pick either of other two/any two */
365     /* Only need to test one edge */
366     if(!test_t0t1 && (e0 & 2))
367     {
368     db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
369     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
370     }
371     if(!test_t1t2 && (e0 & 1))
372     {/* test t1t2 against a */
373     db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
374     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
375     }
376     if(!test_t2t0 && (e0 & 4))
377     {
378     db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
379     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
380     }
381 gwlarson 3.7
382 gwlarson 3.9 return(qt);
383    
384     Lfunc_call:
385     qt = f.func(f.argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
386     s0,s1,s2,sq0,sq1,sq2,0,f,n);
387     return(qt);
388    
389     }
390    
391    
392    
393     /* Leaf node: Do definitive test */
394     QUADTREE
395     qtLeaf_insert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
396     scale, s0,s1,s2,sq0,sq1,sq2,f,n)
397     int root;
398     QUADTREE qt;
399     BCOORD q0[3],q1[3],q2[3];
400     BCOORD t0[3],t1[3],t2[3];
401     BCOORD dt10[3],dt21[3],dt02[3];
402     unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
403     FUNC f;
404     int n;
405     {
406     double db;
407     unsigned int e0,e1,e2;
408     BCOORD a,b,c;
409     double p0[2], p1[2],cp;
410     char al,ah,bl,bh,cl,ch;
411     char testl,testh,test_t0t1,test_t1t2,test_t2t0;
412     unsigned int ls0,ls1,ls2;
413    
414     /* May have gotten passed trivial reject if one/two vertices are ON */
415     a = q1[0];b= q0[1]; c= q0[2];
416     SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
417    
418     /* First check if can trivial accept: if vertex in cell */
419     if(ls0 & ls1 & ls2)
420     goto Lfunc_call;
421    
422     if(ls0==0 || ls1 == 0 || ls2 ==0)
423     return(qt);
424     /* Assumption: Could not trivial reject*/
425     /* IF any coords lie on edge- must be in if couldnt reject */
426    
427     if(t0[0]== a || t0[1] == b || t0[2] == c)
428     goto Lfunc_call;
429     if(t1[0]== a || t1[1] == b || t1[2] == c)
430     goto Lfunc_call;
431     if(t2[0]== a || t2[1] == b || t2[2] == c)
432     goto Lfunc_call;
433    
434     /* Test for edge crossings */
435     /* Test t0t1 against a,b,c */
436     /* First test if t0t1 can be trivially rejected */
437     /* If both edges are outside bounds- dont need to test */
438    
439     /* Test t0t1,t1t2,t2t0 against a */
440     test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
441     ((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
442     e0 = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
443     bl=bh=0;
444     /* Test t0t1,t1t2,t2t0 against a */
445     if(test_t0t1 && (e0 & 2) )
446 gwlarson 3.4 {
447 gwlarson 3.9 db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
448     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
449 gwlarson 3.4 }
450 gwlarson 3.9 test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
451     ((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
452     if(test_t1t2 && (e0 & 1))
453     { /* test t1t2 against a */
454     db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
455     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
456     }
457     test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
458     ((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
459     if(test_t2t0 && (e0 & 4))
460     {
461     db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
462     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
463     }
464     cl = ch = 0;
465     e1 = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
466     if(test_t0t1 && (e1 & 2))
467     {/* test t0t1 against b */
468     db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
469     TEST_INT(testl,testh,db,q1[2],c,cl,ch)
470     }
471     if(test_t1t2 && (e1 & 1))
472     {/* test t1t2 against b */
473     db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
474     TEST_INT(testl,testh,db,q1[2],c,cl,ch)
475     }
476     if(test_t2t0 && (e1 & 4))
477     {/* test t2t0 against b */
478     db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
479     TEST_INT(testl,testh,db,q1[2],c,cl,ch)
480     }
481     al = ah = 0;
482     e2 = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
483     if(test_t0t1 && (e2 & 2))
484     { /* test t0t1 against c */
485     db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
486     TEST_INT(testl,testh,db,q0[0],a,al,ah)
487     }
488     if(test_t1t2 && (e2 & 1))
489     {
490     db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
491     TEST_INT(testl,testh,db,q0[0],a,al,ah)
492     }
493     if(test_t2t0 && (e2 & 4))
494     { /* test t2t0 against c */
495     db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
496     TEST_INT(testl,testh,db,q0[0],a,al,ah)
497     }
498     /* Only remaining case is if some edge trivially rejected */
499     if(!e0 || !e1 || !e2)
500     return(qt);
501 gwlarson 3.4
502 gwlarson 3.9 /* Only one/none got tested - pick either of other two/any two */
503     /* Only need to test one edge */
504     if(!test_t0t1 && (e0 & 2))
505 gwlarson 3.4 {
506 gwlarson 3.9 db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
507     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
508 gwlarson 3.4 }
509 gwlarson 3.9 if(!test_t1t2 && (e0 & 1))
510     {/* test t1t2 against a */
511     db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
512     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
513 gwlarson 3.4 }
514 gwlarson 3.9 if(!test_t2t0 && (e0 & 4))
515 gwlarson 3.7 {
516 gwlarson 3.9 db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
517     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
518 gwlarson 3.7 }
519 gwlarson 3.9
520 gwlarson 3.7 return(qt);
521 gwlarson 3.9 Lfunc_call:
522    
523     qt = f.func(f.argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
524     s0,s1,s2,sq0,sq1,sq2,1,f,n);
525     return(qt);
526     }
527    
528    
529    
530     /* ASSUMPTION: that triangle t0,t1,t2 intersects quad cell q0,q1,q2 */
531    
532     /*-------q2--------- sq2
533     / \
534     s1 /sc \ s0
535     qb_____qa
536     / \ / \
537     \sq0/sa \ /sb \ /sq1
538     \ q0____qc____q1/
539     \ /
540     \ s2 /
541     */
542    
543     int Find_id=0;
544    
545     QUADTREE
546     qtInsert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
547     s0,s1,s2,sq0,sq1,sq2,f,n)
548     int root;
549     QUADTREE qt;
550     BCOORD q0[3],q1[3],q2[3];
551     BCOORD t0[3],t1[3],t2[3];
552     BCOORD dt10[3],dt21[3],dt02[3];
553     BCOORD scale;
554     unsigned int s0,s1,s2,sq0,sq1,sq2;
555     FUNC f;
556     int n;
557     {
558     BCOORD a,b,c;
559     BCOORD va[3],vb[3],vc[3];
560     unsigned int sa,sb,sc;
561    
562     /* If a tree: trivial reject and recurse on potential children */
563     if(QT_IS_TREE(qt))
564     {
565     /* Test against new a edge: if entirely to left: only need
566     to visit child 0
567     */
568     a = q1[0] + scale;
569     b = q0[1] + scale;
570     c = q0[2] + scale;
571    
572     SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
573    
574     if(sa==7)
575     {
576     vb[1] = q0[1];
577     vc[2] = q0[2];
578     vc[1] = b;
579     vb[2] = c;
580     vb[0] = vc[0] = a;
581     QT_NTH_CHILD(qt,0) = qtInsert_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
582     vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
583     return(qt);
584     }
585     if(sb==7)
586     {
587     va[0] = q1[0];
588     vc[2] = q0[2];
589     va[1] = vc[1] = b;
590     va[2] = c;
591     vc[0] = a;
592     QT_NTH_CHILD(qt,1) = qtInsert_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,
593     t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
594     return(qt);
595     }
596     if(sc==7)
597     {
598     va[0] = q1[0];
599     vb[1] = q0[1];
600     va[1] = b;
601     va[2] = vb[2] = c;
602     vb[0] = a;
603     QT_NTH_CHILD(qt,2) = qtInsert_tri(root,QT_NTH_CHILD(qt,2),vb,va,
604     q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
605     return(qt);
606     }
607    
608     va[0] = q1[0];
609     vb[1] = q0[1];
610     vc[2] = q0[2];
611     va[1] = vc[1] = b;
612     va[2] = vb[2] = c;
613     vb[0] = vc[0] = a;
614     /* If outside of all edges: only need to Visit 3 */
615     if(sa==0 && sb==0 && sc==0)
616     {
617     QT_NTH_CHILD(qt,3) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,
618     vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
619     return(qt);
620     }
621    
622     if(sa)
623     QT_NTH_CHILD(qt,0) = qtInsert_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
624     t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
625     if(sb)
626     QT_NTH_CHILD(qt,1) = qtInsert_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
627     t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
628     if(sc)
629     QT_NTH_CHILD(qt,2) = qtInsert_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
630     t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
631     QT_NTH_CHILD(qt,3) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
632     t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
633     return(qt);
634     }
635     /* Leaf node: Do definitive test */
636     else
637     return(qt = qtLeaf_insert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
638     scale,s0,s1,s2,sq0,sq1,sq2,f,n));
639 gwlarson 3.4 }
640    
641    
642 gwlarson 3.7 QUADTREE
643 gwlarson 3.9 qtInsert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
644     s0,s1,s2,sq0,sq1,sq2,f,n)
645     int root;
646     QUADTREE qt;
647     BCOORD q0[3],q1[3],q2[3];
648     BCOORD t0[3],t1[3],t2[3];
649     BCOORD dt10[3],dt21[3],dt02[3];
650     BCOORD scale;
651     unsigned int s0,s1,s2,sq0,sq1,sq2;
652     FUNC f;
653 gwlarson 3.7 {
654 gwlarson 3.9 BCOORD a,b,c;
655     BCOORD va[3],vb[3],vc[3];
656     unsigned int sa,sb,sc;
657 gwlarson 3.4
658 gwlarson 3.9 /* If a tree: trivial reject and recurse on potential children */
659     if(QT_IS_TREE(qt))
660     {
661     /* Test against new a edge: if entirely to left: only need
662     to visit child 0
663     */
664     a = q1[0] - scale;
665     b = q0[1] - scale;
666     c = q0[2] - scale;
667 gwlarson 3.7
668 gwlarson 3.9 SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
669 gwlarson 3.7
670 gwlarson 3.9 if(sa==0)
671     {
672     vb[1] = q0[1];
673     vc[2] = q0[2];
674     vc[1] = b;
675     vb[2] = c;
676     vb[0] = vc[0] = a;
677    
678     QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
679     vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
680     return(qt);
681     }
682     if(sb==0)
683     {
684     va[0] = q1[0];
685     vc[2] = q0[2];
686     va[1] = vc[1] = b;
687     va[2] = c;
688     vc[0] = a;
689     QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
690     t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
691     return(qt);
692     }
693     if(sc==0)
694     {
695     va[0] = q1[0];
696     vb[1] = q0[1];
697     va[1] = b;
698     va[2] = vb[2] = c;
699     vb[0] = a;
700     QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
701     q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
702     return(qt);
703     }
704     va[0] = q1[0];
705     vb[1] = q0[1];
706     vc[2] = q0[2];
707     va[1] = vc[1] = b;
708     va[2] = vb[2] = c;
709     vb[0] = vc[0] = a;
710     /* If outside of all edges: only need to Visit 3 */
711     if(sa==7 && sb==7 && sc==7)
712     {
713     QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,
714     vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
715     return(qt);
716     }
717     if(sa!=7)
718     QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
719     t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
720     if(sb!=7)
721     QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
722     t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
723     if(sc!=7)
724     QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
725     t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
726    
727     QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
728     t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
729     return(qt);
730     }
731     /* Leaf node: Do definitive test */
732 gwlarson 3.7 else
733 gwlarson 3.9 return(qt = qtLeaf_insert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
734     scale,s0,s1,s2,sq0,sq1,sq2,f,n));
735     }
736    
737    
738    
739    
740     /*************************************************************************/
741     /* Visit code for applying functions: NOTE THIS IS THE SAME AS INSERT CODE:
742     except sets flags: wanted insert to be as efficient as possible: so
743     broke into two sets of routines
744     */
745    
746     qtVisit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
747     s0,s1,s2,sq0,sq1,sq2,f)
748     int root;
749     QUADTREE qt;
750     BCOORD q0[3],q1[3],q2[3];
751     BCOORD t0[3],t1[3],t2[3];
752     BCOORD dt10[3],dt21[3],dt02[3];
753     BCOORD scale;
754     unsigned int s0,s1,s2,sq0,sq1,sq2;
755     FUNC f;
756     {
757     BCOORD a,b,c;
758     BCOORD va[3],vb[3],vc[3];
759     unsigned int sa,sb,sc;
760    
761     /* If a tree: trivial reject and recurse on potential children */
762     if(QT_IS_TREE(qt))
763 gwlarson 3.7 {
764 gwlarson 3.9 QT_SET_FLAG(qt);
765    
766     /* Test against new a edge: if entirely to left: only need
767     to visit child 0
768     */
769     a = q1[0] + scale;
770     b = q0[1] + scale;
771     c = q0[2] + scale;
772    
773     SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
774    
775     if(sa==7)
776     {
777     vb[1] = q0[1];
778     vc[2] = q0[2];
779     vc[1] = b;
780     vb[2] = c;
781     vb[0] = vc[0] = a;
782     qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
783     vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f);
784     return;
785     }
786     if(sb==7)
787     {
788     va[0] = q1[0];
789     vc[2] = q0[2];
790     va[1] = vc[1] = b;
791     va[2] = c;
792     vc[0] = a;
793     qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,
794     t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f);
795     return;
796     }
797     if(sc==7)
798     {
799     va[0] = q1[0];
800     vb[1] = q0[1];
801     va[1] = b;
802     va[2] = vb[2] = c;
803     vb[0] = a;
804     qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,
805     q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f);
806     return;
807     }
808    
809     va[0] = q1[0];
810     vb[1] = q0[1];
811     vc[2] = q0[2];
812     va[1] = vc[1] = b;
813     va[2] = vb[2] = c;
814     vb[0] = vc[0] = a;
815     /* If outside of all edges: only need to Visit 3 */
816     if(sa==0 && sb==0 && sc==0)
817     {
818     qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,
819     vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f);
820     return;
821     }
822    
823     if(sa)
824     qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
825     t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f);
826     if(sb)
827     qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
828     t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f);
829     if(sc)
830     qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
831     t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f);
832     qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
833     t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f);
834 gwlarson 3.7 }
835 gwlarson 3.9 /* Leaf node: Do definitive test */
836     else
837     qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
838     scale,s0,s1,s2,sq0,sq1,sq2,f);
839    
840 gwlarson 3.7 }
841 gwlarson 3.4
842 gwlarson 3.3
843 gwlarson 3.9 qtVisit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
844     s0,s1,s2,sq0,sq1,sq2,f)
845     int root;
846     QUADTREE qt;
847     BCOORD q0[3],q1[3],q2[3];
848     BCOORD t0[3],t1[3],t2[3];
849     BCOORD dt10[3],dt21[3],dt02[3];
850     BCOORD scale;
851     unsigned int s0,s1,s2,sq0,sq1,sq2;
852     FUNC f;
853     {
854     BCOORD a,b,c;
855     BCOORD va[3],vb[3],vc[3];
856     unsigned int sa,sb,sc;
857 gwlarson 3.2
858 gwlarson 3.9 /* If a tree: trivial reject and recurse on potential children */
859     if(QT_IS_TREE(qt))
860     {
861     QT_SET_FLAG(qt);
862     /* Test against new a edge: if entirely to left: only need
863     to visit child 0
864     */
865     a = q1[0] - scale;
866     b = q0[1] - scale;
867     c = q0[2] - scale;
868 gwlarson 3.2
869 gwlarson 3.9 SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
870 gwlarson 3.2
871 gwlarson 3.9 if(sa==0)
872     {
873     vb[1] = q0[1];
874     vc[2] = q0[2];
875     vc[1] = b;
876     vb[2] = c;
877     vb[0] = vc[0] = a;
878     qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
879     vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f);
880     return;
881     }
882     if(sb==0)
883     {
884     va[0] = q1[0];
885     vc[2] = q0[2];
886     va[1] = vc[1] = b;
887     va[2] = c;
888     vc[0] = a;
889     qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
890     t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f);
891     return;
892     }
893     if(sc==0)
894     {
895     va[0] = q1[0];
896     vb[1] = q0[1];
897     va[1] = b;
898     va[2] = vb[2] = c;
899     vb[0] = a;
900     qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
901     q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f);
902     return;
903     }
904     va[0] = q1[0];
905     vb[1] = q0[1];
906     vc[2] = q0[2];
907     va[1] = vc[1] = b;
908     va[2] = vb[2] = c;
909     vb[0] = vc[0] = a;
910     /* If outside of all edges: only need to Visit 3 */
911     if(sa==7 && sb==7 && sc==7)
912     {
913     qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,
914     vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f);
915     return;
916     }
917     if(sa!=7)
918     qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
919     t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f);
920     if(sb!=7)
921     qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
922     t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f);
923     if(sc!=7)
924     qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
925     t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f);
926 gwlarson 3.2
927 gwlarson 3.9 qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
928     t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f);
929     return;
930     }
931     /* Leaf node: Do definitive test */
932     else
933     qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
934     scale,s0,s1,s2,sq0,sq1,sq2,f);
935     }
936 gwlarson 3.2
937    
938    
939 gwlarson 3.9 qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
940     scale, s0,s1,s2,sq0,sq1,sq2,f)
941     int root;
942     QUADTREE qt;
943     BCOORD q0[3],q1[3],q2[3];
944     BCOORD t0[3],t1[3],t2[3];
945     BCOORD dt10[3],dt21[3],dt02[3];
946     unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
947     FUNC f;
948     {
949     double db;
950     unsigned int e0,e1,e2;
951     char al,ah,bl,bh,cl,ch,testl,testh;
952     char test_t0t1,test_t1t2,test_t2t0;
953     BCOORD a,b,c;
954 gwlarson 3.2
955 gwlarson 3.9 /* First check if can trivial accept: if vertex in cell */
956     if(s0 & s1 & s2)
957     goto Lfunc_call;
958    
959     /* Assumption: Could not trivial reject*/
960     /* IF any coords lie on edge- must be in if couldnt reject */
961     a = q1[0];b= q0[1]; c= q0[2];
962     if(t0[0]== a || t0[1] == b || t0[2] == c)
963     goto Lfunc_call;
964     if(t1[0]== a || t1[1] == b || t1[2] == c)
965     goto Lfunc_call;
966     if(t2[0]== a || t2[1] == b || t2[2] == c)
967     goto Lfunc_call;
968    
969     /* Test for edge crossings */
970     /* Test t0t1,t1t2,t2t0 against a */
971     /* Calculate edge crossings */
972     e0 = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
973     /* See if edge can be trivially rejected from intersetion testing */
974     test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
975     ((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
976     bl=bh=0;
977     if(test_t0t1 && (e0 & 2) )
978     {
979     /* Must do double calculation to avoid overflow */
980     db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
981     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
982     }
983     test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
984     ((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
985     if(test_t1t2 && (e0 & 1))
986     { /* test t1t2 against a */
987     db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
988     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
989     }
990     test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
991     ((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
992     if(test_t2t0 && (e0 & 4))
993     {
994     db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
995     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
996     }
997     e1 = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
998     cl = ch = 0;
999     if(test_t0t1 && (e1 & 2))
1000     {/* test t0t1 against b */
1001     db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
1002     TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1003     }
1004     if(test_t1t2 && (e1 & 1))
1005     {/* test t1t2 against b */
1006     db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
1007     TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1008     }
1009     if(test_t2t0 && (e1 & 4))
1010     {/* test t2t0 against b */
1011     db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
1012     TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1013     }
1014    
1015     /* test edges against c */
1016     e2 = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
1017     al = ah = 0;
1018     if(test_t0t1 && (e2 & 2))
1019     { /* test t0t1 against c */
1020     db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
1021     TEST_INT(testl,testh,db,a,q0[0],al,ah)
1022     }
1023     if(test_t1t2 && (e2 & 1))
1024     {
1025     db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1026     TEST_INT(testl,testh,db,a,q0[0],al,ah)
1027     }
1028     if(test_t2t0 && (e2 & 4))
1029     { /* test t2t0 against c */
1030     db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1031     TEST_INT(testl,testh,db,a,q0[0],al,ah)
1032     }
1033     /* Only remaining case is if some edge trivially rejected */
1034     if(!e0 || !e1 || !e2)
1035     return;
1036    
1037     /* Only one/none got tested - pick either of other two/any two */
1038     /* Only need to test one edge */
1039     if(!test_t0t1 && (e0 & 2))
1040     {
1041     db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1042     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1043     }
1044     if(!test_t1t2 && (e0 & 1))
1045     {/* test t1t2 against a */
1046     db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1047     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1048     }
1049     if(!test_t2t0 && (e0 & 4))
1050     {
1051     db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1052     TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1053     }
1054    
1055     return;
1056    
1057     Lfunc_call:
1058     f.func(f.argptr,root,qt);
1059     if(!QT_IS_EMPTY(qt))
1060     QT_LEAF_SET_FLAG(qt);
1061     }
1062    
1063    
1064    
1065     /* Leaf node: Do definitive test */
1066     QUADTREE
1067     qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1068     scale, s0,s1,s2,sq0,sq1,sq2,f)
1069     int root;
1070     QUADTREE qt;
1071     BCOORD q0[3],q1[3],q2[3];
1072     BCOORD t0[3],t1[3],t2[3];
1073     BCOORD dt10[3],dt21[3],dt02[3];
1074     unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
1075     FUNC f;
1076     {
1077     double db;
1078     unsigned int e0,e1,e2;
1079     BCOORD a,b,c;
1080     double p0[2], p1[2],cp;
1081     char al,ah,bl,bh,cl,ch;
1082     char testl,testh,test_t0t1,test_t1t2,test_t2t0;
1083     unsigned int ls0,ls1,ls2;
1084    
1085     /* May have gotten passed trivial reject if one/two vertices are ON */
1086     a = q1[0];b= q0[1]; c= q0[2];
1087     SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
1088    
1089     /* First check if can trivial accept: if vertex in cell */
1090     if(ls0 & ls1 & ls2)
1091     goto Lfunc_call;
1092    
1093     if(ls0==0 || ls1 == 0 || ls2 ==0)
1094     return;
1095     /* Assumption: Could not trivial reject*/
1096     /* IF any coords lie on edge- must be in if couldnt reject */
1097    
1098     if(t0[0]== a || t0[1] == b || t0[2] == c)
1099     goto Lfunc_call;
1100     if(t1[0]== a || t1[1] == b || t1[2] == c)
1101     goto Lfunc_call;
1102     if(t2[0]== a || t2[1] == b || t2[2] == c)
1103     goto Lfunc_call;
1104    
1105     /* Test for edge crossings */
1106     /* Test t0t1 against a,b,c */
1107     /* First test if t0t1 can be trivially rejected */
1108     /* If both edges are outside bounds- dont need to test */
1109    
1110     /* Test t0t1,t1t2,t2t0 against a */
1111     test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
1112     ((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
1113     e0 = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
1114     bl=bh=0;
1115     /* Test t0t1,t1t2,t2t0 against a */
1116     if(test_t0t1 && (e0 & 2) )
1117     {
1118     db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
1119     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1120     }
1121     test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
1122     ((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
1123     if(test_t1t2 && (e0 & 1))
1124     { /* test t1t2 against a */
1125     db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1126     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1127     }
1128     test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
1129     ((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
1130     if(test_t2t0 && (e0 & 4))
1131     {
1132     db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1133     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1134     }
1135     cl = ch = 0;
1136     e1 = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
1137     if(test_t0t1 && (e1 & 2))
1138     {/* test t0t1 against b */
1139     db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
1140     TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1141     }
1142     if(test_t1t2 && (e1 & 1))
1143     {/* test t1t2 against b */
1144     db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
1145     TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1146     }
1147     if(test_t2t0 && (e1 & 4))
1148     {/* test t2t0 against b */
1149     db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
1150     TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1151     }
1152     al = ah = 0;
1153     e2 = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
1154     if(test_t0t1 && (e2 & 2))
1155     { /* test t0t1 against c */
1156     db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
1157     TEST_INT(testl,testh,db,q0[0],a,al,ah)
1158     }
1159     if(test_t1t2 && (e2 & 1))
1160     {
1161     db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1162     TEST_INT(testl,testh,db,q0[0],a,al,ah)
1163     }
1164     if(test_t2t0 && (e2 & 4))
1165     { /* test t2t0 against c */
1166     db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1167     TEST_INT(testl,testh,db,q0[0],a,al,ah)
1168     }
1169     /* Only remaining case is if some edge trivially rejected */
1170     if(!e0 || !e1 || !e2)
1171     return;
1172    
1173     /* Only one/none got tested - pick either of other two/any two */
1174     /* Only need to test one edge */
1175     if(!test_t0t1 && (e0 & 2))
1176     {
1177     db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1178     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1179     }
1180     if(!test_t1t2 && (e0 & 1))
1181     {/* test t1t2 against a */
1182     db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1183     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1184     }
1185     if(!test_t2t0 && (e0 & 4))
1186     {
1187     db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1188     TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1189     }
1190    
1191     return;
1192    
1193     Lfunc_call:
1194     f.func(f.argptr,root,qt);
1195     if(!QT_IS_EMPTY(qt))
1196     QT_LEAF_SET_FLAG(qt);
1197     }
1198 gwlarson 3.2
1199    
1200