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.17 by schorsch, Mon Jun 30 14:59:12 2003 UTC

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

Diff Legend

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