ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/radiance/ray/src/hd/sm_qtree.c

Comparing ray/src/hd/sm_qtree.c (file contents):
Revision 3.3 by gwlarson, Tue Aug 25 11:03:27 1998 UTC vs.
Revision 3.11 by gwlarson, Sun Jan 10 10:27:45 1999 UTC

#	Line 14 | Line 14 | static char SCCSid[] = "$SunId$ SGI";
14		*/
15
16		#include "standard.h"
17	<
17	>	#include "sm_flag.h"
18		#include "sm_geom.h"
19		#include "sm_qtree.h"
20	–	#include "object.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	<	int4 *quad_flag;
24	>	int4 *quad_flag= NULL;
25
26	+	#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	+	extern int Pick_q[500];
31	+
32	+	#endif
33	+	int Incnt=0;
34	+
35		QUADTREE
36		qtAlloc()                       /* allocate a quadtree */
37		{
#	Line 41 | Line 49 | qtAlloc()                      /* allocate a quadtree */
49		return(EMPTY);
50		if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
51		QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
52	<	return(EMPTY);
52	>	error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
53	>
54	>	/* Realloc the per/node flags */
55		quad_flag = (int4 *)realloc((char *)quad_flag,
56	<	(QT_BLOCK(freet)+1)*QT_BLOCK_SIZE/(8*sizeof(int4)));
56	>	(QT_BLOCK(freet)+1)*((QT_BLOCK_SIZE+7)>>3));
57		if (quad_flag == NULL)
58	<	return(EMPTY);
58	>	error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
59		}
60		treetop += 4;
61		return(freet);
#	Line 54 | Line 64 | qtAlloc()                      /* allocate a quadtree */
64
65		qtClearAllFlags()               /* clear all quadtree branch flags */
66		{
67	<	if (!treetop) return;
68	<	bzero((char *)quad_flag,
69	<	(QT_BLOCK(treetop-1)+1)*QT_BLOCK_SIZE/(8*sizeof(int4)));
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
62	–
76		qtFree(qt)                      /* free a quadtree */
77		register QUADTREE  qt;
78		{
#	Line 80 | Line 93 | qtFree(qt)                     /* free a quadtree */
93		qtDone()                        /* free EVERYTHING */
94		{
95		register int    i;
96	<
96	>
97	>	qtfreeleaves();
98		for (i = 0; i < QT_MAX_BLK; i++)
99		{
100		if (quad_block[i] == NULL)
#	Line 88 | Line 102 | qtDone()                       /* free EVERYTHING */
102		free((char *)quad_block[i]);
103		quad_block[i] = NULL;
104		}
105	+	/* Free the flags */
106		if (i) free((char *)quad_flag);
107		quad_flag = NULL;
108		quad_free_list = EMPTY;
109		treetop = 0;
110		}
111
97	–
112		QUADTREE
113	<	qtCompress(qt)                  /* recursively combine nodes */
114	<	register QUADTREE  qt;
113	>	qtLocate(qt,bcoord)
114	>	QUADTREE qt;
115	>	BCOORD bcoord[3];
116		{
117	<	register int  i;
103	<	register QUADTREE  qres;
117	>	int i;
118
119	<	if (!QT_IS_TREE(qt))        /* not a tree */
120	<	return(qt);
121	<	qres = QT_NTH_CHILD(qt,0) = qtCompress(QT_NTH_CHILD(qt,0));
122	<	for (i = 1; i < 4; i++)
123	<	if((QT_NTH_CHILD(qt,i) = qtCompress(QT_NTH_CHILD(qt,i))) != qres)
124	<	qres = qt;
125	<	if(!QT_IS_TREE(qres))
126	<	{   /* all were identical leaves */
127	<	QT_NTH_CHILD(qt,0) = quad_free_list;
128	<	quad_free_list = qt;
119	>	if(QT_IS_TREE(qt))
120	>	{
121	>	i = bary_child(bcoord);
122	>
123	>	return(qtLocate(QT_NTH_CHILD(qt,i),bcoord));
124	>	}
125	>	else
126	>	return(qt);
127	>	}
128	>
129	>	int
130	>	move_to_nbr(b,db0,db1,db2,tptr)
131	>	BCOORD b[3];
132	>	BDIR db0,db1,db2;
133	>	double *tptr;
134	>	{
135	>	double t,dt;
136	>	BCOORD bc;
137	>	int nbr;
138	>
139	>	nbr = -1;
140	>	*tptr = 0.0;
141	>	/* Advance to next node */
142	>	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	>	if(db0 < 0)
149	>	{
150	>	t = ((double)b[0])/-db0;
151	>	nbr = 0;
152	>	}
153	>	else
154	>	t = MAXFLOAT;
155	>	if(db1 < 0)
156	>	{
157	>	dt = ((double)b[1])/-db1;
158	>	if( dt < t)
159	>	{
160	>	t = dt;
161	>	nbr = 1;
162		}
163	<	return(qres);
163	>	}
164	>	if(db2 < 0)
165	>	{
166	>	dt = ((double)b[2])/-db2;
167	>	if( dt < t)
168	>	{
169	>	t = dt;
170	>	nbr = 2;
171	>	}
172	>	}
173	>	*tptr = t;
174	>	return(nbr);
175		}
176
177	<
178	<	QUADTREE
179	<	*qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
180	<	QUADTREE *qtptr;
181	<	double bcoord[3];
182	<	FVECT t0,t1,t2;
177	>	qtTrace_ray(qt,b,db0,db1,db2,nextptr,sign,sfactor,func,f)
178	>	QUADTREE qt;
179	>	BCOORD b[3];
180	>	BDIR db0,db1,db2;
181	>	int *nextptr;
182	>	int sign,sfactor;
183	>	FUNC func;
184	>	int *f;
185		{
186	<	int i;
187	<	QUADTREE *child;
188	<	FVECT a,b,c;
186	>	int i,found;
187	>	QUADTREE child;
188	>	int nbr,next,w;
189	>	double t;
190
191	<	if(QT_IS_TREE(*qtptr))
192	<	{
193	<	i = bary2d_child(bcoord);
194	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
195	<	if(t0)
191	>	if(QT_IS_TREE(qt))
192	>	{
193	>	/* Find the appropriate child and reset the coord */
194	>	i = bary_child(b);
195	>
196	>	for(;;)
197		{
198	<	qtSubdivide_tri(t0,t1,t2,a,b,c);
199	<	qtNth_child_tri(t0,t1,t2,a,b,c,i,t0,t1,t2);
198	>	child = QT_NTH_CHILD(qt,i);
199	>	if(i != 3)
200	>	qtTrace_ray(child,b,db0,db1,db2,nextptr,sign,sfactor+1,func,f);
201	>	else
202	>	/* If the center cell- must flip direction signs */
203	>	qtTrace_ray(child,b,-db0,-db1,-db2,nextptr,1-sign,sfactor+1,func,f);
204	>
205	>	if(QT_FLAG_IS_DONE(*f))
206	>	return;
207	>	/* 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	>	bary_parent(b,i);
217	>	/* Else pop up to parent and traverse from there */
218	>	return(qt);
219	>	}
220	>	bary_from_child(b,i,next);
221	>	i = next;
222		}
139	–	return(qtLocate_leaf(child,bcoord,t0,t1,t2));
223		}
224		else
225	<	return(qtptr);
226	<	}
225	>	{
226	>	#ifdef TEST_DRIVER
227	>	if(Pick_cnt < 500)
228	>	Pick_q[Pick_cnt++]=qt;
229	>	#endif;
230	>	F_FUNC(func)(qt,F_ARGS(func),f);
231	>	if(QT_FLAG_IS_DONE(*f))
232	>	return;
233	>	/* Advance to next node */
234	>	nbr = move_to_nbr(b,db0,db1,db2,&t);
235
236	+	if(nbr==-1)
237	+	{
238	+	QT_FLAG_SET_DONE(*f);
239	+	return;
240	+	}
241	+	b[0] += (BCOORD)(t *db0);
242	+	b[1] += (BCOORD)(t *db1);
243	+	b[2] += (BCOORD)(t *db2);
244	+	*nextptr = nbr;
245	+	return;
246	+
247	+	}
248	+	}
249
250
251	<	int
252	<	qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n)
253	<	QUADTREE *qtptr;
254	<	double bcoord[3];
255	<	int id;
256	<	FVECT v0,v1,v2;
251	>
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	>	QUADTREE
265	>	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	<	int i;
277	<	QUADTREE *child;
278	<	OBJECT os[MAXSET+1],*optr;
279	<	int found;
280	<	FVECT r0,r1,r2;
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
296	<	if(QT_IS_TREE(*qtptr))
297	<	{
298	<	QT_SET_FLAG(*qtptr);
299	<	i = bary2d_child(bcoord);
300	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
301	<	return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,++n));
296	>	/* 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	>	{
306	>	/* 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		}
310	<	else
311	<	{
312	<	/* If this leave node emptry- create a new set */
313	<	if(QT_IS_EMPTY(*qtptr))
314	<	*qtptr = qtaddelem(*qtptr,id);
315	<	else
174	<	{
175	<	qtgetset(os,*qtptr);
176	<	/* If the set is too large: subdivide */
177	<	if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
178	<	*qtptr = qtaddelem(*qtptr,id);
179	<	else
180	<	{
181	<	if (n < QT_MAX_LEVELS)
182	<	{
183	<	/* If set size exceeds threshold: subdivide cell and
184	<	reinsert set tris into cell
185	<	*/
186	<	n++;
187	<	qtfreeleaf(*qtptr);
188	<	qtSubdivide(qtptr);
189	<	QT_SET_FLAG(*qtptr);
190	<	found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n);
191	<	#if 0
192	<	if(!found)
193	<	{
194	<	#ifdef TEST_DRIVER
195	<	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
196	<	#else
197	<	eputs("qtAdd_tri():Found in parent but not children\n");
198	<	#endif
199	<	}
200	<	#endif
201	<	for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
202	<	{
203	<	id = QT_SET_NEXT_ELEM(optr);
204	<	qtTri_verts_from_id(id,r0,r1,r2);
205	<	found= qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,++n);
206	<	#ifdef DEBUG
207	<	if(!found)
208	<	eputs("qtAdd_tri():Reinsert-in parent but not children\n");
209	<	#endif
210	<	}
211	<	}
212	<	else
213	<	if(QT_SET_CNT(os) < QT_MAX_SET)
214	<	{
215	<	*qtptr = qtaddelem(*qtptr,id);
216	<	}
217	<	else
218	<	{
219	<	#ifdef DEBUG
220	<	eputs("qtAdd_tri():two many levels\n");
221	<	#endif
222	<	return(FALSE);
223	<	}
224	<	}
225	<	}
310	>	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	<	return(TRUE);
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	>
364	>	/* 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	>
382	>	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
231	–	int
232	–	qtAdd_tri_from_point(qtptr,v0,v1,v2,pt,id)
233	–	QUADTREE *qtptr;
234	–	FVECT v0,v1,v2;
235	–	FVECT pt;
236	–	int id;
237	–	{
238	–	char d,w;
239	–	int i,x,y;
240	–	QUADTREE *child;
241	–	QUADTREE qt;
242	–	FVECT i_pt,n,a,b,c,r0,r1,r2;
243	–	double pd,bcoord[3];
244	–	OBJECT os[MAXSET+1],*optr;
245	–	int found;
246	–
247	–	/* Determine if point lies within pyramid (and therefore
248	–	inside a spherical quadtree cell):GT_INTERIOR, on one of the
249	–	pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
250	–	or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
251	–	For each triangle edge: compare the
252	–	point against the plane formed by the edge and the view center
253	–	*/
254	–	d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);
392
393	<	/* Not in this triangle */
394	<	if(!d)
395	<	return(FALSE);
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	<	/* Will return lowest level triangle containing point: It the
423	<	point is on an edge or vertex: will return first associated
424	<	triangle encountered in the child traversal- the others can
425	<	be derived using triangle adjacency information
264	<	*/
265	<	if(QT_IS_TREE(*qtptr))
266	<	{
267	<	QT_SET_FLAG(*qtptr);
268	<	/* Find the intersection point */
269	<	tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
270	<	intersect_vector_plane(pt,n,pd,NULL,i_pt);
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	<	i = max_index(n);
428	<	x = (i+1)%3;
429	<	y = (i+2)%3;
430	<	/* Calculate barycentric coordinates of i_pt */
431	<	bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
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	<	i = bary2d_child(bcoord);
435	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
436	<	/* NOTE: Optimize: only subdivide for specified child */
437	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
438	<	qtNth_child_tri(v0,v1,v2,a,b,c,i,v0,v1,v2);
439	<	return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,1));
440	<	}
441	<	else
442	<	{
443	<	/* If this leave node emptry- create a new set */
444	<	if(QT_IS_EMPTY(*qtptr))
445	<	*qtptr = qtaddelem(*qtptr,id);
446	<	else
447	<	{
448	<	qtgetset(os,*qtptr);
293	<	/* If the set is too large: subdivide */
294	<	if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
295	<	*qtptr = qtaddelem(*qtptr,id);
296	<	else
297	<	{
298	<	/* If set size exceeds threshold: subdivide cell and
299	<	reinsert set tris into cell
300	<	*/
301	<	qtfreeleaf(*qtptr);
302	<	qtSubdivide(qtptr);
303	<	found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,1);
304	<	#if 0
305	<	if(!found)
306	<	{
307	<	#ifdef TEST_DRIVER
308	<	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
309	<	#else
310	<	eputs("qtAdd_tri():Found in parent but not children\n");
311	<	#endif
312	<	}
313	<	#endif
314	<	for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
315	<	{
316	<	id = QT_SET_NEXT_ELEM(optr);
317	<	qtTri_verts_from_id(id,r0,r1,r2);
318	<	found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,1);
319	<	#ifdef DEBUG
320	<	if(!found)
321	<	eputs("qtAdd_tri():Reinsert-in parent but not children\n");
322	<	#endif
323	<	}
324	<	}
325	<	}
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	>	{
447	>	db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
448	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
449		}
450	<	return(TRUE);
451	<	}
450	>	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
502	+	/* 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	+	{
506	+	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
507	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
508	+	}
509	+	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	+	}
514	+	if(!test_t2t0 && (e0 & 4))
515	+	{
516	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
517	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
518	+	}
519
520	<	QUADTREE
521	<	*qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
333	<	QUADTREE *qtptr;
334	<	FVECT v0,v1,v2;
335	<	FVECT pt;
336	<	FVECT t0,t1,t2;
337	<	{
338	<	char d,w;
339	<	int i,x,y;
340	<	QUADTREE *child;
341	<	QUADTREE qt;
342	<	FVECT n,i_pt,a,b,c;
343	<	double pd,bcoord[3];
520	>	return(qt);
521	>	Lfunc_call:
522
523	<	/* Determine if point lies within pyramid (and therefore
524	<	inside a spherical quadtree cell):GT_INTERIOR, on one of the
525	<	pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
526	<	or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
349	<	For each triangle edge: compare the
350	<	point against the plane formed by the edge and the view center
351	<	*/
352	<	d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);
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
354	–	/* Not in this triangle */
355	–	if(!d)
356	–	{
357	–	return(NULL);
358	–	}
528
360	–	/* Will return lowest level triangle containing point: It the
361	–	point is on an edge or vertex: will return first associated
362	–	triangle encountered in the child traversal- the others can
363	–	be derived using triangle adjacency information
364	–	*/
365	–	if(QT_IS_TREE(*qtptr))
366	–	{
367	–	/* Find the intersection point */
368	–	tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
369	–	intersect_vector_plane(pt,n,pd,NULL,i_pt);
529
530	<	i = max_index(n);
372	<	x = (i+1)%3;
373	<	y = (i+2)%3;
374	<	/* Calculate barycentric coordinates of i_pt */
375	<	bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
530	>	/* ASSUMPTION: that triangle t0,t1,t2 intersects quad cell q0,q1,q2 */
531
532	<	i = bary2d_child(bcoord);
533	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
534	<	if(t0)
535	<	{
536	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
537	<	qtNth_child_tri(v0,v1,v2,a,b,c,i,t0,t1,t2);
538	<	}
539	<	return(qtLocate_leaf(child,bcoord,t0,t1,t2));
540	<	}
541	<	else
387	<	return(qtptr);
388	<	}
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
544	<	qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
545	<	QUADTREE *qtptr;
546	<	FVECT v0,v1,v2;
547	<	FVECT pt;
548	<	char *type,*which;
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	<	OBJECT os[MAXSET+1],*optr;
559	<	char d,w;
560	<	int i,id;
400	<	FVECT p0,p1,p2;
558	>	BCOORD a,b,c;
559	>	BCOORD va[3],vb[3],vc[3];
560	>	unsigned int sa,sb,sc;
561
562	<	qtptr = qtRoot_point_locate(qtptr,v0,v1,v2,pt,NULL,NULL,NULL);
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	<	if(!qtptr || QT_IS_EMPTY(*qtptr))
573	<	return(EMPTY);
574	<
407	<	/* Get the set */
408	<	qtgetset(os,*qtptr);
409	<	for (i = QT_SET_CNT(os),optr = QT_SET_PTR(os); i > 0; i--)
572	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
573	>
574	>	if(sa==7)
575		{
576	<	/* Find the triangle that pt falls in (NOTE:FOR now return first 1) */
577	<	id = QT_SET_NEXT_ELEM(optr);
578	<	qtTri_verts_from_id(id,p0,p1,p2);
579	<	d = test_single_point_against_spherical_tri(p0,p1,p2,pt,&w);
580	<	if(d)
581	<	{
582	<	if(type)
583	<	*type = d;
419	<	if(which)
420	<	*which = w;
421	<	return(id);
422	<	}
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	<	return(EMPTY);
586	<	}
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	<	QUADTREE
609	<	qtSubdivide(qtptr)
610	<	QUADTREE *qtptr;
611	<	{
612	<	QUADTREE node;
613	<	node = qtAlloc();
614	<	QT_CLEAR_CHILDREN(node);
615	<	*qtptr = node;
616	<	return(node);
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		}
640
641
642		QUADTREE
643	<	qtSubdivide_nth_child(qt,n)
643	>	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	<	int n;
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		{
654	<	QUADTREE node;
654	>	BCOORD a,b,c;
655	>	BCOORD va[3],vb[3],vc[3];
656	>	unsigned int sa,sb,sc;
657
658	<	node = qtSubdivide(&(QT_NTH_CHILD(qt,n)));
659	<
660	<	return(node);
661	<	}
658	>	/* 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
668	<	/* for triangle v0-v1-v2- returns a,b,c: children are:
668	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
669
670	<	v2                     0: v0,a,c
671	<	/\                     1: a,v1,b
672	<	/2 \                    2: c,b,v2
673	<	c/____\b                  3: b,c,a
674	<	/\    /\
675	<	/0 \3 /1 \
676	<	v0____\/____\v1
460	<	a
461	<	*/
670	>	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	<	qtSubdivide_tri(v0,v1,v2,a,b,c)
679	<	FVECT v0,v1,v2;
680	<	FVECT a,b,c;
681	<	{
682	<	EDGE_MIDPOINT_VEC3(a,v0,v1);
683	<	EDGE_MIDPOINT_VEC3(b,v1,v2);
684	<	EDGE_MIDPOINT_VEC3(c,v2,v0);
685	<	}
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	<	qtNth_child_tri(v0,v1,v2,a,b,c,i,r0,r1,r2)
728	<	FVECT v0,v1,v2;
729	<	FVECT a,b,c;
475	<	int i;
476	<	FVECT r0,r1,r2;
477	<	{
478	<	switch(i){
479	<	case 0:
480	<	VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
481	<	break;
482	<	case 1:
483	<	VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
484	<	break;
485	<	case 2:
486	<	VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
487	<	break;
488	<	case 3:
489	<	VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
490	<	break;
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	+	else
733	+	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
494	–	/* Add triangle "id" to all leaf level cells that are children of
495	–	quadtree pointed to by "qtptr" with cell vertices "t1,t2,t3"
496	–	that it overlaps (vertex and edge adjacencies do not count
497	–	as overlapping). If the addition of the triangle causes the cell to go over
498	–	threshold- the cell is split- and the triangle must be recursively inserted
499	–	into the new child cells: it is assumed that "v1,v2,v3" are normalized
500	–	*/
737
502	–	int
503	–	qtRoot_add_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
504	–	QUADTREE *qtptr;
505	–	int id;
506	–	FVECT t0,t1,t2;
507	–	FVECT v0,v1,v2;
508	–	int n;
509	–	{
510	–	char test;
511	–	int found;
738
513	–	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
514	–	if(!test)
515	–	return(FALSE);
516	–
517	–	found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
518	–	return(found);
519	–	}
739
740	<	int
741	<	qtRoot_add_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n)
742	<	QUADTREE *qtptr;
743	<	int id;
744	<	FVECT t0,t1,t2;
526	<	FVECT v0,v1,v2;
527	<	int n;
528	<	{
529	<	char test;
530	<	int found;
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	<	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
747	<	if(!test)
748	<	return(FALSE);
749	<
750	<	found = qtAdd_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n);
751	<	return(found);
752	<	}
753	<
754	<	int
755	<	qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
542	<	QUADTREE *qtptr;
543	<	int id;
544	<	FVECT t0,t1,t2;
545	<	FVECT v0,v1,v2;
546	<	int n;
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	<
758	<	char test;
759	<	int i,index;
551	<	FVECT a,b,c;
552	<	OBJECT os[MAXSET+1],*optr;
553	<	QUADTREE qt;
554	<	int found;
555	<	FVECT r0,r1,r2;
757	>	BCOORD a,b,c;
758	>	BCOORD va[3],vb[3],vc[3];
759	>	unsigned int sa,sb,sc;
760
761	<	found = 0;
762	<	/* if this is tree: recurse */
559	<	if(QT_IS_TREE(*qtptr))
761	>	/* If a tree: trivial reject and recurse on potential children */
762	>	if(QT_IS_TREE(qt))
763		{
764	<	n++;
562	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
563	<	test = spherical_tri_intersect(t0,t1,t2,v0,a,c);
564	<	if(test)
565	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c,n);
566	<	test = spherical_tri_intersect(t0,t1,t2,a,v1,b);
567	<	if(test)
568	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b,n);
569	<	test = spherical_tri_intersect(t0,t1,t2,c,b,v2);
570	<	if(test)
571	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2,n);
572	<	test = spherical_tri_intersect(t0,t1,t2,b,c,a);
573	<	if(test)
574	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a,n);
764	>	QT_SET_FLAG(qt);
765
766	<	#if 0
767	<	if(!found)
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	<	#ifdef TEST_DRIVER
778	<	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
779	<	#else
780	<	eputs("qtAdd_tri():Found in parent but not children\n");
781	<	#endif
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	<	#endif
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		}
835	+	/* Leaf node: Do definitive test */
836		else
837	<	{
838	<	/* If this leave node emptry- create a new set */
839	<	if(QT_IS_EMPTY(*qtptr))
591	<	*qtptr = qtaddelem(*qtptr,id);
592	<	else
593	<	{
594	<	qtgetset(os,*qtptr);
595	<	/* If the set is too large: subdivide */
596	<	if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
597	<	*qtptr = qtaddelem(*qtptr,id);
598	<	else
599	<	{
600	<	if (n < QT_MAX_LEVELS)
601	<	{
602	<	/* If set size exceeds threshold: subdivide cell and
603	<	reinsert set tris into cell
604	<	*/
605	<	n++;
606	<	qtfreeleaf(*qtptr);
607	<	qtSubdivide(qtptr);
608	<	found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
609	<	#if 0
610	<	if(!found)
611	<	{
612	<	#ifdef TEST_DRIVER
613	<	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
614	<	#else
615	<	eputs("qtAdd_tri():Found in parent but not children\n");
616	<	#endif
617	<	}
618	<	#endif
619	<	for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
620	<	{
621	<	id = QT_SET_NEXT_ELEM(optr);
622	<	qtTri_verts_from_id(id,r0,r1,r2);
623	<	found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,n);
624	<	#ifdef DEBUG
625	<	if(!found)
626	<	eputs("qtAdd_tri():Reinsert-in parent but not children\n");
627	<	#endif
628	<	}
629	<	}
630	<	else
631	<	if(QT_SET_CNT(os) < QT_MAX_SET)
632	<	{
633	<	*qtptr = qtaddelem(*qtptr,id);
634	<	#if 0
635	<	{
636	<	int k;
637	<	qtgetset(os,*qtptr);
638	<	printf("\n%d:\n",os[0]);
639	<	for(k=1; k <= os[0];k++)
640	<	printf("%d ",os[k]);
641	<	printf("\n");
642	<	}
643	<	#endif
644	<	/*
645	<	insertelem(os,id);
646	<	*qtptr = fullnode(os);
647	<	*/
648	<	}
649	<	else
650	<	{
651	<	#ifdef DEBUG
652	<	eputs("qtAdd_tri():two many levels\n");
653	<	#endif
654	<	return(FALSE);
655	<	}
656	<	}
657	<	}
658	<	}
659	<	return(TRUE);
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		}
841
842
843	<	int
844	<	qtApply_to_tri_cells(qtptr,t0,t1,t2,v0,v1,v2,func,arg)
845	<	QUADTREE *qtptr;
846	<	FVECT t0,t1,t2;
847	<	FVECT v0,v1,v2;
848	<	int (*func)();
849	<	char *arg;
843	>	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	<	char test;
855	<	FVECT a,b,c;
854	>	BCOORD a,b,c;
855	>	BCOORD va[3],vb[3],vc[3];
856	>	unsigned int sa,sb,sc;
857
858	<	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
859	<	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
858	>	/* 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
869	<	/* If triangles do not overlap: done */
678	<	if(!test)
679	<	return(FALSE);
869	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
870
871	<	/* if this is tree: recurse */
872	<	if(QT_IS_TREE(*qtptr))
873	<	{
874	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
875	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t0,t1,t2,v0,a,c,func,arg);
876	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t0,t1,t2,a,v1,b,func,arg);
877	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t0,t1,t2,c,b,v2,func,arg);
878	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t0,t1,t2,b,c,a,func,arg);
871	>	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	>
927	>	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	<	return(func(qtptr,arg));
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
937
695	–	int
696	–	qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
697	–	QUADTREE *qtptr;
698	–	int id;
699	–	FVECT t0,t1,t2;
700	–	FVECT v0,v1,v2;
701	–	{
702	–
703	–	char test;
704	–	int i;
705	–	FVECT a,b,c;
706	–	OBJECT os[MAXSET+1];
707	–
708	–	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
709	–	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
938
939	<	/* If triangles do not overlap: done */
940	<	if(!test)
941	<	return(FALSE);
939	>	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
955	<	/* if this is tree: recurse */
956	<	if(QT_IS_TREE(*qtptr))
955	>	/* 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	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
980	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c);
981	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b);
721	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2);
722	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a);
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	<	else
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	<	if(QT_IS_EMPTY(*qtptr))
995	<	{
728	<	#ifdef DEBUG
729	<	eputs("qtRemove_tri(): triangle not found\n");
730	<	#endif
731	<	}
732	<	/* remove id from set */
733	<	else
734	<	{
735	<	qtgetset(os,*qtptr);
736	<	if(!inset(os,id))
737	<	{
738	<	#ifdef DEBUG
739	<	eputs("qtRemove_tri(): tri not in set\n");
740	<	#endif
741	<	}
742	<	else
743	<	{
744	<	*qtptr = qtdelelem(*qtptr,id);
745	<	}
746	<	}
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	<	return(TRUE);
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
1199
1200

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines