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.4 by gwlarson, Fri Sep 11 11:52:25 1998 UTC vs.
Revision 3.8 by gwlarson, Wed Nov 11 12:05:38 1998 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)smMesh->samples->max_samp+4);
61                 b =  SUMi (zi - zi+1)(xi + xi+1)
62                 c =  SUMi (xi - xi+1)(yi + yi+1)
63   */
# Line 64 | Line 75 | int norm;
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 –
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 125 | Line 134 | FVECT ps,p,c;
134    * t, if parallel, returns t=FHUGE
135    */
136   int
137 < intersect_vector_plane(v,plane_n,plane_d,tptr,r)
138 <   FVECT v,plane_n;
139 <   double plane_d;
137 > intersect_vector_plane(v,peq,tptr,r)
138 >   FVECT v;
139 >   FPEQ peq;
140     double *tptr;
141     FVECT r;
142   {
# Line 141 | Line 150 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
150      /* line is  l = p1 + (p2-p1)t, p1=origin */
151  
152      /* Solve for t: */
153 <  d = -(DOT(plane_n,v));
153 >  d = -(DOT(FP_N(peq),v));
154    if(ZERO(d))
155    {
156        t = FHUGE;
# Line 149 | Line 158 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
158    }
159    else
160    {
161 <      t =  plane_d/d;
161 >      t =  FP_D(peq)/d;
162        if(t < 0 )
163           hit = 0;
164        else
# Line 167 | Line 176 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
176   }
177  
178   int
179 < intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
179 > intersect_ray_plane(orig,dir,peq,pd,r)
180     FVECT orig,dir;
181 <   FVECT plane_n;
173 <   double plane_d;
181 >   FPEQ peq;
182     double *pd;
183     FVECT r;
184   {
185 <  double t;
185 >  double t,d;
186    int hit;
187      /*
188        Plane is Ax + By + Cz +D = 0:
# Line 185 | Line 193 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
193         line is  l = p1 + (p2-p1)t
194       */
195      /* Solve for t: */
196 <    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
196 >  d = DOT(FP_N(peq),dir);
197 >  if(ZERO(d))
198 >     return(0);
199 >  t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
200 >
201 >  if(t < 0)
202 >       hit = 0;
203 >    else
204 >       hit = 1;
205 >
206 >  if(r)
207 >     VSUM(r,orig,dir,t);
208 >
209 >    if(pd)
210 >       *pd = t;
211 >  return(hit);
212 > }
213 >
214 >
215 > int
216 > intersect_ray_oplane(orig,dir,n,pd,r)
217 >   FVECT orig,dir;
218 >   FVECT n;
219 >   double *pd;
220 >   FVECT r;
221 > {
222 >  double t,d;
223 >  int hit;
224 >    /*
225 >      Plane is Ax + By + Cz +D = 0:
226 >      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
227 >    */
228 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
229 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
230 >       line is  l = p1 + (p2-p1)t
231 >     */
232 >    /* Solve for t: */
233 >    d= DOT(n,dir);
234 >    if(ZERO(d))
235 >       return(0);
236 >    t =  -(DOT(n,orig))/d;
237      if(t < 0)
238         hit = 0;
239      else
# Line 201 | Line 249 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
249  
250  
251   int
252 < intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
252 > intersect_edge_plane(e0,e1,peq,pd,r)
253     FVECT e0,e1;
254 <   FVECT plane_n;
207 <   double plane_d;
254 >   FPEQ peq;
255     double *pd;
256     FVECT r;
257   {
# Line 221 | Line 268 | intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
268       */
269      /* Solve for t: */
270    VSUB(d,e1,e0);
271 <  t =  -(DOT(plane_n,e0) + plane_d)/(DOT(plane_n,d));
271 >  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
272      if(t < 0)
273         hit = 0;
274      else
# Line 240 | Line 287 | point_in_cone(p,p0,p1,p2)
287   FVECT p;
288   FVECT p0,p1,p2;
289   {
243 –    FVECT n;
290      FVECT np,x_axis,y_axis;
291 <    double d1,d2,d;
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 p1-p3 */
300 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
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
# Line 258 | Line 305 | FVECT p0,p1,p2;
305       */
306      /* Project p onto the plane */
307      /* NOTE: check this: does sideness check?*/
308 <    if(!intersect_vector_plane(p,n,d,NULL,np))
308 >    if(!intersect_vector_plane(p,peq,NULL,np))
309          return(FALSE);
310  
311 <    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
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,n,x_axis);
315 >    VCROSS(y_axis,FP_N(peq),x_axis);
316      normalize(y_axis);
317  
318      VSUB(p1,p1,p0);
# Line 301 | Line 348 | int sides[3];
348      /* Find the normal to the triangle ORIGIN:v0:v1 */
349      if(!NTH_BIT(*nset,0))
350      {
351 <        VCROSS(n[0],v1,v0);
351 >        VCROSS(n[0],v0,v1);
352          SET_NTH_BIT(*nset,0);
353      }
354      /* Test the point for sidedness */
# Line 319 | Line 366 | int sides[3];
366      /* Test next edge */
367      if(!NTH_BIT(*nset,1))
368      {
369 <        VCROSS(n[1],v2,v1);
369 >        VCROSS(n[1],v1,v2);
370          SET_NTH_BIT(*nset,1);
371      }
372      /* Test the point for sidedness */
# Line 335 | Line 382 | int sides[3];
382      /* Test next edge */
383      if(!NTH_BIT(*nset,2))
384      {
385 <        VCROSS(n[2],v0,v2);
385 >        VCROSS(n[2],v2,v0);
386          SET_NTH_BIT(*nset,2);
387      }
388      /* Test the point for sidedness */
# Line 353 | Line 400 | int sides[3];
400  
401  
402  
403 <
403 >
404   int
405   point_in_stri(v0,v1,v2,p)
406   FVECT v0,v1,v2,p;
# Line 361 | Line 408 | FVECT v0,v1,v2,p;
408      double d;
409      FVECT n;  
410  
411 <    VCROSS(n,v1,v0);
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,v2,v1);
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,v0,v2);
425 >    VCROSS(n,v2,v0);
426      /* Test the point for sidedness */
427      d  = DOT(n,p);
428      if(d > 0.0)
# Line 469 | Line 516 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
516      if(sides[0][0] == GT_INVALID)
517      {
518        if(!NTH_BIT(nset,0))
519 <        VCROSS(n[0],t1,t0);
519 >        VCROSS(n[0],t0,t1);
520        /* Test the point for sidedness */
521        d  = DOT(n[0],p0);
522        if(d >= 0.0)
# Line 482 | Line 529 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
529      if(sides[0][1] == GT_INVALID)
530      {
531        if(!NTH_BIT(nset,1))
532 <        VCROSS(n[1],t2,t1);
532 >        VCROSS(n[1],t1,t2);
533          /* Test the point for sidedness */
534          d  = DOT(n[1],p0);
535          if(d >= 0.0)
# Line 495 | Line 542 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
542      if(sides[0][2] == GT_INVALID)
543      {
544        if(!NTH_BIT(nset,2))
545 <        VCROSS(n[2],t0,t2);
545 >        VCROSS(n[2],t2,t0);
546        /* Test the point for sidedness */
547        d  = DOT(n[2],p0);
548        if(d >= 0.0)
# Line 511 | Line 558 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
558      if(sides[1][0] == GT_INVALID)
559      {
560        if(!NTH_BIT(nset,0))
561 <        VCROSS(n[0],t1,t0);
561 >        VCROSS(n[0],t0,t1);
562        /* Test the point for sidedness */
563        d  = DOT(n[0],p1);
564        if(d >= 0.0)
# Line 525 | Line 572 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
572      if(sides[1][1] == GT_INVALID)
573      {
574        if(!NTH_BIT(nset,1))
575 <        VCROSS(n[1],t2,t1);
575 >        VCROSS(n[1],t1,t2);
576        /* Test the point for sidedness */
577        d  = DOT(n[1],p1);
578        if(d >= 0.0)
# Line 539 | Line 586 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
586      if(sides[1][2] == GT_INVALID)
587      {
588        if(!NTH_BIT(nset,2))
589 <        VCROSS(n[2],t0,t2);
589 >        VCROSS(n[2],t2,t0);
590        /* Test the point for sidedness */
591        d  = DOT(n[2],p1);
592        if(d >= 0.0)
# Line 555 | Line 602 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
602      if(sides[2][0] == GT_INVALID)
603      {
604        if(!NTH_BIT(nset,0))
605 <        VCROSS(n[0],t1,t0);
605 >        VCROSS(n[0],t0,t1);
606        /* Test the point for sidedness */
607        d  = DOT(n[0],p2);
608        if(d >= 0.0)
# Line 568 | Line 615 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
615      if(sides[2][1] == GT_INVALID)
616      {
617        if(!NTH_BIT(nset,1))
618 <        VCROSS(n[1],t2,t1);
618 >        VCROSS(n[1],t1,t2);
619        /* Test the point for sidedness */
620        d  = DOT(n[1],p2);
621        if(d >= 0.0)
# Line 581 | Line 628 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
628      if(sides[2][2] == GT_INVALID)
629      {
630        if(!NTH_BIT(nset,2))
631 <        VCROSS(n[2],t0,t2);
631 >        VCROSS(n[2],t2,t0);
632        /* Test the point for sidedness */
633        d  = DOT(n[2],p2);
634        if(d >= 0.0)
# Line 663 | Line 710 | FVECT orig,dir;
710   FVECT v0,v1,v2;
711   FVECT pt;
712   {
713 <  FVECT p0,p1,p2,p,n;
714 <  double pd;
713 >  FVECT p0,p1,p2,p;
714 >  FPEQ peq;
715    int type;
716  
717 <  point_on_sphere(p0,v0,orig);
718 <  point_on_sphere(p1,v1,orig);
719 <  point_on_sphere(p2,v2,orig);
720 <  
717 >  VSUB(p0,v0,orig);
718 >  VSUB(p1,v1,orig);
719 >  VSUB(p2,v2,orig);
720 >
721    if(point_in_stri(p0,p1,p2,dir))
722    {
723        /* Intersect the ray with the triangle plane */
724 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
725 <      return(intersect_ray_plane(orig,dir,n,pd,NULL,pt));
724 >      tri_plane_equation(v0,v1,v2,&peq,FALSE);
725 >      return(intersect_ray_plane(orig,dir,peq,NULL,pt));
726    }
727    return(FALSE);
728   }
# Line 736 | Line 783 | FVECT fnear[4],ffar[4];
783      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
784   }
785  
786 + int
787 + max_index(v,r)
788 + FVECT v;
789 + double *r;
790 + {
791 +  double p[3];
792 +  int i;
793  
794 +  p[0] = fabs(v[0]);
795 +  p[1] = fabs(v[1]);
796 +  p[2] = fabs(v[2]);
797 +  i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);  
798 +  if(r)
799 +    *r = p[i];
800 +  return(i);
801 + }
802  
803 + int
804 + closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
805 + FVECT p0,p1,p2,p;
806 + int p0id,p1id,p2id;
807 + {
808 +    double d,d1;
809 +    int i;
810 +    
811 +    d =  DIST_SQ(p,p0);
812 +    d1 = DIST_SQ(p,p1);
813 +    if(d < d1)
814 +    {
815 +      d1 = DIST_SQ(p,p2);
816 +      i = (d1 < d)?p2id:p0id;
817 +    }
818 +    else
819 +    {
820 +      d = DIST_SQ(p,p2);
821 +      i = (d < d1)? p2id:p1id;
822 +    }
823 +    return(i);
824 + }
825  
826 +
827   int
828   sedge_intersect(a0,a1,b0,b1)
829   FVECT a0,a1,b0,b1;
# Line 800 | Line 885 | double coord[3];
885    a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
886    coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
887    coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
888 <  coord[2]  = 1.0 - coord[0] - coord[1];
888 >  coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
889  
890   }
891  
# Line 1013 | Line 1098 | int child,next;
1098    }
1099   }
1100  
1101 < int
1102 < max_index(v)
1103 < FVECT v;
1101 >
1102 > baryi_parent(coord,i)
1103 > BCOORD coord[3];
1104 > int i;
1105   {
1020 –  double a,b,c;
1021 –  int i;
1106  
1107 <  a = fabs(v[0]);
1108 <  b = fabs(v[1]);
1109 <  c = fabs(v[2]);
1110 <  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
1111 <  return(i);
1107 >  switch(i) {
1108 >  case 0:
1109 >    /* update bary for child */
1110 >    coord[0] = (coord[0] >> 1) + MAXBCOORD2;
1111 >    coord[1] >>= 1;
1112 >    coord[2] >>= 1;
1113 >    break;
1114 >  case 1:
1115 >    coord[0] >>= 1;
1116 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD2;
1117 >    coord[2] >>= 1;
1118 >    break;
1119 >    
1120 >  case 2:
1121 >    coord[0] >>= 1;
1122 >    coord[1] >>= 1;
1123 >    coord[2] = (coord[2] >> 1) + MAXBCOORD2;
1124 >    break;
1125 >    
1126 >  case 3:
1127 >    coord[0] = MAXBCOORD2 - (coord[0] >> 1);
1128 >    coord[1] = MAXBCOORD2 - (coord[1] >> 1);
1129 >    coord[2] = MAXBCOORD2 - (coord[2] >> 1);
1130 >    break;
1131 > #ifdef DEBUG
1132 >  default:
1133 >    eputs("baryi_parent():Invalid child\n");
1134 >    break;
1135 > #endif
1136 >  }
1137   }
1138  
1139 + baryi_from_child(coord,child,next)
1140 + BCOORD coord[3];
1141 + int child,next;
1142 + {
1143 + #ifdef DEBUG
1144 +  if(child <0 || child > 3)
1145 +  {
1146 +    eputs("baryi_from_child():Invalid child\n");
1147 +    return;
1148 +  }
1149 +  if(next <0 || next > 3)
1150 +  {
1151 +    eputs("baryi_from_child():Invalid next\n");
1152 +    return;
1153 +  }
1154 + #endif
1155 +  if(next == child)
1156 +    return;
1157  
1158 +  switch(child){
1159 +  case 0:
1160 +      coord[0] = 0;
1161 +      coord[1] = MAXBCOORD - coord[1];
1162 +      coord[2] = MAXBCOORD - coord[2];
1163 +      break;
1164 +  case 1:
1165 +      coord[0] = MAXBCOORD - coord[0];
1166 +      coord[1] = 0;
1167 +      coord[2] = MAXBCOORD - coord[2];
1168 +      break;
1169 +  case 2:
1170 +      coord[0] = MAXBCOORD - coord[0];
1171 +      coord[1] = MAXBCOORD - coord[1];
1172 +      coord[2] = 0;
1173 +    break;
1174 +  case 3:
1175 +    switch(next){
1176 +    case 0:
1177 +      coord[0] = 0;
1178 +      coord[1] = MAXBCOORD - coord[1];
1179 +      coord[2] = MAXBCOORD - coord[2];
1180 +      break;
1181 +    case 1:
1182 +      coord[0] = MAXBCOORD - coord[0];
1183 +      coord[1] = 0;
1184 +      coord[2] = MAXBCOORD - coord[2];
1185 +      break;
1186 +    case 2:
1187 +      coord[0] = MAXBCOORD - coord[0];
1188 +      coord[1] = MAXBCOORD - coord[1];
1189 +      coord[2] = 0;
1190 +      break;
1191 +    }
1192 +    break;
1193 +  }
1194 + }
1195  
1032 – /*
1033 – * int
1034 – * traceRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
1035 – *
1036 – *   Intersect the ray with triangle v0v1v2, return intersection point in r
1037 – *
1038 – *    Assumes orig,v0,v1,v2 are in spherical coordinates, and orig is
1039 – *    unit
1040 – */
1196   int
1197 < traceRay(orig,dir,v0,v1,v2,r)
1198 <  FVECT orig,dir;
1044 <  FVECT v0,v1,v2;
1045 <  FVECT r;
1197 > baryi_child(coord)
1198 > BCOORD coord[3];
1199   {
1047 –  FVECT n,p[3],d;
1048 –  double pt[3],r_eps;
1200    int i;
1050 –  int which;
1201  
1202 <  /* Find the plane equation for the triangle defined by the edge v0v1 and
1203 <   the view center*/
1204 <  VCROSS(n,v0,v1);
1205 <  /* Intersect the ray with this plane */
1206 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[0]),p[0]);
1207 <  if(i)
1208 <    which = 0;
1202 >  if(coord[0] > MAXBCOORD2)
1203 >  {
1204 >      /* update bary for child */
1205 >      coord[0] = (coord[0]<< 1) - MAXBCOORD;
1206 >      coord[1] <<= 1;
1207 >      coord[2] <<= 1;
1208 >      return(0);
1209 >  }
1210    else
1211 <    which = -1;
1211 >    if(coord[1] > MAXBCOORD2)
1212 >    {
1213 >      coord[0] <<= 1;
1214 >      coord[1] = (coord[1] << 1) - MAXBCOORD;
1215 >      coord[2] <<= 1;
1216 >      return(1);
1217 >    }
1218 >    else
1219 >      if(coord[2] > MAXBCOORD2)
1220 >      {
1221 >        coord[0] <<= 1;
1222 >        coord[1] <<= 1;
1223 >        coord[2] = (coord[2] << 1) - MAXBCOORD;
1224 >        return(2);
1225 >      }
1226 >      else
1227 >         {
1228 >           coord[0] = MAXBCOORD - (coord[0] << 1);
1229 >           coord[1] = MAXBCOORD - (coord[1] << 1);
1230 >           coord[2] = MAXBCOORD - (coord[2] << 1);
1231 >           return(3);
1232 >         }
1233 > }
1234  
1235 <  VCROSS(n,v1,v2);
1236 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[1]),p[1]);
1237 <  if(i)
1238 <    if((which==-1) || pt[1] < pt[0])
1239 <      which = 1;
1235 > int
1236 > baryi_nth_child(coord,i)
1237 > BCOORD coord[3];
1238 > int i;
1239 > {
1240  
1241 <  VCROSS(n,v2,v0);
1242 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[2]),p[2]);
1243 <  if(i)
1244 <    if((which==-1) || pt[2] < pt[which])
1245 <      which = 2;
1241 >  switch(i){
1242 >  case 0:
1243 >    /* update bary for child */
1244 >    coord[0] = (coord[0]<< 1) - MAXBCOORD;
1245 >    coord[1] <<= 1;
1246 >    coord[2] <<= 1;
1247 >    break;
1248 >  case 1:
1249 >    coord[0] <<= 1;
1250 >    coord[1] = (coord[1] << 1) - MAXBCOORD;
1251 >    coord[2] <<= 1;
1252 >    break;
1253 >  case 2:
1254 >    coord[0] <<= 1;
1255 >    coord[1] <<= 1;
1256 >    coord[2] = (coord[2] << 1) - MAXBCOORD;
1257 >    break;
1258 >  case 3:
1259 >    coord[0] = MAXBCOORD - (coord[0] << 1);
1260 >    coord[1] = MAXBCOORD - (coord[1] << 1);
1261 >    coord[2] = MAXBCOORD - (coord[2] << 1);
1262 >    break;
1263 >  }
1264 > }
1265  
1266 <  if(which != -1)
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 <      /* Return point that is K*eps along projection of the ray on
1279 <         the sphere to push intersection point p[which] into next cell
1280 <      */
1281 <      normalize(p[which]);
1282 <      /* Calculate the ray perpendicular to the intersection point: approx
1283 <       the direction moved along the sphere a distance of K*epsilon*/
1284 <      r_eps = -DOT(p[which],dir);
1285 <      VSUM(n,dir,p[which],r_eps);
1286 <     /* Calculate the length along ray p[which]-dir needed to move to
1287 <         cause a move along the sphere of k*epsilon
1288 <       */
1289 <       r_eps = DOT(n,dir);
1290 <      VSUM(r,p[which],dir,(20*FTINY)/r_eps);
1291 <      normalize(r);
1292 <      return(TRUE);
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 <      /* Unable to find intersection: move along ray and try again */
1357 <      r_eps = -DOT(orig,dir);
1358 <      VSUM(n,dir,orig,r_eps);
1359 <      r_eps = DOT(n,dir);
1360 <      VSUM(r,orig,dir,(20*FTINY)/r_eps);
1361 <      normalize(r);
1362 < #ifdef DEBUG
1363 <      eputs("traceRay:Ray does not intersect triangle\n");
1364 < #endif
1365 <      return(FALSE);
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 < int
1373 < closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
1374 < FVECT p0,p1,p2,p;
1375 < int p0id,p1id,p2id;
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 <    double d,d1;
1382 <    int i;
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 <    d =  DIST_SQ(p,p0);
1413 <    d1 = DIST_SQ(p,p1);
1414 <    if(d < d1)
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 <      d1 = DIST_SQ(p,p2);
1430 <      i = (d1 < d)?p2id:p0id;
1429 >    if(b[j][i] <= 0.0)
1430 >    {
1431 >      bi[j][i]= 0;
1432      }
1433      else
1434 <    {
1435 <      d = DIST_SQ(p,p2);
1436 <      i = (d < d1)? p2id:p1id;
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 <    return(i);
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  

Diff Legend

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