--- ray/src/hd/sm_geom.c 1998/08/19 17:45:24 3.1 +++ ray/src/hd/sm_geom.c 1999/01/05 16:52:38 3.10 @@ -27,33 +27,44 @@ FVECT v1,v2; { return(EQUAL(v1[0],v2[0]) && EQUAL(v1[1],v2[1])&& EQUAL(v1[2],v2[2])); } +#if 0 +extern FVECT Norm[500]; +extern int Ncnt; +#endif - int convex_angle(v0,v1,v2) FVECT v0,v1,v2; { - FVECT cp01,cp12,cp; - + FVECT cp,cp01,cp12,v10,v02; + double dp; + /* 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); + + dp = DOT(cp,v1); +#if 0 + VCOPY(Norm[Ncnt++],cp01); + VCOPY(Norm[Ncnt++],cp12); + VCOPY(Norm[Ncnt++],cp); +#endif + if(ZERO(dp) || dp < 0.0) + return(FALSE); return(TRUE); } /* calculates the normal of a face contour using Newell's formula. e - a = SUMi (yi - yi+1)(zi + zi+1) + 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; @@ -63,8 +74,7 @@ char norm; 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]) + @@ -72,7 +82,6 @@ char norm; if(!norm) return(0); - mag = normalize(n); @@ -80,32 +89,16 @@ char norm; } -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) + FVECT v0,v1,v2; + FPEQ *peqptr; + int norm; { - tri_normal(v0,v1,v2,n,norm); + tri_normal(v0,v1,v2,FP_N(*peqptr),norm); - *nd = -(DOT(n,v0)); + FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0)); } -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) @@ -122,27 +115,34 @@ FVECT p0,p1; det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1); - return (det > 0); + return (det > (1e-10)); } - +double point_on_sphere(ps,p,c) FVECT ps,p,c; { + double d; VSUB(ps,p,c); - normalize(ps); + d= normalize(ps); + return(d); } +/* returns TRUE if ray from origin in direction v intersects plane defined + * by normal plane_n, and plane_d. If plane is not parallel- returns + * intersection point if r != NULL. If tptr!= NULL returns value of + * t, if parallel, returns t=FHUGE + */ int -intersect_vector_plane(v,plane_n,plane_d,pd,r) - FVECT v,plane_n; - double plane_d; - double *pd; +intersect_vector_plane(v,peq,tptr,r) + FVECT v; + FPEQ peq; + double *tptr; FVECT r; { - double t; + double t,d; int hit; /* Plane is Ax + By + Cz +D = 0: @@ -152,44 +152,149 @@ intersect_vector_plane(v,plane_n,plane_d,pd,r) /* line is l = p1 + (p2-p1)t, p1=origin */ /* Solve for t: */ - t = plane_d/-(DOT(plane_n,v)); - if(t >0 || ZERO(t)) - hit = 1; - else + d = -(DOT(FP_N(peq),v)); + if(ZERO(d)) + { + t = FHUGE; + hit = 0; + } + else + { + t = FP_D(peq)/d; + if(t < 0 ) + hit = 0; + else + hit = 1; + if(r) + { + 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,peq,pd,r) + FVECT orig,dir; + FPEQ peq; + double *pd; + FVECT r; +{ + double t,d; + 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: */ + 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; - r[0] = v[0]*t; - r[1] = v[1]*t; - r[2] = v[2]*t; + else + hit = 1; + + if(r) + VSUM(r,orig,dir,t); + if(pd) *pd = t; return(hit); } + int -intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r) +intersect_ray_oplane(orig,dir,n,pd,r) FVECT orig,dir; - FVECT plane_n; - double plane_d; + FVECT n; 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 */ - /* 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 */ + /* 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; + d= DOT(n,dir); + if(ZERO(d)) + return(0); + t = -(DOT(n,orig))/d; + if(t < 0) + hit = 0; else + hit = 1; + + if(r) + VSUM(r,orig,dir,t); + + if(pd) + *pd = t; + return(hit); +} + +/* Assumption: know crosses plane:dont need to check for 'on' case */ +intersect_edge_coord_plane(v0,v1,w,r) +FVECT v0,v1; +int w; +FVECT r; +{ + FVECT dv; + int wnext; + double t; + + VSUB(dv,v1,v0); + t = -v0[w]/dv[w]; + r[w] = 0.0; + wnext = (w+1)%3; + r[wnext] = v0[wnext] + dv[wnext]*t; + wnext = (w+2)%3; + r[wnext] = v0[wnext] + dv[wnext]*t; +} + +int +intersect_edge_plane(e0,e1,peq,pd,r) + FVECT e0,e1; + FPEQ peq; + double *pd; + FVECT r; +{ + double t; + int hit; + FVECT d; + /* + 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: */ + VSUB(d,e1,e0); + t = -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d)); + if(t < 0) hit = 0; + else + hit = 1; - VSUM(r,orig,dir,t); + VSUM(r,e0,d,t); if(pd) *pd = t; @@ -202,38 +307,39 @@ point_in_cone(p,p0,p1,p2) FVECT p; FVECT p0,p1,p2; { - FVECT n; FVECT np,x_axis,y_axis; - double d1,d2,d; + double d1,d2; + FPEQ peq; /* 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); + /* find the equation of the plane defined by p0-p2 */ + tri_plane_equation(p0,p1,p2,&peq,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)) + /* NOTE: check this: does sideness check?*/ + if(!intersect_vector_plane(p,peq,NULL,np)) return(FALSE); - /* create coordinate system on plane: p2-p1 defines the x_axis*/ + /* create coordinate system on plane: p1-p0 defines the x_axis*/ VSUB(x_axis,p1,p0); normalize(x_axis); /* The y axis is */ - VCROSS(y_axis,n,x_axis); + VCROSS(y_axis,FP_N(peq),x_axis); normalize(y_axis); VSUB(p1,p1,p0); VSUB(p2,p2,p0); VSUB(np,np,p0); - p1[0] = VLEN(p1); + p1[0] = DOT(p1,x_axis); p1[1] = 0; d1 = DOT(p2,x_axis); @@ -251,33 +357,28 @@ FVECT p0,p1,p2; } int -test_point_against_spherical_tri(v0,v1,v2,p,n,nset,which,sides) +point_set_in_stri(v0,v1,v2,p,n,nset,sides) FVECT v0,v1,v2,p; FVECT n[3]; -char *nset; -char *which; -char sides[3]; +int *nset; +int sides[3]; { - float d; - + double d; /* Find the normal to the triangle ORIGIN:v0:v1 */ if(!NTH_BIT(*nset,0)) { - VCROSS(n[0],v1,v0); + VCROSS(n[0],v0,v1); 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); + if(d > 0.0) + { + sides[0] = GT_OUT; + sides[1] = sides[2] = GT_INVALID; + return(FALSE); } else sides[0] = GT_INTERIOR; @@ -285,24 +386,13 @@ char sides[3]; /* Test next edge */ if(!NTH_BIT(*nset,1)) { - VCROSS(n[1],v2,v1); + VCROSS(n[1],v1,v2); SET_NTH_BIT(*nset,1); } /* Test the point for sidedness */ d = DOT(n[1],p); - if(ZERO(d)) + if(d > 0.0) { - 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); @@ -312,163 +402,67 @@ char sides[3]; /* Test next edge */ if(!NTH_BIT(*nset,2)) { - - VCROSS(n[2],v0,v2); + VCROSS(n[2],v2,v0); SET_NTH_BIT(*nset,2); } /* Test the point for sidedness */ d = DOT(n[2],p); - if(ZERO(d)) + if(d > 0.0) { - 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); - } + sides[2] = GT_OUT; + return(FALSE); } - 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) +point_in_stri(v0,v1,v2,p) 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(GT_INTERIOR); } int -test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides) +vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides) FVECT t0,t1,t2,p0,p1,p2; -char *nset; +int *nset; FVECT n[3]; FVECT avg; -char pt_sides[3][3]; +int pt_sides[3][3]; { - char below_plane[3],on_edge,test; - char which; + int below_plane[3],test; 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 @@ -476,66 +470,51 @@ char pt_sides[3][3]; */ /* test point 0 */ - if(DOT(avg,p0) < 0) + if(DOT(avg,p0) < 0.0) below_plane[0] = 1; else - below_plane[0]=0; + 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]); + test = point_set_in_stri(t0,t1,t2,p0,n,nset,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) + if(DOT(avg,p1) < 0.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]); + test = point_set_in_stri(t0,t1,t2,p1,n,nset,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) + if(DOT(avg,p2) < 0.0) below_plane[2] = 1; else - below_plane[2]=0; + 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]); + test = point_set_in_stri(t0,t1,t2,p2,n,nset,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); } @@ -543,12 +522,12 @@ char pt_sides[3][3]; 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; + int test[3]; + int sides[3][3]; + int nset; FVECT n[3]; { - char t; + int t; double d; @@ -557,10 +536,10 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, if(sides[0][0] == GT_INVALID) { if(!NTH_BIT(nset,0)) - VCROSS(n[0],t1,t0); + VCROSS(n[0],t0,t1); /* Test the point for sidedness */ d = DOT(n[0],p0); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[0],0); } else @@ -570,10 +549,10 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, if(sides[0][1] == GT_INVALID) { if(!NTH_BIT(nset,1)) - VCROSS(n[1],t2,t1); + VCROSS(n[1],t1,t2); /* Test the point for sidedness */ d = DOT(n[1],p0); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[0],1); } else @@ -583,10 +562,10 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, if(sides[0][2] == GT_INVALID) { if(!NTH_BIT(nset,2)) - VCROSS(n[2],t0,t2); + VCROSS(n[2],t2,t0); /* Test the point for sidedness */ d = DOT(n[2],p0); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[0],2); } else @@ -599,10 +578,10 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, if(sides[1][0] == GT_INVALID) { if(!NTH_BIT(nset,0)) - VCROSS(n[0],t1,t0); + VCROSS(n[0],t0,t1); /* Test the point for sidedness */ d = DOT(n[0],p1); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[1],0); } else @@ -613,10 +592,10 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, if(sides[1][1] == GT_INVALID) { if(!NTH_BIT(nset,1)) - VCROSS(n[1],t2,t1); + VCROSS(n[1],t1,t2); /* Test the point for sidedness */ d = DOT(n[1],p1); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[1],1); } else @@ -627,10 +606,10 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, if(sides[1][2] == GT_INVALID) { if(!NTH_BIT(nset,2)) - VCROSS(n[2],t0,t2); + VCROSS(n[2],t2,t0); /* Test the point for sidedness */ d = DOT(n[2],p1); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[1],2); } else @@ -643,10 +622,10 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, if(sides[2][0] == GT_INVALID) { if(!NTH_BIT(nset,0)) - VCROSS(n[0],t1,t0); + VCROSS(n[0],t0,t1); /* Test the point for sidedness */ d = DOT(n[0],p2); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[2],0); } else @@ -656,10 +635,10 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, if(sides[2][1] == GT_INVALID) { if(!NTH_BIT(nset,1)) - VCROSS(n[1],t2,t1); + VCROSS(n[1],t1,t2); /* Test the point for sidedness */ d = DOT(n[1],p2); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[2],1); } else @@ -669,10 +648,10 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, if(sides[2][2] == GT_INVALID) { if(!NTH_BIT(nset,2)) - VCROSS(n[2],t0,t2); + VCROSS(n[2],t2,t0); /* Test the point for sidedness */ d = DOT(n[2],p2); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[2],2); } else @@ -682,33 +661,47 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, int -spherical_tri_intersect(a1,a2,a3,b1,b2,b3) +stri_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]; + int which,test,n_set[2]; + int sides[2][3][3],i,j,inext,jnext; + int tests[2][3]; + FVECT n[2][3],p,avg[2],t1,t2,t3; +#ifdef DEBUG + tri_normal(b1,b2,b3,p,FALSE); + if(DOT(p,b1) > 0) + { + VCOPY(t3,b1); + VCOPY(t2,b2); + VCOPY(t1,b3); + } + else +#endif + { + VCOPY(t1,b1); + VCOPY(t2,b2); + VCOPY(t3,b3); + } + /* 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])) + if(vertices_in_stri(a1,a2,a3,t1,t2,t3,&(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])) + if(vertices_in_stri(t1,t2,t3,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]); + set_sidedness_tests(t1,t2,t3,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]); + set_sidedness_tests(a1,a2,a3,t1,t2,t3,tests[1],sides[1],n_set[0],n[0]); if(tests[1][0]&tests[1][1]&tests[1][2]) return(FALSE); @@ -748,31 +741,26 @@ FVECT a1,a2,a3,b1,b2,b3; } 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); } @@ -831,11 +819,49 @@ FVECT fnear[4],ffar[4]; ffar[3][2] = width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ; } +int +max_index(v,r) +FVECT v; +double *r; +{ + double p[3]; + int i; + p[0] = fabs(v[0]); + p[1] = fabs(v[1]); + p[2] = fabs(v[2]); + i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2); + if(r) + *r = p[i]; + return(i); +} +int +closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id) +FVECT p0,p1,p2,p; +int p0id,p1id,p2id; +{ + double d,d1; + int i; + + d = DIST_SQ(p,p0); + d1 = DIST_SQ(p,p1); + if(d < d1) + { + d1 = DIST_SQ(p,p2); + i = (d1 < d)?p2id:p0id; + } + else + { + d = DIST_SQ(p,p2); + i = (d < d1)? p2id:p1id; + } + return(i); +} + int -spherical_polygon_edge_intersect(a0,a1,b0,b1) +sedge_intersect(a0,a1,b0,b1) FVECT a0,a1,b0,b1; { FVECT na,nb,avga,avgb,p; @@ -879,4 +905,190 @@ FVECT a0,a1,b0,b1; 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) +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; + +} + + + + +bary_parent(coord,i) +BCOORD coord[3]; +int i; +{ + switch(i) { + case 0: + /* update bary for child */ + coord[0] = (coord[0] >> 1) + MAXBCOORD4; + coord[1] >>= 1; + coord[2] >>= 1; + break; + case 1: + coord[0] >>= 1; + coord[1] = (coord[1] >> 1) + MAXBCOORD4; + coord[2] >>= 1; + break; + + case 2: + coord[0] >>= 1; + coord[1] >>= 1; + coord[2] = (coord[2] >> 1) + MAXBCOORD4; + break; + + case 3: + coord[0] = MAXBCOORD4 - (coord[0] >> 1); + coord[1] = MAXBCOORD4 - (coord[1] >> 1); + coord[2] = MAXBCOORD4 - (coord[2] >> 1); + break; +#ifdef DEBUG + default: + eputs("bary_parent():Invalid child\n"); + break; +#endif + } +} + +bary_from_child(coord,child,next) +BCOORD coord[3]; +int child,next; +{ +#ifdef DEBUG + if(child <0 || child > 3) + { + eputs("bary_from_child():Invalid child\n"); + return; + } + if(next <0 || next > 3) + { + eputs("bary_from_child():Invalid next\n"); + return; + } +#endif + if(next == child) + return; + + switch(child){ + case 0: + coord[0] = 0; + coord[1] = MAXBCOORD2 - coord[1]; + coord[2] = MAXBCOORD2 - coord[2]; + break; + case 1: + coord[0] = MAXBCOORD2 - coord[0]; + coord[1] = 0; + coord[2] = MAXBCOORD2 - coord[2]; + break; + case 2: + coord[0] = MAXBCOORD2 - coord[0]; + coord[1] = MAXBCOORD2 - coord[1]; + coord[2] = 0; + break; + case 3: + switch(next){ + case 0: + coord[0] = 0; + coord[1] = MAXBCOORD2 - coord[1]; + coord[2] = MAXBCOORD2 - coord[2]; + break; + case 1: + coord[0] = MAXBCOORD2 - coord[0]; + coord[1] = 0; + coord[2] = MAXBCOORD2 - coord[2]; + break; + case 2: + coord[0] = MAXBCOORD2 - coord[0]; + coord[1] = MAXBCOORD2 - coord[1]; + coord[2] = 0; + break; + } + break; + } +} + +int +bary_child(coord) +BCOORD coord[3]; +{ + int i; + + if(coord[0] > MAXBCOORD4) + { + /* update bary for child */ + coord[0] = (coord[0]<< 1) - MAXBCOORD2; + coord[1] <<= 1; + coord[2] <<= 1; + return(0); + } + else + if(coord[1] > MAXBCOORD4) + { + coord[0] <<= 1; + coord[1] = (coord[1] << 1) - MAXBCOORD2; + coord[2] <<= 1; + return(1); + } + else + if(coord[2] > MAXBCOORD4) + { + coord[0] <<= 1; + coord[1] <<= 1; + coord[2] = (coord[2] << 1) - MAXBCOORD2; + return(2); + } + else + { + coord[0] = MAXBCOORD2 - (coord[0] << 1); + coord[1] = MAXBCOORD2 - (coord[1] << 1); + coord[2] = MAXBCOORD2 - (coord[2] << 1); + return(3); + } +} + +int +bary_nth_child(coord,i) +BCOORD coord[3]; +int i; +{ + + switch(i){ + case 0: + /* update bary for child */ + coord[0] = (coord[0]<< 1) - MAXBCOORD2; + coord[1] <<= 1; + coord[2] <<= 1; + break; + case 1: + coord[0] <<= 1; + coord[1] = (coord[1] << 1) - MAXBCOORD2; + coord[2] <<= 1; + break; + case 2: + coord[0] <<= 1; + coord[1] <<= 1; + coord[2] = (coord[2] << 1) - MAXBCOORD2; + break; + case 3: + coord[0] = MAXBCOORD2 - (coord[0] << 1); + coord[1] = MAXBCOORD2 - (coord[1] << 1); + coord[2] = MAXBCOORD2 - (coord[2] << 1); + break; + } +} + +