ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/hd/sm_geom.c
(Generate patch)

Comparing ray/src/hd/sm_geom.c (file contents):
Revision 3.6 by gwlarson, Wed Sep 16 18:16:28 1998 UTC vs.
Revision 3.11 by gwlarson, Sun Jan 10 10:27:47 1999 UTC

# Line 27 | Line 27 | FVECT v1,v2;
27   {
28     return(EQUAL(v1[0],v2[0]) && EQUAL(v1[1],v2[1])&& EQUAL(v1[2],v2[2]));
29   }
30 + #if 0
31 + extern FVECT Norm[500];
32 + extern int Ncnt;
33 + #endif
34  
31 –
35   int
36   convex_angle(v0,v1,v2)
37   FVECT v0,v1,v2;
38   {
39 <    FVECT cp01,cp12,cp;
40 <    
39 >    FVECT cp,cp01,cp12,v10,v02;
40 >    double dp;
41 >
42      /* test sign of (v0Xv1)X(v1Xv2). v1 */
43      VCROSS(cp01,v0,v1);
44      VCROSS(cp12,v1,v2);
45      VCROSS(cp,cp01,cp12);
46 <    if(DOT(cp,v1) < 0.0)
47 <       return(FALSE);
46 >        
47 >    dp = DOT(cp,v1);
48 > #if 0
49 >    VCOPY(Norm[Ncnt++],cp01);
50 >    VCOPY(Norm[Ncnt++],cp12);
51 >    VCOPY(Norm[Ncnt++],cp);
52 > #endif
53 >    if(ZERO(dp) || dp < 0.0)
54 >      return(FALSE);
55      return(TRUE);
56   }
57  
58   /* calculates the normal of a face contour using Newell's formula. e
59  
60 <               a =  SUMi (yi - yi+1)(zi + zi+1)
60 >               a =  SUMi (yi - yi+1)(zi + zi+1);
61                 b =  SUMi (zi - zi+1)(xi + xi+1)
62                 c =  SUMi (xi - xi+1)(yi + yi+1)
63   */
# Line 63 | Line 74 | int norm;
74    
75    n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
76             (v1[2] - v2[2]) * (v1[0] + v2[0]) +
77 <           (v2[2] - v0[2]) * (v2[0] + v0[0]);
67 <
77 >            (v2[2] - v0[2]) * (v2[0] + v0[0]);
78    
79    n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
80           (v1[1] + v2[1]) * (v1[0] - v2[0]) +
# Line 72 | Line 82 | int norm;
82  
83    if(!norm)
84       return(0);
75 –
85    
86    mag = normalize(n);
87  
# Line 80 | Line 89 | int norm;
89   }
90  
91  
92 < tri_plane_equation(v0,v1,v2,n,nd,norm)
93 <   FVECT v0,v1,v2,n;
94 <   double *nd;
92 > tri_plane_equation(v0,v1,v2,peqptr,norm)
93 >   FVECT v0,v1,v2;
94 >   FPEQ *peqptr;
95     int norm;
96   {
97 <    tri_normal(v0,v1,v2,n,norm);
97 >    tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
98  
99 <    *nd = -(DOT(n,v0));
99 >    FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
100   }
101  
102   /* From quad_edge-code */
# Line 106 | Line 115 | FVECT p0,p1;
115  
116      det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
117      
118 <    return (det > 0);
118 >    return (det > (1e-10));
119   }
120  
121  
122 <
122 > double
123   point_on_sphere(ps,p,c)
124   FVECT ps,p,c;
125   {
126 +  double d;
127      VSUB(ps,p,c);    
128 <    normalize(ps);
128 >    d= normalize(ps);
129 >    return(d);
130   }
131  
132  
# Line 125 | Line 136 | FVECT ps,p,c;
136    * t, if parallel, returns t=FHUGE
137    */
138   int
139 < intersect_vector_plane(v,plane_n,plane_d,tptr,r)
140 <   FVECT v,plane_n;
141 <   double plane_d;
139 > intersect_vector_plane(v,peq,tptr,r)
140 >   FVECT v;
141 >   FPEQ peq;
142     double *tptr;
143     FVECT r;
144   {
# Line 141 | Line 152 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
152      /* line is  l = p1 + (p2-p1)t, p1=origin */
153  
154      /* Solve for t: */
155 <  d = -(DOT(plane_n,v));
155 >  d = -(DOT(FP_N(peq),v));
156    if(ZERO(d))
157    {
158        t = FHUGE;
# Line 149 | Line 160 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
160    }
161    else
162    {
163 <      t =  plane_d/d;
163 >      t =  FP_D(peq)/d;
164        if(t < 0 )
165           hit = 0;
166        else
# Line 167 | Line 178 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
178   }
179  
180   int
181 < intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
181 > intersect_ray_plane(orig,dir,peq,pd,r)
182     FVECT orig,dir;
183 <   FVECT plane_n;
173 <   double plane_d;
183 >   FPEQ peq;
184     double *pd;
185     FVECT r;
186   {
187 <  double t;
187 >  double t,d;
188    int hit;
189      /*
190        Plane is Ax + By + Cz +D = 0:
# Line 185 | Line 195 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
195         line is  l = p1 + (p2-p1)t
196       */
197      /* Solve for t: */
198 <    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
198 >  d = DOT(FP_N(peq),dir);
199 >  if(ZERO(d))
200 >     return(0);
201 >  t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
202 >
203 >  if(t < 0)
204 >       hit = 0;
205 >    else
206 >       hit = 1;
207 >
208 >  if(r)
209 >     VSUM(r,orig,dir,t);
210 >
211 >    if(pd)
212 >       *pd = t;
213 >  return(hit);
214 > }
215 >
216 >
217 > int
218 > intersect_ray_oplane(orig,dir,n,pd,r)
219 >   FVECT orig,dir;
220 >   FVECT n;
221 >   double *pd;
222 >   FVECT r;
223 > {
224 >  double t,d;
225 >  int hit;
226 >    /*
227 >      Plane is Ax + By + Cz +D = 0:
228 >      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
229 >    */
230 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
231 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
232 >       line is  l = p1 + (p2-p1)t
233 >     */
234 >    /* Solve for t: */
235 >    d= DOT(n,dir);
236 >    if(ZERO(d))
237 >       return(0);
238 >    t =  -(DOT(n,orig))/d;
239      if(t < 0)
240         hit = 0;
241      else
# Line 199 | Line 249 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
249    return(hit);
250   }
251  
252 + /* Assumption: know crosses plane:dont need to check for 'on' case */
253 + intersect_edge_coord_plane(v0,v1,w,r)
254 + FVECT v0,v1;
255 + int w;
256 + FVECT r;
257 + {
258 +  FVECT dv;
259 +  int wnext;
260 +  double t;
261  
262 +  VSUB(dv,v1,v0);
263 +  t = -v0[w]/dv[w];
264 +  r[w] = 0.0;
265 +  wnext = (w+1)%3;
266 +  r[wnext] = v0[wnext] + dv[wnext]*t;
267 +  wnext = (w+2)%3;
268 +  r[wnext] = v0[wnext] + dv[wnext]*t;
269 + }
270 +
271   int
272 < intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
272 > intersect_edge_plane(e0,e1,peq,pd,r)
273     FVECT e0,e1;
274 <   FVECT plane_n;
207 <   double plane_d;
274 >   FPEQ peq;
275     double *pd;
276     FVECT r;
277   {
# Line 221 | Line 288 | intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
288       */
289      /* Solve for t: */
290    VSUB(d,e1,e0);
291 <  t =  -(DOT(plane_n,e0) + plane_d)/(DOT(plane_n,d));
291 >  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
292      if(t < 0)
293         hit = 0;
294      else
# Line 240 | Line 307 | point_in_cone(p,p0,p1,p2)
307   FVECT p;
308   FVECT p0,p1,p2;
309   {
243 –    FVECT n;
310      FVECT np,x_axis,y_axis;
311 <    double d1,d2,d;
311 >    double d1,d2;
312 >    FPEQ peq;
313      
314      /* Find the equation of the circle defined by the intersection
315         of the cone with the plane defined by p1,p2,p3- project p into
316         that plane and do an in-circle test in the plane
317       */
318      
319 <    /* find the equation of the plane defined by p1-p3 */
320 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
319 >    /* find the equation of the plane defined by p0-p2 */
320 >    tri_plane_equation(p0,p1,p2,&peq,FALSE);
321  
322      /* define a coordinate system on the plane: the x axis is in
323         the direction of np2-np1, and the y axis is calculated from
# Line 258 | Line 325 | FVECT p0,p1,p2;
325       */
326      /* Project p onto the plane */
327      /* NOTE: check this: does sideness check?*/
328 <    if(!intersect_vector_plane(p,n,d,NULL,np))
328 >    if(!intersect_vector_plane(p,peq,NULL,np))
329          return(FALSE);
330  
331 <    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
331 >    /* create coordinate system on  plane: p1-p0 defines the x_axis*/
332      VSUB(x_axis,p1,p0);
333      normalize(x_axis);
334      /* The y axis is  */
335 <    VCROSS(y_axis,n,x_axis);
335 >    VCROSS(y_axis,FP_N(peq),x_axis);
336      normalize(y_axis);
337  
338      VSUB(p1,p1,p0);
339      VSUB(p2,p2,p0);
340      VSUB(np,np,p0);
341      
342 <    p1[0] = VLEN(p1);
342 >    p1[0] = DOT(p1,x_axis);
343      p1[1] = 0;
344  
345      d1 = DOT(p2,x_axis);
# Line 301 | Line 368 | int sides[3];
368      /* Find the normal to the triangle ORIGIN:v0:v1 */
369      if(!NTH_BIT(*nset,0))
370      {
371 <        VCROSS(n[0],v1,v0);
371 >        VCROSS(n[0],v0,v1);
372          SET_NTH_BIT(*nset,0);
373      }
374      /* Test the point for sidedness */
# Line 319 | Line 386 | int sides[3];
386      /* Test next edge */
387      if(!NTH_BIT(*nset,1))
388      {
389 <        VCROSS(n[1],v2,v1);
389 >        VCROSS(n[1],v1,v2);
390          SET_NTH_BIT(*nset,1);
391      }
392      /* Test the point for sidedness */
# Line 335 | Line 402 | int sides[3];
402      /* Test next edge */
403      if(!NTH_BIT(*nset,2))
404      {
405 <        VCROSS(n[2],v0,v2);
405 >        VCROSS(n[2],v2,v0);
406          SET_NTH_BIT(*nset,2);
407      }
408      /* Test the point for sidedness */
# Line 353 | Line 420 | int sides[3];
420  
421  
422  
423 <
423 >
424   int
425   point_in_stri(v0,v1,v2,p)
426   FVECT v0,v1,v2,p;
# Line 361 | Line 428 | FVECT v0,v1,v2,p;
428      double d;
429      FVECT n;  
430  
431 <    VCROSS(n,v1,v0);
431 >    VCROSS(n,v0,v1);
432      /* Test the point for sidedness */
433      d  = DOT(n,p);
434      if(d > 0.0)
435        return(FALSE);
436  
437      /* Test next edge */
438 <    VCROSS(n,v2,v1);
438 >    VCROSS(n,v1,v2);
439      /* Test the point for sidedness */
440      d  = DOT(n,p);
441      if(d > 0.0)
442         return(FALSE);
443  
444      /* Test next edge */
445 <    VCROSS(n,v0,v2);
445 >    VCROSS(n,v2,v0);
446      /* Test the point for sidedness */
447      d  = DOT(n,p);
448      if(d > 0.0)
# Line 469 | Line 536 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
536      if(sides[0][0] == GT_INVALID)
537      {
538        if(!NTH_BIT(nset,0))
539 <        VCROSS(n[0],t1,t0);
539 >        VCROSS(n[0],t0,t1);
540        /* Test the point for sidedness */
541        d  = DOT(n[0],p0);
542        if(d >= 0.0)
# Line 482 | Line 549 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
549      if(sides[0][1] == GT_INVALID)
550      {
551        if(!NTH_BIT(nset,1))
552 <        VCROSS(n[1],t2,t1);
552 >        VCROSS(n[1],t1,t2);
553          /* Test the point for sidedness */
554          d  = DOT(n[1],p0);
555          if(d >= 0.0)
# Line 495 | Line 562 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
562      if(sides[0][2] == GT_INVALID)
563      {
564        if(!NTH_BIT(nset,2))
565 <        VCROSS(n[2],t0,t2);
565 >        VCROSS(n[2],t2,t0);
566        /* Test the point for sidedness */
567        d  = DOT(n[2],p0);
568        if(d >= 0.0)
# Line 511 | Line 578 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
578      if(sides[1][0] == GT_INVALID)
579      {
580        if(!NTH_BIT(nset,0))
581 <        VCROSS(n[0],t1,t0);
581 >        VCROSS(n[0],t0,t1);
582        /* Test the point for sidedness */
583        d  = DOT(n[0],p1);
584        if(d >= 0.0)
# Line 525 | Line 592 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
592      if(sides[1][1] == GT_INVALID)
593      {
594        if(!NTH_BIT(nset,1))
595 <        VCROSS(n[1],t2,t1);
595 >        VCROSS(n[1],t1,t2);
596        /* Test the point for sidedness */
597        d  = DOT(n[1],p1);
598        if(d >= 0.0)
# Line 539 | Line 606 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
606      if(sides[1][2] == GT_INVALID)
607      {
608        if(!NTH_BIT(nset,2))
609 <        VCROSS(n[2],t0,t2);
609 >        VCROSS(n[2],t2,t0);
610        /* Test the point for sidedness */
611        d  = DOT(n[2],p1);
612        if(d >= 0.0)
# Line 555 | Line 622 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
622      if(sides[2][0] == GT_INVALID)
623      {
624        if(!NTH_BIT(nset,0))
625 <        VCROSS(n[0],t1,t0);
625 >        VCROSS(n[0],t0,t1);
626        /* Test the point for sidedness */
627        d  = DOT(n[0],p2);
628        if(d >= 0.0)
# Line 568 | Line 635 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
635      if(sides[2][1] == GT_INVALID)
636      {
637        if(!NTH_BIT(nset,1))
638 <        VCROSS(n[1],t2,t1);
638 >        VCROSS(n[1],t1,t2);
639        /* Test the point for sidedness */
640        d  = DOT(n[1],p2);
641        if(d >= 0.0)
# Line 581 | Line 648 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
648      if(sides[2][2] == GT_INVALID)
649      {
650        if(!NTH_BIT(nset,2))
651 <        VCROSS(n[2],t0,t2);
651 >        VCROSS(n[2],t2,t0);
652        /* Test the point for sidedness */
653        d  = DOT(n[2],p2);
654        if(d >= 0.0)
# Line 600 | Line 667 | FVECT a1,a2,a3,b1,b2,b3;
667    int which,test,n_set[2];
668    int sides[2][3][3],i,j,inext,jnext;
669    int tests[2][3];
670 <  FVECT n[2][3],p,avg[2];
670 >  FVECT n[2][3],p,avg[2],t1,t2,t3;
671  
672 + #ifdef DEBUG
673 +    tri_normal(b1,b2,b3,p,FALSE);
674 +    if(DOT(p,b1) > 0)
675 +      {
676 +       VCOPY(t3,b1);
677 +       VCOPY(t2,b2);
678 +       VCOPY(t1,b3);
679 +      }
680 +    else
681 + #endif
682 +      {
683 +       VCOPY(t1,b1);
684 +       VCOPY(t2,b2);
685 +       VCOPY(t3,b3);
686 +      }
687 +
688    /* Test the vertices of triangle a against the pyramid formed by triangle
689       b and the origin. If any vertex of a is interior to triangle b, or
690       if all 3 vertices of a are ON the edges of b,return TRUE. Remember
691       the results of the edge normal and sidedness tests for later.
692     */
693 < if(vertices_in_stri(a1,a2,a3,b1,b2,b3,&(n_set[0]),n[0],avg[0],sides[1]))
693 > if(vertices_in_stri(a1,a2,a3,t1,t2,t3,&(n_set[0]),n[0],avg[0],sides[1]))
694       return(TRUE);
695    
696 < if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
696 > if(vertices_in_stri(t1,t2,t3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
697       return(TRUE);
698  
699  
700 <  set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
700 >  set_sidedness_tests(t1,t2,t3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
701    if(tests[0][0]&tests[0][1]&tests[0][2])
702      return(FALSE);
703  
704 <  set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
704 >  set_sidedness_tests(a1,a2,a3,t1,t2,t3,tests[1],sides[1],n_set[0],n[0]);
705    if(tests[1][0]&tests[1][1]&tests[1][2])
706      return(FALSE);
707  
# Line 663 | Line 746 | FVECT orig,dir;
746   FVECT v0,v1,v2;
747   FVECT pt;
748   {
749 <  FVECT p0,p1,p2,p,n;
750 <  double pd;
749 >  FVECT p0,p1,p2,p;
750 >  FPEQ peq;
751    int type;
752  
753    VSUB(p0,v0,orig);
# Line 674 | Line 757 | FVECT pt;
757    if(point_in_stri(p0,p1,p2,dir))
758    {
759        /* Intersect the ray with the triangle plane */
760 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
761 <      return(intersect_ray_plane(orig,dir,n,pd,NULL,pt));
760 >      tri_plane_equation(v0,v1,v2,&peq,FALSE);
761 >      return(intersect_ray_plane(orig,dir,peq,NULL,pt));
762    }
763    return(FALSE);
764   }
# Line 824 | Line 907 | FVECT a0,a1,b0,b1;
907   }
908  
909  
827 –
910   /* Find the normalized barycentric coordinates of p relative to
911   * triangle v0,v1,v2. Return result in coord
912   */
# Line 842 | Line 924 | double coord[3];
924  
925   }
926  
845 – bary_ith_child(coord,i)
846 – double coord[3];
847 – int i;
848 – {
849 –
850 –  switch(i){
851 –  case 0:
852 –      /* update bary for child */
853 –      coord[0] = 2.0*coord[0]- 1.0;
854 –      coord[1] *= 2.0;
855 –      coord[2] *= 2.0;
856 –      break;
857 –  case 1:
858 –    coord[0] *= 2.0;
859 –    coord[1] = 2.0*coord[1]- 1.0;
860 –    coord[2] *= 2.0;
861 –    break;
862 –  case 2:
863 –    coord[0] *= 2.0;
864 –    coord[1] *= 2.0;
865 –    coord[2] = 2.0*coord[2]- 1.0;
866 –    break;
867 –  case 3:
868 –    coord[0] = 1.0 - 2.0*coord[0];
869 –    coord[1] = 1.0 - 2.0*coord[1];
870 –    coord[2] = 1.0 - 2.0*coord[2];
871 –    break;
872 – #ifdef DEBUG
873 –  default:
874 –    eputs("bary_ith_child():Invalid child\n");
875 –    break;
876 – #endif
877 –  }
878 – }
927  
928  
881 – int
882 – bary_child(coord)
883 – double coord[3];
884 – {
885 –  int i;
929  
887 –  if(coord[0] > 0.5)
888 –  {
889 –      /* update bary for child */
890 –      coord[0] = 2.0*coord[0]- 1.0;
891 –      coord[1] *= 2.0;
892 –      coord[2] *= 2.0;
893 –      return(0);
894 –  }
895 –  else
896 –    if(coord[1] > 0.5)
897 –    {
898 –      coord[0] *= 2.0;
899 –      coord[1] = 2.0*coord[1]- 1.0;
900 –      coord[2] *= 2.0;
901 –      return(1);
902 –    }
903 –    else
904 –      if(coord[2] > 0.5)
905 –      {
906 –        coord[0] *= 2.0;
907 –        coord[1] *= 2.0;
908 –        coord[2] = 2.0*coord[2]- 1.0;
909 –        return(2);
910 –      }
911 –      else
912 –         {
913 –           coord[0] = 1.0 - 2.0*coord[0];
914 –           coord[1] = 1.0 - 2.0*coord[1];
915 –           coord[2] = 1.0 - 2.0*coord[2];
916 –           return(3);
917 –         }
918 – }
919 –
920 – /* Coord was the ith child of the parent: set the coordinate
921 –   relative to the parent
922 – */
930   bary_parent(coord,i)
924 – double coord[3];
925 – int i;
926 – {
927 –
928 –  switch(i) {
929 –  case 0:
930 –    /* update bary for child */
931 –    coord[0] = coord[0]*0.5 + 0.5;
932 –    coord[1] *= 0.5;
933 –    coord[2] *= 0.5;
934 –    break;
935 –  case 1:
936 –    coord[0] *= 0.5;
937 –    coord[1]  = coord[1]*0.5 + 0.5;
938 –    coord[2] *= 0.5;
939 –    break;
940 –    
941 –  case 2:
942 –    coord[0] *= 0.5;
943 –    coord[1] *= 0.5;
944 –    coord[2] = coord[2]*0.5 + 0.5;
945 –    break;
946 –    
947 –  case 3:
948 –    coord[0] = 0.5 - 0.5*coord[0];
949 –    coord[1] = 0.5 - 0.5*coord[1];
950 –    coord[2] = 0.5 - 0.5*coord[2];
951 –    break;
952 – #ifdef DEBUG
953 –  default:
954 –    eputs("bary_parent():Invalid child\n");
955 –    break;
956 – #endif
957 –  }
958 – }
959 –
960 – bary_from_child(coord,child,next)
961 – double coord[3];
962 – int child,next;
963 – {
964 – #ifdef DEBUG
965 –  if(child <0 || child > 3)
966 –  {
967 –    eputs("bary_from_child():Invalid child\n");
968 –    return;
969 –  }
970 –  if(next <0 || next > 3)
971 –  {
972 –    eputs("bary_from_child():Invalid next\n");
973 –    return;
974 –  }
975 – #endif
976 –  if(next == child)
977 –    return;
978 –
979 –  switch(child){
980 –  case 0:
981 –    switch(next){
982 –    case 1:
983 –      coord[0] += 1.0;
984 –      coord[1] -= 1.0;
985 –      break;
986 –    case 2:
987 –      coord[0] += 1.0;
988 –      coord[2] -= 1.0;
989 –      break;
990 –    case 3:
991 –      coord[0] *= -1.0;
992 –      coord[1] = 1 - coord[1];
993 –      coord[2] = 1 - coord[2];
994 –      break;
995 –
996 –    }
997 –    break;
998 –  case 1:
999 –    switch(next){
1000 –    case 0:
1001 –      coord[0] -= 1.0;
1002 –      coord[1] += 1.0;
1003 –      break;
1004 –    case 2:
1005 –      coord[1] += 1.0;
1006 –      coord[2] -= 1.0;
1007 –      break;
1008 –    case 3:
1009 –      coord[0] = 1 - coord[0];
1010 –      coord[1] *= -1.0;
1011 –      coord[2] = 1 - coord[2];
1012 –      break;
1013 –    }
1014 –    break;
1015 –  case 2:
1016 –    switch(next){
1017 –    case 0:
1018 –      coord[0] -= 1.0;
1019 –      coord[2] += 1.0;
1020 –      break;
1021 –    case 1:
1022 –      coord[1] -= 1.0;
1023 –      coord[2] += 1.0;
1024 –      break;
1025 –    case 3:
1026 –      coord[0] = 1 - coord[0];
1027 –      coord[1] = 1 - coord[1];
1028 –      coord[2] *= -1.0;
1029 –      break;
1030 –    }
1031 –    break;
1032 –  case 3:
1033 –    switch(next){
1034 –    case 0:
1035 –      coord[0] *= -1.0;
1036 –      coord[1] = 1 - coord[1];
1037 –      coord[2] = 1 - coord[2];
1038 –      break;
1039 –    case 1:
1040 –      coord[0] = 1 - coord[0];
1041 –      coord[1] *= -1.0;
1042 –      coord[2] = 1 - coord[2];
1043 –      break;
1044 –    case 2:
1045 –      coord[0] = 1 - coord[0];
1046 –      coord[1] = 1 - coord[1];
1047 –      coord[2] *= -1.0;
1048 –      break;
1049 –    }
1050 –    break;
1051 –  }
1052 – }
1053 –
1054 –
1055 – baryi_parent(coord,i)
931   BCOORD coord[3];
932   int i;
933   {
1059 –
934    switch(i) {
935    case 0:
936      /* update bary for child */
937 <    coord[0] = (coord[0] >> 1) + MAXBCOORD2;
937 >    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
938      coord[1] >>= 1;
939      coord[2] >>= 1;
940      break;
941    case 1:
942      coord[0] >>= 1;
943 <    coord[1]  = (coord[1] >> 1) + MAXBCOORD2;
943 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
944      coord[2] >>= 1;
945      break;
946      
947    case 2:
948      coord[0] >>= 1;
949      coord[1] >>= 1;
950 <    coord[2] = (coord[2] >> 1) + MAXBCOORD2;
950 >    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
951      break;
952      
953    case 3:
954 <    coord[0] = MAXBCOORD2 - (coord[0] >> 1);
955 <    coord[1] = MAXBCOORD2 - (coord[1] >> 1);
956 <    coord[2] = MAXBCOORD2 - (coord[2] >> 1);
954 >    coord[0] = MAXBCOORD4 - (coord[0] >> 1);
955 >    coord[1] = MAXBCOORD4 - (coord[1] >> 1);
956 >    coord[2] = MAXBCOORD4 - (coord[2] >> 1);
957      break;
958   #ifdef DEBUG
959    default:
960 <    eputs("baryi_parent():Invalid child\n");
960 >    eputs("bary_parent():Invalid child\n");
961      break;
962   #endif
963    }
964   }
965  
966 < baryi_from_child(coord,child,next)
966 > bary_from_child(coord,child,next)
967   BCOORD coord[3];
968   int child,next;
969   {
970   #ifdef DEBUG
971    if(child <0 || child > 3)
972    {
973 <    eputs("baryi_from_child():Invalid child\n");
973 >    eputs("bary_from_child():Invalid child\n");
974      return;
975    }
976    if(next <0 || next > 3)
977    {
978 <    eputs("baryi_from_child():Invalid next\n");
978 >    eputs("bary_from_child():Invalid next\n");
979      return;
980    }
981   #endif
# Line 1111 | Line 985 | int child,next;
985    switch(child){
986    case 0:
987        coord[0] = 0;
988 <      coord[1] = MAXBCOORD - coord[1];
989 <      coord[2] = MAXBCOORD - coord[2];
988 >      coord[1] = MAXBCOORD2 - coord[1];
989 >      coord[2] = MAXBCOORD2 - coord[2];
990        break;
991    case 1:
992 <      coord[0] = MAXBCOORD - coord[0];
992 >      coord[0] = MAXBCOORD2 - coord[0];
993        coord[1] = 0;
994 <      coord[2] = MAXBCOORD - coord[2];
994 >      coord[2] = MAXBCOORD2 - coord[2];
995        break;
996    case 2:
997 <      coord[0] = MAXBCOORD - coord[0];
998 <      coord[1] = MAXBCOORD - coord[1];
997 >      coord[0] = MAXBCOORD2 - coord[0];
998 >      coord[1] = MAXBCOORD2 - coord[1];
999        coord[2] = 0;
1000      break;
1001    case 3:
1002      switch(next){
1003      case 0:
1004        coord[0] = 0;
1005 <      coord[1] = MAXBCOORD - coord[1];
1006 <      coord[2] = MAXBCOORD - coord[2];
1005 >      coord[1] = MAXBCOORD2 - coord[1];
1006 >      coord[2] = MAXBCOORD2 - coord[2];
1007        break;
1008      case 1:
1009 <      coord[0] = MAXBCOORD - coord[0];
1009 >      coord[0] = MAXBCOORD2 - coord[0];
1010        coord[1] = 0;
1011 <      coord[2] = MAXBCOORD - coord[2];
1011 >      coord[2] = MAXBCOORD2 - coord[2];
1012        break;
1013      case 2:
1014 <      coord[0] = MAXBCOORD - coord[0];
1015 <      coord[1] = MAXBCOORD - coord[1];
1014 >      coord[0] = MAXBCOORD2 - coord[0];
1015 >      coord[1] = MAXBCOORD2 - coord[1];
1016        coord[2] = 0;
1017        break;
1018      }
# Line 1147 | Line 1021 | int child,next;
1021   }
1022  
1023   int
1024 < baryi_child(coord)
1024 > bary_child(coord)
1025   BCOORD coord[3];
1026   {
1027    int i;
1028  
1029 <  if(coord[0] > MAXBCOORD2)
1029 >  if(coord[0] > MAXBCOORD4)
1030    {
1031        /* update bary for child */
1032 <      coord[0] = (coord[0]<< 1) - MAXBCOORD;
1032 >      coord[0] = (coord[0]<< 1) - MAXBCOORD2;
1033        coord[1] <<= 1;
1034        coord[2] <<= 1;
1035        return(0);
1036    }
1037    else
1038 <    if(coord[1] > MAXBCOORD2)
1038 >    if(coord[1] > MAXBCOORD4)
1039      {
1040        coord[0] <<= 1;
1041 <      coord[1] = (coord[1] << 1) - MAXBCOORD;
1041 >      coord[1] = (coord[1] << 1) - MAXBCOORD2;
1042        coord[2] <<= 1;
1043        return(1);
1044      }
1045      else
1046 <      if(coord[2] > MAXBCOORD2)
1046 >      if(coord[2] > MAXBCOORD4)
1047        {
1048          coord[0] <<= 1;
1049          coord[1] <<= 1;
1050 <        coord[2] = (coord[2] << 1) - MAXBCOORD;
1050 >        coord[2] = (coord[2] << 1) - MAXBCOORD2;
1051          return(2);
1052        }
1053        else
1054           {
1055 <           coord[0] = MAXBCOORD - (coord[0] << 1);
1056 <           coord[1] = MAXBCOORD - (coord[1] << 1);
1057 <           coord[2] = MAXBCOORD - (coord[2] << 1);
1055 >           coord[0] = MAXBCOORD2 - (coord[0] << 1);
1056 >           coord[1] = MAXBCOORD2 - (coord[1] << 1);
1057 >           coord[2] = MAXBCOORD2 - (coord[2] << 1);
1058             return(3);
1059           }
1060   }
1061  
1062 < /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1063 <   dir unbounded to [-MAXLONG,MAXLONG]
1064 < */
1065 < bary_dtol(b,db,bi,dbi,t)
1192 < double b[3],db[3][3];
1193 < BCOORD bi[3];
1194 < BDIR dbi[3][3];
1195 < TINT t[3];
1062 > int
1063 > bary_nth_child(coord,i)
1064 > BCOORD coord[3];
1065 > int i;
1066   {
1197 –  int i,id,j;
1198 –  double d;
1067  
1068 <  for(i=0; i < 2;i++)
1069 <  {
1070 <    if(b[i] <= 0.0)
1071 <    {
1072 <      bi[i]= 0;
1073 <    }
1074 <    else
1075 <      if(b[i] >= 1.0)
1076 <      {
1077 <        bi[i]= MAXBCOORD;
1078 <      }
1079 <      else
1080 <        bi[i] = (BCOORD)(b[i]*MAXBCOORD);
1068 >  switch(i){
1069 >  case 0:
1070 >    /* update bary for child */
1071 >    coord[0] = (coord[0]<< 1) - MAXBCOORD2;
1072 >    coord[1] <<= 1;
1073 >    coord[2] <<= 1;
1074 >    break;
1075 >  case 1:
1076 >    coord[0] <<= 1;
1077 >    coord[1] = (coord[1] << 1) - MAXBCOORD2;
1078 >    coord[2] <<= 1;
1079 >    break;
1080 >  case 2:
1081 >    coord[0] <<= 1;
1082 >    coord[1] <<= 1;
1083 >    coord[2] = (coord[2] << 1) - MAXBCOORD2;
1084 >    break;
1085 >  case 3:
1086 >    coord[0] = MAXBCOORD2 - (coord[0] << 1);
1087 >    coord[1] = MAXBCOORD2 - (coord[1] << 1);
1088 >    coord[2] = MAXBCOORD2 - (coord[2] << 1);
1089 >    break;
1090    }
1214 –  bi[2] = MAXBCOORD - bi[0] - bi[1];
1215 –
1216 –  if(bi[2] < 0)
1217 –  {
1218 –      bi[2] = 0;
1219 –      bi[1] = MAXBCOORD - bi[0];
1220 –  }
1221 –  for(j=0; j< 3; j++)
1222 –  {
1223 –    if(t[j]==0)
1224 –      continue;
1225 –    if(t[j] == HUGET)
1226 –       max_index(db[j],&d);
1227 –    for(i=0; i< 3; i++)
1228 –       if(t[j] != HUGET)
1229 –          dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1230 –       else
1231 –          dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1232 –  }
1091   }
1234 –
1235 –
1236 –
1237 –
1238 –
1239 –
1092  
1093  
1094  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines