ViewVC Help

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

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

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

Diff Legend

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