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

Diff Legend

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