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.8 by gwlarson, Wed Nov 11 12:05:38 1998 UTC vs.
Revision 3.12 by gwlarson, Thu Jun 10 15:22:22 1999 UTC

# Line 36 | Line 36 | int
36   convex_angle(v0,v1,v2)
37   FVECT v0,v1,v2;
38   {
39 <    FVECT cp,cp01,cp12,v10,v02;
40 <    double dp;
39 >    double dp,dp1;
40 >    int test,test1;
41 >    FVECT nv0,nv1,nv2;
42 >    FVECT cp01,cp12,cp;
43  
44      /* test sign of (v0Xv1)X(v1Xv2). v1 */
45 <    VCROSS(cp01,v0,v1);
46 <    VCROSS(cp12,v1,v2);
45 >    /* un-Simplified: */
46 >    VCOPY(nv0,v0);
47 >    normalize(nv0);
48 >    VCOPY(nv1,v1);
49 >    normalize(nv1);
50 >    VCOPY(nv2,v2);
51 >    normalize(nv2);
52 >
53 >    VCROSS(cp01,nv0,nv1);
54 >    VCROSS(cp12,nv1,nv2);
55      VCROSS(cp,cp01,cp12);
56 <        
57 <    dp = DOT(cp,v1);
58 < #if 0
59 <    VCOPY(Norm[Ncnt++],cp01);
60 <    VCOPY(Norm[Ncnt++],cp12);
61 <    VCOPY(Norm[Ncnt++],cp);
62 < #endif
63 <    if(ZERO(dp) || dp < 0.0)
64 <      return(FALSE);
65 <    return(TRUE);
56 >    normalize(cp);
57 >    dp1 = DOT(cp,nv1);
58 >    if(dp1 <= 1e-8 || dp1 >= (1-1e-8)) /* Test if on other side,or colinear*/
59 >      test1 = FALSE;
60 >    else
61 >      test1 = TRUE;
62 >
63 >    dp =   nv0[2]*nv1[0]*nv2[1] -  nv0[2]*nv1[1]*nv2[0] - nv0[0]*nv1[2]*nv2[1]
64 >         + nv0[0]*nv1[1]*nv2[2] +  nv0[1]*nv1[2]*nv2[0] - nv0[1]*nv1[0]*nv2[2];
65 >    
66 >    if(dp <= 1e-8 || dp1 >= (1-1e-8))
67 >      test = FALSE;
68 >    else
69 >      test = TRUE;
70 >
71 >    if(test != test1)
72 >      fprintf(stderr,"test %f simplified  %f\n",dp1,dp);
73 >    return(test1);
74   }
75  
76   /* calculates the normal of a face contour using Newell's formula. e
77  
78 <               a =  SUMi (yi - yi+1)(zi + zi+1)smMesh->samples->max_samp+4);
78 >               a =  SUMi (yi - yi+1)(zi + zi+1);
79                 b =  SUMi (zi - zi+1)(xi + xi+1)
80                 c =  SUMi (xi - xi+1)(yi + yi+1)
81   */
# Line 74 | Line 92 | int norm;
92    
93    n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
94             (v1[2] - v2[2]) * (v1[0] + v2[0]) +
95 <           (v2[2] - v0[2]) * (v2[0] + v0[0]);
95 >            (v2[2] - v0[2]) * (v2[0] + v0[0]);
96    
97    n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
98           (v1[1] + v2[1]) * (v1[0] - v2[0]) +
# Line 90 | Line 108 | int norm;
108  
109  
110   tri_plane_equation(v0,v1,v2,peqptr,norm)
111 <   FVECT v0,v1,v2;
111 >   FVECT v0,v1,v2;
112     FPEQ *peqptr;
113     int norm;
114   {
# Line 99 | Line 117 | tri_plane_equation(v0,v1,v2,peqptr,norm)
117      FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
118   }
119  
102 – /* From quad_edge-code */
103 – int
104 – point_in_circle_thru_origin(p,p0,p1)
105 – FVECT p;
106 – FVECT p0,p1;
107 – {
120  
109 –    double dp0,dp1;
110 –    double dp,det;
111 –    
112 –    dp0 = DOT_VEC2(p0,p0);
113 –    dp1 = DOT_VEC2(p1,p1);
114 –    dp  = DOT_VEC2(p,p);    
115 –
116 –    det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
117 –    
118 –    return (det > 0);
119 – }
120 –
121 –
122 –
123 – point_on_sphere(ps,p,c)
124 – FVECT ps,p,c;
125 – {
126 –    VSUB(ps,p,c);    
127 –    normalize(ps);
128 – }
129 –
130 –
121   /* returns TRUE if ray from origin in direction v intersects plane defined
122    * by normal plane_n, and plane_d. If plane is not parallel- returns
123    * intersection point if r != NULL. If tptr!= NULL returns value of
# Line 247 | Line 237 | intersect_ray_oplane(orig,dir,n,pd,r)
237    return(hit);
238   }
239  
240 + /* Assumption: know crosses plane:dont need to check for 'on' case */
241 + intersect_edge_coord_plane(v0,v1,w,r)
242 + FVECT v0,v1;
243 + int w;
244 + FVECT r;
245 + {
246 +  FVECT dv;
247 +  int wnext;
248 +  double t;
249  
250 +  VSUB(dv,v1,v0);
251 +  t = -v0[w]/dv[w];
252 +  r[w] = 0.0;
253 +  wnext = (w+1)%3;
254 +  r[wnext] = v0[wnext] + dv[wnext]*t;
255 +  wnext = (w+2)%3;
256 +  r[wnext] = v0[wnext] + dv[wnext]*t;
257 + }
258 +
259   int
260   intersect_edge_plane(e0,e1,peq,pd,r)
261     FVECT e0,e1;
# Line 281 | Line 289 | intersect_edge_plane(e0,e1,peq,pd,r)
289    return(hit);
290   }
291  
284 –
292   int
286 – point_in_cone(p,p0,p1,p2)
287 – FVECT p;
288 – FVECT p0,p1,p2;
289 – {
290 –    FVECT np,x_axis,y_axis;
291 –    double d1,d2;
292 –    FPEQ peq;
293 –    
294 –    /* Find the equation of the circle defined by the intersection
295 –       of the cone with the plane defined by p1,p2,p3- project p into
296 –       that plane and do an in-circle test in the plane
297 –     */
298 –    
299 –    /* find the equation of the plane defined by p0-p2 */
300 –    tri_plane_equation(p0,p1,p2,&peq,FALSE);
301 –
302 –    /* define a coordinate system on the plane: the x axis is in
303 –       the direction of np2-np1, and the y axis is calculated from
304 –       n cross x-axis
305 –     */
306 –    /* Project p onto the plane */
307 –    /* NOTE: check this: does sideness check?*/
308 –    if(!intersect_vector_plane(p,peq,NULL,np))
309 –        return(FALSE);
310 –
311 –    /* create coordinate system on  plane: p1-p0 defines the x_axis*/
312 –    VSUB(x_axis,p1,p0);
313 –    normalize(x_axis);
314 –    /* The y axis is  */
315 –    VCROSS(y_axis,FP_N(peq),x_axis);
316 –    normalize(y_axis);
317 –
318 –    VSUB(p1,p1,p0);
319 –    VSUB(p2,p2,p0);
320 –    VSUB(np,np,p0);
321 –    
322 –    p1[0] = VLEN(p1);
323 –    p1[1] = 0;
324 –
325 –    d1 = DOT(p2,x_axis);
326 –    d2 = DOT(p2,y_axis);
327 –    p2[0] = d1;
328 –    p2[1] = d2;
329 –    
330 –    d1 = DOT(np,x_axis);
331 –    d2 = DOT(np,y_axis);
332 –    np[0] = d1;
333 –    np[1] = d2;
334 –
335 –    /* perform the in-circle test in the new coordinate system */
336 –    return(point_in_circle_thru_origin(np,p1,p2));
337 – }
338 –
339 – int
293   point_set_in_stri(v0,v1,v2,p,n,nset,sides)
294   FVECT v0,v1,v2,p;
295   FVECT n[3];
# Line 400 | Line 353 | int sides[3];
353  
354  
355  
403 –
404 – int
405 – point_in_stri(v0,v1,v2,p)
406 – FVECT v0,v1,v2,p;
407 – {
408 –    double d;
409 –    FVECT n;  
410 –
411 –    VCROSS(n,v0,v1);
412 –    /* Test the point for sidedness */
413 –    d  = DOT(n,p);
414 –    if(d > 0.0)
415 –      return(FALSE);
416 –
417 –    /* Test next edge */
418 –    VCROSS(n,v1,v2);
419 –    /* Test the point for sidedness */
420 –    d  = DOT(n,p);
421 –    if(d > 0.0)
422 –       return(FALSE);
423 –
424 –    /* Test next edge */
425 –    VCROSS(n,v2,v0);
426 –    /* Test the point for sidedness */
427 –    d  = DOT(n,p);
428 –    if(d > 0.0)
429 –       return(FALSE);
430 –    /* Must be interior to the pyramid */
431 –    return(GT_INTERIOR);
432 – }
433 –
434 – int
435 – vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
436 – FVECT t0,t1,t2,p0,p1,p2;
437 – int *nset;
438 – FVECT n[3];
439 – FVECT avg;
440 – int pt_sides[3][3];
441 –
442 – {
443 –    int below_plane[3],test;
444 –
445 –    SUM_3VEC3(avg,t0,t1,t2);
446 –    *nset = 0;
447 –    /* Test vertex v[i] against triangle j*/
448 –    /* Check if v[i] lies below plane defined by avg of 3 vectors
449 –       defining triangle
450 –       */
451 –
452 –    /* test point 0 */
453 –    if(DOT(avg,p0) < 0.0)
454 –      below_plane[0] = 1;
455 –    else
456 –      below_plane[0] = 0;
457 –    /* Test if b[i] lies in or on triangle a */
458 –    test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
459 –    /* If pts[i] is interior: done */
460 –    if(!below_plane[0])
461 –      {
462 –        if(test == GT_INTERIOR)
463 –          return(TRUE);
464 –      }
465 –    /* Now test point 1*/
466 –
467 –    if(DOT(avg,p1) < 0.0)
468 –      below_plane[1] = 1;
469 –    else
470 –      below_plane[1]=0;
471 –    /* Test if b[i] lies in or on triangle a */
472 –    test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
473 –    /* If pts[i] is interior: done */
474 –    if(!below_plane[1])
475 –    {
476 –      if(test == GT_INTERIOR)
477 –        return(TRUE);
478 –    }
479 –    
480 –    /* Now test point 2 */
481 –    if(DOT(avg,p2) < 0.0)
482 –      below_plane[2] = 1;
483 –    else
484 –      below_plane[2] = 0;
485 –        /* Test if b[i] lies in or on triangle a */
486 –    test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
487 –
488 –    /* If pts[i] is interior: done */
489 –    if(!below_plane[2])
490 –      {
491 –        if(test == GT_INTERIOR)
492 –          return(TRUE);
493 –      }
494 –
495 –    /* If all three points below separating plane: trivial reject */
496 –    if(below_plane[0] && below_plane[1] && below_plane[2])
497 –       return(FALSE);
498 –    /* Now check vertices in a against triangle b */
499 –    return(FALSE);
500 – }
501 –
502 –
356   set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
357     FVECT t0,t1,t2,p0,p1,p2;
358     int test[3];
# Line 639 | Line 492 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
492          SET_NTH_BIT(test[2],2);
493   }
494  
495 + double
496 + point_on_sphere(ps,p,c)
497 + FVECT ps,p,c;
498 + {
499 +  double d;
500 +    VSUB(ps,p,c);    
501 +    d= normalize(ps);
502 +    return(d);
503 + }
504  
505   int
506 < stri_intersect(a1,a2,a3,b1,b2,b3)
507 < FVECT a1,a2,a3,b1,b2,b3;
506 > point_in_stri(v0,v1,v2,p)
507 > FVECT v0,v1,v2,p;
508   {
509 <  int which,test,n_set[2];
510 <  int sides[2][3][3],i,j,inext,jnext;
649 <  int tests[2][3];
650 <  FVECT n[2][3],p,avg[2];
651 <
652 <  /* Test the vertices of triangle a against the pyramid formed by triangle
653 <     b and the origin. If any vertex of a is interior to triangle b, or
654 <     if all 3 vertices of a are ON the edges of b,return TRUE. Remember
655 <     the results of the edge normal and sidedness tests for later.
656 <   */
657 < if(vertices_in_stri(a1,a2,a3,b1,b2,b3,&(n_set[0]),n[0],avg[0],sides[1]))
658 <     return(TRUE);
659 <  
660 < if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
661 <     return(TRUE);
509 >    double d;
510 >    FVECT n;  
511  
512 +    VCROSS(n,v0,v1);
513 +    /* Test the point for sidedness */
514 +    d  = DOT(n,p);
515 +    if(d > 0.0)
516 +      return(FALSE);
517  
518 <  set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
519 <  if(tests[0][0]&tests[0][1]&tests[0][2])
520 <    return(FALSE);
518 >    /* Test next edge */
519 >    VCROSS(n,v1,v2);
520 >    /* Test the point for sidedness */
521 >    d  = DOT(n,p);
522 >    if(d > 0.0)
523 >       return(FALSE);
524  
525 <  set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
526 <  if(tests[1][0]&tests[1][1]&tests[1][2])
527 <    return(FALSE);
528 <
529 <  for(j=0; j < 3;j++)
530 <  {
531 <    jnext = (j+1)%3;
532 <    /* IF edge b doesnt cross any great circles of a, punt */
676 <    if(tests[1][j] & tests[1][jnext])
677 <      continue;
678 <    for(i=0;i<3;i++)
679 <    {
680 <      inext = (i+1)%3;
681 <      /* IF edge a doesnt cross any great circles of b, punt */
682 <      if(tests[0][i] & tests[0][inext])
683 <        continue;
684 <      /* Now find the great circles that cross and test */
685 <      if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
686 <          && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
687 <      {
688 <        VCROSS(p,n[0][i],n[1][j]);
689 <                    
690 <        /* If zero cp= done */
691 <        if(ZERO_VEC3(p))
692 <          continue;
693 <        /* check above both planes */
694 <        if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
695 <          {
696 <            NEGATE_VEC3(p);
697 <            if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
698 <              continue;
699 <          }
700 <        return(TRUE);
701 <      }
702 <    }
703 <  }
704 <  return(FALSE);
525 >    /* Test next edge */
526 >    VCROSS(n,v2,v0);
527 >    /* Test the point for sidedness */
528 >    d  = DOT(n,p);
529 >    if(d > 0.0)
530 >       return(FALSE);
531 >    /* Must be interior to the pyramid */
532 >    return(GT_INTERIOR);
533   }
534  
535 +
536   int
537   ray_intersect_tri(orig,dir,v0,v1,v2,pt)
538   FVECT orig,dir;
# Line 823 | Line 652 | int p0id,p1id,p2id;
652      return(i);
653   }
654  
826 –
827 – int
828 – sedge_intersect(a0,a1,b0,b1)
829 – FVECT a0,a1,b0,b1;
830 – {
831 –    FVECT na,nb,avga,avgb,p;
832 –    double d;
833 –    int sb0,sb1,sa0,sa1;
834 –
835 –    /* First test if edge b straddles great circle of a */
836 –    VCROSS(na,a0,a1);
837 –    d = DOT(na,b0);
838 –    sb0 = ZERO(d)?0:(d<0)? -1: 1;
839 –    d = DOT(na,b1);
840 –    sb1 = ZERO(d)?0:(d<0)? -1: 1;
841 –    /* edge b entirely on one side of great circle a: edges cannot intersect*/
842 –    if(sb0*sb1 > 0)
843 –       return(FALSE);
844 –    /* test if edge a straddles great circle of b */
845 –    VCROSS(nb,b0,b1);
846 –    d = DOT(nb,a0);
847 –    sa0 = ZERO(d)?0:(d<0)? -1: 1;
848 –    d = DOT(nb,a1);
849 –    sa1 = ZERO(d)?0:(d<0)? -1: 1;
850 –    /* edge a entirely on one side of great circle b: edges cannot intersect*/
851 –    if(sa0*sa1 > 0)
852 –       return(FALSE);
853 –
854 –    /* Find one of intersection points of the great circles */
855 –    VCROSS(p,na,nb);
856 –    /* If they lie on same great circle: call an intersection */
857 –    if(ZERO_VEC3(p))
858 –       return(TRUE);
859 –
860 –    VADD(avga,a0,a1);
861 –    VADD(avgb,b0,b1);
862 –    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
863 –    {
864 –      NEGATE_VEC3(p);
865 –      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
866 –        return(FALSE);
867 –    }
868 –    if((!sb0 || !sb1) && (!sa0 || !sa1))
869 –      return(FALSE);
870 –    return(TRUE);
871 – }
872 –
873 –
874 –
655   /* Find the normalized barycentric coordinates of p relative to
656   * triangle v0,v1,v2. Return result in coord
657   */
# Line 889 | Line 669 | double coord[3];
669  
670   }
671  
892 – bary_ith_child(coord,i)
893 – double coord[3];
894 – int i;
895 – {
896 –
897 –  switch(i){
898 –  case 0:
899 –      /* update bary for child */
900 –      coord[0] = 2.0*coord[0]- 1.0;
901 –      coord[1] *= 2.0;
902 –      coord[2] *= 2.0;
903 –      break;
904 –  case 1:
905 –    coord[0] *= 2.0;
906 –    coord[1] = 2.0*coord[1]- 1.0;
907 –    coord[2] *= 2.0;
908 –    break;
909 –  case 2:
910 –    coord[0] *= 2.0;
911 –    coord[1] *= 2.0;
912 –    coord[2] = 2.0*coord[2]- 1.0;
913 –    break;
914 –  case 3:
915 –    coord[0] = 1.0 - 2.0*coord[0];
916 –    coord[1] = 1.0 - 2.0*coord[1];
917 –    coord[2] = 1.0 - 2.0*coord[2];
918 –    break;
919 – #ifdef DEBUG
920 –  default:
921 –    eputs("bary_ith_child():Invalid child\n");
922 –    break;
923 – #endif
924 –  }
925 – }
672  
673  
928 – int
929 – bary_child(coord)
930 – double coord[3];
931 – {
932 –  int i;
674  
934 –  if(coord[0] > 0.5)
935 –  {
936 –      /* update bary for child */
937 –      coord[0] = 2.0*coord[0]- 1.0;
938 –      coord[1] *= 2.0;
939 –      coord[2] *= 2.0;
940 –      return(0);
941 –  }
942 –  else
943 –    if(coord[1] > 0.5)
944 –    {
945 –      coord[0] *= 2.0;
946 –      coord[1] = 2.0*coord[1]- 1.0;
947 –      coord[2] *= 2.0;
948 –      return(1);
949 –    }
950 –    else
951 –      if(coord[2] > 0.5)
952 –      {
953 –        coord[0] *= 2.0;
954 –        coord[1] *= 2.0;
955 –        coord[2] = 2.0*coord[2]- 1.0;
956 –        return(2);
957 –      }
958 –      else
959 –         {
960 –           coord[0] = 1.0 - 2.0*coord[0];
961 –           coord[1] = 1.0 - 2.0*coord[1];
962 –           coord[2] = 1.0 - 2.0*coord[2];
963 –           return(3);
964 –         }
965 – }
966 –
967 – /* Coord was the ith child of the parent: set the coordinate
968 –   relative to the parent
969 – */
675   bary_parent(coord,i)
971 – double coord[3];
972 – int i;
973 – {
974 –
975 –  switch(i) {
976 –  case 0:
977 –    /* update bary for child */
978 –    coord[0] = coord[0]*0.5 + 0.5;
979 –    coord[1] *= 0.5;
980 –    coord[2] *= 0.5;
981 –    break;
982 –  case 1:
983 –    coord[0] *= 0.5;
984 –    coord[1]  = coord[1]*0.5 + 0.5;
985 –    coord[2] *= 0.5;
986 –    break;
987 –    
988 –  case 2:
989 –    coord[0] *= 0.5;
990 –    coord[1] *= 0.5;
991 –    coord[2] = coord[2]*0.5 + 0.5;
992 –    break;
993 –    
994 –  case 3:
995 –    coord[0] = 0.5 - 0.5*coord[0];
996 –    coord[1] = 0.5 - 0.5*coord[1];
997 –    coord[2] = 0.5 - 0.5*coord[2];
998 –    break;
999 – #ifdef DEBUG
1000 –  default:
1001 –    eputs("bary_parent():Invalid child\n");
1002 –    break;
1003 – #endif
1004 –  }
1005 – }
1006 –
1007 – bary_from_child(coord,child,next)
1008 – double coord[3];
1009 – int child,next;
1010 – {
1011 – #ifdef DEBUG
1012 –  if(child <0 || child > 3)
1013 –  {
1014 –    eputs("bary_from_child():Invalid child\n");
1015 –    return;
1016 –  }
1017 –  if(next <0 || next > 3)
1018 –  {
1019 –    eputs("bary_from_child():Invalid next\n");
1020 –    return;
1021 –  }
1022 – #endif
1023 –  if(next == child)
1024 –    return;
1025 –
1026 –  switch(child){
1027 –  case 0:
1028 –    switch(next){
1029 –    case 1:
1030 –      coord[0] += 1.0;
1031 –      coord[1] -= 1.0;
1032 –      break;
1033 –    case 2:
1034 –      coord[0] += 1.0;
1035 –      coord[2] -= 1.0;
1036 –      break;
1037 –    case 3:
1038 –      coord[0] *= -1.0;
1039 –      coord[1] = 1 - coord[1];
1040 –      coord[2] = 1 - coord[2];
1041 –      break;
1042 –
1043 –    }
1044 –    break;
1045 –  case 1:
1046 –    switch(next){
1047 –    case 0:
1048 –      coord[0] -= 1.0;
1049 –      coord[1] += 1.0;
1050 –      break;
1051 –    case 2:
1052 –      coord[1] += 1.0;
1053 –      coord[2] -= 1.0;
1054 –      break;
1055 –    case 3:
1056 –      coord[0] = 1 - coord[0];
1057 –      coord[1] *= -1.0;
1058 –      coord[2] = 1 - coord[2];
1059 –      break;
1060 –    }
1061 –    break;
1062 –  case 2:
1063 –    switch(next){
1064 –    case 0:
1065 –      coord[0] -= 1.0;
1066 –      coord[2] += 1.0;
1067 –      break;
1068 –    case 1:
1069 –      coord[1] -= 1.0;
1070 –      coord[2] += 1.0;
1071 –      break;
1072 –    case 3:
1073 –      coord[0] = 1 - coord[0];
1074 –      coord[1] = 1 - coord[1];
1075 –      coord[2] *= -1.0;
1076 –      break;
1077 –    }
1078 –    break;
1079 –  case 3:
1080 –    switch(next){
1081 –    case 0:
1082 –      coord[0] *= -1.0;
1083 –      coord[1] = 1 - coord[1];
1084 –      coord[2] = 1 - coord[2];
1085 –      break;
1086 –    case 1:
1087 –      coord[0] = 1 - coord[0];
1088 –      coord[1] *= -1.0;
1089 –      coord[2] = 1 - coord[2];
1090 –      break;
1091 –    case 2:
1092 –      coord[0] = 1 - coord[0];
1093 –      coord[1] = 1 - coord[1];
1094 –      coord[2] *= -1.0;
1095 –      break;
1096 –    }
1097 –    break;
1098 –  }
1099 – }
1100 –
1101 –
1102 – baryi_parent(coord,i)
676   BCOORD coord[3];
677   int i;
678   {
1106 –
679    switch(i) {
680    case 0:
681      /* update bary for child */
682 <    coord[0] = (coord[0] >> 1) + MAXBCOORD2;
682 >    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
683      coord[1] >>= 1;
684      coord[2] >>= 1;
685      break;
686    case 1:
687      coord[0] >>= 1;
688 <    coord[1]  = (coord[1] >> 1) + MAXBCOORD2;
688 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
689      coord[2] >>= 1;
690      break;
691      
692    case 2:
693      coord[0] >>= 1;
694      coord[1] >>= 1;
695 <    coord[2] = (coord[2] >> 1) + MAXBCOORD2;
695 >    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
696      break;
697      
698    case 3:
699 <    coord[0] = MAXBCOORD2 - (coord[0] >> 1);
700 <    coord[1] = MAXBCOORD2 - (coord[1] >> 1);
701 <    coord[2] = MAXBCOORD2 - (coord[2] >> 1);
699 >    coord[0] = MAXBCOORD4 - (coord[0] >> 1);
700 >    coord[1] = MAXBCOORD4 - (coord[1] >> 1);
701 >    coord[2] = MAXBCOORD4 - (coord[2] >> 1);
702      break;
703   #ifdef DEBUG
704    default:
705 <    eputs("baryi_parent():Invalid child\n");
705 >    eputs("bary_parent():Invalid child\n");
706      break;
707   #endif
708    }
709   }
710  
711 < baryi_from_child(coord,child,next)
711 > bary_from_child(coord,child,next)
712   BCOORD coord[3];
713   int child,next;
714   {
715   #ifdef DEBUG
716    if(child <0 || child > 3)
717    {
718 <    eputs("baryi_from_child():Invalid child\n");
718 >    eputs("bary_from_child():Invalid child\n");
719      return;
720    }
721    if(next <0 || next > 3)
722    {
723 <    eputs("baryi_from_child():Invalid next\n");
723 >    eputs("bary_from_child():Invalid next\n");
724      return;
725    }
726   #endif
# Line 1158 | Line 730 | int child,next;
730    switch(child){
731    case 0:
732        coord[0] = 0;
733 <      coord[1] = MAXBCOORD - coord[1];
734 <      coord[2] = MAXBCOORD - coord[2];
733 >      coord[1] = MAXBCOORD2 - coord[1];
734 >      coord[2] = MAXBCOORD2 - coord[2];
735        break;
736    case 1:
737 <      coord[0] = MAXBCOORD - coord[0];
737 >      coord[0] = MAXBCOORD2 - coord[0];
738        coord[1] = 0;
739 <      coord[2] = MAXBCOORD - coord[2];
739 >      coord[2] = MAXBCOORD2 - coord[2];
740        break;
741    case 2:
742 <      coord[0] = MAXBCOORD - coord[0];
743 <      coord[1] = MAXBCOORD - coord[1];
742 >      coord[0] = MAXBCOORD2 - coord[0];
743 >      coord[1] = MAXBCOORD2 - coord[1];
744        coord[2] = 0;
745      break;
746    case 3:
747      switch(next){
748      case 0:
749        coord[0] = 0;
750 <      coord[1] = MAXBCOORD - coord[1];
751 <      coord[2] = MAXBCOORD - coord[2];
750 >      coord[1] = MAXBCOORD2 - coord[1];
751 >      coord[2] = MAXBCOORD2 - coord[2];
752        break;
753      case 1:
754 <      coord[0] = MAXBCOORD - coord[0];
754 >      coord[0] = MAXBCOORD2 - coord[0];
755        coord[1] = 0;
756 <      coord[2] = MAXBCOORD - coord[2];
756 >      coord[2] = MAXBCOORD2 - coord[2];
757        break;
758      case 2:
759 <      coord[0] = MAXBCOORD - coord[0];
760 <      coord[1] = MAXBCOORD - coord[1];
759 >      coord[0] = MAXBCOORD2 - coord[0];
760 >      coord[1] = MAXBCOORD2 - coord[1];
761        coord[2] = 0;
762        break;
763      }
# Line 1194 | Line 766 | int child,next;
766   }
767  
768   int
769 < baryi_child(coord)
769 > bary_child(coord)
770   BCOORD coord[3];
771   {
772    int i;
773  
774 <  if(coord[0] > MAXBCOORD2)
774 >  if(coord[0] > MAXBCOORD4)
775    {
776        /* update bary for child */
777 <      coord[0] = (coord[0]<< 1) - MAXBCOORD;
777 >      coord[0] = (coord[0]<< 1) - MAXBCOORD2;
778        coord[1] <<= 1;
779        coord[2] <<= 1;
780        return(0);
781    }
782    else
783 <    if(coord[1] > MAXBCOORD2)
783 >    if(coord[1] > MAXBCOORD4)
784      {
785        coord[0] <<= 1;
786 <      coord[1] = (coord[1] << 1) - MAXBCOORD;
786 >      coord[1] = (coord[1] << 1) - MAXBCOORD2;
787        coord[2] <<= 1;
788        return(1);
789      }
790      else
791 <      if(coord[2] > MAXBCOORD2)
791 >      if(coord[2] > MAXBCOORD4)
792        {
793          coord[0] <<= 1;
794          coord[1] <<= 1;
795 <        coord[2] = (coord[2] << 1) - MAXBCOORD;
795 >        coord[2] = (coord[2] << 1) - MAXBCOORD2;
796          return(2);
797        }
798        else
799           {
800 <           coord[0] = MAXBCOORD - (coord[0] << 1);
801 <           coord[1] = MAXBCOORD - (coord[1] << 1);
802 <           coord[2] = MAXBCOORD - (coord[2] << 1);
800 >           coord[0] = MAXBCOORD2 - (coord[0] << 1);
801 >           coord[1] = MAXBCOORD2 - (coord[1] << 1);
802 >           coord[2] = MAXBCOORD2 - (coord[2] << 1);
803             return(3);
804           }
805   }
806  
807   int
808 < baryi_nth_child(coord,i)
808 > bary_nth_child(coord,i)
809   BCOORD coord[3];
810   int i;
811   {
# Line 1241 | Line 813 | int i;
813    switch(i){
814    case 0:
815      /* update bary for child */
816 <    coord[0] = (coord[0]<< 1) - MAXBCOORD;
816 >    coord[0] = (coord[0]<< 1) - MAXBCOORD2;
817      coord[1] <<= 1;
818      coord[2] <<= 1;
819      break;
820    case 1:
821      coord[0] <<= 1;
822 <    coord[1] = (coord[1] << 1) - MAXBCOORD;
822 >    coord[1] = (coord[1] << 1) - MAXBCOORD2;
823      coord[2] <<= 1;
824      break;
825    case 2:
826      coord[0] <<= 1;
827      coord[1] <<= 1;
828 <    coord[2] = (coord[2] << 1) - MAXBCOORD;
828 >    coord[2] = (coord[2] << 1) - MAXBCOORD2;
829      break;
830    case 3:
831 <    coord[0] = MAXBCOORD - (coord[0] << 1);
832 <    coord[1] = MAXBCOORD - (coord[1] << 1);
833 <    coord[2] = MAXBCOORD - (coord[2] << 1);
831 >    coord[0] = MAXBCOORD2 - (coord[0] << 1);
832 >    coord[1] = MAXBCOORD2 - (coord[1] << 1);
833 >    coord[2] = MAXBCOORD2 - (coord[2] << 1);
834      break;
835    }
836   }
1265 –
1266 –
1267 – baryi_children(coord,i,in_tri,rcoord,rin_tri)
1268 – BCOORD coord[3][3];
1269 – int i;
1270 – int in_tri[3];
1271 – BCOORD rcoord[3][3];
1272 – int rin_tri[3];
1273 – {
1274 –  int j;
1275 –
1276 –  for(j=0; j< 3; j++)
1277 –  {
1278 –    if(!in_tri[j])
1279 –    {
1280 –      rin_tri[j]=0;
1281 –      continue;
1282 –    }
1283 –    
1284 –    if(i != 3)
1285 –    {
1286 –      if(coord[j][i] < MAXBCOORD2)
1287 –        {
1288 –          rin_tri[j] = 0;
1289 –          continue;
1290 –        }
1291 –    }
1292 –    else
1293 –      if( !(coord[j][0] <= MAXBCOORD2 && coord[j][1] <= MAXBCOORD2 &&
1294 –            coord[j][2] <= MAXBCOORD2))
1295 –        {
1296 –          rin_tri[j] = 0;
1297 –          continue;
1298 –        }
1299 –      
1300 –    rin_tri[j]=1;
1301 –    VCOPY(rcoord[j],coord[j]);
1302 –    baryi_nth_child(rcoord[j],i);
1303 –  }
1304 –
1305 – }
1306 –
1307 – convert_dtol(b,bi)
1308 – double b[3];
1309 – BCOORD bi[3];
1310 – {
1311 –  int i;
1312 –
1313 –  for(i=0; i < 2;i++)
1314 –  {
1315 –
1316 –    if(b[i] <= 0.0)
1317 –    {
1318 –      bi[i]= 0;
1319 –    }
1320 –    else
1321 –      if(b[i] >= 1.0)
1322 –      {
1323 –        bi[i]= MAXBCOORD;
1324 –      }
1325 –      else
1326 –        bi[i] = (BCOORD)(b[i]*MAXBCOORD);
1327 –  }
1328 –  bi[2] = bi[0] +  bi[1];
1329 –  if(bi[2] > MAXBCOORD)
1330 –  {
1331 –      bi[2] = 0;
1332 –      bi[1] = MAXBCOORD - bi[0];
1333 –  }
1334 –  else
1335 –    bi[2] = MAXBCOORD - bi[2];
1336 –
1337 – }
1338 –
1339 – /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1340 –   dir unbounded to [-MAXLONG,MAXLONG]
1341 – */
1342 – bary_dtol(b,db,bi,dbi,t,w)
1343 – double b[3],db[3][3];
1344 – BCOORD bi[3];
1345 – BDIR dbi[3][3];
1346 – TINT t[3];
1347 – int w;
1348 – {
1349 –  int i,id,j,k;
1350 –  double d;
1351 –
1352 –  convert_dtol(b,bi);
1353 –
1354 –  for(j=w; j< 3; j++)
1355 –  {
1356 –    if(t[j] == HUGET)
1357 –    {
1358 –      max_index(db[j],&d);
1359 –      for(i=0; i< 3; i++)
1360 –        dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1361 –      break;
1362 –    }
1363 –    else
1364 –    {
1365 –      for(i=0; i< 3; i++)
1366 –          dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1367 –    }
1368 –  }
1369 – }
1370 –
1371 –
1372 – /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1373 –   dir unbounded to [-MAXLONG,MAXLONG]
1374 – */
1375 – bary_dtol_new(b,db,bi,boi,dbi,t)
1376 – double b[4][3],db[3][3];
1377 – BCOORD bi[3],boi[3][3];
1378 – BDIR dbi[3][3];
1379 – int t[3];
1380 – {
1381 –  int i,id,j,k;
1382 –  double d;
1383 –
1384 –  convert_dtol(b[3],bi);
1385 –
1386 –  for(j=0; j<3;j++)
1387 –  {
1388 –    if(t[j] != 1)
1389 –      continue;
1390 –    convert_dtol(b[j],boi[j]);
1391 –  }
1392 –  for(j=0; j< 3; j++)
1393 –  {
1394 –    k = (j+1)%3;
1395 –    if(t[k]==0)
1396 –      continue;
1397 –    if(t[k] == -1)
1398 –      {
1399 –        max_index(db[j],&d);
1400 –        for(i=0; i< 3; i++)
1401 –          dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1402 –        t[k] = 0;
1403 –      }
1404 –    else
1405 –      if(t[j] != 1)
1406 –        for(i=0; i< 3; i++)
1407 –          dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1408 –    else
1409 –      for(i=0; i< 3; i++)
1410 –        dbi[j][i] = boi[k][i] - boi[j][i];
1411 –    
1412 –  }
1413 – }
1414 –
1415 –
1416 – bary_dtolb(b,bi,in_tri)
1417 – double b[3][3];
1418 – BCOORD bi[3][3];
1419 – int in_tri[3];
1420 – {
1421 –  int i,j;
1422 –
1423 –  for(j=0; j<3;j++)
1424 –  {
1425 –    if(!in_tri[j])
1426 –      continue;
1427 –    for(i=0; i < 2;i++)
1428 –    {
1429 –    if(b[j][i] <= 0.0)
1430 –    {
1431 –      bi[j][i]= 0;
1432 –    }
1433 –    else
1434 –      if(b[j][i] >= 1.0)
1435 –      {
1436 –        bi[j][i]= MAXBCOORD;
1437 –      }
1438 –      else
1439 –        bi[j][i] = (BCOORD)(b[j][i]*MAXBCOORD);
1440 –    }
1441 –    bi[j][2] = MAXBCOORD - bi[j][0] - bi[j][1];
1442 –    if(bi[j][2] < 0)
1443 –      {
1444 –        bi[j][2] = 0;
1445 –        bi[j][1] = MAXBCOORD - bi[j][0];
1446 –      }
1447 –  }
1448 – }
1449 –
1450 –
1451 –
1452 –
1453 –
1454 –
1455 –
1456 –
837  
838  
839  

Diff Legend

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