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.9 by gwlarson, Mon Dec 28 18:07:36 1998 UTC

#	Line 110 | Line 110 | qtDone()                       /* free EVERYTHING */
110		}
111
112		QUADTREE
113	<	qtLocate_leaf(qt,bcoord)
113	>	qtLocate(qt,bcoord)
114		QUADTREE qt;
115		BCOORD bcoord[3];
116		{
#	Line 118 | Line 118 | BCOORD bcoord[3];
118
119		if(QT_IS_TREE(qt))
120		{
121	<	i = baryi_child(bcoord);
121	>	i = bary_child(bcoord);
122
123	<	return(qtLocate_leaf(QT_NTH_CHILD(qt,i),bcoord));
123	>	return(qtLocate(QT_NTH_CHILD(qt,i),bcoord));
124		}
125		else
126		return(qt);
127		}
128
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	–
129		int
130	<	move_to_nbri(b,db0,db1,db2,tptr)
130	>	move_to_nbr(b,db0,db1,db2,tptr)
131		BCOORD b[3];
132		BDIR db0,db1,db2;
133	<	TINT *tptr;
133	>	double *tptr;
134		{
135	<	TINT t,dt;
135	>	double t,dt;
136		BCOORD bc;
137		int nbr;
138
139		nbr = -1;
140	<	*tptr = 0;
140	>	*tptr = 0.0;
141		/* Advance to next node */
142		if(b[0]==0 && db0 < 0)
143		return(0);
#	Line 480 | Line 145 | TINT *tptr;
145		return(1);
146		if(b[2]==0 && db2 < 0)
147		return(2);
483	–
148		if(db0 < 0)
149		{
150	<	bc = b[0]<<SHIFT_MAXBCOORD;
487	<	t = bc/-db0;
150	>	t = ((double)b[0])/-db0;
151		nbr = 0;
152		}
153		else
154	<	t = HUGET;
154	>	t = MAXFLOAT;
155		if(db1 < 0)
156		{
157	<	bc = b[1] <<SHIFT_MAXBCOORD;
495	<	dt = bc/-db1;
157	>	dt = ((double)b[1])/-db1;
158		if( dt < t)
159		{
160		t = dt;
#	Line 501 | Line 163 | TINT *tptr;
163		}
164		if(db2 < 0)
165		{
166	<	bc = b[2] << SHIFT_MAXBCOORD;
505	<	dt = bc/-db2;
166	>	dt = ((double)b[2])/-db2;
167		if( dt < t)
168		{
169		t = dt;
#	Line 513 | Line 174 | TINT *tptr;
174		return(nbr);
175		}
176
177	<	QUADTREE
517	<	qtVisit_tri_edges(qt,b,db0,db1,db2,db,wptr,nextptr,t,sign,sfactor,func,f,argptr)
177	>	qtTrace_ray(qt,b,db0,db1,db2,nextptr,sign,sfactor,func,f)
178		QUADTREE qt;
179		BCOORD b[3];
180	<	BDIR db0,db1,db2,db[3][3];
181	<	int *wptr,*nextptr;
522	<	TINT t[3];
180	>	BDIR db0,db1,db2;
181	>	int *nextptr;
182		int sign,sfactor;
183	<	int (*func)();
184	<	int *f,*argptr;
183	>	FUNC func;
184	>	int *f;
185		{
186		int i,found;
187		QUADTREE child;
188		int nbr,next,w;
189	<	TINT t_g,t_l,t_i,l;
189	>	double t;
190
191		if(QT_IS_TREE(qt))
192		{
193		/* Find the appropriate child and reset the coord */
194	<	i = baryi_child(b);
194	>	i = bary_child(b);
195
537	–	QT_SET_FLAG(qt);
538	–
196		for(;;)
197		{
541	–	w = *wptr;
198		child = QT_NTH_CHILD(qt,i);
199		if(i != 3)
200	<	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);
200	>	qtTrace_ray(child,b,db0,db1,db2,nextptr,sign,sfactor+1,func,f);
201		else
202		/* If the center cell- must flip direction signs */
203	<	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);
203	>	qtTrace_ray(child,b,-db0,-db1,-db2,nextptr,1-sign,sfactor+1,func,f);
204
205		if(QT_FLAG_IS_DONE(*f))
206	<	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	<	}
206	>	return;
207		/* If in same block: traverse */
208		if(i==3)
209		next = *nextptr;
#	Line 568 | Line 213 | qtVisit_tri_edges(qt,b,db0,db1,db2,db,wptr,nextptr,t,s
213		else
214		{
215		/* reset the barycentric coordinates in the parents*/
216	<	baryi_parent(b,i);
216	>	bary_parent(b,i);
217		/* Else pop up to parent and traverse from there */
218		return(qt);
219		}
220	<	baryi_from_child(b,i,next);
220	>	bary_from_child(b,i,next);
221		i = next;
222		}
223		}
#	Line 582 | Line 227 | qtVisit_tri_edges(qt,b,db0,db1,db2,db,wptr,nextptr,t,s
227		if(Pick_cnt < 500)
228		Pick_q[Pick_cnt++]=qt;
229		#endif;
230	<	func(&qt,f,argptr);
230	>	F_FUNC(func)(qt,F_ARGS(func),f);
231		if(QT_FLAG_IS_DONE(*f))
232	+	return;
233	+	/* Advance to next node */
234	+	nbr = move_to_nbr(b,db0,db1,db2,&t);
235	+
236	+	if(nbr==-1)
237		{
238	<	if(!QT_IS_EMPTY(qt))
239	<	QT_LEAF_SET_FLAG(qt);
590	<	return(qt);
238	>	QT_FLAG_SET_DONE(*f);
239	>	return;
240		}
241	+	b[0] += (BCOORD)(t *db0);
242	+	b[1] += (BCOORD)(t *db1);
243	+	b[2] += (BCOORD)(t *db2);
244	+	*nextptr = nbr;
245	+	return;
246
247	<	if(!QT_IS_EMPTY(qt))
248	<	QT_LEAF_SET_FLAG(qt);
595	<	/* Advance to next node */
596	<	w = *wptr;
597	<	while(1)
598	<	{
599	<	nbr = move_to_nbri(b,db0,db1,db2,&t_i);
600	<
601	<	t_g = t_i >> sfactor;
602	<
603	<	if(t_g >= t[w])
604	<	{
605	<	if(w == 2)
606	<	{
607	<	QT_FLAG_SET_DONE(*f);
608	<	return(qt);
609	<	}
610	<	/* The edge fits in the cell- therefore the result of shifting
611	<	db by sfactor is guaranteed to be less than MAXBCOORD
612	<	*/
613	<	/* Caution: (t[w]*db) must occur before divide by MAXBCOORD
614	<	since t[w] will always be < MAXBCOORD
615	<	*/
616	<	l = t[w] << sfactor;
617	<	/* NOTE: Change divides to Shift and multiply by sign*/
618	<	b[0] += (l*db0)/MAXBCOORD;
619	<	b[1] += (l*db1)/MAXBCOORD;
620	<	b[2] += (l*db2)/MAXBCOORD;
621	<	w++;
622	<	db0 = db[w][0]; db1 = db[w][1]; db2 = db[w][2];
623	<	if(sign)
624	<	{  db0 *= -1;db1 *= -1; db2 *= -1;}
625	<	}
626	<	else
627	<	{
628	<	/* Caution: (t_i*db) must occur before divide by MAXBCOORD
629	<	since t_i will always be < MAXBCOORD*/
630	<	/* NOTE: Change divides to Shift and  by sign*/
631	<	b[0] += (t_i *db0) / MAXBCOORD;
632	<	b[1] += (t_i *db1) / MAXBCOORD;
633	<	b[2] += (t_i *db2) / MAXBCOORD;
634	<
635	<	t[w] -= t_g;
636	<	*wptr = w;
637	<	*nextptr = nbr;
638	<	return(qt);
639	<	}
640	<	}
641	<	}
642	<	}
247	>	}
248	>	}
249
250
251	+
252	+
253	+
254	+
255	+
256	+
257	+	#define TEST_INT(tl,th,d,q0,q1,h,l) \
258	+	tl=d>q0;th=d<q1; if(tl&&th) goto Lfunc_call; \
259	+	if(tl) if(l) goto Lfunc_call; else h=TRUE; \
260	+	if(th) if(h) goto Lfunc_call; else l = TRUE; \
261	+	if(th) if(h) goto Lfunc_call; else l = TRUE;
262	+
263	+	/* Leaf node: Do definitive test */
264		QUADTREE
265	<	qtRoot_visit_tri_edges(qt,q0,q1,q2,peq,tri,i_pt,wptr,nextptr,func,f,argptr)
266	<	QUADTREE qt;
267	<	FVECT q0,q1,q2;
268	<	FPEQ  peq;
269	<	FVECT tri[3],i_pt;
270	<	int *wptr,*nextptr;
271	<	int (*func)();
272	<	int *f,*argptr;
273	<	{
274	<	int x,y,z,w,i,j,first;
275	<	QUADTREE child;
276	<	FVECT c,d,v[3];
277	<	double b[4][3],db[3][3],et[3],exit_pt;
278	<	BCOORD bi[3];
279	<	TINT t[3];
280	<	BDIR dbi[3][3];
265	>	qtLeaf_insert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
266	>	scale, s0,s1,s2,sq0,sq1,sq2,f,n)
267	>	int root;
268	>	QUADTREE qt;
269	>	BCOORD q0[3],q1[3],q2[3];
270	>	BCOORD t0[3],t1[3],t2[3];
271	>	BCOORD dt10[3],dt21[3],dt02[3];
272	>	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
273	>	FUNC f;
274	>	int n;
275	>	{
276	>	double db;
277	>	unsigned int e0,e1,e2;
278	>	char al,ah,bl,bh,cl,ch,testl,testh;
279	>	char test_t0t1,test_t1t2,test_t2t0;
280	>	BCOORD a,b,c;
281	>
282	>	/* First check if can trivial accept: if vertex in cell */
283	>	if(s0 & s1 & s2)
284	>	goto Lfunc_call;
285	>
286	>	/* Assumption: Could not trivial reject*/
287	>	/* IF any coords lie on edge- must be in if couldnt reject */
288	>	a = q1[0];b= q0[1]; c= q0[2];
289	>	if(t0[0]== a || t0[1] == b || t0[2] == c)
290	>	goto Lfunc_call;
291	>	if(t1[0]== a || t1[1] == b || t1[2] == c)
292	>	goto Lfunc_call;
293	>	if(t2[0]== a || t2[1] == b || t2[2] == c)
294	>	goto Lfunc_call;
295
296	<	first =0;
297	<	w = *wptr;
298	<	if(w==-1)
296	>	/* Test for edge crossings */
297	>	/* Test t0t1,t1t2,t2t0 against a */
298	>	/* Calculate edge crossings */
299	>	e0  = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
300	>	/* See if edge can be trivially rejected from intersetion testing */
301	>	test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
302	>	((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
303	>	bl=bh=0;
304	>	if(test_t0t1 && (e0 & 2) )
305		{
306	<	first = 1;
307	<	*wptr = w = 0;
306	>	/* Must do double calculation to avoid overflow */
307	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
308	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
309		}
310	<	/* Project the origin onto the root node plane */
310	>	test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
311	>	((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
312	>	if(test_t1t2 && (e0 & 1))
313	>	{    /* test t1t2 against a */
314	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
315	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
316	>	}
317	>	test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
318	>	((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
319	>	if(test_t2t0 && (e0 & 4))
320	>	{
321	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
322	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
323	>	}
324	>	e1  = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
325	>	cl = ch = 0;
326	>	if(test_t0t1 && (e1 & 2))
327	>	{/* test t0t1 against b */
328	>	db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
329	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
330	>	}
331	>	if(test_t1t2 && (e1 & 1))
332	>	{/* test t1t2 against b */
333	>	db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
334	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
335	>	}
336	>	if(test_t2t0 && (e1 & 4))
337	>	{/* test t2t0 against b */
338	>	db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
339	>	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
340	>	}
341	>
342	>	/* test edges against c */
343	>	e2  = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
344	>	al = ah = 0;
345	>	if(test_t0t1 && (e2 & 2))
346	>	{ /* test t0t1 against c */
347	>	db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
348	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
349	>	}
350	>	if(test_t1t2 && (e2 & 1))
351	>	{
352	>	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
353	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
354	>	}
355	>	if(test_t2t0 && (e2 & 4))
356	>	{ /* test t2t0 against c */
357	>	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
358	>	TEST_INT(testl,testh,db,a,q0[0],al,ah)
359	>	}
360	>	/* Only remaining case is if some edge trivially rejected */
361	>	if(!e0 || !e1 || !e2)
362	>	return(qt);
363
364	<	t[0] = t[1] = t[2] = 0;
365	<	/* Find the intersection point of the origin */
366	<	/* map to 2d by dropping maximum magnitude component of normal */
364	>	/* Only one/none got tested - pick either of other two/any two */
365	>	/* Only need to test one edge */
366	>	if(!test_t0t1 && (e0 & 2))
367	>	{
368	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
369	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
370	>	}
371	>	if(!test_t1t2 && (e0 & 1))
372	>	{/* test t1t2 against a */
373	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
374	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
375	>	}
376	>	if(!test_t2t0 && (e0 & 4))
377	>	{
378	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
379	>	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
380	>	}
381
382	<	x = FP_X(peq);
383	<	y = FP_Y(peq);
384	<	z = FP_Z(peq);
385	<	/* Calculate barycentric coordinates for current vertex */
386	<	if(!first)
387	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[3]);
388	<	else
389	<	/* Just starting: b[0] is the origin point: guaranteed to be valid b*/
382	>	return(qt);
383	>
384	>	Lfunc_call:
385	>	qt = f.func(f.argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
386	>	s0,s1,s2,sq0,sq1,sq2,0,f,n);
387	>	return(qt);
388	>
389	>	}
390	>
391	>
392	>
393	>	/* Leaf node: Do definitive test */
394	>	QUADTREE
395	>	qtLeaf_insert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
396	>	scale, s0,s1,s2,sq0,sq1,sq2,f,n)
397	>	int root;
398	>	QUADTREE qt;
399	>	BCOORD q0[3],q1[3],q2[3];
400	>	BCOORD t0[3],t1[3],t2[3];
401	>	BCOORD dt10[3],dt21[3],dt02[3];
402	>	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
403	>	FUNC f;
404	>	int n;
405	>	{
406	>	double db;
407	>	unsigned int e0,e1,e2;
408	>	BCOORD a,b,c;
409	>	double p0[2], p1[2],cp;
410	>	char al,ah,bl,bh,cl,ch;
411	>	char testl,testh,test_t0t1,test_t1t2,test_t2t0;
412	>	unsigned int ls0,ls1,ls2;
413	>
414	>	/* May have gotten passed trivial reject if one/two vertices are ON */
415	>	a = q1[0];b= q0[1]; c= q0[2];
416	>	SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
417	>
418	>	/* First check if can trivial accept: if vertex in cell */
419	>	if(ls0 & ls1 & ls2)
420	>	goto Lfunc_call;
421	>
422	>	if(ls0==0 || ls1 == 0 || ls2 ==0)
423	>	return(qt);
424	>	/* Assumption: Could not trivial reject*/
425	>	/* IF any coords lie on edge- must be in if couldnt reject */
426	>
427	>	if(t0[0]== a || t0[1] == b || t0[2] == c)
428	>	goto Lfunc_call;
429	>	if(t1[0]== a || t1[1] == b || t1[2] == c)
430	>	goto Lfunc_call;
431	>	if(t2[0]== a || t2[1] == b || t2[2] == c)
432	>	goto Lfunc_call;
433	>
434	>	/* Test for edge crossings */
435	>	/* Test t0t1 against a,b,c */
436	>	/* First test if t0t1 can be trivially rejected */
437	>	/* If both edges are outside bounds- dont need to test */
438	>
439	>	/* Test t0t1,t1t2,t2t0 against a */
440	>	test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
441	>	((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
442	>	e0  = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
443	>	bl=bh=0;
444	>	/* Test t0t1,t1t2,t2t0 against a */
445	>	if(test_t0t1 && (e0 & 2) )
446		{
447	<	intersect_vector_plane(tri[0],peq,&(et[0]),v[0]);
448	<	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]);
447	>	db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
448	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
449		}
450	+	test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
451	+	((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
452	+	if(test_t1t2 && (e0 & 1))
453	+	{    /* test t1t2 against a */
454	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
455	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
456	+	}
457	+	test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
458	+	((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
459	+	if(test_t2t0 && (e0 & 4))
460	+	{
461	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
462	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
463	+	}
464	+	cl = ch = 0;
465	+	e1  = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
466	+	if(test_t0t1 && (e1 & 2))
467	+	{/* test t0t1 against b */
468	+	db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
469	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
470	+	}
471	+	if(test_t1t2 && (e1 & 1))
472	+	{/* test t1t2 against b */
473	+	db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
474	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
475	+	}
476	+	if(test_t2t0 && (e1 & 4))
477	+	{/* test t2t0 against b */
478	+	db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
479	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
480	+	}
481	+	al = ah = 0;
482	+	e2  = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
483	+	if(test_t0t1 && (e2 & 2))
484	+	{ /* test t0t1 against c */
485	+	db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
486	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
487	+	}
488	+	if(test_t1t2 && (e2 & 1))
489	+	{
490	+	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
491	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
492	+	}
493	+	if(test_t2t0 && (e2 & 4))
494	+	{ /* test t2t0 against c */
495	+	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
496	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
497	+	}
498	+	/* Only remaining case is if some edge trivially rejected */
499	+	if(!e0 || !e1 || !e2)
500	+	return(qt);
501
502	<	j = (w+1)%3;
503	<	intersect_vector_plane(tri[j],peq,&(et[j]),v[j]);
504	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[j][x],v[j][y],b[j]);
693	<	if(et[j] < 0.0)
502	>	/* Only one/none got tested - pick either of other two/any two */
503	>	/* Only need to test one edge */
504	>	if(!test_t0t1 && (e0 & 2))
505		{
506	<	VSUB(db[w],b[3],b[j]);
507	<	t[w] = HUGET;
506	>	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
507	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
508		}
509	<	else
510	<	{
511	<	/* NOTE: for stability: do not increment with ipt- use full dir and
512	<	calculate t: but for wrap around case: could get same problem?
702	<	*/
703	<	VSUB(db[w],b[j],b[3]);
704	<	/* Check if the point is out side of the triangle: if so t[w] =HUGET */
705	<	if((fabs(b[j][0])>(FTINY+1.0)) ||(fabs(b[j][1])>(FTINY+1.0)) ||
706	<	(fabs(b[j][2])>(FTINY+1.0))||(b[j][0] <-FTINY) ||
707	<	(b[j][1]<-FTINY) ||(b[j][2]<-FTINY))
708	<	t[w] = HUGET;
709	<	else
710	<	{
711	<	/* The next point is in the triangle- so db values will be in valid
712	<	range and t[w]= MAXT
713	<	*/
714	<	t[w] = MAXT;
715	<	if(j != 0)
716	<	for(;j < 3;j++)
717	<	{
718	<	i = (j+1)%3;
719	<	if(!first || i != 0)
720	<	{
721	<	intersect_vector_plane(tri[i],peq,&(et[i]),v[i]);
722	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],v[i][x],
723	<	v[i][y],b[i]);
724	<	}
725	<	if(et[i] < 0.0)
726	<	{
727	<	VSUB(db[j],b[j],b[i]);
728	<	t[j] = HUGET;
729	<	break;
730	<	}
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	<	}
509	>	if(!test_t1t2 && (e0 & 1))
510	>	{/* test t1t2 against a */
511	>	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
512	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
513		}
514	<	}
747	<	bary_dtol(b[3],db,bi,dbi,t,w);
748	<
749	<	/* trace the ray starting with this node */
750	<	qt = qtVisit_tri_edges(qt,bi,dbi[w][0],dbi[w][1],dbi[w][2],
751	<	dbi,wptr,nextptr,t,0,0,func,f,argptr);
752	<	if(!QT_FLAG_IS_DONE(*f))
514	>	if(!test_t2t0 && (e0 & 4))
515		{
516	<	b[3][0] = (double)bi[0]/(double)MAXBCOORD;
517	<	b[3][1] = (double)bi[1]/(double)MAXBCOORD;
756	<	b[3][2] = (double)bi[2]/(double)MAXBCOORD;
757	<	i_pt[x] = b[3][0]*q0[x] + b[3][1]*q1[x] + b[3][2]*q2[x];
758	<	i_pt[y] = b[3][0]*q0[y] + b[3][1]*q1[y] + b[3][2]*q2[y];
759	<	i_pt[z] = (-FP_N(peq)[x]*i_pt[x] - FP_N(peq)[y]*i_pt[y]-FP_D(peq))/FP_N(peq)[z];
516	>	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
517	>	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
518		}
519	+
520		return(qt);
521	<
521	>	Lfunc_call:
522	>
523	>	qt = f.func(f.argptr,root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
524	>	s0,s1,s2,sq0,sq1,sq2,1,f,n);
525	>	return(qt);
526	>	}
527	>
528	>
529	>
530	>	/* ASSUMPTION: that triangle t0,t1,t2 intersects quad cell q0,q1,q2 */
531	>
532	>	/*-------q2--------- sq2
533	>	/ \
534	>	s1     /sc \ s0
535	>	qb_____qa
536	>	/ \   / \
537	>	\sq0/sa \ /sb \   /sq1
538	>	\ q0____qc____q1/
539	>	\             /
540	>	\     s2    /
541	>	*/
542	>
543	>	int Find_id=0;
544	>
545	>	QUADTREE
546	>	qtInsert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
547	>	s0,s1,s2,sq0,sq1,sq2,f,n)
548	>	int root;
549	>	QUADTREE qt;
550	>	BCOORD q0[3],q1[3],q2[3];
551	>	BCOORD t0[3],t1[3],t2[3];
552	>	BCOORD dt10[3],dt21[3],dt02[3];
553	>	BCOORD scale;
554	>	unsigned int s0,s1,s2,sq0,sq1,sq2;
555	>	FUNC f;
556	>	int n;
557	>	{
558	>	BCOORD a,b,c;
559	>	BCOORD va[3],vb[3],vc[3];
560	>	unsigned int sa,sb,sc;
561	>
562	>	/* If a tree: trivial reject and recurse on potential children */
563	>	if(QT_IS_TREE(qt))
564	>	{
565	>	/* Test against new a edge: if entirely to left: only need
566	>	to visit child 0
567	>	*/
568	>	a = q1[0] + scale;
569	>	b = q0[1] + scale;
570	>	c = q0[2] + scale;
571	>
572	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
573	>
574	>	if(sa==7)
575	>	{
576	>	vb[1] = q0[1];
577	>	vc[2] = q0[2];
578	>	vc[1] = b;
579	>	vb[2] = c;
580	>	vb[0] = vc[0] = a;
581	>	QT_NTH_CHILD(qt,0) = qtInsert_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
582	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
583	>	return(qt);
584	>	}
585	>	if(sb==7)
586	>	{
587	>	va[0] = q1[0];
588	>	vc[2] = q0[2];
589	>	va[1] = vc[1] = b;
590	>	va[2] = c;
591	>	vc[0] = a;
592	>	QT_NTH_CHILD(qt,1) = qtInsert_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,
593	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
594	>	return(qt);
595	>	}
596	>	if(sc==7)
597	>	{
598	>	va[0] = q1[0];
599	>	vb[1] = q0[1];
600	>	va[1] = b;
601	>	va[2] = vb[2] = c;
602	>	vb[0] = a;
603	>	QT_NTH_CHILD(qt,2) = qtInsert_tri(root,QT_NTH_CHILD(qt,2),vb,va,
604	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
605	>	return(qt);
606	>	}
607	>
608	>	va[0] = q1[0];
609	>	vb[1] = q0[1];
610	>	vc[2] = q0[2];
611	>	va[1] = vc[1] = b;
612	>	va[2] = vb[2] = c;
613	>	vb[0] = vc[0] = a;
614	>	/* If outside of all edges: only need to Visit 3  */
615	>	if(sa==0 && sb==0 && sc==0)
616	>	{
617	>	QT_NTH_CHILD(qt,3) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,
618	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
619	>	return(qt);
620	>	}
621	>
622	>	if(sa)
623	>	QT_NTH_CHILD(qt,0) = qtInsert_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
624	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
625	>	if(sb)
626	>	QT_NTH_CHILD(qt,1) = qtInsert_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
627	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
628	>	if(sc)
629	>	QT_NTH_CHILD(qt,2) = qtInsert_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
630	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
631	>	QT_NTH_CHILD(qt,3) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
632	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
633	>	return(qt);
634	>	}
635	>	/* Leaf node: Do definitive test */
636	>	else
637	>	return(qt = qtLeaf_insert_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
638	>	scale,s0,s1,s2,sq0,sq1,sq2,f,n));
639		}
640
641
642		QUADTREE
643	<	qtRoot_trace_ray(qt,q0,q1,q2,peq,orig,dir,nextptr,func,f,argptr)
644	<	QUADTREE qt;
645	<	FVECT q0,q1,q2;
646	<	FPEQ  peq;
647	<	FVECT orig,dir;
648	<	int *nextptr;
649	<	int (*func)();
650	<	int *f,*argptr;
643	>	qtInsert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
644	>	s0,s1,s2,sq0,sq1,sq2,f,n)
645	>	int root;
646	>	QUADTREE qt;
647	>	BCOORD q0[3],q1[3],q2[3];
648	>	BCOORD t0[3],t1[3],t2[3];
649	>	BCOORD dt10[3],dt21[3],dt02[3];
650	>	BCOORD scale;
651	>	unsigned int s0,s1,s2,sq0,sq1,sq2;
652	>	FUNC f;
653		{
654	<	int x,y,z,nbr,w,i;
655	<	QUADTREE child;
656	<	FVECT v,i_pt;
779	<	double b[2][3],db[3],et[2],d,br[3];
780	<	BCOORD bi[3];
781	<	TINT t[3];
782	<	BDIR dbi[3][3];
783	<
784	<	/* Project the origin onto the root node plane */
785	<	t[0] = t[1] = t[2] = 0;
654	>	BCOORD a,b,c;
655	>	BCOORD va[3],vb[3],vc[3];
656	>	unsigned int sa,sb,sc;
657
658	<	VADD(v,orig,dir);
659	<	/* Find the intersection point of the origin */
660	<	/* map to 2d by dropping maximum magnitude component of normal */
661	<	x = FP_X(peq);
662	<	y = FP_Y(peq);
663	<	z = FP_Z(peq);
658	>	/* If a tree: trivial reject and recurse on potential children */
659	>	if(QT_IS_TREE(qt))
660	>	{
661	>	/* Test against new a edge: if entirely to left: only need
662	>	to visit child 0
663	>	*/
664	>	a = q1[0] - scale;
665	>	b = q0[1] - scale;
666	>	c = q0[2] - scale;
667
668	<	/* Calculate barycentric coordinates for current vertex */
795	<	intersect_vector_plane(orig,peq,&(et[0]),i_pt);
796	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[0]);
668	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
669
670	<	intersect_vector_plane(v,peq,&(et[1]),i_pt);
671	<	bary2d(q0[x],q0[y],q1[x],q1[y],q2[x],q2[y],i_pt[x],i_pt[y],b[1]);
672	<	if(et[1] < 0.0)
673	<	VSUB(db,b[0],b[1]);
670	>	if(sa==0)
671	>	{
672	>	vb[1] = q0[1];
673	>	vc[2] = q0[2];
674	>	vc[1] = b;
675	>	vb[2] = c;
676	>	vb[0] = vc[0] = a;
677	>
678	>	QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
679	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
680	>	return(qt);
681	>	}
682	>	if(sb==0)
683	>	{
684	>	va[0] = q1[0];
685	>	vc[2] = q0[2];
686	>	va[1] = vc[1] = b;
687	>	va[2] = c;
688	>	vc[0] = a;
689	>	QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
690	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
691	>	return(qt);
692	>	}
693	>	if(sc==0)
694	>	{
695	>	va[0] = q1[0];
696	>	vb[1] = q0[1];
697	>	va[1] = b;
698	>	va[2] = vb[2] = c;
699	>	vb[0] = a;
700	>	QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
701	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
702	>	return(qt);
703	>	}
704	>	va[0] = q1[0];
705	>	vb[1] = q0[1];
706	>	vc[2] = q0[2];
707	>	va[1] = vc[1] = b;
708	>	va[2] = vb[2] = c;
709	>	vb[0] = vc[0] = a;
710	>	/* If outside of all edges: only need to Visit 3  */
711	>	if(sa==7 && sb==7 && sc==7)
712	>	{
713	>	QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,
714	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
715	>	return(qt);
716	>	}
717	>	if(sa!=7)
718	>	QT_NTH_CHILD(qt,0) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
719	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f,n+1);
720	>	if(sb!=7)
721	>	QT_NTH_CHILD(qt,1) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
722	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f,n+1);
723	>	if(sc!=7)
724	>	QT_NTH_CHILD(qt,2) = qtInsert_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
725	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f,n+1);
726	>
727	>	QT_NTH_CHILD(qt,3) = qtInsert_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
728	>	t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f,n+1);
729	>	return(qt);
730	>	}
731	>	/* Leaf node: Do definitive test */
732		else
733	<	VSUB(db,b[1],b[0]);
734	<	t[0] = HUGET;
735	<	convert_dtol(b[0],bi);
736	<	if(et[1]<0.0 ||(fabs(b[1][0])>(FTINY+1.0))||(fabs(b[1][1])>(FTINY+1.0)) ||
737	<	(fabs(b[1][2])>(FTINY+1.0))||(b[1][0] <-FTINY) ||
738	<	(b[1][1]<-FTINY) ||(b[1][2]<-FTINY))
739	<	{
740	<	max_index(db,&d);
741	<	for(i=0; i< 3; i++)
742	<	dbi[0][i] = (BDIR)(db[i]/d*MAXBDIR);
743	<	}
744	<	else
745	<	for(i=0; i< 3; i++)
746	<	dbi[0][i] = (BDIR)(db[i]*MAXBDIR);
747	<	w=0;
748	<	/* trace the ray starting with this node */
749	<	qt = qtVisit_tri_edges(qt,bi,dbi[0][0],dbi[0][1],dbi[0][2], dbi,&w,
750	<	nextptr,t,0,0,func,f,argptr);
751	<	if(!QT_FLAG_IS_DONE(*f))
733	>	return(qt = qtLeaf_insert_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
734	>	scale,s0,s1,s2,sq0,sq1,sq2,f,n));
735	>	}
736	>
737	>
738	>
739	>
740	>	/*************************************************************************/
741	>	/* Visit code for applying functions: NOTE THIS IS THE SAME AS INSERT CODE:
742	>	except sets flags: wanted insert to be as efficient as possible: so
743	>	broke into two sets of routines
744	>	*/
745	>
746	>	qtVisit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
747	>	s0,s1,s2,sq0,sq1,sq2,f)
748	>	int root;
749	>	QUADTREE qt;
750	>	BCOORD q0[3],q1[3],q2[3];
751	>	BCOORD t0[3],t1[3],t2[3];
752	>	BCOORD dt10[3],dt21[3],dt02[3];
753	>	BCOORD scale;
754	>	unsigned int s0,s1,s2,sq0,sq1,sq2;
755	>	FUNC f;
756	>	{
757	>	BCOORD a,b,c;
758	>	BCOORD va[3],vb[3],vc[3];
759	>	unsigned int sa,sb,sc;
760	>
761	>	/* If a tree: trivial reject and recurse on potential children */
762	>	if(QT_IS_TREE(qt))
763		{
764	<	br[0] = (double)bi[0]/(double)MAXBCOORD;
765	<	br[1] = (double)bi[1]/(double)MAXBCOORD;
766	<	br[2] = (double)bi[2]/(double)MAXBCOORD;
767	<	orig[x] = br[0]*q0[x] + br[1]*q1[x] + br[2]*q2[x];
768	<	orig[y] = br[0]*q0[y] + br[1]*q1[y] + br[2]*q2[y];
769	<	orig[z]=(-FP_N(peq)[x]*orig[x] -
770	<	FP_N(peq)[y]*orig[y]-FP_D(peq))/FP_N(peq)[z];
764	>	QT_SET_FLAG(qt);
765	>
766	>	/* Test against new a edge: if entirely to left: only need
767	>	to visit child 0
768	>	*/
769	>	a = q1[0] + scale;
770	>	b = q0[1] + scale;
771	>	c = q0[2] + scale;
772	>
773	>	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
774	>
775	>	if(sa==7)
776	>	{
777	>	vb[1] = q0[1];
778	>	vc[2] = q0[2];
779	>	vc[1] = b;
780	>	vb[2] = c;
781	>	vb[0] = vc[0] = a;
782	>	qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,
783	>	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f);
784	>	return;
785	>	}
786	>	if(sb==7)
787	>	{
788	>	va[0] = q1[0];
789	>	vc[2] = q0[2];
790	>	va[1] = vc[1] = b;
791	>	va[2] = c;
792	>	vc[0] = a;
793	>	qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,
794	>	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f);
795	>	return;
796	>	}
797	>	if(sc==7)
798	>	{
799	>	va[0] = q1[0];
800	>	vb[1] = q0[1];
801	>	va[1] = b;
802	>	va[2] = vb[2] = c;
803	>	vb[0] = a;
804	>	qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,
805	>	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f);
806	>	return;
807	>	}
808	>
809	>	va[0] = q1[0];
810	>	vb[1] = q0[1];
811	>	vc[2] = q0[2];
812	>	va[1] = vc[1] = b;
813	>	va[2] = vb[2] = c;
814	>	vb[0] = vc[0] = a;
815	>	/* If outside of all edges: only need to Visit 3  */
816	>	if(sa==0 && sb==0 && sc==0)
817	>	{
818	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,
819	>	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f);
820	>	return;
821	>	}
822	>
823	>	if(sa)
824	>	qtVisit_tri(root,QT_NTH_CHILD(qt,0),q0,vc,vb,t0,
825	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f);
826	>	if(sb)
827	>	qtVisit_tri(root,QT_NTH_CHILD(qt,1),vc,q1,va,t0,
828	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f);
829	>	if(sc)
830	>	qtVisit_tri(root,QT_NTH_CHILD(qt,2),vb,va,q2,t0,
831	>	t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f);
832	>	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,
833	>	t1,t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f);
834		}
835	<	return(qt);
836	<
835	>	/* Leaf node: Do definitive test */
836	>	else
837	>	qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
838	>	scale,s0,s1,s2,sq0,sq1,sq2,f);
839	>
840		}
841
842
843	+	qtVisit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,scale,
844	+	s0,s1,s2,sq0,sq1,sq2,f)
845	+	int root;
846	+	QUADTREE qt;
847	+	BCOORD q0[3],q1[3],q2[3];
848	+	BCOORD t0[3],t1[3],t2[3];
849	+	BCOORD dt10[3],dt21[3],dt02[3];
850	+	BCOORD scale;
851	+	unsigned int s0,s1,s2,sq0,sq1,sq2;
852	+	FUNC f;
853	+	{
854	+	BCOORD a,b,c;
855	+	BCOORD va[3],vb[3],vc[3];
856	+	unsigned int sa,sb,sc;
857
858	+	/* If a tree: trivial reject and recurse on potential children */
859	+	if(QT_IS_TREE(qt))
860	+	{
861	+	QT_SET_FLAG(qt);
862	+	/* Test against new a edge: if entirely to left: only need
863	+	to visit child 0
864	+	*/
865	+	a = q1[0] - scale;
866	+	b = q0[1] - scale;
867	+	c = q0[2] - scale;
868
869	+	SIDES_GTR(t0,t1,t2,sa,sb,sc,a,b,c);
870
871	+	if(sa==0)
872	+	{
873	+	vb[1] = q0[1];
874	+	vc[2] = q0[2];
875	+	vc[1] = b;
876	+	vb[2] = c;
877	+	vb[0] = vc[0] = a;
878	+	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,
879	+	vb,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,s1,s2,sq0,sb,sc,f);
880	+	return;
881	+	}
882	+	if(sb==0)
883	+	{
884	+	va[0] = q1[0];
885	+	vc[2] = q0[2];
886	+	va[1] = vc[1] = b;
887	+	va[2] = c;
888	+	vc[0] = a;
889	+	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
890	+	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f);
891	+	return;
892	+	}
893	+	if(sc==0)
894	+	{
895	+	va[0] = q1[0];
896	+	vb[1] = q0[1];
897	+	va[1] = b;
898	+	va[2] = vb[2] = c;
899	+	vb[0] = a;
900	+	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,
901	+	q2,t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f);
902	+	return;
903	+	}
904	+	va[0] = q1[0];
905	+	vb[1] = q0[1];
906	+	vc[2] = q0[2];
907	+	va[1] = vc[1] = b;
908	+	va[2] = vb[2] = c;
909	+	vb[0] = vc[0] = a;
910	+	/* If outside of all edges: only need to Visit 3  */
911	+	if(sa==7 && sb==7 && sc==7)
912	+	{
913	+	qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,
914	+	vc,t0,t1,t2,dt10,dt21,dt02, scale>>1,sa,sb,sc,s0,s1,s2,f);
915	+	return;
916	+	}
917	+	if(sa!=7)
918	+	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,0),q0,vc,vb,
919	+	t0,t1,t2,dt10,dt21,dt02,scale>>1,sa,s1,s2,sq0,sb,sc,f);
920	+	if(sb!=7)
921	+	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,1),vc,q1,va,
922	+	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,sb,s2,sa,sq1,sc,f);
923	+	if(sc!=7)
924	+	qtVisit_tri_rev(root,QT_NTH_CHILD(qt,2),vb,va,q2,
925	+	t0,t1,t2,dt10,dt21,dt02,scale>>1,s0,s1,sc,sa,sb,sq2,f);
926
927	+	qtVisit_tri(root,QT_NTH_CHILD(qt,3),va,vb,vc,t0,t1,
928	+	t2,dt10,dt21,dt02,scale>>1,sa,sb,sc,s0,s1,s2,f);
929	+	return;
930	+	}
931	+	/* Leaf node: Do definitive test */
932	+	else
933	+	qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
934	+	scale,s0,s1,s2,sq0,sq1,sq2,f);
935	+	}
936
937
938
939	+	qtLeaf_visit_tri(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
940	+	scale, s0,s1,s2,sq0,sq1,sq2,f)
941	+	int root;
942	+	QUADTREE qt;
943	+	BCOORD q0[3],q1[3],q2[3];
944	+	BCOORD t0[3],t1[3],t2[3];
945	+	BCOORD dt10[3],dt21[3],dt02[3];
946	+	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
947	+	FUNC f;
948	+	{
949	+	double db;
950	+	unsigned int e0,e1,e2;
951	+	char al,ah,bl,bh,cl,ch,testl,testh;
952	+	char test_t0t1,test_t1t2,test_t2t0;
953	+	BCOORD a,b,c;
954
955	+	/* First check if can trivial accept: if vertex in cell */
956	+	if(s0 & s1 & s2)
957	+	goto Lfunc_call;
958	+
959	+	/* Assumption: Could not trivial reject*/
960	+	/* IF any coords lie on edge- must be in if couldnt reject */
961	+	a = q1[0];b= q0[1]; c= q0[2];
962	+	if(t0[0]== a || t0[1] == b || t0[2] == c)
963	+	goto Lfunc_call;
964	+	if(t1[0]== a || t1[1] == b || t1[2] == c)
965	+	goto Lfunc_call;
966	+	if(t2[0]== a || t2[1] == b || t2[2] == c)
967	+	goto Lfunc_call;
968	+
969	+	/* Test for edge crossings */
970	+	/* Test t0t1,t1t2,t2t0 against a */
971	+	/* Calculate edge crossings */
972	+	e0  = (s0 ^ ((s0 >> 1) | (s0 << 2 & 4)));
973	+	/* See if edge can be trivially rejected from intersetion testing */
974	+	test_t0t1 = !(((s0 & 6)==0) || ((s1 & 6)==0)|| ((s2 & 6)==0) ||
975	+	((sq0 & 6)==6) ||((sq1 & 6)==6)|| ((sq2 & 6)==6));
976	+	bl=bh=0;
977	+	if(test_t0t1 && (e0 & 2) )
978	+	{
979	+	/* Must do double calculation to avoid overflow */
980	+	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
981	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
982	+	}
983	+	test_t1t2= !(((s0 & 3)==0) || ((s1 & 3)==0)|| ((s2 & 3)==0) ||
984	+	((sq0 & 3)==3) ||((sq1 & 3)==3)|| ((sq2 & 3)==3));
985	+	if(test_t1t2 && (e0 & 1))
986	+	{    /* test t1t2 against a */
987	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
988	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
989	+	}
990	+	test_t2t0 = !(((s0 & 5)==0) || ((s1 & 5)==0)|| ((s2 & 5)==0) ||
991	+	((sq0 & 5)==5) ||((sq1 & 5)==5)|| ((sq2 & 5)==5));
992	+	if(test_t2t0 && (e0 & 4))
993	+	{
994	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
995	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
996	+	}
997	+	e1  = (s1 ^ ((s1 >> 1) | (s1 << 2 & 4)));
998	+	cl = ch = 0;
999	+	if(test_t0t1 && (e1 & 2))
1000	+	{/* test t0t1 against b */
1001	+	db = t0[2] + dt10[2]*((double)(b-t0[1]))/dt10[1];
1002	+	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1003	+	}
1004	+	if(test_t1t2 && (e1 & 1))
1005	+	{/* test t1t2 against b */
1006	+	db = t1[2] + dt21[2]*((double)(q0[1] - t1[1]))/dt21[1];
1007	+	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1008	+	}
1009	+	if(test_t2t0 && (e1 & 4))
1010	+	{/* test t2t0 against b */
1011	+	db = t2[2] + dt02[2]*((double)(q0[1] - t2[1]))/dt02[1];
1012	+	TEST_INT(testl,testh,db,c,q2[2],cl,ch)
1013	+	}
1014	+
1015	+	/* test edges against c */
1016	+	e2  = (s2 ^ ((s2 >> 1) | (s2 << 2 & 4)));
1017	+	al = ah = 0;
1018	+	if(test_t0t1 && (e2 & 2))
1019	+	{ /* test t0t1 against c */
1020	+	db = t0[0] + dt10[0]*((double)(c-t0[2]))/dt10[2];
1021	+	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1022	+	}
1023	+	if(test_t1t2 && (e2 & 1))
1024	+	{
1025	+	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1026	+	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1027	+	}
1028	+	if(test_t2t0 && (e2 & 4))
1029	+	{ /* test t2t0 against c */
1030	+	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1031	+	TEST_INT(testl,testh,db,a,q0[0],al,ah)
1032	+	}
1033	+	/* Only remaining case is if some edge trivially rejected */
1034	+	if(!e0 || !e1 || !e2)
1035	+	return;
1036	+
1037	+	/* Only one/none got tested - pick either of other two/any two */
1038	+	/* Only need to test one edge */
1039	+	if(!test_t0t1 && (e0 & 2))
1040	+	{
1041	+	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1042	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1043	+	}
1044	+	if(!test_t1t2 && (e0 & 1))
1045	+	{/* test t1t2 against a */
1046	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1047	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1048	+	}
1049	+	if(!test_t2t0 && (e0 & 4))
1050	+	{
1051	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1052	+	TEST_INT(testl,testh,db,b,q1[1],bl,bh)
1053	+	}
1054	+
1055	+	return;
1056	+
1057	+	Lfunc_call:
1058	+	f.func(f.argptr,root,qt);
1059	+	if(!QT_IS_EMPTY(qt))
1060	+	QT_LEAF_SET_FLAG(qt);
1061	+	}
1062	+
1063	+
1064	+
1065	+	/* Leaf node: Do definitive test */
1066	+	QUADTREE
1067	+	qtLeaf_visit_tri_rev(root,qt,q0,q1,q2,t0,t1,t2,dt10,dt21,dt02,
1068	+	scale, s0,s1,s2,sq0,sq1,sq2,f)
1069	+	int root;
1070	+	QUADTREE qt;
1071	+	BCOORD q0[3],q1[3],q2[3];
1072	+	BCOORD t0[3],t1[3],t2[3];
1073	+	BCOORD dt10[3],dt21[3],dt02[3];
1074	+	unsigned int scale,s0,s1,s2,sq0,sq1,sq2;
1075	+	FUNC f;
1076	+	{
1077	+	double db;
1078	+	unsigned int e0,e1,e2;
1079	+	BCOORD a,b,c;
1080	+	double p0[2], p1[2],cp;
1081	+	char al,ah,bl,bh,cl,ch;
1082	+	char testl,testh,test_t0t1,test_t1t2,test_t2t0;
1083	+	unsigned int ls0,ls1,ls2;
1084	+
1085	+	/* May have gotten passed trivial reject if one/two vertices are ON */
1086	+	a = q1[0];b= q0[1]; c= q0[2];
1087	+	SIDES_LESS(t0,t1,t2,ls0,ls1,ls2,a,b,c);
1088	+
1089	+	/* First check if can trivial accept: if vertex in cell */
1090	+	if(ls0 & ls1 & ls2)
1091	+	goto Lfunc_call;
1092	+
1093	+	if(ls0==0 || ls1 == 0 || ls2 ==0)
1094	+	return;
1095	+	/* Assumption: Could not trivial reject*/
1096	+	/* IF any coords lie on edge- must be in if couldnt reject */
1097	+
1098	+	if(t0[0]== a || t0[1] == b || t0[2] == c)
1099	+	goto Lfunc_call;
1100	+	if(t1[0]== a || t1[1] == b || t1[2] == c)
1101	+	goto Lfunc_call;
1102	+	if(t2[0]== a || t2[1] == b || t2[2] == c)
1103	+	goto Lfunc_call;
1104	+
1105	+	/* Test for edge crossings */
1106	+	/* Test t0t1 against a,b,c */
1107	+	/* First test if t0t1 can be trivially rejected */
1108	+	/* If both edges are outside bounds- dont need to test */
1109	+
1110	+	/* Test t0t1,t1t2,t2t0 against a */
1111	+	test_t0t1 = !(((ls0 & 6)==0) || ((ls1 & 6)==0)|| ((ls2 & 6)==0) ||
1112	+	((sq0 & 6)==0) ||((sq1 & 6)==0)|| ((sq2 & 6)==0));
1113	+	e0  = (ls0 ^ ((ls0 >> 1) | (ls0 << 2 & 4)));
1114	+	bl=bh=0;
1115	+	/* Test t0t1,t1t2,t2t0 against a */
1116	+	if(test_t0t1 && (e0 & 2) )
1117	+	{
1118	+	db = t0[1] + dt10[1]*((double)(a-t0[0])/dt10[0]);
1119	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1120	+	}
1121	+	test_t1t2= !(((ls0 & 3)==0) || ((ls1 & 3)==0)|| ((ls2 & 3)==0) ||
1122	+	((sq0 & 3)==0) ||((sq1 & 3)==0)|| ((sq2 & 3)==0));
1123	+	if(test_t1t2 && (e0 & 1))
1124	+	{    /* test t1t2 against a */
1125	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1126	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1127	+	}
1128	+	test_t2t0 = !(((ls0 & 5)==0) || ((ls1 & 5)==0)|| ((ls2 & 5)==0) ||
1129	+	((sq0 & 5)==0) ||((sq1 & 5)==0)|| ((sq2 & 5)==0));
1130	+	if(test_t2t0 && (e0 & 4))
1131	+	{
1132	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1133	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1134	+	}
1135	+	cl = ch = 0;
1136	+	e1  = (ls1 ^ ((ls1 >> 1) | (ls1 << 2 & 4)));
1137	+	if(test_t0t1 && (e1 & 2))
1138	+	{/* test t0t1 against b */
1139	+	db = t0[2] + dt10[2]*(((double)(b-t0[1]))/dt10[1]);
1140	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1141	+	}
1142	+	if(test_t1t2 && (e1 & 1))
1143	+	{/* test t1t2 against b */
1144	+	db = t1[2] + dt21[2]*((double)(b - t1[1]))/dt21[1];
1145	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1146	+	}
1147	+	if(test_t2t0 && (e1 & 4))
1148	+	{/* test t2t0 against b */
1149	+	db = t2[2] + dt02[2]*((double)(b - t2[1]))/dt02[1];
1150	+	TEST_INT(testl,testh,db,q1[2],c,cl,ch)
1151	+	}
1152	+	al = ah = 0;
1153	+	e2  = (ls2 ^ ((ls2 >> 1) | (ls2 << 2 & 4)));
1154	+	if(test_t0t1 && (e2 & 2))
1155	+	{ /* test t0t1 against c */
1156	+	db = t0[0] + dt10[0]*(((double)(c-t0[2]))/dt10[2]);
1157	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1158	+	}
1159	+	if(test_t1t2 && (e2 & 1))
1160	+	{
1161	+	db = t1[0] + dt21[0]*((double)(c - t1[2]))/dt21[2];
1162	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1163	+	}
1164	+	if(test_t2t0 && (e2 & 4))
1165	+	{ /* test t2t0 against c */
1166	+	db = t2[0] + dt02[0]*((double)(c - t2[2]))/dt02[2];
1167	+	TEST_INT(testl,testh,db,q0[0],a,al,ah)
1168	+	}
1169	+	/* Only remaining case is if some edge trivially rejected */
1170	+	if(!e0 || !e1 || !e2)
1171	+	return;
1172	+
1173	+	/* Only one/none got tested - pick either of other two/any two */
1174	+	/* Only need to test one edge */
1175	+	if(!test_t0t1 && (e0 & 2))
1176	+	{
1177	+	db = t0[1] + dt10[1]*((double)(a-t0[0]))/dt10[0];
1178	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1179	+	}
1180	+	if(!test_t1t2 && (e0 & 1))
1181	+	{/* test t1t2 against a */
1182	+	db = t1[1] + dt21[1]*((double)(a - t1[0]))/dt21[0];
1183	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1184	+	}
1185	+	if(!test_t2t0 && (e0 & 4))
1186	+	{
1187	+	db = t2[1] + dt02[1]*((double)(a - t2[0]))/dt02[0];
1188	+	TEST_INT(testl,testh,db,q1[1],b,bl,bh)
1189	+	}
1190	+
1191	+	return;
1192	+
1193	+	Lfunc_call:
1194	+	f.func(f.argptr,root,qt);
1195	+	if(!QT_IS_EMPTY(qt))
1196	+	QT_LEAF_SET_FLAG(qt);
1197	+	}
1198
1199
1200

Diff Legend

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