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.1 by gwlarson, Wed Aug 19 17:45:24 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,pd,r)
138	<	FVECT v,plane_n;
139	<	double plane_d;
140	<	double *pd;
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,pd,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
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,peq,pd,r)
180	>	FVECT orig,dir;
181	>	FPEQ peq;
182	>	double *pd;
183	>	FVECT r;
184	>	{
185	>	double t,d;
186	>	int hit;
187	>	/*
188	>	Plane is Ax + By + Cz +D = 0:
189	>	plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
190	>	*/
191	>	/*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
192	>	t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
193	>	line is  l = p1 + (p2-p1)t
194	>	*/
195	>	/* Solve for t: */
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	<	r[0] = v[0]*t;
204	<	r[1] = v[1]*t;
205	<	r[2] = v[2]*t;
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_plane(orig,dir,plane_n,plane_d,pd,r)
216	>	intersect_ray_oplane(orig,dir,n,pd,r)
217		FVECT orig,dir;
218	<	FVECT plane_n;
172	<	double plane_d;
218	>	FVECT n;
219		double *pd;
220		FVECT r;
221		{
222	<	double t;
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 */
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	<	t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
234	<	if(ZERO(t) || t >0)
235	<	hit = 1;
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 202 | Line 287 | point_in_cone(p,p0,p1,p2)
287		FVECT p;
288		FVECT p0,p1,p2;
289		{
205	–	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 251 | 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;
259	<	char sides[3];
343	>	int *nset;
344	>	int sides[3];
345
346		{
347	<	float d;
263	<
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	<	{
278	<	sides[0] =  GT_OUT;
279	<	sides[1] = sides[2] = GT_INVALID;
280	<	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 285 | 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		{
295	–	sides[1] = GT_EDGE;
296	–	/* If on plane 0-and on plane 1: lies on edge */
297	–	if(sides[0] == GT_EDGE)
298	–	{
299	–	*which = 1;
300	–	sides[2] = GT_INVALID;
301	–	return(GT_EDGE);
302	–	}
303	–	}
304	–	else if(d > 0)
305	–	{
376		sides[1] = GT_OUT;
377		sides[2] = GT_INVALID;
378		return(FALSE);
#	Line 312 | Line 382 | char sides[3];
382		/* Test next edge */
383		if(!NTH_BIT(*nset,2))
384		{
385	<
316	<	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	<
325	<	/* If on plane 0 and 2: lies on edge 0*/
326	<	if(sides[0] == GT_EDGE)
327	<	{
328	<	*which = 0;
329	<	return(GT_EDGE);
330	<	}
331	<	/* If on plane 1 and 2: lies on edge  2*/
332	<	if(sides[1] == GT_EDGE)
333	<	{
334	<	*which = 2;
335	<	return(GT_EDGE);
336	<	}
337	<	/* otherwise: on face 2 */
338	<	else
339	<	{
340	<	*which = 2;
341	<	return(GT_FACE);
342	<	}
392	>	sides[2] = GT_OUT;
393	>	return(FALSE);
394		}
344	–	else if(d > 0)
345	–	{
346	–	sides[2] = GT_OUT;
347	–	return(FALSE);
348	–	}
349	–	/* If on edge */
395		else
396		sides[2] = GT_INTERIOR;
352	–
353	–	/* If on plane 0 only: on face 0 */
354	–	if(sides[0] == GT_EDGE)
355	–	{
356	–	*which = 0;
357	–	return(GT_FACE);
358	–	}
359	–	/* If on plane 1 only: on face 1 */
360	–	if(sides[1] == GT_EDGE)
361	–	{
362	–	*which = 1;
363	–	return(GT_FACE);
364	–	}
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;
375	–	char *which;
407		{
408	<	float d;
408	>	double d;
409		FVECT n;
379	–	char sides[3];
410
411	<	/* First test if point coincides with any of the vertices */
382	<	if(EQUAL_VEC3(p,v0))
383	<	{
384	<	*which = 0;
385	<	return(GT_VERTEX);
386	<	}
387	<	if(EQUAL_VEC3(p,v1))
388	<	{
389	<	*which = 1;
390	<	return(GT_VERTEX);
391	<	}
392	<	if(EQUAL_VEC3(p,v2))
393	<	{
394	<	*which = 2;
395	<	return(GT_VERTEX);
396	<	}
397	<	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
403	<	if(d > 0)
404	<	return(FALSE);
405	<	else
406	<	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))
412	<	{
413	<	sides[1] = GT_EDGE;
414	<	/* If on plane 0-and on plane 1: lies on edge */
415	<	if(sides[0] == GT_EDGE)
416	<	{
417	<	*which = 1;
418	<	return(GT_VERTEX);
419	<	}
420	<	}
421	<	else if(d > 0)
421	>	if(d > 0.0)
422		return(FALSE);
423	–	else
424	–	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))
431	<	{
432	<	sides[2] = GT_EDGE;
433	<
434	<	/* If on plane 0 and 2: lies on edge 0*/
435	<	if(sides[0] == GT_EDGE)
436	<	{
437	<	*which = 0;
438	<	return(GT_VERTEX);
439	<	}
440	<	/* If on plane 1 and 2: lies on edge  2*/
441	<	if(sides[1] == GT_EDGE)
442	<	{
443	<	*which = 2;
444	<	return(GT_VERTEX);
445	<	}
446	<	/* otherwise: on face 2 */
447	<	else
448	<	{
449	<	return(GT_FACE);
450	<	}
451	<	}
452	<	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;
468	<	char which;
443	>	int below_plane[3],test;
444
445		SUM_3VEC3(avg,t0,t1,t2);
471	–	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 476 | 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,
485	<	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);
491	–	/* Remember if b[i] fell on one of the 3 defining planes */
492	–	if(test)
493	–	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,
503	<	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);
509	–	/* Remember if b[i] fell on one of the 3 defining planes */
510	–	if(test)
511	–	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,
521	<	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);
528	–	/* Remember if b[i] fell on one of the 3 defining planes */
529	–	if(test)
530	–	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);
536	–	/* Accept if all points lie on a triangle vertex/edge edge- accept*/
537	–	if(on_edge == 3)
538	–	return(TRUE);
498		/* Now check vertices in a against triangle b */
499		return(FALSE);
500		}
#	Line 543 | 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 557 | 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 570 | 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 583 | 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 599 | 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 613 | 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 627 | 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 643 | 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 656 | 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 669 | 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 682 | 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 695 | 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,
699	<	&(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,
703	<	&(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 748 | 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;
755	–	char *wptr;
712		{
713	<	FVECT p0,p1,p2,p,n;
714	<	char type,which;
715	<	double pd;
760	<
761	<	point_on_sphere(p0,v0,orig);
762	<	point_on_sphere(p1,v1,orig);
763	<	point_on_sphere(p2,v2,orig);
764	<	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)
773	<	*wptr = which;
774	<
775	<	return(type);
727	>	return(FALSE);
728		}
729
730
#	Line 831 | 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 879 | Line 869 | FVECT a0,a1,b0,b1;
869		return(FALSE);
870		return(TRUE);
871		}
872	+
873	+
874	+
875	+	/* Find the normalized barycentric coordinates of p relative to
876	+	* triangle v0,v1,v2. Return result in coord
877	+	*/
878	+	bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
879	+	double x1,y1,x2,y2,x3,y3;
880	+	double px,py;
881	+	double coord[3];
882	+	{
883	+	double a;
884	+
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] = ((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	+	bary_child(coord)
930	+	double coord[3];
931	+	{
932	+	int i;
933	+
934	+	if(coord[0] > 0.5)
935	+	{
936	+	/* update bary for child */
937	+	coord[0] = 2.0*coord[0]- 1.0;
938	+	coord[1] *= 2.0;
939	+	coord[2] *= 2.0;
940	+	return(0);
941	+	}
942	+	else
943	+	if(coord[1] > 0.5)
944	+	{
945	+	coord[0] *= 2.0;
946	+	coord[1] = 2.0*coord[1]- 1.0;
947	+	coord[2] *= 2.0;
948	+	return(1);
949	+	}
950	+	else
951	+	if(coord[2] > 0.5)
952	+	{
953	+	coord[0] *= 2.0;
954	+	coord[1] *= 2.0;
955	+	coord[2] = 2.0*coord[2]- 1.0;
956	+	return(2);
957	+	}
958	+	else
959	+	{
960	+	coord[0] = 1.0 - 2.0*coord[0];
961	+	coord[1] = 1.0 - 2.0*coord[1];
962	+	coord[2] = 1.0 - 2.0*coord[2];
963	+	return(3);
964	+	}
965	+	}
966	+
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	+	{
974	+
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	+	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	+	baryi_child(coord)
1198	+	BCOORD coord[3];
1199	+	{
1200	+	int i;
1201	+
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	+	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	+	int
1236	+	baryi_nth_child(coord,i)
1237	+	BCOORD coord[3];
1238	+	int i;
1239	+	{
1240	+
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	+
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	+	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	+	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	+
1454	+
1455	+
1456	+
1457	+
1458	+
1459

Diff Legend

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