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.4 by gwlarson, Fri Sep 11 11:52:26 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
#	Line 27 | Line 27 | int4 *quad_flag= NULL;
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 46 | 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));
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 60 | Line 64 | qtAlloc()                      /* allocate a quadtree */
64
65		qtClearAllFlags()               /* clear all quadtree branch flags */
66		{
67	<	if (!treetop) return;
68	<	bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*(QT_BLOCK_SIZE/8));
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
67	–
76		qtFree(qt)                      /* free a quadtree */
77		register QUADTREE  qt;
78		{
#	Line 85 | 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		{
#	Line 94 | 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;
#	Line 101 | Line 110 | qtDone()                       /* free EVERYTHING */
110		}
111
112		QUADTREE
113	<	*qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
105	<	QUADTREE *qtptr;
106	<	double bcoord[3];
107	<	FVECT t0,t1,t2;
108	<	{
109	<	int i;
110	<	QUADTREE *child;
111	<	FVECT a,b,c;
112	<
113	<	if(QT_IS_TREE(*qtptr))
114	<	{
115	<	i = bary_child(bcoord);
116	<	#ifdef DEBUG_TEST_DRIVER
117	<	qtSubdivide_tri(Pick_v0[Pick_cnt-1],Pick_v1[Pick_cnt-1],
118	<	Pick_v2[Pick_cnt-1],a,b,c);
119	<	qtNth_child_tri(Pick_v0[Pick_cnt-1],Pick_v1[Pick_cnt-1],
120	<	Pick_v2[Pick_cnt-1],a,b,c,i,
121	<	Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
122	<	Pick_v2[Pick_cnt]);
123	<	Pick_cnt++;
124	<	#endif
125	<
126	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
127	<	if(t0)
128	<	{
129	<	qtSubdivide_tri(t0,t1,t2,a,b,c);
130	<	qtNth_child_tri(t0,t1,t2,a,b,c,i,t0,t1,t2);
131	<	}
132	<	return(qtLocate_leaf(child,bcoord,t0,t1,t2));
133	<	}
134	<	else
135	<	return(qtptr);
136	<	}
137	<
138	<
139	<	QUADTREE
140	<	*qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
141	<	QUADTREE *qtptr;
142	<	FVECT v0,v1,v2;
143	<	FVECT pt;
144	<	FVECT t0,t1,t2;
145	<	{
146	<	int d;
147	<	int i,x,y;
148	<	QUADTREE *child;
149	<	FVECT n,i_pt,a,b,c;
150	<	double pd,bcoord[3];
151	<
152	<	/* Determine if point lies within pyramid (and therefore
153	<	inside a spherical quadtree cell):GT_INTERIOR, on one of the
154	<	pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
155	<	or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
156	<	For each triangle edge: compare the
157	<	point against the plane formed by the edge and the view center
158	<	*/
159	<	d = point_in_stri(v0,v1,v2,pt);
160	<
161	<
162	<	/* Not in this triangle */
163	<	if(!d)
164	<	return(NULL);
165	<
166	<	/* Will return lowest level triangle containing point: It the
167	<	point is on an edge or vertex: will return first associated
168	<	triangle encountered in the child traversal- the others can
169	<	be derived using triangle adjacency information
170	<	*/
171	<	if(QT_IS_TREE(*qtptr))
172	<	{
173	<	/* Find the intersection point */
174	<	tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
175	<	intersect_vector_plane(pt,n,pd,NULL,i_pt);
176	<
177	<	i = max_index(n);
178	<	x = (i+1)%3;
179	<	y = (i+2)%3;
180	<	/* Calculate barycentric coordinates of i_pt */
181	<	bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
182	<
183	<	i = bary_child(bcoord);
184	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
185	<	#ifdef DEBUG_TEST_DRIVER
186	<	Pick_cnt = 0;
187	<	VCOPY(Pick_v0[0],v0);
188	<	VCOPY(Pick_v1[0],v1);
189	<	VCOPY(Pick_v2[0],v2);
190	<	Pick_cnt++;
191	<	qtSubdivide_tri(Pick_v0[Pick_cnt-1],Pick_v1[Pick_cnt-1],
192	<	Pick_v2[Pick_cnt-1],a,b,c);
193	<	qtNth_child_tri(Pick_v0[Pick_cnt-1],Pick_v1[Pick_cnt-1],
194	<	Pick_v2[Pick_cnt-1],a,b,c,i,
195	<	Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
196	<	Pick_v2[Pick_cnt]);
197	<	Pick_cnt++;
198	<	#endif
199	<	if(t0)
200	<	{
201	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
202	<	qtNth_child_tri(v0,v1,v2,a,b,c,i,t0,t1,t2);
203	<	}
204	<	return(qtLocate_leaf(child,bcoord,t0,t1,t2));
205	<	}
206	<	else
207	<	{
208	<	if(t0)
209	<	{
210	<	VCOPY(t0,v0);
211	<	VCOPY(t1,v1);
212	<	VCOPY(t2,v2);
213	<	}
214	<	return(qtptr);
215	<	}
216	<	}
217	<
218	<
219	<	QUADTREE
220	<	qtSubdivide(qtptr)
221	<	QUADTREE *qtptr;
222	<	{
223	<	QUADTREE node;
224	<	node = qtAlloc();
225	<	QT_CLEAR_CHILDREN(node);
226	<	*qtptr = node;
227	<	return(node);
228	<	}
229	<
230	<
231	<	QUADTREE
232	<	qtSubdivide_nth_child(qt,n)
113	>	qtLocate(qt,bcoord)
114		QUADTREE qt;
115	<	int n;
115	>	BCOORD bcoord[3];
116		{
236	–	QUADTREE node;
237	–
238	–	node = qtSubdivide(&(QT_NTH_CHILD(qt,n)));
239	–
240	–	return(node);
241	–	}
242	–
243	–	/* for triangle v0-v1-v2- returns a,b,c: children are:
244	–
245	–	v2                     0: v0,a,c
246	–	/\                     1: a,v1,b
247	–	/2 \                    2: c,b,v2
248	–	c/____\b                  3: b,c,a
249	–	/\    /\
250	–	/0 \3 /1 \
251	–	v0____\/____\v1
252	–	a
253	–	*/
254	–
255	–	qtSubdivide_tri(v0,v1,v2,a,b,c)
256	–	FVECT v0,v1,v2;
257	–	FVECT a,b,c;
258	–	{
259	–	EDGE_MIDPOINT_VEC3(a,v0,v1);
260	–	EDGE_MIDPOINT_VEC3(b,v1,v2);
261	–	EDGE_MIDPOINT_VEC3(c,v2,v0);
262	–	}
263	–
264	–	qtNth_child_tri(v0,v1,v2,a,b,c,i,r0,r1,r2)
265	–	FVECT v0,v1,v2;
266	–	FVECT a,b,c;
267	–	int i;
268	–	FVECT r0,r1,r2;
269	–	{
270	–
271	–	switch(i){
272	–	case 0:
273	–	VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
274	–	break;
275	–	case 1:
276	–	VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
277	–	break;
278	–	case 2:
279	–	VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
280	–	break;
281	–	case 3:
282	–	VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
283	–	break;
284	–	}
285	–	}
286	–
287	–	/* Add triangle "id" to all leaf level cells that are children of
288	–	quadtree pointed to by "qtptr" with cell vertices "t1,t2,t3"
289	–	that it overlaps (vertex and edge adjacencies do not count
290	–	as overlapping). If the addition of the triangle causes the cell to go over
291	–	threshold- the cell is split- and the triangle must be recursively inserted
292	–	into the new child cells: it is assumed that "v1,v2,v3" are normalized
293	–	*/
294	–
295	–	int
296	–	qtRoot_add_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
297	–	QUADTREE *qtptr;
298	–	FVECT q0,q1,q2;
299	–	FVECT t0,t1,t2;
300	–	int id;
301	–	int n;
302	–	{
303	–	int test;
304	–	int found;
305	–
306	–	test = stri_intersect(q0,q1,q2,t0,t1,t2);
307	–	if(!test)
308	–	return(FALSE);
309	–
310	–	found = qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n);
311	–
312	–	return(found);
313	–	}
314	–
315	–	int
316	–	qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
317	–	QUADTREE *qtptr;
318	–	FVECT q0,q1,q2;
319	–	FVECT t0,t1,t2;
320	–	int id;
321	–	int n;
322	–	{
323	–	int i,index,test,found;
324	–	FVECT a,b,c;
325	–	OBJECT os[QT_MAXSET+1],*optr;
326	–	QUADTREE qt;
327	–	FVECT r0,r1,r2;
328	–
329	–	found = 0;
330	–	/* if this is tree: recurse */
331	–	if(QT_IS_TREE(*qtptr))
332	–	{
333	–	n++;
334	–	qtSubdivide_tri(q0,q1,q2,a,b,c);
335	–	test = stri_intersect(t0,t1,t2,q0,a,c);
336	–	if(test)
337	–	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),q0,a,c,t0,t1,t2,id,n);
338	–	test = stri_intersect(t0,t1,t2,a,q1,b);
339	–	if(test)
340	–	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),a,q1,b,t0,t1,t2,id,n);
341	–	test = stri_intersect(t0,t1,t2,c,b,q2);
342	–	if(test)
343	–	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),c,b,q2,t0,t1,t2,id,n);
344	–	test = stri_intersect(t0,t1,t2,b,c,a);
345	–	if(test)
346	–	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),b,c,a,t0,t1,t2,id,n);
347	–	}
348	–	else
349	–	{
350	–	/* If this leave node emptry- create a new set */
351	–	if(QT_IS_EMPTY(*qtptr))
352	–	*qtptr = qtaddelem(*qtptr,id);
353	–	else
354	–	{
355	–	/* If the set is too large: subdivide */
356	–	optr = qtqueryset(*qtptr);
357	–
358	–	if(QT_SET_CNT(optr) < QT_SET_THRESHOLD)
359	–	*qtptr = qtaddelem(*qtptr,id);
360	–	else
361	–	{
362	–	if (n < QT_MAX_LEVELS)
363	–	{
364	–	/* If set size exceeds threshold: subdivide cell and
365	–	reinsert set tris into cell
366	–	*/
367	–	qtgetset(os,*qtptr);
368	–
369	–	n++;
370	–	qtfreeleaf(*qtptr);
371	–	qtSubdivide(qtptr);
372	–	found = qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n);
373	–
374	–	for(optr = QT_SET_PTR(os),i = QT_SET_CNT(os); i > 0; i--)
375	–	{
376	–	id = QT_SET_NEXT_ELEM(optr);
377	–	qtTri_from_id(id,NULL,NULL,NULL,r0,r1,r2,NULL,NULL,NULL);
378	–	found=qtAdd_tri(qtptr,q0,q1,q2,r0,r1,r2,id,n);
379	–	#ifdef DEBUG
380	–	if(!found)
381	–	eputs("qtAdd_tri():Reinsert-in parent but not children\n");
382	–	#endif
383	–	}
384	–	}
385	–	else
386	–	if(QT_SET_CNT(optr) < QT_MAXSET)
387	–	*qtptr = qtaddelem(*qtptr,id);
388	–	else
389	–	{
390	–	#ifdef DEBUG
391	–	eputs("qtAdd_tri():two many levels\n");
392	–	#endif
393	–	return(FALSE);
394	–	}
395	–	}
396	–	}
397	–	}
398	–	return(TRUE);
399	–	}
400	–
401	–
402	–	int
403	–	qtApply_to_tri_cells(qtptr,t0,t1,t2,v0,v1,v2,func,arg)
404	–	QUADTREE *qtptr;
405	–	FVECT t0,t1,t2;
406	–	FVECT v0,v1,v2;
407	–	int (*func)();
408	–	int *arg;
409	–	{
410	–	int test;
411	–	FVECT a,b,c;
412	–
413	–	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
414	–	test = stri_intersect(t0,t1,t2,v0,v1,v2);
415	–
416	–	/* If triangles do not overlap: done */
417	–	if(!test)
418	–	return(FALSE);
419	–
420	–	/* if this is tree: recurse */
421	–	func(qtptr,arg);
422	–
423	–	if(QT_IS_TREE(*qtptr))
424	–	{
425	–	QT_SET_FLAG(*qtptr);
426	–	qtSubdivide_tri(v0,v1,v2,a,b,c);
427	–	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t0,t1,t2,v0,a,c,func,arg);
428	–	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t0,t1,t2,a,v1,b,func,arg);
429	–	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t0,t1,t2,c,b,v2,func,arg);
430	–	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t0,t1,t2,b,c,a,func,arg);
431	–	}
432	–	}
433	–
434	–	int
435	–	qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
436	–	QUADTREE *qtptr;
437	–	int id;
438	–	FVECT t0,t1,t2;
439	–	FVECT v0,v1,v2;
440	–	{
441	–
442	–	int test;
117		int i;
444	–	FVECT a,b,c;
445	–	OBJECT os[QT_MAXSET+1];
446	–
447	–	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
448	–	test = stri_intersect(t0,t1,t2,v0,v1,v2);
118
119	<	/* If triangles do not overlap: done */
120	<	if(!test)
121	<	return(FALSE);
119	>	if(QT_IS_TREE(qt))
120	>	{
121	>	i = bary_child(bcoord);
122
123	<	/* if this is tree: recurse */
455	<	if(QT_IS_TREE(*qtptr))
456	<	{
457	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
458	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c);
459	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b);
460	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2);
461	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a);
123	>	return(qtLocate(QT_NTH_CHILD(qt,i),bcoord));
124		}
125		else
126	<	{
465	<	if(QT_IS_EMPTY(*qtptr))
466	<	{
467	<	#ifdef DEBUG
468	<	eputs("qtRemove_tri(): triangle not found\n");
469	<	#endif
470	<	}
471	<	/* remove id from set */
472	<	else
473	<	{
474	<	if(!qtinset(*qtptr,id))
475	<	{
476	<	#ifdef DEBUG
477	<	eputs("qtRemove_tri(): tri not in set\n");
478	<	#endif
479	<	}
480	<	else
481	<	{
482	<	*qtptr = qtdelelem(*qtptr,id);
483	<	}
484	<	}
485	<	}
486	<	return(TRUE);
126	>	return(qt);
127		}
128
489	–
129		int
130		move_to_nbr(b,db0,db1,db2,tptr)
131	<	double b[3],db0,db1,db2;
131	>	BCOORD b[3];
132	>	BDIR db0,db1,db2;
133		double *tptr;
134		{
135		double t,dt;
136	+	BCOORD bc;
137		int nbr;
138	<
138	>
139		nbr = -1;
140	+	*tptr = 0.0;
141		/* Advance to next node */
142	<	if(!ZERO(db0) && db0 < 0.0)
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 = -b[0]/db0;
150	>	t = ((double)b[0])/-db0;
151		nbr = 0;
152		}
153		else
154	<	t = FHUGE;
155	<	if(!ZERO(db1) && db1 < 0.0 )
154	>	t = MAXFLOAT;
155	>	if(db1 < 0)
156		{
157	<	dt = -b[1]/db1;
157	>	dt = ((double)b[1])/-db1;
158		if( dt < t)
159		{
160		t = dt;
161		nbr = 1;
162		}
163		}
164	<	if(!ZERO(db2) && db2 < 0.0 )
165	<	{
166	<	dt = -b[2]/db2;
167	<	if( dt < t)
164	>	if(db2 < 0)
165	>	{
166	>	dt = ((double)b[2])/-db2;
167	>	if( dt < t)
168		{
169		t = dt;
170		nbr = 2;
#	Line 526 | Line 174 | double *tptr;
174		return(nbr);
175		}
176
177	<	int
178	<	qtTrace_edge(qtptr,b,db,orig,dir,max_t,func,arg1,arg2)
179	<	QUADTREE *qtptr;
180	<	double b[3],db[3];
181	<	FVECT orig,dir;
182	<	double max_t;
183	<	int (*func)();
184	<	int *arg1,arg2;
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		{
538	–
186		int i,found;
187	<	QUADTREE *child;
188	<	int nbr,next;
187	>	QUADTREE child;
188	>	int nbr,next,w;
189		double t;
543	–	#ifdef DEBUG_TEST_DRIVER
544	–
545	–	FVECT a1,b1,c1;
546	–	int Pick_parent = Pick_cnt-1;
547	–	qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
548	–	Pick_v2[Pick_parent],a1,b1,c1);
190
191	<	#endif
551	<	if(QT_IS_TREE(*qtptr))
191	>	if(QT_IS_TREE(qt))
192		{
193	<	/* Find the appropriate child and reset the coord */
194	<	i = bary_child(b);
193	>	/* Find the appropriate child and reset the coord */
194	>	i = bary_child(b);
195
196	<	QT_SET_FLAG(*qtptr);
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	<	for(;;)
206	<	{
207	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
208	<
209	<	if(i != 3)
210	<	{
211	<
212	<	db[0] *= 2.0;db[1] *= 2.0; db[2] *= 2.0;
213	<	nbr = qtTrace_edge(child,b,db,orig,dir,max_t,func,arg1,arg2);
214	<	db[0] *= 0.5;db[1] *= 0.5; db[2] *= 0.5;
215	<	}
216	<	else
217	<	{
218	<	db[0] *=-2.0;db[1] *= -2.0; db[2] *= -2.0;
572	<	/* If the center cell- must flip direction signs */
573	<	nbr = qtTrace_edge(child,b,db,orig,dir,max_t,func,arg1,arg2);
574	<	db[0] *=-0.5;db[1] *= -0.5; db[2] *= -0.5;
575	<	}
576	<	if(nbr == QT_DONE)
577	<	return(nbr);
578	<
579	<	/* If in same block: traverse */
580	<	if(i==3)
581	<	next = nbr;
582	<	else
583	<	if(nbr == i)
584	<	next = 3;
585	<	else
586	<	{
587	<	/* reset the barycentric coordinates in the parents*/
588	<	bary_parent(b,i);
589	<	/* Else pop up to parent and traverse from there */
590	<	return(nbr);
591	<	}
592	<	bary_from_child(b,i,next);
593	<	i = next;
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		{
226	<	#ifdef DEBUG_TEST_DRIVER
227	<	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
228	<	Pick_v2[Pick_parent],a1,b1,c1,i,
229	<	Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
230	<	Pick_v2[Pick_cnt]);
231	<	Pick_cnt++;
232	<	#endif
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(func(qtptr,arg1,arg2) == QT_DONE)
237	<	return(QT_DONE);
238	<
239	<	/* Advance to next node */
240	<	/* NOTE: Optimize: should only have to check 1/2 */
241	<	nbr = move_to_nbr(b,db[0],db[1],db[2],&t);
242	<
243	<	if(t >= max_t)
244	<	return(QT_DONE);
245	<	if(nbr != -1)
246	<	{
617	<	b[0] += t * db[0];
618	<	b[1] += t * db[1];
619	<	b[2] += t * db[2];
620	<	db[0] *= (1.0 - t);
621	<	db[1] *= (1.0 - t);
622	<	db[2] *= (1.0 - t);
623	<	}
624	<	return(nbr);
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	<
627	<	}
248	>	}
249
250
630	–	int
631	–	qtTrace_ray(qtptr,b,db0,db1,db2,orig,dir,func,arg1,arg2)
632	–	QUADTREE *qtptr;
633	–	double b[3],db0,db1,db2;
634	–	FVECT orig,dir;
635	–	int (*func)();
636	–	int *arg1,arg2;
637	–	{
251
639	–	int i,found;
640	–	QUADTREE *child;
641	–	int nbr,next;
642	–	double t;
643	–	#ifdef DEBUG_TEST_DRIVER
644	–
645	–	FVECT a1,b1,c1;
646	–	int Pick_parent = Pick_cnt-1;
647	–	qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
648	–	Pick_v2[Pick_parent],a1,b1,c1);
252
650	–	#endif
651	–	if(QT_IS_TREE(*qtptr))
652	–	{
653	–	/* Find the appropriate child and reset the coord */
654	–	i = bary_child(b);
253
656	–	QT_SET_FLAG(*qtptr);
254
658	–	for(;;)
659	–	{
660	–	child = QT_NTH_CHILD_PTR(*qtptr,i);
255
662	–	if(i != 3)
663	–	nbr = qtTrace_ray(child,b,db0,db1,db2,orig,dir,func,arg1,arg2);
664	–	else
665	–	/* If the center cell- must flip direction signs */
666	–	nbr =qtTrace_ray(child,b,-db0,-db1,-db2,orig,dir,func,arg1,arg2);
667	–	if(nbr == QT_DONE)
668	–	return(nbr);
256
257	<	/* If in same block: traverse */
258	<	if(i==3)
259	<	next = nbr;
260	<	else
261	<	if(nbr == i)
675	<	next = 3;
676	<	else
677	<	{
678	<	/* reset the barycentric coordinates in the parents*/
679	<	bary_parent(b,i);
680	<	/* Else pop up to parent and traverse from there */
681	<	return(nbr);
682	<	}
683	<	bary_from_child(b,i,next);
684	<	i = next;
685	<	}
686	<	}
687	<	else
688	<	{
689	<	#ifdef DEBUG_TEST_DRIVER
690	<	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
691	<	Pick_v2[Pick_parent],a1,b1,c1,i,
692	<	Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
693	<	Pick_v2[Pick_cnt]);
694	<	Pick_cnt++;
695	<	#endif
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	<	if(func(qtptr,orig,dir,arg1,arg2) == QT_DONE)
264	<	return(QT_DONE);
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	>	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	<	/* Advance to next node */
283	<	/* NOTE: Optimize: should only have to check 1/2 */
284	<	nbr = move_to_nbr(b,db0,db1,db2,&t);
282	>	/* First check if can trivial accept: if vertex in cell */
283	>	if(s0 & s1 & s2)
284	>	goto Lfunc_call;
285
286	<	if(nbr != -1)
287	<	{
288	<	b[0] += t * db0;
289	<	b[1] += t * db1;
290	<	b[2] += t * db2;
291	<	}
292	<	return(nbr);
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	<
351	<	}
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	<	int
365	<	qtRoot_trace_ray(qtptr,q0,q1,q2,orig,dir,func,arg1,arg2)
366	<	QUADTREE *qtptr;
367	<	FVECT q0,q1,q2;
368	<	FVECT orig,dir;
369	<	int (*func)();
370	<	int *arg1,arg2;
371	<	{
372	<	int i,x,y,nbr;
373	<	QUADTREE *child;
374	<	FVECT n,c,i_pt,d;
375	<	double pd,b[3],db[3],t;
376	<	/* Determine if point lies within pyramid (and therefore
377	<	inside a spherical quadtree cell):GT_INTERIOR, on one of the
378	<	pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
379	<	or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
380	<	For each triangle edge: compare the
732	<	point against the plane formed by the edge and the view center
733	<	*/
734	<	i = point_in_stri(q0,q1,q2,orig);
735	<
736	<	/* Not in this triangle */
737	<	if(!i)
738	<	return(INVALID);
739	<	/* Project the origin onto the root node plane */
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	<	/* Find the intersection point of the origin */
742	<	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
743	<	intersect_vector_plane(orig,n,pd,NULL,i_pt);
744	<	/* project the dir as well */
745	<	VADD(c,orig,dir);
746	<	intersect_vector_plane(c,n,pd,&t,c);
382	>	return(qt);
383
384	<	/* map to 2d by dropping maximum magnitude component of normal */
385	<	i = max_index(n);
386	<	x = (i+1)%3;
387	<	y = (i+2)%3;
752	<	/* Calculate barycentric coordinates of orig */
753	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b);
754	<	/* Calculate barycentric coordinates of dir */
755	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db);
756	<	if(t < 0.0)
757	<	VSUB(db,b,db);
758	<	else
759	<	VSUB(db,db,b);
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
761	–
762	–	#ifdef DEBUG_TEST_DRIVER
763	–	VCOPY(Pick_v0[Pick_cnt],q0);
764	–	VCOPY(Pick_v1[Pick_cnt],q1);
765	–	VCOPY(Pick_v2[Pick_cnt],q2);
766	–	Pick_cnt++;
767	–	#endif
768	–
769	–	/* trace the ray starting with this node */
770	–	nbr = qtTrace_ray(qtptr,b,db[0],db[1],db[2],orig,dir,func,arg1,arg2);
771	–	return(nbr);
772	–
389		}
390
391
776	–	int
777	–	qtRoot_trace_edge(qtptr,q0,q1,q2,orig,dir,max_t,func,arg1,arg2)
778	–	QUADTREE *qtptr;
779	–	FVECT q0,q1,q2;
780	–	FVECT orig,dir;
781	–	double max_t;
782	–	int (*func)();
783	–	int *arg1,arg2;
784	–	{
785	–	int i,x,y,nbr;
786	–	QUADTREE *child;
787	–	FVECT n,c,i_pt,d;
788	–	double pd,b[3],db[3],t;
789	–	/* Determine if point lies within pyramid (and therefore
790	–	inside a spherical quadtree cell):GT_INTERIOR, on one of the
791	–	pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
792	–	or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
793	–	For each triangle edge: compare the
794	–	point against the plane formed by the edge and the view center
795	–	*/
796	–	i = point_in_stri(q0,q1,q2,orig);
797	–
798	–	/* Not in this triangle */
799	–	if(!i)
800	–	return(-1);
801	–	/* Project the origin onto the root node plane */
392
393	<	/* Find the intersection point of the origin */
394	<	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
395	<	intersect_vector_plane(orig,n,pd,NULL,i_pt);
396	<	/* project the dir as well */
397	<	VADD(c,orig,dir);
398	<	intersect_vector_plane(c,n,pd,&t,c);
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	<	/* map to 2d by dropping maximum magnitude component of normal */
423	<	i = max_index(n);
424	<	x = (i+1)%3;
425	<	y = (i+2)%3;
814	<	/* Calculate barycentric coordinates of orig */
815	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b);
816	<	/* Calculate barycentric coordinates of dir */
817	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db);
818	<	if(t < 0.0)
819	<	VSUB(db,b,db);
820	<	else
821	<	VSUB(db,db,b);
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	<	#ifdef DEBUG_TEST_DRIVER
435	<	VCOPY(Pick_v0[Pick_cnt],q0);
436	<	VCOPY(Pick_v1[Pick_cnt],q1);
437	<	VCOPY(Pick_v2[Pick_cnt],q2);
438	<	Pick_cnt++;
439	<	#endif
440	<	/* trace the ray starting with this node */
441	<	nbr = qtTrace_edge(qtptr,b,db,orig,d,max_t,func,arg1,arg2);
442	<	return(nbr);
443	<
444	<	}
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	<	qtVisit_tri_interior(qtptr,q0,q1,q2,t0,t1,t2,n,func,arg1,arg2)
521	<	QUADTREE *qtptr;
839	<	FVECT q0,q1,q2;
840	<	FVECT t0,t1,t2;
841	<	int n;
842	<	int (*func)();
843	<	int *arg1,arg2;
844	<	{
845	<	int i,found,test;
846	<	QUADTREE *child;
847	<	FVECT c0,c1,c2,a,b,c;
848	<	OBJECT os[QT_MAXSET+1],*optr;
849	<	int w;
850	<
851	<	/* If qt Flag set, or qt vertices interior to t0t1t2-descend */
852	<	tree_modified:
520	>	return(qt);
521	>	Lfunc_call:
522
523	<	if(QT_IS_TREE(*qtptr))
524	<	{
525	<	if(QT_IS_FLAG(*qtptr) ||  point_in_stri(t0,t1,t2,q0))
526	<	{
858	<	QT_SET_FLAG(*qtptr);
859	<	qtSubdivide_tri(q0,q1,q2,a,b,c);
860	<	/* descend to children */
861	<	for(i=0;i < 4; i++)
862	<	{
863	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
864	<	qtNth_child_tri(q0,q1,q2,a,b,c,i,c0,c1,c2);
865	<	qtVisit_tri_interior(child,c0,c1,c2,t0,t1,t2,n+1,
866	<	func,arg1,arg2);
867	<	}
868	<	}
869	<	}
870	<	else
871	<	{
872	<	/* NOTE THIS IN TRI TEST Could be replaced by a flag */
873	<	if(!QT_IS_EMPTY(*qtptr))
874	<	{
875	<	if(qtinset(*qtptr,arg2))
876	<	if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2)==QT_MODIFIED)
877	<	goto tree_modified;
878	<	else
879	<	return;
880	<	}
881	<	if(point_in_stri(t0,t1,t2,q0) )
882	<	if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2)==QT_MODIFIED)
883	<	goto tree_modified;
884	<	}
885	<	}
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	<	int
546	<	qtVisit_tri_edges(qtptr,b,db,wptr,sfactor,func,arg1,arg2)
547	<	QUADTREE *qtptr;
548	<	double b[3],db[3][3];
549	<	int *wptr;
550	<	double sfactor;
551	<	int (*func)();
552	<	int *arg1,arg2;
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	<	int i,found;
559	<	QUADTREE *child;
560	<	int nbr,next,w;
904	<	double t;
905	<	#ifdef DEBUG_TEST_DRIVER
906	<	FVECT a1,b1,c1;
907	<	int Pick_parent = Pick_cnt-1;
908	<	qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
909	<	Pick_v2[Pick_parent],a1,b1,c1);
910	<	#endif
558	>	BCOORD a,b,c;
559	>	BCOORD va[3],vb[3],vc[3];
560	>	unsigned int sa,sb,sc;
561
562	<	if(QT_IS_TREE(*qtptr))
563	<	{
564	<	/* Find the appropriate child and reset the coord */
565	<	i = bary_child(b);
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	<	QT_SET_FLAG(*qtptr);
572	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
573
574	<	for(;;)
575	<	{
576	<	w = *wptr;
577	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
578	<
579	<	if(i != 3)
580	<	{
581	<
582	<	db[w][0] *= 2.0;db[w][1] *= 2.0; db[w][2] *= 2.0;
583	<	nbr = qtVisit_tri_edges(child,b,db,wptr,sfactor*2.0,
929	<	func,arg1,arg2);
930	<	w = *wptr;
931	<	db[w][0] *= 0.5;db[w][1] *= 0.5; db[w][2] *= 0.5;
932	<	}
933	<	else
934	<	{
935	<	db[w][0] *=-2.0;db[w][1] *= -2.0; db[w][2] *= -2.0;
936	<	/* If the center cell- must flip direction signs */
937	<	nbr = qtVisit_tri_edges(child,b,db,wptr,sfactor*(-2.0),
938	<	func,arg1,arg2);
939	<	w = *wptr;
940	<	db[w][0] *=-0.5;db[w][1] *= -0.5; db[w][2] *= -0.5;
941	<	}
942	<	if(nbr == QT_DONE)
943	<	return(nbr);
944	<
945	<	/* If in same block: traverse */
946	<	if(i==3)
947	<	next = nbr;
948	<	else
949	<	if(nbr == i)
950	<	next = 3;
951	<	else
952	<	{
953	<	/* reset the barycentric coordinates in the parents*/
954	<	bary_parent(b,i);
955	<	/* Else pop up to parent and traverse from there */
956	<	return(nbr);
957	<	}
958	<	bary_from_child(b,i,next);
959	<	i = next;
960	<	}
574	>	if(sa==7)
575	>	{
576	>	vb[1] = q0[1];
577	>	vc[2] = q0[2];
578	>	vc[1] = b;
579	>	vb[2] = c;
580	>	vb[0] = vc[0] = a;
581	>	QT_NTH_CHILD(qt,0) = qtInsert_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
582	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
583	>	return(qt);
584		}
585	<	else
585	>	if(sb==7)
586		{
587	<	#ifdef DEBUG_TEST_DRIVER
588	<	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
589	<	Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
590	<	Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
591	<	Pick_cnt++;
592	<	#endif
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	<	if(func(qtptr,arg1,arg2) == QT_DONE)
609	<	return(QT_DONE);
610	<
611	<	/* Advance to next node */
612	<	/* NOTE: Optimize: should only have to check 1/2 */
613	<	w = *wptr;
614	<	while(1)
615	<	{
616	<	nbr = move_to_nbr(b,db[w][0],db[w][1],db[w][2],&t);
617	<
618	<	if(t >= 1.0)
619	<	{
983	<	if(w == 2)
984	<	return(QT_DONE);
985	<	b[0] += db[w][0];
986	<	b[1] += db[w][1];
987	<	b[2] += db[w][2];
988	<	w++;
989	<	db[w][0] *= sfactor;
990	<	db[w][1] *= sfactor;
991	<	db[w][2] *= sfactor;
992	<	}
993	<	else
994	<	if(nbr != INVALID)
995	<	{
996	<	b[0] += t * db[w][0];
997	<	b[1] += t * db[w][1];
998	<	b[2] += t * db[w][2];
999	<	db[w][0] *= (1.0 - t);
1000	<	db[w][1] *= (1.0 - t);
1001	<	db[w][2] *= (1.0 - t);
1002	<	*wptr = w;
1003	<	return(nbr);
1004	<	}
1005	<	else
1006	<	return(INVALID);
1007	<	}
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	<
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	<	int
643	<	qtRoot_visit_tri_edges(qtptr,q0,q1,q2,tri,dir,wptr,func,arg1,arg2)
644	<	QUADTREE *qtptr;
645	<	FVECT q0,q1,q2;
646	<	FVECT tri[3],dir[3];
647	<	int *wptr;
648	<	int (*func)();
649	<	int *arg1,arg2;
642	>	QUADTREE
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	>	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	<	int i,x,y,nbr,w;
655	<	QUADTREE *child;
656	<	FVECT n,c,i_pt,d;
1025	<	double pd,b[3][3],db[3][3],t;
1026	<
1027	<	w = *wptr;
654	>	BCOORD a,b,c;
655	>	BCOORD va[3],vb[3],vc[3];
656	>	unsigned int sa,sb,sc;
657
658	<	/* Project the origin onto the root node plane */
659	<
1031	<	/* Find the intersection point of the origin */
1032	<	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
1033	<	/* map to 2d by dropping maximum magnitude component of normal */
1034	<	i = max_index(n);
1035	<	x = (i+1)%3;
1036	<	y = (i+2)%3;
1037	<	/* Calculate barycentric coordinates for current vertex */
1038	<
1039	<	for(i=0;i < 3; i++)
658	>	/* If a tree: trivial reject and recurse on potential children */
659	>	if(QT_IS_TREE(qt))
660		{
661	<	/* If processing 3rd edge-dont need info for t1 */
662	<	if(i==1 && w==2)
663	<	continue;
664	<	/* project the dir as well */
665	<	intersect_vector_plane(tri[i],n,pd,NULL,i_pt);
666	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[i]);
667	<	VADD(c,tri[i],dir[i]);
668	<	intersect_vector_plane(c,n,pd,&t,c);
669	<	/* Calculate barycentric coordinates of dir */
670	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db[i]);
671	<	if(t < 0.0)
672	<	VSUB(db[i],b[i],db[i]);
673	<	else
674	<	VSUB(db[i],db[i],b[i]);
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	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
669	>
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	>	QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
679	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
680	>	return(qt);
681	>	}
682	>	if(sb==0)
683	>	{
684	>	va[0] = q1[0];
685	>	vc[2] = q0[2];
686	>	va[1] = vc[1] = b;
687	>	va[2] = c;
688	>	vc[0] = a;
689	>	QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
690	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
691	>	return(qt);
692	>	}
693	>	if(sc==0)
694	>	{
695	>	va[0] = q1[0];
696	>	vb[1] = q0[1];
697	>	va[1] = b;
698	>	va[2] = vb[2] = c;
699	>	vb[0] = a;
700	>	QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
701	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
702	>	return(qt);
703	>	}
704	>	va[0] = q1[0];
705	>	vb[1] = q0[1];
706	>	vc[2] = q0[2];
707	>	va[1] = vc[1] = b;
708	>	va[2] = vb[2] = c;
709	>	vb[0] = vc[0] = a;
710	>	/* If outside of all edges: only need to Visit 3  */
711	>	if(sa==7 && sb==7 && sc==7)
712	>	{
713	>	QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,
714	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
715	>	return(qt);
716	>	}
717	>	if(sa!=7)
718	>	QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
719	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
720	>	if(sb!=7)
721	>	QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
722	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
723	>	if(sc!=7)
724	>	QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
725	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
726	>
727	>	QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
728	>	t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
729	>	return(qt);
730		}
731	<	#ifdef DEBUG_TEST_DRIVER
732	<	VCOPY(Pick_v0[Pick_cnt],q0);
733	<	VCOPY(Pick_v1[Pick_cnt],q1);
734	<	VCOPY(Pick_v2[Pick_cnt],q2);
1060	<	Pick_cnt++;
1061	<	#endif
1062	<	/* trace the ray starting with this node */
1063	<	nbr = qtVisit_tri_edges(qtptr,b[w],db,wptr,1.0,func,arg1,arg2);
1064	<	return(nbr);
1065	<
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
738
739
740	<	/* NOTE: SINCE DIR could be unit: then we could use integer math */
741	<	int
742	<	qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
743	<	db,wptr,t,sign,sfactor,func,arg1,arg2)
744	<	QUADTREE *qtptr;
745	<	double b[3],db0,db1,db2,db[3][3];
746	<	int *wptr;
747	<	double t[3];
748	<	int sign;
749	<	double sfactor;
750	<	int (*func)();
751	<	int *arg1,arg2;
740	>	/*************************************************************************/
741	>	/* Visit code for applying functions: NOTE THIS IS THE SAME AS INSERT CODE:
742	>	except sets flags: wanted insert to be as efficient as possible: so
743	>	broke into two sets of routines
744	>	*/
745	>
746	>	qtVisit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
747	>	s0,s1,s2,sq0,sq1,sq2,f)
748	>	int root;
749	>	QUADTREE qt;
750	>	BCOORD q0[3],q1[3],q2[3];
751	>	BCOORD t0[3],t1[3],t2[3];
752	>	BCOORD dt10[3],dt21[3],dt02[3];
753	>	BCOORD scale;
754	>	unsigned int s0,s1,s2,sq0,sq1,sq2;
755	>	FUNC f;
756		{
757	<	int i,found;
758	<	QUADTREE *child;
759	<	int nbr,next,w;
1087	<	double t_l,t_g;
1088	<	#ifdef DEBUG_TEST_DRIVER
1089	<	FVECT a1,b1,c1;
1090	<	int Pick_parent = Pick_cnt-1;
1091	<	qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
1092	<	Pick_v2[Pick_parent],a1,b1,c1);
1093	<	#endif
1094	<	if(QT_IS_TREE(*qtptr))
1095	<	{
1096	<	/* Find the appropriate child and reset the coord */
1097	<	i = bary_child(b);
757	>	BCOORD a,b,c;
758	>	BCOORD va[3],vb[3],vc[3];
759	>	unsigned int sa,sb,sc;
760
761	<	QT_SET_FLAG(*qtptr);
761	>	/* If a tree: trivial reject and recurse on potential children */
762	>	if(QT_IS_TREE(qt))
763	>	{
764	>	QT_SET_FLAG(qt);
765
766	<	for(;;)
767	<	{
768	<	w = *wptr;
769	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
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	<	if(i != 3)
1107	<	nbr = qtVisit_tri_edges2(child,b,db0,db1,db2,
1108	<	db,wptr,t,sign,
1109	<	sfactor*2.0,func,arg1,arg2);
1110	<	else
1111	<	/* If the center cell- must flip direction signs */
1112	<	nbr = qtVisit_tri_edges2(child,b,-db0,-db1,-db2,
1113	<	db,wptr,t,1-sign,
1114	<	sfactor*2.0,func,arg1,arg2);
773	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
774
775	<	if(nbr == QT_DONE)
776	<	return(nbr);
777	<	if(*wptr != w)
778	<	{
779	<	w = *wptr;
780	<	db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
781	<	if(sign)
782	<	{  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
783	<	}
784	<	/* If in same block: traverse */
1126	<	if(i==3)
1127	<	next = nbr;
1128	<	else
1129	<	if(nbr == i)
1130	<	next = 3;
1131	<	else
1132	<	{
1133	<	/* reset the barycentric coordinates in the parents*/
1134	<	bary_parent(b,i);
1135	<	/* Else pop up to parent and traverse from there */
1136	<	return(nbr);
1137	<	}
1138	<	bary_from_child(b,i,next);
1139	<	i = next;
1140	<	}
775	>	if(sa==7)
776	>	{
777	>	vb[1] = q0[1];
778	>	vc[2] = q0[2];
779	>	vc[1] = b;
780	>	vb[2] = c;
781	>	vb[0] = vc[0] = a;
782	>	qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
783	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f);
784	>	return;
785		}
786	<	else
786	>	if(sb==7)
787		{
788	<	#ifdef DEBUG_TEST_DRIVER
789	<	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
790	<	Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
791	<	Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
792	<	Pick_cnt++;
793	<	#endif
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	<	if(func(qtptr,arg1,arg2) == QT_DONE)
810	<	return(QT_DONE);
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	<	/* Advance to next node */
824	<	w = *wptr;
825	<	while(1)
826	<	{
827	<	nbr = move_to_nbr(b,db0,db1,db2,&t_l);
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	>	qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
838	>	scale,s0,s1,s2,sq0,sq1,sq2,f);
839
1160	–	t_g = t_l/sfactor;
1161	–	#ifdef DEBUG
1162	–	if(t[w] <= 0.0)
1163	–	eputs("qtVisit_tri_edges2():negative t\n");
1164	–	#endif
1165	–	if(t_g >= t[w])
1166	–	{
1167	–	if(w == 2)
1168	–	return(QT_DONE);
1169	–
1170	–	b[0] += (t[w])*sfactor*db0;
1171	–	b[1] += (t[w])*sfactor*db1;
1172	–	b[2] += (t[w])*sfactor*db2;
1173	–	w++;
1174	–	db0 = db[w][0];
1175	–	db1 = db[w][1];
1176	–	db2 = db[w][2];
1177	–	if(sign)
1178	–	{  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
1179	–	}
1180	–	else
1181	–	if(nbr != INVALID)
1182	–	{
1183	–	b[0] += t_l * db0;
1184	–	b[1] += t_l * db1;
1185	–	b[2] += t_l * db2;
1186	–
1187	–	t[w] -= t_g;
1188	–	*wptr = w;
1189	–	return(nbr);
1190	–	}
1191	–	else
1192	–	return(INVALID);
1193	–	}
1194	–	}
1195	–
840		}
841
842
843	<	int
844	<	qtRoot_visit_tri_edges2(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2)
845	<	QUADTREE *qtptr;
846	<	FVECT q0,q1,q2;
847	<	FVECT tri[3],i_pt;
848	<	int *wptr;
849	<	int (*func)();
850	<	int *arg1,arg2;
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	<	int x,y,z,nbr,w,i,j;
855	<	QUADTREE *child;
856	<	FVECT n,c,d,v[3];
1211	<	double pd,b[4][3],db[3][3],et[3],t[3],exit_pt;
1212	<
1213	<	w = *wptr;
854	>	BCOORD a,b,c;
855	>	BCOORD va[3],vb[3],vc[3];
856	>	unsigned int sa,sb,sc;
857
858	<	/* Project the origin onto the root node plane */
859	<
1217	<	/* Find the intersection point of the origin */
1218	<	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
1219	<	/* map to 2d by dropping maximum magnitude component of normal */
1220	<	z = max_index(n);
1221	<	x = (z+1)%3;
1222	<	y = (z+2)%3;
1223	<	/* Calculate barycentric coordinates for current vertex */
1224	<	if(w != -1)
858	>	/* If a tree: trivial reject and recurse on potential children */
859	>	if(QT_IS_TREE(qt))
860		{
861	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);
862	<	intersect_vector_plane(tri[w],n,pd,&(et[w]),v[w]);
863	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[w][x],v[w][y],b[w]);
864	<	}
865	<	else
866	<	/* Just starting: b[0] is the origin point: guaranteed to be valid b*/
867	<	{
1233	<	w = 0;
1234	<	intersect_vector_plane(tri[0],n,pd,&(et[0]),v[0]);
1235	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);
1236	<	VCOPY(b[3],b[0]);
1237	<	}
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	+	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
870
871	<	j = (w+1)%3;
1241	<	intersect_vector_plane(tri[j],n,pd,&(et[j]),v[j]);
1242	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);
1243	<	if(et[j] < 0.0)
1244	<	{
1245	<	VSUB(db[w],b[3],b[j]);
1246	<	t[w] = FHUGE;
1247	<	}
1248	<	else
1249	<	{
1250	<	VSUB(db[w],b[j],b[3]);
1251	<	t[w] = 1.0;
1252	<	move_to_nbr(b[3],db[w][0],db[w][1],db[w][2],&exit_pt);
1253	<	if(exit_pt >= 1.0)
1254	<	{
1255	<	for(;j < 3;j++)
1256	<	{
1257	<	i = (j+1)%3;
1258	<	if(i!= w)
1259	<	{
1260	<	intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
1261	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
1262	<	v[i][y],b[i]);
1263	<	}
1264	<	if(et[i] < 0.0)
1265	<	{
1266	<	VSUB(db[j],b[j],b[i]);
1267	<	t[j] = FHUGE;
1268	<	break;
1269	<	}
1270	<	else
1271	<	{
1272	<	VSUB(db[j],b[i],b[j]);
1273	<	t[j] = 1.0;
1274	<	}
1275	<	move_to_nbr(b[j],db[j][0],db[j][1],db[j][2],&exit_pt);
1276	<	if(exit_pt < 1.0)
1277	<	break;
1278	<	}
1279	<	}
1280	<	}
1281	<	*wptr = w;
1282	<	/* trace the ray starting with this node */
1283	<	nbr = qtVisit_tri_edges2(qtptr,b[3],db[w][0],db[w][1],db[w][2],
1284	<	db,wptr,t,0,1.0,func,arg1,arg2);
1285	<	if(nbr != INVALID && nbr != QT_DONE)
871	>	if(sa==0)
872		{
873	<	i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
874	<	i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
875	<	i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
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	<	return(nbr);
883	<
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	>	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
938
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	+	/* First check if can trivial accept: if vertex in cell */
956	+	if(s0 & s1 & s2)
957	+	goto Lfunc_call;
958
959	+	/* Assumption: Could not trivial reject*/
960	+	/* IF any coords lie on edge- must be in if couldnt reject */
961	+	a = q1[0];b= q0[1]; c= q0[2];
962	+	if(t0[0]== a || t0[1] == b || t0[2] == c)
963	+	goto Lfunc_call;
964	+	if(t1[0]== a || t1[1] == b || t1[2] == c)
965	+	goto Lfunc_call;
966	+	if(t2[0]== a || t2[1] == b || t2[2] == c)
967	+	goto Lfunc_call;
968	+
969	+	/* Test for edge crossings */
970	+	/* Test t0t1,t1t2,t2t0 against a */
971	+	/* Calculate edge crossings */
972	+	e0  = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
973	+	/* See if edge can be trivially rejected from intersetion testing */
974	+	test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
975	+	((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
976	+	bl=bh=0;
977	+	if(test_t0t1 && (e0 & 2) )
978	+	{
979	+	/* Must do double calculation to avoid overflow */
980	+	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
981	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
982	+	}
983	+	test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
984	+	((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
985	+	if(test_t1t2 && (e0 & 1))
986	+	{    /* test t1t2 against a */
987	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
988	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
989	+	}
990	+	test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
991	+	((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
992	+	if(test_t2t0 && (e0 & 4))
993	+	{
994	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
995	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
996	+	}
997	+	e1  = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
998	+	cl = ch = 0;
999	+	if(test_t0t1 && (e1 & 2))
1000	+	{/* test t0t1 against b */
1001	+	db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
1002	+	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1003	+	}
1004	+	if(test_t1t2 && (e1 & 1))
1005	+	{/* test t1t2 against b */
1006	+	db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
1007	+	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1008	+	}
1009	+	if(test_t2t0 && (e1 & 4))
1010	+	{/* test t2t0 against b */
1011	+	db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
1012	+	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1013	+	}
1014	+
1015	+	/* test edges against c */
1016	+	e2  = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
1017	+	al = ah = 0;
1018	+	if(test_t0t1 && (e2 & 2))
1019	+	{ /* test t0t1 against c */
1020	+	db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
1021	+	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1022	+	}
1023	+	if(test_t1t2 && (e2 & 1))
1024	+	{
1025	+	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1026	+	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1027	+	}
1028	+	if(test_t2t0 && (e2 & 4))
1029	+	{ /* test t2t0 against c */
1030	+	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1031	+	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1032	+	}
1033	+	/* Only remaining case is if some edge trivially rejected */
1034	+	if(!e0 || !e1 || !e2)
1035	+	return;
1036
1037	+	/* Only one/none got tested - pick either of other two/any two */
1038	+	/* Only need to test one edge */
1039	+	if(!test_t0t1 && (e0 & 2))
1040	+	{
1041	+	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1042	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1043	+	}
1044	+	if(!test_t1t2 && (e0 & 1))
1045	+	{/* test t1t2 against a */
1046	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1047	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1048	+	}
1049	+	if(!test_t2t0 && (e0 & 4))
1050	+	{
1051	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1052	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1053	+	}
1054
1055	+	return;
1056
1057	+	Lfunc_call:
1058	+	f.func(f.argptr,root,qt);
1059	+	if(!QT_IS_EMPTY(qt))
1060	+	QT_LEAF_SET_FLAG(qt);
1061	+	}
1062
1063
1064
1065	+	/* Leaf node: Do definitive test */
1066	+	QUADTREE
1067	+	qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1068	+	scale, s0,s1,s2,sq0,sq1,sq2,f)
1069	+	int root;
1070	+	QUADTREE qt;
1071	+	BCOORD q0[3],q1[3],q2[3];
1072	+	BCOORD t0[3],t1[3],t2[3];
1073	+	BCOORD dt10[3],dt21[3],dt02[3];
1074	+	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
1075	+	FUNC f;
1076	+	{
1077	+	double db;
1078	+	unsigned int e0,e1,e2;
1079	+	BCOORD a,b,c;
1080	+	double p0[2], p1[2],cp;
1081	+	char al,ah,bl,bh,cl,ch;
1082	+	char testl,testh,test_t0t1,test_t1t2,test_t2t0;
1083	+	unsigned int ls0,ls1,ls2;
1084	+
1085	+	/* May have gotten passed trivial reject if one/two vertices are ON */
1086	+	a = q1[0];b= q0[1]; c= q0[2];
1087	+	SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
1088	+
1089	+	/* First check if can trivial accept: if vertex in cell */
1090	+	if(ls0 & ls1 & ls2)
1091	+	goto Lfunc_call;
1092
1093	+	if(ls0==0 || ls1 == 0 || ls2 ==0)
1094	+	return;
1095	+	/* Assumption: Could not trivial reject*/
1096	+	/* IF any coords lie on edge- must be in if couldnt reject */
1097
1098	+	if(t0[0]== a || t0[1] == b || t0[2] == c)
1099	+	goto Lfunc_call;
1100	+	if(t1[0]== a || t1[1] == b || t1[2] == c)
1101	+	goto Lfunc_call;
1102	+	if(t2[0]== a || t2[1] == b || t2[2] == c)
1103	+	goto Lfunc_call;
1104
1105	+	/* Test for edge crossings */
1106	+	/* Test t0t1 against a,b,c */
1107	+	/* First test if t0t1 can be trivially rejected */
1108	+	/* If both edges are outside bounds- dont need to test */
1109	+
1110	+	/* Test t0t1,t1t2,t2t0 against a */
1111	+	test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
1112	+	((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
1113	+	e0  = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
1114	+	bl=bh=0;
1115	+	/* Test t0t1,t1t2,t2t0 against a */
1116	+	if(test_t0t1 && (e0 & 2) )
1117	+	{
1118	+	db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
1119	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1120	+	}
1121	+	test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
1122	+	((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
1123	+	if(test_t1t2 && (e0 & 1))
1124	+	{    /* test t1t2 against a */
1125	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1126	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1127	+	}
1128	+	test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
1129	+	((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
1130	+	if(test_t2t0 && (e0 & 4))
1131	+	{
1132	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1133	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1134	+	}
1135	+	cl = ch = 0;
1136	+	e1  = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
1137	+	if(test_t0t1 && (e1 & 2))
1138	+	{/* test t0t1 against b */
1139	+	db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
1140	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1141	+	}
1142	+	if(test_t1t2 && (e1 & 1))
1143	+	{/* test t1t2 against b */
1144	+	db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
1145	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1146	+	}
1147	+	if(test_t2t0 && (e1 & 4))
1148	+	{/* test t2t0 against b */
1149	+	db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
1150	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1151	+	}
1152	+	al = ah = 0;
1153	+	e2  = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
1154	+	if(test_t0t1 && (e2 & 2))
1155	+	{ /* test t0t1 against c */
1156	+	db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
1157	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1158	+	}
1159	+	if(test_t1t2 && (e2 & 1))
1160	+	{
1161	+	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1162	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1163	+	}
1164	+	if(test_t2t0 && (e2 & 4))
1165	+	{ /* test t2t0 against c */
1166	+	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1167	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1168	+	}
1169	+	/* Only remaining case is if some edge trivially rejected */
1170	+	if(!e0 || !e1 || !e2)
1171	+	return;
1172
1173	+	/* Only one/none got tested - pick either of other two/any two */
1174	+	/* Only need to test one edge */
1175	+	if(!test_t0t1 && (e0 & 2))
1176	+	{
1177	+	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1178	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1179	+	}
1180	+	if(!test_t1t2 && (e0 & 1))
1181	+	{/* test t1t2 against a */
1182	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1183	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1184	+	}
1185	+	if(!test_t2t0 && (e0 & 4))
1186	+	{
1187	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1188	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1189	+	}
1190	+
1191	+	return;
1192	+
1193	+	Lfunc_call:
1194	+	f.func(f.argptr,root,qt);
1195	+	if(!QT_IS_EMPTY(qt))
1196	+	QT_LEAF_SET_FLAG(qt);
1197	+	}
1198
1199
1200

Diff Legend

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