src/hd/sm_geom.c

#ifndef lint
static const char       RCSid[] = "$Id: sm_geom.c,v 3.13 2003/02/22 02:07:25 greg Exp $";
#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.14
Committed:	Fri Nov 5 03:31:37 2004 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 3.13:	+1 -1 lines
State:	*FILE REMOVED*
Log Message:	Removed unused programs and files from distribution (sources to CVS attic)
#	User	Rev	Content
1	gwlarson	3.1	#ifndef lint
2	greg	3.14	static const char RCSid[] = "$Id: sm_geom.c,v 3.13 2003/02/22 02:07:25 greg Exp $";
3	gwlarson	3.1	#endif
4			/*
5			* sm_geom.c
6	greg	3.13	* some geometric utility routines
7	gwlarson	3.1	*/
8
9			#include "standard.h"
10			#include "sm_geom.h"
11
12	greg	3.13	/*
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	gwlarson	3.1	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	greg	3.13	/*
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	gwlarson	3.1	*/
97			double
98			tri_normal(v0,v1,v2,n,norm)
99			FVECT v0,v1,v2,n;
100	gwlarson	3.4	int norm;
101	gwlarson	3.1	{
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	gwlarson	3.9	(v2[2] - v0[2]) * (v2[0] + v0[0]);
110	gwlarson	3.1	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	greg	3.13	/*
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	gwlarson	3.7	tri_plane_equation(v0,v1,v2,peqptr,norm)
128	gwlarson	3.12	FVECT v0,v1,v2;
129	gwlarson	3.7	FPEQ *peqptr;
130	gwlarson	3.4	int norm;
131	gwlarson	3.1	{
132	gwlarson	3.7	tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
133			FP_D(peqptr) = -(DOT(FP_N(peqptr),v0));
134	gwlarson	3.1	}
135
136	greg	3.13	/*
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	gwlarson	3.1	int
155	gwlarson	3.7	intersect_ray_plane(orig,dir,peq,pd,r)
156	gwlarson	3.1	FVECT orig,dir;
157	gwlarson	3.7	FPEQ peq;
158	gwlarson	3.1	double *pd;
159			FVECT r;
160			{
161	gwlarson	3.8	double t,d;
162	gwlarson	3.1	int hit;
163	greg	3.13
164	gwlarson	3.8	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	gwlarson	3.4	hit = 0;
171			else
172	gwlarson	3.1	hit = 1;
173	gwlarson	3.4	if(r)
174			VSUM(r,orig,dir,t);
175
176			if(pd)
177			*pd = t;
178			return(hit);
179			}
180
181	greg	3.13	/*
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	gwlarson	3.12	double
188			point_on_sphere(ps,p,c)
189			FVECT ps,p,c;
190			{
191			double d;
192	greg	3.13
193			VSUB(ps,p,c);
194			d = normalize(ps);
195			return(d);
196	gwlarson	3.12	}
197	gwlarson	3.1
198	greg	3.13	/*
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	gwlarson	3.1	int
208	gwlarson	3.12	point_in_stri(v0,v1,v2,p)
209			FVECT v0,v1,v2,p;
210	gwlarson	3.1	{
211	gwlarson	3.12	double d;
212			FVECT n;
213	gwlarson	3.9
214	gwlarson	3.12	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	greg	3.13	return(TRUE);
233	gwlarson	3.12	}
234	gwlarson	3.1
235	greg	3.13	/*
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	gwlarson	3.1	int
244	gwlarson	3.4	ray_intersect_tri(orig,dir,v0,v1,v2,pt)
245	gwlarson	3.1	FVECT orig,dir;
246			FVECT v0,v1,v2;
247			FVECT pt;
248			{
249	gwlarson	3.7	FVECT p0,p1,p2,p;
250			FPEQ peq;
251	gwlarson	3.4	int type;
252
253	gwlarson	3.5	VSUB(p0,v0,orig);
254			VSUB(p1,v1,orig);
255			VSUB(p2,v2,orig);
256
257	gwlarson	3.4	if(point_in_stri(p0,p1,p2,dir))
258	gwlarson	3.1	{
259			/* Intersect the ray with the triangle plane */
260	gwlarson	3.7	tri_plane_equation(v0,v1,v2,&peq,FALSE);
261			return(intersect_ray_plane(orig,dir,peq,NULL,pt));
262	gwlarson	3.1	}
263	gwlarson	3.4	return(FALSE);
264	gwlarson	3.1	}
265
266	greg	3.13	/*
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	gwlarson	3.1	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	gwlarson	3.2
331	greg	3.13	/*
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	gwlarson	3.2	*/
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	gwlarson	3.6	coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
350	gwlarson	3.2
351			}
352
353	gwlarson	3.4
354
355	gwlarson	3.2
356
357	gwlarson	3.6
358
359
360
361	greg	3.13
362
363
364	gwlarson	3.2
365
366	gwlarson	3.1