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.7 by gwlarson, Tue Oct 6 18:16:53 1998 UTC

#	Line 33 | Line 33 | int
33		convex_angle(v0,v1,v2)
34		FVECT v0,v1,v2;
35		{
36	<	FVECT cp01,cp12,cp;
37	<
38	<	/* test sign of (v0Xv1)X(v1Xv2). v1 */
36	>	FVECT cp,cp01,cp12,v10,v02;
37	>	double dp;
38	>	/*
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	<	if(DOT(cp,v1) < 0)
47	>
48	>	dp = DOT(cp,v1);
49	>	if(ZERO(dp) || dp < 0.0)
50		return(FALSE);
51		return(TRUE);
52		}
#	Line 53 | Line 60 | FVECT v0,v1,v2;
60		double
61		tri_normal(v0,v1,v2,n,norm)
62		FVECT v0,v1,v2,n;
63	<	char norm;
63	>	int norm;
64		{
65		double mag;
66
#	Line 64 | Line 71 | char norm;
71		n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
72		(v1[2] - v2[2]) * (v1[0] + v2[0]) +
73		(v2[2] - v0[2]) * (v2[0] + v0[0]);
67	–
74
75		n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
76		(v1[1] + v2[1]) * (v1[0] - v2[0]) +
#	Line 80 | Line 86 | char norm;
86		}
87
88
89	<	tri_plane_equation(v0,v1,v2,n,nd,norm)
90	<	FVECT v0,v1,v2,n;
91	<	double *nd;
92	<	char norm;
89	>	tri_plane_equation(v0,v1,v2,peqptr,norm)
90	>	FVECT v0,v1,v2;
91	>	FPEQ *peqptr;
92	>	int norm;
93		{
94	<	tri_normal(v0,v1,v2,n,norm);
94	>	tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
95
96	<	*nd = -(DOT(n,v0));
96	>	FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
97		}
98
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	–
99		/* From quad_edge-code */
100		int
101		point_in_circle_thru_origin(p,p0,p1)
#	Line 135 | Line 125 | FVECT ps,p,c;
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,plane_n,plane_d,tptr,r)
135	<	FVECT v,plane_n;
136	<	double plane_d;
134	>	intersect_vector_plane(v,peq,tptr,r)
135	>	FVECT v;
136	>	FPEQ peq;
137		double *tptr;
138		FVECT r;
139		{
140	<	double t;
140	>	double t,d;
141		int hit;
142		/*
143		Plane is Ax + By + Cz +D = 0:
#	Line 152 | Line 147 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
147		/* line is  l = p1 + (p2-p1)t, p1=origin */
148
149		/* Solve for t: */
150	<	t =  plane_d/-(DOT(plane_n,v));
151	<	if(t >0 || ZERO(t))
152	<	hit = 1;
153	<	else
154	<	hit = 0;
155	<	r[0] = v[0]*t;
156	<	r[1] = v[1]*t;
157	<	r[2] = v[2]*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
176	<	intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
176	>	intersect_ray_plane(orig,dir,peq,pd,r)
177		FVECT orig,dir;
178	<	FVECT plane_n;
172	<	double plane_d;
178	>	FPEQ peq;
179		double *pd;
180		FVECT r;
181		{
#	Line 184 | Line 190 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
190		line is  l = p1 + (p2-p1)t
191		*/
192		/* Solve for t: */
193	<	t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
194	<	if(ZERO(t) || t >0)
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;
198	+
199	+	if(r)
200	+	VSUM(r,orig,dir,t);
201	+
202	+	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)
228	+	hit = 0;
229		else
230	+	hit = 1;
231	+
232	+	if(r)
233	+	VSUM(r,orig,dir,t);
234	+
235	+	if(pd)
236	+	*pd = t;
237	+	return(hit);
238	+	}
239	+
240	+
241	+	int
242	+	intersect_edge_plane(e0,e1,peq,pd,r)
243	+	FVECT e0,e1;
244	+	FPEQ peq;
245	+	double *pd;
246	+	FVECT r;
247	+	{
248	+	double t;
249	+	int hit;
250	+	FVECT d;
251	+	/*
252	+	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,orig,dir,t);
267	>	VSUM(r,e0,d,t);
268
269		if(pd)
270		*pd = t;
#	Line 203 | Line 277 | point_in_cone(p,p0,p1,p2)
277		FVECT p;
278		FVECT p0,p1,p2;
279		{
206	–	FVECT n;
280		FVECT np,x_axis,y_axis;
281	<	double d1,d2,d;
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
#	Line 213 | Line 287 | FVECT p0,p1,p2;
287		*/
288
289		/* find the equation of the plane defined by p1-p3 */
290	<	tri_plane_equation(p0,p1,p2,n,&d,FALSE);
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	<	if(!intersect_vector_plane(p,n,d,NULL,np))
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,n,x_axis);
305	>	VCROSS(y_axis,FP_N(peq),x_axis);
306		normalize(y_axis);
307
308		VSUB(p1,p1,p0);
#	Line 252 | Line 327 | FVECT p0,p1,p2;
327		}
328
329		int
330	<	test_point_against_spherical_tri(v0,v1,v2,p,n,nset,which,sides)
330	>	point_set_in_stri(v0,v1,v2,p,n,nset,sides)
331		FVECT v0,v1,v2,p;
332		FVECT n[3];
333	<	char *nset;
334	<	char *which;
260	<	char sides[3];
333	>	int *nset;
334	>	int sides[3];
335
336		{
337	<	float d;
264	<
337	>	double d;
338		/* Find the normal to the triangle ORIGIN:v0:v1 */
339		if(!NTH_BIT(*nset,0))
340		{
#	Line 271 | Line 344 | char sides[3];
344		/* Test the point for sidedness */
345		d  = DOT(n[0],p);
346
347	<	if(ZERO(d))
348	<	sides[0] = GT_EDGE;
349	<	else
350	<	if(d > 0)
351	<	{
279	<	sides[0] =  GT_OUT;
280	<	sides[1] = sides[2] = GT_INVALID;
281	<	return(FALSE);
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;
#	Line 291 | Line 361 | char sides[3];
361		}
362		/* Test the point for sidedness */
363		d  = DOT(n[1],p);
364	<	if(ZERO(d))
364	>	if(d > 0.0)
365		{
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	–	{
366		sides[1] = GT_OUT;
367		sides[2] = GT_INVALID;
368		return(FALSE);
#	Line 313 | Line 372 | char sides[3];
372		/* Test next edge */
373		if(!NTH_BIT(*nset,2))
374		{
316	–
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(ZERO(d))
380	>	if(d > 0.0)
381		{
382	<	sides[2] = GT_EDGE;
383	<
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	<	}
382	>	sides[2] = GT_OUT;
383	>	return(FALSE);
384		}
345	–	else if(d > 0)
346	–	{
347	–	sides[2] = GT_OUT;
348	–	return(FALSE);
349	–	}
350	–	/* If on edge */
385		else
386		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	–	}
387		/* Must be interior to the pyramid */
388		return(GT_INTERIOR);
389		}
390
391
392
393	<
393	>
394		int
395	<	test_single_point_against_spherical_tri(v0,v1,v2,p,which)
395	>	point_in_stri(v0,v1,v2,p)
396		FVECT v0,v1,v2,p;
376	–	char *which;
397		{
398	<	float d;
398	>	double d;
399		FVECT n;
380	–	char sides[3];
400
382	–	/* 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	–	}
401		VCROSS(n,v1,v0);
402		/* Test the point for sidedness */
403		d  = DOT(n,p);
404	<	if(ZERO(d))
405	<	sides[0] = GT_EDGE;
406	<	else
404	<	if(d > 0)
405	<	return(FALSE);
406	<	else
407	<	sides[0] = GT_INTERIOR;
404	>	if(d > 0.0)
405	>	return(FALSE);
406	>
407		/* Test next edge */
408		VCROSS(n,v2,v1);
409		/* Test the point for sidedness */
410		d  = DOT(n,p);
411	<	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)
411	>	if(d > 0.0)
412		return(FALSE);
424	–	else
425	–	sides[1] = GT_INTERIOR;
413
414		/* Test next edge */
415		VCROSS(n,v0,v2);
416		/* Test the point for sidedness */
417		d  = DOT(n,p);
418	<	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)
418	>	if(d > 0.0)
419		return(FALSE);
420		/* Must be interior to the pyramid */
421	<	return(GT_FACE);
421	>	return(GT_INTERIOR);
422		}
423
424		int
425	<	test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
425	>	vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
426		FVECT t0,t1,t2,p0,p1,p2;
427	<	char *nset;
427	>	int *nset;
428		FVECT n[3];
429		FVECT avg;
430	<	char pt_sides[3][3];
430	>	int pt_sides[3][3];
431
432		{
433	<	char below_plane[3],on_edge,test;
469	<	char which;
433	>	int below_plane[3],test;
434
435		SUM_3VEC3(avg,t0,t1,t2);
472	–	on_edge = 0;
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
#	Line 477 | Line 440 | char pt_sides[3][3];
440		*/
441
442		/* test point 0 */
443	<	if(DOT(avg,p0) < 0)
443	>	if(DOT(avg,p0) < 0.0)
444		below_plane[0] = 1;
445		else
446	<	below_plane[0]=0;
446	>	below_plane[0] = 0;
447		/* Test if b[i] lies in or on triangle a */
448	<	test = test_point_against_spherical_tri(t0,t1,t2,p0,
486	<	n,nset,&which,pt_sides[0]);
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);
492	–	/* Remember if b[i] fell on one of the 3 defining planes */
493	–	if(test)
494	–	on_edge++;
454		}
455		/* Now test point 1*/
456
457	<	if(DOT(avg,p1) < 0)
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 = test_point_against_spherical_tri(t0,t1,t2,p1,
504	<	n,nset,&which,pt_sides[1]);
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);
510	–	/* Remember if b[i] fell on one of the 3 defining planes */
511	–	if(test)
512	–	on_edge++;
468		}
469
470		/* Now test point 2 */
471	<	if(DOT(avg,p2) < 0)
471	>	if(DOT(avg,p2) < 0.0)
472		below_plane[2] = 1;
473		else
474	<	below_plane[2]=0;
474	>	below_plane[2] = 0;
475		/* Test if b[i] lies in or on triangle a */
476	<	test = test_point_against_spherical_tri(t0,t1,t2,p2,
522	<	n,nset,&which,pt_sides[2]);
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);
529	–	/* Remember if b[i] fell on one of the 3 defining planes */
530	–	if(test)
531	–	on_edge++;
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);
537	–	/* Accept if all points lie on a triangle vertex/edge edge- accept*/
538	–	if(on_edge == 3)
539	–	return(TRUE);
488		/* Now check vertices in a against triangle b */
489		return(FALSE);
490		}
#	Line 544 | Line 492 | char pt_sides[3][3];
492
493		set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
494		FVECT t0,t1,t2,p0,p1,p2;
495	<	char test[3];
496	<	char sides[3][3];
497	<	char nset;
495	>	int test[3];
496	>	int sides[3][3];
497	>	int nset;
498		FVECT n[3];
499		{
500	<	char t;
500	>	int t;
501		double d;
502
503
#	Line 561 | Line 509 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
509		VCROSS(n[0],t1,t0);
510		/* Test the point for sidedness */
511		d  = DOT(n[0],p0);
512	<	if(d >= 0)
512	>	if(d >= 0.0)
513		SET_NTH_BIT(test[0],0);
514		}
515		else
#	Line 574 | Line 522 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
522		VCROSS(n[1],t2,t1);
523		/* Test the point for sidedness */
524		d  = DOT(n[1],p0);
525	<	if(d >= 0)
525	>	if(d >= 0.0)
526		SET_NTH_BIT(test[0],1);
527		}
528		else
#	Line 587 | Line 535 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
535		VCROSS(n[2],t0,t2);
536		/* Test the point for sidedness */
537		d  = DOT(n[2],p0);
538	<	if(d >= 0)
538	>	if(d >= 0.0)
539		SET_NTH_BIT(test[0],2);
540		}
541		else
#	Line 603 | Line 551 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
551		VCROSS(n[0],t1,t0);
552		/* Test the point for sidedness */
553		d  = DOT(n[0],p1);
554	<	if(d >= 0)
554	>	if(d >= 0.0)
555		SET_NTH_BIT(test[1],0);
556		}
557		else
#	Line 617 | Line 565 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
565		VCROSS(n[1],t2,t1);
566		/* Test the point for sidedness */
567		d  = DOT(n[1],p1);
568	<	if(d >= 0)
568	>	if(d >= 0.0)
569		SET_NTH_BIT(test[1],1);
570		}
571		else
#	Line 631 | Line 579 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
579		VCROSS(n[2],t0,t2);
580		/* Test the point for sidedness */
581		d  = DOT(n[2],p1);
582	<	if(d >= 0)
582	>	if(d >= 0.0)
583		SET_NTH_BIT(test[1],2);
584		}
585		else
#	Line 647 | Line 595 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
595		VCROSS(n[0],t1,t0);
596		/* Test the point for sidedness */
597		d  = DOT(n[0],p2);
598	<	if(d >= 0)
598	>	if(d >= 0.0)
599		SET_NTH_BIT(test[2],0);
600		}
601		else
#	Line 660 | Line 608 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
608		VCROSS(n[1],t2,t1);
609		/* Test the point for sidedness */
610		d  = DOT(n[1],p2);
611	<	if(d >= 0)
611	>	if(d >= 0.0)
612		SET_NTH_BIT(test[2],1);
613		}
614		else
#	Line 673 | Line 621 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
621		VCROSS(n[2],t0,t2);
622		/* Test the point for sidedness */
623		d  = DOT(n[2],p2);
624	<	if(d >= 0)
624	>	if(d >= 0.0)
625		SET_NTH_BIT(test[2],2);
626		}
627		else
#	Line 683 | Line 631 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
631
632
633		int
634	<	spherical_tri_intersect(a1,a2,a3,b1,b2,b3)
634	>	stri_intersect(a1,a2,a3,b1,b2,b3)
635		FVECT a1,a2,a3,b1,b2,b3;
636		{
637	<	char which,test,n_set[2];
638	<	char sides[2][3][3],i,j,inext,jnext;
639	<	char tests[2][3];
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
#	Line 696 | Line 644 | FVECT a1,a2,a3,b1,b2,b3;
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(test_vertices_for_tri_inclusion(a1,a2,a3,b1,b2,b3,
700	<	&(n_set[0]),n[0],avg[0],sides[1]))
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(test_vertices_for_tri_inclusion(b1,b2,b3,a1,a2,a3,
704	<	&(n_set[1]),n[1],avg[1],sides[0]))
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
#	Line 749 | Line 695 | FVECT a1,a2,a3,b1,b2,b3;
695		}
696
697		int
698	<	ray_intersect_tri(orig,dir,v0,v1,v2,pt,wptr)
698	>	ray_intersect_tri(orig,dir,v0,v1,v2,pt)
699		FVECT orig,dir;
700		FVECT v0,v1,v2;
701		FVECT pt;
756	–	char *wptr;
702		{
703	<	FVECT p0,p1,p2,p,n;
704	<	char type,which;
705	<	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);
703	>	FVECT p0,p1,p2,p;
704	>	FPEQ peq;
705	>	int type;
706
707	<	if(type)
707	>	VSUB(p0,v0,orig);
708	>	VSUB(p1,v1,orig);
709	>	VSUB(p2,v2,orig);
710	>
711	>	if(point_in_stri(p0,p1,p2,dir))
712		{
713		/* Intersect the ray with the triangle plane */
714	<	tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
715	<	intersect_ray_plane(orig,dir,n,pd,NULL,pt);
714	>	tri_plane_equation(v0,v1,v2,&peq,FALSE);
715	>	return(intersect_ray_plane(orig,dir,peq,NULL,pt));
716		}
717	<	if(wptr)
774	<	*wptr = which;
775	<
776	<	return(type);
717	>	return(FALSE);
718		}
719
720
#	Line 832 | Line 773 | FVECT fnear[4],ffar[4];
773		ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
774		}
775
776	+	int
777	+	max_index(v,r)
778	+	FVECT v;
779	+	double *r;
780	+	{
781	+	double p[3];
782	+	int i;
783
784	+	p[0] = fabs(v[0]);
785	+	p[1] = fabs(v[1]);
786	+	p[2] = fabs(v[2]);
787	+	i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);
788	+	if(r)
789	+	*r = p[i];
790	+	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	<	spherical_polygon_edge_intersect(a0,a1,b0,b1)
818	>	sedge_intersect(a0,a1,b0,b1)
819		FVECT a0,a1,b0,b1;
820		{
821		FVECT na,nb,avga,avgb,p;
#	Line 896 | Line 875 | double coord[3];
875		a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
876		coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
877		coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
878	<	coord[2]  = 1.0 - coord[0] - coord[1];
878	>	coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
879
880		}
881
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	+	}
916	+
917	+
918		int
919	<	bary2d_child(coord)
919	>	bary_child(coord)
920		double coord[3];
921		{
922		int i;
923
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 */
924		if(coord[0] > 0.5)
925		{
926		/* update bary for child */
#	Line 968 | Line 954 | double coord[3];
954		}
955		}
956
957	<	int
958	<	max_index(v)
959	<	FVECT v;
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		{
975	–	double a,b,c;
976	–	int i;
964
965	<	a = fabs(v[0]);
966	<	b = fabs(v[1]);
967	<	c = fabs(v[2]);
968	<	i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);
969	<	return(i);
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	<	* int
1035	<	* smRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
1036	<	*
1037	<	*   Intersect the ray with triangle v0v1v2, return intersection point in r
1038	<	*
1039	<	*    Assumes orig,v0,v1,v2 are in spherical coordinates, and orig is
1040	<	*    unit
1041	<	*/
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	<	traceRay(orig,dir,v0,v1,v2,r)
1188	<	FVECT orig,dir;
999	<	FVECT v0,v1,v2;
1000	<	FVECT r;
1187	>	baryi_child(coord)
1188	>	BCOORD coord[3];
1189		{
1190	<	FVECT n,p[3],d;
1003	<	double pt[3],r_eps;
1004	<	char i;
1005	<	int which;
1190	>	int i;
1191
1192	<	/* Find the plane equation for the triangle defined by the edge v0v1 and
1193	<	the view center*/
1194	<	VCROSS(n,v0,v1);
1195	<	/* Intersect the ray with this plane */
1196	<	i = intersect_ray_plane(orig,dir,n,0.0,&(pt[0]),p[0]);
1197	<	if(i)
1198	<	which = 0;
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	<	which = -1;
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	<	VCROSS(n,v1,v2);
1226	<	i = intersect_ray_plane(orig,dir,n,0.0,&(pt[1]),p[1]);
1227	<	if(i)
1228	<	if((which==-1) || pt[1] < pt[0])
1229	<	which = 1;
1225	>	int
1226	>	baryi_nth_child(coord,i)
1227	>	BCOORD coord[3];
1228	>	int i;
1229	>	{
1230
1231	<	VCROSS(n,v2,v0);
1232	<	i = intersect_ray_plane(orig,dir,n,0.0,&(pt[2]),p[2]);
1233	<	if(i)
1234	<	if((which==-1) || pt[2] < pt[which])
1235	<	which = 2;
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	<	if(which != -1)
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	<	/* Return point that is K*eps along projection of the ray on
1269	<	the sphere to push intersection point p[which] into next cell
1270	<	*/
1271	<	normalize(p[which]);
1272	<	/* Calculate the ray perpendicular to the intersection point: approx
1273	<	the direction moved along the sphere a distance of K*epsilon*/
1274	<	r_eps = -DOT(p[which],dir);
1275	<	VSUM(n,dir,p[which],r_eps);
1276	<	/* Calculate the length along ray p[which]-dir needed to move to
1277	<	cause a move along the sphere of k*epsilon
1278	<	*/
1279	<	r_eps = DOT(n,dir);
1280	<	VSUM(r,p[which],dir,(20*FTINY)/r_eps);
1281	<	normalize(r);
1282	<	return(TRUE);
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	<	/* Unable to find intersection: move along ray and try again */
1358	<	r_eps = -DOT(orig,dir);
1359	<	VSUM(n,dir,orig,r_eps);
1360	<	r_eps = DOT(n,dir);
1361	<	VSUM(r,orig,dir,(20*FTINY)/r_eps);
1362	<	normalize(r);
1363	<	#ifdef DEBUG
1364	<	eputs("traceRay:Ray does not intersect triangle");
1365	<	#endif
1366	<	return(FALSE);
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	+	}
1450	+
1451	+
1452
1453
1454

Diff Legend

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