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.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,pd,r)
135	<	FVECT v,plane_n;
136	<	double plane_d;
137	<	double *pd;
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,pd,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
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,peq,pd,r)
177	>	FVECT orig,dir;
178	>	FPEQ peq;
179	>	double *pd;
180	>	FVECT r;
181	>	{
182	>	double t;
183	>	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	<	r[0] = v[0]*t;
197	<	r[1] = v[1]*t;
198	<	r[2] = v[2]*t;
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_plane(orig,dir,plane_n,plane_d,pd,r)
209	>	intersect_ray_oplane(orig,dir,n,pd,r)
210		FVECT orig,dir;
211	<	FVECT plane_n;
172	<	double plane_d;
211	>	FVECT n;
212		double *pd;
213		FVECT r;
214		{
#	Line 179 | Line 218 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
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 */
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(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
227	<	if(ZERO(t) || t >0)
228	<	hit = 1;
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 202 | Line 277 | point_in_cone(p,p0,p1,p2)
277		FVECT p;
278		FVECT p0,p1,p2;
279		{
205	–	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 212 | 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 251 | 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;
259	<	char sides[3];
333	>	int *nset;
334	>	int sides[3];
335
336		{
337	<	float d;
263	<
337	>	double d;
338		/* Find the normal to the triangle ORIGIN:v0:v1 */
339		if(!NTH_BIT(*nset,0))
340		{
#	Line 270 | 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	<	{
278	<	sides[0] =  GT_OUT;
279	<	sides[1] = sides[2] = GT_INVALID;
280	<	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 290 | 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		{
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	–	{
366		sides[1] = GT_OUT;
367		sides[2] = GT_INVALID;
368		return(FALSE);
#	Line 312 | Line 372 | char sides[3];
372		/* Test next edge */
373		if(!NTH_BIT(*nset,2))
374		{
315	–
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	<
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	<	}
382	>	sides[2] = GT_OUT;
383	>	return(FALSE);
384		}
344	–	else if(d > 0)
345	–	{
346	–	sides[2] = GT_OUT;
347	–	return(FALSE);
348	–	}
349	–	/* If on edge */
385		else
386		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	–	}
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;
375	–	char *which;
397		{
398	<	float d;
398	>	double d;
399		FVECT n;
379	–	char sides[3];
400
381	–	/* 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	–	}
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
403	<	if(d > 0)
404	<	return(FALSE);
405	<	else
406	<	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))
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)
411	>	if(d > 0.0)
412		return(FALSE);
423	–	else
424	–	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))
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)
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;
468	<	char which;
433	>	int below_plane[3],test;
434
435		SUM_3VEC3(avg,t0,t1,t2);
471	–	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 476 | 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,
485	<	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);
491	–	/* Remember if b[i] fell on one of the 3 defining planes */
492	–	if(test)
493	–	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,
503	<	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);
509	–	/* Remember if b[i] fell on one of the 3 defining planes */
510	–	if(test)
511	–	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,
521	<	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);
528	–	/* Remember if b[i] fell on one of the 3 defining planes */
529	–	if(test)
530	–	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);
536	–	/* Accept if all points lie on a triangle vertex/edge edge- accept*/
537	–	if(on_edge == 3)
538	–	return(TRUE);
488		/* Now check vertices in a against triangle b */
489		return(FALSE);
490		}
#	Line 543 | 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 560 | 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 573 | 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 586 | 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 602 | 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 616 | 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 630 | 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 646 | 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 659 | 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 672 | 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 682 | 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 695 | 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,
699	<	&(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,
703	<	&(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 748 | 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;
755	–	char *wptr;
702		{
703	<	FVECT p0,p1,p2,p,n;
704	<	char type,which;
705	<	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);
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)
773	<	*wptr = which;
774	<
775	<	return(type);
717	>	return(FALSE);
718		}
719
720
#	Line 831 | 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 879 | Line 859 | FVECT a0,a1,b0,b1;
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
867	+	*/
868	+	bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
869	+	double x1,y1,x2,y2,x3,y3;
870	+	double px,py;
871	+	double coord[3];
872	+	{
873	+	double a;
874	+
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] = ((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	+	bary_child(coord)
920	+	double coord[3];
921	+	{
922	+	int i;
923	+
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	+	}
1450	+
1451	+
1452	+
1453	+
1454	+
1455	+
1456	+
1457	+
1458	+
1459	+
1460

Diff Legend

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