--- ray/src/hd/sm_geom.c 1998/08/20 16:47:21 3.2 +++ ray/src/hd/sm_geom.c 2003/02/22 02:07:25 3.13 @@ -1,782 +1,281 @@ -/* Copyright (c) 1998 Silicon Graphics, Inc. */ - #ifndef lint -static char SCCSid[] = "$SunId$ SGI"; +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" -tri_centroid(v0,v1,v2,c) -FVECT v0,v1,v2,c; +/* + * 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]; { -/* Average three triangle vertices to give centroid: return in 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 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); -int -vec3_equal(v1,v2) -FVECT v1,v2; -{ - return(EQUAL(v1[0],v2[0]) && EQUAL(v1[1],v2[1])&& EQUAL(v1[2],v2[2])); + 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); } - -int -convex_angle(v0,v1,v2) -FVECT v0,v1,v2; +/* + * 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; { - FVECT cp01,cp12,cp; - - /* test sign of (v0Xv1)X(v1Xv2). v1 */ - VCROSS(cp01,v0,v1); - VCROSS(cp12,v1,v2); - VCROSS(cp,cp01,cp12); - if(DOT(cp,v1) < 0) - return(FALSE); - return(TRUE); + 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; } -/* calculates the normal of a face contour using Newell's formula. e - 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) : 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; -char norm; +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]); - - + (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,n,nd,norm) - FVECT v0,v1,v2,n; - double *nd; - char norm; +/* + * 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,n,norm); - - *nd = -(DOT(n,v0)); + 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 -point_relative_to_plane(p,n,nd) - FVECT p,n; - double nd; -{ - double d; - - d = p[0]*n[0] + p[1]*n[1] + p[2]*n[2] + nd; - if(d < 0) - return(-1); - if(ZERO(d)) - return(0); - else - return(1); -} - -/* From quad_edge-code */ -int -point_in_circle_thru_origin(p,p0,p1) -FVECT p; -FVECT p0,p1; -{ - - double dp0,dp1; - double dp,det; - - dp0 = DOT_VEC2(p0,p0); - dp1 = DOT_VEC2(p1,p1); - dp = DOT_VEC2(p,p); - - det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1); - - return (det > 0); -} - - - -point_on_sphere(ps,p,c) -FVECT ps,p,c; -{ - VSUB(ps,p,c); - normalize(ps); -} - - -int -intersect_vector_plane(v,plane_n,plane_d,tptr,r) - FVECT v,plane_n; - double plane_d; - double *tptr; +intersect_ray_plane(orig,dir,peq,pd,r) + FVECT orig,dir; + FPEQ peq; + double *pd; FVECT r; { - double t; + double t,d; int hit; - /* - Plane is Ax + By + Cz +D = 0: - plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D - */ - /* line is l = p1 + (p2-p1)t, p1=origin */ + d = DOT(FP_N(peq),dir); + if(ZERO(d)) + return(0); + t = -(DOT(FP_N(peq),orig) + FP_D(peq))/d; - /* Solve for t: */ - t = plane_d/-(DOT(plane_n,v)); - if(t >0 || ZERO(t)) - hit = 1; - else + if(t < 0) hit = 0; - r[0] = v[0]*t; - r[1] = v[1]*t; - r[2] = v[2]*t; - if(tptr) - *tptr = t; - return(hit); -} - -int -intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r) - FVECT orig,dir; - FVECT plane_n; - double plane_d; - double *pd; - FVECT r; -{ - double t; - int hit; - /* - Plane is Ax + By + Cz +D = 0: - plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D - */ - /* 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: */ - t = -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir)); - if(ZERO(t) || t >0) - hit = 1; else - hit = 0; + hit = 1; + if(r) + VSUM(r,orig,dir,t); - VSUM(r,orig,dir,t); - if(pd) *pd = t; return(hit); } - -int -point_in_cone(p,p0,p1,p2) -FVECT p; -FVECT p0,p1,p2; +/* + * 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; { - FVECT n; - FVECT np,x_axis,y_axis; - double d1,d2,d; - - /* Find the equation of the circle defined by the intersection - of the cone with the plane defined by p1,p2,p3- project p into - that plane and do an in-circle test in the plane - */ - - /* find the equation of the plane defined by p1-p3 */ - tri_plane_equation(p0,p1,p2,n,&d,FALSE); - - /* define a coordinate system on the plane: the x axis is in - the direction of np2-np1, and the y axis is calculated from - n cross x-axis - */ - /* Project p onto the plane */ - if(!intersect_vector_plane(p,n,d,NULL,np)) - return(FALSE); - - /* create coordinate system on plane: p2-p1 defines the x_axis*/ - VSUB(x_axis,p1,p0); - normalize(x_axis); - /* The y axis is */ - VCROSS(y_axis,n,x_axis); - normalize(y_axis); - - VSUB(p1,p1,p0); - VSUB(p2,p2,p0); - VSUB(np,np,p0); - - p1[0] = VLEN(p1); - p1[1] = 0; - - d1 = DOT(p2,x_axis); - d2 = DOT(p2,y_axis); - p2[0] = d1; - p2[1] = d2; - - d1 = DOT(np,x_axis); - d2 = DOT(np,y_axis); - np[0] = d1; - np[1] = d2; - - /* perform the in-circle test in the new coordinate system */ - return(point_in_circle_thru_origin(np,p1,p2)); + 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 -test_point_against_spherical_tri(v0,v1,v2,p,n,nset,which,sides) +point_in_stri(v0,v1,v2,p) FVECT v0,v1,v2,p; -FVECT n[3]; -char *nset; -char *which; -char sides[3]; - { - float d; - - /* Find the normal to the triangle ORIGIN:v0:v1 */ - if(!NTH_BIT(*nset,0)) - { - VCROSS(n[0],v1,v0); - SET_NTH_BIT(*nset,0); - } - /* Test the point for sidedness */ - d = DOT(n[0],p); - - if(ZERO(d)) - sides[0] = GT_EDGE; - else - if(d > 0) - { - sides[0] = GT_OUT; - sides[1] = sides[2] = GT_INVALID; - return(FALSE); - } - else - sides[0] = GT_INTERIOR; - - /* Test next edge */ - if(!NTH_BIT(*nset,1)) - { - VCROSS(n[1],v2,v1); - SET_NTH_BIT(*nset,1); - } - /* Test the point for sidedness */ - d = DOT(n[1],p); - if(ZERO(d)) - { - sides[1] = GT_EDGE; - /* If on plane 0-and on plane 1: lies on edge */ - if(sides[0] == GT_EDGE) - { - *which = 1; - sides[2] = GT_INVALID; - return(GT_EDGE); - } - } - else if(d > 0) - { - sides[1] = GT_OUT; - sides[2] = GT_INVALID; - return(FALSE); - } - else - sides[1] = GT_INTERIOR; - /* Test next edge */ - if(!NTH_BIT(*nset,2)) - { - - VCROSS(n[2],v0,v2); - SET_NTH_BIT(*nset,2); - } - /* Test the point for sidedness */ - d = DOT(n[2],p); - if(ZERO(d)) - { - sides[2] = GT_EDGE; - - /* If on plane 0 and 2: lies on edge 0*/ - if(sides[0] == GT_EDGE) - { - *which = 0; - return(GT_EDGE); - } - /* If on plane 1 and 2: lies on edge 2*/ - if(sides[1] == GT_EDGE) - { - *which = 2; - return(GT_EDGE); - } - /* otherwise: on face 2 */ - else - { - *which = 2; - return(GT_FACE); - } - } - else if(d > 0) - { - sides[2] = GT_OUT; - return(FALSE); - } - /* If on edge */ - else - sides[2] = GT_INTERIOR; - - /* If on plane 0 only: on face 0 */ - if(sides[0] == GT_EDGE) - { - *which = 0; - return(GT_FACE); - } - /* If on plane 1 only: on face 1 */ - if(sides[1] == GT_EDGE) - { - *which = 1; - return(GT_FACE); - } - /* Must be interior to the pyramid */ - return(GT_INTERIOR); -} - - - - -int -test_single_point_against_spherical_tri(v0,v1,v2,p,which) -FVECT v0,v1,v2,p; -char *which; -{ - float d; + double d; FVECT n; - char sides[3]; - /* First test if point coincides with any of the vertices */ - if(EQUAL_VEC3(p,v0)) - { - *which = 0; - return(GT_VERTEX); - } - if(EQUAL_VEC3(p,v1)) - { - *which = 1; - return(GT_VERTEX); - } - if(EQUAL_VEC3(p,v2)) - { - *which = 2; - return(GT_VERTEX); - } - VCROSS(n,v1,v0); + VCROSS(n,v0,v1); /* Test the point for sidedness */ d = DOT(n,p); - if(ZERO(d)) - sides[0] = GT_EDGE; - else - if(d > 0) - return(FALSE); - else - sides[0] = GT_INTERIOR; + if(d > 0.0) + return(FALSE); /* Test next edge */ - VCROSS(n,v2,v1); + VCROSS(n,v1,v2); /* Test the point for sidedness */ d = DOT(n,p); - if(ZERO(d)) - { - sides[1] = GT_EDGE; - /* If on plane 0-and on plane 1: lies on edge */ - if(sides[0] == GT_EDGE) - { - *which = 1; - return(GT_VERTEX); - } - } - else if(d > 0) + if(d > 0.0) return(FALSE); - else - sides[1] = GT_INTERIOR; - /* Test next edge */ - VCROSS(n,v0,v2); + VCROSS(n,v2,v0); /* Test the point for sidedness */ d = DOT(n,p); - if(ZERO(d)) - { - sides[2] = GT_EDGE; - - /* If on plane 0 and 2: lies on edge 0*/ - if(sides[0] == GT_EDGE) - { - *which = 0; - return(GT_VERTEX); - } - /* If on plane 1 and 2: lies on edge 2*/ - if(sides[1] == GT_EDGE) - { - *which = 2; - return(GT_VERTEX); - } - /* otherwise: on face 2 */ - else - { - return(GT_FACE); - } - } - else if(d > 0) + if(d > 0.0) return(FALSE); /* Must be interior to the pyramid */ - return(GT_FACE); + 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 -test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides) -FVECT t0,t1,t2,p0,p1,p2; -char *nset; -FVECT n[3]; -FVECT avg; -char pt_sides[3][3]; - -{ - char below_plane[3],on_edge,test; - char which; - - SUM_3VEC3(avg,t0,t1,t2); - on_edge = 0; - *nset = 0; - /* Test vertex v[i] against triangle j*/ - /* Check if v[i] lies below plane defined by avg of 3 vectors - defining triangle - */ - - /* test point 0 */ - if(DOT(avg,p0) < 0) - below_plane[0] = 1; - else - below_plane[0]=0; - /* Test if b[i] lies in or on triangle a */ - test = test_point_against_spherical_tri(t0,t1,t2,p0, - n,nset,&which,pt_sides[0]); - /* If pts[i] is interior: done */ - if(!below_plane[0]) - { - if(test == GT_INTERIOR) - return(TRUE); - /* Remember if b[i] fell on one of the 3 defining planes */ - if(test) - on_edge++; - } - /* Now test point 1*/ - - if(DOT(avg,p1) < 0) - below_plane[1] = 1; - else - below_plane[1]=0; - /* Test if b[i] lies in or on triangle a */ - test = test_point_against_spherical_tri(t0,t1,t2,p1, - n,nset,&which,pt_sides[1]); - /* If pts[i] is interior: done */ - if(!below_plane[1]) - { - if(test == GT_INTERIOR) - return(TRUE); - /* Remember if b[i] fell on one of the 3 defining planes */ - if(test) - on_edge++; - } - - /* Now test point 2 */ - if(DOT(avg,p2) < 0) - below_plane[2] = 1; - else - below_plane[2]=0; - /* Test if b[i] lies in or on triangle a */ - test = test_point_against_spherical_tri(t0,t1,t2,p2, - n,nset,&which,pt_sides[2]); - - /* If pts[i] is interior: done */ - if(!below_plane[2]) - { - if(test == GT_INTERIOR) - return(TRUE); - /* Remember if b[i] fell on one of the 3 defining planes */ - if(test) - on_edge++; - } - - /* If all three points below separating plane: trivial reject */ - if(below_plane[0] && below_plane[1] && below_plane[2]) - return(FALSE); - /* Accept if all points lie on a triangle vertex/edge edge- accept*/ - if(on_edge == 3) - return(TRUE); - /* Now check vertices in a against triangle b */ - return(FALSE); -} - - -set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n) - FVECT t0,t1,t2,p0,p1,p2; - char test[3]; - char sides[3][3]; - char nset; - FVECT n[3]; -{ - char t; - double d; - - - /* p=0 */ - test[0] = 0; - if(sides[0][0] == GT_INVALID) - { - if(!NTH_BIT(nset,0)) - VCROSS(n[0],t1,t0); - /* Test the point for sidedness */ - d = DOT(n[0],p0); - if(d >= 0) - SET_NTH_BIT(test[0],0); - } - else - if(sides[0][0] != GT_INTERIOR) - SET_NTH_BIT(test[0],0); - - if(sides[0][1] == GT_INVALID) - { - if(!NTH_BIT(nset,1)) - VCROSS(n[1],t2,t1); - /* Test the point for sidedness */ - d = DOT(n[1],p0); - if(d >= 0) - SET_NTH_BIT(test[0],1); - } - else - if(sides[0][1] != GT_INTERIOR) - SET_NTH_BIT(test[0],1); - - if(sides[0][2] == GT_INVALID) - { - if(!NTH_BIT(nset,2)) - VCROSS(n[2],t0,t2); - /* Test the point for sidedness */ - d = DOT(n[2],p0); - if(d >= 0) - SET_NTH_BIT(test[0],2); - } - else - if(sides[0][2] != GT_INTERIOR) - SET_NTH_BIT(test[0],2); - - /* p=1 */ - test[1] = 0; - /* t=0*/ - if(sides[1][0] == GT_INVALID) - { - if(!NTH_BIT(nset,0)) - VCROSS(n[0],t1,t0); - /* Test the point for sidedness */ - d = DOT(n[0],p1); - if(d >= 0) - SET_NTH_BIT(test[1],0); - } - else - if(sides[1][0] != GT_INTERIOR) - SET_NTH_BIT(test[1],0); - - /* t=1 */ - if(sides[1][1] == GT_INVALID) - { - if(!NTH_BIT(nset,1)) - VCROSS(n[1],t2,t1); - /* Test the point for sidedness */ - d = DOT(n[1],p1); - if(d >= 0) - SET_NTH_BIT(test[1],1); - } - else - if(sides[1][1] != GT_INTERIOR) - SET_NTH_BIT(test[1],1); - - /* t=2 */ - if(sides[1][2] == GT_INVALID) - { - if(!NTH_BIT(nset,2)) - VCROSS(n[2],t0,t2); - /* Test the point for sidedness */ - d = DOT(n[2],p1); - if(d >= 0) - SET_NTH_BIT(test[1],2); - } - else - if(sides[1][2] != GT_INTERIOR) - SET_NTH_BIT(test[1],2); - - /* p=2 */ - test[2] = 0; - /* t = 0 */ - if(sides[2][0] == GT_INVALID) - { - if(!NTH_BIT(nset,0)) - VCROSS(n[0],t1,t0); - /* Test the point for sidedness */ - d = DOT(n[0],p2); - if(d >= 0) - SET_NTH_BIT(test[2],0); - } - else - if(sides[2][0] != GT_INTERIOR) - SET_NTH_BIT(test[2],0); - /* t=1 */ - if(sides[2][1] == GT_INVALID) - { - if(!NTH_BIT(nset,1)) - VCROSS(n[1],t2,t1); - /* Test the point for sidedness */ - d = DOT(n[1],p2); - if(d >= 0) - SET_NTH_BIT(test[2],1); - } - else - if(sides[2][1] != GT_INTERIOR) - SET_NTH_BIT(test[2],1); - /* t=2 */ - if(sides[2][2] == GT_INVALID) - { - if(!NTH_BIT(nset,2)) - VCROSS(n[2],t0,t2); - /* Test the point for sidedness */ - d = DOT(n[2],p2); - if(d >= 0) - SET_NTH_BIT(test[2],2); - } - else - if(sides[2][2] != GT_INTERIOR) - SET_NTH_BIT(test[2],2); -} - - -int -spherical_tri_intersect(a1,a2,a3,b1,b2,b3) -FVECT a1,a2,a3,b1,b2,b3; -{ - char which,test,n_set[2]; - char sides[2][3][3],i,j,inext,jnext; - char tests[2][3]; - FVECT n[2][3],p,avg[2]; - - /* Test the vertices of triangle a against the pyramid formed by triangle - b and the origin. If any vertex of a is interior to triangle b, or - if all 3 vertices of a are ON the edges of b,return TRUE. Remember - the results of the edge normal and sidedness tests for later. - */ - if(test_vertices_for_tri_inclusion(a1,a2,a3,b1,b2,b3, - &(n_set[0]),n[0],avg[0],sides[1])) - return(TRUE); - - if(test_vertices_for_tri_inclusion(b1,b2,b3,a1,a2,a3, - &(n_set[1]),n[1],avg[1],sides[0])) - return(TRUE); - - - set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]); - if(tests[0][0]&tests[0][1]&tests[0][2]) - return(FALSE); - - set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]); - if(tests[1][0]&tests[1][1]&tests[1][2]) - return(FALSE); - - for(j=0; j < 3;j++) - { - jnext = (j+1)%3; - /* IF edge b doesnt cross any great circles of a, punt */ - if(tests[1][j] & tests[1][jnext]) - continue; - for(i=0;i<3;i++) - { - inext = (i+1)%3; - /* IF edge a doesnt cross any great circles of b, punt */ - if(tests[0][i] & tests[0][inext]) - continue; - /* Now find the great circles that cross and test */ - if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j))) - && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i))) - { - VCROSS(p,n[0][i],n[1][j]); - - /* If zero cp= done */ - if(ZERO_VEC3(p)) - continue; - /* check above both planes */ - if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0) - { - NEGATE_VEC3(p); - if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0) - continue; - } - return(TRUE); - } - } - } - return(FALSE); -} - -int -ray_intersect_tri(orig,dir,v0,v1,v2,pt,wptr) +ray_intersect_tri(orig,dir,v0,v1,v2,pt) FVECT orig,dir; FVECT v0,v1,v2; FVECT pt; -char *wptr; { - FVECT p0,p1,p2,p,n; - char type,which; - double pd; - - point_on_sphere(p0,v0,orig); - point_on_sphere(p1,v1,orig); - point_on_sphere(p2,v2,orig); - type = test_single_point_against_spherical_tri(p0,p1,p2,dir,&which); + FVECT p0,p1,p2,p; + FPEQ peq; + int type; - if(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,n,&pd,FALSE); - intersect_ray_plane(orig,dir,n,pd,NULL,pt); + tri_plane_equation(v0,v1,v2,&peq,FALSE); + return(intersect_ray_plane(orig,dir,peq,NULL,pt)); } - if(wptr) - *wptr = which; - - return(type); + 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; @@ -786,10 +285,6 @@ FVECT fnear[4],ffar[4]; FVECT t,nhv,nvv,ndv; double w2,h2; /* Calculate the x and y dimensions of the near face */ - /* hv and vv are the horizontal and vertical vectors in the - view frame-the magnitude is the dimension of the front frustum - face at z =1 - */ VCOPY(nhv,hv); VCOPY(nvv,vv); w2 = normalize(nhv); @@ -833,58 +328,13 @@ FVECT fnear[4],ffar[4]; } - - -int -spherical_polygon_edge_intersect(a0,a1,b0,b1) -FVECT a0,a1,b0,b1; -{ - FVECT na,nb,avga,avgb,p; - double d; - int sb0,sb1,sa0,sa1; - - /* First test if edge b straddles great circle of a */ - VCROSS(na,a0,a1); - d = DOT(na,b0); - sb0 = ZERO(d)?0:(d<0)? -1: 1; - d = DOT(na,b1); - sb1 = ZERO(d)?0:(d<0)? -1: 1; - /* edge b entirely on one side of great circle a: edges cannot intersect*/ - if(sb0*sb1 > 0) - return(FALSE); - /* test if edge a straddles great circle of b */ - VCROSS(nb,b0,b1); - d = DOT(nb,a0); - sa0 = ZERO(d)?0:(d<0)? -1: 1; - d = DOT(nb,a1); - sa1 = ZERO(d)?0:(d<0)? -1: 1; - /* edge a entirely on one side of great circle b: edges cannot intersect*/ - if(sa0*sa1 > 0) - return(FALSE); - - /* Find one of intersection points of the great circles */ - VCROSS(p,na,nb); - /* If they lie on same great circle: call an intersection */ - if(ZERO_VEC3(p)) - return(TRUE); - - VADD(avga,a0,a1); - VADD(avgb,b0,b1); - if(DOT(avga,p) < 0 || DOT(avgb,p) < 0) - { - NEGATE_VEC3(p); - if(DOT(avga,p) < 0 || DOT(avgb,p) < 0) - return(FALSE); - } - if((!sb0 || !sb1) && (!sa0 || !sa1)) - return(FALSE); - return(TRUE); -} - - - -/* Find the normalized barycentric coordinates of p relative to - * triangle v0,v1,v2. Return result in coord +/* + * 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; @@ -896,91 +346,13 @@ double coord[3]; 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] = 1.0 - coord[0] - coord[1]; + coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a; } -int -bary2d_child(coord) -double coord[3]; -{ - int i; - /* First check if one of the original vertices */ - for(i=0;i<3;i++) - if(EQUAL(coord[i],1.0)) - return(i); - /* Check if one of the new vertices: for all return center child */ - if(ZERO(coord[0]) && EQUAL(coord[1],0.5)) - { - coord[0] = 1.0f; - coord[1] = 0.0f; - coord[2] = 0.0f; - return(3); - } - if(ZERO(coord[1]) && EQUAL(coord[0],0.5)) - { - coord[0] = 0.0f; - coord[1] = 1.0f; - coord[2] = 0.0f; - return(3); - } - if(ZERO(coord[2]) && EQUAL(coord[0],0.5)) - { - coord[0] = 0.0f; - coord[1] = 0.0f; - coord[2] = 1.0f; - return(3); - } - /* Otherwise return child */ - if(coord[0] > 0.5) - { - /* update bary for child */ - coord[0] = 2.0*coord[0]- 1.0; - coord[1] *= 2.0; - coord[2] *= 2.0; - return(0); - } - else - if(coord[1] > 0.5) - { - coord[0] *= 2.0; - coord[1] = 2.0*coord[1]- 1.0; - coord[2] *= 2.0; - return(1); - } - else - if(coord[2] > 0.5) - { - coord[0] *= 2.0; - coord[1] *= 2.0; - coord[2] = 2.0*coord[2]- 1.0; - return(2); - } - else - { - coord[0] = 1.0 - 2.0*coord[0]; - coord[1] = 1.0 - 2.0*coord[1]; - coord[2] = 1.0 - 2.0*coord[2]; - return(3); - } -} - -int -max_index(v) -FVECT v; -{ - double a,b,c; - int i; - - a = fabs(v[0]); - b = fabs(v[1]); - c = fabs(v[2]); - i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2); - return(i); -}