ViewVC Help

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

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

Comparing ray/src/hd/sm_qtree.c (file contents):
Revision 3.3 by gwlarson, Tue Aug 25 11:03:27 1998 UTC vs.
Revision 3.6 by gwlarson, Wed Sep 16 18:16:29 1998 UTC

#	Line 17 | Line 17 | static char SCCSid[] = "$SunId$ SGI";
17
18		#include "sm_geom.h"
19		#include "sm_qtree.h"
20	–	#include "object.h"
20
21		QUADTREE  *quad_block[QT_MAX_BLK];              /* our quadtree */
22		static QUADTREE  quad_free_list = EMPTY;  /* freed octree nodes */
23		static QUADTREE  treetop = 0;           /* next free node */
24	<	int4 *quad_flag;
24	>	int4 *quad_flag= NULL;
25
26	+	#ifdef TEST_DRIVER
27	+	extern FVECT Pick_v0[500],Pick_v1[500],Pick_v2[500];
28	+	extern int Pick_cnt,Pick_tri,Pick_samp;
29	+	extern FVECT Pick_point[500];
30	+	#endif
31	+
32		QUADTREE
33		qtAlloc()                       /* allocate a quadtree */
34		{
#	Line 42 | Line 47 | qtAlloc()                      /* allocate a quadtree */
47		if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
48		QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
49		return(EMPTY);
50	+
51		quad_flag = (int4 *)realloc((char *)quad_flag,
52	<	(QT_BLOCK(freet)+1)*QT_BLOCK_SIZE/(8*sizeof(int4)));
52	>	(QT_BLOCK(freet)+1)*(QT_BLOCK_SIZE/8));
53		if (quad_flag == NULL)
54		return(EMPTY);
55		}
#	Line 55 | Line 61 | qtAlloc()                      /* allocate a quadtree */
61		qtClearAllFlags()               /* clear all quadtree branch flags */
62		{
63		if (!treetop) return;
64	<	bzero((char *)quad_flag,
59	<	(QT_BLOCK(treetop-1)+1)*QT_BLOCK_SIZE/(8*sizeof(int4)));
64	>	bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*(QT_BLOCK_SIZE/8));
65		}
66
67
#	Line 81 | Line 86 | qtDone()                       /* free EVERYTHING */
86		{
87		register int    i;
88
89	+	qtfreeleaves();
90		for (i = 0; i < QT_MAX_BLK; i++)
91		{
92		if (quad_block[i] == NULL)
#	Line 94 | Line 100 | qtDone()                       /* free EVERYTHING */
100		treetop = 0;
101		}
102
97	–
103		QUADTREE
99	–	qtCompress(qt)                  /* recursively combine nodes */
100	–	register QUADTREE  qt;
101	–	{
102	–	register int  i;
103	–	register QUADTREE  qres;
104	–
105	–	if (!QT_IS_TREE(qt))        /* not a tree */
106	–	return(qt);
107	–	qres = QT_NTH_CHILD(qt,0) = qtCompress(QT_NTH_CHILD(qt,0));
108	–	for (i = 1; i < 4; i++)
109	–	if((QT_NTH_CHILD(qt,i) = qtCompress(QT_NTH_CHILD(qt,i))) != qres)
110	–	qres = qt;
111	–	if(!QT_IS_TREE(qres))
112	–	{   /* all were identical leaves */
113	–	QT_NTH_CHILD(qt,0) = quad_free_list;
114	–	quad_free_list = qt;
115	–	}
116	–	return(qres);
117	–	}
118	–
119	–
120	–	QUADTREE
104		*qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
105		QUADTREE *qtptr;
106		double bcoord[3];
#	Line 129 | Line 112 | FVECT t0,t1,t2;
112
113		if(QT_IS_TREE(*qtptr))
114		{
115	<	i = bary2d_child(bcoord);
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		{
#	Line 143 | Line 136 | FVECT t0,t1,t2;
136		}
137
138
146	–
147	–	int
148	–	qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n)
149	–	QUADTREE *qtptr;
150	–	double bcoord[3];
151	–	int id;
152	–	FVECT v0,v1,v2;
153	–	int n;
154	–	{
155	–	int i;
156	–	QUADTREE *child;
157	–	OBJECT os[MAXSET+1],*optr;
158	–	int found;
159	–	FVECT r0,r1,r2;
160	–
161	–	if(QT_IS_TREE(*qtptr))
162	–	{
163	–	QT_SET_FLAG(*qtptr);
164	–	i = bary2d_child(bcoord);
165	–	child = QT_NTH_CHILD_PTR(*qtptr,i);
166	–	return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,++n));
167	–	}
168	–	else
169	–	{
170	–	/* If this leave node emptry- create a new set */
171	–	if(QT_IS_EMPTY(*qtptr))
172	–	*qtptr = qtaddelem(*qtptr,id);
173	–	else
174	–	{
175	–	qtgetset(os,*qtptr);
176	–	/* If the set is too large: subdivide */
177	–	if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
178	–	*qtptr = qtaddelem(*qtptr,id);
179	–	else
180	–	{
181	–	if (n < QT_MAX_LEVELS)
182	–	{
183	–	/* If set size exceeds threshold: subdivide cell and
184	–	reinsert set tris into cell
185	–	*/
186	–	n++;
187	–	qtfreeleaf(*qtptr);
188	–	qtSubdivide(qtptr);
189	–	QT_SET_FLAG(*qtptr);
190	–	found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,n);
191	–	#if 0
192	–	if(!found)
193	–	{
194	–	#ifdef TEST_DRIVER
195	–	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
196	–	#else
197	–	eputs("qtAdd_tri():Found in parent but not children\n");
198	–	#endif
199	–	}
200	–	#endif
201	–	for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
202	–	{
203	–	id = QT_SET_NEXT_ELEM(optr);
204	–	qtTri_verts_from_id(id,r0,r1,r2);
205	–	found= qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,++n);
206	–	#ifdef DEBUG
207	–	if(!found)
208	–	eputs("qtAdd_tri():Reinsert-in parent but not children\n");
209	–	#endif
210	–	}
211	–	}
212	–	else
213	–	if(QT_SET_CNT(os) < QT_MAX_SET)
214	–	{
215	–	*qtptr = qtaddelem(*qtptr,id);
216	–	}
217	–	else
218	–	{
219	–	#ifdef DEBUG
220	–	eputs("qtAdd_tri():two many levels\n");
221	–	#endif
222	–	return(FALSE);
223	–	}
224	–	}
225	–	}
226	–	}
227	–	return(TRUE);
228	–	}
229	–
230	–
231	–	int
232	–	qtAdd_tri_from_point(qtptr,v0,v1,v2,pt,id)
233	–	QUADTREE *qtptr;
234	–	FVECT v0,v1,v2;
235	–	FVECT pt;
236	–	int id;
237	–	{
238	–	char d,w;
239	–	int i,x,y;
240	–	QUADTREE *child;
241	–	QUADTREE qt;
242	–	FVECT i_pt,n,a,b,c,r0,r1,r2;
243	–	double pd,bcoord[3];
244	–	OBJECT os[MAXSET+1],*optr;
245	–	int found;
246	–
247	–	/* Determine if point lies within pyramid (and therefore
248	–	inside a spherical quadtree cell):GT_INTERIOR, on one of the
249	–	pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
250	–	or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
251	–	For each triangle edge: compare the
252	–	point against the plane formed by the edge and the view center
253	–	*/
254	–	d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);
255	–
256	–	/* Not in this triangle */
257	–	if(!d)
258	–	return(FALSE);
259	–
260	–	/* Will return lowest level triangle containing point: It the
261	–	point is on an edge or vertex: will return first associated
262	–	triangle encountered in the child traversal- the others can
263	–	be derived using triangle adjacency information
264	–	*/
265	–	if(QT_IS_TREE(*qtptr))
266	–	{
267	–	QT_SET_FLAG(*qtptr);
268	–	/* Find the intersection point */
269	–	tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
270	–	intersect_vector_plane(pt,n,pd,NULL,i_pt);
271	–
272	–	i = max_index(n);
273	–	x = (i+1)%3;
274	–	y = (i+2)%3;
275	–	/* Calculate barycentric coordinates of i_pt */
276	–	bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
277	–
278	–	i = bary2d_child(bcoord);
279	–	child = QT_NTH_CHILD_PTR(*qtptr,i);
280	–	/* NOTE: Optimize: only subdivide for specified child */
281	–	qtSubdivide_tri(v0,v1,v2,a,b,c);
282	–	qtNth_child_tri(v0,v1,v2,a,b,c,i,v0,v1,v2);
283	–	return(qtLeaf_add_tri_from_pt(child,bcoord,id,v0,v1,v2,1));
284	–	}
285	–	else
286	–	{
287	–	/* If this leave node emptry- create a new set */
288	–	if(QT_IS_EMPTY(*qtptr))
289	–	*qtptr = qtaddelem(*qtptr,id);
290	–	else
291	–	{
292	–	qtgetset(os,*qtptr);
293	–	/* If the set is too large: subdivide */
294	–	if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
295	–	*qtptr = qtaddelem(*qtptr,id);
296	–	else
297	–	{
298	–	/* If set size exceeds threshold: subdivide cell and
299	–	reinsert set tris into cell
300	–	*/
301	–	qtfreeleaf(*qtptr);
302	–	qtSubdivide(qtptr);
303	–	found = qtLeaf_add_tri_from_pt(qtptr,bcoord,id,v0,v1,v2,1);
304	–	#if 0
305	–	if(!found)
306	–	{
307	–	#ifdef TEST_DRIVER
308	–	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
309	–	#else
310	–	eputs("qtAdd_tri():Found in parent but not children\n");
311	–	#endif
312	–	}
313	–	#endif
314	–	for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
315	–	{
316	–	id = QT_SET_NEXT_ELEM(optr);
317	–	qtTri_verts_from_id(id,r0,r1,r2);
318	–	found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,1);
319	–	#ifdef DEBUG
320	–	if(!found)
321	–	eputs("qtAdd_tri():Reinsert-in parent but not children\n");
322	–	#endif
323	–	}
324	–	}
325	–	}
326	–	}
327	–	return(TRUE);
328	–	}
329	–
330	–
139		QUADTREE
140		*qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
141		QUADTREE *qtptr;
#	Line 335 | Line 143 | FVECT v0,v1,v2;
143		FVECT pt;
144		FVECT t0,t1,t2;
145		{
146	<	char d,w;
146	>	int d;
147		int i,x,y;
148		QUADTREE *child;
341	–	QUADTREE qt;
149		FVECT n,i_pt,a,b,c;
150		double pd,bcoord[3];
151
#	Line 349 | Line 156 | FVECT t0,t1,t2;
156		For each triangle edge: compare the
157		point against the plane formed by the edge and the view center
158		*/
159	<	d = test_single_point_against_spherical_tri(v0,v1,v2,pt,&w);
159	>	d = point_in_stri(v0,v1,v2,pt);
160
161	+
162		/* Not in this triangle */
163		if(!d)
356	–	{
164		return(NULL);
358	–	}
165
166		/* Will return lowest level triangle containing point: It the
167		point is on an edge or vertex: will return first associated
#	Line 368 | Line 174 | FVECT t0,t1,t2;
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);
177	>	i = max_index(n,NULL);
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 = bary2d_child(bcoord);
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);
#	Line 384 | Line 204 | FVECT t0,t1,t2;
204		return(qtLocate_leaf(child,bcoord,t0,t1,t2));
205		}
206		else
387	–	return(qtptr);
388	–	}
389	–
390	–	int
391	–	qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
392	–	QUADTREE *qtptr;
393	–	FVECT v0,v1,v2;
394	–	FVECT pt;
395	–	char *type,*which;
396	–	{
397	–	OBJECT os[MAXSET+1],*optr;
398	–	char d,w;
399	–	int i,id;
400	–	FVECT p0,p1,p2;
401	–
402	–	qtptr = qtRoot_point_locate(qtptr,v0,v1,v2,pt,NULL,NULL,NULL);
403	–
404	–	if(!qtptr || QT_IS_EMPTY(*qtptr))
405	–	return(EMPTY);
406	–
407	–	/* Get the set */
408	–	qtgetset(os,*qtptr);
409	–	for (i = QT_SET_CNT(os),optr = QT_SET_PTR(os); i > 0; i--)
207		{
208	<	/* Find the triangle that pt falls in (NOTE:FOR now return first 1) */
209	<	id = QT_SET_NEXT_ELEM(optr);
210	<	qtTri_verts_from_id(id,p0,p1,p2);
211	<	d = test_single_point_against_spherical_tri(p0,p1,p2,pt,&w);
212	<	if(d)
416	<	{
417	<	if(type)
418	<	*type = d;
419	<	if(which)
420	<	*which = w;
421	<	return(id);
208	>	if(t0)
209	>	{
210	>	VCOPY(t0,v0);
211	>	VCOPY(t1,v1);
212	>	VCOPY(t2,v2);
213		}
214	+	return(qtptr);
215		}
424	–	return(EMPTY);
216		}
217
218	+
219		QUADTREE
220		qtSubdivide(qtptr)
221		QUADTREE *qtptr;
#	Line 475 | Line 267 | 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);
#	Line 500 | Line 293 | into the new child cells: it is assumed that "v1,v2,v3
293		*/
294
295		int
296	<	qtRoot_add_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n)
296	>	qtRoot_add_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
297		QUADTREE *qtptr;
298	<	int id;
298	>	FVECT q0,q1,q2;
299		FVECT t0,t1,t2;
507	–	FVECT v0,v1,v2;
508	–	int n;
509	–	{
510	–	char test;
511	–	int found;
512	–
513	–	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
514	–	if(!test)
515	–	return(FALSE);
516	–
517	–	found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
518	–	return(found);
519	–	}
520	–
521	–	int
522	–	qtRoot_add_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n)
523	–	QUADTREE *qtptr;
300		int id;
525	–	FVECT t0,t1,t2;
526	–	FVECT v0,v1,v2;
301		int n;
302		{
303	<	char test;
303	>	int test;
304		int found;
305
306	<	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
306	>	test = stri_intersect(q0,q1,q2,t0,t1,t2);
307		if(!test)
308		return(FALSE);
309
310	<	found = qtAdd_tri_from_point(qtptr,id,t0,t1,t2,v0,v1,v2,n);
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,id,t0,t1,t2,v0,v1,v2,n)
316	>	qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
317		QUADTREE *qtptr;
318	<	int id;
318	>	FVECT q0,q1,q2;
319		FVECT t0,t1,t2;
320	<	FVECT v0,v1,v2;
320	>	int id;
321		int n;
322		{
323	<
549	<	char test;
550	<	int i,index;
323	>	int i,index,test,found;
324		FVECT a,b,c;
325	<	OBJECT os[MAXSET+1],*optr;
325	>	OBJECT os[QT_MAXSET+1],*optr;
326		QUADTREE qt;
554	–	int found;
327		FVECT r0,r1,r2;
328
329		found = 0;
#	Line 559 | Line 331 | int n;
331		if(QT_IS_TREE(*qtptr))
332		{
333		n++;
334	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
335	<	test = spherical_tri_intersect(t0,t1,t2,v0,a,c);
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),id,t0,t1,t2,v0,a,c,n);
338	<	test = spherical_tri_intersect(t0,t1,t2,a,v1,b);
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),id,t0,t1,t2,a,v1,b,n);
341	<	test = spherical_tri_intersect(t0,t1,t2,c,b,v2);
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),id,t0,t1,t2,c,b,v2,n);
344	<	test = spherical_tri_intersect(t0,t1,t2,b,c,a);
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),id,t0,t1,t2,b,c,a,n);
575	<
576	<	#if 0
577	<	if(!found)
578	<	{
579	<	#ifdef TEST_DRIVER
580	<	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
581	<	#else
582	<	eputs("qtAdd_tri():Found in parent but not children\n");
583	<	#endif
584	<	}
585	<	#endif
346	>	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),b,c,a,t0,t1,t2,id,n);
347		}
348		else
349		{
#	Line 591 | Line 352 | int n;
352		*qtptr = qtaddelem(*qtptr,id);
353		else
354		{
594	–	qtgetset(os,*qtptr);
355		/* If the set is too large: subdivide */
356	<	if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
357	<	*qtptr = qtaddelem(*qtptr,id);
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)
#	Line 602 | Line 364 | int n;
364		/* If set size exceeds threshold: subdivide cell and
365		reinsert set tris into cell
366		*/
367	<	n++;
368	<	qtfreeleaf(*qtptr);
369	<	qtSubdivide(qtptr);
370	<	found = qtAdd_tri(qtptr,id,t0,t1,t2,v0,v1,v2,n);
371	<	#if 0
372	<	if(!found)
373	<	{
374	<	#ifdef TEST_DRIVER
375	<	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
376	<	#else
377	<	eputs("qtAdd_tri():Found in parent but not children\n");
378	<	#endif
617	<	}
618	<	#endif
619	<	for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
620	<	{
621	<	id = QT_SET_NEXT_ELEM(optr);
622	<	qtTri_verts_from_id(id,r0,r1,r2);
623	<	found=qtAdd_tri(qtptr,id,r0,r1,r2,v0,v1,v2,n);
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,r0,r1,r2,NULL,NULL,NULL,NULL,NULL,NULL);
378	>	found=qtAdd_tri(qtptr,q0,q1,q2,r0,r1,r2,id,n);
379		#ifdef DEBUG
380	<	if(!found)
380	>	if(!found)
381		eputs("qtAdd_tri():Reinsert-in parent but not children\n");
382		#endif
383	<	}
384	<	}
383	>	}
384	>	}
385		else
386	<	if(QT_SET_CNT(os) < QT_MAX_SET)
632	<	{
386	>	if(QT_SET_CNT(optr) < QT_MAXSET)
387		*qtptr = qtaddelem(*qtptr,id);
634	–	#if 0
635	–	{
636	–	int k;
637	–	qtgetset(os,*qtptr);
638	–	printf("\n%d:\n",os[0]);
639	–	for(k=1; k <= os[0];k++)
640	–	printf("%d ",os[k]);
641	–	printf("\n");
642	–	}
643	–	#endif
644	–	/*
645	–	insertelem(os,id);
646	–	*qtptr = fullnode(os);
647	–	*/
648	–	}
388		else
389		{
390		#ifdef DEBUG
#	Line 666 | Line 405 | QUADTREE *qtptr;
405		FVECT t0,t1,t2;
406		FVECT v0,v1,v2;
407		int (*func)();
408	<	char *arg;
408	>	int *arg;
409		{
410	<	char test;
410	>	int test;
411		FVECT a,b,c;
412
413		/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
414	<	test = spherical_tri_intersect(t0,t1,t2,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	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
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		}
690	–	else
691	–	return(func(qtptr,arg));
432		}
433
694	–
434		int
435		qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
436		QUADTREE *qtptr;
#	Line 700 | Line 439 | FVECT t0,t1,t2;
439		FVECT v0,v1,v2;
440		{
441
442	<	char test;
442	>	int test;
443		int i;
444		FVECT a,b,c;
445	<	OBJECT os[MAXSET+1];
445	>	OBJECT os[QT_MAXSET+1];
446
447		/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
448	<	test = spherical_tri_intersect(t0,t1,t2,v0,v1,v2);
448	>	test = stri_intersect(t0,t1,t2,v0,v1,v2);
449
450		/* If triangles do not overlap: done */
451		if(!test)
#	Line 732 | Line 471 | FVECT v0,v1,v2;
471		/* remove id from set */
472		else
473		{
474	<	qtgetset(os,*qtptr);
736	<	if(!inset(os,id))
474	>	if(!qtinset(*qtptr,id))
475		{
476		#ifdef DEBUG
477		eputs("qtRemove_tri(): tri not in set\n");
#	Line 747 | Line 485 | FVECT v0,v1,v2;
485		}
486		return(TRUE);
487		}
488	+
489	+
490	+	int
491	+	move_to_nbr(b,db0,db1,db2,tptr)
492	+	double b[3],db0,db1,db2;
493	+	double *tptr;
494	+	{
495	+	double t,dt;
496	+	int nbr;
497	+
498	+	nbr = -1;
499	+	/* Advance to next node */
500	+	if(!ZERO(db0) && db0 < 0.0)
501	+	{
502	+	t = -b[0]/db0;
503	+	nbr = 0;
504	+	}
505	+	else
506	+	t = FHUGE;
507	+	if(!ZERO(db1) && db1 < 0.0 )
508	+	{
509	+	dt = -b[1]/db1;
510	+	if( dt < t)
511	+	{
512	+	t = dt;
513	+	nbr = 1;
514	+	}
515	+	}
516	+	if(!ZERO(db2) && db2 < 0.0 )
517	+	{
518	+	dt = -b[2]/db2;
519	+	if( dt < t)
520	+	{
521	+	t = dt;
522	+	nbr = 2;
523	+	}
524	+	}
525	+	*tptr = t;
526	+	return(nbr);
527	+	}
528	+
529	+	int
530	+	qtTrace_ray(qtptr,b,db0,db1,db2,orig,dir,func,arg1,arg2)
531	+	QUADTREE *qtptr;
532	+	double b[3],db0,db1,db2;
533	+	FVECT orig,dir;
534	+	int (*func)();
535	+	int *arg1,arg2;
536	+	{
537	+
538	+	int i,found;
539	+	QUADTREE *child;
540	+	int nbr,next;
541	+	double t;
542	+	#ifdef DEBUG_TEST_DRIVER
543	+
544	+	FVECT a1,b1,c1;
545	+	int Pick_parent = Pick_cnt-1;
546	+	qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
547	+	Pick_v2[Pick_parent],a1,b1,c1);
548	+
549	+	#endif
550	+	if(QT_IS_TREE(*qtptr))
551	+	{
552	+	/* Find the appropriate child and reset the coord */
553	+	i = bary_child(b);
554	+
555	+	QT_SET_FLAG(*qtptr);
556	+
557	+	for(;;)
558	+	{
559	+	child = QT_NTH_CHILD_PTR(*qtptr,i);
560	+
561	+	if(i != 3)
562	+	nbr = qtTrace_ray(child,b,db0,db1,db2,orig,dir,func,arg1,arg2);
563	+	else
564	+	/* If the center cell- must flip direction signs */
565	+	nbr =qtTrace_ray(child,b,-db0,-db1,-db2,orig,dir,func,arg1,arg2);
566	+	if(nbr == QT_DONE)
567	+	return(nbr);
568	+
569	+	/* If in same block: traverse */
570	+	if(i==3)
571	+	next = nbr;
572	+	else
573	+	if(nbr == i)
574	+	next = 3;
575	+	else
576	+	{
577	+	/* reset the barycentric coordinates in the parents*/
578	+	bary_parent(b,i);
579	+	/* Else pop up to parent and traverse from there */
580	+	return(nbr);
581	+	}
582	+	bary_from_child(b,i,next);
583	+	i = next;
584	+	}
585	+	}
586	+	else
587	+	{
588	+	#ifdef DEBUG_TEST_DRIVER
589	+	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
590	+	Pick_v2[Pick_parent],a1,b1,c1,i,
591	+	Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
592	+	Pick_v2[Pick_cnt]);
593	+	Pick_cnt++;
594	+	#endif
595	+
596	+	if(func(qtptr,orig,dir,arg1,arg2) == QT_DONE)
597	+	return(QT_DONE);
598	+
599	+	/* Advance to next node */
600	+	/* NOTE: Optimize: should only have to check 1/2 */
601	+	nbr = move_to_nbr(b,db0,db1,db2,&t);
602	+
603	+	if(nbr != -1)
604	+	{
605	+	b[0] += t * db0;
606	+	b[1] += t * db1;
607	+	b[2] += t * db2;
608	+	}
609	+	return(nbr);
610	+	}
611	+
612	+	}
613	+
614	+	int
615	+	qtRoot_trace_ray(qtptr,q0,q1,q2,orig,dir,func,arg1,arg2)
616	+	QUADTREE *qtptr;
617	+	FVECT q0,q1,q2;
618	+	FVECT orig,dir;
619	+	int (*func)();
620	+	int *arg1,arg2;
621	+	{
622	+	int i,x,y,nbr;
623	+	QUADTREE *child;
624	+	FVECT n,c,i_pt,d;
625	+	double pd,b[3],db[3],t;
626	+	/* Determine if point lies within pyramid (and therefore
627	+	inside a spherical quadtree cell):GT_INTERIOR, on one of the
628	+	pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
629	+	or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
630	+	For each triangle edge: compare the
631	+	point against the plane formed by the edge and the view center
632	+	*/
633	+	i = point_in_stri(q0,q1,q2,orig);
634	+
635	+	/* Not in this triangle */
636	+	if(!i)
637	+	return(INVALID);
638	+	/* Project the origin onto the root node plane */
639	+
640	+	/* Find the intersection point of the origin */
641	+	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
642	+	intersect_vector_plane(orig,n,pd,NULL,i_pt);
643	+	/* project the dir as well */
644	+	VADD(c,orig,dir);
645	+	intersect_vector_plane(c,n,pd,&t,c);
646	+
647	+	/* map to 2d by dropping maximum magnitude component of normal */
648	+	i = max_index(n,NULL);
649	+	x = (i+1)%3;
650	+	y = (i+2)%3;
651	+	/* Calculate barycentric coordinates of orig */
652	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b);
653	+	/* Calculate barycentric coordinates of dir */
654	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db);
655	+	if(t < 0.0)
656	+	VSUB(db,b,db);
657	+	else
658	+	VSUB(db,db,b);
659	+
660	+
661	+	#ifdef DEBUG_TEST_DRIVER
662	+	VCOPY(Pick_v0[Pick_cnt],q0);
663	+	VCOPY(Pick_v1[Pick_cnt],q1);
664	+	VCOPY(Pick_v2[Pick_cnt],q2);
665	+	Pick_cnt++;
666	+	#endif
667	+
668	+	/* trace the ray starting with this node */
669	+	nbr = qtTrace_ray(qtptr,b,db[0],db[1],db[2],orig,dir,func,arg1,arg2);
670	+	return(nbr);
671	+
672	+	}
673	+
674	+	qtVisit_tri_interior(qtptr,q0,q1,q2,t0,t1,t2,n,func,arg1,arg2,arg3)
675	+	QUADTREE *qtptr;
676	+	FVECT q0,q1,q2;
677	+	FVECT t0,t1,t2;
678	+	int n;
679	+	int (*func)();
680	+	int *arg1,arg2,*arg3;
681	+	{
682	+	int i,found,test;
683	+	QUADTREE *child;
684	+	FVECT c0,c1,c2,a,b,c;
685	+	OBJECT os[QT_MAXSET+1],*optr;
686	+	int w;
687	+
688	+	/* If qt Flag set, or qt vertices interior to t0t1t2-descend */
689	+	tree_modified:
690	+
691	+	if(QT_IS_TREE(*qtptr))
692	+	{
693	+	if(QT_IS_FLAG(*qtptr) ||  point_in_stri(t0,t1,t2,q0))
694	+	{
695	+	QT_SET_FLAG(*qtptr);
696	+	qtSubdivide_tri(q0,q1,q2,a,b,c);
697	+	/* descend to children */
698	+	for(i=0;i < 4; i++)
699	+	{
700	+	child = QT_NTH_CHILD_PTR(*qtptr,i);
701	+	qtNth_child_tri(q0,q1,q2,a,b,c,i,c0,c1,c2);
702	+	qtVisit_tri_interior(child,c0,c1,c2,t0,t1,t2,n+1,
703	+	func,arg1,arg2,arg3);
704	+	}
705	+	}
706	+	}
707	+	else
708	+	{
709	+	/* NOTE THIS IN TRI TEST Could be replaced by a flag */
710	+	if(!QT_IS_EMPTY(*qtptr))
711	+	{
712	+	if(qtinset(*qtptr,arg2))
713	+	if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2,arg3)==QT_MODIFIED)
714	+	goto tree_modified;
715	+	else
716	+	return;
717	+	}
718	+	if(point_in_stri(t0,t1,t2,q0) )
719	+	if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2,arg3)==QT_MODIFIED)
720	+	goto tree_modified;
721	+	}
722	+	}
723	+
724	+
725	+
726	+
727	+
728	+
729	+	/* NOTE: SINCE DIR could be unit: then we could use integer math */
730	+	int
731	+	qtVisit_tri_edges(qtptr,b,db0,db1,db2,
732	+	db,wptr,t,sign,sfactor,func,arg1,arg2,arg3)
733	+	QUADTREE *qtptr;
734	+	double b[3],db0,db1,db2,db[3][3];
735	+	int *wptr;
736	+	double t[3];
737	+	int sign;
738	+	double sfactor;
739	+	int (*func)();
740	+	int *arg1,arg2,*arg3;
741	+	{
742	+	int i,found;
743	+	QUADTREE *child;
744	+	int nbr,next,w;
745	+	double t_l,t_g;
746	+	#ifdef DEBUG_TEST_DRIVER
747	+	FVECT a1,b1,c1;
748	+	int Pick_parent = Pick_cnt-1;
749	+	qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
750	+	Pick_v2[Pick_parent],a1,b1,c1);
751	+	#endif
752	+	if(QT_IS_TREE(*qtptr))
753	+	{
754	+	/* Find the appropriate child and reset the coord */
755	+	i = bary_child(b);
756	+
757	+	QT_SET_FLAG(*qtptr);
758	+
759	+	for(;;)
760	+	{
761	+	w = *wptr;
762	+	child = QT_NTH_CHILD_PTR(*qtptr,i);
763	+
764	+	if(i != 3)
765	+	nbr = qtVisit_tri_edges(child,b,db0,db1,db2,
766	+	db,wptr,t,sign,
767	+	sfactor*2.0,func,arg1,arg2,arg3);
768	+	else
769	+	/* If the center cell- must flip direction signs */
770	+	nbr = qtVisit_tri_edges(child,b,-db0,-db1,-db2,
771	+	db,wptr,t,1-sign,
772	+	sfactor*2.0,func,arg1,arg2,arg3);
773	+
774	+	if(nbr == QT_DONE)
775	+	return(nbr);
776	+	if(*wptr != w)
777	+	{
778	+	w = *wptr;
779	+	db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
780	+	if(sign)
781	+	{  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
782	+	}
783	+	/* If in same block: traverse */
784	+	if(i==3)
785	+	next = nbr;
786	+	else
787	+	if(nbr == i)
788	+	next = 3;
789	+	else
790	+	{
791	+	/* reset the barycentric coordinates in the parents*/
792	+	bary_parent(b,i);
793	+	/* Else pop up to parent and traverse from there */
794	+	return(nbr);
795	+	}
796	+	bary_from_child(b,i,next);
797	+	i = next;
798	+	}
799	+	}
800	+	else
801	+	{
802	+	#ifdef DEBUG_TEST_DRIVER
803	+	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
804	+	Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
805	+	Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
806	+	Pick_cnt++;
807	+	#endif
808	+
809	+	if(func(qtptr,arg1,arg2,arg3) == QT_DONE)
810	+	return(QT_DONE);
811	+
812	+	/* Advance to next node */
813	+	w = *wptr;
814	+	while(1)
815	+	{
816	+	nbr = move_to_nbr(b,db0,db1,db2,&t_l);
817	+
818	+	t_g = t_l/sfactor;
819	+	#ifdef DEBUG
820	+	if(t[w] <= 0.0)
821	+	eputs("qtVisit_tri_edges():negative t\n");
822	+	#endif
823	+	if(t_g >= t[w])
824	+	{
825	+	if(w == 2)
826	+	return(QT_DONE);
827	+
828	+	b[0] += (t[w])*sfactor*db0;
829	+	b[1] += (t[w])*sfactor*db1;
830	+	b[2] += (t[w])*sfactor*db2;
831	+	w++;
832	+	db0 = db[w][0];
833	+	db1 = db[w][1];
834	+	db2 = db[w][2];
835	+	if(sign)
836	+	{  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
837	+	}
838	+	else
839	+	if(nbr != INVALID)
840	+	{
841	+	b[0] += t_l * db0;
842	+	b[1] += t_l * db1;
843	+	b[2] += t_l * db2;
844	+
845	+	t[w] -= t_g;
846	+	*wptr = w;
847	+	return(nbr);
848	+	}
849	+	else
850	+	return(INVALID);
851	+	}
852	+	}
853	+
854	+	}
855	+
856	+
857	+	int
858	+	qtRoot_visit_tri_edges(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2,arg3)
859	+	QUADTREE *qtptr;
860	+	FVECT q0,q1,q2;
861	+	FVECT tri[3],i_pt;
862	+	int *wptr;
863	+	int (*func)();
864	+	int *arg1,arg2,*arg3;
865	+	{
866	+	int x,y,z,nbr,w,i,j;
867	+	QUADTREE *child;
868	+	FVECT n,c,d,v[3];
869	+	double pd,b[4][3],db[3][3],et[3],t[3],exit_pt;
870	+
871	+	w = *wptr;
872	+
873	+	/* Project the origin onto the root node plane */
874	+
875	+	/* Find the intersection point of the origin */
876	+	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
877	+	/* map to 2d by dropping maximum magnitude component of normal */
878	+	z = max_index(n,NULL);
879	+	x = (z+1)%3;
880	+	y = (z+2)%3;
881	+	/* Calculate barycentric coordinates for current vertex */
882	+	if(w != -1)
883	+	{
884	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);
885	+	intersect_vector_plane(tri[w],n,pd,&(et[w]),v[w]);
886	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[w][x],v[w][y],b[w]);
887	+	}
888	+	else
889	+	/* Just starting: b[0] is the origin point: guaranteed to be valid b*/
890	+	{
891	+	w = 0;
892	+	intersect_vector_plane(tri[0],n,pd,&(et[0]),v[0]);
893	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);
894	+	VCOPY(b[3],b[0]);
895	+	}
896	+
897	+
898	+	j = (w+1)%3;
899	+	intersect_vector_plane(tri[j],n,pd,&(et[j]),v[j]);
900	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);
901	+	if(et[j] < 0.0)
902	+	{
903	+	VSUB(db[w],b[3],b[j]);
904	+	t[w] = FHUGE;
905	+	}
906	+	else
907	+	{
908	+	/* NOTE: for stability: do not increment with ipt- use full dir and
909	+	calculate t: but for wrap around case: could get same problem?
910	+	*/
911	+	VSUB(db[w],b[j],b[3]);
912	+	t[w] = 1.0;
913	+	move_to_nbr(b[3],db[w][0],db[w][1],db[w][2],&exit_pt);
914	+	if(exit_pt >= 1.0)
915	+	{
916	+	for(;j < 3;j++)
917	+	{
918	+	i = (j+1)%3;
919	+	if(i!= w)
920	+	{
921	+	intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
922	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
923	+	v[i][y],b[i]);
924	+	}
925	+	if(et[i] < 0.0)
926	+	{
927	+	VSUB(db[j],b[j],b[i]);
928	+	t[j] = FHUGE;
929	+	break;
930	+	}
931	+	else
932	+	{
933	+	VSUB(db[j],b[i],b[j]);
934	+	t[j] = 1.0;
935	+	}
936	+	move_to_nbr(b[j],db[j][0],db[j][1],db[j][2],&exit_pt);
937	+	if(exit_pt < 1.0)
938	+	break;
939	+	}
940	+	}
941	+	}
942	+	*wptr = w;
943	+	/* trace the ray starting with this node */
944	+	nbr = qtVisit_tri_edges(qtptr,b[3],db[w][0],db[w][1],db[w][2],
945	+	db,wptr,t,0,1.0,func,arg1,arg2,arg3);
946	+	if(nbr != INVALID && nbr != QT_DONE)
947	+	{
948	+	i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
949	+	i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
950	+	i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
951	+	}
952	+	return(nbr);
953	+
954	+	}
955	+
956	+	int
957	+	move_to_nbri(b,db0,db1,db2,tptr)
958	+	BCOORD b[3];
959	+	BDIR db0,db1,db2;
960	+	TINT *tptr;
961	+	{
962	+	TINT t,dt;
963	+	BCOORD bc;
964	+	int nbr;
965	+
966	+	nbr = -1;
967	+	/* Advance to next node */
968	+	if(db0 < 0)
969	+	{
970	+	bc = MAXBCOORD*b[0];
971	+	t = bc/-db0;
972	+	nbr = 0;
973	+	}
974	+	else
975	+	t = HUGET;
976	+	if(db1 < 0)
977	+	{
978	+	bc = MAXBCOORD*b[1];
979	+	dt = bc/-db1;
980	+	if( dt < t)
981	+	{
982	+	t = dt;
983	+	nbr = 1;
984	+	}
985	+	}
986	+	if(db2 < 0)
987	+	{
988	+	bc = MAXBCOORD*b[2];
989	+	dt = bc/-db2;
990	+	if( dt < t)
991	+	{
992	+	t = dt;
993	+	nbr = 2;
994	+	}
995	+	}
996	+	*tptr = t;
997	+	return(nbr);
998	+	}
999	+	/* NOTE: SINCE DIR could be unit: then we could use integer math */
1000	+	int
1001	+	qtVisit_tri_edgesi(qtptr,b,db0,db1,db2,
1002	+	db,wptr,t,sign,sfactor,func,arg1,arg2,arg3)
1003	+	QUADTREE *qtptr;
1004	+	BCOORD b[3];
1005	+	BDIR db0,db1,db2,db[3][3];
1006	+	int *wptr;
1007	+	TINT t[3];
1008	+	int sign;
1009	+	int sfactor;
1010	+	int (*func)();
1011	+	int *arg1,arg2,*arg3;
1012	+	{
1013	+	int i,found;
1014	+	QUADTREE *child;
1015	+	int nbr,next,w;
1016	+	TINT t_g,t_l,t_i;
1017	+	#ifdef DEBUG_TEST_DRIVER
1018	+	FVECT a1,b1,c1;
1019	+	int Pick_parent = Pick_cnt-1;
1020	+	qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
1021	+	Pick_v2[Pick_parent],a1,b1,c1);
1022	+	#endif
1023	+	if(QT_IS_TREE(*qtptr))
1024	+	{
1025	+	/* Find the appropriate child and reset the coord */
1026	+	i = baryi_child(b);
1027	+
1028	+	QT_SET_FLAG(*qtptr);
1029	+
1030	+	for(;;)
1031	+	{
1032	+	w = *wptr;
1033	+	child = QT_NTH_CHILD_PTR(*qtptr,i);
1034	+
1035	+	if(i != 3)
1036	+	nbr = qtVisit_tri_edgesi(child,b,db0,db1,db2,
1037	+	db,wptr,t,sign,
1038	+	sfactor+1,func,arg1,arg2,arg3);
1039	+	else
1040	+	/* If the center cell- must flip direction signs */
1041	+	nbr = qtVisit_tri_edgesi(child,b,-db0,-db1,-db2,
1042	+	db,wptr,t,1-sign,
1043	+	sfactor+1,func,arg1,arg2,arg3);
1044	+
1045	+	if(nbr == QT_DONE)
1046	+	return(nbr);
1047	+	if(*wptr != w)
1048	+	{
1049	+	w = *wptr;
1050	+	db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
1051	+	if(sign)
1052	+	{  db0 *= -1;db1 *= -1; db2 *= -1;}
1053	+	}
1054	+	/* If in same block: traverse */
1055	+	if(i==3)
1056	+	next = nbr;
1057	+	else
1058	+	if(nbr == i)
1059	+	next = 3;
1060	+	else
1061	+	{
1062	+	/* reset the barycentric coordinates in the parents*/
1063	+	baryi_parent(b,i);
1064	+	/* Else pop up to parent and traverse from there */
1065	+	return(nbr);
1066	+	}
1067	+	baryi_from_child(b,i,next);
1068	+	i = next;
1069	+	}
1070	+	}
1071	+	else
1072	+	{
1073	+	#ifdef DEBUG_TEST_DRIVER
1074	+	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
1075	+	Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
1076	+	Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
1077	+	Pick_cnt++;
1078	+	#endif
1079	+
1080	+	if(func(qtptr,arg1,arg2,arg3) == QT_DONE)
1081	+	return(QT_DONE);
1082	+
1083	+	/* Advance to next node */
1084	+	w = *wptr;
1085	+	while(1)
1086	+	{
1087	+	nbr = move_to_nbri(b,db0,db1,db2,&t_i);
1088	+
1089	+	t_g = t_i >> sfactor;
1090	+
1091	+	if(t_g >= t[w])
1092	+	{
1093	+	if(w == 2)
1094	+	return(QT_DONE);
1095	+
1096	+	/* The edge fits in the cell- therefore the result of shifting db by
1097	+	sfactor is guaranteed to be less than MAXBCOORD
1098	+	*/
1099	+	/* Caution: (t[w]*db) must occur before divide by MAXBCOORD
1100	+	since t[w] will always be < MAXBCOORD
1101	+	*/
1102	+	b[0] = b[0] + (t[w]*db0*(1 << sfactor))/MAXBCOORD;
1103	+	b[1] = b[1] + (t[w]*db1*(1 << sfactor))/MAXBCOORD;
1104	+	b[2] = b[2] + (t[w]*db2*(1 << sfactor))/MAXBCOORD;
1105	+	w++;
1106	+	db0 = db[w][0]; db1 = db[w][1]; db2 = db[w][2];
1107	+	if(sign)
1108	+	{  db0 *= -1;db1 *= -1; db2 *= -1;}
1109	+	}
1110	+	else
1111	+	if(nbr != INVALID)
1112	+	{
1113	+	/* Caution: (t_i*db) must occur before divide by MAXBCOORD
1114	+	since t_i will always be < MAXBCOORD
1115	+	*/
1116	+	b[0] = b[0] + (t_i *db0) / MAXBCOORD;
1117	+	b[1] = b[1] + (t_i *db1) / MAXBCOORD;
1118	+	b[2] = b[2] + (t_i *db2) / MAXBCOORD;
1119	+
1120	+	t[w] -= t_g;
1121	+	*wptr = w;
1122	+	return(nbr);
1123	+	}
1124	+	else
1125	+	return(INVALID);
1126	+	}
1127	+	}
1128	+	}
1129	+
1130	+
1131	+	int
1132	+	qtRoot_visit_tri_edgesi(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2,arg3)
1133	+	QUADTREE *qtptr;
1134	+	FVECT q0,q1,q2;
1135	+	FVECT tri[3],i_pt;
1136	+	int *wptr;
1137	+	int (*func)();
1138	+	int *arg1,arg2,*arg3;
1139	+	{
1140	+	int x,y,z,nbr,w,i,j;
1141	+	QUADTREE *child;
1142	+	FVECT n,c,d,v[3];
1143	+	double pd,b[4][3],db[3][3],et[3],exit_pt;
1144	+	BCOORD bi[3];
1145	+	TINT t[3];
1146	+	BDIR dbi[3][3];
1147	+	w = *wptr;
1148	+
1149	+	/* Project the origin onto the root node plane */
1150	+
1151	+	t[0] = t[1] = t[2] = 0;
1152	+	/* Find the intersection point of the origin */
1153	+	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
1154	+	/* map to 2d by dropping maximum magnitude component of normal */
1155	+	z = max_index(n,NULL);
1156	+	x = (z+1)%3;
1157	+	y = (z+2)%3;
1158	+	/* Calculate barycentric coordinates for current vertex */
1159	+	if(w != -1)
1160	+	{
1161	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);
1162	+	intersect_vector_plane(tri[w],n,pd,&(et[w]),v[w]);
1163	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[w][x],v[w][y],b[w]);
1164	+	}
1165	+	else
1166	+	/* Just starting: b[0] is the origin point: guaranteed to be valid b*/
1167	+	{
1168	+	w = 0;
1169	+	intersect_vector_plane(tri[0],n,pd,&(et[0]),v[0]);
1170	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);
1171	+	VCOPY(b[3],b[0]);
1172	+	}
1173	+
1174	+
1175	+	j = (w+1)%3;
1176	+	intersect_vector_plane(tri[j],n,pd,&(et[j]),v[j]);
1177	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);
1178	+	if(et[j] < 0.0)
1179	+	{
1180	+	VSUB(db[w],b[3],b[j]);
1181	+	t[w] = HUGET;
1182	+	}
1183	+	else
1184	+	{
1185	+	/* NOTE: for stability: do not increment with ipt- use full dir and
1186	+	calculate t: but for wrap around case: could get same problem?
1187	+	*/
1188	+	VSUB(db[w],b[j],b[3]);
1189	+	/* Check if the point is out side of the triangle: if so t[w] =HUGET */
1190	+	if((fabs(b[j][0])>(FTINY+1.0)) ||(fabs(b[j][1])>(FTINY+1.0)) ||
1191	+	(fabs(b[j][2])>(FTINY+1.0)))
1192	+	t[w] = HUGET;
1193	+	else
1194	+	{
1195	+	/* The next point is in the triangle- so db values will be in valid
1196	+	range and t[w]= MAXT
1197	+	*/
1198	+	t[w] = MAXT;
1199	+	if(j != 0)
1200	+	for(;j < 3;j++)
1201	+	{
1202	+	i = (j+1)%3;
1203	+	if(i!= w)
1204	+	{
1205	+	intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
1206	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
1207	+	v[i][y],b[i]);
1208	+	}
1209	+	if(et[i] < 0.0)
1210	+	{
1211	+	VSUB(db[j],b[j],b[i]);
1212	+	t[j] = HUGET;
1213	+	break;
1214	+	}
1215	+	else
1216	+	{
1217	+	VSUB(db[j],b[i],b[j]);
1218	+	if((fabs(b[j][0])>(FTINY+1.0)) ||
1219	+	(fabs(b[j][1])>(FTINY+1.0)) ||
1220	+	(fabs(b[j][2])>(FTINY+1.0)))
1221	+	{
1222	+	t[j] = HUGET;
1223	+	break;
1224	+	}
1225	+	else
1226	+	t[j] = MAXT;
1227	+	}
1228	+	}
1229	+	}
1230	+	}
1231	+	*wptr = w;
1232	+	bary_dtol(b[3],db,bi,dbi,t);
1233	+
1234	+	/* trace the ray starting with this node */
1235	+	nbr = qtVisit_tri_edgesi(qtptr,bi,dbi[w][0],dbi[w][1],dbi[w][2],
1236	+	dbi,wptr,t,0,0,func,arg1,arg2,arg3);
1237	+	if(nbr != INVALID && nbr != QT_DONE)
1238	+	{
1239	+	b[3][0] = (double)bi[0]/(double)MAXBCOORD;
1240	+	b[3][1] = (double)bi[1]/(double)MAXBCOORD;
1241	+	b[3][2] = (double)bi[2]/(double)MAXBCOORD;
1242	+	i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
1243	+	i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
1244	+	i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
1245	+	}
1246	+	return(nbr);
1247	+
1248	+	}
1249	+
1250	+
1251	+
1252	+
1253
1254
1255

Diff Legend

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