--- ray/src/hd/sm_geom.c 1998/08/25 11:03:28 3.3 +++ ray/src/hd/sm_geom.c 1998/10/06 18:16:53 3.7 @@ -33,13 +33,20 @@ int convex_angle(v0,v1,v2) FVECT v0,v1,v2; { - FVECT cp01,cp12,cp; - - /* test sign of (v0Xv1)X(v1Xv2). v1 */ + FVECT cp,cp01,cp12,v10,v02; + double dp; + /* + VSUB(v10,v1,v0); + VSUB(v02,v0,v2); + VCROSS(cp,v10,v02); + */ + /* test sign of (v0Xv1)X(v1Xv2). v1 */ VCROSS(cp01,v0,v1); VCROSS(cp12,v1,v2); VCROSS(cp,cp01,cp12); - if(DOT(cp,v1) < 0) + + dp = DOT(cp,v1); + if(ZERO(dp) || dp < 0.0) return(FALSE); return(TRUE); } @@ -53,7 +60,7 @@ FVECT v0,v1,v2; double tri_normal(v0,v1,v2,n,norm) FVECT v0,v1,v2,n; -char norm; +int norm; { double mag; @@ -64,7 +71,6 @@ 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]); - n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) + (v1[1] + v2[1]) * (v1[0] - v2[0]) + @@ -80,32 +86,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) @@ -135,14 +125,19 @@ FVECT ps,p,c; } +/* 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,tptr,r) - FVECT v,plane_n; - double plane_d; +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,24 +147,35 @@ intersect_vector_plane(v,plane_n,plane_d,tptr,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 - hit = 0; - r[0] = v[0]*t; - r[1] = v[1]*t; - r[2] = v[2]*t; + 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,plane_n,plane_d,pd,r) +intersect_ray_plane(orig,dir,peq,pd,r) FVECT orig,dir; - FVECT plane_n; - double plane_d; + FPEQ peq; double *pd; FVECT r; { @@ -184,13 +190,81 @@ intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r) 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) + t = -(DOT(FP_N(peq),orig) + FP_D(peq))/(DOT(FP_N(peq),dir)); + if(t < 0) + hit = 0; + else hit = 1; + + if(r) + VSUM(r,orig,dir,t); + + if(pd) + *pd = t; + return(hit); +} + + +int +intersect_ray_oplane(orig,dir,n,pd,r) + FVECT orig,dir; + FVECT n; + 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(n,orig))/(DOT(n,dir)); + if(t < 0) + hit = 0; else + hit = 1; + + if(r) + VSUM(r,orig,dir,t); + + if(pd) + *pd = t; + return(hit); +} + + +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; @@ -203,9 +277,9 @@ 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 @@ -213,21 +287,22 @@ FVECT p0,p1,p2; */ /* find the equation of the plane defined by p1-p3 */ - tri_plane_equation(p0,p1,p2,n,&d,FALSE); + 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*/ 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); @@ -252,16 +327,14 @@ 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)) { @@ -271,14 +344,11 @@ char sides[3]; /* 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; @@ -291,19 +361,8 @@ char sides[3]; } /* 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); @@ -313,163 +372,67 @@ char sides[3]; /* 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)) + 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); /* 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); /* 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); /* 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 @@ -477,66 +440,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); } @@ -544,12 +492,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; @@ -561,7 +509,7 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, VCROSS(n[0],t1,t0); /* Test the point for sidedness */ d = DOT(n[0],p0); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[0],0); } else @@ -574,7 +522,7 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, VCROSS(n[1],t2,t1); /* Test the point for sidedness */ d = DOT(n[1],p0); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[0],1); } else @@ -587,7 +535,7 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, VCROSS(n[2],t0,t2); /* Test the point for sidedness */ d = DOT(n[2],p0); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[0],2); } else @@ -603,7 +551,7 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, VCROSS(n[0],t1,t0); /* Test the point for sidedness */ d = DOT(n[0],p1); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[1],0); } else @@ -617,7 +565,7 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, VCROSS(n[1],t2,t1); /* Test the point for sidedness */ d = DOT(n[1],p1); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[1],1); } else @@ -631,7 +579,7 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, VCROSS(n[2],t0,t2); /* Test the point for sidedness */ d = DOT(n[2],p1); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[1],2); } else @@ -647,7 +595,7 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, VCROSS(n[0],t1,t0); /* Test the point for sidedness */ d = DOT(n[0],p2); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[2],0); } else @@ -660,7 +608,7 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, VCROSS(n[1],t2,t1); /* Test the point for sidedness */ d = DOT(n[1],p2); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[2],1); } else @@ -673,7 +621,7 @@ set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset, VCROSS(n[2],t0,t2); /* Test the point for sidedness */ d = DOT(n[2],p2); - if(d >= 0) + if(d >= 0.0) SET_NTH_BIT(test[2],2); } else @@ -683,12 +631,12 @@ 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]; + 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]; /* Test the vertices of triangle a against the pyramid formed by triangle @@ -696,12 +644,10 @@ FVECT a1,a2,a3,b1,b2,b3; 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,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])) + if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0])) return(TRUE); @@ -749,31 +695,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); } @@ -832,11 +773,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; @@ -896,45 +875,52 @@ 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; } +bary_ith_child(coord,i) +double coord[3]; + int i; +{ + + switch(i){ + case 0: + /* update bary for child */ + coord[0] = 2.0*coord[0]- 1.0; + coord[1] *= 2.0; + coord[2] *= 2.0; + break; + case 1: + coord[0] *= 2.0; + coord[1] = 2.0*coord[1]- 1.0; + coord[2] *= 2.0; + break; + case 2: + coord[0] *= 2.0; + coord[1] *= 2.0; + coord[2] = 2.0*coord[2]- 1.0; + break; + case 3: + coord[0] = 1.0 - 2.0*coord[0]; + coord[1] = 1.0 - 2.0*coord[1]; + coord[2] = 1.0 - 2.0*coord[2]; + break; +#ifdef DEBUG + default: + eputs("bary_ith_child():Invalid child\n"); + break; +#endif + } +} + + int -bary2d_child(coord) +bary_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 */ @@ -968,96 +954,501 @@ double coord[3]; } } -int -max_index(v) -FVECT v; +/* Coord was the ith child of the parent: set the coordinate + relative to the parent +*/ +bary_parent(coord,i) +double coord[3]; +int i; { - 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); + switch(i) { + case 0: + /* update bary for child */ + coord[0] = coord[0]*0.5 + 0.5; + coord[1] *= 0.5; + coord[2] *= 0.5; + break; + case 1: + coord[0] *= 0.5; + coord[1] = coord[1]*0.5 + 0.5; + coord[2] *= 0.5; + break; + + case 2: + coord[0] *= 0.5; + coord[1] *= 0.5; + coord[2] = coord[2]*0.5 + 0.5; + break; + + case 3: + coord[0] = 0.5 - 0.5*coord[0]; + coord[1] = 0.5 - 0.5*coord[1]; + coord[2] = 0.5 - 0.5*coord[2]; + break; +#ifdef DEBUG + default: + eputs("bary_parent():Invalid child\n"); + break; +#endif + } } +bary_from_child(coord,child,next) +double 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: + switch(next){ + case 1: + coord[0] += 1.0; + coord[1] -= 1.0; + break; + case 2: + coord[0] += 1.0; + coord[2] -= 1.0; + break; + case 3: + coord[0] *= -1.0; + coord[1] = 1 - coord[1]; + coord[2] = 1 - coord[2]; + break; -/* - * int - * smRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r) - * - * Intersect the ray with triangle v0v1v2, return intersection point in r - * - * Assumes orig,v0,v1,v2 are in spherical coordinates, and orig is - * unit - */ + } + break; + case 1: + switch(next){ + case 0: + coord[0] -= 1.0; + coord[1] += 1.0; + break; + case 2: + coord[1] += 1.0; + coord[2] -= 1.0; + break; + case 3: + coord[0] = 1 - coord[0]; + coord[1] *= -1.0; + coord[2] = 1 - coord[2]; + break; + } + break; + case 2: + switch(next){ + case 0: + coord[0] -= 1.0; + coord[2] += 1.0; + break; + case 1: + coord[1] -= 1.0; + coord[2] += 1.0; + break; + case 3: + coord[0] = 1 - coord[0]; + coord[1] = 1 - coord[1]; + coord[2] *= -1.0; + break; + } + break; + case 3: + switch(next){ + case 0: + coord[0] *= -1.0; + coord[1] = 1 - coord[1]; + coord[2] = 1 - coord[2]; + break; + case 1: + coord[0] = 1 - coord[0]; + coord[1] *= -1.0; + coord[2] = 1 - coord[2]; + break; + case 2: + coord[0] = 1 - coord[0]; + coord[1] = 1 - coord[1]; + coord[2] *= -1.0; + break; + } + break; + } +} + + +baryi_parent(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; +#ifdef DEBUG + default: + eputs("baryi_parent():Invalid child\n"); + break; +#endif + } +} + +baryi_from_child(coord,child,next) +BCOORD coord[3]; +int child,next; +{ +#ifdef DEBUG + if(child <0 || child > 3) + { + eputs("baryi_from_child():Invalid child\n"); + return; + } + if(next <0 || next > 3) + { + eputs("baryi_from_child():Invalid next\n"); + return; + } +#endif + if(next == child) + return; + + switch(child){ + case 0: + coord[0] = 0; + coord[1] = MAXBCOORD - coord[1]; + coord[2] = MAXBCOORD - coord[2]; + break; + case 1: + coord[0] = MAXBCOORD - coord[0]; + coord[1] = 0; + coord[2] = MAXBCOORD - coord[2]; + break; + case 2: + coord[0] = MAXBCOORD - coord[0]; + coord[1] = MAXBCOORD - coord[1]; + coord[2] = 0; + break; + case 3: + switch(next){ + case 0: + coord[0] = 0; + coord[1] = MAXBCOORD - coord[1]; + coord[2] = MAXBCOORD - coord[2]; + break; + case 1: + coord[0] = MAXBCOORD - coord[0]; + coord[1] = 0; + coord[2] = MAXBCOORD - coord[2]; + break; + case 2: + coord[0] = MAXBCOORD - coord[0]; + coord[1] = MAXBCOORD - coord[1]; + coord[2] = 0; + break; + } + break; + } +} + int -traceRay(orig,dir,v0,v1,v2,r) - FVECT orig,dir; - FVECT v0,v1,v2; - FVECT r; +baryi_child(coord) +BCOORD coord[3]; { - FVECT n,p[3],d; - double pt[3],r_eps; - char i; - int which; + int i; - /* Find the plane equation for the triangle defined by the edge v0v1 and - the view center*/ - VCROSS(n,v0,v1); - /* Intersect the ray with this plane */ - i = intersect_ray_plane(orig,dir,n,0.0,&(pt[0]),p[0]); - if(i) - which = 0; + if(coord[0] > MAXBCOORD2) + { + /* update bary for child */ + coord[0] = (coord[0]<< 1) - MAXBCOORD; + coord[1] <<= 1; + coord[2] <<= 1; + return(0); + } else - which = -1; + if(coord[1] > MAXBCOORD2) + { + coord[0] <<= 1; + coord[1] = (coord[1] << 1) - MAXBCOORD; + coord[2] <<= 1; + return(1); + } + else + if(coord[2] > MAXBCOORD2) + { + coord[0] <<= 1; + coord[1] <<= 1; + coord[2] = (coord[2] << 1) - MAXBCOORD; + return(2); + } + else + { + coord[0] = MAXBCOORD - (coord[0] << 1); + coord[1] = MAXBCOORD - (coord[1] << 1); + coord[2] = MAXBCOORD - (coord[2] << 1); + return(3); + } +} - VCROSS(n,v1,v2); - i = intersect_ray_plane(orig,dir,n,0.0,&(pt[1]),p[1]); - if(i) - if((which==-1) || pt[1] < pt[0]) - which = 1; +int +baryi_nth_child(coord,i) +BCOORD coord[3]; +int i; +{ - VCROSS(n,v2,v0); - i = intersect_ray_plane(orig,dir,n,0.0,&(pt[2]),p[2]); - if(i) - if((which==-1) || pt[2] < pt[which]) - which = 2; + switch(i){ + case 0: + /* update bary for child */ + coord[0] = (coord[0]<< 1) - MAXBCOORD; + coord[1] <<= 1; + coord[2] <<= 1; + break; + case 1: + coord[0] <<= 1; + coord[1] = (coord[1] << 1) - MAXBCOORD; + coord[2] <<= 1; + break; + case 2: + coord[0] <<= 1; + coord[1] <<= 1; + coord[2] = (coord[2] << 1) - MAXBCOORD; + break; + case 3: + coord[0] = MAXBCOORD - (coord[0] << 1); + coord[1] = MAXBCOORD - (coord[1] << 1); + coord[2] = MAXBCOORD - (coord[2] << 1); + break; + } +} - if(which != -1) + +baryi_children(coord,i,in_tri,rcoord,rin_tri) +BCOORD coord[3][3]; +int i; +int in_tri[3]; +BCOORD rcoord[3][3]; +int rin_tri[3]; +{ + int j; + + for(j=0; j< 3; j++) { - /* Return point that is K*eps along projection of the ray on - the sphere to push intersection point p[which] into next cell - */ - normalize(p[which]); - /* Calculate the ray perpendicular to the intersection point: approx - the direction moved along the sphere a distance of K*epsilon*/ - r_eps = -DOT(p[which],dir); - VSUM(n,dir,p[which],r_eps); - /* Calculate the length along ray p[which]-dir needed to move to - cause a move along the sphere of k*epsilon - */ - r_eps = DOT(n,dir); - VSUM(r,p[which],dir,(20*FTINY)/r_eps); - normalize(r); - return(TRUE); + if(!in_tri[j]) + { + rin_tri[j]=0; + continue; + } + + if(i != 3) + { + if(coord[j][i] < MAXBCOORD2) + { + rin_tri[j] = 0; + continue; + } + } + else + if( !(coord[j][0] <= MAXBCOORD2 && coord[j][1] <= MAXBCOORD2 && + coord[j][2] <= MAXBCOORD2)) + { + rin_tri[j] = 0; + continue; + } + + rin_tri[j]=1; + VCOPY(rcoord[j],coord[j]); + baryi_nth_child(rcoord[j],i); } + +} + +convert_dtol(b,bi) +double b[3]; +BCOORD bi[3]; +{ + int i; + + for(i=0; i < 2;i++) + { + + if(b[i] <= 0.0) + { +#ifdef EXTRA_DEBUG + if(b[i] < 0.0) + printf("under %f\n",b[i]); +#endif + bi[i]= 0; + } + else + if(b[i] >= 1.0) + { +#ifdef EXTRA_DEBUG + if(b[i] > 1.0) + printf("over %f\n",b[i]); +#endif + bi[i]= MAXBCOORD; + } + else + bi[i] = (BCOORD)(b[i]*MAXBCOORD); + } + bi[2] = bi[0] + bi[1]; + if(bi[2] > MAXBCOORD) + { +#ifdef EXTRA_DEBUG + printf("sum over %f\n",b[0]+b[1]); +#endif + bi[2] = 0; + bi[1] = MAXBCOORD - bi[0]; + } else + bi[2] = MAXBCOORD - bi[2]; + +} + +/* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG], + dir unbounded to [-MAXLONG,MAXLONG] + */ +bary_dtol(b,db,bi,dbi,t,w) +double b[3],db[3][3]; +BCOORD bi[3]; +BDIR dbi[3][3]; +TINT t[3]; +int w; +{ + int i,id,j,k; + double d; + + convert_dtol(b,bi); + + for(j=w; j< 3; j++) { - /* Unable to find intersection: move along ray and try again */ - r_eps = -DOT(orig,dir); - VSUM(n,dir,orig,r_eps); - r_eps = DOT(n,dir); - VSUM(r,orig,dir,(20*FTINY)/r_eps); - normalize(r); -#ifdef DEBUG - eputs("traceRay:Ray does not intersect triangle"); -#endif - return(FALSE); + if(t[j] == HUGET) + { + max_index(db[j],&d); + for(i=0; i< 3; i++) + dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR); + break; + } + else + { + for(i=0; i< 3; i++) + dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR); + } } } + + +/* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG], + dir unbounded to [-MAXLONG,MAXLONG] + */ +bary_dtol_new(b,db,bi,boi,dbi,t) +double b[4][3],db[3][3]; +BCOORD bi[3],boi[3][3]; +BDIR dbi[3][3]; +int t[3]; +{ + int i,id,j,k; + double d; + + convert_dtol(b[3],bi); + + for(j=0; j<3;j++) + { + if(t[j] != 1) + continue; + convert_dtol(b[j],boi[j]); + } + for(j=0; j< 3; j++) + { + k = (j+1)%3; + if(t[k]==0) + continue; + if(t[k] == -1) + { + max_index(db[j],&d); + for(i=0; i< 3; i++) + dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR); + t[k] = 0; + } + else + if(t[j] != 1) + for(i=0; i< 3; i++) + dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR); + else + for(i=0; i< 3; i++) + dbi[j][i] = boi[k][i] - boi[j][i]; + + } +} + + +bary_dtolb(b,bi,in_tri) +double b[3][3]; +BCOORD bi[3][3]; +int in_tri[3]; +{ + int i,j; + + for(j=0; j<3;j++) + { + if(!in_tri[j]) + continue; + for(i=0; i < 2;i++) + { + if(b[j][i] <= 0.0) + { + bi[j][i]= 0; + } + else + if(b[j][i] >= 1.0) + { + bi[j][i]= MAXBCOORD; + } + else + bi[j][i] = (BCOORD)(b[j][i]*MAXBCOORD); + } + bi[j][2] = MAXBCOORD - bi[j][0] - bi[j][1]; + if(bi[j][2] < 0) + { + bi[j][2] = 0; + bi[j][1] = MAXBCOORD - bi[j][0]; + } + } +} + +