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.2 by gwlarson, Thu Aug 20 16:47:21 1998 UTC vs.
Revision 3.12 by gwlarson, Fri Mar 5 16:32:49 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= 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 */
#	Line 41 | Line 48 | qtAlloc()                      /* allocate a quadtree */
48		if (QT_BLOCK(freet) >= QT_MAX_BLK)
49		return(EMPTY);
50		if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
51	<	(unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
52	<	return(EMPTY);
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((char *)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		{
#	Line 69 | 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	<	free((char *)quad_block[i],
101	<	(unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE));
100	>	if (quad_block[i] == NULL)
101	>	break;
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
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;
88	<	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	<	QUADTREE
178	<	qtLocate_leaf(qtptr,bcoord,type,which)
179	<	QUADTREE *qtptr;
180	<	double bcoord[3];
181	<	char *type,*which;
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;
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	<	return(qtLocate_leaf(child,bcoord,type,which));
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	>	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	>	}
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
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	<	qtRoot_point_locate(qtptr,v0,v1,v2,pt,type,which)
266	<	QUADTREE *qtptr;
267	<	FVECT v0,v1,v2;
268	<	FVECT pt;
269	<	char *type,*which;
270	<	{
271	<	char d,w;
272	<	int i,x,y;
273	<	QUADTREE *child;
274	<	QUADTREE qt;
275	<	FVECT n,i_pt;
276	<	double pd,bcoord[3];
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	>	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	<	/* Determine if point lies within pyramid (and therefore
283	<	inside a spherical quadtree cell):GT_INTERIOR, on one of the
284	<	pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
143	<	or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
144	<	For each triangle edge: compare the
145	<	point against the plane formed by the edge and the view center
146	<	*/
147	<	d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);
282	>	/* First check if can trivial accept: if vertex in cell */
283	>	if(s0 & s1 & s2)
284	>	goto Lfunc_call;
285
286	<	/* Not in this triangle */
287	<	if(!d)
288	<	{
289	<	if(which)
290	<	*which = 0;
291	<	return(EMPTY);
292	<	}
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	>	/* 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	>	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
364	<	/* Will return lowest level triangle containing point: It the
365	<	point is on an edge or vertex: will return first associated
366	<	triangle encountered in the child traversal- the others can
367	<	be derived using triangle adjacency information
368	<	*/
369	<	if(QT_IS_TREE(*qtptr))
370	<	{
371	<	/* Find the intersection point */
372	<	tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
373	<	intersect_vector_plane(pt,n,pd,NULL,i_pt);
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	<	i = max_index(n);
169	<	x = (i+1)%3;
170	<	y = (i+2)%3;
171	<	/* Calculate barycentric coordinates of i_pt */
172	<	bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
382	>	return(qt);
383
384	<	i = bary2d_child(bcoord);
385	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
386	<	return(qtLocate_leaf(child,bcoord,type,which));
387	<	}
388	<	else
179	<	if(!QT_IS_EMPTY(*qtptr))
180	<	{
181	<	/* map GT_VERTEX,GT_EDGE,GT_FACE GT_INTERIOR of pyramid to
182	<	spherical triangle primitives
183	<	*/
184	<	if(type)
185	<	*type = d;
186	<	if(which)
187	<	*which = w;
188	<	return(*qtptr);
189	<	}
190	<	return(EMPTY);
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	<	int
392	<	qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
393	<	QUADTREE *qtptr;
394	<	FVECT v0,v1,v2;
395	<	FVECT pt;
396	<	char *type,*which;
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	>	{
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	>	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	>	return(qt);
521	>	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	<	QUADTREE qt;
559	<	OBJECT os[MAXSET+1],*optr;
560	<	char d,w;
203	<	int i,id;
204	<	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	<	qt = qtRoot_point_locate(qtptr,v0,v1,v2,pt,type,which);
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(QT_IS_EMPTY(qt))
573	<	return(EMPTY);
574	<
211	<	/* Get the set */
212	<	qtgetset(os,qt);
213	<	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;
223	<	if(which)
224	<	*which = w;
225	<	return(id);
226	<	}
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
264	<	a
265	<	*/
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;
279	<	int i;
280	<	FVECT r0,r1,r2;
281	<	{
282	<	switch(i){
283	<	case 0:
284	<	VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
285	<	break;
286	<	case 1:
287	<	VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
288	<	break;
289	<	case 2:
290	<	VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
291	<	break;
292	<	case 3:
293	<	VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
294	<	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
737	<	/* Add triangle "id" to all leaf level cells that are children of
738	<	quadtree pointed to by "qtptr" with cell vertices "t1,t2,t3"
739	<	that it overlaps (vertex and edge adjacencies do not count
740	<	as overlapping). If the addition of the triangle causes the cell to go over
741	<	threshold- the cell is split- and the triangle must be recursively inserted
742	<	into the new child cells: it is assumed that "v1,v2,v3" are normalized
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	<	int
747	<	qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
748	<	QUADTREE *qtptr;
749	<	int id;
750	<	FVECT t0,t1,t2;
751	<	FVECT v0,v1,v2;
746	>	qtVisit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
747	>	s0,s1,s2,sq0,sq1,sq2,f,n)
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		int n;
757		{
758	<
759	<	char test;
760	<	int i,index;
317	<	FVECT a,b,c;
318	<	OBJECT os[MAXSET+1],*optr;
319	<	QUADTREE qt;
320	<	int found;
321	<	FVECT r0,r1,r2;
758	>	BCOORD a,b,c;
759	>	BCOORD va[3],vb[3],vc[3];
760	>	unsigned int sa,sb,sc;
761
762	<	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
763	<	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
762	>	if(n == 0)
763	>	return;
764	>	/* If a tree: trivial reject and recurse on potential children */
765	>	if(QT_IS_TREE(qt))
766	>	{
767	>	QT_SET_FLAG(qt);
768
769	<	/* If triangles do not overlap: done */
770	<	if(!test)
771	<	return(FALSE);
772	<	found = 0;
769	>	/* Test against new a edge: if entirely to left: only need
770	>	to visit child 0
771	>	*/
772	>	a = q1[0] + scale;
773	>	b = q0[1] + scale;
774	>	c = q0[2] + scale;
775
776	<	/* if this is tree: recurse */
332	<	if(QT_IS_TREE(*qtptr))
333	<	{
334	<	n++;
335	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
336	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c,n);
337	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b,n);
338	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2,n);
339	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a,n);
776	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
777
778	<	#if 0
342	<	if(!found)
778	>	if(sa==7)
779		{
780	<	#ifdef TEST_DRIVER
781	<	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
782	<	#else
783	<	eputs("qtAdd_tri():Found in parent but not children\n");
784	<	#endif
780	>	vb[1] = q0[1];
781	>	vc[2] = q0[2];
782	>	vc[1] = b;
783	>	vb[2] = c;
784	>	vb[0] = vc[0] = a;
785	>	qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
786	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
787	>	return;
788		}
789	<	#endif
789	>	if(sb==7)
790	>	{
791	>	va[0] = q1[0];
792	>	vc[2] = q0[2];
793	>	va[1] = vc[1] = b;
794	>	va[2] = c;
795	>	vc[0] = a;
796	>	qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,
797	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
798	>	return;
799	>	}
800	>	if(sc==7)
801	>	{
802	>	va[0] = q1[0];
803	>	vb[1] = q0[1];
804	>	va[1] = b;
805	>	va[2] = vb[2] = c;
806	>	vb[0] = a;
807	>	qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,
808	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
809	>	return;
810	>	}
811	>
812	>	va[0] = q1[0];
813	>	vb[1] = q0[1];
814	>	vc[2] = q0[2];
815	>	va[1] = vc[1] = b;
816	>	va[2] = vb[2] = c;
817	>	vb[0] = vc[0] = a;
818	>	/* If outside of all edges: only need to Visit 3  */
819	>	if(sa==0 && sb==0 && sc==0)
820	>	{
821	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,
822	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
823	>	return;
824	>	}
825	>
826	>	if(sa)
827	>	qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
828	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
829	>	if(sb)
830	>	qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
831	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
832	>	if(sc)
833	>	qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
834	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
835	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
836	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
837		}
838	+	/* Leaf node: Do definitive test */
839		else
840	<	{
841	<	/* If this leave node emptry- create a new set */
842	<	if(QT_IS_EMPTY(*qtptr))
356	<	{
357	<	*qtptr = qtaddelem(*qtptr,id);
358	<	#if 0
359	<	{
360	<	int k;
361	<	qtgetset(os,*qtptr);
362	<	printf("\n%d:\n",os[0]);
363	<	for(k=1; k <= os[0];k++)
364	<	printf("%d ",os[k]);
365	<	printf("\n");
366	<	}
367	<	#endif
368	<	/*
369	<	os[0] = 0;
370	<	insertelem(os,id);
371	<	qt = fullnode(os);
372	<	*qtptr = qt;
373	<	*/
374	<	}
375	<	else
376	<	{
377	<	qtgetset(os,*qtptr);
378	<	/* If the set is too large: subdivide */
379	<	if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
380	<	{
381	<	*qtptr = qtaddelem(*qtptr,id);
382	<	#if 0
383	<	{
384	<	int k;
385	<	qtgetset(os,*qtptr);
386	<	printf("\n%d:\n",os[0]);
387	<	for(k=1; k <= os[0];k++)
388	<	printf("%d ",os[k]);
389	<	printf("\n");
390	<	}
391	<	#endif
392	<	/*
393	<	insertelem(os,id);
394	<	*qtptr = fullnode(os);
395	<	*/
396	<	}
397	<	else
398	<	{
399	<	if (n < QT_MAX_LEVELS)
400	<	{
401	<	/* If set size exceeds threshold: subdivide cell and
402	<	reinsert set tris into cell
403	<	*/
404	<	n++;
405	<	qtfreeleaf(*qtptr);
406	<	qtSubdivide(qtptr);
407	<	found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
408	<	#if 0
409	<	if(!found)
410	<	{
411	<	#ifdef TEST_DRIVER
412	<	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
413	<	#else
414	<	eputs("qtAdd_tri():Found in parent but not children\n");
415	<	#endif
416	<	}
417	<	#endif
418	<	for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
419	<	{
420	<	id = QT_SET_NEXT_ELEM(optr);
421	<	qtTri_verts_from_id(id,r0,r1,r2);
422	<	found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,n);
423	<	#ifdef DEBUG
424	<	if(!found)
425	<	eputs("qtAdd_tri():Reinsert-in parent but not children\n");
426	<	#endif
427	<	}
428	<	}
429	<	else
430	<	if(QT_SET_CNT(os) < QT_MAX_SET)
431	<	{
432	<	*qtptr = qtaddelem(*qtptr,id);
433	<	#if 0
434	<	{
435	<	int k;
436	<	qtgetset(os,*qtptr);
437	<	printf("\n%d:\n",os[0]);
438	<	for(k=1; k <= os[0];k++)
439	<	printf("%d ",os[k]);
440	<	printf("\n");
441	<	}
442	<	#endif
443	<	/*
444	<	insertelem(os,id);
445	<	*qtptr = fullnode(os);
446	<	*/
447	<	}
448	<	else
449	<	{
450	<	#ifdef DEBUG
451	<	eputs("qtAdd_tri():two many levels\n");
452	<	#endif
453	<	return(FALSE);
454	<	}
455	<	}
456	<	}
457	<	}
458	<	return(TRUE);
840	>	qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
841	>	scale,s0,s1,s2,sq0,sq1,sq2,f,n);
842	>
843		}
844
845
846	<	int
847	<	qtApply_to_tri_cells(qtptr,t0,t1,t2,v0,v1,v2,func,arg)
848	<	QUADTREE *qtptr;
849	<	FVECT t0,t1,t2;
850	<	FVECT v0,v1,v2;
851	<	int (*func)();
852	<	char *arg;
846	>	qtVisit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
847	>	s0,s1,s2,sq0,sq1,sq2,f,n)
848	>	int root;
849	>	QUADTREE qt;
850	>	BCOORD q0[3],q1[3],q2[3];
851	>	BCOORD t0[3],t1[3],t2[3];
852	>	BCOORD dt10[3],dt21[3],dt02[3];
853	>	BCOORD scale;
854	>	unsigned int s0,s1,s2,sq0,sq1,sq2;
855	>	FUNC f;
856	>	int n;
857		{
858	<	char test;
859	<	FVECT a,b,c;
858	>	BCOORD a,b,c;
859	>	BCOORD va[3],vb[3],vc[3];
860	>	unsigned int sa,sb,sc;
861
862	<	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
863	<	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
862	>	if(n==0)
863	>	return;
864	>	/* If a tree: trivial reject and recurse on potential children */
865	>	if(QT_IS_TREE(qt))
866	>	{
867	>	QT_SET_FLAG(qt);
868	>	/* Test against new a edge: if entirely to left: only need
869	>	to visit child 0
870	>	*/
871	>	a = q1[0] - scale;
872	>	b = q0[1] - scale;
873	>	c = q0[2] - scale;
874
875	<	/* If triangles do not overlap: done */
477	<	if(!test)
478	<	return(FALSE);
875	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
876
877	<	/* if this is tree: recurse */
878	<	if(QT_IS_TREE(*qtptr))
879	<	{
880	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
881	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t0,t1,t2,v0,a,c,func,arg);
882	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t0,t1,t2,a,v1,b,func,arg);
883	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t0,t1,t2,c,b,v2,func,arg);
884	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t0,t1,t2,b,c,a,func,arg);
877	>	if(sa==0)
878	>	{
879	>	vb[1] = q0[1];
880	>	vc[2] = q0[2];
881	>	vc[1] = b;
882	>	vb[2] = c;
883	>	vb[0] = vc[0] = a;
884	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
885	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
886	>	return;
887	>	}
888	>	if(sb==0)
889	>	{
890	>	va[0] = q1[0];
891	>	vc[2] = q0[2];
892	>	va[1] = vc[1] = b;
893	>	va[2] = c;
894	>	vc[0] = a;
895	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
896	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
897	>	return;
898	>	}
899	>	if(sc==0)
900	>	{
901	>	va[0] = q1[0];
902	>	vb[1] = q0[1];
903	>	va[1] = b;
904	>	va[2] = vb[2] = c;
905	>	vb[0] = a;
906	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
907	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
908	>	return;
909	>	}
910	>	va[0] = q1[0];
911	>	vb[1] = q0[1];
912	>	vc[2] = q0[2];
913	>	va[1] = vc[1] = b;
914	>	va[2] = vb[2] = c;
915	>	vb[0] = vc[0] = a;
916	>	/* If outside of all edges: only need to Visit 3  */
917	>	if(sa==7 && sb==7 && sc==7)
918	>	{
919	>	qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,
920	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
921	>	return;
922	>	}
923	>	if(sa!=7)
924	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
925	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
926	>	if(sb!=7)
927	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
928	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
929	>	if(sc!=7)
930	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
931	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
932	>
933	>	qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
934	>	t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
935	>	return;
936		}
937	+	/* Leaf node: Do definitive test */
938		else
939	<	return(func(qtptr,arg));
939	>	qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
940	>	scale,s0,s1,s2,sq0,sq1,sq2,f,n);
941		}
942
943
944	<	int
945	<	qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
946	<	QUADTREE *qtptr;
947	<	int id;
948	<	FVECT t0,t1,t2;
949	<	FVECT v0,v1,v2;
950	<	{
944	>
945	>	qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
946	>	scale, s0,s1,s2,sq0,sq1,sq2,f,n)
947	>	int root;
948	>	QUADTREE qt;
949	>	BCOORD q0[3],q1[3],q2[3];
950	>	BCOORD t0[3],t1[3],t2[3];
951	>	BCOORD dt10[3],dt21[3],dt02[3];
952	>	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
953	>	FUNC f;
954	>	int n;
955	>	{
956	>	double db;
957	>	unsigned int e0,e1,e2;
958	>	char al,ah,bl,bh,cl,ch,testl,testh;
959	>	char test_t0t1,test_t1t2,test_t2t0;
960	>	BCOORD a,b,c;
961
962	<	char test;
963	<	int i;
964	<	FVECT a,b,c;
505	<	OBJECT os[MAXSET+1];
506	<
507	<	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
508	<	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
962	>	/* First check if can trivial accept: if vertex in cell */
963	>	if(s0 & s1 & s2)
964	>	goto Lfunc_call;
965
966	<	/* If triangles do not overlap: done */
967	<	if(!test)
968	<	return(FALSE);
969	<
970	<	/* if this is tree: recurse */
971	<	if(QT_IS_TREE(*qtptr))
966	>	/* Assumption: Could not trivial reject*/
967	>	/* IF any coords lie on edge- must be in if couldnt reject */
968	>	a = q1[0];b= q0[1]; c= q0[2];
969	>	if(t0[0]== a || t0[1] == b || t0[2] == c)
970	>	goto Lfunc_call;
971	>	if(t1[0]== a || t1[1] == b || t1[2] == c)
972	>	goto Lfunc_call;
973	>	if(t2[0]== a || t2[1] == b || t2[2] == c)
974	>	goto Lfunc_call;
975	>
976	>	/* Test for edge crossings */
977	>	/* Test t0t1,t1t2,t2t0 against a */
978	>	/* Calculate edge crossings */
979	>	e0  = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
980	>	/* See if edge can be trivially rejected from intersetion testing */
981	>	test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
982	>	((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
983	>	bl=bh=0;
984	>	if(test_t0t1 && (e0 & 2) )
985		{
986	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
987	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c);
988	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b);
520	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2);
521	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a);
986	>	/* Must do double calculation to avoid overflow */
987	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
988	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
989		}
990	<	else
990	>	test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
991	>	((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
992	>	if(test_t1t2 && (e0 & 1))
993	>	{    /* test t1t2 against a */
994	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
995	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
996	>	}
997	>	test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
998	>	((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
999	>	if(test_t2t0 && (e0 & 4))
1000		{
1001	<	if(QT_IS_EMPTY(*qtptr))
1002	<	{
1003	<	#ifdef DEBUG
1004	<	eputs("qtRemove_tri(): triangle not found\n");
1005	<	#endif
1006	<	}
1007	<	/* remove id from set */
1008	<	else
1009	<	{
1010	<	qtgetset(os,*qtptr);
1011	<	if(!inset(os,id))
1012	<	{
1013	<	#ifdef DEBUG
1014	<	eputs("qtRemove_tri(): tri not in set\n");
1015	<	#endif
1016	<	}
1017	<	else
1018	<	{
1019	<	*qtptr = qtdelelem(*qtptr,id);
1020	<	#if 0
1021	<	{
1022	<	int k;
1023	<	if(!QT_IS_EMPTY(*qtptr))
1024	<	{qtgetset(os,*qtptr);
1025	<	printf("\n%d:\n",os[0]);
1026	<	for(k=1; k <= os[0];k++)
1027	<	printf("%d ",os[k]);
1028	<	printf("\n");
1029	<	}
1001	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1002	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1003	>	}
1004	>	e1  = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
1005	>	cl = ch = 0;
1006	>	if(test_t0t1 && (e1 & 2))
1007	>	{/* test t0t1 against b */
1008	>	db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
1009	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1010	>	}
1011	>	if(test_t1t2 && (e1 & 1))
1012	>	{/* test t1t2 against b */
1013	>	db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
1014	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1015	>	}
1016	>	if(test_t2t0 && (e1 & 4))
1017	>	{/* test t2t0 against b */
1018	>	db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
1019	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1020	>	}
1021	>
1022	>	/* test edges against c */
1023	>	e2  = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
1024	>	al = ah = 0;
1025	>	if(test_t0t1 && (e2 & 2))
1026	>	{ /* test t0t1 against c */
1027	>	db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
1028	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1029	>	}
1030	>	if(test_t1t2 && (e2 & 1))
1031	>	{
1032	>	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1033	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1034	>	}
1035	>	if(test_t2t0 && (e2 & 4))
1036	>	{ /* test t2t0 against c */
1037	>	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1038	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1039	>	}
1040	>	/* Only remaining case is if some edge trivially rejected */
1041	>	if(!e0 || !e1 || !e2)
1042	>	return;
1043
1044	<	}
1045	<	#endif
1046	<	}
1047	<	}
1044	>	/* Only one/none got tested - pick either of other two/any two */
1045	>	/* Only need to test one edge */
1046	>	if(!test_t0t1 && (e0 & 2))
1047	>	{
1048	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1049	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1050		}
1051	<	return(TRUE);
1051	>	if(!test_t1t2 && (e0 & 1))
1052	>	{/* test t1t2 against a */
1053	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1054	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1055	>	}
1056	>	if(!test_t2t0 && (e0 & 4))
1057	>	{
1058	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1059	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1060	>	}
1061	>
1062	>	return;
1063	>
1064	>	Lfunc_call:
1065	>
1066	>	f.func(f.argptr,root,qt,n);
1067	>	if(!QT_IS_EMPTY(qt))
1068	>	QT_LEAF_SET_FLAG(qt);
1069	>
1070		}
1071
1072
1073
1074	+	/* Leaf node: Do definitive test */
1075	+	QUADTREE
1076	+	qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1077	+	scale, s0,s1,s2,sq0,sq1,sq2,f,n)
1078	+	int root;
1079	+	QUADTREE qt;
1080	+	BCOORD q0[3],q1[3],q2[3];
1081	+	BCOORD t0[3],t1[3],t2[3];
1082	+	BCOORD dt10[3],dt21[3],dt02[3];
1083	+	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
1084	+	FUNC f;
1085	+	int n;
1086	+	{
1087	+	double db;
1088	+	unsigned int e0,e1,e2;
1089	+	BCOORD a,b,c;
1090	+	double p0[2], p1[2],cp;
1091	+	char al,ah,bl,bh,cl,ch;
1092	+	char testl,testh,test_t0t1,test_t1t2,test_t2t0;
1093	+	unsigned int ls0,ls1,ls2;
1094	+
1095	+	/* May have gotten passed trivial reject if one/two vertices are ON */
1096	+	a = q1[0];b= q0[1]; c= q0[2];
1097	+	SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
1098	+
1099	+	/* First check if can trivial accept: if vertex in cell */
1100	+	if(ls0 & ls1 & ls2)
1101	+	goto Lfunc_call;
1102
1103	+	if(ls0==0 || ls1 == 0 || ls2 ==0)
1104	+	return;
1105	+	/* Assumption: Could not trivial reject*/
1106	+	/* IF any coords lie on edge- must be in if couldnt reject */
1107
1108	+	if(t0[0]== a || t0[1] == b || t0[2] == c)
1109	+	goto Lfunc_call;
1110	+	if(t1[0]== a || t1[1] == b || t1[2] == c)
1111	+	goto Lfunc_call;
1112	+	if(t2[0]== a || t2[1] == b || t2[2] == c)
1113	+	goto Lfunc_call;
1114
1115	+	/* Test for edge crossings */
1116	+	/* Test t0t1 against a,b,c */
1117	+	/* First test if t0t1 can be trivially rejected */
1118	+	/* If both edges are outside bounds- dont need to test */
1119	+
1120	+	/* Test t0t1,t1t2,t2t0 against a */
1121	+	test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
1122	+	((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
1123	+	e0  = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
1124	+	bl=bh=0;
1125	+	/* Test t0t1,t1t2,t2t0 against a */
1126	+	if(test_t0t1 && (e0 & 2) )
1127	+	{
1128	+	db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
1129	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1130	+	}
1131	+	test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
1132	+	((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
1133	+	if(test_t1t2 && (e0 & 1))
1134	+	{    /* test t1t2 against a */
1135	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1136	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1137	+	}
1138	+	test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
1139	+	((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
1140	+	if(test_t2t0 && (e0 & 4))
1141	+	{
1142	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1143	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1144	+	}
1145	+	cl = ch = 0;
1146	+	e1  = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
1147	+	if(test_t0t1 && (e1 & 2))
1148	+	{/* test t0t1 against b */
1149	+	db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
1150	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1151	+	}
1152	+	if(test_t1t2 && (e1 & 1))
1153	+	{/* test t1t2 against b */
1154	+	db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
1155	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1156	+	}
1157	+	if(test_t2t0 && (e1 & 4))
1158	+	{/* test t2t0 against b */
1159	+	db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
1160	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1161	+	}
1162	+	al = ah = 0;
1163	+	e2  = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
1164	+	if(test_t0t1 && (e2 & 2))
1165	+	{ /* test t0t1 against c */
1166	+	db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
1167	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1168	+	}
1169	+	if(test_t1t2 && (e2 & 1))
1170	+	{
1171	+	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1172	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1173	+	}
1174	+	if(test_t2t0 && (e2 & 4))
1175	+	{ /* test t2t0 against c */
1176	+	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1177	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1178	+	}
1179	+	/* Only remaining case is if some edge trivially rejected */
1180	+	if(!e0 || !e1 || !e2)
1181	+	return;
1182
1183	+	/* Only one/none got tested - pick either of other two/any two */
1184	+	/* Only need to test one edge */
1185	+	if(!test_t0t1 && (e0 & 2))
1186	+	{
1187	+	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1188	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1189	+	}
1190	+	if(!test_t1t2 && (e0 & 1))
1191	+	{/* test t1t2 against a */
1192	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1193	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1194	+	}
1195	+	if(!test_t2t0 && (e0 & 4))
1196	+	{
1197	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1198	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1199	+	}
1200
1201	+	return;
1202	+
1203	+	Lfunc_call:
1204	+	f.func(f.argptr,root,qt,n);
1205	+	if(!QT_IS_EMPTY(qt))
1206	+	QT_LEAF_SET_FLAG(qt);
1207	+	}
1208
1209
1210

Diff Legend

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