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.1 by gwlarson, Wed Aug 19 17:45:24 1998 UTC vs.
Revision 3.4 by gwlarson, Fri Sep 11 11:52:26 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= 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
27	–
32		QUADTREE
33		qtAlloc()                       /* allocate a quadtree */
34		{
#	Line 41 | Line 45 | qtAlloc()                      /* allocate a quadtree */
45		if (QT_BLOCK(freet) >= QT_MAX_BLK)
46		return(EMPTY);
47		if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
48	<	(unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
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));
53	+	if (quad_flag == NULL)
54	+	return(EMPTY);
55		}
56		treetop += 4;
57		return(freet);
58		}
59
60
61	+	qtClearAllFlags()               /* clear all quadtree branch flags */
62	+	{
63	+	if (!treetop) return;
64	+	bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*(QT_BLOCK_SIZE/8));
65	+	}
66	+
67	+
68		qtFree(qt)                      /* free a quadtree */
69		register QUADTREE  qt;
70		{
#	Line 70 | Line 86 | qtDone()                       /* free EVERYTHING */
86		{
87		register int    i;
88
89	+	qtfreeleaves();
90		for (i = 0; i < QT_MAX_BLK; i++)
91		{
92	<	free((char *)quad_block[i],
93	<	(unsigned)QT_BLOCK_SIZE*4*sizeof(QUADTREE));
92	>	if (quad_block[i] == NULL)
93	>	break;
94	>	free((char *)quad_block[i]);
95		quad_block[i] = NULL;
96		}
97	+	if (i) free((char *)quad_flag);
98	+	quad_flag = NULL;
99		quad_free_list = EMPTY;
100		treetop = 0;
101		}
102
103		QUADTREE
104	<	qtCompress(qt)                  /* recursively combine nodes */
105	<	register QUADTREE  qt;
104	>	*qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
105	>	QUADTREE *qtptr;
106	>	double bcoord[3];
107	>	FVECT t0,t1,t2;
108		{
109	<	register int  i;
110	<	register QUADTREE  qres;
109	>	int i;
110	>	QUADTREE *child;
111	>	FVECT a,b,c;
112
113	<	if (!QT_IS_TREE(qt))        /* not a tree */
114	<	return(qt);
115	<	qres = QT_NTH_CHILD(qt,0) = qtCompress(QT_NTH_CHILD(qt,0));
116	<	for (i = 1; i < 4; i++)
117	<	if((QT_NTH_CHILD(qt,i) = qtCompress(QT_NTH_CHILD(qt,i))) != qres)
118	<	qres = qt;
119	<	if(!QT_IS_TREE(qres))
120	<	{   /* all were identical leaves */
121	<	QT_NTH_CHILD(qt,0) = quad_free_list;
122	<	quad_free_list = qt;
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	<	return(qres);
134	>	else
135	>	return(qtptr);
136		}
137
138	+
139		QUADTREE
140	<	qtPoint_locate(qtptr,v1,v2,v3,pt,type,which,p0,p1,p2)
140	>	*qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
141		QUADTREE *qtptr;
142	<	FVECT v1,v2,v3;
142	>	FVECT v0,v1,v2;
143		FVECT pt;
144	<	char *type,*which;
110	<	FVECT p0,p1,p2;
144	>	FVECT t0,t1,t2;
145		{
146	<	char d,w;
147	<	int i;
146	>	int d;
147	>	int i,x,y;
148		QUADTREE *child;
149	<	QUADTREE qt;
150	<	FVECT a,b,c;
117	<	FVECT t1,t2,t3;
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
#	Line 123 | Line 156 | FVECT p0,p1,p2;
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(v1,v2,v3,pt,&w);
159	>	d = point_in_stri(v0,v1,v2,pt);
160
161	+
162		/* Not in this triangle */
163		if(!d)
164	<	{
131	<	if(which)
132	<	*which = 0;
133	<	return(EMPTY);
134	<	}
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		*/
141	–
171		if(QT_IS_TREE(*qtptr))
172		{
173	<	qtSubdivide_tri(v1,v2,v3,a,b,c);
174	<	child = QT_NTH_CHILD_PTR(*qtptr,0);
175	<	if(!QT_IS_EMPTY(*child))
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	<	qt = qtPoint_locate(child,v1,a,c,pt,type,which,p0,p1,p2);
202	<	if(!QT_IS_EMPTY(qt))
150	<	return(qt);
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	<	child = QT_NTH_CHILD_PTR(*qtptr,1);
153	<	if(!QT_IS_EMPTY(*child))
154	<	{
155	<	qt = qtPoint_locate(child,a,b,c,pt,type,which,p0,p1,p2);
156	<	if(!QT_IS_EMPTY(qt))
157	<	return(qt);
158	<	}
159	<	child = QT_NTH_CHILD_PTR(*qtptr,2);
160	<	if(!QT_IS_EMPTY(*child))
161	<	{
162	<	qt = qtPoint_locate(child,a,v2,b,pt,type,which,p0,p1,p2);
163	<	if(!QT_IS_EMPTY(qt))
164	<	return(qt);
165	<	}
166	<	child = QT_NTH_CHILD_PTR(*qtptr,3);
167	<	if(!QT_IS_EMPTY(*child))
168	<	{
169	<	qt = qtPoint_locate(child,c,b,v3,pt,type,which,p0,p1,p2);
170	<	if(!QT_IS_EMPTY(qt))
171	<	return(qt);
172	<	}
204	>	return(qtLocate_leaf(child,bcoord,t0,t1,t2));
205		}
206		else
175	–	if(!QT_IS_EMPTY(*qtptr))
176	–	{
177	–	/* map GT_VERTEX,GT_EDGE,GT_FACE GT_INTERIOR of pyramid to
178	–	spherical triangle primitives
179	–	*/
180	–	if(type)
181	–	*type = d;
182	–	if(which)
183	–	*which = w;
184	–	VCOPY(p0,v1);
185	–	VCOPY(p1,v2);
186	–	VCOPY(p2,v3);
187	–	return(*qtptr);
188	–	}
189	–	return(EMPTY);
190	–	}
191	–
192	–	int
193	–	qtPoint_in_tri(qtptr,v0,v1,v2,pt,type,which)
194	–	QUADTREE *qtptr;
195	–	FVECT v0,v1,v2;
196	–	FVECT pt;
197	–	char *type,*which;
198	–	{
199	–	QUADTREE qt;
200	–	OBJECT os[MAXSET+1],*optr;
201	–	char d,w;
202	–	int i,id;
203	–	FVECT p0,p1,p2;
204	–
205	–	qt = qtPoint_locate(qtptr,v0,v1,v2,pt,type,which,p0,p1,p2);
206	–	if(QT_IS_EMPTY(qt))
207	–	return(EMPTY);
208	–
209	–	/* Get the set */
210	–	qtgetset(os,qt);
211	–	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)
218	<	{
219	<	if(type)
220	<	*type = d;
221	<	if(which)
222	<	*which = w;
223	<	return(id);
208	>	if(t0)
209	>	{
210	>	VCOPY(t0,v0);
211	>	VCOPY(t1,v1);
212	>	VCOPY(t2,v2);
213		}
214	+	return(qtptr);
215		}
226	–	return(EMPTY);
216		}
217
218	+
219		QUADTREE
220		qtSubdivide(qtptr)
221		QUADTREE *qtptr;
#	Line 250 | Line 240 | int n;
240		return(node);
241		}
242
243	<	/* for triangle v1-v2-v3- returns a,b,c: children are:
243	>	/* for triangle v0-v1-v2- returns a,b,c: children are:
244
245	<	v3                     0: v1,a,c
246	<	/\                     1: a,b,c
247	<	/3 \                    2: a,v2,b
248	<	c/____\b                  3: c,b,v3
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 \1 /2 \
251	<	v1/____\/____\v2
250	>	/0 \3 /1 \
251	>	v0____\/____\v1
252		a
253		*/
254
255	<	qtSubdivide_tri(v1,v2,v3,a,b,c)
256	<	FVECT v1,v2,v3;
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,v1,v2);
260	<	normalize(a);
261	<	EDGE_MIDPOINT_VEC3(b,v2,v3);
272	<	normalize(b);
273	<	EDGE_MIDPOINT_VEC3(c,v3,v1);
274	<	normalize(c);
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(v1,v2,v3,a,b,c,i,r1,r2,r3)
265	<	FVECT v1,v2,v3;
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 r1,r2,r3;
268	>	FVECT r0,r1,r2;
269		{
270	<	VCOPY(r1,a); VCOPY(r2,b);    VCOPY(r3,c);
270	>
271		switch(i){
272		case 0:
273	<	VCOPY(r2,r1);
287	<	VCOPY(r1,v1);
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(r3,r2);
293	<	VCOPY(r2,v2);
279	>	VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
280		break;
281		case 3:
282	<	VCOPY(r1,r3);
297	<	VCOPY(r3,v3);
282	>	VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
283		break;
284		}
285		}
#	Line 308 | Line 293 | into the new child cells: it is assumed that "v1,v2,v3
293		*/
294
295		int
296	<	qtAdd_tri(qtptr,id,t1,t2,t3,v1,v2,v3,n)
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;
314	–	FVECT t1,t2,t3;
315	–	FVECT v1,v2,v3;
301		int n;
302		{
303	<
319	<	char test;
320	<	int i,index;
321	<	FVECT a,b,c;
322	<	OBJECT os[MAXSET+1],*optr;
323	<	QUADTREE qt;
303	>	int test;
304		int found;
325	–	FVECT r1,r2,r3;
305
306	<	/* test if triangle (t1,t2,t3) overlaps cell triangle (v1,v2,v3) */
328	<	test = spherical_tri_intersect(t1,t2,t3,v1,v2,v3);
329	<
330	<	/* If triangles do not overlap: done */
306	>	test = stri_intersect(q0,q1,q2,t0,t1,t2);
307		if(!test)
308		return(FALSE);
309	<	found = 0;
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(v1,v2,v3,a,b,c);
335	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t1,t2,t3,v1,a,c,n);
336	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t1,t2,t3,a,b,c,n);
337	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t1,t2,t3,a,v2,b,n);
338	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t1,t2,t3,c,b,v3,n);
339	<
340	<	#if 0
341	<	if(!found)
342	<	{
343	<	#ifdef TEST_DRIVER
344	<	HANDLE_ERROR("qtAdd_tri():Found in parent but not children\n");
345	<	#else
346	<	eputs("qtAdd_tri():Found in parent but not children\n");
352	<	#endif
353	<	}
354	<	#endif
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	<	{
361	<	*qtptr = qtaddelem(*qtptr,id);
362	<	#if 0
363	<	{
364	<	int k;
365	<	qtgetset(os,*qtptr);
366	<	printf("\n%d:\n",os[0]);
367	<	for(k=1; k <= os[0];k++)
368	<	printf("%d ",os[k]);
369	<	printf("\n");
370	<	}
371	<	#endif
372	<	/*
373	<	os[0] = 0;
374	<	insertelem(os,id);
375	<	qt = fullnode(os);
376	<	*qtptr = qt;
377	<	*/
378	<	}
352	>	*qtptr = qtaddelem(*qtptr,id);
353		else
354		{
381	–	qtgetset(os,*qtptr);
355		/* If the set is too large: subdivide */
356	<	if(QT_SET_CNT(os) < QT_SET_THRESHOLD)
357	<	{
358	<	*qtptr = qtaddelem(*qtptr,id);
359	<	#if 0
387	<	{
388	<	int k;
389	<	qtgetset(os,*qtptr);
390	<	printf("\n%d:\n",os[0]);
391	<	for(k=1; k <= os[0];k++)
392	<	printf("%d ",os[k]);
393	<	printf("\n");
394	<	}
395	<	#endif
396	<	/*
397	<	insertelem(os,id);
398	<	*qtptr = fullnode(os);
399	<	*/
400	<	}
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 405 | 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,t1,t2,t3,v1,v2,v3,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
420	<	}
421	<	#endif
422	<	for(optr = &(os[1]),i = QT_SET_CNT(os); i > 0; i--)
423	<	{
424	<	id = QT_SET_NEXT_ELEM(optr);
425	<	qtTri_verts_from_id(id,r1,r2,r3);
426	<	found=qtAdd_tri(qtptr,id,r1,r2,r3,v1,v2,v3,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,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)
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)
435	<	{
386	>	if(QT_SET_CNT(optr) < QT_MAXSET)
387		*qtptr = qtaddelem(*qtptr,id);
437	–	#if 0
438	–	{
439	–	int k;
440	–	qtgetset(os,*qtptr);
441	–	printf("\n%d:\n",os[0]);
442	–	for(k=1; k <= os[0];k++)
443	–	printf("%d ",os[k]);
444	–	printf("\n");
445	–	}
446	–	#endif
447	–	/*
448	–	insertelem(os,id);
449	–	*qtptr = fullnode(os);
450	–	*/
451	–	}
388		else
389		{
390		#ifdef DEBUG
#	Line 464 | Line 400 | int n;
400
401
402		int
403	<	qtApply_to_tri_cells(qtptr,t1,t2,t3,v1,v2,v3,func,arg)
403	>	qtApply_to_tri_cells(qtptr,t0,t1,t2,v0,v1,v2,func,arg)
404		QUADTREE *qtptr;
405	<	FVECT t1,t2,t3;
406	<	FVECT v1,v2,v3;
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 (t1,t2,t3) overlaps cell triangle (v1,v2,v3) */
414	<	test = spherical_tri_intersect(t1,t2,t3,v1,v2,v3);
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	<	qtSubdivide_tri(v1,v2,v3,a,b,c);
426	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,0),t1,t2,t3,v1,a,c,func,arg);
427	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,1),t1,t2,t3,a,b,c,func,arg);
428	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,2),t1,t2,t3,a,v2,b,func,arg);
429	<	qtApply_to_tri_cells(QT_NTH_CHILD_PTR(*qtptr,3),t1,t2,t3,c,b,v3,func,arg);
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		}
493	–	else
494	–	return(func(qtptr,arg));
432		}
433
497	–
434		int
435	<	qtRemove_tri(qtptr,id,t1,t2,t3,v1,v2,v3)
435	>	qtRemove_tri(qtptr,id,t0,t1,t2,v0,v1,v2)
436		QUADTREE *qtptr;
437		int id;
438	<	FVECT t1,t2,t3;
439	<	FVECT v1,v2,v3;
438	>	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 (t1,t2,t3) overlaps cell triangle (v1,v2,v3) */
448	<	test = spherical_tri_intersect(t1,t2,t3,v1,v2,v3);
447	>	/* test if triangle (t0,t1,t2) overlaps cell triangle (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 518 | Line 454 | FVECT v1,v2,v3;
454		/* if this is tree: recurse */
455		if(QT_IS_TREE(*qtptr))
456		{
457	<	qtSubdivide_tri(v1,v2,v3,a,b,c);
458	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t1,t2,t3,v1,a,c);
459	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t1,t2,t3,a,b,c);
460	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t1,t2,t3,a,v2,b);
461	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t1,t2,t3,c,b,v3);
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);
462		}
463		else
464		{
#	Line 535 | Line 471 | FVECT v1,v2,v3;
471		/* remove id from set */
472		else
473		{
474	<	qtgetset(os,*qtptr);
539	<	if(!inset(os,id))
474	>	if(!qtinset(*qtptr,id))
475		{
476		#ifdef DEBUG
477		eputs("qtRemove_tri(): tri not in set\n");
#	Line 545 | Line 480 | FVECT v1,v2,v3;
480		else
481		{
482		*qtptr = qtdelelem(*qtptr,id);
483	<	#if 0
484	<	{
485	<	int k;
486	<	if(!QT_IS_EMPTY(*qtptr))
487	<	{qtgetset(os,*qtptr);
553	<	printf("\n%d:\n",os[0]);
554	<	for(k=1; k <= os[0];k++)
555	<	printf("%d ",os[k]);
556	<	printf("\n");
557	<	}
483	>	}
484	>	}
485	>	}
486	>	return(TRUE);
487	>	}
488
489	<	}
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_edge(qtptr,b,db,orig,dir,max_t,func,arg1,arg2)
531	>	QUADTREE *qtptr;
532	>	double b[3],db[3];
533	>	FVECT orig,dir;
534	>	double max_t;
535	>	int (*func)();
536	>	int *arg1,arg2;
537	>	{
538	>
539	>	int i,found;
540	>	QUADTREE *child;
541	>	int nbr,next;
542	>	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);
549	>
550		#endif
551	+	if(QT_IS_TREE(*qtptr))
552	+	{
553	+	/* Find the appropriate child and reset the coord */
554	+	i = bary_child(b);
555	+
556	+	QT_SET_FLAG(*qtptr);
557	+
558	+	for(;;)
559	+	{
560	+	child = QT_NTH_CHILD_PTR(*qtptr,i);
561	+
562	+	if(i != 3)
563	+	{
564	+
565	+	db[0] *= 2.0;db[1] *= 2.0; db[2] *= 2.0;
566	+	nbr = qtTrace_edge(child,b,db,orig,dir,max_t,func,arg1,arg2);
567	+	db[0] *= 0.5;db[1] *= 0.5; db[2] *= 0.5;
568	+	}
569	+	else
570	+	{
571	+	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;
594	+	}
595	+	}
596	+	else
597	+	{
598	+	#ifdef DEBUG_TEST_DRIVER
599	+	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
600	+	Pick_v2[Pick_parent],a1,b1,c1,i,
601	+	Pick_v0[Pick_cnt],Pick_v1[Pick_cnt],
602	+	Pick_v2[Pick_cnt]);
603	+	Pick_cnt++;
604	+	#endif
605	+
606	+	if(func(qtptr,arg1,arg2) == QT_DONE)
607	+	return(QT_DONE);
608	+
609	+	/* Advance to next node */
610	+	/* NOTE: Optimize: should only have to check 1/2 */
611	+	nbr = move_to_nbr(b,db[0],db[1],db[2],&t);
612	+
613	+	if(t >= max_t)
614	+	return(QT_DONE);
615	+	if(nbr != -1)
616	+	{
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);
625	+	}
626	+
627	+	}
628	+
629	+
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	+	{
638	+
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);
649	+
650	+	#endif
651	+	if(QT_IS_TREE(*qtptr))
652	+	{
653	+	/* Find the appropriate child and reset the coord */
654	+	i = bary_child(b);
655	+
656	+	QT_SET_FLAG(*qtptr);
657	+
658	+	for(;;)
659	+	{
660	+	child = QT_NTH_CHILD_PTR(*qtptr,i);
661	+
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);
669	+
670	+	/* If in same block: traverse */
671	+	if(i==3)
672	+	next = nbr;
673	+	else
674	+	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
696	+
697	+	if(func(qtptr,orig,dir,arg1,arg2) == QT_DONE)
698	+	return(QT_DONE);
699	+
700	+	/* Advance to next node */
701	+	/* NOTE: Optimize: should only have to check 1/2 */
702	+	nbr = move_to_nbr(b,db0,db1,db2,&t);
703	+
704	+	if(nbr != -1)
705	+	{
706	+	b[0] += t * db0;
707	+	b[1] += t * db1;
708	+	b[2] += t * db2;
709	+	}
710	+	return(nbr);
711	+	}
712	+
713	+	}
714	+
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 */
740	+
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);
747	+
748	+	/* map to 2d by dropping maximum magnitude component of normal */
749	+	i = max_index(n);
750	+	x = (i+1)%3;
751	+	y = (i+2)%3;
752	+	/* Calculate barycentric coordinates of orig */
753	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b);
754	+	/* Calculate barycentric coordinates of dir */
755	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db);
756	+	if(t < 0.0)
757	+	VSUB(db,b,db);
758	+	else
759	+	VSUB(db,db,b);
760	+
761	+
762	+	#ifdef DEBUG_TEST_DRIVER
763	+	VCOPY(Pick_v0[Pick_cnt],q0);
764	+	VCOPY(Pick_v1[Pick_cnt],q1);
765	+	VCOPY(Pick_v2[Pick_cnt],q2);
766	+	Pick_cnt++;
767	+	#endif
768	+
769	+	/* trace the ray starting with this node */
770	+	nbr = qtTrace_ray(qtptr,b,db[0],db[1],db[2],orig,dir,func,arg1,arg2);
771	+	return(nbr);
772	+
773	+	}
774	+
775	+
776	+	int
777	+	qtRoot_trace_edge(qtptr,q0,q1,q2,orig,dir,max_t,func,arg1,arg2)
778	+	QUADTREE *qtptr;
779	+	FVECT q0,q1,q2;
780	+	FVECT orig,dir;
781	+	double max_t;
782	+	int (*func)();
783	+	int *arg1,arg2;
784	+	{
785	+	int i,x,y,nbr;
786	+	QUADTREE *child;
787	+	FVECT n,c,i_pt,d;
788	+	double pd,b[3],db[3],t;
789	+	/* Determine if point lies within pyramid (and therefore
790	+	inside a spherical quadtree cell):GT_INTERIOR, on one of the
791	+	pyramid sides (and on cell edge):GT_EDGE(1,2 or 3),
792	+	or on pyramid vertex (and on cell vertex):GT_VERTEX(1,2, or 3).
793	+	For each triangle edge: compare the
794	+	point against the plane formed by the edge and the view center
795	+	*/
796	+	i = point_in_stri(q0,q1,q2,orig);
797	+
798	+	/* Not in this triangle */
799	+	if(!i)
800	+	return(-1);
801	+	/* Project the origin onto the root node plane */
802	+
803	+	/* 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);
809	+
810	+	/* map to 2d by dropping maximum magnitude component of normal */
811	+	i = max_index(n);
812	+	x = (i+1)%3;
813	+	y = (i+2)%3;
814	+	/* Calculate barycentric coordinates of orig */
815	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b);
816	+	/* Calculate barycentric coordinates of dir */
817	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],c[x],c[y],db);
818	+	if(t < 0.0)
819	+	VSUB(db,b,db);
820	+	else
821	+	VSUB(db,db,b);
822	+
823	+
824	+	#ifdef DEBUG_TEST_DRIVER
825	+	VCOPY(Pick_v0[Pick_cnt],q0);
826	+	VCOPY(Pick_v1[Pick_cnt],q1);
827	+	VCOPY(Pick_v2[Pick_cnt],q2);
828	+	Pick_cnt++;
829	+	#endif
830	+	/* trace the ray starting with this node */
831	+	nbr = qtTrace_edge(qtptr,b,db,orig,d,max_t,func,arg1,arg2);
832	+	return(nbr);
833	+
834	+	}
835	+
836	+
837	+	qtVisit_tri_interior(qtptr,q0,q1,q2,t0,t1,t2,n,func,arg1,arg2)
838	+	QUADTREE *qtptr;
839	+	FVECT q0,q1,q2;
840	+	FVECT t0,t1,t2;
841	+	int n;
842	+	int (*func)();
843	+	int *arg1,arg2;
844	+	{
845	+	int i,found,test;
846	+	QUADTREE *child;
847	+	FVECT c0,c1,c2,a,b,c;
848	+	OBJECT os[QT_MAXSET+1],*optr;
849	+	int w;
850	+
851	+	/* If qt Flag set, or qt vertices interior to t0t1t2-descend */
852	+	tree_modified:
853	+
854	+	if(QT_IS_TREE(*qtptr))
855	+	{
856	+	if(QT_IS_FLAG(*qtptr) ||  point_in_stri(t0,t1,t2,q0))
857	+	{
858	+	QT_SET_FLAG(*qtptr);
859	+	qtSubdivide_tri(q0,q1,q2,a,b,c);
860	+	/* descend to children */
861	+	for(i=0;i < 4; i++)
862	+	{
863	+	child = QT_NTH_CHILD_PTR(*qtptr,i);
864	+	qtNth_child_tri(q0,q1,q2,a,b,c,i,c0,c1,c2);
865	+	qtVisit_tri_interior(child,c0,c1,c2,t0,t1,t2,n+1,
866	+	func,arg1,arg2);
867	+	}
868		}
869	+	}
870	+	else
871	+	{
872	+	/* NOTE THIS IN TRI TEST Could be replaced by a flag */
873	+	if(!QT_IS_EMPTY(*qtptr))
874	+	{
875	+	if(qtinset(*qtptr,arg2))
876	+	if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2)==QT_MODIFIED)
877	+	goto tree_modified;
878	+	else
879	+	return;
880		}
881	+	if(point_in_stri(t0,t1,t2,q0) )
882	+	if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2)==QT_MODIFIED)
883	+	goto tree_modified;
884	+	}
885	+	}
886	+
887	+
888	+
889	+
890	+
891	+
892	+	int
893	+	qtVisit_tri_edges(qtptr,b,db,wptr,sfactor,func,arg1,arg2)
894	+	QUADTREE *qtptr;
895	+	double b[3],db[3][3];
896	+	int *wptr;
897	+	double sfactor;
898	+	int (*func)();
899	+	int *arg1,arg2;
900	+	{
901	+	int i,found;
902	+	QUADTREE *child;
903	+	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
911	+
912	+	if(QT_IS_TREE(*qtptr))
913	+	{
914	+	/* Find the appropriate child and reset the coord */
915	+	i = bary_child(b);
916	+
917	+	QT_SET_FLAG(*qtptr);
918	+
919	+	for(;;)
920	+	{
921	+	w = *wptr;
922	+	child = QT_NTH_CHILD_PTR(*qtptr,i);
923	+
924	+	if(i != 3)
925	+	{
926	+
927	+	db[w][0] *= 2.0;db[w][1] *= 2.0; db[w][2] *= 2.0;
928	+	nbr = qtVisit_tri_edges(child,b,db,wptr,sfactor*2.0,
929	+	func,arg1,arg2);
930	+	w = *wptr;
931	+	db[w][0] *= 0.5;db[w][1] *= 0.5; db[w][2] *= 0.5;
932	+	}
933	+	else
934	+	{
935	+	db[w][0] *=-2.0;db[w][1] *= -2.0; db[w][2] *= -2.0;
936	+	/* If the center cell- must flip direction signs */
937	+	nbr = qtVisit_tri_edges(child,b,db,wptr,sfactor*(-2.0),
938	+	func,arg1,arg2);
939	+	w = *wptr;
940	+	db[w][0] *=-0.5;db[w][1] *= -0.5; db[w][2] *= -0.5;
941	+	}
942	+	if(nbr == QT_DONE)
943	+	return(nbr);
944	+
945	+	/* If in same block: traverse */
946	+	if(i==3)
947	+	next = nbr;
948	+	else
949	+	if(nbr == i)
950	+	next = 3;
951	+	else
952	+	{
953	+	/* reset the barycentric coordinates in the parents*/
954	+	bary_parent(b,i);
955	+	/* Else pop up to parent and traverse from there */
956	+	return(nbr);
957	+	}
958	+	bary_from_child(b,i,next);
959	+	i = next;
960	+	}
961	+	}
962	+	else
963	+	{
964	+	#ifdef DEBUG_TEST_DRIVER
965	+	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
966	+	Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
967	+	Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
968	+	Pick_cnt++;
969	+	#endif
970	+
971	+	if(func(qtptr,arg1,arg2) == QT_DONE)
972	+	return(QT_DONE);
973	+
974	+	/* Advance to next node */
975	+	/* NOTE: Optimize: should only have to check 1/2 */
976	+	w = *wptr;
977	+	while(1)
978	+	{
979	+	nbr = move_to_nbr(b,db[w][0],db[w][1],db[w][2],&t);
980	+
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	+
1010	+	}
1011	+
1012	+
1013	+	int
1014	+	qtRoot_visit_tri_edges(qtptr,q0,q1,q2,tri,dir,wptr,func,arg1,arg2)
1015	+	QUADTREE *qtptr;
1016	+	FVECT q0,q1,q2;
1017	+	FVECT tri[3],dir[3];
1018	+	int *wptr;
1019	+	int (*func)();
1020	+	int *arg1,arg2;
1021	+	{
1022	+	int i,x,y,nbr,w;
1023	+	QUADTREE *child;
1024	+	FVECT n,c,i_pt,d;
1025	+	double pd,b[3][3],db[3][3],t;
1026	+
1027	+	w = *wptr;
1028	+
1029	+	/* Project the origin onto the root node plane */
1030	+
1031	+	/* Find the intersection point of the origin */
1032	+	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
1033	+	/* map to 2d by dropping maximum magnitude component of normal */
1034	+	i = max_index(n);
1035	+	x = (i+1)%3;
1036	+	y = (i+2)%3;
1037	+	/* Calculate barycentric coordinates for current vertex */
1038	+
1039	+	for(i=0;i < 3; i++)
1040	+	{
1041	+	/* If processing 3rd edge-dont need info for t1 */
1042	+	if(i==1 && w==2)
1043	+	continue;
1044	+	/* project the dir as well */
1045	+	intersect_vector_plane(tri[i],n,pd,NULL,i_pt);
1046	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[i]);
1047	+	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	<	return(TRUE);
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		}
1067	+
1068	+
1069	+
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))
1095	+	{
1096	+	/* Find the appropriate child and reset the coord */
1097	+	i = bary_child(b);
1098	+
1099	+	QT_SET_FLAG(*qtptr);
1100	+
1101	+	for(;;)
1102	+	{
1103	+	w = *wptr;
1104	+	child = QT_NTH_CHILD_PTR(*qtptr,i);
1105	+
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);
1115	+
1116	+	if(nbr == QT_DONE)
1117	+	return(nbr);
1118	+	if(*wptr != w)
1119	+	{
1120	+	w = *wptr;
1121	+	db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
1122	+	if(sign)
1123	+	{  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
1124	+	}
1125	+	/* If in same block: traverse */
1126	+	if(i==3)
1127	+	next = nbr;
1128	+	else
1129	+	if(nbr == i)
1130	+	next = 3;
1131	+	else
1132	+	{
1133	+	/* reset the barycentric coordinates in the parents*/
1134	+	bary_parent(b,i);
1135	+	/* Else pop up to parent and traverse from there */
1136	+	return(nbr);
1137	+	}
1138	+	bary_from_child(b,i,next);
1139	+	i = next;
1140	+	}
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
1150	+
1151	+	if(func(qtptr,arg1,arg2) == QT_DONE)
1152	+	return(QT_DONE);
1153	+
1154	+	/* Advance to next node */
1155	+	w = *wptr;
1156	+	while(1)
1157	+	{
1158	+	nbr = move_to_nbr(b,db0,db1,db2,&t_l);
1159	+
1160	+	t_g = t_l/sfactor;
1161	+	#ifdef DEBUG
1162	+	if(t[w] <= 0.0)
1163	+	eputs("qtVisit_tri_edges2():negative t\n");
1164	+	#endif
1165	+	if(t_g >= t[w])
1166	+	{
1167	+	if(w == 2)
1168	+	return(QT_DONE);
1169	+
1170	+	b[0] += (t[w])*sfactor*db0;
1171	+	b[1] += (t[w])*sfactor*db1;
1172	+	b[2] += (t[w])*sfactor*db2;
1173	+	w++;
1174	+	db0 = db[w][0];
1175	+	db1 = db[w][1];
1176	+	db2 = db[w][2];
1177	+	if(sign)
1178	+	{  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
1179	+	}
1180	+	else
1181	+	if(nbr != INVALID)
1182	+	{
1183	+	b[0] += t_l * db0;
1184	+	b[1] += t_l * db1;
1185	+	b[2] += t_l * db2;
1186	+
1187	+	t[w] -= t_g;
1188	+	*wptr = w;
1189	+	return(nbr);
1190	+	}
1191	+	else
1192	+	return(INVALID);
1193	+	}
1194	+	}
1195	+
1196	+	}
1197	+
1198	+
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;
1214	+
1215	+	/* Project the origin onto the root node plane */
1216	+
1217	+	/* Find the intersection point of the origin */
1218	+	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
1219	+	/* map to 2d by dropping maximum magnitude component of normal */
1220	+	z = max_index(n);
1221	+	x = (z+1)%3;
1222	+	y = (z+2)%3;
1223	+	/* Calculate barycentric coordinates for current vertex */
1224	+	if(w != -1)
1225	+	{
1226	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);
1227	+	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]);
1229	+	}
1230	+	else
1231	+	/* Just starting: b[0] is the origin point: guaranteed to be valid b*/
1232	+	{
1233	+	w = 0;
1234	+	intersect_vector_plane(tri[0],n,pd,&(et[0]),v[0]);
1235	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);
1236	+	VCOPY(b[3],b[0]);
1237	+	}
1238	+
1239	+
1240	+	j = (w+1)%3;
1241	+	intersect_vector_plane(tri[j],n,pd,&(et[j]),v[j]);
1242	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);
1243	+	if(et[j] < 0.0)
1244	+	{
1245	+	VSUB(db[w],b[3],b[j]);
1246	+	t[w] = FHUGE;
1247	+	}
1248	+	else
1249	+	{
1250	+	VSUB(db[w],b[j],b[3]);
1251	+	t[w] = 1.0;
1252	+	move_to_nbr(b[3],db[w][0],db[w][1],db[w][2],&exit_pt);
1253	+	if(exit_pt >= 1.0)
1254	+	{
1255	+	for(;j < 3;j++)
1256	+	{
1257	+	i = (j+1)%3;
1258	+	if(i!= w)
1259	+	{
1260	+	intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
1261	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
1262	+	v[i][y],b[i]);
1263	+	}
1264	+	if(et[i] < 0.0)
1265	+	{
1266	+	VSUB(db[j],b[j],b[i]);
1267	+	t[j] = FHUGE;
1268	+	break;
1269	+	}
1270	+	else
1271	+	{
1272	+	VSUB(db[j],b[i],b[j]);
1273	+	t[j] = 1.0;
1274	+	}
1275	+	move_to_nbr(b[j],db[j][0],db[j][1],db[j][2],&exit_pt);
1276	+	if(exit_pt < 1.0)
1277	+	break;
1278	+	}
1279	+	}
1280	+	}
1281	+	*wptr = w;
1282	+	/* trace the ray starting with this node */
1283	+	nbr = qtVisit_tri_edges2(qtptr,b[3],db[w][0],db[w][1],db[w][2],
1284	+	db,wptr,t,0,1.0,func,arg1,arg2);
1285	+	if(nbr != INVALID && nbr != QT_DONE)
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	+	}
1294	+
1295	+
1296	+
1297	+
1298	+
1299	+
1300	+
1301	+
1302	+
1303	+
1304	+
1305	+
1306	+
1307	+
1308	+
1309	+
1310	+
1311	+

Diff Legend

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