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.8 by gwlarson, Wed Nov 11 12:05:39 1998 UTC

#	Line 14 | Line 14 | static char SCCSid[] = "$SunId$ SGI";
14		*/
15
16		#include "standard.h"
17	<
17	>	#include "sm_flag.h"
18		#include "sm_geom.h"
19		#include "sm_qtree.h"
20
#	Line 27 | Line 27 | int4 *quad_flag= NULL;
27		extern FVECT Pick_v0[500],Pick_v1[500],Pick_v2[500];
28		extern int Pick_cnt,Pick_tri,Pick_samp;
29		extern FVECT Pick_point[500];
30	+	extern int Pick_q[500];
31	+
32		#endif
33	+	int Incnt=0;
34
35		QUADTREE
36		qtAlloc()                       /* allocate a quadtree */
#	Line 46 | Line 49 | qtAlloc()                      /* allocate a quadtree */
49		return(EMPTY);
50		if ((quad_block[QT_BLOCK(freet)] = (QUADTREE *)malloc(
51		QT_BLOCK_SIZE*4*sizeof(QUADTREE))) == NULL)
52	<	return(EMPTY);
52	>	error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
53
54	+	/* Realloc the per/node flags */
55		quad_flag = (int4 *)realloc((char *)quad_flag,
56	<	(QT_BLOCK(freet)+1)*(QT_BLOCK_SIZE/8));
56	>	(QT_BLOCK(freet)+1)*((QT_BLOCK_SIZE+7)>>3));
57		if (quad_flag == NULL)
58	<	return(EMPTY);
58	>	error(SYSTEM,"qtAlloc(): Unable to allocate memory\n");
59		}
60		treetop += 4;
61		return(freet);
#	Line 60 | Line 64 | qtAlloc()                      /* allocate a quadtree */
64
65		qtClearAllFlags()               /* clear all quadtree branch flags */
66		{
67	<	if (!treetop) return;
68	<	bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*(QT_BLOCK_SIZE/8));
67	>	if (!treetop)
68	>	return;
69	>
70	>	/* Clear the node flags*/
71	>	bzero((char *)quad_flag, (QT_BLOCK(treetop-4)+1)*((QT_BLOCK_SIZE+7)>>3));
72	>	/* Clear set flags */
73	>	qtclearsetflags();
74		}
75
67	–
76		qtFree(qt)                      /* free a quadtree */
77		register QUADTREE  qt;
78		{
#	Line 85 | Line 93 | qtFree(qt)                     /* free a quadtree */
93		qtDone()                        /* free EVERYTHING */
94		{
95		register int    i;
96	<
96	>
97		qtfreeleaves();
98		for (i = 0; i < QT_MAX_BLK; i++)
99		{
#	Line 94 | Line 102 | qtDone()                       /* free EVERYTHING */
102		free((char *)quad_block[i]);
103		quad_block[i] = NULL;
104		}
105	+	/* Free the flags */
106		if (i) free((char *)quad_flag);
107		quad_flag = NULL;
108		quad_free_list = EMPTY;
#	Line 101 | Line 110 | qtDone()                       /* free EVERYTHING */
110		}
111
112		QUADTREE
113	<	*qtLocate_leaf(qtptr,bcoord,t0,t1,t2)
114	<	QUADTREE *qtptr;
115	<	double bcoord[3];
107	<	FVECT t0,t1,t2;
113	>	qtLocate_leaf(qt,bcoord)
114	>	QUADTREE qt;
115	>	BCOORD bcoord[3];
116		{
117		int i;
110	–	QUADTREE *child;
111	–	FVECT a,b,c;
118
119	<	if(QT_IS_TREE(*qtptr))
120	<	{
121	<	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
119	>	if(QT_IS_TREE(qt))
120	>	{
121	>	i = baryi_child(bcoord);
122
123	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
124	<	if(t0)
125	<	{
126	<	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);
123	>	return(qtLocate_leaf(QT_NTH_CHILD(qt,i),bcoord));
124	>	}
125	>	else
126	>	return(qt);
127		}
128
129	<
129	>	/*
130	>	Return the quadtree node containing pt. It is assumed that pt is in
131	>	the root node qt with ws vertices q0,q1,q2 and plane equation peq.
132	>	*/
133		QUADTREE
134	<	*qtRoot_point_locate(qtptr,v0,v1,v2,pt,t0,t1,t2)
135	<	QUADTREE *qtptr;
136	<	FVECT v0,v1,v2;
134	>	qtRoot_point_locate(qt,q0,q1,q2,peq,pt)
135	>	QUADTREE qt;
136	>	FVECT q0,q1,q2;
137	>	FPEQ peq;
138		FVECT pt;
144	–	FVECT t0,t1,t2;
139		{
146	–	int d;
140		int i,x,y;
141	<	QUADTREE *child;
142	<	FVECT n,i_pt,a,b,c;
143	<	double pd,bcoord[3];
141	>	FVECT i_pt;
142	>	double bcoord[3];
143	>	BCOORD bcoordi[3];
144
145	<	/* 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
145	>	/* Will return lowest level triangle containing point: It the
146		point is on an edge or vertex: will return first associated
147		triangle encountered in the child traversal- the others can
148		be derived using triangle adjacency information
149		*/
150	<	if(QT_IS_TREE(*qtptr))
150	>	if(QT_IS_TREE(qt))
151		{
152		/* Find the intersection point */
153	<	tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
175	<	intersect_vector_plane(pt,n,pd,NULL,i_pt);
153	>	intersect_vector_plane(pt,peq,NULL,i_pt);
154
155	<	i = max_index(n);
156	<	x = (i+1)%3;
179	<	y = (i+2)%3;
155	>	x = FP_X(peq);
156	>	y = FP_Y(peq);
157		/* Calculate barycentric coordinates of i_pt */
158	<	bary2d(v0[x],v0[y],v1[x],v1[y],v2[x],v2[y],i_pt[x],i_pt[y],bcoord);
158	>	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],bcoord);
159	>
160	>	/* convert to integer coordinate */
161	>	convert_dtol(bcoord,bcoordi);
162
163	<	i = bary_child(bcoord);
164	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
165	<	#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));
163	>	i = baryi_child(bcoordi);
164	>
165	>	return(qtLocate_leaf(QT_NTH_CHILD(qt,i),bcoordi));
166		}
167		else
168	<	{
208	<	if(t0)
209	<	{
210	<	VCOPY(t0,v0);
211	<	VCOPY(t1,v1);
212	<	VCOPY(t2,v2);
213	<	}
214	<	return(qtptr);
215	<	}
168	>	return(qt);
169		}
170
171
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	–	}
172
173
231	–	QUADTREE
232	–	qtSubdivide_nth_child(qt,n)
233	–	QUADTREE qt;
234	–	int n;
235	–	{
236	–	QUADTREE node;
237	–
238	–	node = qtSubdivide(&(QT_NTH_CHILD(qt,n)));
239	–
240	–	return(node);
241	–	}
242	–
174		/* for triangle v0-v1-v2- returns a,b,c: children are:
175
176		v2                     0: v0,a,c
#	Line 252 | Line 183 | int n;
183		a
184		*/
185
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	–	}
186
187		qtNth_child_tri(v0,v1,v2,a,b,c,i,r0,r1,r2)
188		FVECT v0,v1,v2;
#	Line 268 | Line 191 | int i;
191		FVECT r0,r1,r2;
192		{
193
194	<	switch(i){
195	<	case 0:
196	<	VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
197	<	break;
198	<	case 1:
199	<	VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
200	<	break;
201	<	case 2:
202	<	VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
203	<	break;
204	<	case 3:
205	<	VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
206	<	break;
194	>	if(!a)
195	>	{
196	>	/* Caution: r's must not be equal to v's:will be incorrect */
197	>	switch(i){
198	>	case 0:
199	>	VCOPY(r0,v0);
200	>	EDGE_MIDPOINT_VEC3(r1,v0,v1);
201	>	EDGE_MIDPOINT_VEC3(r2,v2,v0);
202	>	break;
203	>	case 1:
204	>	EDGE_MIDPOINT_VEC3(r0,v0,v1);
205	>	VCOPY(r1,v1);
206	>	EDGE_MIDPOINT_VEC3(r2,v1,v2);
207	>	break;
208	>	case 2:
209	>	EDGE_MIDPOINT_VEC3(r0,v2,v0);
210	>	EDGE_MIDPOINT_VEC3(r1,v1,v2);
211	>	VCOPY(r2,v2);
212	>	break;
213	>	case 3:
214	>	EDGE_MIDPOINT_VEC3(r0,v1,v2);
215	>	EDGE_MIDPOINT_VEC3(r1,v2,v0);
216	>	EDGE_MIDPOINT_VEC3(r2,v0,v1);
217	>	break;
218	>	}
219		}
220	+	else
221	+	{
222	+	switch(i){
223	+	case 0:
224	+	VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
225	+	break;
226	+	case 1:
227	+	VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
228	+	break;
229	+	case 2:
230	+	VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
231	+	break;
232	+	case 3:
233	+	VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
234	+	break;
235	+	}
236	+	}
237		}
238
239		/* Add triangle "id" to all leaf level cells that are children of
#	Line 292 | Line 244 | threshold- the cell is split- and the triangle must be
244		into the new child cells: it is assumed that "v1,v2,v3" are normalized
245		*/
246
247	<	int
248	<	qtRoot_add_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
249	<	QUADTREE *qtptr;
247	>	QUADTREE
248	>	qtRoot_add_tri(qt,q0,q1,q2,t0,t1,t2,id,n)
249	>	QUADTREE qt;
250		FVECT q0,q1,q2;
251		FVECT t0,t1,t2;
252	<	int id;
301	<	int n;
252	>	int id,n;
253		{
254	<	int test;
255	<	int found;
254	>	if(stri_intersect(q0,q1,q2,t0,t1,t2))
255	>	qt = qtAdd_tri(qt,q0,q1,q2,t0,t1,t2,id,n);
256
257	<	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);
257	>	return(qt);
258		}
259
260	<	int
261	<	qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n)
262	<	QUADTREE *qtptr;
260	>	QUADTREE
261	>	qtRoot_remove_tri(qt,q0,q1,q2,t0,t1,t2,id,n)
262	>	QUADTREE qt;
263		FVECT q0,q1,q2;
264		FVECT t0,t1,t2;
265	+	int id,n;
266	+	{
267	+
268	+	if(stri_intersect(q0,q1,q2,t0,t1,t2))
269	+	qt = qtRemove_tri(qt,q0,q1,q2,t0,t1,t2,id,n);
270	+	return(qt);
271	+	}
272	+
273	+
274	+	QUADTREE
275	+	qtAdd_tri(qt,q0,q1,q2,t0,t1,t2,id,n)
276	+	QUADTREE qt;
277	+	FVECT q0,q1,q2;
278	+	FVECT t0,t1,t2;
279		int id;
280		int n;
281		{
323	–	int i,index,test,found;
282		FVECT a,b,c;
283	<	OBJECT os[QT_MAXSET+1],*optr;
326	<	QUADTREE qt;
283	>	OBJECT tset[QT_MAXSET+1],*optr,*tptr;
284		FVECT r0,r1,r2;
285	+	int i;
286
329	–	found = 0;
287		/* if this is tree: recurse */
288	<	if(QT_IS_TREE(*qtptr))
288	>	if(QT_IS_TREE(qt))
289		{
290	+	QT_SET_FLAG(qt);
291		n++;
292		qtSubdivide_tri(q0,q1,q2,a,b,c);
293	<	test = stri_intersect(t0,t1,t2,q0,a,c);
294	<	if(test)
295	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,0),q0,a,c,t0,t1,t2,id,n);
296	<	test = stri_intersect(t0,t1,t2,a,q1,b);
297	<	if(test)
298	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,1),a,q1,b,t0,t1,t2,id,n);
299	<	test = stri_intersect(t0,t1,t2,c,b,q2);
300	<	if(test)
301	<	found |= qtAdd_tri(QT_NTH_CHILD_PTR(*qtptr,2),c,b,q2,t0,t1,t2,id,n);
302	<	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);
293	>
294	>	if(stri_intersect(t0,t1,t2,q0,a,c))
295	>	QT_NTH_CHILD(qt,0) = qtAdd_tri(QT_NTH_CHILD(qt,0),q0,a,c,t0,t1,t2,id,n);
296	>	if(stri_intersect(t0,t1,t2,a,q1,b))
297	>	QT_NTH_CHILD(qt,1) = qtAdd_tri(QT_NTH_CHILD(qt,1),a,q1,b,t0,t1,t2,id,n);
298	>	if(stri_intersect(t0,t1,t2,c,b,q2))
299	>	QT_NTH_CHILD(qt,2) = qtAdd_tri(QT_NTH_CHILD(qt,2),c,b,q2,t0,t1,t2,id,n);
300	>	if(stri_intersect(t0,t1,t2,b,c,a))
301	>	QT_NTH_CHILD(qt,3) = qtAdd_tri(QT_NTH_CHILD(qt,3),b,c,a,t0,t1,t2,id,n);
302	>	return(qt);
303		}
304		else
305		{
306		/* If this leave node emptry- create a new set */
307	<	if(QT_IS_EMPTY(*qtptr))
308	<	*qtptr = qtaddelem(*qtptr,id);
307	>	if(QT_IS_EMPTY(qt))
308	>	qt = qtaddelem(qt,id);
309		else
310		{
311		/* If the set is too large: subdivide */
312	<	optr = qtqueryset(*qtptr);
312	>	optr = qtqueryset(qt);
313
314		if(QT_SET_CNT(optr) < QT_SET_THRESHOLD)
315	<	*qtptr = qtaddelem(*qtptr,id);
316	<	else
317	<	{
318	<	if (n < QT_MAX_LEVELS)
319	<	{
320	<	/* If set size exceeds threshold: subdivide cell and
321	<	reinsert set tris into cell
322	<	*/
323	<	qtgetset(os,*qtptr);
315	>	qt = qtaddelem(qt,id);
316	>	else
317	>	{
318	>	if (n < QT_MAX_LEVELS)
319	>	{
320	>	/* If set size exceeds threshold: subdivide cell and
321	>	reinsert set tris into cell
322	>	*/
323	>	/* CAUTION:If QT_SET_THRESHOLD << QT_MAXSET, and dont add
324	>	more than a few triangles before expanding: then are safe here
325	>	otherwise must check to make sure set size is < MAXSET,
326	>	or  qtgetset could overflow os.
327	>	*/
328	>	tptr = qtqueryset(qt);
329	>	if(QT_SET_CNT(tptr) > QT_MAXSET)
330	>	tptr = (OBJECT *)malloc((QT_SET_CNT(tptr)+1)*sizeof(OBJECT));
331	>	else
332	>	tptr = tset;
333	>	if(!tptr)
334	>	goto memerr;
335
336	<	n++;
337	<	qtfreeleaf(*qtptr);
338	<	qtSubdivide(qtptr);
339	<	found = qtAdd_tri(qtptr,q0,q1,q2,t0,t1,t2,id,n);
336	>	qtgetset(tptr,qt);
337	>	n++;
338	>	qtfreeleaf(qt);
339	>	qtSubdivide(qt);
340	>	qt = qtAdd_tri(qt,q0,q1,q2,t0,t1,t2,id,n);
341
342	<	for(optr = QT_SET_PTR(os),i = QT_SET_CNT(os); i > 0; i--)
343	<	{
344	<	id = QT_SET_NEXT_ELEM(optr);
345	<	qtTri_from_id(id,NULL,NULL,NULL,r0,r1,r2,NULL,NULL,NULL);
346	<	found=qtAdd_tri(qtptr,q0,q1,q2,r0,r1,r2,id,n);
347	<	#ifdef DEBUG
348	<	if(!found)
349	<	eputs("qtAdd_tri():Reinsert-in parent but not children\n");
350	<	#endif
383	<	}
342	>	for(optr = QT_SET_PTR(tptr),i = QT_SET_CNT(tptr); i > 0; i--)
343	>	{
344	>	id = QT_SET_NEXT_ELEM(optr);
345	>	if(!qtTri_from_id(id,r0,r1,r2))
346	>	continue;
347	>	qt  = qtAdd_tri(qt,q0,q1,q2,r0,r1,r2,id,n);
348	>	}
349	>	if(tptr != tset)
350	>	free(tptr);
351		}
352		else
353	<	if(QT_SET_CNT(optr) < QT_MAXSET)
354	<	*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	<	}
353	>	qt = qtaddelem(qt,id);
354	>	}
355		}
356		}
357	<	return(TRUE);
357	>	return(qt);
358	>	memerr:
359	>	error(SYSTEM,"qtAdd_tri():Unable to allocate memory");
360		}
361
362
363	<	int
364	<	qtApply_to_tri_cells(qtptr,t0,t1,t2,v0,v1,v2,func,arg)
365	<	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;
363	>	QUADTREE
364	>	qtRemove_tri(qt,id,q0,q1,q2,t0,t1,t2)
365	>	QUADTREE qt;
366		int id;
367	+	FVECT q0,q1,q2;
368		FVECT t0,t1,t2;
439	–	FVECT v0,v1,v2;
369		{
441	–
442	–	int test;
443	–	int i;
370		FVECT a,b,c;
445	–	OBJECT os[QT_MAXSET+1];
371
372		/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
373	<	test = stri_intersect(t0,t1,t2,v0,v1,v2);
373	>	if(!stri_intersect(t0,t1,t2,q0,q1,q2))
374	>	return(qt);
375
450	–	/* If triangles do not overlap: done */
451	–	if(!test)
452	–	return(FALSE);
453	–
376		/* if this is tree: recurse */
377	<	if(QT_IS_TREE(*qtptr))
377	>	if(QT_IS_TREE(qt))
378		{
379	<	qtSubdivide_tri(v0,v1,v2,a,b,c);
380	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,0),id,t0,t1,t2,v0,a,c);
381	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,1),id,t0,t1,t2,a,v1,b);
382	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,2),id,t0,t1,t2,c,b,v2);
383	<	qtRemove_tri(QT_NTH_CHILD_PTR(*qtptr,3),id,t0,t1,t2,b,c,a);
379	>	qtSubdivide_tri(q0,q1,q2,a,b,c);
380	>	QT_NTH_CHILD(qt,0) = qtRemove_tri(QT_NTH_CHILD(qt,0),id,t0,t1,t2,q0,a,c);
381	>	QT_NTH_CHILD(qt,1) = qtRemove_tri(QT_NTH_CHILD(qt,1),id,t0,t1,t2,a,q1,b);
382	>	QT_NTH_CHILD(qt,2) = qtRemove_tri(QT_NTH_CHILD(qt,2),id,t0,t1,t2,c,b,q2);
383	>	QT_NTH_CHILD(qt,3) = qtRemove_tri(QT_NTH_CHILD(qt,3),id,t0,t1,t2,b,c,a);
384	>	return(qt);
385		}
386		else
387		{
388	<	if(QT_IS_EMPTY(*qtptr))
388	>	if(QT_IS_EMPTY(qt))
389		{
390		#ifdef DEBUG
391		eputs("qtRemove_tri(): triangle not found\n");
#	Line 471 | Line 394 | FVECT v0,v1,v2;
394		/* remove id from set */
395		else
396		{
397	<	if(!qtinset(*qtptr,id))
397	>	if(!qtinset(qt,id))
398		{
399		#ifdef DEBUG
400		eputs("qtRemove_tri(): tri not in set\n");
401		#endif
402		}
403		else
404	<	{
482	<	*qtptr = qtdelelem(*qtptr,id);
483	<	}
404	>	qt = qtdelelem(qt,id);
405		}
406		}
407	<	return(TRUE);
407	>	return(qt);
408		}
409
410
411	<	int
412	<	move_to_nbr(b,db0,db1,db2,tptr)
413	<	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;
411	>	QUADTREE
412	>	qtVisit_tri_interior(qt,q0,q1,q2,t0,t1,t2,n0,n1,n2,n,func,f,argptr)
413	>	QUADTREE qt;
414		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;
415		FVECT t0,t1,t2;
416	+	FVECT n0,n1,n2;
417		int n;
418	<	int (*func)();
419	<	int *arg1,arg2;
418	>	int (*func)(),*f;
419	>	int *argptr;
420		{
421	<	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;
421	>	FVECT a,b,c;
422
423		/* If qt Flag set, or qt vertices interior to t0t1t2-descend */
424		tree_modified:
425
426	<	if(QT_IS_TREE(*qtptr))
426	>	if(QT_IS_TREE(qt))
427		{
428	<	if(QT_IS_FLAG(*qtptr) ||  point_in_stri(t0,t1,t2,q0))
428	>	if(QT_IS_FLAG(qt) ||  point_in_stri_n(n0,n1,n2,q0))
429		{
430	<	QT_SET_FLAG(*qtptr);
430	>	QT_SET_FLAG(qt);
431		qtSubdivide_tri(q0,q1,q2,a,b,c);
432		/* descend to children */
433	<	for(i=0;i < 4; i++)
434	<	{
435	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
436	<	qtNth_child_tri(q0,q1,q2,a,b,c,i,c0,c1,c2);
437	<	qtVisit_tri_interior(child,c0,c1,c2,t0,t1,t2,n+1,
438	<	func,arg1,arg2);
439	<	}
433	>	QT_NTH_CHILD(qt,0) = qtVisit_tri_interior(QT_NTH_CHILD(qt,0),
434	>	q0,a,c,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
435	>	QT_NTH_CHILD(qt,1) = qtVisit_tri_interior(QT_NTH_CHILD(qt,1),
436	>	a,q1,b,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
437	>	QT_NTH_CHILD(qt,2) = qtVisit_tri_interior(QT_NTH_CHILD(qt,2),
438	>	c,b,q2,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
439	>	QT_NTH_CHILD(qt,3) = qtVisit_tri_interior(QT_NTH_CHILD(qt,3),
440	>	b,c,a,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
441		}
442		}
443		else
444	<	{
445	<	/* NOTE THIS IN TRI TEST Could be replaced by a flag */
873	<	if(!QT_IS_EMPTY(*qtptr))
444	>	if((!QT_IS_EMPTY(qt) && QT_LEAF_IS_FLAG(qt)) ||
445	>	point_in_stri_n(n0,n1,n2,q0))
446		{
447	<	if(qtinset(*qtptr,arg2))
448	<	if(func(qtptr,q0,q1,q2,t0,t1,t2,n,arg1,arg2)==QT_MODIFIED)
449	<	goto tree_modified;
450	<	else
451	<	return;
447	>	func(&qt,f,argptr,q0,q1,q2,t0,t1,t2,n);
448	>	if(QT_FLAG_IS_MODIFIED(*f))
449	>	{
450	>	QT_SET_FLAG(qt);
451	>	goto tree_modified;
452	>	}
453	>	if(QT_IS_LEAF(qt))
454	>	QT_LEAF_SET_FLAG(qt);
455	>	else
456	>	if(QT_IS_TREE(qt))
457	>	QT_SET_FLAG(qt);
458		}
459	<	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	<	}
459	>	return(qt);
460		}
461
462
463
889	–
890	–
891	–
464		int
465	<	qtVisit_tri_edges(qtptr,b,db,wptr,sfactor,func,arg1,arg2)
466	<	QUADTREE *qtptr;
467	<	double b[3],db[3][3];
468	<	int *wptr;
897	<	double sfactor;
898	<	int (*func)();
899	<	int *arg1,arg2;
465	>	move_to_nbri(b,db0,db1,db2,tptr)
466	>	BCOORD b[3];
467	>	BDIR db0,db1,db2;
468	>	TINT *tptr;
469		{
470	<	int i,found;
471	<	QUADTREE *child;
472	<	int nbr,next,w;
473	<	double t;
474	<	#ifdef DEBUG_TEST_DRIVER
475	<	FVECT a1,b1,c1;
476	<	int Pick_parent = Pick_cnt-1;
477	<	qtSubdivide_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
478	<	Pick_v2[Pick_parent],a1,b1,c1);
479	<	#endif
470	>	TINT t,dt;
471	>	BCOORD bc;
472	>	int nbr;
473	>
474	>	nbr = -1;
475	>	*tptr = 0;
476	>	/* Advance to next node */
477	>	if(b[0]==0 && db0 < 0)
478	>	return(0);
479	>	if(b[1]==0 && db1 < 0)
480	>	return(1);
481	>	if(b[2]==0 && db2 < 0)
482	>	return(2);
483
484	<	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
484	>	if(db0 < 0)
485		{
486	<	#ifdef DEBUG_TEST_DRIVER
487	<	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
488	<	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	<	}
486	>	bc = b[0]<<SHIFT_MAXBCOORD;
487	>	t = bc/-db0;
488	>	nbr = 0;
489		}
490	<
491	<	}
492	<
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++)
490	>	else
491	>	t = HUGET;
492	>	if(db1 < 0)
493		{
494	<	/* If processing 3rd edge-dont need info for t1 */
495	<	if(i==1 && w==2)
496	<	continue;
497	<	/* project the dir as well */
498	<	intersect_vector_plane(tri[i],n,pd,NULL,i_pt);
499	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[i]);
500	<	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]);
494	>	bc = b[1] <<SHIFT_MAXBCOORD;
495	>	dt = bc/-db1;
496	>	if( dt < t)
497	>	{
498	>	t = dt;
499	>	nbr = 1;
500	>	}
501		}
502	<	#ifdef DEBUG_TEST_DRIVER
503	<	VCOPY(Pick_v0[Pick_cnt],q0);
504	<	VCOPY(Pick_v1[Pick_cnt],q1);
505	<	VCOPY(Pick_v2[Pick_cnt],q2);
506	<	Pick_cnt++;
507	<	#endif
508	<	/* trace the ray starting with this node */
509	<	nbr = qtVisit_tri_edges(qtptr,b[w],db,wptr,1.0,func,arg1,arg2);
510	<	return(nbr);
511	<
502	>	if(db2 < 0)
503	>	{
504	>	bc = b[2] << SHIFT_MAXBCOORD;
505	>	dt = bc/-db2;
506	>	if( dt < t)
507	>	{
508	>	t = dt;
509	>	nbr = 2;
510	>	}
511	>	}
512	>	*tptr = t;
513	>	return(nbr);
514		}
515
516	<
517	<
518	<
519	<	/* NOTE: SINCE DIR could be unit: then we could use integer math */
520	<	int
521	<	qtVisit_tri_edges2(qtptr,b,db0,db1,db2,
522	<	db,wptr,t,sign,sfactor,func,arg1,arg2)
523	<	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;
516	>	QUADTREE
517	>	qtVisit_tri_edges(qt,b,db0,db1,db2,db,wptr,nextptr,t,sign,sfactor,func,f,argptr)
518	>	QUADTREE qt;
519	>	BCOORD b[3];
520	>	BDIR db0,db1,db2,db[3][3];
521	>	int *wptr,*nextptr;
522	>	TINT t[3];
523	>	int sign,sfactor;
524		int (*func)();
525	<	int *arg1,arg2;
525	>	int *f,*argptr;
526		{
527		int i,found;
528	<	QUADTREE *child;
528	>	QUADTREE child;
529		int nbr,next,w;
530	<	double t_l,t_g;
531	<	#ifdef DEBUG_TEST_DRIVER
532	<	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))
530	>	TINT t_g,t_l,t_i,l;
531	>
532	>	if(QT_IS_TREE(qt))
533		{
534	<	/* Find the appropriate child and reset the coord */
535	<	i = bary_child(b);
534	>	/* Find the appropriate child and reset the coord */
535	>	i = baryi_child(b);
536
537	<	QT_SET_FLAG(*qtptr);
537	>	QT_SET_FLAG(qt);
538
539	<	for(;;)
540	<	{
541	<	w = *wptr;
542	<	child = QT_NTH_CHILD_PTR(*qtptr,i);
539	>	for(;;)
540	>	{
541	>	w = *wptr;
542	>	child = QT_NTH_CHILD(qt,i);
543	>	if(i != 3)
544	>	QT_NTH_CHILD(qt,i) =
545	>	qtVisit_tri_edges(child,b,db0,db1,db2,db,wptr,nextptr,t,sign,
546	>	sfactor+1,func,f,argptr);
547	>	else
548	>	/* If the center cell- must flip direction signs */
549	>	QT_NTH_CHILD(qt,i) =
550	>	qtVisit_tri_edges(child,b,-db0,-db1,-db2,db,wptr,nextptr,t,1-sign,
551	>	sfactor+1,func,f,argptr);
552
553	<	if(i != 3)
554	<	nbr = qtVisit_tri_edges2(child,b,db0,db1,db2,
555	<	db,wptr,t,sign,
556	<	sfactor*2.0,func,arg1,arg2);
557	<	else
558	<	/* If the center cell- must flip direction signs */
559	<	nbr = qtVisit_tri_edges2(child,b,-db0,-db1,-db2,
560	<	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;
553	>	if(QT_FLAG_IS_DONE(*f))
554	>	return(qt);
555	>	if(*wptr != w)
556	>	{
557	>	w = *wptr;
558	>	db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
559	>	if(sign)
560	>	{  db0 *= -1;db1 *= -1; db2 *= -1;}
561		}
562	+	/* If in same block: traverse */
563	+	if(i==3)
564	+	next = *nextptr;
565	+	else
566	+	if(*nextptr == i)
567	+	next = 3;
568	+	else
569	+	{
570	+	/* reset the barycentric coordinates in the parents*/
571	+	baryi_parent(b,i);
572	+	/* Else pop up to parent and traverse from there */
573	+	return(qt);
574	+	}
575	+	baryi_from_child(b,i,next);
576	+	i = next;
577	+	}
578		}
579		else
580		{
581	<	#ifdef DEBUG_TEST_DRIVER
582	<	qtNth_child_tri(Pick_v0[Pick_parent],Pick_v1[Pick_parent],
583	<	Pick_v2[Pick_parent],a1,b1,c1,i,Pick_v0[Pick_cnt],
584	<	Pick_v1[Pick_cnt],Pick_v2[Pick_cnt]);
585	<	Pick_cnt++;
586	<	#endif
587	<
588	<	if(func(qtptr,arg1,arg2) == QT_DONE)
589	<	return(QT_DONE);
590	<
591	<	/* Advance to next node */
592	<	w = *wptr;
593	<	while(1)
581	>	#ifdef TEST_DRIVER
582	>	if(Pick_cnt < 500)
583	>	Pick_q[Pick_cnt++]=qt;
584	>	#endif;
585	>	func(&qt,f,argptr);
586	>	if(QT_FLAG_IS_DONE(*f))
587	>	{
588	>	if(!QT_IS_EMPTY(qt))
589	>	QT_LEAF_SET_FLAG(qt);
590	>	return(qt);
591	>	}
592	>
593	>	if(!QT_IS_EMPTY(qt))
594	>	QT_LEAF_SET_FLAG(qt);
595	>	/* Advance to next node */
596	>	w = *wptr;
597	>	while(1)
598		{
599	<	nbr = move_to_nbr(b,db0,db1,db2,&t_l);
600	<
601	<	t_g = t_l/sfactor;
602	<	#ifdef DEBUG
1162	<	if(t[w] <= 0.0)
1163	<	eputs("qtVisit_tri_edges2():negative t\n");
1164	<	#endif
599	>	nbr = move_to_nbri(b,db0,db1,db2,&t_i);
600	>
601	>	t_g = t_i >> sfactor;
602	>
603		if(t_g >= t[w])
604		{
605		if(w == 2)
606	<	return(QT_DONE);
607	<
608	<	b[0] += (t[w])*sfactor*db0;
609	<	b[1] += (t[w])*sfactor*db1;
610	<	b[2] += (t[w])*sfactor*db2;
606	>	{
607	>	QT_FLAG_SET_DONE(*f);
608	>	return(qt);
609	>	}
610	>	/* The edge fits in the cell- therefore the result of shifting
611	>	db by sfactor is guaranteed to be less than MAXBCOORD
612	>	*/
613	>	/* Caution: (t[w]*db) must occur before divide by MAXBCOORD
614	>	since t[w] will always be < MAXBCOORD
615	>	*/
616	>	l = t[w] << sfactor;
617	>	/* NOTE: Change divides to Shift and multiply by sign*/
618	>	b[0] += (l*db0)/MAXBCOORD;
619	>	b[1] += (l*db1)/MAXBCOORD;
620	>	b[2] += (l*db2)/MAXBCOORD;
621		w++;
622	<	db0 = db[w][0];
1175	<	db1 = db[w][1];
1176	<	db2 = db[w][2];
622	>	db0 = db[w][0]; db1 = db[w][1]; db2 = db[w][2];
623		if(sign)
624	<	{  db0 *= -1.0;db1 *= -1.0; db2 *= -1.0;}
624	>	{  db0 *= -1;db1 *= -1; db2 *= -1;}
625		}
626	<	else
627	<	if(nbr != INVALID)
628	<	{
629	<	b[0] += t_l * db0;
630	<	b[1] += t_l * db1;
631	<	b[2] += t_l * db2;
632	<
626	>	else
627	>	{
628	>	/* Caution: (t_i*db) must occur before divide by MAXBCOORD
629	>	since t_i will always be < MAXBCOORD*/
630	>	/* NOTE: Change divides to Shift and  by sign*/
631	>	b[0] += (t_i *db0) / MAXBCOORD;
632	>	b[1] += (t_i *db1) / MAXBCOORD;
633	>	b[2] += (t_i *db2) / MAXBCOORD;
634	>
635		t[w] -= t_g;
636		*wptr = w;
637	<	return(nbr);
637	>	*nextptr = nbr;
638	>	return(qt);
639		}
1191	–	else
1192	–	return(INVALID);
640		}
641	<	}
1195	<
641	>	}
642		}
643
644
645	<	int
646	<	qtRoot_visit_tri_edges2(qtptr,q0,q1,q2,tri,i_pt,wptr,func,arg1,arg2)
647	<	QUADTREE *qtptr;
645	>	QUADTREE
646	>	qtRoot_visit_tri_edges(qt,q0,q1,q2,peq,tri,i_pt,wptr,nextptr,func,f,argptr)
647	>	QUADTREE qt;
648		FVECT q0,q1,q2;
649	+	FPEQ  peq;
650		FVECT tri[3],i_pt;
651	<	int *wptr;
651	>	int *wptr,*nextptr;
652		int (*func)();
653	<	int *arg1,arg2;
653	>	int *f,*argptr;
654		{
655	<	int x,y,z,nbr,w,i,j;
656	<	QUADTREE *child;
657	<	FVECT n,c,d,v[3];
658	<	double pd,b[4][3],db[3][3],et[3],t[3],exit_pt;
659	<
655	>	int x,y,z,w,i,j,first;
656	>	QUADTREE child;
657	>	FVECT c,d,v[3];
658	>	double b[4][3],db[3][3],et[3],exit_pt;
659	>	BCOORD bi[3];
660	>	TINT t[3];
661	>	BDIR dbi[3][3];
662	>
663	>	first =0;
664		w = *wptr;
665	<
665	>	if(w==-1)
666	>	{
667	>	first = 1;
668	>	*wptr = w = 0;
669	>	}
670		/* Project the origin onto the root node plane */
671
672	+	t[0] = t[1] = t[2] = 0;
673		/* Find the intersection point of the origin */
1218	–	tri_plane_equation(q0,q1,q2,n,&pd,FALSE);
674		/* map to 2d by dropping maximum magnitude component of normal */
675	<	z = max_index(n);
676	<	x = (z+1)%3;
677	<	y = (z+2)%3;
675	>
676	>	x = FP_X(peq);
677	>	y = FP_Y(peq);
678	>	z = FP_Z(peq);
679		/* Calculate barycentric coordinates for current vertex */
680	<	if(w != -1)
1225	<	{
680	>	if(!first)
681		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	–	}
682		else
683		/* Just starting: b[0] is the origin point: guaranteed to be valid b*/
684		{
685	<	w = 0;
1234	<	intersect_vector_plane(tri[0],n,pd,&(et[0]),v[0]);
685	>	intersect_vector_plane(tri[0],peq,&(et[0]),v[0]);
686		bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);
687		VCOPY(b[3],b[0]);
688		}
689
1239	–
690		j = (w+1)%3;
691	<	intersect_vector_plane(tri[j],n,pd,&(et[j]),v[j]);
691	>	intersect_vector_plane(tri[j],peq,&(et[j]),v[j]);
692		bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);
693		if(et[j] < 0.0)
694		{
695		VSUB(db[w],b[3],b[j]);
696	<	t[w] = FHUGE;
696	>	t[w] = HUGET;
697		}
698		else
699		{
700	+	/* NOTE: for stability: do not increment with ipt- use full dir and
701	+	calculate t: but for wrap around case: could get same problem?
702	+	*/
703		VSUB(db[w],b[j],b[3]);
704	<	t[w] = 1.0;
705	<	move_to_nbr(b[3],db[w][0],db[w][1],db[w][2],&exit_pt);
706	<	if(exit_pt >= 1.0)
704	>	/* Check if the point is out side of the triangle: if so t[w] =HUGET */
705	>	if((fabs(b[j][0])>(FTINY+1.0)) ||(fabs(b[j][1])>(FTINY+1.0)) ||
706	>	(fabs(b[j][2])>(FTINY+1.0))||(b[j][0] <-FTINY) ||
707	>	(b[j][1]<-FTINY) ||(b[j][2]<-FTINY))
708	>	t[w] = HUGET;
709	>	else
710		{
711	<	for(;j < 3;j++)
712	<	{
713	<	i = (j+1)%3;
714	<	if(i!= w)
715	<	{
716	<	intersect_vector_plane(tri[i],n,pd,&(et[i]),v[i]);
717	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
718	<	v[i][y],b[i]);
719	<	}
720	<	if(et[i] < 0.0)
721	<	{
722	<	VSUB(db[j],b[j],b[i]);
723	<	t[j] = FHUGE;
724	<	break;
725	<	}
726	<	else
727	<	{
728	<	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)
711	>	/* The next point is in the triangle- so db values will be in valid
712	>	range and t[w]= MAXT
713	>	*/
714	>	t[w] = MAXT;
715	>	if(j != 0)
716	>	for(;j < 3;j++)
717	>	{
718	>	i = (j+1)%3;
719	>	if(!first || i != 0)
720	>	{
721	>	intersect_vector_plane(tri[i],peq,&(et[i]),v[i]);
722	>	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
723	>	v[i][y],b[i]);
724	>	}
725	>	if(et[i] < 0.0)
726	>	{
727	>	VSUB(db[j],b[j],b[i]);
728	>	t[j] = HUGET;
729		break;
730	<	}
730	>	}
731	>	else
732	>	{
733	>	VSUB(db[j],b[i],b[j]);
734	>	if((fabs(b[j][0])>(FTINY+1.0))||(fabs(b[j][1])>(FTINY+1.0)) ||
735	>	(fabs(b[j][2])>(FTINY+1.0))||(b[i][0] <-FTINY) ||
736	>	(b[i][1]<-FTINY) ||(b[i][2]<-FTINY))
737	>	{
738	>	t[j] = HUGET;
739	>	break;
740	>	}
741	>	else
742	>	t[j] = MAXT;
743	>	}
744	>	}
745		}
746		}
747	<	*wptr = w;
747	>	bary_dtol(b[3],db,bi,dbi,t,w);
748	>
749		/* trace the ray starting with this node */
750	<	nbr = qtVisit_tri_edges2(qtptr,b[3],db[w][0],db[w][1],db[w][2],
751	<	db,wptr,t,0,1.0,func,arg1,arg2);
752	<	if(nbr != INVALID && nbr != QT_DONE)
753	<	{
754	<	i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
755	<	i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
756	<	i_pt[z] = (-n[x]*i_pt[x] - n[y]*i_pt[y] -pd)/n[z];
757	<	}
758	<	return(nbr);
750	>	qt = qtVisit_tri_edges(qt,bi,dbi[w][0],dbi[w][1],dbi[w][2],
751	>	dbi,wptr,nextptr,t,0,0,func,f,argptr);
752	>	if(!QT_FLAG_IS_DONE(*f))
753	>	{
754	>	b[3][0] = (double)bi[0]/(double)MAXBCOORD;
755	>	b[3][1] = (double)bi[1]/(double)MAXBCOORD;
756	>	b[3][2] = (double)bi[2]/(double)MAXBCOORD;
757	>	i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
758	>	i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
759	>	i_pt[z] = (-FP_N(peq)[x]*i_pt[x] - FP_N(peq)[y]*i_pt[y]-FP_D(peq))/FP_N(peq)[z];
760	>	}
761	>	return(qt);
762
763		}
764
765
766	+	QUADTREE
767	+	qtRoot_trace_ray(qt,q0,q1,q2,peq,orig,dir,nextptr,func,f,argptr)
768	+	QUADTREE qt;
769	+	FVECT q0,q1,q2;
770	+	FPEQ  peq;
771	+	FVECT orig,dir;
772	+	int *nextptr;
773	+	int (*func)();
774	+	int *f,*argptr;
775	+	{
776	+	int x,y,z,nbr,w,i;
777	+	QUADTREE child;
778	+	FVECT v,i_pt;
779	+	double b[2][3],db[3],et[2],d,br[3];
780	+	BCOORD bi[3];
781	+	TINT t[3];
782	+	BDIR dbi[3][3];
783	+
784	+	/* Project the origin onto the root node plane */
785	+	t[0] = t[1] = t[2] = 0;
786
787	+	VADD(v,orig,dir);
788	+	/* Find the intersection point of the origin */
789	+	/* map to 2d by dropping maximum magnitude component of normal */
790	+	x = FP_X(peq);
791	+	y = FP_Y(peq);
792	+	z = FP_Z(peq);
793
794	+	/* Calculate barycentric coordinates for current vertex */
795	+	intersect_vector_plane(orig,peq,&(et[0]),i_pt);
796	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[0]);
797
798	+	intersect_vector_plane(v,peq,&(et[1]),i_pt);
799	+	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[1]);
800	+	if(et[1] < 0.0)
801	+	VSUB(db,b[0],b[1]);
802	+	else
803	+	VSUB(db,b[1],b[0]);
804	+	t[0] = HUGET;
805	+	convert_dtol(b[0],bi);
806	+	if(et[1]<0.0 ||(fabs(b[1][0])>(FTINY+1.0))||(fabs(b[1][1])>(FTINY+1.0)) ||
807	+	(fabs(b[1][2])>(FTINY+1.0))||(b[1][0] <-FTINY) ||
808	+	(b[1][1]<-FTINY) ||(b[1][2]<-FTINY))
809	+	{
810	+	max_index(db,&d);
811	+	for(i=0; i< 3; i++)
812	+	dbi[0][i] = (BDIR)(db[i]/d*MAXBDIR);
813	+	}
814	+	else
815	+	for(i=0; i< 3; i++)
816	+	dbi[0][i] = (BDIR)(db[i]*MAXBDIR);
817	+	w=0;
818	+	/* trace the ray starting with this node */
819	+	qt = qtVisit_tri_edges(qt,bi,dbi[0][0],dbi[0][1],dbi[0][2], dbi,&w,
820	+	nextptr,t,0,0,func,f,argptr);
821	+	if(!QT_FLAG_IS_DONE(*f))
822	+	{
823	+	br[0] = (double)bi[0]/(double)MAXBCOORD;
824	+	br[1] = (double)bi[1]/(double)MAXBCOORD;
825	+	br[2] = (double)bi[2]/(double)MAXBCOORD;
826	+	orig[x] = br[0]*q0[x] + br[1]*q1[x] + br[2]*q2[x];
827	+	orig[y] = br[0]*q0[y] + br[1]*q1[y] + br[2]*q2[y];
828	+	orig[z]=(-FP_N(peq)[x]*orig[x] -
829	+	FP_N(peq)[y]*orig[y]-FP_D(peq))/FP_N(peq)[z];
830	+	}
831	+	return(qt);
832	+
833	+	}
834
835
836

Diff Legend

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