ViewVC Help

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

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

Comparing ray/src/hd/sm_geom.c (file contents):
Revision 3.3 by gwlarson, Tue Aug 25 11:03:28 1998 UTC vs.
Revision 3.8 by gwlarson, Wed Nov 11 12:05:38 1998 UTC

#	Line 27 | Line 27 | FVECT v1,v2;
27		{
28		return(EQUAL(v1[0],v2[0]) && EQUAL(v1[1],v2[1])&& EQUAL(v1[2],v2[2]));
29		}
30	+	#if 0
31	+	extern FVECT Norm[500];
32	+	extern int Ncnt;
33	+	#endif
34
31	–
35		int
36		convex_angle(v0,v1,v2)
37		FVECT v0,v1,v2;
38		{
39	<	FVECT cp01,cp12,cp;
40	<
39	>	FVECT cp,cp01,cp12,v10,v02;
40	>	double dp;
41	>
42		/* test sign of (v0Xv1)X(v1Xv2). v1 */
43		VCROSS(cp01,v0,v1);
44		VCROSS(cp12,v1,v2);
45		VCROSS(cp,cp01,cp12);
46	<	if(DOT(cp,v1) < 0)
47	<	return(FALSE);
46	>
47	>	dp = DOT(cp,v1);
48	>	#if 0
49	>	VCOPY(Norm[Ncnt++],cp01);
50	>	VCOPY(Norm[Ncnt++],cp12);
51	>	VCOPY(Norm[Ncnt++],cp);
52	>	#endif
53	>	if(ZERO(dp) || dp < 0.0)
54	>	return(FALSE);
55		return(TRUE);
56		}
57
58		/* calculates the normal of a face contour using Newell's formula. e
59
60	<	a =  SUMi (yi - yi+1)(zi + zi+1)
60	>	a =  SUMi (yi - yi+1)(zi + zi+1)smMesh->samples->max_samp+4);
61		b =  SUMi (zi - zi+1)(xi + xi+1)
62		c =  SUMi (xi - xi+1)(yi + yi+1)
63		*/
64		double
65		tri_normal(v0,v1,v2,n,norm)
66		FVECT v0,v1,v2,n;
67	<	char norm;
67	>	int norm;
68		{
69		double mag;
70
#	Line 64 | Line 75 | char norm;
75		n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
76		(v1[2] - v2[2]) * (v1[0] + v2[0]) +
77		(v2[2] - v0[2]) * (v2[0] + v0[0]);
67	–
78
79		n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
80		(v1[1] + v2[1]) * (v1[0] - v2[0]) +
#	Line 72 | Line 82 | char norm;
82
83		if(!norm)
84		return(0);
75	–
85
86		mag = normalize(n);
87
#	Line 80 | Line 89 | char norm;
89		}
90
91
92	<	tri_plane_equation(v0,v1,v2,n,nd,norm)
93	<	FVECT v0,v1,v2,n;
94	<	double *nd;
95	<	char norm;
92	>	tri_plane_equation(v0,v1,v2,peqptr,norm)
93	>	FVECT v0,v1,v2;
94	>	FPEQ *peqptr;
95	>	int norm;
96		{
97	<	tri_normal(v0,v1,v2,n,norm);
97	>	tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
98
99	<	*nd = -(DOT(n,v0));
99	>	FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
100		}
101
93	–	int
94	–	point_relative_to_plane(p,n,nd)
95	–	FVECT p,n;
96	–	double nd;
97	–	{
98	–	double d;
99	–
100	–	d = p[0]*n[0] + p[1]*n[1] + p[2]*n[2] + nd;
101	–	if(d < 0)
102	–	return(-1);
103	–	if(ZERO(d))
104	–	return(0);
105	–	else
106	–	return(1);
107	–	}
108	–
102		/* From quad_edge-code */
103		int
104		point_in_circle_thru_origin(p,p0,p1)
#	Line 135 | Line 128 | FVECT ps,p,c;
128		}
129
130
131	+	/* returns TRUE if ray from origin in direction v intersects plane defined
132	+	* by normal plane_n, and plane_d. If plane is not parallel- returns
133	+	* intersection point if r != NULL. If tptr!= NULL returns value of
134	+	* t, if parallel, returns t=FHUGE
135	+	*/
136		int
137	<	intersect_vector_plane(v,plane_n,plane_d,tptr,r)
138	<	FVECT v,plane_n;
139	<	double plane_d;
137	>	intersect_vector_plane(v,peq,tptr,r)
138	>	FVECT v;
139	>	FPEQ peq;
140		double *tptr;
141		FVECT r;
142		{
143	<	double t;
143	>	double t,d;
144		int hit;
145		/*
146		Plane is Ax + By + Cz +D = 0:
#	Line 152 | Line 150 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
150		/* line is  l = p1 + (p2-p1)t, p1=origin */
151
152		/* Solve for t: */
153	<	t =  plane_d/-(DOT(plane_n,v));
154	<	if(t >0 || ZERO(t))
155	<	hit = 1;
156	<	else
157	<	hit = 0;
158	<	r[0] = v[0]*t;
159	<	r[1] = v[1]*t;
160	<	r[2] = v[2]*t;
153	>	d = -(DOT(FP_N(peq),v));
154	>	if(ZERO(d))
155	>	{
156	>	t = FHUGE;
157	>	hit = 0;
158	>	}
159	>	else
160	>	{
161	>	t =  FP_D(peq)/d;
162	>	if(t < 0 )
163	>	hit = 0;
164	>	else
165	>	hit = 1;
166	>	if(r)
167	>	{
168	>	r[0] = v[0]*t;
169	>	r[1] = v[1]*t;
170	>	r[2] = v[2]*t;
171	>	}
172	>	}
173		if(tptr)
174		*tptr = t;
175		return(hit);
176		}
177
178		int
179	<	intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
179	>	intersect_ray_plane(orig,dir,peq,pd,r)
180		FVECT orig,dir;
181	<	FVECT plane_n;
172	<	double plane_d;
181	>	FPEQ peq;
182		double *pd;
183		FVECT r;
184		{
185	<	double t;
185	>	double t,d;
186		int hit;
187		/*
188		Plane is Ax + By + Cz +D = 0:
#	Line 184 | Line 193 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
193		line is  l = p1 + (p2-p1)t
194		*/
195		/* Solve for t: */
196	<	t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
197	<	if(ZERO(t) || t >0)
196	>	d = DOT(FP_N(peq),dir);
197	>	if(ZERO(d))
198	>	return(0);
199	>	t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
200	>
201	>	if(t < 0)
202	>	hit = 0;
203	>	else
204		hit = 1;
205	+
206	+	if(r)
207	+	VSUM(r,orig,dir,t);
208	+
209	+	if(pd)
210	+	*pd = t;
211	+	return(hit);
212	+	}
213	+
214	+
215	+	int
216	+	intersect_ray_oplane(orig,dir,n,pd,r)
217	+	FVECT orig,dir;
218	+	FVECT n;
219	+	double *pd;
220	+	FVECT r;
221	+	{
222	+	double t,d;
223	+	int hit;
224	+	/*
225	+	Plane is Ax + By + Cz +D = 0:
226	+	plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
227	+	*/
228	+	/*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
229	+	t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
230	+	line is  l = p1 + (p2-p1)t
231	+	*/
232	+	/* Solve for t: */
233	+	d= DOT(n,dir);
234	+	if(ZERO(d))
235	+	return(0);
236	+	t =  -(DOT(n,orig))/d;
237	+	if(t < 0)
238	+	hit = 0;
239		else
240	+	hit = 1;
241	+
242	+	if(r)
243	+	VSUM(r,orig,dir,t);
244	+
245	+	if(pd)
246	+	*pd = t;
247	+	return(hit);
248	+	}
249	+
250	+
251	+	int
252	+	intersect_edge_plane(e0,e1,peq,pd,r)
253	+	FVECT e0,e1;
254	+	FPEQ peq;
255	+	double *pd;
256	+	FVECT r;
257	+	{
258	+	double t;
259	+	int hit;
260	+	FVECT d;
261	+	/*
262	+	Plane is Ax + By + Cz +D = 0:
263	+	plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
264	+	*/
265	+	/*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
266	+	t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
267	+	line is  l = p1 + (p2-p1)t
268	+	*/
269	+	/* Solve for t: */
270	+	VSUB(d,e1,e0);
271	+	t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
272	+	if(t < 0)
273		hit = 0;
274	+	else
275	+	hit = 1;
276
277	<	VSUM(r,orig,dir,t);
277	>	VSUM(r,e0,d,t);
278
279		if(pd)
280		*pd = t;
#	Line 203 | Line 287 | point_in_cone(p,p0,p1,p2)
287		FVECT p;
288		FVECT p0,p1,p2;
289		{
206	–	FVECT n;
290		FVECT np,x_axis,y_axis;
291	<	double d1,d2,d;
291	>	double d1,d2;
292	>	FPEQ peq;
293
294		/* Find the equation of the circle defined by the intersection
295		of the cone with the plane defined by p1,p2,p3- project p into
296		that plane and do an in-circle test in the plane
297		*/
298
299	<	/* find the equation of the plane defined by p1-p3 */
300	<	tri_plane_equation(p0,p1,p2,n,&d,FALSE);
299	>	/* find the equation of the plane defined by p0-p2 */
300	>	tri_plane_equation(p0,p1,p2,&peq,FALSE);
301
302		/* define a coordinate system on the plane: the x axis is in
303		the direction of np2-np1, and the y axis is calculated from
304		n cross x-axis
305		*/
306		/* Project p onto the plane */
307	<	if(!intersect_vector_plane(p,n,d,NULL,np))
307	>	/* NOTE: check this: does sideness check?*/
308	>	if(!intersect_vector_plane(p,peq,NULL,np))
309		return(FALSE);
310
311	<	/* create coordinate system on  plane: p2-p1 defines the x_axis*/
311	>	/* create coordinate system on  plane: p1-p0 defines the x_axis*/
312		VSUB(x_axis,p1,p0);
313		normalize(x_axis);
314		/* The y axis is  */
315	<	VCROSS(y_axis,n,x_axis);
315	>	VCROSS(y_axis,FP_N(peq),x_axis);
316		normalize(y_axis);
317
318		VSUB(p1,p1,p0);
#	Line 252 | Line 337 | FVECT p0,p1,p2;
337		}
338
339		int
340	<	test_point_against_spherical_tri(v0,v1,v2,p,n,nset,which,sides)
340	>	point_set_in_stri(v0,v1,v2,p,n,nset,sides)
341		FVECT v0,v1,v2,p;
342		FVECT n[3];
343	<	char *nset;
344	<	char *which;
260	<	char sides[3];
343	>	int *nset;
344	>	int sides[3];
345
346		{
347	<	float d;
264	<
347	>	double d;
348		/* Find the normal to the triangle ORIGIN:v0:v1 */
349		if(!NTH_BIT(*nset,0))
350		{
351	<	VCROSS(n[0],v1,v0);
351	>	VCROSS(n[0],v0,v1);
352		SET_NTH_BIT(*nset,0);
353		}
354		/* Test the point for sidedness */
355		d  = DOT(n[0],p);
356
357	<	if(ZERO(d))
358	<	sides[0] = GT_EDGE;
359	<	else
360	<	if(d > 0)
361	<	{
279	<	sides[0] =  GT_OUT;
280	<	sides[1] = sides[2] = GT_INVALID;
281	<	return(FALSE);
357	>	if(d > 0.0)
358	>	{
359	>	sides[0] =  GT_OUT;
360	>	sides[1] = sides[2] = GT_INVALID;
361	>	return(FALSE);
362		}
363		else
364		sides[0] = GT_INTERIOR;
#	Line 286 | Line 366 | char sides[3];
366		/* Test next edge */
367		if(!NTH_BIT(*nset,1))
368		{
369	<	VCROSS(n[1],v2,v1);
369	>	VCROSS(n[1],v1,v2);
370		SET_NTH_BIT(*nset,1);
371		}
372		/* Test the point for sidedness */
373		d  = DOT(n[1],p);
374	<	if(ZERO(d))
374	>	if(d > 0.0)
375		{
296	–	sides[1] = GT_EDGE;
297	–	/* If on plane 0-and on plane 1: lies on edge */
298	–	if(sides[0] == GT_EDGE)
299	–	{
300	–	*which = 1;
301	–	sides[2] = GT_INVALID;
302	–	return(GT_EDGE);
303	–	}
304	–	}
305	–	else if(d > 0)
306	–	{
376		sides[1] = GT_OUT;
377		sides[2] = GT_INVALID;
378		return(FALSE);
#	Line 313 | Line 382 | char sides[3];
382		/* Test next edge */
383		if(!NTH_BIT(*nset,2))
384		{
385	<
317	<	VCROSS(n[2],v0,v2);
385	>	VCROSS(n[2],v2,v0);
386		SET_NTH_BIT(*nset,2);
387		}
388		/* Test the point for sidedness */
389		d  = DOT(n[2],p);
390	<	if(ZERO(d))
390	>	if(d > 0.0)
391		{
392	<	sides[2] = GT_EDGE;
393	<
326	<	/* If on plane 0 and 2: lies on edge 0*/
327	<	if(sides[0] == GT_EDGE)
328	<	{
329	<	*which = 0;
330	<	return(GT_EDGE);
331	<	}
332	<	/* If on plane 1 and 2: lies on edge  2*/
333	<	if(sides[1] == GT_EDGE)
334	<	{
335	<	*which = 2;
336	<	return(GT_EDGE);
337	<	}
338	<	/* otherwise: on face 2 */
339	<	else
340	<	{
341	<	*which = 2;
342	<	return(GT_FACE);
343	<	}
392	>	sides[2] = GT_OUT;
393	>	return(FALSE);
394		}
345	–	else if(d > 0)
346	–	{
347	–	sides[2] = GT_OUT;
348	–	return(FALSE);
349	–	}
350	–	/* If on edge */
395		else
396		sides[2] = GT_INTERIOR;
353	–
354	–	/* If on plane 0 only: on face 0 */
355	–	if(sides[0] == GT_EDGE)
356	–	{
357	–	*which = 0;
358	–	return(GT_FACE);
359	–	}
360	–	/* If on plane 1 only: on face 1 */
361	–	if(sides[1] == GT_EDGE)
362	–	{
363	–	*which = 1;
364	–	return(GT_FACE);
365	–	}
397		/* Must be interior to the pyramid */
398		return(GT_INTERIOR);
399		}
400
401
402
403	<
403	>
404		int
405	<	test_single_point_against_spherical_tri(v0,v1,v2,p,which)
405	>	point_in_stri(v0,v1,v2,p)
406		FVECT v0,v1,v2,p;
376	–	char *which;
407		{
408	<	float d;
408	>	double d;
409		FVECT n;
380	–	char sides[3];
410
411	<	/* First test if point coincides with any of the vertices */
383	<	if(EQUAL_VEC3(p,v0))
384	<	{
385	<	*which = 0;
386	<	return(GT_VERTEX);
387	<	}
388	<	if(EQUAL_VEC3(p,v1))
389	<	{
390	<	*which = 1;
391	<	return(GT_VERTEX);
392	<	}
393	<	if(EQUAL_VEC3(p,v2))
394	<	{
395	<	*which = 2;
396	<	return(GT_VERTEX);
397	<	}
398	<	VCROSS(n,v1,v0);
411	>	VCROSS(n,v0,v1);
412		/* Test the point for sidedness */
413		d  = DOT(n,p);
414	<	if(ZERO(d))
415	<	sides[0] = GT_EDGE;
416	<	else
404	<	if(d > 0)
405	<	return(FALSE);
406	<	else
407	<	sides[0] = GT_INTERIOR;
414	>	if(d > 0.0)
415	>	return(FALSE);
416	>
417		/* Test next edge */
418	<	VCROSS(n,v2,v1);
418	>	VCROSS(n,v1,v2);
419		/* Test the point for sidedness */
420		d  = DOT(n,p);
421	<	if(ZERO(d))
413	<	{
414	<	sides[1] = GT_EDGE;
415	<	/* If on plane 0-and on plane 1: lies on edge */
416	<	if(sides[0] == GT_EDGE)
417	<	{
418	<	*which = 1;
419	<	return(GT_VERTEX);
420	<	}
421	<	}
422	<	else if(d > 0)
421	>	if(d > 0.0)
422		return(FALSE);
424	–	else
425	–	sides[1] = GT_INTERIOR;
423
424		/* Test next edge */
425	<	VCROSS(n,v0,v2);
425	>	VCROSS(n,v2,v0);
426		/* Test the point for sidedness */
427		d  = DOT(n,p);
428	<	if(ZERO(d))
432	<	{
433	<	sides[2] = GT_EDGE;
434	<
435	<	/* If on plane 0 and 2: lies on edge 0*/
436	<	if(sides[0] == GT_EDGE)
437	<	{
438	<	*which = 0;
439	<	return(GT_VERTEX);
440	<	}
441	<	/* If on plane 1 and 2: lies on edge  2*/
442	<	if(sides[1] == GT_EDGE)
443	<	{
444	<	*which = 2;
445	<	return(GT_VERTEX);
446	<	}
447	<	/* otherwise: on face 2 */
448	<	else
449	<	{
450	<	return(GT_FACE);
451	<	}
452	<	}
453	<	else if(d > 0)
428	>	if(d > 0.0)
429		return(FALSE);
430		/* Must be interior to the pyramid */
431	<	return(GT_FACE);
431	>	return(GT_INTERIOR);
432		}
433
434		int
435	<	test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
435	>	vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
436		FVECT t0,t1,t2,p0,p1,p2;
437	<	char *nset;
437	>	int *nset;
438		FVECT n[3];
439		FVECT avg;
440	<	char pt_sides[3][3];
440	>	int pt_sides[3][3];
441
442		{
443	<	char below_plane[3],on_edge,test;
469	<	char which;
443	>	int below_plane[3],test;
444
445		SUM_3VEC3(avg,t0,t1,t2);
472	–	on_edge = 0;
446		*nset = 0;
447		/* Test vertex v[i] against triangle j*/
448		/* Check if v[i] lies below plane defined by avg of 3 vectors
#	Line 477 | Line 450 | char pt_sides[3][3];
450		*/
451
452		/* test point 0 */
453	<	if(DOT(avg,p0) < 0)
453	>	if(DOT(avg,p0) < 0.0)
454		below_plane[0] = 1;
455		else
456	<	below_plane[0]=0;
456	>	below_plane[0] = 0;
457		/* Test if b[i] lies in or on triangle a */
458	<	test = test_point_against_spherical_tri(t0,t1,t2,p0,
486	<	n,nset,&which,pt_sides[0]);
458	>	test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
459		/* If pts[i] is interior: done */
460		if(!below_plane[0])
461		{
462		if(test == GT_INTERIOR)
463		return(TRUE);
492	–	/* Remember if b[i] fell on one of the 3 defining planes */
493	–	if(test)
494	–	on_edge++;
464		}
465		/* Now test point 1*/
466
467	<	if(DOT(avg,p1) < 0)
467	>	if(DOT(avg,p1) < 0.0)
468		below_plane[1] = 1;
469		else
470		below_plane[1]=0;
471		/* Test if b[i] lies in or on triangle a */
472	<	test = test_point_against_spherical_tri(t0,t1,t2,p1,
504	<	n,nset,&which,pt_sides[1]);
472	>	test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
473		/* If pts[i] is interior: done */
474		if(!below_plane[1])
475		{
476		if(test == GT_INTERIOR)
477		return(TRUE);
510	–	/* Remember if b[i] fell on one of the 3 defining planes */
511	–	if(test)
512	–	on_edge++;
478		}
479
480		/* Now test point 2 */
481	<	if(DOT(avg,p2) < 0)
481	>	if(DOT(avg,p2) < 0.0)
482		below_plane[2] = 1;
483		else
484	<	below_plane[2]=0;
484	>	below_plane[2] = 0;
485		/* Test if b[i] lies in or on triangle a */
486	<	test = test_point_against_spherical_tri(t0,t1,t2,p2,
522	<	n,nset,&which,pt_sides[2]);
486	>	test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
487
488		/* If pts[i] is interior: done */
489		if(!below_plane[2])
490		{
491		if(test == GT_INTERIOR)
492		return(TRUE);
529	–	/* Remember if b[i] fell on one of the 3 defining planes */
530	–	if(test)
531	–	on_edge++;
493		}
494
495		/* If all three points below separating plane: trivial reject */
496		if(below_plane[0] && below_plane[1] && below_plane[2])
497		return(FALSE);
537	–	/* Accept if all points lie on a triangle vertex/edge edge- accept*/
538	–	if(on_edge == 3)
539	–	return(TRUE);
498		/* Now check vertices in a against triangle b */
499		return(FALSE);
500		}
#	Line 544 | Line 502 | char pt_sides[3][3];
502
503		set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
504		FVECT t0,t1,t2,p0,p1,p2;
505	<	char test[3];
506	<	char sides[3][3];
507	<	char nset;
505	>	int test[3];
506	>	int sides[3][3];
507	>	int nset;
508		FVECT n[3];
509		{
510	<	char t;
510	>	int t;
511		double d;
512
513
#	Line 558 | Line 516 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
516		if(sides[0][0] == GT_INVALID)
517		{
518		if(!NTH_BIT(nset,0))
519	<	VCROSS(n[0],t1,t0);
519	>	VCROSS(n[0],t0,t1);
520		/* Test the point for sidedness */
521		d  = DOT(n[0],p0);
522	<	if(d >= 0)
522	>	if(d >= 0.0)
523		SET_NTH_BIT(test[0],0);
524		}
525		else
#	Line 571 | Line 529 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
529		if(sides[0][1] == GT_INVALID)
530		{
531		if(!NTH_BIT(nset,1))
532	<	VCROSS(n[1],t2,t1);
532	>	VCROSS(n[1],t1,t2);
533		/* Test the point for sidedness */
534		d  = DOT(n[1],p0);
535	<	if(d >= 0)
535	>	if(d >= 0.0)
536		SET_NTH_BIT(test[0],1);
537		}
538		else
#	Line 584 | Line 542 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
542		if(sides[0][2] == GT_INVALID)
543		{
544		if(!NTH_BIT(nset,2))
545	<	VCROSS(n[2],t0,t2);
545	>	VCROSS(n[2],t2,t0);
546		/* Test the point for sidedness */
547		d  = DOT(n[2],p0);
548	<	if(d >= 0)
548	>	if(d >= 0.0)
549		SET_NTH_BIT(test[0],2);
550		}
551		else
#	Line 600 | Line 558 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
558		if(sides[1][0] == GT_INVALID)
559		{
560		if(!NTH_BIT(nset,0))
561	<	VCROSS(n[0],t1,t0);
561	>	VCROSS(n[0],t0,t1);
562		/* Test the point for sidedness */
563		d  = DOT(n[0],p1);
564	<	if(d >= 0)
564	>	if(d >= 0.0)
565		SET_NTH_BIT(test[1],0);
566		}
567		else
#	Line 614 | Line 572 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
572		if(sides[1][1] == GT_INVALID)
573		{
574		if(!NTH_BIT(nset,1))
575	<	VCROSS(n[1],t2,t1);
575	>	VCROSS(n[1],t1,t2);
576		/* Test the point for sidedness */
577		d  = DOT(n[1],p1);
578	<	if(d >= 0)
578	>	if(d >= 0.0)
579		SET_NTH_BIT(test[1],1);
580		}
581		else
#	Line 628 | Line 586 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
586		if(sides[1][2] == GT_INVALID)
587		{
588		if(!NTH_BIT(nset,2))
589	<	VCROSS(n[2],t0,t2);
589	>	VCROSS(n[2],t2,t0);
590		/* Test the point for sidedness */
591		d  = DOT(n[2],p1);
592	<	if(d >= 0)
592	>	if(d >= 0.0)
593		SET_NTH_BIT(test[1],2);
594		}
595		else
#	Line 644 | Line 602 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
602		if(sides[2][0] == GT_INVALID)
603		{
604		if(!NTH_BIT(nset,0))
605	<	VCROSS(n[0],t1,t0);
605	>	VCROSS(n[0],t0,t1);
606		/* Test the point for sidedness */
607		d  = DOT(n[0],p2);
608	<	if(d >= 0)
608	>	if(d >= 0.0)
609		SET_NTH_BIT(test[2],0);
610		}
611		else
#	Line 657 | Line 615 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
615		if(sides[2][1] == GT_INVALID)
616		{
617		if(!NTH_BIT(nset,1))
618	<	VCROSS(n[1],t2,t1);
618	>	VCROSS(n[1],t1,t2);
619		/* Test the point for sidedness */
620		d  = DOT(n[1],p2);
621	<	if(d >= 0)
621	>	if(d >= 0.0)
622		SET_NTH_BIT(test[2],1);
623		}
624		else
#	Line 670 | Line 628 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
628		if(sides[2][2] == GT_INVALID)
629		{
630		if(!NTH_BIT(nset,2))
631	<	VCROSS(n[2],t0,t2);
631	>	VCROSS(n[2],t2,t0);
632		/* Test the point for sidedness */
633		d  = DOT(n[2],p2);
634	<	if(d >= 0)
634	>	if(d >= 0.0)
635		SET_NTH_BIT(test[2],2);
636		}
637		else
#	Line 683 | Line 641 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
641
642
643		int
644	<	spherical_tri_intersect(a1,a2,a3,b1,b2,b3)
644	>	stri_intersect(a1,a2,a3,b1,b2,b3)
645		FVECT a1,a2,a3,b1,b2,b3;
646		{
647	<	char which,test,n_set[2];
648	<	char sides[2][3][3],i,j,inext,jnext;
649	<	char tests[2][3];
647	>	int which,test,n_set[2];
648	>	int sides[2][3][3],i,j,inext,jnext;
649	>	int tests[2][3];
650		FVECT n[2][3],p,avg[2];
651
652		/* Test the vertices of triangle a against the pyramid formed by triangle
#	Line 696 | Line 654 | FVECT a1,a2,a3,b1,b2,b3;
654		if all 3 vertices of a are ON the edges of b,return TRUE. Remember
655		the results of the edge normal and sidedness tests for later.
656		*/
657	<	if(test_vertices_for_tri_inclusion(a1,a2,a3,b1,b2,b3,
700	<	&(n_set[0]),n[0],avg[0],sides[1]))
657	>	if(vertices_in_stri(a1,a2,a3,b1,b2,b3,&(n_set[0]),n[0],avg[0],sides[1]))
658		return(TRUE);
659
660	<	if(test_vertices_for_tri_inclusion(b1,b2,b3,a1,a2,a3,
704	<	&(n_set[1]),n[1],avg[1],sides[0]))
660	>	if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
661		return(TRUE);
662
663
#	Line 749 | Line 705 | FVECT a1,a2,a3,b1,b2,b3;
705		}
706
707		int
708	<	ray_intersect_tri(orig,dir,v0,v1,v2,pt,wptr)
708	>	ray_intersect_tri(orig,dir,v0,v1,v2,pt)
709		FVECT orig,dir;
710		FVECT v0,v1,v2;
711		FVECT pt;
756	–	char *wptr;
712		{
713	<	FVECT p0,p1,p2,p,n;
714	<	char type,which;
715	<	double pd;
761	<
762	<	point_on_sphere(p0,v0,orig);
763	<	point_on_sphere(p1,v1,orig);
764	<	point_on_sphere(p2,v2,orig);
765	<	type = test_single_point_against_spherical_tri(p0,p1,p2,dir,&which);
713	>	FVECT p0,p1,p2,p;
714	>	FPEQ peq;
715	>	int type;
716
717	<	if(type)
717	>	VSUB(p0,v0,orig);
718	>	VSUB(p1,v1,orig);
719	>	VSUB(p2,v2,orig);
720	>
721	>	if(point_in_stri(p0,p1,p2,dir))
722		{
723		/* Intersect the ray with the triangle plane */
724	<	tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
725	<	intersect_ray_plane(orig,dir,n,pd,NULL,pt);
724	>	tri_plane_equation(v0,v1,v2,&peq,FALSE);
725	>	return(intersect_ray_plane(orig,dir,peq,NULL,pt));
726		}
727	<	if(wptr)
774	<	*wptr = which;
775	<
776	<	return(type);
727	>	return(FALSE);
728		}
729
730
#	Line 832 | Line 783 | FVECT fnear[4],ffar[4];
783		ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
784		}
785
786	+	int
787	+	max_index(v,r)
788	+	FVECT v;
789	+	double *r;
790	+	{
791	+	double p[3];
792	+	int i;
793
794	+	p[0] = fabs(v[0]);
795	+	p[1] = fabs(v[1]);
796	+	p[2] = fabs(v[2]);
797	+	i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);
798	+	if(r)
799	+	*r = p[i];
800	+	return(i);
801	+	}
802
803	+	int
804	+	closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
805	+	FVECT p0,p1,p2,p;
806	+	int p0id,p1id,p2id;
807	+	{
808	+	double d,d1;
809	+	int i;
810	+
811	+	d =  DIST_SQ(p,p0);
812	+	d1 = DIST_SQ(p,p1);
813	+	if(d < d1)
814	+	{
815	+	d1 = DIST_SQ(p,p2);
816	+	i = (d1 < d)?p2id:p0id;
817	+	}
818	+	else
819	+	{
820	+	d = DIST_SQ(p,p2);
821	+	i = (d < d1)? p2id:p1id;
822	+	}
823	+	return(i);
824	+	}
825
826	+
827		int
828	<	spherical_polygon_edge_intersect(a0,a1,b0,b1)
828	>	sedge_intersect(a0,a1,b0,b1)
829		FVECT a0,a1,b0,b1;
830		{
831		FVECT na,nb,avga,avgb,p;
#	Line 896 | Line 885 | double coord[3];
885		a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
886		coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
887		coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
888	<	coord[2]  = 1.0 - coord[0] - coord[1];
888	>	coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
889
890		}
891
892	+	bary_ith_child(coord,i)
893	+	double coord[3];
894	+	int i;
895	+	{
896	+
897	+	switch(i){
898	+	case 0:
899	+	/* update bary for child */
900	+	coord[0] = 2.0*coord[0]- 1.0;
901	+	coord[1] *= 2.0;
902	+	coord[2] *= 2.0;
903	+	break;
904	+	case 1:
905	+	coord[0] *= 2.0;
906	+	coord[1] = 2.0*coord[1]- 1.0;
907	+	coord[2] *= 2.0;
908	+	break;
909	+	case 2:
910	+	coord[0] *= 2.0;
911	+	coord[1] *= 2.0;
912	+	coord[2] = 2.0*coord[2]- 1.0;
913	+	break;
914	+	case 3:
915	+	coord[0] = 1.0 - 2.0*coord[0];
916	+	coord[1] = 1.0 - 2.0*coord[1];
917	+	coord[2] = 1.0 - 2.0*coord[2];
918	+	break;
919	+	#ifdef DEBUG
920	+	default:
921	+	eputs("bary_ith_child():Invalid child\n");
922	+	break;
923	+	#endif
924	+	}
925	+	}
926	+
927	+
928		int
929	<	bary2d_child(coord)
929	>	bary_child(coord)
930		double coord[3];
931		{
932		int i;
933
909	–	/* First check if one of the original vertices */
910	–	for(i=0;i<3;i++)
911	–	if(EQUAL(coord[i],1.0))
912	–	return(i);
913	–
914	–	/* Check if one of the new vertices: for all return center child */
915	–	if(ZERO(coord[0]) && EQUAL(coord[1],0.5))
916	–	{
917	–	coord[0] = 1.0f;
918	–	coord[1] = 0.0f;
919	–	coord[2] = 0.0f;
920	–	return(3);
921	–	}
922	–	if(ZERO(coord[1]) && EQUAL(coord[0],0.5))
923	–	{
924	–	coord[0] = 0.0f;
925	–	coord[1] = 1.0f;
926	–	coord[2] = 0.0f;
927	–	return(3);
928	–	}
929	–	if(ZERO(coord[2]) && EQUAL(coord[0],0.5))
930	–	{
931	–	coord[0] = 0.0f;
932	–	coord[1] = 0.0f;
933	–	coord[2] = 1.0f;
934	–	return(3);
935	–	}
936	–
937	–	/* Otherwise return child */
934		if(coord[0] > 0.5)
935		{
936		/* update bary for child */
#	Line 968 | Line 964 | double coord[3];
964		}
965		}
966
967	<	int
968	<	max_index(v)
969	<	FVECT v;
967	>	/* Coord was the ith child of the parent: set the coordinate
968	>	relative to the parent
969	>	*/
970	>	bary_parent(coord,i)
971	>	double coord[3];
972	>	int i;
973		{
975	–	double a,b,c;
976	–	int i;
974
975	<	a = fabs(v[0]);
976	<	b = fabs(v[1]);
977	<	c = fabs(v[2]);
978	<	i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);
979	<	return(i);
975	>	switch(i) {
976	>	case 0:
977	>	/* update bary for child */
978	>	coord[0] = coord[0]*0.5 + 0.5;
979	>	coord[1] *= 0.5;
980	>	coord[2] *= 0.5;
981	>	break;
982	>	case 1:
983	>	coord[0] *= 0.5;
984	>	coord[1]  = coord[1]*0.5 + 0.5;
985	>	coord[2] *= 0.5;
986	>	break;
987	>
988	>	case 2:
989	>	coord[0] *= 0.5;
990	>	coord[1] *= 0.5;
991	>	coord[2] = coord[2]*0.5 + 0.5;
992	>	break;
993	>
994	>	case 3:
995	>	coord[0] = 0.5 - 0.5*coord[0];
996	>	coord[1] = 0.5 - 0.5*coord[1];
997	>	coord[2] = 0.5 - 0.5*coord[2];
998	>	break;
999	>	#ifdef DEBUG
1000	>	default:
1001	>	eputs("bary_parent():Invalid child\n");
1002	>	break;
1003	>	#endif
1004	>	}
1005		}
1006
1007	+	bary_from_child(coord,child,next)
1008	+	double coord[3];
1009	+	int child,next;
1010	+	{
1011	+	#ifdef DEBUG
1012	+	if(child <0 || child > 3)
1013	+	{
1014	+	eputs("bary_from_child():Invalid child\n");
1015	+	return;
1016	+	}
1017	+	if(next <0 || next > 3)
1018	+	{
1019	+	eputs("bary_from_child():Invalid next\n");
1020	+	return;
1021	+	}
1022	+	#endif
1023	+	if(next == child)
1024	+	return;
1025
1026	+	switch(child){
1027	+	case 0:
1028	+	switch(next){
1029	+	case 1:
1030	+	coord[0] += 1.0;
1031	+	coord[1] -= 1.0;
1032	+	break;
1033	+	case 2:
1034	+	coord[0] += 1.0;
1035	+	coord[2] -= 1.0;
1036	+	break;
1037	+	case 3:
1038	+	coord[0] *= -1.0;
1039	+	coord[1] = 1 - coord[1];
1040	+	coord[2] = 1 - coord[2];
1041	+	break;
1042
1043	<	/*
1044	<	* int
1045	<	* smRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
1046	<	*
1047	<	*   Intersect the ray with triangle v0v1v2, return intersection point in r
1048	<	*
1049	<	*    Assumes orig,v0,v1,v2 are in spherical coordinates, and orig is
1050	<	*    unit
1051	<	*/
1043	>	}
1044	>	break;
1045	>	case 1:
1046	>	switch(next){
1047	>	case 0:
1048	>	coord[0] -= 1.0;
1049	>	coord[1] += 1.0;
1050	>	break;
1051	>	case 2:
1052	>	coord[1] += 1.0;
1053	>	coord[2] -= 1.0;
1054	>	break;
1055	>	case 3:
1056	>	coord[0] = 1 - coord[0];
1057	>	coord[1] *= -1.0;
1058	>	coord[2] = 1 - coord[2];
1059	>	break;
1060	>	}
1061	>	break;
1062	>	case 2:
1063	>	switch(next){
1064	>	case 0:
1065	>	coord[0] -= 1.0;
1066	>	coord[2] += 1.0;
1067	>	break;
1068	>	case 1:
1069	>	coord[1] -= 1.0;
1070	>	coord[2] += 1.0;
1071	>	break;
1072	>	case 3:
1073	>	coord[0] = 1 - coord[0];
1074	>	coord[1] = 1 - coord[1];
1075	>	coord[2] *= -1.0;
1076	>	break;
1077	>	}
1078	>	break;
1079	>	case 3:
1080	>	switch(next){
1081	>	case 0:
1082	>	coord[0] *= -1.0;
1083	>	coord[1] = 1 - coord[1];
1084	>	coord[2] = 1 - coord[2];
1085	>	break;
1086	>	case 1:
1087	>	coord[0] = 1 - coord[0];
1088	>	coord[1] *= -1.0;
1089	>	coord[2] = 1 - coord[2];
1090	>	break;
1091	>	case 2:
1092	>	coord[0] = 1 - coord[0];
1093	>	coord[1] = 1 - coord[1];
1094	>	coord[2] *= -1.0;
1095	>	break;
1096	>	}
1097	>	break;
1098	>	}
1099	>	}
1100	>
1101	>
1102	>	baryi_parent(coord,i)
1103	>	BCOORD coord[3];
1104	>	int i;
1105	>	{
1106	>
1107	>	switch(i) {
1108	>	case 0:
1109	>	/* update bary for child */
1110	>	coord[0] = (coord[0] >> 1) + MAXBCOORD2;
1111	>	coord[1] >>= 1;
1112	>	coord[2] >>= 1;
1113	>	break;
1114	>	case 1:
1115	>	coord[0] >>= 1;
1116	>	coord[1]  = (coord[1] >> 1) + MAXBCOORD2;
1117	>	coord[2] >>= 1;
1118	>	break;
1119	>
1120	>	case 2:
1121	>	coord[0] >>= 1;
1122	>	coord[1] >>= 1;
1123	>	coord[2] = (coord[2] >> 1) + MAXBCOORD2;
1124	>	break;
1125	>
1126	>	case 3:
1127	>	coord[0] = MAXBCOORD2 - (coord[0] >> 1);
1128	>	coord[1] = MAXBCOORD2 - (coord[1] >> 1);
1129	>	coord[2] = MAXBCOORD2 - (coord[2] >> 1);
1130	>	break;
1131	>	#ifdef DEBUG
1132	>	default:
1133	>	eputs("baryi_parent():Invalid child\n");
1134	>	break;
1135	>	#endif
1136	>	}
1137	>	}
1138	>
1139	>	baryi_from_child(coord,child,next)
1140	>	BCOORD coord[3];
1141	>	int child,next;
1142	>	{
1143	>	#ifdef DEBUG
1144	>	if(child <0 || child > 3)
1145	>	{
1146	>	eputs("baryi_from_child():Invalid child\n");
1147	>	return;
1148	>	}
1149	>	if(next <0 || next > 3)
1150	>	{
1151	>	eputs("baryi_from_child():Invalid next\n");
1152	>	return;
1153	>	}
1154	>	#endif
1155	>	if(next == child)
1156	>	return;
1157	>
1158	>	switch(child){
1159	>	case 0:
1160	>	coord[0] = 0;
1161	>	coord[1] = MAXBCOORD - coord[1];
1162	>	coord[2] = MAXBCOORD - coord[2];
1163	>	break;
1164	>	case 1:
1165	>	coord[0] = MAXBCOORD - coord[0];
1166	>	coord[1] = 0;
1167	>	coord[2] = MAXBCOORD - coord[2];
1168	>	break;
1169	>	case 2:
1170	>	coord[0] = MAXBCOORD - coord[0];
1171	>	coord[1] = MAXBCOORD - coord[1];
1172	>	coord[2] = 0;
1173	>	break;
1174	>	case 3:
1175	>	switch(next){
1176	>	case 0:
1177	>	coord[0] = 0;
1178	>	coord[1] = MAXBCOORD - coord[1];
1179	>	coord[2] = MAXBCOORD - coord[2];
1180	>	break;
1181	>	case 1:
1182	>	coord[0] = MAXBCOORD - coord[0];
1183	>	coord[1] = 0;
1184	>	coord[2] = MAXBCOORD - coord[2];
1185	>	break;
1186	>	case 2:
1187	>	coord[0] = MAXBCOORD - coord[0];
1188	>	coord[1] = MAXBCOORD - coord[1];
1189	>	coord[2] = 0;
1190	>	break;
1191	>	}
1192	>	break;
1193	>	}
1194	>	}
1195	>
1196		int
1197	<	traceRay(orig,dir,v0,v1,v2,r)
1198	<	FVECT orig,dir;
999	<	FVECT v0,v1,v2;
1000	<	FVECT r;
1197	>	baryi_child(coord)
1198	>	BCOORD coord[3];
1199		{
1200	<	FVECT n,p[3],d;
1003	<	double pt[3],r_eps;
1004	<	char i;
1005	<	int which;
1200	>	int i;
1201
1202	<	/* Find the plane equation for the triangle defined by the edge v0v1 and
1203	<	the view center*/
1204	<	VCROSS(n,v0,v1);
1205	<	/* Intersect the ray with this plane */
1206	<	i = intersect_ray_plane(orig,dir,n,0.0,&(pt[0]),p[0]);
1207	<	if(i)
1208	<	which = 0;
1202	>	if(coord[0] > MAXBCOORD2)
1203	>	{
1204	>	/* update bary for child */
1205	>	coord[0] = (coord[0]<< 1) - MAXBCOORD;
1206	>	coord[1] <<= 1;
1207	>	coord[2] <<= 1;
1208	>	return(0);
1209	>	}
1210		else
1211	<	which = -1;
1211	>	if(coord[1] > MAXBCOORD2)
1212	>	{
1213	>	coord[0] <<= 1;
1214	>	coord[1] = (coord[1] << 1) - MAXBCOORD;
1215	>	coord[2] <<= 1;
1216	>	return(1);
1217	>	}
1218	>	else
1219	>	if(coord[2] > MAXBCOORD2)
1220	>	{
1221	>	coord[0] <<= 1;
1222	>	coord[1] <<= 1;
1223	>	coord[2] = (coord[2] << 1) - MAXBCOORD;
1224	>	return(2);
1225	>	}
1226	>	else
1227	>	{
1228	>	coord[0] = MAXBCOORD - (coord[0] << 1);
1229	>	coord[1] = MAXBCOORD - (coord[1] << 1);
1230	>	coord[2] = MAXBCOORD - (coord[2] << 1);
1231	>	return(3);
1232	>	}
1233	>	}
1234
1235	<	VCROSS(n,v1,v2);
1236	<	i = intersect_ray_plane(orig,dir,n,0.0,&(pt[1]),p[1]);
1237	<	if(i)
1238	<	if((which==-1) || pt[1] < pt[0])
1239	<	which = 1;
1235	>	int
1236	>	baryi_nth_child(coord,i)
1237	>	BCOORD coord[3];
1238	>	int i;
1239	>	{
1240
1241	<	VCROSS(n,v2,v0);
1242	<	i = intersect_ray_plane(orig,dir,n,0.0,&(pt[2]),p[2]);
1243	<	if(i)
1244	<	if((which==-1) || pt[2] < pt[which])
1245	<	which = 2;
1241	>	switch(i){
1242	>	case 0:
1243	>	/* update bary for child */
1244	>	coord[0] = (coord[0]<< 1) - MAXBCOORD;
1245	>	coord[1] <<= 1;
1246	>	coord[2] <<= 1;
1247	>	break;
1248	>	case 1:
1249	>	coord[0] <<= 1;
1250	>	coord[1] = (coord[1] << 1) - MAXBCOORD;
1251	>	coord[2] <<= 1;
1252	>	break;
1253	>	case 2:
1254	>	coord[0] <<= 1;
1255	>	coord[1] <<= 1;
1256	>	coord[2] = (coord[2] << 1) - MAXBCOORD;
1257	>	break;
1258	>	case 3:
1259	>	coord[0] = MAXBCOORD - (coord[0] << 1);
1260	>	coord[1] = MAXBCOORD - (coord[1] << 1);
1261	>	coord[2] = MAXBCOORD - (coord[2] << 1);
1262	>	break;
1263	>	}
1264	>	}
1265
1266	<	if(which != -1)
1266	>
1267	>	baryi_children(coord,i,in_tri,rcoord,rin_tri)
1268	>	BCOORD coord[3][3];
1269	>	int i;
1270	>	int in_tri[3];
1271	>	BCOORD rcoord[3][3];
1272	>	int rin_tri[3];
1273	>	{
1274	>	int j;
1275	>
1276	>	for(j=0; j< 3; j++)
1277		{
1278	<	/* Return point that is K*eps along projection of the ray on
1279	<	the sphere to push intersection point p[which] into next cell
1280	<	*/
1281	<	normalize(p[which]);
1282	<	/* Calculate the ray perpendicular to the intersection point: approx
1283	<	the direction moved along the sphere a distance of K*epsilon*/
1284	<	r_eps = -DOT(p[which],dir);
1285	<	VSUM(n,dir,p[which],r_eps);
1286	<	/* Calculate the length along ray p[which]-dir needed to move to
1287	<	cause a move along the sphere of k*epsilon
1288	<	*/
1289	<	r_eps = DOT(n,dir);
1290	<	VSUM(r,p[which],dir,(20*FTINY)/r_eps);
1291	<	normalize(r);
1292	<	return(TRUE);
1278	>	if(!in_tri[j])
1279	>	{
1280	>	rin_tri[j]=0;
1281	>	continue;
1282	>	}
1283	>
1284	>	if(i != 3)
1285	>	{
1286	>	if(coord[j][i] < MAXBCOORD2)
1287	>	{
1288	>	rin_tri[j] = 0;
1289	>	continue;
1290	>	}
1291	>	}
1292	>	else
1293	>	if( !(coord[j][0] <= MAXBCOORD2 && coord[j][1] <= MAXBCOORD2 &&
1294	>	coord[j][2] <= MAXBCOORD2))
1295	>	{
1296	>	rin_tri[j] = 0;
1297	>	continue;
1298	>	}
1299	>
1300	>	rin_tri[j]=1;
1301	>	VCOPY(rcoord[j],coord[j]);
1302	>	baryi_nth_child(rcoord[j],i);
1303		}
1304	+
1305	+	}
1306	+
1307	+	convert_dtol(b,bi)
1308	+	double b[3];
1309	+	BCOORD bi[3];
1310	+	{
1311	+	int i;
1312	+
1313	+	for(i=0; i < 2;i++)
1314	+	{
1315	+
1316	+	if(b[i] <= 0.0)
1317	+	{
1318	+	bi[i]= 0;
1319	+	}
1320	+	else
1321	+	if(b[i] >= 1.0)
1322	+	{
1323	+	bi[i]= MAXBCOORD;
1324	+	}
1325	+	else
1326	+	bi[i] = (BCOORD)(b[i]*MAXBCOORD);
1327	+	}
1328	+	bi[2] = bi[0] +  bi[1];
1329	+	if(bi[2] > MAXBCOORD)
1330	+	{
1331	+	bi[2] = 0;
1332	+	bi[1] = MAXBCOORD - bi[0];
1333	+	}
1334		else
1335	+	bi[2] = MAXBCOORD - bi[2];
1336	+
1337	+	}
1338	+
1339	+	/* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1340	+	dir unbounded to [-MAXLONG,MAXLONG]
1341	+	*/
1342	+	bary_dtol(b,db,bi,dbi,t,w)
1343	+	double b[3],db[3][3];
1344	+	BCOORD bi[3];
1345	+	BDIR dbi[3][3];
1346	+	TINT t[3];
1347	+	int w;
1348	+	{
1349	+	int i,id,j,k;
1350	+	double d;
1351	+
1352	+	convert_dtol(b,bi);
1353	+
1354	+	for(j=w; j< 3; j++)
1355		{
1356	<	/* Unable to find intersection: move along ray and try again */
1357	<	r_eps = -DOT(orig,dir);
1358	<	VSUM(n,dir,orig,r_eps);
1359	<	r_eps = DOT(n,dir);
1360	<	VSUM(r,orig,dir,(20*FTINY)/r_eps);
1361	<	normalize(r);
1362	<	#ifdef DEBUG
1363	<	eputs("traceRay:Ray does not intersect triangle");
1364	<	#endif
1365	<	return(FALSE);
1356	>	if(t[j] == HUGET)
1357	>	{
1358	>	max_index(db[j],&d);
1359	>	for(i=0; i< 3; i++)
1360	>	dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1361	>	break;
1362	>	}
1363	>	else
1364	>	{
1365	>	for(i=0; i< 3; i++)
1366	>	dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1367	>	}
1368		}
1369		}
1370	+
1371	+
1372	+	/* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1373	+	dir unbounded to [-MAXLONG,MAXLONG]
1374	+	*/
1375	+	bary_dtol_new(b,db,bi,boi,dbi,t)
1376	+	double b[4][3],db[3][3];
1377	+	BCOORD bi[3],boi[3][3];
1378	+	BDIR dbi[3][3];
1379	+	int t[3];
1380	+	{
1381	+	int i,id,j,k;
1382	+	double d;
1383	+
1384	+	convert_dtol(b[3],bi);
1385	+
1386	+	for(j=0; j<3;j++)
1387	+	{
1388	+	if(t[j] != 1)
1389	+	continue;
1390	+	convert_dtol(b[j],boi[j]);
1391	+	}
1392	+	for(j=0; j< 3; j++)
1393	+	{
1394	+	k = (j+1)%3;
1395	+	if(t[k]==0)
1396	+	continue;
1397	+	if(t[k] == -1)
1398	+	{
1399	+	max_index(db[j],&d);
1400	+	for(i=0; i< 3; i++)
1401	+	dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1402	+	t[k] = 0;
1403	+	}
1404	+	else
1405	+	if(t[j] != 1)
1406	+	for(i=0; i< 3; i++)
1407	+	dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1408	+	else
1409	+	for(i=0; i< 3; i++)
1410	+	dbi[j][i] = boi[k][i] - boi[j][i];
1411	+
1412	+	}
1413	+	}
1414	+
1415	+
1416	+	bary_dtolb(b,bi,in_tri)
1417	+	double b[3][3];
1418	+	BCOORD bi[3][3];
1419	+	int in_tri[3];
1420	+	{
1421	+	int i,j;
1422	+
1423	+	for(j=0; j<3;j++)
1424	+	{
1425	+	if(!in_tri[j])
1426	+	continue;
1427	+	for(i=0; i < 2;i++)
1428	+	{
1429	+	if(b[j][i] <= 0.0)
1430	+	{
1431	+	bi[j][i]= 0;
1432	+	}
1433	+	else
1434	+	if(b[j][i] >= 1.0)
1435	+	{
1436	+	bi[j][i]= MAXBCOORD;
1437	+	}
1438	+	else
1439	+	bi[j][i] = (BCOORD)(b[j][i]*MAXBCOORD);
1440	+	}
1441	+	bi[j][2] = MAXBCOORD - bi[j][0] - bi[j][1];
1442	+	if(bi[j][2] < 0)
1443	+	{
1444	+	bi[j][2] = 0;
1445	+	bi[j][1] = MAXBCOORD - bi[j][0];
1446	+	}
1447	+	}
1448	+	}
1449	+
1450	+
1451
1452
1453

Diff Legend

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