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.8 by gwlarson, Wed Nov 11 12:05:39 1998 UTC vs.
Revision 3.13 by gwlarson, Thu Jun 10 15:22:23 1999 UTC

#	Line 30 | Line 30 | extern FVECT Pick_point[500];
30		extern int Pick_q[500];
31
32		#endif
33	–	int Incnt=0;
33
34	+	qtremovelast(qt,id)
35	+	QUADTREE qt;
36	+	int id;
37	+	{
38	+	if(qtqueryset(qt)[(*(qtqueryset(qt)))] != id)
39	+	error(CONSISTENCY,"not removed\n");
40	+	((*(qtqueryset(qt)))--);
41	+	}
42		QUADTREE
43		qtAlloc()                       /* allocate a quadtree */
44		{
#	Line 110 | Line 117 | qtDone()                       /* free EVERYTHING */
117		}
118
119		QUADTREE
120	<	qtLocate_leaf(qt,bcoord)
120	>	qtLocate(qt,bcoord)
121		QUADTREE qt;
122		BCOORD bcoord[3];
123		{
#	Line 118 | Line 125 | BCOORD bcoord[3];
125
126		if(QT_IS_TREE(qt))
127		{
128	<	i = baryi_child(bcoord);
128	>	i = bary_child(bcoord);
129
130	<	return(qtLocate_leaf(QT_NTH_CHILD(qt,i),bcoord));
130	>	return(qtLocate(QT_NTH_CHILD(qt,i),bcoord));
131		}
132		else
133		return(qt);
134		}
135
129	–	/*
130	–	Return the quadtree node containing pt. It is assumed that pt is in
131	–	the root node qt with ws vertices q0,q1,q2 and plane equation peq.
132	–	*/
133	–	QUADTREE
134	–	qtRoot_point_locate(qt,q0,q1,q2,peq,pt)
135	–	QUADTREE qt;
136	–	FVECT q0,q1,q2;
137	–	FPEQ peq;
138	–	FVECT pt;
139	–	{
140	–	int i,x,y;
141	–	FVECT i_pt;
142	–	double bcoord[3];
143	–	BCOORD bcoordi[3];
144	–
145	–	/* Will return lowest level triangle containing point: It the
146	–	point is on an edge or vertex: will return first associated
147	–	triangle encountered in the child traversal- the others can
148	–	be derived using triangle adjacency information
149	–	*/
150	–	if(QT_IS_TREE(qt))
151	–	{
152	–	/* Find the intersection point */
153	–	intersect_vector_plane(pt,peq,NULL,i_pt);
154	–
155	–	x = FP_X(peq);
156	–	y = FP_Y(peq);
157	–	/* Calculate barycentric coordinates of i_pt */
158	–	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],bcoord);
159	–
160	–	/* convert to integer coordinate */
161	–	convert_dtol(bcoord,bcoordi);
162	–
163	–	i = baryi_child(bcoordi);
164	–
165	–	return(qtLocate_leaf(QT_NTH_CHILD(qt,i),bcoordi));
166	–	}
167	–	else
168	–	return(qt);
169	–	}
170	–
171	–
172	–
173	–
174	–	/* for triangle v0-v1-v2- returns a,b,c: children are:
175	–
176	–	v2                     0: v0,a,c
177	–	/\                     1: a,v1,b
178	–	/2 \                    2: c,b,v2
179	–	c/____\b                  3: b,c,a
180	–	/\    /\
181	–	/0 \3 /1 \
182	–	v0____\/____\v1
183	–	a
184	–	*/
185	–
186	–
187	–	qtNth_child_tri(v0,v1,v2,a,b,c,i,r0,r1,r2)
188	–	FVECT v0,v1,v2;
189	–	FVECT a,b,c;
190	–	int i;
191	–	FVECT r0,r1,r2;
192	–	{
193	–
194	–	if(!a)
195	–	{
196	–	/* Caution: r's must not be equal to v's:will be incorrect */
197	–	switch(i){
198	–	case 0:
199	–	VCOPY(r0,v0);
200	–	EDGE_MIDPOINT_VEC3(r1,v0,v1);
201	–	EDGE_MIDPOINT_VEC3(r2,v2,v0);
202	–	break;
203	–	case 1:
204	–	EDGE_MIDPOINT_VEC3(r0,v0,v1);
205	–	VCOPY(r1,v1);
206	–	EDGE_MIDPOINT_VEC3(r2,v1,v2);
207	–	break;
208	–	case 2:
209	–	EDGE_MIDPOINT_VEC3(r0,v2,v0);
210	–	EDGE_MIDPOINT_VEC3(r1,v1,v2);
211	–	VCOPY(r2,v2);
212	–	break;
213	–	case 3:
214	–	EDGE_MIDPOINT_VEC3(r0,v1,v2);
215	–	EDGE_MIDPOINT_VEC3(r1,v2,v0);
216	–	EDGE_MIDPOINT_VEC3(r2,v0,v1);
217	–	break;
218	–	}
219	–	}
220	–	else
221	–	{
222	–	switch(i){
223	–	case 0:
224	–	VCOPY(r0,v0); VCOPY(r1,a);    VCOPY(r2,c);
225	–	break;
226	–	case 1:
227	–	VCOPY(r0,a); VCOPY(r1,v1);    VCOPY(r2,b);
228	–	break;
229	–	case 2:
230	–	VCOPY(r0,c); VCOPY(r1,b);    VCOPY(r2,v2);
231	–	break;
232	–	case 3:
233	–	VCOPY(r0,b); VCOPY(r1,c);    VCOPY(r2,a);
234	–	break;
235	–	}
236	–	}
237	–	}
238	–
239	–	/* Add triangle "id" to all leaf level cells that are children of
240	–	quadtree pointed to by "qtptr" with cell vertices "t1,t2,t3"
241	–	that it overlaps (vertex and edge adjacencies do not count
242	–	as overlapping). If the addition of the triangle causes the cell to go over
243	–	threshold- the cell is split- and the triangle must be recursively inserted
244	–	into the new child cells: it is assumed that "v1,v2,v3" are normalized
245	–	*/
246	–
247	–	QUADTREE
248	–	qtRoot_add_tri(qt,q0,q1,q2,t0,t1,t2,id,n)
249	–	QUADTREE qt;
250	–	FVECT q0,q1,q2;
251	–	FVECT t0,t1,t2;
252	–	int id,n;
253	–	{
254	–	if(stri_intersect(q0,q1,q2,t0,t1,t2))
255	–	qt = qtAdd_tri(qt,q0,q1,q2,t0,t1,t2,id,n);
256	–
257	–	return(qt);
258	–	}
259	–
260	–	QUADTREE
261	–	qtRoot_remove_tri(qt,q0,q1,q2,t0,t1,t2,id,n)
262	–	QUADTREE qt;
263	–	FVECT q0,q1,q2;
264	–	FVECT t0,t1,t2;
265	–	int id,n;
266	–	{
267	–
268	–	if(stri_intersect(q0,q1,q2,t0,t1,t2))
269	–	qt = qtRemove_tri(qt,q0,q1,q2,t0,t1,t2,id,n);
270	–	return(qt);
271	–	}
272	–
273	–
274	–	QUADTREE
275	–	qtAdd_tri(qt,q0,q1,q2,t0,t1,t2,id,n)
276	–	QUADTREE qt;
277	–	FVECT q0,q1,q2;
278	–	FVECT t0,t1,t2;
279	–	int id;
280	–	int n;
281	–	{
282	–	FVECT a,b,c;
283	–	OBJECT tset[QT_MAXSET+1],*optr,*tptr;
284	–	FVECT r0,r1,r2;
285	–	int i;
286	–
287	–	/* if this is tree: recurse */
288	–	if(QT_IS_TREE(qt))
289	–	{
290	–	QT_SET_FLAG(qt);
291	–	n++;
292	–	qtSubdivide_tri(q0,q1,q2,a,b,c);
293	–
294	–	if(stri_intersect(t0,t1,t2,q0,a,c))
295	–	QT_NTH_CHILD(qt,0) = qtAdd_tri(QT_NTH_CHILD(qt,0),q0,a,c,t0,t1,t2,id,n);
296	–	if(stri_intersect(t0,t1,t2,a,q1,b))
297	–	QT_NTH_CHILD(qt,1) = qtAdd_tri(QT_NTH_CHILD(qt,1),a,q1,b,t0,t1,t2,id,n);
298	–	if(stri_intersect(t0,t1,t2,c,b,q2))
299	–	QT_NTH_CHILD(qt,2) = qtAdd_tri(QT_NTH_CHILD(qt,2),c,b,q2,t0,t1,t2,id,n);
300	–	if(stri_intersect(t0,t1,t2,b,c,a))
301	–	QT_NTH_CHILD(qt,3) = qtAdd_tri(QT_NTH_CHILD(qt,3),b,c,a,t0,t1,t2,id,n);
302	–	return(qt);
303	–	}
304	–	else
305	–	{
306	–	/* If this leave node emptry- create a new set */
307	–	if(QT_IS_EMPTY(qt))
308	–	qt = qtaddelem(qt,id);
309	–	else
310	–	{
311	–	/* If the set is too large: subdivide */
312	–	optr = qtqueryset(qt);
313	–
314	–	if(QT_SET_CNT(optr) < QT_SET_THRESHOLD)
315	–	qt = qtaddelem(qt,id);
316	–	else
317	–	{
318	–	if (n < QT_MAX_LEVELS)
319	–	{
320	–	/* If set size exceeds threshold: subdivide cell and
321	–	reinsert set tris into cell
322	–	*/
323	–	/* CAUTION:If QT_SET_THRESHOLD << QT_MAXSET, and dont add
324	–	more than a few triangles before expanding: then are safe here
325	–	otherwise must check to make sure set size is < MAXSET,
326	–	or  qtgetset could overflow os.
327	–	*/
328	–	tptr = qtqueryset(qt);
329	–	if(QT_SET_CNT(tptr) > QT_MAXSET)
330	–	tptr = (OBJECT *)malloc((QT_SET_CNT(tptr)+1)*sizeof(OBJECT));
331	–	else
332	–	tptr = tset;
333	–	if(!tptr)
334	–	goto memerr;
335	–
336	–	qtgetset(tptr,qt);
337	–	n++;
338	–	qtfreeleaf(qt);
339	–	qtSubdivide(qt);
340	–	qt = qtAdd_tri(qt,q0,q1,q2,t0,t1,t2,id,n);
341	–
342	–	for(optr = QT_SET_PTR(tptr),i = QT_SET_CNT(tptr); i > 0; i--)
343	–	{
344	–	id = QT_SET_NEXT_ELEM(optr);
345	–	if(!qtTri_from_id(id,r0,r1,r2))
346	–	continue;
347	–	qt  = qtAdd_tri(qt,q0,q1,q2,r0,r1,r2,id,n);
348	–	}
349	–	if(tptr != tset)
350	–	free(tptr);
351	–	}
352	–	else
353	–	qt = qtaddelem(qt,id);
354	–	}
355	–	}
356	–	}
357	–	return(qt);
358	–	memerr:
359	–	error(SYSTEM,"qtAdd_tri():Unable to allocate memory");
360	–	}
361	–
362	–
363	–	QUADTREE
364	–	qtRemove_tri(qt,id,q0,q1,q2,t0,t1,t2)
365	–	QUADTREE qt;
366	–	int id;
367	–	FVECT q0,q1,q2;
368	–	FVECT t0,t1,t2;
369	–	{
370	–	FVECT a,b,c;
371	–
372	–	/* test if triangle (t0,t1,t2) overlaps cell triangle (v0,v1,v2) */
373	–	if(!stri_intersect(t0,t1,t2,q0,q1,q2))
374	–	return(qt);
375	–
376	–	/* if this is tree: recurse */
377	–	if(QT_IS_TREE(qt))
378	–	{
379	–	qtSubdivide_tri(q0,q1,q2,a,b,c);
380	–	QT_NTH_CHILD(qt,0) = qtRemove_tri(QT_NTH_CHILD(qt,0),id,t0,t1,t2,q0,a,c);
381	–	QT_NTH_CHILD(qt,1) = qtRemove_tri(QT_NTH_CHILD(qt,1),id,t0,t1,t2,a,q1,b);
382	–	QT_NTH_CHILD(qt,2) = qtRemove_tri(QT_NTH_CHILD(qt,2),id,t0,t1,t2,c,b,q2);
383	–	QT_NTH_CHILD(qt,3) = qtRemove_tri(QT_NTH_CHILD(qt,3),id,t0,t1,t2,b,c,a);
384	–	return(qt);
385	–	}
386	–	else
387	–	{
388	–	if(QT_IS_EMPTY(qt))
389	–	{
390	–	#ifdef DEBUG
391	–	eputs("qtRemove_tri(): triangle not found\n");
392	–	#endif
393	–	}
394	–	/* remove id from set */
395	–	else
396	–	{
397	–	if(!qtinset(qt,id))
398	–	{
399	–	#ifdef DEBUG
400	–	eputs("qtRemove_tri(): tri not in set\n");
401	–	#endif
402	–	}
403	–	else
404	–	qt = qtdelelem(qt,id);
405	–	}
406	–	}
407	–	return(qt);
408	–	}
409	–
410	–
411	–	QUADTREE
412	–	qtVisit_tri_interior(qt,q0,q1,q2,t0,t1,t2,n0,n1,n2,n,func,f,argptr)
413	–	QUADTREE qt;
414	–	FVECT q0,q1,q2;
415	–	FVECT t0,t1,t2;
416	–	FVECT n0,n1,n2;
417	–	int n;
418	–	int (*func)(),*f;
419	–	int *argptr;
420	–	{
421	–	FVECT a,b,c;
422	–
423	–	/* If qt Flag set, or qt vertices interior to t0t1t2-descend */
424	–	tree_modified:
425	–
426	–	if(QT_IS_TREE(qt))
427	–	{
428	–	if(QT_IS_FLAG(qt) ||  point_in_stri_n(n0,n1,n2,q0))
429	–	{
430	–	QT_SET_FLAG(qt);
431	–	qtSubdivide_tri(q0,q1,q2,a,b,c);
432	–	/* descend to children */
433	–	QT_NTH_CHILD(qt,0) = qtVisit_tri_interior(QT_NTH_CHILD(qt,0),
434	–	q0,a,c,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
435	–	QT_NTH_CHILD(qt,1) = qtVisit_tri_interior(QT_NTH_CHILD(qt,1),
436	–	a,q1,b,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
437	–	QT_NTH_CHILD(qt,2) = qtVisit_tri_interior(QT_NTH_CHILD(qt,2),
438	–	c,b,q2,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
439	–	QT_NTH_CHILD(qt,3) = qtVisit_tri_interior(QT_NTH_CHILD(qt,3),
440	–	b,c,a,t0,t1,t2,n0,n1,n2,n+1,func,f,argptr);
441	–	}
442	–	}
443	–	else
444	–	if((!QT_IS_EMPTY(qt) && QT_LEAF_IS_FLAG(qt)) ||
445	–	point_in_stri_n(n0,n1,n2,q0))
446	–	{
447	–	func(&qt,f,argptr,q0,q1,q2,t0,t1,t2,n);
448	–	if(QT_FLAG_IS_MODIFIED(*f))
449	–	{
450	–	QT_SET_FLAG(qt);
451	–	goto tree_modified;
452	–	}
453	–	if(QT_IS_LEAF(qt))
454	–	QT_LEAF_SET_FLAG(qt);
455	–	else
456	–	if(QT_IS_TREE(qt))
457	–	QT_SET_FLAG(qt);
458	–	}
459	–	return(qt);
460	–	}
461	–
462	–
463	–
136		int
137	<	move_to_nbri(b,db0,db1,db2,tptr)
137	>	move_to_nbr(b,db0,db1,db2,tptr)
138		BCOORD b[3];
139		BDIR db0,db1,db2;
140	<	TINT *tptr;
140	>	double *tptr;
141		{
142	<	TINT t,dt;
142	>	double t,dt;
143		BCOORD bc;
144		int nbr;
145
146		nbr = -1;
147	<	*tptr = 0;
147	>	*tptr = 0.0;
148		/* Advance to next node */
149		if(b[0]==0 && db0 < 0)
150		return(0);
#	Line 480 | Line 152 | TINT *tptr;
152		return(1);
153		if(b[2]==0 && db2 < 0)
154		return(2);
483	–
155		if(db0 < 0)
156		{
157	<	bc = b[0]<<SHIFT_MAXBCOORD;
487	<	t = bc/-db0;
157	>	t = ((double)b[0])/-db0;
158		nbr = 0;
159		}
160		else
161	<	t = HUGET;
161	>	t = MAXFLOAT;
162		if(db1 < 0)
163		{
164	<	bc = b[1] <<SHIFT_MAXBCOORD;
495	<	dt = bc/-db1;
164	>	dt = ((double)b[1])/-db1;
165		if( dt < t)
166		{
167		t = dt;
#	Line 501 | Line 170 | TINT *tptr;
170		}
171		if(db2 < 0)
172		{
173	<	bc = b[2] << SHIFT_MAXBCOORD;
505	<	dt = bc/-db2;
173	>	dt = ((double)b[2])/-db2;
174		if( dt < t)
175		{
176		t = dt;
#	Line 513 | Line 181 | TINT *tptr;
181		return(nbr);
182		}
183
184	<	QUADTREE
517	<	qtVisit_tri_edges(qt,b,db0,db1,db2,db,wptr,nextptr,t,sign,sfactor,func,f,argptr)
184	>	qtTrace_ray(qt,b,db0,db1,db2,nextptr,sign,sfactor,func,f)
185		QUADTREE qt;
186		BCOORD b[3];
187	<	BDIR db0,db1,db2,db[3][3];
188	<	int *wptr,*nextptr;
522	<	TINT t[3];
187	>	BDIR db0,db1,db2;
188	>	int *nextptr;
189		int sign,sfactor;
190	<	int (*func)();
191	<	int *f,*argptr;
190	>	FUNC func;
191	>	int *f;
192		{
193		int i,found;
194		QUADTREE child;
195		int nbr,next,w;
196	<	TINT t_g,t_l,t_i,l;
196	>	double t;
197
198		if(QT_IS_TREE(qt))
199		{
200		/* Find the appropriate child and reset the coord */
201	<	i = baryi_child(b);
201	>	i = bary_child(b);
202
537	–	QT_SET_FLAG(qt);
538	–
203		for(;;)
204		{
541	–	w = *wptr;
205		child = QT_NTH_CHILD(qt,i);
206		if(i != 3)
207	<	QT_NTH_CHILD(qt,i) =
545	<	qtVisit_tri_edges(child,b,db0,db1,db2,db,wptr,nextptr,t,sign,
546	<	sfactor+1,func,f,argptr);
207	>	qtTrace_ray(child,b,db0,db1,db2,nextptr,sign,sfactor+1,func,f);
208		else
209		/* If the center cell- must flip direction signs */
210	<	QT_NTH_CHILD(qt,i) =
550	<	qtVisit_tri_edges(child,b,-db0,-db1,-db2,db,wptr,nextptr,t,1-sign,
551	<	sfactor+1,func,f,argptr);
210	>	qtTrace_ray(child,b,-db0,-db1,-db2,nextptr,1-sign,sfactor+1,func,f);
211
212		if(QT_FLAG_IS_DONE(*f))
213	<	return(qt);
555	<	if(*wptr != w)
556	<	{
557	<	w = *wptr;
558	<	db0 = db[w][0];db1 = db[w][1];db2 = db[w][2];
559	<	if(sign)
560	<	{  db0 *= -1;db1 *= -1; db2 *= -1;}
561	<	}
213	>	return;
214		/* If in same block: traverse */
215		if(i==3)
216		next = *nextptr;
#	Line 568 | Line 220 | qtVisit_tri_edges(qt,b,db0,db1,db2,db,wptr,nextptr,t,s
220		else
221		{
222		/* reset the barycentric coordinates in the parents*/
223	<	baryi_parent(b,i);
223	>	bary_parent(b,i);
224		/* Else pop up to parent and traverse from there */
225		return(qt);
226		}
227	<	baryi_from_child(b,i,next);
227	>	bary_from_child(b,i,next);
228		i = next;
229		}
230		}
#	Line 582 | Line 234 | qtVisit_tri_edges(qt,b,db0,db1,db2,db,wptr,nextptr,t,s
234		if(Pick_cnt < 500)
235		Pick_q[Pick_cnt++]=qt;
236		#endif;
237	<	func(&qt,f,argptr);
237	>	F_FUNC(func)(qt,F_ARGS(func),f);
238		if(QT_FLAG_IS_DONE(*f))
239	+	return;
240	+	/* Advance to next node */
241	+	nbr = move_to_nbr(b,db0,db1,db2,&t);
242	+
243	+	if(nbr==-1)
244		{
245	<	if(!QT_IS_EMPTY(qt))
246	<	QT_LEAF_SET_FLAG(qt);
590	<	return(qt);
245	>	QT_FLAG_SET_DONE(*f);
246	>	return;
247		}
248	+	b[0] += (BCOORD)(t *db0);
249	+	b[1] += (BCOORD)(t *db1);
250	+	b[2] += (BCOORD)(t *db2);
251	+	*nextptr = nbr;
252	+	return;
253
254	<	if(!QT_IS_EMPTY(qt))
255	<	QT_LEAF_SET_FLAG(qt);
256	<	/* Advance to next node */
257	<	w = *wptr;
258	<	while(1)
259	<	{
260	<	nbr = move_to_nbri(b,db0,db1,db2,&t_i);
261	<
262	<	t_g = t_i >> sfactor;
263	<
264	<	if(t_g >= t[w])
265	<	{
266	<	if(w == 2)
267	<	{
268	<	QT_FLAG_SET_DONE(*f);
269	<	return(qt);
270	<	}
271	<	/* The edge fits in the cell- therefore the result of shifting
272	<	db by sfactor is guaranteed to be less than MAXBCOORD
273	<	*/
274	<	/* Caution: (t[w]*db) must occur before divide by MAXBCOORD
275	<	since t[w] will always be < MAXBCOORD
276	<	*/
277	<	l = t[w] << sfactor;
278	<	/* NOTE: Change divides to Shift and multiply by sign*/
279	<	b[0] += (l*db0)/MAXBCOORD;
280	<	b[1] += (l*db1)/MAXBCOORD;
281	<	b[2] += (l*db2)/MAXBCOORD;
282	<	w++;
283	<	db0 = db[w][0]; db1 = db[w][1]; db2 = db[w][2];
284	<	if(sign)
285	<	{  db0 *= -1;db1 *= -1; db2 *= -1;}
286	<	}
287	<	else
288	<	{
289	<	/* Caution: (t_i*db) must occur before divide by MAXBCOORD
290	<	since t_i will always be < MAXBCOORD*/
291	<	/* NOTE: Change divides to Shift and  by sign*/
292	<	b[0] += (t_i *db0) / MAXBCOORD;
293	<	b[1] += (t_i *db1) / MAXBCOORD;
294	<	b[2] += (t_i *db2) / MAXBCOORD;
295	<
296	<	t[w] -= t_g;
297	<	*wptr = w;
298	<	*nextptr = nbr;
299	<	return(qt);
300	<	}
301	<	}
302	<	}
254	>	}
255	>	}
256	>
257	>
258	>
259	>
260	>
261	>
262	>
263	>
264	>	#define TEST_INT(tl,th,d,q0,q1,h,l) \
265	>	tl=d>q0;th=d<q1; if(tl&&th) goto Lfunc_call; \
266	>	if(tl) if(l) goto Lfunc_call; else h=TRUE; \
267	>	if(th) if(h) goto Lfunc_call; else l = TRUE; \
268	>	if(th) if(h) goto Lfunc_call; else l = TRUE;
269	>
270	>	/* Leaf node: Do definitive test */
271	>	QUADTREE
272	>	qtLeaf_insert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
273	>	scale, s0,s1,s2,sq0,sq1,sq2,f,n)
274	>	int root;
275	>	QUADTREE qt;
276	>	BCOORD q0[3],q1[3],q2[3];
277	>	BCOORD t0[3],t1[3],t2[3];
278	>	BCOORD dt10[3],dt21[3],dt02[3];
279	>	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
280	>	FUNC f;
281	>	int n;
282	>	{
283	>	double db;
284	>	unsigned int e0,e1,e2;
285	>	char al,ah,bl,bh,cl,ch,testl,testh;
286	>	char test_t0t1,test_t1t2,test_t2t0;
287	>	BCOORD a,b,c;
288	>
289	>	/* First check if can trivial accept: if vertex in cell */
290	>	if(s0 & s1 & s2)
291	>	{
292	>	goto Lfunc_call;
293	>	}
294	>	/* Assumption: Could not trivial reject*/
295	>	/* IF any coords lie on edge- must be in if couldnt reject */
296	>	a = q1[0];b= q0[1]; c= q0[2];
297	>	if(t0[0]== a || t0[1] == b || t0[2] == c)
298	>	{
299	>	goto Lfunc_call;
300	>	}
301	>	if(t1[0]== a || t1[1] == b || t1[2] == c)
302	>	{
303	>	goto Lfunc_call;
304	>	}
305	>	if(t2[0]== a || t2[1] == b || t2[2] == c)
306	>	{
307	>	goto Lfunc_call;
308	>	}
309	>	/* Test for edge crossings */
310	>	/* Test t0t1,t1t2,t2t0 against a */
311	>	/* Calculate edge crossings */
312	>	e0  = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
313	>	/* See if edge can be trivially rejected from intersetion testing */
314	>	test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
315	>	((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
316	>	bl=bh=0;
317	>	if(test_t0t1 && (e0 & 2) )
318	>	{
319	>	/* Must do double calculation to avoid overflow */
320	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
321	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
322	>	}
323	>	test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
324	>	((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
325	>	if(test_t1t2 && (e0 & 1))
326	>	{    /* test t1t2 against a */
327	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
328	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
329	>	}
330	>	test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
331	>	((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
332	>	if(test_t2t0 && (e0 & 4))
333	>	{
334	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
335	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
336	>	}
337	>	e1  = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
338	>	cl = ch = 0;
339	>	if(test_t0t1 && (e1 & 2))
340	>	{/* test t0t1 against b */
341	>	db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
342	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
343	>	}
344	>	if(test_t1t2 && (e1 & 1))
345	>	{/* test t1t2 against b */
346	>	db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
347	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
348	>	}
349	>	if(test_t2t0 && (e1 & 4))
350	>	{/* test t2t0 against b */
351	>	db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
352	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
353	>	}
354	>
355	>	/* test edges against c */
356	>	e2  = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
357	>	al = ah = 0;
358	>	if(test_t0t1 && (e2 & 2))
359	>	{ /* test t0t1 against c */
360	>	db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
361	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
362	>	}
363	>	if(test_t1t2 && (e2 & 1))
364	>	{
365	>	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
366	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
367	>	}
368	>	if(test_t2t0 && (e2 & 4))
369	>	{ /* test t2t0 against c */
370	>	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
371	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
372	>	}
373	>	/* Only remaining case is if some edge trivially rejected */
374	>	if(!e0 || !e1 || !e2)
375	>	return(qt);
376	>
377	>	/* Only one/none got tested - pick either of other two/any two */
378	>	/* Only need to test one edge */
379	>	if(!test_t0t1 && (e0 & 2))
380	>	{
381	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
382	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
383	>	}
384	>	if(!test_t1t2 && (e0 & 1))
385	>	{/* test t1t2 against a */
386	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
387	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
388	>	}
389	>	if(!test_t2t0 && (e0 & 4))
390	>	{
391	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
392	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
393	>	}
394	>
395	>	return(qt);
396	>
397	>	Lfunc_call:
398	>	qt = f.func(f.argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
399	>	s0,s1,s2,sq0,sq1,sq2,0,f,n);
400	>	return(qt);
401	>
402		}
403
404
405	+
406	+	/* Leaf node: Do definitive test */
407		QUADTREE
408	<	qtRoot_visit_tri_edges(qt,q0,q1,q2,peq,tri,i_pt,wptr,nextptr,func,f,argptr)
409	<	QUADTREE qt;
410	<	FVECT q0,q1,q2;
411	<	FPEQ  peq;
412	<	FVECT tri[3],i_pt;
413	<	int *wptr,*nextptr;
414	<	int (*func)();
415	<	int *f,*argptr;
416	<	{
417	<	int x,y,z,w,i,j,first;
418	<	QUADTREE child;
419	<	FVECT c,d,v[3];
420	<	double b[4][3],db[3][3],et[3],exit_pt;
421	<	BCOORD bi[3];
422	<	TINT t[3];
423	<	BDIR dbi[3][3];
408	>	qtLeaf_insert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
409	>	scale, s0,s1,s2,sq0,sq1,sq2,f,n)
410	>	int root;
411	>	QUADTREE qt;
412	>	BCOORD q0[3],q1[3],q2[3];
413	>	BCOORD t0[3],t1[3],t2[3];
414	>	BCOORD dt10[3],dt21[3],dt02[3];
415	>	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
416	>	FUNC f;
417	>	int n;
418	>	{
419	>	double db;
420	>	unsigned int e0,e1,e2;
421	>	BCOORD a,b,c;
422	>	double p0[2], p1[2],cp;
423	>	char al,ah,bl,bh,cl,ch;
424	>	char testl,testh,test_t0t1,test_t1t2,test_t2t0;
425	>	unsigned int ls0,ls1,ls2;
426
427	<	first =0;
428	<	w = *wptr;
429	<	if(w==-1)
427	>
428	>	/* May have gotten passed trivial reject if one/two vertices are ON */
429	>	a = q1[0];b= q0[1]; c= q0[2];
430	>	SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
431	>
432	>	/* First check if can trivial accept: if vertex in cell */
433	>	if(ls0 & ls1 & ls2)
434		{
435	<	first = 1;
668	<	*wptr = w = 0;
435	>	goto Lfunc_call;
436		}
437	<	/* Project the origin onto the root node plane */
437	>	if(ls0==0 || ls1 == 0 || ls2 ==0)
438	>	return(qt);
439	>	/* Assumption: Could not trivial reject*/
440	>	/* IF any coords lie on edge- must be in if couldnt reject */
441
442	<	t[0] = t[1] = t[2] = 0;
443	<	/* Find the intersection point of the origin */
444	<	/* map to 2d by dropping maximum magnitude component of normal */
442	>	if(t0[0]== a || t0[1] == b || t0[2] == c)
443	>	{
444	>	goto Lfunc_call;
445	>	}
446	>	if(t1[0]== a || t1[1] == b || t1[2] == c)
447	>	{
448	>	goto Lfunc_call;
449	>	}
450	>	if(t2[0]== a || t2[1] == b || t2[2] == c)
451	>	{
452	>	goto Lfunc_call;
453	>	}
454	>	/* Test for edge crossings */
455	>	/* Test t0t1 against a,b,c */
456	>	/* First test if t0t1 can be trivially rejected */
457	>	/* If both edges are outside bounds- dont need to test */
458	>
459	>	/* Test t0t1,t1t2,t2t0 against a */
460	>	test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
461	>	((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
462	>	e0  = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
463	>	bl=bh=0;
464	>	/* Test t0t1,t1t2,t2t0 against a */
465	>	if(test_t0t1 && (e0 & 2) )
466	>	{
467	>	db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
468	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
469	>	}
470	>	test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
471	>	((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
472	>	if(test_t1t2 && (e0 & 1))
473	>	{    /* test t1t2 against a */
474	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
475	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
476	>	}
477	>	test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
478	>	((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
479	>	if(test_t2t0 && (e0 & 4))
480	>	{
481	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
482	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
483	>	}
484	>	cl = ch = 0;
485	>	e1  = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
486	>	if(test_t0t1 && (e1 & 2))
487	>	{/* test t0t1 against b */
488	>	db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
489	>	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
490	>	}
491	>	if(test_t1t2 && (e1 & 1))
492	>	{/* test t1t2 against b */
493	>	db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
494	>	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
495	>	}
496	>	if(test_t2t0 && (e1 & 4))
497	>	{/* test t2t0 against b */
498	>	db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
499	>	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
500	>	}
501	>	al = ah = 0;
502	>	e2  = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
503	>	if(test_t0t1 && (e2 & 2))
504	>	{ /* test t0t1 against c */
505	>	db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
506	>	TEST_INT(testl,testh,db,q0[0],a,al,ah)
507	>	}
508	>	if(test_t1t2 && (e2 & 1))
509	>	{
510	>	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
511	>	TEST_INT(testl,testh,db,q0[0],a,al,ah)
512	>	}
513	>	if(test_t2t0 && (e2 & 4))
514	>	{ /* test t2t0 against c */
515	>	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
516	>	TEST_INT(testl,testh,db,q0[0],a,al,ah)
517	>	}
518	>	/* Only remaining case is if some edge trivially rejected */
519	>	if(!e0 || !e1 || !e2)
520	>	return(qt);
521
522	<	x = FP_X(peq);
523	<	y = FP_Y(peq);
524	<	z = FP_Z(peq);
679	<	/* Calculate barycentric coordinates for current vertex */
680	<	if(!first)
681	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);
682	<	else
683	<	/* Just starting: b[0] is the origin point: guaranteed to be valid b*/
522	>	/* Only one/none got tested - pick either of other two/any two */
523	>	/* Only need to test one edge */
524	>	if(!test_t0t1 && (e0 & 2))
525		{
526	<	intersect_vector_plane(tri[0],peq,&(et[0]),v[0]);
527	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[0][x],v[0][y],b[0]);
687	<	VCOPY(b[3],b[0]);
526	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
527	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
528		}
529	+	if(!test_t1t2 && (e0 & 1))
530	+	{/* test t1t2 against a */
531	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
532	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
533	+	}
534	+	if(!test_t2t0 && (e0 & 4))
535	+	{
536	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
537	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
538	+	}
539
540	<	j = (w+1)%3;
541	<	intersect_vector_plane(tri[j],peq,&(et[j]),v[j]);
542	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);
543	<	if(et[j] < 0.0)
540	>	return(qt);
541	>	Lfunc_call:
542	>	qt = f.func(f.argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
543	>	s0,s1,s2,sq0,sq1,sq2,1,f,n);
544	>	return(qt);
545	>	}
546	>
547	>
548	>
549	>	/* ASSUMPTION: that triangle t0,t1,t2 intersects quad cell q0,q1,q2 */
550	>
551	>	/*-------q2--------- sq2
552	>	/ \
553	>	s1     /sc \ s0
554	>	qb_____qa
555	>	/ \   / \
556	>	\sq0/sa \ /sb \   /sq1
557	>	\ q0____qc____q1/
558	>	\             /
559	>	\     s2    /
560	>	*/
561	>
562	>	int Find_id=0;
563	>
564	>	QUADTREE
565	>	qtInsert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
566	>	s0,s1,s2,sq0,sq1,sq2,f,n)
567	>	int root;
568	>	QUADTREE qt;
569	>	BCOORD q0[3],q1[3],q2[3];
570	>	BCOORD t0[3],t1[3],t2[3];
571	>	BCOORD dt10[3],dt21[3],dt02[3];
572	>	BCOORD scale;
573	>	unsigned int s0,s1,s2,sq0,sq1,sq2;
574	>	FUNC f;
575	>	int n;
576	>	{
577	>	BCOORD a,b,c;
578	>	BCOORD va[3],vb[3],vc[3];
579	>	unsigned int sa,sb,sc;
580	>
581	>	/* If a tree: trivial reject and recurse on potential children */
582	>	if(QT_IS_TREE(qt))
583		{
584	<	VSUB(db[w],b[3],b[j]);
585	<	t[w] = HUGET;
584	>	/* Test against new a edge: if entirely to left: only need
585	>	to visit child 0
586	>	*/
587	>	a = q1[0] + scale;
588	>	b = q0[1] + scale;
589	>	c = q0[2] + scale;
590	>
591	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
592	>
593	>	if(sa==7)
594	>	{
595	>	vb[1] = q0[1];
596	>	vc[2] = q0[2];
597	>	vc[1] = b;
598	>	vb[2] = c;
599	>	vb[0] = vc[0] = a;
600	>	QT_NTH_CHILD(qt,0) = qtInsert_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
601	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
602	>	return(qt);
603	>	}
604	>	if(sb==7)
605	>	{
606	>	va[0] = q1[0];
607	>	vc[2] = q0[2];
608	>	va[1] = vc[1] = b;
609	>	va[2] = c;
610	>	vc[0] = a;
611	>	QT_NTH_CHILD(qt,1) = qtInsert_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,
612	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
613	>	return(qt);
614	>	}
615	>	if(sc==7)
616	>	{
617	>	va[0] = q1[0];
618	>	vb[1] = q0[1];
619	>	va[1] = b;
620	>	va[2] = vb[2] = c;
621	>	vb[0] = a;
622	>	QT_NTH_CHILD(qt,2) = qtInsert_tri(root,QT_NTH_CHILD(qt,2),vb,va,
623	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
624	>	return(qt);
625	>	}
626	>
627	>	va[0] = q1[0];
628	>	vb[1] = q0[1];
629	>	vc[2] = q0[2];
630	>	va[1] = vc[1] = b;
631	>	va[2] = vb[2] = c;
632	>	vb[0] = vc[0] = a;
633	>	/* If outside of all edges: only need to Visit 3  */
634	>	if(sa==0 && sb==0 && sc==0)
635	>	{
636	>	QT_NTH_CHILD(qt,3) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,
637	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
638	>	return(qt);
639	>	}
640	>
641	>	if(sa)
642	>	QT_NTH_CHILD(qt,0) = qtInsert_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
643	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
644	>	if(sb)
645	>	QT_NTH_CHILD(qt,1) = qtInsert_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
646	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
647	>	if(sc)
648	>	QT_NTH_CHILD(qt,2) = qtInsert_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
649	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
650	>	QT_NTH_CHILD(qt,3) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
651	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
652	>	return(qt);
653		}
654	+	/* Leaf node: Do definitive test */
655		else
656	<	{
657	<	/* NOTE: for stability: do not increment with ipt- use full dir and
658	<	calculate t: but for wrap around case: could get same problem?
659	<	*/
660	<	VSUB(db[w],b[j],b[3]);
661	<	/* Check if the point is out side of the triangle: if so t[w] =HUGET */
662	<	if((fabs(b[j][0])>(FTINY+1.0)) ||(fabs(b[j][1])>(FTINY+1.0)) ||
663	<	(fabs(b[j][2])>(FTINY+1.0))||(b[j][0] <-FTINY) ||
664	<	(b[j][1]<-FTINY) ||(b[j][2]<-FTINY))
665	<	t[w] = HUGET;
666	<	else
656	>	return(qt = qtLeaf_insert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
657	>	scale,s0,s1,s2,sq0,sq1,sq2,f,n));
658	>	}
659	>
660	>
661	>	QUADTREE
662	>	qtInsert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
663	>	s0,s1,s2,sq0,sq1,sq2,f,n)
664	>	int root;
665	>	QUADTREE qt;
666	>	BCOORD q0[3],q1[3],q2[3];
667	>	BCOORD t0[3],t1[3],t2[3];
668	>	BCOORD dt10[3],dt21[3],dt02[3];
669	>	BCOORD scale;
670	>	unsigned int s0,s1,s2,sq0,sq1,sq2;
671	>	FUNC f;
672	>	int n;
673	>	{
674	>	BCOORD a,b,c;
675	>	BCOORD va[3],vb[3],vc[3];
676	>	unsigned int sa,sb,sc;
677	>
678	>	/* If a tree: trivial reject and recurse on potential children */
679	>	if(QT_IS_TREE(qt))
680	>	{
681	>	/* Test against new a edge: if entirely to left: only need
682	>	to visit child 0
683	>	*/
684	>	a = q1[0] - scale;
685	>	b = q0[1] - scale;
686	>	c = q0[2] - scale;
687	>
688	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
689	>
690	>	if(sa==0)
691	>	{
692	>	vb[1] = q0[1];
693	>	vc[2] = q0[2];
694	>	vc[1] = b;
695	>	vb[2] = c;
696	>	vb[0] = vc[0] = a;
697	>
698	>	QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
699	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
700	>	return(qt);
701	>	}
702	>	if(sb==0)
703	>	{
704	>	va[0] = q1[0];
705	>	vc[2] = q0[2];
706	>	va[1] = vc[1] = b;
707	>	va[2] = c;
708	>	vc[0] = a;
709	>	QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
710	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
711	>	return(qt);
712	>	}
713	>	if(sc==0)
714	>	{
715	>	va[0] = q1[0];
716	>	vb[1] = q0[1];
717	>	va[1] = b;
718	>	va[2] = vb[2] = c;
719	>	vb[0] = a;
720	>	QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
721	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
722	>	return(qt);
723	>	}
724	>	va[0] = q1[0];
725	>	vb[1] = q0[1];
726	>	vc[2] = q0[2];
727	>	va[1] = vc[1] = b;
728	>	va[2] = vb[2] = c;
729	>	vb[0] = vc[0] = a;
730	>	/* If outside of all edges: only need to Visit 3  */
731	>	if(sa==7 && sb==7 && sc==7)
732	>	{
733	>	QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,
734	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
735	>	return(qt);
736	>	}
737	>	if(sa!=7)
738	>	QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
739	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
740	>	if(sb!=7)
741	>	QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
742	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
743	>	if(sc!=7)
744	>	QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
745	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
746	>
747	>	QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
748	>	t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
749	>	return(qt);
750	>	}
751	>	/* Leaf node: Do definitive test */
752	>	else
753	>	return(qt = qtLeaf_insert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
754	>	scale,s0,s1,s2,sq0,sq1,sq2,f,n));
755	>	}
756	>
757	>
758	>
759	>
760	>	/*************************************************************************/
761	>	/* Visit code for applying functions: NOTE THIS IS THE SAME AS INSERT CODE:
762	>	except sets flags: wanted insert to be as efficient as possible: so
763	>	broke into two sets of routines. Also traverses only to n levels.
764	>	*/
765	>
766	>	qtVisit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
767	>	s0,s1,s2,sq0,sq1,sq2,f,n)
768	>	int root;
769	>	QUADTREE qt;
770	>	BCOORD q0[3],q1[3],q2[3];
771	>	BCOORD t0[3],t1[3],t2[3];
772	>	BCOORD dt10[3],dt21[3],dt02[3];
773	>	BCOORD scale;
774	>	unsigned int s0,s1,s2,sq0,sq1,sq2;
775	>	FUNC f;
776	>	int n;
777	>	{
778	>	BCOORD a,b,c;
779	>	BCOORD va[3],vb[3],vc[3];
780	>	unsigned int sa,sb,sc;
781	>
782	>	if(n == 0)
783	>	return;
784	>	/* If a tree: trivial reject and recurse on potential children */
785	>	if(QT_IS_TREE(qt))
786	>	{
787	>	QT_SET_FLAG(qt);
788	>
789	>	/* Test against new a edge: if entirely to left: only need
790	>	to visit child 0
791	>	*/
792	>	a = q1[0] + scale;
793	>	b = q0[1] + scale;
794	>	c = q0[2] + scale;
795	>
796	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
797	>
798	>	if(sa==7)
799	>	{
800	>	vb[1] = q0[1];
801	>	vc[2] = q0[2];
802	>	vc[1] = b;
803	>	vb[2] = c;
804	>	vb[0] = vc[0] = a;
805	>	qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
806	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
807	>	return;
808	>	}
809	>	if(sb==7)
810		{
811	<	/* The next point is in the triangle- so db values will be in valid
812	<	range and t[w]= MAXT
813	<	*/
814	<	t[w] = MAXT;
815	<	if(j != 0)
816	<	for(;j < 3;j++)
817	<	{
818	<	i = (j+1)%3;
719	<	if(!first || i != 0)
720	<	{
721	<	intersect_vector_plane(tri[i],peq,&(et[i]),v[i]);
722	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
723	<	v[i][y],b[i]);
724	<	}
725	<	if(et[i] < 0.0)
726	<	{
727	<	VSUB(db[j],b[j],b[i]);
728	<	t[j] = HUGET;
729	<	break;
730	<	}
731	<	else
732	<	{
733	<	VSUB(db[j],b[i],b[j]);
734	<	if((fabs(b[j][0])>(FTINY+1.0))||(fabs(b[j][1])>(FTINY+1.0)) ||
735	<	(fabs(b[j][2])>(FTINY+1.0))||(b[i][0] <-FTINY) ||
736	<	(b[i][1]<-FTINY) ||(b[i][2]<-FTINY))
737	<	{
738	<	t[j] = HUGET;
739	<	break;
740	<	}
741	<	else
742	<	t[j] = MAXT;
743	<	}
744	<	}
811	>	va[0] = q1[0];
812	>	vc[2] = q0[2];
813	>	va[1] = vc[1] = b;
814	>	va[2] = c;
815	>	vc[0] = a;
816	>	qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,
817	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
818	>	return;
819		}
820	<	}
821	<	bary_dtol(b[3],db,bi,dbi,t,w);
820	>	if(sc==7)
821	>	{
822	>	va[0] = q1[0];
823	>	vb[1] = q0[1];
824	>	va[1] = b;
825	>	va[2] = vb[2] = c;
826	>	vb[0] = a;
827	>	qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,
828	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
829	>	return;
830	>	}
831	>
832	>	va[0] = q1[0];
833	>	vb[1] = q0[1];
834	>	vc[2] = q0[2];
835	>	va[1] = vc[1] = b;
836	>	va[2] = vb[2] = c;
837	>	vb[0] = vc[0] = a;
838	>	/* If outside of all edges: only need to Visit 3  */
839	>	if(sa==0 && sb==0 && sc==0)
840	>	{
841	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,
842	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
843	>	return;
844	>	}
845	>
846	>	if(sa)
847	>	qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
848	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
849	>	if(sb)
850	>	qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
851	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
852	>	if(sc)
853	>	qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
854	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
855	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
856	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
857	>	}
858	>	/* Leaf node: Do definitive test */
859	>	else
860	>	qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
861	>	scale,s0,s1,s2,sq0,sq1,sq2,f,n);
862	>
863	>	}
864	>
865	>
866	>	qtVisit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
867	>	s0,s1,s2,sq0,sq1,sq2,f,n)
868	>	int root;
869	>	QUADTREE qt;
870	>	BCOORD q0[3],q1[3],q2[3];
871	>	BCOORD t0[3],t1[3],t2[3];
872	>	BCOORD dt10[3],dt21[3],dt02[3];
873	>	BCOORD scale;
874	>	unsigned int s0,s1,s2,sq0,sq1,sq2;
875	>	FUNC f;
876	>	int n;
877	>	{
878	>	BCOORD a,b,c;
879	>	BCOORD va[3],vb[3],vc[3];
880	>	unsigned int sa,sb,sc;
881	>
882	>	if(n==0)
883	>	return;
884	>	/* If a tree: trivial reject and recurse on potential children */
885	>	if(QT_IS_TREE(qt))
886	>	{
887	>	QT_SET_FLAG(qt);
888	>	/* Test against new a edge: if entirely to left: only need
889	>	to visit child 0
890	>	*/
891	>	a = q1[0] - scale;
892	>	b = q0[1] - scale;
893	>	c = q0[2] - scale;
894	>
895	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
896	>
897	>	if(sa==0)
898	>	{
899	>	vb[1] = q0[1];
900	>	vc[2] = q0[2];
901	>	vc[1] = b;
902	>	vb[2] = c;
903	>	vb[0] = vc[0] = a;
904	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
905	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
906	>	return;
907	>	}
908	>	if(sb==0)
909	>	{
910	>	va[0] = q1[0];
911	>	vc[2] = q0[2];
912	>	va[1] = vc[1] = b;
913	>	va[2] = c;
914	>	vc[0] = a;
915	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
916	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
917	>	return;
918	>	}
919	>	if(sc==0)
920	>	{
921	>	va[0] = q1[0];
922	>	vb[1] = q0[1];
923	>	va[1] = b;
924	>	va[2] = vb[2] = c;
925	>	vb[0] = a;
926	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
927	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
928	>	return;
929	>	}
930	>	va[0] = q1[0];
931	>	vb[1] = q0[1];
932	>	vc[2] = q0[2];
933	>	va[1] = vc[1] = b;
934	>	va[2] = vb[2] = c;
935	>	vb[0] = vc[0] = a;
936	>	/* If outside of all edges: only need to Visit 3  */
937	>	if(sa==7 && sb==7 && sc==7)
938	>	{
939	>	qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,
940	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
941	>	return;
942	>	}
943	>	if(sa!=7)
944	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
945	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n-1);
946	>	if(sb!=7)
947	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
948	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n-1);
949	>	if(sc!=7)
950	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
951	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n-1);
952	>
953	>	qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
954	>	t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n-1);
955	>	return;
956	>	}
957	>	/* Leaf node: Do definitive test */
958	>	else
959	>	qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
960	>	scale,s0,s1,s2,sq0,sq1,sq2,f,n);
961	>	}
962	>
963	>
964	>
965	>	qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
966	>	scale, s0,s1,s2,sq0,sq1,sq2,f,n)
967	>	int root;
968	>	QUADTREE qt;
969	>	BCOORD q0[3],q1[3],q2[3];
970	>	BCOORD t0[3],t1[3],t2[3];
971	>	BCOORD dt10[3],dt21[3],dt02[3];
972	>	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
973	>	FUNC f;
974	>	int n;
975	>	{
976	>	double db;
977	>	unsigned int e0,e1,e2;
978	>	char al,ah,bl,bh,cl,ch,testl,testh;
979	>	char test_t0t1,test_t1t2,test_t2t0;
980	>	BCOORD a,b,c;
981
982	<	/* trace the ray starting with this node */
983	<	qt = qtVisit_tri_edges(qt,bi,dbi[w][0],dbi[w][1],dbi[w][2],
984	<	dbi,wptr,nextptr,t,0,0,func,f,argptr);
985	<	if(!QT_FLAG_IS_DONE(*f))
982	>	/* First check if can trivial accept: if vertex in cell */
983	>	if(s0 & s1 & s2)
984	>	goto Lfunc_call;
985	>
986	>	/* Assumption: Could not trivial reject*/
987	>	/* IF any coords lie on edge- must be in if couldnt reject */
988	>	a = q1[0];b= q0[1]; c= q0[2];
989	>	if(t0[0]== a || t0[1] == b || t0[2] == c)
990	>	goto Lfunc_call;
991	>	if(t1[0]== a || t1[1] == b || t1[2] == c)
992	>	goto Lfunc_call;
993	>	if(t2[0]== a || t2[1] == b || t2[2] == c)
994	>	goto Lfunc_call;
995	>
996	>	/* Test for edge crossings */
997	>	/* Test t0t1,t1t2,t2t0 against a */
998	>	/* Calculate edge crossings */
999	>	e0  = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
1000	>	/* See if edge can be trivially rejected from intersetion testing */
1001	>	test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
1002	>	((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
1003	>	bl=bh=0;
1004	>	if(test_t0t1 && (e0 & 2) )
1005		{
1006	<	b[3][0] = (double)bi[0]/(double)MAXBCOORD;
1007	<	b[3][1] = (double)bi[1]/(double)MAXBCOORD;
1008	<	b[3][2] = (double)bi[2]/(double)MAXBCOORD;
757	<	i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
758	<	i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
759	<	i_pt[z] = (-FP_N(peq)[x]*i_pt[x] - FP_N(peq)[y]*i_pt[y]-FP_D(peq))/FP_N(peq)[z];
1006	>	/* Must do double calculation to avoid overflow */
1007	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1008	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1009		}
1010	<	return(qt);
1011	<
1010	>	test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
1011	>	((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
1012	>	if(test_t1t2 && (e0 & 1))
1013	>	{    /* test t1t2 against a */
1014	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1015	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1016	>	}
1017	>	test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
1018	>	((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
1019	>	if(test_t2t0 && (e0 & 4))
1020	>	{
1021	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1022	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1023	>	}
1024	>	e1  = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
1025	>	cl = ch = 0;
1026	>	if(test_t0t1 && (e1 & 2))
1027	>	{/* test t0t1 against b */
1028	>	db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
1029	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1030	>	}
1031	>	if(test_t1t2 && (e1 & 1))
1032	>	{/* test t1t2 against b */
1033	>	db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
1034	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1035	>	}
1036	>	if(test_t2t0 && (e1 & 4))
1037	>	{/* test t2t0 against b */
1038	>	db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
1039	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1040	>	}
1041	>
1042	>	/* test edges against c */
1043	>	e2  = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
1044	>	al = ah = 0;
1045	>	if(test_t0t1 && (e2 & 2))
1046	>	{ /* test t0t1 against c */
1047	>	db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
1048	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1049	>	}
1050	>	if(test_t1t2 && (e2 & 1))
1051	>	{
1052	>	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1053	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1054	>	}
1055	>	if(test_t2t0 && (e2 & 4))
1056	>	{ /* test t2t0 against c */
1057	>	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1058	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1059	>	}
1060	>	/* Only remaining case is if some edge trivially rejected */
1061	>	if(!e0 || !e1 || !e2)
1062	>	return;
1063	>
1064	>	/* Only one/none got tested - pick either of other two/any two */
1065	>	/* Only need to test one edge */
1066	>	if(!test_t0t1 && (e0 & 2))
1067	>	{
1068	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1069	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1070	>	}
1071	>	if(!test_t1t2 && (e0 & 1))
1072	>	{/* test t1t2 against a */
1073	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1074	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1075	>	}
1076	>	if(!test_t2t0 && (e0 & 4))
1077	>	{
1078	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1079	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1080	>	}
1081	>
1082	>	return;
1083	>
1084	>	Lfunc_call:
1085	>
1086	>	f.func(f.argptr,root,qt,n);
1087	>	if(!QT_IS_EMPTY(qt))
1088	>	QT_LEAF_SET_FLAG(qt);
1089	>
1090		}
1091
1092
1093	+
1094	+	/* Leaf node: Do definitive test */
1095		QUADTREE
1096	<	qtRoot_trace_ray(qt,q0,q1,q2,peq,orig,dir,nextptr,func,f,argptr)
1097	<	QUADTREE qt;
1098	<	FVECT q0,q1,q2;
1099	<	FPEQ  peq;
1100	<	FVECT orig,dir;
1101	<	int *nextptr;
1102	<	int (*func)();
1103	<	int *f,*argptr;
1104	<	{
1105	<	int x,y,z,nbr,w,i;
1106	<	QUADTREE child;
1107	<	FVECT v,i_pt;
1108	<	double b[2][3],db[3],et[2],d,br[3];
1109	<	BCOORD bi[3];
1110	<	TINT t[3];
1111	<	BDIR dbi[3][3];
1096	>	qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1097	>	scale, s0,s1,s2,sq0,sq1,sq2,f,n)
1098	>	int root;
1099	>	QUADTREE qt;
1100	>	BCOORD q0[3],q1[3],q2[3];
1101	>	BCOORD t0[3],t1[3],t2[3];
1102	>	BCOORD dt10[3],dt21[3],dt02[3];
1103	>	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
1104	>	FUNC f;
1105	>	int n;
1106	>	{
1107	>	double db;
1108	>	unsigned int e0,e1,e2;
1109	>	BCOORD a,b,c;
1110	>	double p0[2], p1[2],cp;
1111	>	char al,ah,bl,bh,cl,ch;
1112	>	char testl,testh,test_t0t1,test_t1t2,test_t2t0;
1113	>	unsigned int ls0,ls1,ls2;
1114
1115	<	/* Project the origin onto the root node plane */
1116	<	t[0] = t[1] = t[2] = 0;
1115	>	/* May have gotten passed trivial reject if one/two vertices are ON */
1116	>	a = q1[0];b= q0[1]; c= q0[2];
1117	>	SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
1118	>
1119	>	/* First check if can trivial accept: if vertex in cell */
1120	>	if(ls0 & ls1 & ls2)
1121	>	goto Lfunc_call;
1122
1123	<	VADD(v,orig,dir);
1124	<	/* Find the intersection point of the origin */
1125	<	/* map to 2d by dropping maximum magnitude component of normal */
1126	<	x = FP_X(peq);
791	<	y = FP_Y(peq);
792	<	z = FP_Z(peq);
1123	>	if(ls0==0 || ls1 == 0 || ls2 ==0)
1124	>	return;
1125	>	/* Assumption: Could not trivial reject*/
1126	>	/* IF any coords lie on edge- must be in if couldnt reject */
1127
1128	<	/* Calculate barycentric coordinates for current vertex */
1129	<	intersect_vector_plane(orig,peq,&(et[0]),i_pt);
1130	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[0]);
1128	>	if(t0[0]== a || t0[1] == b || t0[2] == c)
1129	>	goto Lfunc_call;
1130	>	if(t1[0]== a || t1[1] == b || t1[2] == c)
1131	>	goto Lfunc_call;
1132	>	if(t2[0]== a || t2[1] == b || t2[2] == c)
1133	>	goto Lfunc_call;
1134
1135	<	intersect_vector_plane(v,peq,&(et[1]),i_pt);
1136	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[1]);
1137	<	if(et[1] < 0.0)
1138	<	VSUB(db,b[0],b[1]);
1135	>	/* Test for edge crossings */
1136	>	/* Test t0t1 against a,b,c */
1137	>	/* First test if t0t1 can be trivially rejected */
1138	>	/* If both edges are outside bounds- dont need to test */
1139	>
1140	>	/* Test t0t1,t1t2,t2t0 against a */
1141	>	test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
1142	>	((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
1143	>	e0  = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
1144	>	bl=bh=0;
1145	>	/* Test t0t1,t1t2,t2t0 against a */
1146	>	if(test_t0t1 && (e0 & 2) )
1147	>	{
1148	>	db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
1149	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1150	>	}
1151	>	test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
1152	>	((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
1153	>	if(test_t1t2 && (e0 & 1))
1154	>	{    /* test t1t2 against a */
1155	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1156	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1157	>	}
1158	>	test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
1159	>	((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
1160	>	if(test_t2t0 && (e0 & 4))
1161	>	{
1162	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1163	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1164	>	}
1165	>	cl = ch = 0;
1166	>	e1  = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
1167	>	if(test_t0t1 && (e1 & 2))
1168	>	{/* test t0t1 against b */
1169	>	db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
1170	>	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1171	>	}
1172	>	if(test_t1t2 && (e1 & 1))
1173	>	{/* test t1t2 against b */
1174	>	db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
1175	>	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1176	>	}
1177	>	if(test_t2t0 && (e1 & 4))
1178	>	{/* test t2t0 against b */
1179	>	db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
1180	>	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1181	>	}
1182	>	al = ah = 0;
1183	>	e2  = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
1184	>	if(test_t0t1 && (e2 & 2))
1185	>	{ /* test t0t1 against c */
1186	>	db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
1187	>	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1188	>	}
1189	>	if(test_t1t2 && (e2 & 1))
1190	>	{
1191	>	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1192	>	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1193	>	}
1194	>	if(test_t2t0 && (e2 & 4))
1195	>	{ /* test t2t0 against c */
1196	>	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1197	>	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1198	>	}
1199	>	/* Only remaining case is if some edge trivially rejected */
1200	>	if(!e0 || !e1 || !e2)
1201	>	return;
1202	>
1203	>	/* Only one/none got tested - pick either of other two/any two */
1204	>	/* Only need to test one edge */
1205	>	if(!test_t0t1 && (e0 & 2))
1206	>	{
1207	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1208	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1209	>	}
1210	>	if(!test_t1t2 && (e0 & 1))
1211	>	{/* test t1t2 against a */
1212	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1213	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1214	>	}
1215	>	if(!test_t2t0 && (e0 & 4))
1216	>	{
1217	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1218	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1219	>	}
1220	>
1221	>	return;
1222	>
1223	>	Lfunc_call:
1224	>	f.func(f.argptr,root,qt,n);
1225	>	if(!QT_IS_EMPTY(qt))
1226	>	QT_LEAF_SET_FLAG(qt);
1227	>	}
1228	>
1229	>
1230	>	QUADTREE
1231	>	qtInsert_point(root,qt,parent,q0,q1,q2,bc,scale,f,n,doneptr)
1232	>	int root;
1233	>	QUADTREE qt,parent;
1234	>	BCOORD q0[3],q1[3],q2[3],bc[3],scale;
1235	>	FUNC f;
1236	>	int n,*doneptr;
1237	>	{
1238	>	BCOORD a,b,c;
1239	>	BCOORD va[3],vb[3],vc[3];
1240	>
1241	>	if(QT_IS_TREE(qt))
1242	>	{
1243	>	a = q1[0] + scale;
1244	>	b = q0[1] + scale;
1245	>	c = q0[2] + scale;
1246	>	if(bc[0] > a)
1247	>	{
1248	>	vc[0] = a; vc[1] = b; vc[2] = q0[2];
1249	>	vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1250	>	QT_NTH_CHILD(qt,0) = qtInsert_point(root,QT_NTH_CHILD(qt,0),
1251	>	qt,q0,vc,vb,bc,scale>>1,f,n+1,doneptr);
1252	>	if(!*doneptr)
1253	>	{
1254	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1255	>	if(!*doneptr)
1256	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1257	>	if(!*doneptr)
1258	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1259	>	}
1260	>	return(qt);
1261	>	}
1262	>	if(bc[1] > b)
1263	>	{
1264	>	vc[0] = a; vc[1] = b; vc[2] = q0[2];
1265	>	va[0] = q1[0]; va[1] = b; va[2] = c;
1266	>	QT_NTH_CHILD(qt,1) =qtInsert_point(root,QT_NTH_CHILD(qt,1),
1267	>	qt,vc,q1,va,bc,scale >>1,f,n+1,doneptr);
1268	>	if(!*doneptr)
1269	>	{
1270	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1271	>	if(!*doneptr)
1272	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1273	>	if(!*doneptr)
1274	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1275	>	}
1276	>	return(qt);
1277	>	}
1278	>	if(bc[2] > c)
1279	>	{
1280	>	va[0] = q1[0]; va[1] = b; va[2] = c;
1281	>	vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1282	>	QT_NTH_CHILD(qt,2) = qtInsert_point(root,QT_NTH_CHILD(qt,2),
1283	>	qt,vb,va,q2,bc,scale>>1,f,n+1,doneptr);
1284	>	if(!*doneptr)
1285	>	{
1286	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1287	>	if(!*doneptr)
1288	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1289	>	if(!*doneptr)
1290	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1291	>	}
1292	>	return(qt);
1293	>	}
1294	>	else
1295	>	{
1296	>	va[0] = q1[0]; va[1] = b; va[2] = c;
1297	>	vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1298	>	vc[0] = a; vc[1] = b; vc[2] = q0[2];
1299	>	QT_NTH_CHILD(qt,3) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,3),
1300	>	qt,va,vb,vc,bc,scale>>1,f,n+1,doneptr);
1301	>	if(!*doneptr)
1302	>	{
1303	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1304	>	if(!*doneptr)
1305	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1306	>	if(!*doneptr)
1307	>	f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1308	>	}
1309	>	return(qt);
1310	>	}
1311	>	}
1312		else
803	–	VSUB(db,b[1],b[0]);
804	–	t[0] = HUGET;
805	–	convert_dtol(b[0],bi);
806	–	if(et[1]<0.0 ||(fabs(b[1][0])>(FTINY+1.0))||(fabs(b[1][1])>(FTINY+1.0)) ||
807	–	(fabs(b[1][2])>(FTINY+1.0))||(b[1][0] <-FTINY) ||
808	–	(b[1][1]<-FTINY) ||(b[1][2]<-FTINY))
809	–	{
810	–	max_index(db,&d);
811	–	for(i=0; i< 3; i++)
812	–	dbi[0][i] = (BDIR)(db[i]/d*MAXBDIR);
813	–	}
814	–	else
815	–	for(i=0; i< 3; i++)
816	–	dbi[0][i] = (BDIR)(db[i]*MAXBDIR);
817	–	w=0;
818	–	/* trace the ray starting with this node */
819	–	qt = qtVisit_tri_edges(qt,bi,dbi[0][0],dbi[0][1],dbi[0][2], dbi,&w,
820	–	nextptr,t,0,0,func,f,argptr);
821	–	if(!QT_FLAG_IS_DONE(*f))
1313		{
1314	<	br[0] = (double)bi[0]/(double)MAXBCOORD;
1315	<	br[1] = (double)bi[1]/(double)MAXBCOORD;
825	<	br[2] = (double)bi[2]/(double)MAXBCOORD;
826	<	orig[x] = br[0]*q0[x] + br[1]*q1[x] + br[2]*q2[x];
827	<	orig[y] = br[0]*q0[y] + br[1]*q1[y] + br[2]*q2[y];
828	<	orig[z]=(-FP_N(peq)[x]*orig[x] -
829	<	FP_N(peq)[y]*orig[y]-FP_D(peq))/FP_N(peq)[z];
1314	>	qt = f.func(f.argptr,root,qt,parent,q0,q1,q2,bc,scale,0,f,n,doneptr);
1315	>	return(qt);
1316		}
831	–	return(qt);
832	–
1317		}
1318
1319
1320	+	QUADTREE
1321	+	qtInsert_point_rev(root,qt,parent,q0,q1,q2,bc,scale,f,n,doneptr)
1322	+	int root;
1323	+	QUADTREE qt,parent;
1324	+	BCOORD q0[3],q1[3],q2[3],bc[3];
1325	+	BCOORD scale;
1326	+	FUNC f;
1327	+	int n,*doneptr;
1328	+	{
1329	+	BCOORD a,b,c;
1330	+	BCOORD va[3],vb[3],vc[3];
1331
1332	+	if(QT_IS_TREE(qt))
1333	+	{
1334	+	a = q1[0] - scale;
1335	+	b = q0[1] - scale;
1336	+	c = q0[2] - scale;
1337	+	if(bc[0] < a)
1338	+	{
1339	+	vc[0] = a; vc[1] = b; vc[2] = q0[2];
1340	+	vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1341	+	QT_NTH_CHILD(qt,0) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,0),
1342	+	qt,q0,vc,vb,bc,scale>>1,f,n+1,doneptr);
1343	+	if(!*doneptr)
1344	+	{
1345	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1346	+	if(!*doneptr)
1347	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1348	+	if(!*doneptr)
1349	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1350	+	}
1351	+	return(qt);
1352	+	}
1353	+	if(bc[1] < b)
1354	+	{
1355	+	vc[0] = a; vc[1] = b; vc[2] = q0[2];
1356	+	va[0] = q1[0]; va[1] = b; va[2] = c;
1357	+	QT_NTH_CHILD(qt,1) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,1),
1358	+	qt,vc,q1,va,bc,scale>>1,f,n+1,doneptr);
1359	+	if(!*doneptr)
1360	+	{
1361	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1362	+	if(!*doneptr)
1363	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1364	+	if(!*doneptr)
1365	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1366	+	}
1367	+	return(qt);
1368	+	}
1369	+	if(bc[2] < c)
1370	+	{
1371	+	va[0] = q1[0]; va[1] = b; va[2] = c;
1372	+	vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1373	+	QT_NTH_CHILD(qt,2) = qtInsert_point_rev(root,QT_NTH_CHILD(qt,2),
1374	+	qt,vb,va,q2,bc,scale>>1,f,n+1,doneptr);
1375	+	if(!*doneptr)
1376	+	{
1377	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1378	+	if(!*doneptr)
1379	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1380	+	if(!*doneptr)
1381	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,3),f,doneptr);
1382	+	}
1383	+	return(qt);
1384	+	}
1385	+	else
1386	+	{
1387	+	va[0] = q1[0]; va[1] = b; va[2] = c;
1388	+	vb[0] = a; vb[1] = q0[1]; vb[2] = c;
1389	+	vc[0] = a; vc[1] = b; vc[2] = q0[2];
1390	+	QT_NTH_CHILD(qt,3) = qtInsert_point(root,QT_NTH_CHILD(qt,3),
1391	+	qt,va,vb,vc,bc,scale>>1,f,n+1,doneptr);
1392	+	if(!*doneptr)
1393	+	{
1394	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,0),f,doneptr);
1395	+	if(!*doneptr)
1396	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,1),f,doneptr);
1397	+	if(!*doneptr)
1398	+	f.func_after(f.argptr,QT_NTH_CHILD(qt,2),f,doneptr);
1399	+	}
1400	+	return(qt);
1401	+	}
1402	+	}
1403	+	else
1404	+	{
1405	+	qt = f.func(f.argptr,root,qt,parent,q0,q1,q2,bc,scale,1,f,n,doneptr);
1406	+	return(qt);
1407	+	}
1408	+	}
1409
1410
1411

Diff Legend

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