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.7 by gwlarson, Tue Oct 6 18:16:53 1998 UTC vs.
Revision 3.13 by greg, Sat Feb 22 02:07:25 2003 UTC

#	Line 1 | Line 1
1	–	/* Copyright (c) 1998 Silicon Graphics, Inc. */
2	–
1		#ifndef lint
2	<	static char SCCSid[] = "$SunId$ SGI";
2	>	static const char       RCSid[] = "$Id$";
3		#endif
6	–
4		/*
5		*  sm_geom.c
6	+	*    some geometric utility routines
7		*/
8
9		#include "standard.h"
10		#include "sm_geom.h"
11
12	<	tri_centroid(v0,v1,v2,c)
13	<	FVECT v0,v1,v2,c;
12	>	/*
13	>	* int
14	>	* pt_in_cone(p,a,b,c)
15	>	*             : test if point p lies in cone defined by a,b,c and origin
16	>	* double p[3];              : point to test
17	>	* double a[3],b[3],c[3];    : points forming triangle
18	>	*
19	>	*  Assumes apex at origin, a,b,c are unit vectors defining the
20	>	*  triangle which the cone circumscribes. Assumes p is also normalized
21	>	*  Test is implemented as:
22	>	*             r =  (b-a)X(c-a)
23	>	*             in = (p.r) > (a.r)
24	>	*  The center of the cone is r, and corresponds to the triangle normal.
25	>	*  p.r is the proportional to the cosine of the angle between p and the
26	>	*  the cone center ray, and a.r to the radius of the cone. If the cosine
27	>	*  of the angle for p is greater than that for a, the angle between p
28	>	*  and r is smaller, and p lies in the cone.
29	>	*/
30	>	int
31	>	pt_in_cone(p,a,b,c)
32	>	double p[3],a[3],b[3],c[3];
33		{
34	<	/* Average three triangle vertices to give centroid: return in c */
35	<	c[0] = (v0[0] + v1[0] + v2[0])/3.0;
36	<	c[1] = (v0[1] + v1[1] + v2[1])/3.0;
37	<	c[2] = (v0[2] + v1[2] + v2[2])/3.0;
34	>	double r[3];
35	>	double pr,ar;
36	>	double ab[3],ac[3];
37	>
38	>	#ifdef DEBUG
39	>	#if DEBUG > 1
40	>	{
41	>	double l;
42	>	VSUB(ab,b,a);
43	>	normalize(ab);
44	>	VSUB(ac,c,a);
45	>	normalize(ac);
46	>	VCROSS(r,ab,ac);
47	>	l = normalize(r);
48	>	/* l = sin@ between ab,ac - if 0 vectors are colinear */
49	>	if( l <= COLINEAR_EPS)
50	>	{
51	>	eputs("pt in cone: null triangle:returning FALSE\n");
52	>	return(FALSE);
53	>	}
54		}
55	+	#endif
56	+	#endif
57
58	+	VSUB(ab,b,a);
59	+	VSUB(ac,c,a);
60	+	VCROSS(r,ab,ac);
61
62	<	int
63	<	vec3_equal(v1,v2)
64	<	FVECT v1,v2;
65	<	{
66	<	return(EQUAL(v1[0],v2[0]) && EQUAL(v1[1],v2[1])&& EQUAL(v1[2],v2[2]));
62	>	pr = DOT(p,r);
63	>	ar = DOT(a,r);
64	>	/* Need to check for equality for degeneracy of 4 points on circle */
65	>	if( pr > ar *( 1.0 + EQUALITY_EPS))
66	>	return(TRUE);
67	>	else
68	>	return(FALSE);
69		}
70
71	<
72	<	int
73	<	convex_angle(v0,v1,v2)
74	<	FVECT v0,v1,v2;
71	>	/*
72	>	* tri_centroid(v0,v1,v2,c)
73	>	*             : Average triangle vertices to give centroid: return in c
74	>	*FVECT v0,v1,v2,c; : triangle vertices(v0,v1,v2) and vector to hold result(c)
75	>	*/
76	>	tri_centroid(v0,v1,v2,c)
77	>	FVECT v0,v1,v2,c;
78		{
79	<	FVECT cp,cp01,cp12,v10,v02;
80	<	double dp;
81	<	/*
39	<	VSUB(v10,v1,v0);
40	<	VSUB(v02,v0,v2);
41	<	VCROSS(cp,v10,v02);
42	<	*/
43	<	/* test sign of (v0Xv1)X(v1Xv2). v1 */
44	<	VCROSS(cp01,v0,v1);
45	<	VCROSS(cp12,v1,v2);
46	<	VCROSS(cp,cp01,cp12);
47	<
48	<	dp = DOT(cp,v1);
49	<	if(ZERO(dp) || dp < 0.0)
50	<	return(FALSE);
51	<	return(TRUE);
79	>	c[0] = (v0[0] + v1[0] + v2[0])/3.0;
80	>	c[1] = (v0[1] + v1[1] + v2[1])/3.0;
81	>	c[2] = (v0[2] + v1[2] + v2[2])/3.0;
82		}
83
54	–	/* calculates the normal of a face contour using Newell's formula. e
84
85	<	a =  SUMi (yi - yi+1)(zi + zi+1)
86	<	b =  SUMi (zi - zi+1)(xi + xi+1)
87	<	c =  SUMi (xi - xi+1)(yi + yi+1)
85	>	/*
86	>	* double
87	>	* tri_normal(v0,v1,v2,n,norm) : Calculates the normal of a face contour using
88	>	*                            Newell's formula.
89	>	* FVECT v0,v1,v2,n;  : Triangle vertices(v0,v1,v2) and vector for result(n)
90	>	* int norm;          : If true result is normalized
91	>	*
92	>	* Triangle normal is calculated using the following:
93	>	*    A =  SUMi (yi - yi+1)(zi + zi+1);
94	>	*    B =  SUMi (zi - zi+1)(xi + xi+1)
95	>	*    C =  SUMi (xi - xi+1)(yi + yi+1)
96		*/
97		double
98		tri_normal(v0,v1,v2,n,norm)
#	Line 67 | Line 104 | int norm;
104		n[0] = (v0[2] + v1[2]) * (v0[1] - v1[1]) +
105		(v1[2] + v2[2]) * (v1[1] - v2[1])  +
106		(v2[2] + v0[2]) * (v2[1] - v0[1]);
70	–
107		n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
108		(v1[2] - v2[2]) * (v1[0] + v2[0]) +
109	<	(v2[2] - v0[2]) * (v2[0] + v0[0]);
74	<
109	>	(v2[2] - v0[2]) * (v2[0] + v0[0]);
110		n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
111		(v1[1] + v2[1]) * (v1[0] - v2[0]) +
112		(v2[1] + v0[1]) * (v2[0] - v0[0]);
78	–
113		if(!norm)
114		return(0);
81	–
82	–
115		mag = normalize(n);
84	–
116		return(mag);
117		}
118
119	<
119	>	/*
120	>	* tri_plane_equation(v0,v1,v2,peqptr,norm)
121	>	*             : Calculates the plane equation (A,B,C,D) for triangle
122	>	*               v0,v1,v2 ( Ax + By + Cz = D)
123	>	* FVECT v0,v1,v2;   : Triangle vertices
124	>	* FPEQ *peqptr;     : ptr to structure to hold result
125	>	*  int norm;        : if TRUE, return unit normal
126	>	*/
127		tri_plane_equation(v0,v1,v2,peqptr,norm)
128	<	FVECT v0,v1,v2;
128	>	FVECT v0,v1,v2;
129		FPEQ *peqptr;
130		int norm;
131		{
132		tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
95	–
133		FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
134		}
135
136	<	/* From quad_edge-code */
136	>	/*
137	>	* int
138	>	* intersect_ray_plane(orig,dir,peq,pd,r)
139	>	*             : Intersects ray (orig,dir) with plane (peq). Returns TRUE
140	>	*             if intersection occurs. If r!=NULL, sets with resulting i
141	>	*             intersection point, and pd is set with parametric value of the
142	>	*             intersection.
143	>	* FVECT orig,dir;      : vectors defining the ray
144	>	* FPEQ peq;            : plane equation
145	>	* double *pd;          : holds resulting parametric intersection point
146	>	* FVECT r;             : holds resulting intersection point
147	>	*
148	>	*    Plane is Ax + By + Cz +D = 0:
149	>	*    A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
150	>	*     t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
151	>	*      line is  l = p1 + (p2-p1)t
152	>	*    Solve for t
153	>	*/
154		int
101	–	point_in_circle_thru_origin(p,p0,p1)
102	–	FVECT p;
103	–	FVECT p0,p1;
104	–	{
105	–
106	–	double dp0,dp1;
107	–	double dp,det;
108	–
109	–	dp0 = DOT_VEC2(p0,p0);
110	–	dp1 = DOT_VEC2(p1,p1);
111	–	dp  = DOT_VEC2(p,p);
112	–
113	–	det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
114	–
115	–	return (det > 0);
116	–	}
117	–
118	–
119	–
120	–	point_on_sphere(ps,p,c)
121	–	FVECT ps,p,c;
122	–	{
123	–	VSUB(ps,p,c);
124	–	normalize(ps);
125	–	}
126	–
127	–
128	–	/* returns TRUE if ray from origin in direction v intersects plane defined
129	–	* by normal plane_n, and plane_d. If plane is not parallel- returns
130	–	* intersection point if r != NULL. If tptr!= NULL returns value of
131	–	* t, if parallel, returns t=FHUGE
132	–	*/
133	–	int
134	–	intersect_vector_plane(v,peq,tptr,r)
135	–	FVECT v;
136	–	FPEQ peq;
137	–	double *tptr;
138	–	FVECT r;
139	–	{
140	–	double t,d;
141	–	int hit;
142	–	/*
143	–	Plane is Ax + By + Cz +D = 0:
144	–	plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
145	–	*/
146	–
147	–	/* line is  l = p1 + (p2-p1)t, p1=origin */
148	–
149	–	/* Solve for t: */
150	–	d = -(DOT(FP_N(peq),v));
151	–	if(ZERO(d))
152	–	{
153	–	t = FHUGE;
154	–	hit = 0;
155	–	}
156	–	else
157	–	{
158	–	t =  FP_D(peq)/d;
159	–	if(t < 0 )
160	–	hit = 0;
161	–	else
162	–	hit = 1;
163	–	if(r)
164	–	{
165	–	r[0] = v[0]*t;
166	–	r[1] = v[1]*t;
167	–	r[2] = v[2]*t;
168	–	}
169	–	}
170	–	if(tptr)
171	–	*tptr = t;
172	–	return(hit);
173	–	}
174	–
175	–	int
155		intersect_ray_plane(orig,dir,peq,pd,r)
156		FVECT orig,dir;
157		FPEQ peq;
158		double *pd;
159		FVECT r;
160		{
161	<	double t;
161	>	double t,d;
162		int hit;
184	–	/*
185	–	Plane is Ax + By + Cz +D = 0:
186	–	plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
187	–	*/
188	–	/*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
189	–	t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
190	–	line is  l = p1 + (p2-p1)t
191	–	*/
192	–	/* Solve for t: */
193	–	t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/(DOT(FP_N(peq),dir));
194	–	if(t < 0)
195	–	hit = 0;
196	–	else
197	–	hit = 1;
163
164	<	if(r)
165	<	VSUM(r,orig,dir,t);
164	>	d = DOT(FP_N(peq),dir);
165	>	if(ZERO(d))
166	>	return(0);
167	>	t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
168
169	<	if(pd)
203	<	*pd = t;
204	<	return(hit);
205	<	}
206	<
207	<
208	<	int
209	<	intersect_ray_oplane(orig,dir,n,pd,r)
210	<	FVECT orig,dir;
211	<	FVECT n;
212	<	double *pd;
213	<	FVECT r;
214	<	{
215	<	double t;
216	<	int hit;
217	<	/*
218	<	Plane is Ax + By + Cz +D = 0:
219	<	plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
220	<	*/
221	<	/*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
222	<	t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
223	<	line is  l = p1 + (p2-p1)t
224	<	*/
225	<	/* Solve for t: */
226	<	t =  -(DOT(n,orig))/(DOT(n,dir));
227	<	if(t < 0)
169	>	if(t < 0)
170		hit = 0;
171		else
172		hit = 1;
231	–
173		if(r)
174		VSUM(r,orig,dir,t);
175
#	Line 237 | Line 178 | intersect_ray_oplane(orig,dir,n,pd,r)
178		return(hit);
179		}
180
181	<
182	<	int
183	<	intersect_edge_plane(e0,e1,peq,pd,r)
184	<	FVECT e0,e1;
185	<	FPEQ peq;
186	<	double *pd;
187	<	FVECT r;
181	>	/*
182	>	* double
183	>	* point_on_sphere(ps,p,c) : normalize p relative to sphere with center c
184	>	* FVECT ps,p,c;       : ps Holds result vector,p is the original point,
185	>	*                       and c is the sphere center
186	>	*/
187	>	double
188	>	point_on_sphere(ps,p,c)
189	>	FVECT ps,p,c;
190		{
191	<	double t;
192	<	int hit;
193	<	FVECT d;
194	<	/*
195	<	Plane is Ax + By + Cz +D = 0:
253	<	plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
254	<	*/
255	<	/*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
256	<	t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
257	<	line is  l = p1 + (p2-p1)t
258	<	*/
259	<	/* Solve for t: */
260	<	VSUB(d,e1,e0);
261	<	t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
262	<	if(t < 0)
263	<	hit = 0;
264	<	else
265	<	hit = 1;
266	<
267	<	VSUM(r,e0,d,t);
268	<
269	<	if(pd)
270	<	*pd = t;
271	<	return(hit);
191	>	double d;
192	>
193	>	VSUB(ps,p,c);
194	>	d = normalize(ps);
195	>	return(d);
196		}
197
198	<
198	>	/*
199	>	* int
200	>	* point_in_stri(v0,v1,v2,p) : Return TRUE if p is in pyramid defined by
201	>	*                            tri v0,v1,v2 and origin
202	>	* FVECT v0,v1,v2,p;     :Triangle vertices(v0,v1,v2) and point in question(p)
203	>	*
204	>	*   Tests orientation of p relative to each edge (v0v1,v1v2,v2v0), if it is
205	>	*   inside of all 3 edges, returns TRUE, else FALSE.
206	>	*/
207		int
276	–	point_in_cone(p,p0,p1,p2)
277	–	FVECT p;
278	–	FVECT p0,p1,p2;
279	–	{
280	–	FVECT np,x_axis,y_axis;
281	–	double d1,d2;
282	–	FPEQ peq;
283	–
284	–	/* Find the equation of the circle defined by the intersection
285	–	of the cone with the plane defined by p1,p2,p3- project p into
286	–	that plane and do an in-circle test in the plane
287	–	*/
288	–
289	–	/* find the equation of the plane defined by p1-p3 */
290	–	tri_plane_equation(p0,p1,p2,&peq,FALSE);
291	–
292	–	/* define a coordinate system on the plane: the x axis is in
293	–	the direction of np2-np1, and the y axis is calculated from
294	–	n cross x-axis
295	–	*/
296	–	/* Project p onto the plane */
297	–	/* NOTE: check this: does sideness check?*/
298	–	if(!intersect_vector_plane(p,peq,NULL,np))
299	–	return(FALSE);
300	–
301	–	/* create coordinate system on  plane: p2-p1 defines the x_axis*/
302	–	VSUB(x_axis,p1,p0);
303	–	normalize(x_axis);
304	–	/* The y axis is  */
305	–	VCROSS(y_axis,FP_N(peq),x_axis);
306	–	normalize(y_axis);
307	–
308	–	VSUB(p1,p1,p0);
309	–	VSUB(p2,p2,p0);
310	–	VSUB(np,np,p0);
311	–
312	–	p1[0] = VLEN(p1);
313	–	p1[1] = 0;
314	–
315	–	d1 = DOT(p2,x_axis);
316	–	d2 = DOT(p2,y_axis);
317	–	p2[0] = d1;
318	–	p2[1] = d2;
319	–
320	–	d1 = DOT(np,x_axis);
321	–	d2 = DOT(np,y_axis);
322	–	np[0] = d1;
323	–	np[1] = d2;
324	–
325	–	/* perform the in-circle test in the new coordinate system */
326	–	return(point_in_circle_thru_origin(np,p1,p2));
327	–	}
328	–
329	–	int
330	–	point_set_in_stri(v0,v1,v2,p,n,nset,sides)
331	–	FVECT v0,v1,v2,p;
332	–	FVECT n[3];
333	–	int *nset;
334	–	int sides[3];
335	–
336	–	{
337	–	double d;
338	–	/* Find the normal to the triangle ORIGIN:v0:v1 */
339	–	if(!NTH_BIT(*nset,0))
340	–	{
341	–	VCROSS(n[0],v1,v0);
342	–	SET_NTH_BIT(*nset,0);
343	–	}
344	–	/* Test the point for sidedness */
345	–	d  = DOT(n[0],p);
346	–
347	–	if(d > 0.0)
348	–	{
349	–	sides[0] =  GT_OUT;
350	–	sides[1] = sides[2] = GT_INVALID;
351	–	return(FALSE);
352	–	}
353	–	else
354	–	sides[0] = GT_INTERIOR;
355	–
356	–	/* Test next edge */
357	–	if(!NTH_BIT(*nset,1))
358	–	{
359	–	VCROSS(n[1],v2,v1);
360	–	SET_NTH_BIT(*nset,1);
361	–	}
362	–	/* Test the point for sidedness */
363	–	d  = DOT(n[1],p);
364	–	if(d > 0.0)
365	–	{
366	–	sides[1] = GT_OUT;
367	–	sides[2] = GT_INVALID;
368	–	return(FALSE);
369	–	}
370	–	else
371	–	sides[1] = GT_INTERIOR;
372	–	/* Test next edge */
373	–	if(!NTH_BIT(*nset,2))
374	–	{
375	–	VCROSS(n[2],v0,v2);
376	–	SET_NTH_BIT(*nset,2);
377	–	}
378	–	/* Test the point for sidedness */
379	–	d  = DOT(n[2],p);
380	–	if(d > 0.0)
381	–	{
382	–	sides[2] = GT_OUT;
383	–	return(FALSE);
384	–	}
385	–	else
386	–	sides[2] = GT_INTERIOR;
387	–	/* Must be interior to the pyramid */
388	–	return(GT_INTERIOR);
389	–	}
390	–
391	–
392	–
393	–
394	–	int
208		point_in_stri(v0,v1,v2,p)
209		FVECT v0,v1,v2,p;
210		{
211		double d;
212		FVECT n;
213
214	<	VCROSS(n,v1,v0);
214	>	VCROSS(n,v0,v1);
215		/* Test the point for sidedness */
216		d  = DOT(n,p);
217		if(d > 0.0)
218		return(FALSE);
406	–
219		/* Test next edge */
220	<	VCROSS(n,v2,v1);
220	>	VCROSS(n,v1,v2);
221		/* Test the point for sidedness */
222		d  = DOT(n,p);
223		if(d > 0.0)
224		return(FALSE);
413	–
225		/* Test next edge */
226	<	VCROSS(n,v0,v2);
226	>	VCROSS(n,v2,v0);
227		/* Test the point for sidedness */
228		d  = DOT(n,p);
229		if(d > 0.0)
230		return(FALSE);
231		/* Must be interior to the pyramid */
232	<	return(GT_INTERIOR);
232	>	return(TRUE);
233		}
234
235	+	/*
236	+	* int
237	+	* ray_intersect_tri(orig,dir,v0,v1,v2,pt)
238	+	*    : test if ray orig-dir intersects triangle v0v1v2, result in pt
239	+	* FVECT orig,dir;   : Vectors defining ray origin and direction
240	+	* FVECT v0,v1,v2;   : Triangle vertices
241	+	* FVECT pt;         : Intersection point (if any)
242	+	*/
243		int
425	–	vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
426	–	FVECT t0,t1,t2,p0,p1,p2;
427	–	int *nset;
428	–	FVECT n[3];
429	–	FVECT avg;
430	–	int pt_sides[3][3];
431	–
432	–	{
433	–	int below_plane[3],test;
434	–
435	–	SUM_3VEC3(avg,t0,t1,t2);
436	–	*nset = 0;
437	–	/* Test vertex v[i] against triangle j*/
438	–	/* Check if v[i] lies below plane defined by avg of 3 vectors
439	–	defining triangle
440	–	*/
441	–
442	–	/* test point 0 */
443	–	if(DOT(avg,p0) < 0.0)
444	–	below_plane[0] = 1;
445	–	else
446	–	below_plane[0] = 0;
447	–	/* Test if b[i] lies in or on triangle a */
448	–	test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
449	–	/* If pts[i] is interior: done */
450	–	if(!below_plane[0])
451	–	{
452	–	if(test == GT_INTERIOR)
453	–	return(TRUE);
454	–	}
455	–	/* Now test point 1*/
456	–
457	–	if(DOT(avg,p1) < 0.0)
458	–	below_plane[1] = 1;
459	–	else
460	–	below_plane[1]=0;
461	–	/* Test if b[i] lies in or on triangle a */
462	–	test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
463	–	/* If pts[i] is interior: done */
464	–	if(!below_plane[1])
465	–	{
466	–	if(test == GT_INTERIOR)
467	–	return(TRUE);
468	–	}
469	–
470	–	/* Now test point 2 */
471	–	if(DOT(avg,p2) < 0.0)
472	–	below_plane[2] = 1;
473	–	else
474	–	below_plane[2] = 0;
475	–	/* Test if b[i] lies in or on triangle a */
476	–	test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
477	–
478	–	/* If pts[i] is interior: done */
479	–	if(!below_plane[2])
480	–	{
481	–	if(test == GT_INTERIOR)
482	–	return(TRUE);
483	–	}
484	–
485	–	/* If all three points below separating plane: trivial reject */
486	–	if(below_plane[0] && below_plane[1] && below_plane[2])
487	–	return(FALSE);
488	–	/* Now check vertices in a against triangle b */
489	–	return(FALSE);
490	–	}
491	–
492	–
493	–	set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
494	–	FVECT t0,t1,t2,p0,p1,p2;
495	–	int test[3];
496	–	int sides[3][3];
497	–	int nset;
498	–	FVECT n[3];
499	–	{
500	–	int t;
501	–	double d;
502	–
503	–
504	–	/* p=0 */
505	–	test[0] = 0;
506	–	if(sides[0][0] == GT_INVALID)
507	–	{
508	–	if(!NTH_BIT(nset,0))
509	–	VCROSS(n[0],t1,t0);
510	–	/* Test the point for sidedness */
511	–	d  = DOT(n[0],p0);
512	–	if(d >= 0.0)
513	–	SET_NTH_BIT(test[0],0);
514	–	}
515	–	else
516	–	if(sides[0][0] != GT_INTERIOR)
517	–	SET_NTH_BIT(test[0],0);
518	–
519	–	if(sides[0][1] == GT_INVALID)
520	–	{
521	–	if(!NTH_BIT(nset,1))
522	–	VCROSS(n[1],t2,t1);
523	–	/* Test the point for sidedness */
524	–	d  = DOT(n[1],p0);
525	–	if(d >= 0.0)
526	–	SET_NTH_BIT(test[0],1);
527	–	}
528	–	else
529	–	if(sides[0][1] != GT_INTERIOR)
530	–	SET_NTH_BIT(test[0],1);
531	–
532	–	if(sides[0][2] == GT_INVALID)
533	–	{
534	–	if(!NTH_BIT(nset,2))
535	–	VCROSS(n[2],t0,t2);
536	–	/* Test the point for sidedness */
537	–	d  = DOT(n[2],p0);
538	–	if(d >= 0.0)
539	–	SET_NTH_BIT(test[0],2);
540	–	}
541	–	else
542	–	if(sides[0][2] != GT_INTERIOR)
543	–	SET_NTH_BIT(test[0],2);
544	–
545	–	/* p=1 */
546	–	test[1] = 0;
547	–	/* t=0*/
548	–	if(sides[1][0] == GT_INVALID)
549	–	{
550	–	if(!NTH_BIT(nset,0))
551	–	VCROSS(n[0],t1,t0);
552	–	/* Test the point for sidedness */
553	–	d  = DOT(n[0],p1);
554	–	if(d >= 0.0)
555	–	SET_NTH_BIT(test[1],0);
556	–	}
557	–	else
558	–	if(sides[1][0] != GT_INTERIOR)
559	–	SET_NTH_BIT(test[1],0);
560	–
561	–	/* t=1 */
562	–	if(sides[1][1] == GT_INVALID)
563	–	{
564	–	if(!NTH_BIT(nset,1))
565	–	VCROSS(n[1],t2,t1);
566	–	/* Test the point for sidedness */
567	–	d  = DOT(n[1],p1);
568	–	if(d >= 0.0)
569	–	SET_NTH_BIT(test[1],1);
570	–	}
571	–	else
572	–	if(sides[1][1] != GT_INTERIOR)
573	–	SET_NTH_BIT(test[1],1);
574	–
575	–	/* t=2 */
576	–	if(sides[1][2] == GT_INVALID)
577	–	{
578	–	if(!NTH_BIT(nset,2))
579	–	VCROSS(n[2],t0,t2);
580	–	/* Test the point for sidedness */
581	–	d  = DOT(n[2],p1);
582	–	if(d >= 0.0)
583	–	SET_NTH_BIT(test[1],2);
584	–	}
585	–	else
586	–	if(sides[1][2] != GT_INTERIOR)
587	–	SET_NTH_BIT(test[1],2);
588	–
589	–	/* p=2 */
590	–	test[2] = 0;
591	–	/* t = 0 */
592	–	if(sides[2][0] == GT_INVALID)
593	–	{
594	–	if(!NTH_BIT(nset,0))
595	–	VCROSS(n[0],t1,t0);
596	–	/* Test the point for sidedness */
597	–	d  = DOT(n[0],p2);
598	–	if(d >= 0.0)
599	–	SET_NTH_BIT(test[2],0);
600	–	}
601	–	else
602	–	if(sides[2][0] != GT_INTERIOR)
603	–	SET_NTH_BIT(test[2],0);
604	–	/* t=1 */
605	–	if(sides[2][1] == GT_INVALID)
606	–	{
607	–	if(!NTH_BIT(nset,1))
608	–	VCROSS(n[1],t2,t1);
609	–	/* Test the point for sidedness */
610	–	d  = DOT(n[1],p2);
611	–	if(d >= 0.0)
612	–	SET_NTH_BIT(test[2],1);
613	–	}
614	–	else
615	–	if(sides[2][1] != GT_INTERIOR)
616	–	SET_NTH_BIT(test[2],1);
617	–	/* t=2 */
618	–	if(sides[2][2] == GT_INVALID)
619	–	{
620	–	if(!NTH_BIT(nset,2))
621	–	VCROSS(n[2],t0,t2);
622	–	/* Test the point for sidedness */
623	–	d  = DOT(n[2],p2);
624	–	if(d >= 0.0)
625	–	SET_NTH_BIT(test[2],2);
626	–	}
627	–	else
628	–	if(sides[2][2] != GT_INTERIOR)
629	–	SET_NTH_BIT(test[2],2);
630	–	}
631	–
632	–
633	–	int
634	–	stri_intersect(a1,a2,a3,b1,b2,b3)
635	–	FVECT a1,a2,a3,b1,b2,b3;
636	–	{
637	–	int which,test,n_set[2];
638	–	int sides[2][3][3],i,j,inext,jnext;
639	–	int tests[2][3];
640	–	FVECT n[2][3],p,avg[2];
641	–
642	–	/* Test the vertices of triangle a against the pyramid formed by triangle
643	–	b and the origin. If any vertex of a is interior to triangle b, or
644	–	if all 3 vertices of a are ON the edges of b,return TRUE. Remember
645	–	the results of the edge normal and sidedness tests for later.
646	–	*/
647	–	if(vertices_in_stri(a1,a2,a3,b1,b2,b3,&(n_set[0]),n[0],avg[0],sides[1]))
648	–	return(TRUE);
649	–
650	–	if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
651	–	return(TRUE);
652	–
653	–
654	–	set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
655	–	if(tests[0][0]&tests[0][1]&tests[0][2])
656	–	return(FALSE);
657	–
658	–	set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
659	–	if(tests[1][0]&tests[1][1]&tests[1][2])
660	–	return(FALSE);
661	–
662	–	for(j=0; j < 3;j++)
663	–	{
664	–	jnext = (j+1)%3;
665	–	/* IF edge b doesnt cross any great circles of a, punt */
666	–	if(tests[1][j] & tests[1][jnext])
667	–	continue;
668	–	for(i=0;i<3;i++)
669	–	{
670	–	inext = (i+1)%3;
671	–	/* IF edge a doesnt cross any great circles of b, punt */
672	–	if(tests[0][i] & tests[0][inext])
673	–	continue;
674	–	/* Now find the great circles that cross and test */
675	–	if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
676	–	&& (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
677	–	{
678	–	VCROSS(p,n[0][i],n[1][j]);
679	–
680	–	/* If zero cp= done */
681	–	if(ZERO_VEC3(p))
682	–	continue;
683	–	/* check above both planes */
684	–	if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
685	–	{
686	–	NEGATE_VEC3(p);
687	–	if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
688	–	continue;
689	–	}
690	–	return(TRUE);
691	–	}
692	–	}
693	–	}
694	–	return(FALSE);
695	–	}
696	–
697	–	int
244		ray_intersect_tri(orig,dir,v0,v1,v2,pt)
245		FVECT orig,dir;
246		FVECT v0,v1,v2;
#	Line 717 | Line 263 | FVECT pt;
263		return(FALSE);
264		}
265
266	<
266	>	/*
267	>	* calculate_view_frustum(vp,hv,vv,horiz,vert,near,far,fnear,ffar)
268	>	*    : Calculate vertices defining front and rear clip rectangles of
269	>	*      view frustum defined by vp,hv,vv,horiz,vert,near, and far and
270	>	*      return in fnear and ffar.
271	>	* FVECT vp,hv,vv;        : Viewpoint(vp),hv and vv are the horizontal and
272	>	*                          vertical vectors in the view frame-magnitude is
273	>	*                          the dimension of the front frustum face at z =1
274	>	* double horiz,vert,near,far; : View angle horizontal and vertical(horiz,vert)
275	>	*                              and distance to the near,far clipping planes
276	>	* FVECT fnear[4],ffar[4];     : holds results
277	>	*
278	>	*/
279		calculate_view_frustum(vp,hv,vv,horiz,vert,near,far,fnear,ffar)
280		FVECT vp,hv,vv;
281		double horiz,vert,near,far;
#	Line 727 | Line 285 | FVECT fnear[4],ffar[4];
285		FVECT t,nhv,nvv,ndv;
286		double w2,h2;
287		/* Calculate the x and y dimensions of the near face */
730	–	/* hv and vv are the horizontal and vertical vectors in the
731	–	view frame-the magnitude is the dimension of the front frustum
732	–	face at z =1
733	–	*/
288		VCOPY(nhv,hv);
289		VCOPY(nvv,vv);
290		w2 = normalize(nhv);
#	Line 773 | Line 327 | FVECT fnear[4],ffar[4];
327		ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
328		}
329
776	–	int
777	–	max_index(v,r)
778	–	FVECT v;
779	–	double *r;
780	–	{
781	–	double p[3];
782	–	int i;
330
331	<	p[0] = fabs(v[0]);
332	<	p[1] = fabs(v[1]);
333	<	p[2] = fabs(v[2]);
334	<	i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);
335	<	if(r)
336	<	*r = p[i];
337	<	return(i);
791	<	}
792	<
793	<	int
794	<	closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
795	<	FVECT p0,p1,p2,p;
796	<	int p0id,p1id,p2id;
797	<	{
798	<	double d,d1;
799	<	int i;
800	<
801	<	d =  DIST_SQ(p,p0);
802	<	d1 = DIST_SQ(p,p1);
803	<	if(d < d1)
804	<	{
805	<	d1 = DIST_SQ(p,p2);
806	<	i = (d1 < d)?p2id:p0id;
807	<	}
808	<	else
809	<	{
810	<	d = DIST_SQ(p,p2);
811	<	i = (d < d1)? p2id:p1id;
812	<	}
813	<	return(i);
814	<	}
815	<
816	<
817	<	int
818	<	sedge_intersect(a0,a1,b0,b1)
819	<	FVECT a0,a1,b0,b1;
820	<	{
821	<	FVECT na,nb,avga,avgb,p;
822	<	double d;
823	<	int sb0,sb1,sa0,sa1;
824	<
825	<	/* First test if edge b straddles great circle of a */
826	<	VCROSS(na,a0,a1);
827	<	d = DOT(na,b0);
828	<	sb0 = ZERO(d)?0:(d<0)? -1: 1;
829	<	d = DOT(na,b1);
830	<	sb1 = ZERO(d)?0:(d<0)? -1: 1;
831	<	/* edge b entirely on one side of great circle a: edges cannot intersect*/
832	<	if(sb0*sb1 > 0)
833	<	return(FALSE);
834	<	/* test if edge a straddles great circle of b */
835	<	VCROSS(nb,b0,b1);
836	<	d = DOT(nb,a0);
837	<	sa0 = ZERO(d)?0:(d<0)? -1: 1;
838	<	d = DOT(nb,a1);
839	<	sa1 = ZERO(d)?0:(d<0)? -1: 1;
840	<	/* edge a entirely on one side of great circle b: edges cannot intersect*/
841	<	if(sa0*sa1 > 0)
842	<	return(FALSE);
843	<
844	<	/* Find one of intersection points of the great circles */
845	<	VCROSS(p,na,nb);
846	<	/* If they lie on same great circle: call an intersection */
847	<	if(ZERO_VEC3(p))
848	<	return(TRUE);
849	<
850	<	VADD(avga,a0,a1);
851	<	VADD(avgb,b0,b1);
852	<	if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
853	<	{
854	<	NEGATE_VEC3(p);
855	<	if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
856	<	return(FALSE);
857	<	}
858	<	if((!sb0 || !sb1) && (!sa0 || !sa1))
859	<	return(FALSE);
860	<	return(TRUE);
861	<	}
862	<
863	<
864	<
865	<	/* Find the normalized barycentric coordinates of p relative to
866	<	* triangle v0,v1,v2. Return result in coord
331	>	/*
332	>	* bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
333	>	*  : Find the normalized barycentric coordinates of p relative to
334	>	*    triangle v0,v1,v2. Return result in coord
335	>	* double x1,y1,x2,y2,x3,y3; : defines triangle vertices 1,2,3
336	>	* double px,py;             : coordinates of pt
337	>	* double coord[3];          : result
338		*/
339		bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
340		double x1,y1,x2,y2,x3,y3;
#	Line 879 | Line 350 | double coord[3];
350
351		}
352
882	–	bary_ith_child(coord,i)
883	–	double coord[3];
884	–	int i;
885	–	{
886	–
887	–	switch(i){
888	–	case 0:
889	–	/* update bary for child */
890	–	coord[0] = 2.0*coord[0]- 1.0;
891	–	coord[1] *= 2.0;
892	–	coord[2] *= 2.0;
893	–	break;
894	–	case 1:
895	–	coord[0] *= 2.0;
896	–	coord[1] = 2.0*coord[1]- 1.0;
897	–	coord[2] *= 2.0;
898	–	break;
899	–	case 2:
900	–	coord[0] *= 2.0;
901	–	coord[1] *= 2.0;
902	–	coord[2] = 2.0*coord[2]- 1.0;
903	–	break;
904	–	case 3:
905	–	coord[0] = 1.0 - 2.0*coord[0];
906	–	coord[1] = 1.0 - 2.0*coord[1];
907	–	coord[2] = 1.0 - 2.0*coord[2];
908	–	break;
909	–	#ifdef DEBUG
910	–	default:
911	–	eputs("bary_ith_child():Invalid child\n");
912	–	break;
913	–	#endif
914	–	}
915	–	}
353
354
918	–	int
919	–	bary_child(coord)
920	–	double coord[3];
921	–	{
922	–	int i;
355
924	–	if(coord[0] > 0.5)
925	–	{
926	–	/* update bary for child */
927	–	coord[0] = 2.0*coord[0]- 1.0;
928	–	coord[1] *= 2.0;
929	–	coord[2] *= 2.0;
930	–	return(0);
931	–	}
932	–	else
933	–	if(coord[1] > 0.5)
934	–	{
935	–	coord[0] *= 2.0;
936	–	coord[1] = 2.0*coord[1]- 1.0;
937	–	coord[2] *= 2.0;
938	–	return(1);
939	–	}
940	–	else
941	–	if(coord[2] > 0.5)
942	–	{
943	–	coord[0] *= 2.0;
944	–	coord[1] *= 2.0;
945	–	coord[2] = 2.0*coord[2]- 1.0;
946	–	return(2);
947	–	}
948	–	else
949	–	{
950	–	coord[0] = 1.0 - 2.0*coord[0];
951	–	coord[1] = 1.0 - 2.0*coord[1];
952	–	coord[2] = 1.0 - 2.0*coord[2];
953	–	return(3);
954	–	}
955	–	}
956	–
957	–	/* Coord was the ith child of the parent: set the coordinate
958	–	relative to the parent
959	–	*/
960	–	bary_parent(coord,i)
961	–	double coord[3];
962	–	int i;
963	–	{
964	–
965	–	switch(i) {
966	–	case 0:
967	–	/* update bary for child */
968	–	coord[0] = coord[0]*0.5 + 0.5;
969	–	coord[1] *= 0.5;
970	–	coord[2] *= 0.5;
971	–	break;
972	–	case 1:
973	–	coord[0] *= 0.5;
974	–	coord[1]  = coord[1]*0.5 + 0.5;
975	–	coord[2] *= 0.5;
976	–	break;
977	–
978	–	case 2:
979	–	coord[0] *= 0.5;
980	–	coord[1] *= 0.5;
981	–	coord[2] = coord[2]*0.5 + 0.5;
982	–	break;
983	–
984	–	case 3:
985	–	coord[0] = 0.5 - 0.5*coord[0];
986	–	coord[1] = 0.5 - 0.5*coord[1];
987	–	coord[2] = 0.5 - 0.5*coord[2];
988	–	break;
989	–	#ifdef DEBUG
990	–	default:
991	–	eputs("bary_parent():Invalid child\n");
992	–	break;
993	–	#endif
994	–	}
995	–	}
996	–
997	–	bary_from_child(coord,child,next)
998	–	double coord[3];
999	–	int child,next;
1000	–	{
1001	–	#ifdef DEBUG
1002	–	if(child <0 || child > 3)
1003	–	{
1004	–	eputs("bary_from_child():Invalid child\n");
1005	–	return;
1006	–	}
1007	–	if(next <0 || next > 3)
1008	–	{
1009	–	eputs("bary_from_child():Invalid next\n");
1010	–	return;
1011	–	}
1012	–	#endif
1013	–	if(next == child)
1014	–	return;
1015	–
1016	–	switch(child){
1017	–	case 0:
1018	–	switch(next){
1019	–	case 1:
1020	–	coord[0] += 1.0;
1021	–	coord[1] -= 1.0;
1022	–	break;
1023	–	case 2:
1024	–	coord[0] += 1.0;
1025	–	coord[2] -= 1.0;
1026	–	break;
1027	–	case 3:
1028	–	coord[0] *= -1.0;
1029	–	coord[1] = 1 - coord[1];
1030	–	coord[2] = 1 - coord[2];
1031	–	break;
1032	–
1033	–	}
1034	–	break;
1035	–	case 1:
1036	–	switch(next){
1037	–	case 0:
1038	–	coord[0] -= 1.0;
1039	–	coord[1] += 1.0;
1040	–	break;
1041	–	case 2:
1042	–	coord[1] += 1.0;
1043	–	coord[2] -= 1.0;
1044	–	break;
1045	–	case 3:
1046	–	coord[0] = 1 - coord[0];
1047	–	coord[1] *= -1.0;
1048	–	coord[2] = 1 - coord[2];
1049	–	break;
1050	–	}
1051	–	break;
1052	–	case 2:
1053	–	switch(next){
1054	–	case 0:
1055	–	coord[0] -= 1.0;
1056	–	coord[2] += 1.0;
1057	–	break;
1058	–	case 1:
1059	–	coord[1] -= 1.0;
1060	–	coord[2] += 1.0;
1061	–	break;
1062	–	case 3:
1063	–	coord[0] = 1 - coord[0];
1064	–	coord[1] = 1 - coord[1];
1065	–	coord[2] *= -1.0;
1066	–	break;
1067	–	}
1068	–	break;
1069	–	case 3:
1070	–	switch(next){
1071	–	case 0:
1072	–	coord[0] *= -1.0;
1073	–	coord[1] = 1 - coord[1];
1074	–	coord[2] = 1 - coord[2];
1075	–	break;
1076	–	case 1:
1077	–	coord[0] = 1 - coord[0];
1078	–	coord[1] *= -1.0;
1079	–	coord[2] = 1 - coord[2];
1080	–	break;
1081	–	case 2:
1082	–	coord[0] = 1 - coord[0];
1083	–	coord[1] = 1 - coord[1];
1084	–	coord[2] *= -1.0;
1085	–	break;
1086	–	}
1087	–	break;
1088	–	}
1089	–	}
1090	–
1091	–
1092	–	baryi_parent(coord,i)
1093	–	BCOORD coord[3];
1094	–	int i;
1095	–	{
1096	–
1097	–	switch(i) {
1098	–	case 0:
1099	–	/* update bary for child */
1100	–	coord[0] = (coord[0] >> 1) + MAXBCOORD2;
1101	–	coord[1] >>= 1;
1102	–	coord[2] >>= 1;
1103	–	break;
1104	–	case 1:
1105	–	coord[0] >>= 1;
1106	–	coord[1]  = (coord[1] >> 1) + MAXBCOORD2;
1107	–	coord[2] >>= 1;
1108	–	break;
1109	–
1110	–	case 2:
1111	–	coord[0] >>= 1;
1112	–	coord[1] >>= 1;
1113	–	coord[2] = (coord[2] >> 1) + MAXBCOORD2;
1114	–	break;
1115	–
1116	–	case 3:
1117	–	coord[0] = MAXBCOORD2 - (coord[0] >> 1);
1118	–	coord[1] = MAXBCOORD2 - (coord[1] >> 1);
1119	–	coord[2] = MAXBCOORD2 - (coord[2] >> 1);
1120	–	break;
1121	–	#ifdef DEBUG
1122	–	default:
1123	–	eputs("baryi_parent():Invalid child\n");
1124	–	break;
1125	–	#endif
1126	–	}
1127	–	}
1128	–
1129	–	baryi_from_child(coord,child,next)
1130	–	BCOORD coord[3];
1131	–	int child,next;
1132	–	{
1133	–	#ifdef DEBUG
1134	–	if(child <0 || child > 3)
1135	–	{
1136	–	eputs("baryi_from_child():Invalid child\n");
1137	–	return;
1138	–	}
1139	–	if(next <0 || next > 3)
1140	–	{
1141	–	eputs("baryi_from_child():Invalid next\n");
1142	–	return;
1143	–	}
1144	–	#endif
1145	–	if(next == child)
1146	–	return;
1147	–
1148	–	switch(child){
1149	–	case 0:
1150	–	coord[0] = 0;
1151	–	coord[1] = MAXBCOORD - coord[1];
1152	–	coord[2] = MAXBCOORD - coord[2];
1153	–	break;
1154	–	case 1:
1155	–	coord[0] = MAXBCOORD - coord[0];
1156	–	coord[1] = 0;
1157	–	coord[2] = MAXBCOORD - coord[2];
1158	–	break;
1159	–	case 2:
1160	–	coord[0] = MAXBCOORD - coord[0];
1161	–	coord[1] = MAXBCOORD - coord[1];
1162	–	coord[2] = 0;
1163	–	break;
1164	–	case 3:
1165	–	switch(next){
1166	–	case 0:
1167	–	coord[0] = 0;
1168	–	coord[1] = MAXBCOORD - coord[1];
1169	–	coord[2] = MAXBCOORD - coord[2];
1170	–	break;
1171	–	case 1:
1172	–	coord[0] = MAXBCOORD - coord[0];
1173	–	coord[1] = 0;
1174	–	coord[2] = MAXBCOORD - coord[2];
1175	–	break;
1176	–	case 2:
1177	–	coord[0] = MAXBCOORD - coord[0];
1178	–	coord[1] = MAXBCOORD - coord[1];
1179	–	coord[2] = 0;
1180	–	break;
1181	–	}
1182	–	break;
1183	–	}
1184	–	}
1185	–
1186	–	int
1187	–	baryi_child(coord)
1188	–	BCOORD coord[3];
1189	–	{
1190	–	int i;
1191	–
1192	–	if(coord[0] > MAXBCOORD2)
1193	–	{
1194	–	/* update bary for child */
1195	–	coord[0] = (coord[0]<< 1) - MAXBCOORD;
1196	–	coord[1] <<= 1;
1197	–	coord[2] <<= 1;
1198	–	return(0);
1199	–	}
1200	–	else
1201	–	if(coord[1] > MAXBCOORD2)
1202	–	{
1203	–	coord[0] <<= 1;
1204	–	coord[1] = (coord[1] << 1) - MAXBCOORD;
1205	–	coord[2] <<= 1;
1206	–	return(1);
1207	–	}
1208	–	else
1209	–	if(coord[2] > MAXBCOORD2)
1210	–	{
1211	–	coord[0] <<= 1;
1212	–	coord[1] <<= 1;
1213	–	coord[2] = (coord[2] << 1) - MAXBCOORD;
1214	–	return(2);
1215	–	}
1216	–	else
1217	–	{
1218	–	coord[0] = MAXBCOORD - (coord[0] << 1);
1219	–	coord[1] = MAXBCOORD - (coord[1] << 1);
1220	–	coord[2] = MAXBCOORD - (coord[2] << 1);
1221	–	return(3);
1222	–	}
1223	–	}
1224	–
1225	–	int
1226	–	baryi_nth_child(coord,i)
1227	–	BCOORD coord[3];
1228	–	int i;
1229	–	{
1230	–
1231	–	switch(i){
1232	–	case 0:
1233	–	/* update bary for child */
1234	–	coord[0] = (coord[0]<< 1) - MAXBCOORD;
1235	–	coord[1] <<= 1;
1236	–	coord[2] <<= 1;
1237	–	break;
1238	–	case 1:
1239	–	coord[0] <<= 1;
1240	–	coord[1] = (coord[1] << 1) - MAXBCOORD;
1241	–	coord[2] <<= 1;
1242	–	break;
1243	–	case 2:
1244	–	coord[0] <<= 1;
1245	–	coord[1] <<= 1;
1246	–	coord[2] = (coord[2] << 1) - MAXBCOORD;
1247	–	break;
1248	–	case 3:
1249	–	coord[0] = MAXBCOORD - (coord[0] << 1);
1250	–	coord[1] = MAXBCOORD - (coord[1] << 1);
1251	–	coord[2] = MAXBCOORD - (coord[2] << 1);
1252	–	break;
1253	–	}
1254	–	}
1255	–
1256	–
1257	–	baryi_children(coord,i,in_tri,rcoord,rin_tri)
1258	–	BCOORD coord[3][3];
1259	–	int i;
1260	–	int in_tri[3];
1261	–	BCOORD rcoord[3][3];
1262	–	int rin_tri[3];
1263	–	{
1264	–	int j;
1265	–
1266	–	for(j=0; j< 3; j++)
1267	–	{
1268	–	if(!in_tri[j])
1269	–	{
1270	–	rin_tri[j]=0;
1271	–	continue;
1272	–	}
1273	–
1274	–	if(i != 3)
1275	–	{
1276	–	if(coord[j][i] < MAXBCOORD2)
1277	–	{
1278	–	rin_tri[j] = 0;
1279	–	continue;
1280	–	}
1281	–	}
1282	–	else
1283	–	if( !(coord[j][0] <= MAXBCOORD2 && coord[j][1] <= MAXBCOORD2 &&
1284	–	coord[j][2] <= MAXBCOORD2))
1285	–	{
1286	–	rin_tri[j] = 0;
1287	–	continue;
1288	–	}
1289	–
1290	–	rin_tri[j]=1;
1291	–	VCOPY(rcoord[j],coord[j]);
1292	–	baryi_nth_child(rcoord[j],i);
1293	–	}
1294	–
1295	–	}
1296	–
1297	–	convert_dtol(b,bi)
1298	–	double b[3];
1299	–	BCOORD bi[3];
1300	–	{
1301	–	int i;
1302	–
1303	–	for(i=0; i < 2;i++)
1304	–	{
1305	–
1306	–	if(b[i] <= 0.0)
1307	–	{
1308	–	#ifdef EXTRA_DEBUG
1309	–	if(b[i] < 0.0)
1310	–	printf("under %f\n",b[i]);
1311	–	#endif
1312	–	bi[i]= 0;
1313	–	}
1314	–	else
1315	–	if(b[i] >= 1.0)
1316	–	{
1317	–	#ifdef EXTRA_DEBUG
1318	–	if(b[i] > 1.0)
1319	–	printf("over %f\n",b[i]);
1320	–	#endif
1321	–	bi[i]= MAXBCOORD;
1322	–	}
1323	–	else
1324	–	bi[i] = (BCOORD)(b[i]*MAXBCOORD);
1325	–	}
1326	–	bi[2] = bi[0] +  bi[1];
1327	–	if(bi[2] > MAXBCOORD)
1328	–	{
1329	–	#ifdef EXTRA_DEBUG
1330	–	printf("sum over %f\n",b[0]+b[1]);
1331	–	#endif
1332	–	bi[2] = 0;
1333	–	bi[1] = MAXBCOORD - bi[0];
1334	–	}
1335	–	else
1336	–	bi[2] = MAXBCOORD - bi[2];
1337	–
1338	–	}
1339	–
1340	–	/* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1341	–	dir unbounded to [-MAXLONG,MAXLONG]
1342	–	*/
1343	–	bary_dtol(b,db,bi,dbi,t,w)
1344	–	double b[3],db[3][3];
1345	–	BCOORD bi[3];
1346	–	BDIR dbi[3][3];
1347	–	TINT t[3];
1348	–	int w;
1349	–	{
1350	–	int i,id,j,k;
1351	–	double d;
1352	–
1353	–	convert_dtol(b,bi);
1354	–
1355	–	for(j=w; j< 3; j++)
1356	–	{
1357	–	if(t[j] == HUGET)
1358	–	{
1359	–	max_index(db[j],&d);
1360	–	for(i=0; i< 3; i++)
1361	–	dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1362	–	break;
1363	–	}
1364	–	else
1365	–	{
1366	–	for(i=0; i< 3; i++)
1367	–	dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1368	–	}
1369	–	}
1370	–	}
1371	–
1372	–
1373	–	/* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1374	–	dir unbounded to [-MAXLONG,MAXLONG]
1375	–	*/
1376	–	bary_dtol_new(b,db,bi,boi,dbi,t)
1377	–	double b[4][3],db[3][3];
1378	–	BCOORD bi[3],boi[3][3];
1379	–	BDIR dbi[3][3];
1380	–	int t[3];
1381	–	{
1382	–	int i,id,j,k;
1383	–	double d;
1384	–
1385	–	convert_dtol(b[3],bi);
1386	–
1387	–	for(j=0; j<3;j++)
1388	–	{
1389	–	if(t[j] != 1)
1390	–	continue;
1391	–	convert_dtol(b[j],boi[j]);
1392	–	}
1393	–	for(j=0; j< 3; j++)
1394	–	{
1395	–	k = (j+1)%3;
1396	–	if(t[k]==0)
1397	–	continue;
1398	–	if(t[k] == -1)
1399	–	{
1400	–	max_index(db[j],&d);
1401	–	for(i=0; i< 3; i++)
1402	–	dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1403	–	t[k] = 0;
1404	–	}
1405	–	else
1406	–	if(t[j] != 1)
1407	–	for(i=0; i< 3; i++)
1408	–	dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1409	–	else
1410	–	for(i=0; i< 3; i++)
1411	–	dbi[j][i] = boi[k][i] - boi[j][i];
1412	–
1413	–	}
1414	–	}
1415	–
1416	–
1417	–	bary_dtolb(b,bi,in_tri)
1418	–	double b[3][3];
1419	–	BCOORD bi[3][3];
1420	–	int in_tri[3];
1421	–	{
1422	–	int i,j;
1423	–
1424	–	for(j=0; j<3;j++)
1425	–	{
1426	–	if(!in_tri[j])
1427	–	continue;
1428	–	for(i=0; i < 2;i++)
1429	–	{
1430	–	if(b[j][i] <= 0.0)
1431	–	{
1432	–	bi[j][i]= 0;
1433	–	}
1434	–	else
1435	–	if(b[j][i] >= 1.0)
1436	–	{
1437	–	bi[j][i]= MAXBCOORD;
1438	–	}
1439	–	else
1440	–	bi[j][i] = (BCOORD)(b[j][i]*MAXBCOORD);
1441	–	}
1442	–	bi[j][2] = MAXBCOORD - bi[j][0] - bi[j][1];
1443	–	if(bi[j][2] < 0)
1444	–	{
1445	–	bi[j][2] = 0;
1446	–	bi[j][1] = MAXBCOORD - bi[j][0];
1447	–	}
1448	–	}
1449	–	}
356
357
358

Diff Legend

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