src/hd/sm_geom.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 *  sm_geom.c
 *    some geometric utility routines
 */

#include "standard.h"
#include "sm_geom.h"

/*
 * int
 * pt_in_cone(p,a,b,c)
 *             : test if point p lies in cone defined by a,b,c and origin
 * double p[3];              : point to test
 * double a[3],b[3],c[3];    : points forming triangle
 * 
 *  Assumes apex at origin, a,b,c are unit vectors defining the
 *  triangle which the cone circumscribes. Assumes p is also normalized
 *  Test is implemented as:
 *             r =  (b-a)X(c-a) 
 *             in = (p.r) > (a.r) 
 *  The center of the cone is r, and corresponds to the triangle normal.
 *  p.r is the proportional to the cosine of the angle between p and the
 *  the cone center ray, and a.r to the radius of the cone. If the cosine
 *  of the angle for p is greater than that for a, the angle between p
 *  and r is smaller, and p lies in the cone.
 */
int
pt_in_cone(p,a,b,c)
double p[3],a[3],b[3],c[3];
{
  double r[3];
  double pr,ar;
  double ab[3],ac[3];
 
#ifdef DEBUG
#if DEBUG > 1
{
  double l;
  VSUB(ab,b,a);
  normalize(ab);
  VSUB(ac,c,a);
  normalize(ac);
  VCROSS(r,ab,ac);
  l = normalize(r);
  /* l = sin@ between ab,ac - if 0 vectors are colinear */
  if( l <= COLINEAR_EPS)
  {
    eputs("pt in cone: null triangle:returning FALSE\n");
    return(FALSE);
  }
}
#endif
#endif

  VSUB(ab,b,a);
  VSUB(ac,c,a);
  VCROSS(r,ab,ac);

  pr = DOT(p,r);        
  ar = DOT(a,r);
  /* Need to check for equality for degeneracy of 4 points on circle */
  if( pr > ar *( 1.0 + EQUALITY_EPS))
    return(TRUE); 
  else
    return(FALSE);
}

/*
 * tri_centroid(v0,v1,v2,c)
 *             : Average triangle vertices to give centroid: return in c
 *FVECT v0,v1,v2,c; : triangle vertices(v0,v1,v2) and vector to hold result(c)
 */                   
tri_centroid(v0,v1,v2,c)
FVECT v0,v1,v2,c;
{
  c[0] = (v0[0] + v1[0] + v2[0])/3.0;
  c[1] = (v0[1] + v1[1] + v2[1])/3.0;
  c[2] = (v0[2] + v1[2] + v2[2])/3.0;
}


/*
 * double  
 * tri_normal(v0,v1,v2,n,norm) : Calculates the normal of a face contour using
 *                            Newell's formula.
 * FVECT v0,v1,v2,n;  : Triangle vertices(v0,v1,v2) and vector for result(n)
 * int norm;          : If true result is normalized
 *
 * Triangle normal is calculated using the following:
 *    A =  SUMi (yi - yi+1)(zi + zi+1);
 *    B =  SUMi (zi - zi+1)(xi + xi+1)
 *    C =  SUMi (xi - xi+1)(yi + yi+1)
*/
double
tri_normal(v0,v1,v2,n,norm)
FVECT v0,v1,v2,n;
int norm;
{
  double mag;

  n[0] = (v0[2] + v1[2]) * (v0[1] - v1[1]) + 
         (v1[2] + v2[2]) * (v1[1] - v2[1])  +
         (v2[2] + v0[2]) * (v2[1] - v0[1]);
  n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
           (v1[2] - v2[2]) * (v1[0] + v2[0]) +
            (v2[2] - v0[2]) * (v2[0] + v0[0]);
  n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
         (v1[1] + v2[1]) * (v1[0] - v2[0]) +
         (v2[1] + v0[1]) * (v2[0] - v0[0]);
  if(!norm)
     return(0);
  mag = normalize(n);
  return(mag);
}

/*
 * tri_plane_equation(v0,v1,v2,peqptr,norm)
 *             : Calculates the plane equation (A,B,C,D) for triangle
 *               v0,v1,v2 ( Ax + By + Cz = D)
 * FVECT v0,v1,v2;   : Triangle vertices 
 * FPEQ *peqptr;     : ptr to structure to hold result
 *  int norm;        : if TRUE, return unit normal
 */
tri_plane_equation(v0,v1,v2,peqptr,norm)
   FVECT v0,v1,v2; 
   FPEQ *peqptr;
   int norm;
{
    tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
    FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
}

/*
 * int 
 * intersect_ray_plane(orig,dir,peq,pd,r)
 *             : Intersects ray (orig,dir) with plane (peq). Returns TRUE
 *             if intersection occurs. If r!=NULL, sets with resulting i
 *             intersection point, and pd is set with parametric value of the
 *             intersection.
 * FVECT orig,dir;      : vectors defining the ray
 * FPEQ peq;            : plane equation
 * double *pd;          : holds resulting parametric intersection point
 * FVECT r;             : holds resulting intersection point
 *
 *    Plane is Ax + By + Cz +D = 0:
 *    A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
 *     t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
 *      line is  l = p1 + (p2-p1)t 
 *    Solve for t
 */
int
intersect_ray_plane(orig,dir,peq,pd,r)
   FVECT orig,dir;
   FPEQ peq;
   double *pd;
   FVECT r;
{
  double t,d;
  int hit;

  d = DOT(FP_N(peq),dir);
  if(ZERO(d))
     return(0);
  t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d; 

  if(t < 0)
       hit = 0;
    else
       hit = 1;
  if(r)
     VSUM(r,orig,dir,t);

    if(pd)
       *pd = t;
  return(hit);
}

/*
 * double
 * point_on_sphere(ps,p,c) : normalize p relative to sphere with center c
 * FVECT ps,p,c;       : ps Holds result vector,p is the original point,
 *                       and c is the sphere center
 */
double
point_on_sphere(ps,p,c)
FVECT ps,p,c;
{
  double d;
  
  VSUB(ps,p,c);    
  d = normalize(ps);
  return(d);
}

/*
 * int
 * point_in_stri(v0,v1,v2,p) : Return TRUE if p is in pyramid defined by
 *                            tri v0,v1,v2 and origin
 * FVECT v0,v1,v2,p;     :Triangle vertices(v0,v1,v2) and point in question(p)
 *
 *   Tests orientation of p relative to each edge (v0v1,v1v2,v2v0), if it is
 *   inside of all 3 edges, returns TRUE, else FALSE.
 */
int
point_in_stri(v0,v1,v2,p)
FVECT v0,v1,v2,p;
{
    double d;
    FVECT n;  

    VCROSS(n,v0,v1);
    /* Test the point for sidedness */
    d  = DOT(n,p);
    if(d > 0.0)
      return(FALSE);
    /* Test next edge */
    VCROSS(n,v1,v2);
    /* Test the point for sidedness */
    d  = DOT(n,p);
    if(d > 0.0)
       return(FALSE);
    /* Test next edge */
    VCROSS(n,v2,v0);
    /* Test the point for sidedness */
    d  = DOT(n,p);
    if(d > 0.0)
       return(FALSE);
    /* Must be interior to the pyramid */
    return(TRUE);
}

/*
 * int
 * ray_intersect_tri(orig,dir,v0,v1,v2,pt)
 *    : test if ray orig-dir intersects triangle v0v1v2, result in pt
 * FVECT orig,dir;   : Vectors defining ray origin and direction
 * FVECT v0,v1,v2;   : Triangle vertices
 * FVECT pt;         : Intersection point (if any)
 */
int
ray_intersect_tri(orig,dir,v0,v1,v2,pt)
FVECT orig,dir;
FVECT v0,v1,v2;
FVECT pt;
{
  FVECT p0,p1,p2,p;
  FPEQ peq;
  int type;

  VSUB(p0,v0,orig);
  VSUB(p1,v1,orig);
  VSUB(p2,v2,orig);

  if(point_in_stri(p0,p1,p2,dir))
  {
      /* Intersect the ray with the triangle plane */
      tri_plane_equation(v0,v1,v2,&peq,FALSE);
      return(intersect_ray_plane(orig,dir,peq,NULL,pt));
  }
  return(FALSE);
}

/*
 * calculate_view_frustum(vp,hv,vv,horiz,vert,near,far,fnear,ffar)
 *    : Calculate vertices defining front and rear clip rectangles of
 *      view frustum defined by vp,hv,vv,horiz,vert,near, and far and
 *      return in fnear and ffar.
 * FVECT vp,hv,vv;        : Viewpoint(vp),hv and vv are the horizontal and
 *                          vertical vectors in the view frame-magnitude is
 *                          the dimension of the front frustum face at z =1
 * double horiz,vert,near,far; : View angle horizontal and vertical(horiz,vert)
 *                              and distance to the near,far clipping planes
 * FVECT fnear[4],ffar[4];     : holds results
 * 
 */
calculate_view_frustum(vp,hv,vv,horiz,vert,near,far,fnear,ffar)
FVECT vp,hv,vv;
double horiz,vert,near,far;
FVECT fnear[4],ffar[4];
{
    double height,width;
    FVECT t,nhv,nvv,ndv;
    double w2,h2;
    /* Calculate the x and y dimensions of the near face */
    VCOPY(nhv,hv);
    VCOPY(nvv,vv);
    w2 = normalize(nhv);
    h2 = normalize(nvv);
    /* Use similar triangles to calculate the dimensions at z=near */
    width  = near*0.5*w2;
    height = near*0.5*h2;

    VCROSS(ndv,nvv,nhv);
    /* Calculate the world space points corresponding to the 4 corners
       of the front face of the view frustum
     */
    fnear[0][0] =  width*nhv[0] + height*nvv[0] + near*ndv[0] + vp[0] ;
    fnear[0][1] =  width*nhv[1] + height*nvv[1] + near*ndv[1] + vp[1];
    fnear[0][2] =  width*nhv[2] + height*nvv[2] + near*ndv[2] + vp[2];
    fnear[1][0] = -width*nhv[0] + height*nvv[0] + near*ndv[0] + vp[0];
    fnear[1][1] = -width*nhv[1] + height*nvv[1] + near*ndv[1] + vp[1];
    fnear[1][2] = -width*nhv[2] + height*nvv[2] + near*ndv[2] + vp[2];
    fnear[2][0] = -width*nhv[0] - height*nvv[0] + near*ndv[0] + vp[0];
    fnear[2][1] = -width*nhv[1] - height*nvv[1] + near*ndv[1] + vp[1];
    fnear[2][2] = -width*nhv[2] - height*nvv[2] + near*ndv[2] + vp[2];
    fnear[3][0] =  width*nhv[0] - height*nvv[0] + near*ndv[0] + vp[0];
    fnear[3][1] =  width*nhv[1] - height*nvv[1] + near*ndv[1] + vp[1];
    fnear[3][2] =  width*nhv[2] - height*nvv[2] + near*ndv[2] + vp[2];

    /* Now do the far face */
    width  = far*0.5*w2;
    height = far*0.5*h2;
    ffar[0][0] =  width*nhv[0] + height*nvv[0] + far*ndv[0] + vp[0] ;
    ffar[0][1] =  width*nhv[1] + height*nvv[1] + far*ndv[1] + vp[1] ;
    ffar[0][2] =  width*nhv[2] + height*nvv[2] + far*ndv[2] + vp[2] ;
    ffar[1][0] = -width*nhv[0] + height*nvv[0] + far*ndv[0] + vp[0] ;
    ffar[1][1] = -width*nhv[1] + height*nvv[1] + far*ndv[1] + vp[1] ;
    ffar[1][2] = -width*nhv[2] + height*nvv[2] + far*ndv[2] + vp[2] ;
    ffar[2][0] = -width*nhv[0] - height*nvv[0] + far*ndv[0] + vp[0] ;
    ffar[2][1] = -width*nhv[1] - height*nvv[1] + far*ndv[1] + vp[1] ;
    ffar[2][2] = -width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
    ffar[3][0] =  width*nhv[0] - height*nvv[0] + far*ndv[0] + vp[0] ;
    ffar[3][1] =  width*nhv[1] - height*nvv[1] + far*ndv[1] + vp[1] ;
    ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
}


/*
 * bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
 *  : Find the normalized barycentric coordinates of p relative to
 *    triangle v0,v1,v2. Return result in coord
 * double x1,y1,x2,y2,x3,y3; : defines triangle vertices 1,2,3
 * double px,py;             : coordinates of pt
 * double coord[3];          : result
 */
bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
double x1,y1,x2,y2,x3,y3;
double px,py;
double coord[3];
{
  double a;

  a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
  coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a; 
  coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
  coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
 
}


Revision:	3.13
Committed:	Sat Feb 22 02:07:25 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad3R6P1, rad3R5, rad3R6
Changes since 3.12:	+155 -628 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	Content
1	#ifndef lint
2	static const char RCSid[] = "$Id$";
3	#endif
4	/*
5	* sm_geom.c
6	* some geometric utility routines
7	*/
8
9	#include "standard.h"
10	#include "sm_geom.h"
11
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	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	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	* 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	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
84
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)
99	FVECT v0,v1,v2,n;
100	int norm;
101	{
102	double mag;
103
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]);
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]);
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]);
113	if(!norm)
114	return(0);
115	mag = normalize(n);
116	return(mag);
117	}
118
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;
129	FPEQ *peqptr;
130	int norm;
131	{
132	tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
133	FP_D(peqptr) = -(DOT(FP_N(peqptr),v0));
134	}
135
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
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,d;
162	int hit;
163
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(t < 0)
170	hit = 0;
171	else
172	hit = 1;
173	if(r)
174	VSUM(r,orig,dir,t);
175
176	if(pd)
177	*pd = t;
178	return(hit);
179	}
180
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 d;
192
193	VSUB(ps,p,c);
194	d = normalize(ps);
195	return(d);
196	}
197
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
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,v0,v1);
215	/* Test the point for sidedness */
216	d = DOT(n,p);
217	if(d > 0.0)
218	return(FALSE);
219	/* Test next edge */
220	VCROSS(n,v1,v2);
221	/* Test the point for sidedness */
222	d = DOT(n,p);
223	if(d > 0.0)
224	return(FALSE);
225	/* Test next edge */
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(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
244	ray_intersect_tri(orig,dir,v0,v1,v2,pt)
245	FVECT orig,dir;
246	FVECT v0,v1,v2;
247	FVECT pt;
248	{
249	FVECT p0,p1,p2,p;
250	FPEQ peq;
251	int type;
252
253	VSUB(p0,v0,orig);
254	VSUB(p1,v1,orig);
255	VSUB(p2,v2,orig);
256
257	if(point_in_stri(p0,p1,p2,dir))
258	{
259	/* Intersect the ray with the triangle plane */
260	tri_plane_equation(v0,v1,v2,&peq,FALSE);
261	return(intersect_ray_plane(orig,dir,peq,NULL,pt));
262	}
263	return(FALSE);
264	}
265
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;
282	FVECT fnear[4],ffar[4];
283	{
284	double height,width;
285	FVECT t,nhv,nvv,ndv;
286	double w2,h2;
287	/* Calculate the x and y dimensions of the near face */
288	VCOPY(nhv,hv);
289	VCOPY(nvv,vv);
290	w2 = normalize(nhv);
291	h2 = normalize(nvv);
292	/* Use similar triangles to calculate the dimensions at z=near */
293	width = near0.5w2;
294	height = near0.5h2;
295
296	VCROSS(ndv,nvv,nhv);
297	/* Calculate the world space points corresponding to the 4 corners
298	of the front face of the view frustum
299	*/
300	fnear[0][0] = widthnhv[0] + heightnvv[0] + near*ndv[0] + vp[0] ;
301	fnear[0][1] = widthnhv[1] + heightnvv[1] + near*ndv[1] + vp[1];
302	fnear[0][2] = widthnhv[2] + heightnvv[2] + near*ndv[2] + vp[2];
303	fnear[1][0] = -widthnhv[0] + heightnvv[0] + near*ndv[0] + vp[0];
304	fnear[1][1] = -widthnhv[1] + heightnvv[1] + near*ndv[1] + vp[1];
305	fnear[1][2] = -widthnhv[2] + heightnvv[2] + near*ndv[2] + vp[2];
306	fnear[2][0] = -widthnhv[0] - heightnvv[0] + near*ndv[0] + vp[0];
307	fnear[2][1] = -widthnhv[1] - heightnvv[1] + near*ndv[1] + vp[1];
308	fnear[2][2] = -widthnhv[2] - heightnvv[2] + near*ndv[2] + vp[2];
309	fnear[3][0] = widthnhv[0] - heightnvv[0] + near*ndv[0] + vp[0];
310	fnear[3][1] = widthnhv[1] - heightnvv[1] + near*ndv[1] + vp[1];
311	fnear[3][2] = widthnhv[2] - heightnvv[2] + near*ndv[2] + vp[2];
312
313	/* Now do the far face */
314	width = far0.5w2;
315	height = far0.5h2;
316	ffar[0][0] = widthnhv[0] + heightnvv[0] + far*ndv[0] + vp[0] ;
317	ffar[0][1] = widthnhv[1] + heightnvv[1] + far*ndv[1] + vp[1] ;
318	ffar[0][2] = widthnhv[2] + heightnvv[2] + far*ndv[2] + vp[2] ;
319	ffar[1][0] = -widthnhv[0] + heightnvv[0] + far*ndv[0] + vp[0] ;
320	ffar[1][1] = -widthnhv[1] + heightnvv[1] + far*ndv[1] + vp[1] ;
321	ffar[1][2] = -widthnhv[2] + heightnvv[2] + far*ndv[2] + vp[2] ;
322	ffar[2][0] = -widthnhv[0] - heightnvv[0] + far*ndv[0] + vp[0] ;
323	ffar[2][1] = -widthnhv[1] - heightnvv[1] + far*ndv[1] + vp[1] ;
324	ffar[2][2] = -widthnhv[2] - heightnvv[2] + far*ndv[2] + vp[2] ;
325	ffar[3][0] = widthnhv[0] - heightnvv[0] + far*ndv[0] + vp[0] ;
326	ffar[3][1] = widthnhv[1] - heightnvv[1] + far*ndv[1] + vp[1] ;
327	ffar[3][2] = widthnhv[2] - heightnvv[2] + far*ndv[2] + vp[2] ;
328	}
329
330
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;
341	double px,py;
342	double coord[3];
343	{
344	double a;
345
346	a = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
347	coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
348	coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
349	coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
350
351	}
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366