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.7 by gwlarson, Tue Oct 6 18:16:53 1998 UTC

# Line 33 | Line 33 | int
33   convex_angle(v0,v1,v2)
34   FVECT v0,v1,v2;
35   {
36 <    FVECT cp01,cp12,cp;
37 <    
38 <    /* test sign of (v0Xv1)X(v1Xv2). v1 */
36 >    FVECT cp,cp01,cp12,v10,v02;
37 >    double dp;
38 >    /*
39 >      VSUB(v10,v1,v0);
40 >      VSUB(v02,v0,v2);
41 >      VCROSS(cp,v10,v02);
42 >   */
43 >      /* test sign of (v0Xv1)X(v1Xv2). v1 */
44      VCROSS(cp01,v0,v1);
45      VCROSS(cp12,v1,v2);
46      VCROSS(cp,cp01,cp12);
47 <    if(DOT(cp,v1) < 0.0)
47 >        
48 >    dp = DOT(cp,v1);
49 >    if(ZERO(dp) || dp < 0.0)
50         return(FALSE);
51      return(TRUE);
52   }
# Line 64 | Line 71 | int norm;
71    n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
72             (v1[2] - v2[2]) * (v1[0] + v2[0]) +
73             (v2[2] - v0[2]) * (v2[0] + v0[0]);
67 –
74    
75    n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
76           (v1[1] + v2[1]) * (v1[0] - v2[0]) +
# Line 80 | Line 86 | int norm;
86   }
87  
88  
89 < tri_plane_equation(v0,v1,v2,n,nd,norm)
90 <   FVECT v0,v1,v2,n;
91 <   double *nd;
89 > tri_plane_equation(v0,v1,v2,peqptr,norm)
90 >   FVECT v0,v1,v2;
91 >   FPEQ *peqptr;
92     int norm;
93   {
94 <    tri_normal(v0,v1,v2,n,norm);
94 >    tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
95  
96 <    *nd = -(DOT(n,v0));
96 >    FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
97   }
98  
99   /* From quad_edge-code */
# Line 125 | Line 131 | FVECT ps,p,c;
131    * t, if parallel, returns t=FHUGE
132    */
133   int
134 < intersect_vector_plane(v,plane_n,plane_d,tptr,r)
135 <   FVECT v,plane_n;
136 <   double plane_d;
134 > intersect_vector_plane(v,peq,tptr,r)
135 >   FVECT v;
136 >   FPEQ peq;
137     double *tptr;
138     FVECT r;
139   {
# Line 141 | Line 147 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
147      /* line is  l = p1 + (p2-p1)t, p1=origin */
148  
149      /* Solve for t: */
150 <  d = -(DOT(plane_n,v));
150 >  d = -(DOT(FP_N(peq),v));
151    if(ZERO(d))
152    {
153        t = FHUGE;
# Line 149 | Line 155 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
155    }
156    else
157    {
158 <      t =  plane_d/d;
158 >      t =  FP_D(peq)/d;
159        if(t < 0 )
160           hit = 0;
161        else
# Line 167 | Line 173 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
173   }
174  
175   int
176 < intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
176 > intersect_ray_plane(orig,dir,peq,pd,r)
177     FVECT orig,dir;
178 <   FVECT plane_n;
173 <   double plane_d;
178 >   FPEQ peq;
179     double *pd;
180     FVECT r;
181   {
# Line 185 | Line 190 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
190         line is  l = p1 + (p2-p1)t
191       */
192      /* Solve for t: */
193 <    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
193 >    t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/(DOT(FP_N(peq),dir));
194      if(t < 0)
195         hit = 0;
196      else
# Line 201 | Line 206 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
206  
207  
208   int
209 < intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
209 > intersect_ray_oplane(orig,dir,n,pd,r)
210 >   FVECT orig,dir;
211 >   FVECT n;
212 >   double *pd;
213 >   FVECT r;
214 > {
215 >  double t;
216 >  int hit;
217 >    /*
218 >      Plane is Ax + By + Cz +D = 0:
219 >      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
220 >    */
221 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
222 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
223 >       line is  l = p1 + (p2-p1)t
224 >     */
225 >    /* Solve for t: */
226 >    t =  -(DOT(n,orig))/(DOT(n,dir));
227 >    if(t < 0)
228 >       hit = 0;
229 >    else
230 >       hit = 1;
231 >
232 >  if(r)
233 >     VSUM(r,orig,dir,t);
234 >
235 >    if(pd)
236 >       *pd = t;
237 >  return(hit);
238 > }
239 >
240 >
241 > int
242 > intersect_edge_plane(e0,e1,peq,pd,r)
243     FVECT e0,e1;
244 <   FVECT plane_n;
207 <   double plane_d;
244 >   FPEQ peq;
245     double *pd;
246     FVECT r;
247   {
# Line 221 | Line 258 | intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
258       */
259      /* Solve for t: */
260    VSUB(d,e1,e0);
261 <  t =  -(DOT(plane_n,e0) + plane_d)/(DOT(plane_n,d));
261 >  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
262      if(t < 0)
263         hit = 0;
264      else
# Line 240 | Line 277 | point_in_cone(p,p0,p1,p2)
277   FVECT p;
278   FVECT p0,p1,p2;
279   {
243 –    FVECT n;
280      FVECT np,x_axis,y_axis;
281 <    double d1,d2,d;
281 >    double d1,d2;
282 >    FPEQ peq;
283      
284      /* Find the equation of the circle defined by the intersection
285         of the cone with the plane defined by p1,p2,p3- project p into
# Line 250 | Line 287 | FVECT p0,p1,p2;
287       */
288      
289      /* find the equation of the plane defined by p1-p3 */
290 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
290 >    tri_plane_equation(p0,p1,p2,&peq,FALSE);
291  
292      /* define a coordinate system on the plane: the x axis is in
293         the direction of np2-np1, and the y axis is calculated from
# Line 258 | Line 295 | FVECT p0,p1,p2;
295       */
296      /* Project p onto the plane */
297      /* NOTE: check this: does sideness check?*/
298 <    if(!intersect_vector_plane(p,n,d,NULL,np))
298 >    if(!intersect_vector_plane(p,peq,NULL,np))
299          return(FALSE);
300  
301      /* create coordinate system on  plane: p2-p1 defines the x_axis*/
302      VSUB(x_axis,p1,p0);
303      normalize(x_axis);
304      /* The y axis is  */
305 <    VCROSS(y_axis,n,x_axis);
305 >    VCROSS(y_axis,FP_N(peq),x_axis);
306      normalize(y_axis);
307  
308      VSUB(p1,p1,p0);
# Line 353 | Line 390 | int sides[3];
390  
391  
392  
393 <
393 >
394   int
395   point_in_stri(v0,v1,v2,p)
396   FVECT v0,v1,v2,p;
# Line 663 | Line 700 | FVECT orig,dir;
700   FVECT v0,v1,v2;
701   FVECT pt;
702   {
703 <  FVECT p0,p1,p2,p,n;
704 <  double pd;
703 >  FVECT p0,p1,p2,p;
704 >  FPEQ peq;
705    int type;
706  
707 <  point_on_sphere(p0,v0,orig);
708 <  point_on_sphere(p1,v1,orig);
709 <  point_on_sphere(p2,v2,orig);
710 <  
707 >  VSUB(p0,v0,orig);
708 >  VSUB(p1,v1,orig);
709 >  VSUB(p2,v2,orig);
710 >
711    if(point_in_stri(p0,p1,p2,dir))
712    {
713        /* Intersect the ray with the triangle plane */
714 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
715 <      return(intersect_ray_plane(orig,dir,n,pd,NULL,pt));
714 >      tri_plane_equation(v0,v1,v2,&peq,FALSE);
715 >      return(intersect_ray_plane(orig,dir,peq,NULL,pt));
716    }
717    return(FALSE);
718   }
# Line 736 | Line 773 | FVECT fnear[4],ffar[4];
773      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
774   }
775  
776 + int
777 + max_index(v,r)
778 + FVECT v;
779 + double *r;
780 + {
781 +  double p[3];
782 +  int i;
783  
784 +  p[0] = fabs(v[0]);
785 +  p[1] = fabs(v[1]);
786 +  p[2] = fabs(v[2]);
787 +  i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);  
788 +  if(r)
789 +    *r = p[i];
790 +  return(i);
791 + }
792  
793 + int
794 + closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
795 + FVECT p0,p1,p2,p;
796 + int p0id,p1id,p2id;
797 + {
798 +    double d,d1;
799 +    int i;
800 +    
801 +    d =  DIST_SQ(p,p0);
802 +    d1 = DIST_SQ(p,p1);
803 +    if(d < d1)
804 +    {
805 +      d1 = DIST_SQ(p,p2);
806 +      i = (d1 < d)?p2id:p0id;
807 +    }
808 +    else
809 +    {
810 +      d = DIST_SQ(p,p2);
811 +      i = (d < d1)? p2id:p1id;
812 +    }
813 +    return(i);
814 + }
815  
816 +
817   int
818   sedge_intersect(a0,a1,b0,b1)
819   FVECT a0,a1,b0,b1;
# Line 800 | Line 875 | double coord[3];
875    a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
876    coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
877    coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
878 <  coord[2]  = 1.0 - coord[0] - coord[1];
878 >  coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
879  
880   }
881  
# Line 1013 | Line 1088 | int child,next;
1088    }
1089   }
1090  
1091 < int
1092 < max_index(v)
1093 < FVECT v;
1091 >
1092 > baryi_parent(coord,i)
1093 > BCOORD coord[3];
1094 > int i;
1095   {
1020 –  double a,b,c;
1021 –  int i;
1096  
1097 <  a = fabs(v[0]);
1098 <  b = fabs(v[1]);
1099 <  c = fabs(v[2]);
1100 <  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
1101 <  return(i);
1097 >  switch(i) {
1098 >  case 0:
1099 >    /* update bary for child */
1100 >    coord[0] = (coord[0] >> 1) + MAXBCOORD2;
1101 >    coord[1] >>= 1;
1102 >    coord[2] >>= 1;
1103 >    break;
1104 >  case 1:
1105 >    coord[0] >>= 1;
1106 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD2;
1107 >    coord[2] >>= 1;
1108 >    break;
1109 >    
1110 >  case 2:
1111 >    coord[0] >>= 1;
1112 >    coord[1] >>= 1;
1113 >    coord[2] = (coord[2] >> 1) + MAXBCOORD2;
1114 >    break;
1115 >    
1116 >  case 3:
1117 >    coord[0] = MAXBCOORD2 - (coord[0] >> 1);
1118 >    coord[1] = MAXBCOORD2 - (coord[1] >> 1);
1119 >    coord[2] = MAXBCOORD2 - (coord[2] >> 1);
1120 >    break;
1121 > #ifdef DEBUG
1122 >  default:
1123 >    eputs("baryi_parent():Invalid child\n");
1124 >    break;
1125 > #endif
1126 >  }
1127   }
1128  
1129 + baryi_from_child(coord,child,next)
1130 + BCOORD coord[3];
1131 + int child,next;
1132 + {
1133 + #ifdef DEBUG
1134 +  if(child <0 || child > 3)
1135 +  {
1136 +    eputs("baryi_from_child():Invalid child\n");
1137 +    return;
1138 +  }
1139 +  if(next <0 || next > 3)
1140 +  {
1141 +    eputs("baryi_from_child():Invalid next\n");
1142 +    return;
1143 +  }
1144 + #endif
1145 +  if(next == child)
1146 +    return;
1147  
1148 +  switch(child){
1149 +  case 0:
1150 +      coord[0] = 0;
1151 +      coord[1] = MAXBCOORD - coord[1];
1152 +      coord[2] = MAXBCOORD - coord[2];
1153 +      break;
1154 +  case 1:
1155 +      coord[0] = MAXBCOORD - coord[0];
1156 +      coord[1] = 0;
1157 +      coord[2] = MAXBCOORD - coord[2];
1158 +      break;
1159 +  case 2:
1160 +      coord[0] = MAXBCOORD - coord[0];
1161 +      coord[1] = MAXBCOORD - coord[1];
1162 +      coord[2] = 0;
1163 +    break;
1164 +  case 3:
1165 +    switch(next){
1166 +    case 0:
1167 +      coord[0] = 0;
1168 +      coord[1] = MAXBCOORD - coord[1];
1169 +      coord[2] = MAXBCOORD - coord[2];
1170 +      break;
1171 +    case 1:
1172 +      coord[0] = MAXBCOORD - coord[0];
1173 +      coord[1] = 0;
1174 +      coord[2] = MAXBCOORD - coord[2];
1175 +      break;
1176 +    case 2:
1177 +      coord[0] = MAXBCOORD - coord[0];
1178 +      coord[1] = MAXBCOORD - coord[1];
1179 +      coord[2] = 0;
1180 +      break;
1181 +    }
1182 +    break;
1183 +  }
1184 + }
1185  
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 – */
1186   int
1187 < traceRay(orig,dir,v0,v1,v2,r)
1188 <  FVECT orig,dir;
1044 <  FVECT v0,v1,v2;
1045 <  FVECT r;
1187 > baryi_child(coord)
1188 > BCOORD coord[3];
1189   {
1047 –  FVECT n,p[3],d;
1048 –  double pt[3],r_eps;
1190    int i;
1050 –  int which;
1191  
1192 <  /* Find the plane equation for the triangle defined by the edge v0v1 and
1193 <   the view center*/
1194 <  VCROSS(n,v0,v1);
1195 <  /* Intersect the ray with this plane */
1196 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[0]),p[0]);
1197 <  if(i)
1198 <    which = 0;
1192 >  if(coord[0] > MAXBCOORD2)
1193 >  {
1194 >      /* update bary for child */
1195 >      coord[0] = (coord[0]<< 1) - MAXBCOORD;
1196 >      coord[1] <<= 1;
1197 >      coord[2] <<= 1;
1198 >      return(0);
1199 >  }
1200    else
1201 <    which = -1;
1201 >    if(coord[1] > MAXBCOORD2)
1202 >    {
1203 >      coord[0] <<= 1;
1204 >      coord[1] = (coord[1] << 1) - MAXBCOORD;
1205 >      coord[2] <<= 1;
1206 >      return(1);
1207 >    }
1208 >    else
1209 >      if(coord[2] > MAXBCOORD2)
1210 >      {
1211 >        coord[0] <<= 1;
1212 >        coord[1] <<= 1;
1213 >        coord[2] = (coord[2] << 1) - MAXBCOORD;
1214 >        return(2);
1215 >      }
1216 >      else
1217 >         {
1218 >           coord[0] = MAXBCOORD - (coord[0] << 1);
1219 >           coord[1] = MAXBCOORD - (coord[1] << 1);
1220 >           coord[2] = MAXBCOORD - (coord[2] << 1);
1221 >           return(3);
1222 >         }
1223 > }
1224  
1225 <  VCROSS(n,v1,v2);
1226 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[1]),p[1]);
1227 <  if(i)
1228 <    if((which==-1) || pt[1] < pt[0])
1229 <      which = 1;
1225 > int
1226 > baryi_nth_child(coord,i)
1227 > BCOORD coord[3];
1228 > int i;
1229 > {
1230  
1231 <  VCROSS(n,v2,v0);
1232 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[2]),p[2]);
1233 <  if(i)
1234 <    if((which==-1) || pt[2] < pt[which])
1235 <      which = 2;
1231 >  switch(i){
1232 >  case 0:
1233 >    /* update bary for child */
1234 >    coord[0] = (coord[0]<< 1) - MAXBCOORD;
1235 >    coord[1] <<= 1;
1236 >    coord[2] <<= 1;
1237 >    break;
1238 >  case 1:
1239 >    coord[0] <<= 1;
1240 >    coord[1] = (coord[1] << 1) - MAXBCOORD;
1241 >    coord[2] <<= 1;
1242 >    break;
1243 >  case 2:
1244 >    coord[0] <<= 1;
1245 >    coord[1] <<= 1;
1246 >    coord[2] = (coord[2] << 1) - MAXBCOORD;
1247 >    break;
1248 >  case 3:
1249 >    coord[0] = MAXBCOORD - (coord[0] << 1);
1250 >    coord[1] = MAXBCOORD - (coord[1] << 1);
1251 >    coord[2] = MAXBCOORD - (coord[2] << 1);
1252 >    break;
1253 >  }
1254 > }
1255  
1256 <  if(which != -1)
1256 >
1257 > baryi_children(coord,i,in_tri,rcoord,rin_tri)
1258 > BCOORD coord[3][3];
1259 > int i;
1260 > int in_tri[3];
1261 > BCOORD rcoord[3][3];
1262 > int rin_tri[3];
1263 > {
1264 >  int j;
1265 >
1266 >  for(j=0; j< 3; j++)
1267    {
1268 <      /* Return point that is K*eps along projection of the ray on
1269 <         the sphere to push intersection point p[which] into next cell
1270 <      */
1271 <      normalize(p[which]);
1272 <      /* Calculate the ray perpendicular to the intersection point: approx
1273 <       the direction moved along the sphere a distance of K*epsilon*/
1274 <      r_eps = -DOT(p[which],dir);
1275 <      VSUM(n,dir,p[which],r_eps);
1276 <     /* Calculate the length along ray p[which]-dir needed to move to
1277 <         cause a move along the sphere of k*epsilon
1278 <       */
1279 <       r_eps = DOT(n,dir);
1280 <      VSUM(r,p[which],dir,(20*FTINY)/r_eps);
1281 <      normalize(r);
1282 <      return(TRUE);
1268 >    if(!in_tri[j])
1269 >    {
1270 >      rin_tri[j]=0;
1271 >      continue;
1272 >    }
1273 >    
1274 >    if(i != 3)
1275 >    {
1276 >      if(coord[j][i] < MAXBCOORD2)
1277 >        {
1278 >          rin_tri[j] = 0;
1279 >          continue;
1280 >        }
1281 >    }
1282 >    else
1283 >      if( !(coord[j][0] <= MAXBCOORD2 && coord[j][1] <= MAXBCOORD2 &&
1284 >            coord[j][2] <= MAXBCOORD2))
1285 >        {
1286 >          rin_tri[j] = 0;
1287 >          continue;
1288 >        }
1289 >      
1290 >    rin_tri[j]=1;
1291 >    VCOPY(rcoord[j],coord[j]);
1292 >    baryi_nth_child(rcoord[j],i);
1293    }
1294 +
1295 + }
1296 +
1297 + convert_dtol(b,bi)
1298 + double b[3];
1299 + BCOORD bi[3];
1300 + {
1301 +  int i;
1302 +
1303 +  for(i=0; i < 2;i++)
1304 +  {
1305 +
1306 +    if(b[i] <= 0.0)
1307 +    {
1308 + #ifdef EXTRA_DEBUG
1309 +      if(b[i] < 0.0)
1310 +        printf("under %f\n",b[i]);
1311 + #endif
1312 +      bi[i]= 0;
1313 +    }
1314 +    else
1315 +      if(b[i] >= 1.0)
1316 +      {
1317 + #ifdef EXTRA_DEBUG
1318 +        if(b[i] > 1.0)
1319 +          printf("over %f\n",b[i]);
1320 + #endif
1321 +        bi[i]= MAXBCOORD;
1322 +      }
1323 +      else
1324 +        bi[i] = (BCOORD)(b[i]*MAXBCOORD);
1325 +  }
1326 +  bi[2] = bi[0] +  bi[1];
1327 +  if(bi[2] > MAXBCOORD)
1328 +  {
1329 + #ifdef EXTRA_DEBUG
1330 +      printf("sum over %f\n",b[0]+b[1]);
1331 + #endif
1332 +      bi[2] = 0;
1333 +      bi[1] = MAXBCOORD - bi[0];
1334 +  }
1335    else
1336 +    bi[2] = MAXBCOORD - bi[2];
1337 +
1338 + }
1339 +
1340 + /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1341 +   dir unbounded to [-MAXLONG,MAXLONG]
1342 + */
1343 + bary_dtol(b,db,bi,dbi,t,w)
1344 + double b[3],db[3][3];
1345 + BCOORD bi[3];
1346 + BDIR dbi[3][3];
1347 + TINT t[3];
1348 + int w;
1349 + {
1350 +  int i,id,j,k;
1351 +  double d;
1352 +
1353 +  convert_dtol(b,bi);
1354 +
1355 +  for(j=w; j< 3; j++)
1356    {
1357 <      /* Unable to find intersection: move along ray and try again */
1358 <      r_eps = -DOT(orig,dir);
1359 <      VSUM(n,dir,orig,r_eps);
1360 <      r_eps = DOT(n,dir);
1361 <      VSUM(r,orig,dir,(20*FTINY)/r_eps);
1362 <      normalize(r);
1363 < #ifdef DEBUG
1364 <      eputs("traceRay:Ray does not intersect triangle\n");
1365 < #endif
1366 <      return(FALSE);
1357 >    if(t[j] == HUGET)
1358 >    {
1359 >      max_index(db[j],&d);
1360 >      for(i=0; i< 3; i++)
1361 >        dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1362 >      break;
1363 >    }
1364 >    else
1365 >    {
1366 >      for(i=0; i< 3; i++)
1367 >          dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1368 >    }
1369    }
1370   }
1371  
1372  
1373 < int
1374 < closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
1375 < FVECT p0,p1,p2,p;
1376 < int p0id,p1id,p2id;
1373 > /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1374 >   dir unbounded to [-MAXLONG,MAXLONG]
1375 > */
1376 > bary_dtol_new(b,db,bi,boi,dbi,t)
1377 > double b[4][3],db[3][3];
1378 > BCOORD bi[3],boi[3][3];
1379 > BDIR dbi[3][3];
1380 > int t[3];
1381   {
1382 <    double d,d1;
1383 <    int i;
1382 >  int i,id,j,k;
1383 >  double d;
1384 >
1385 >  convert_dtol(b[3],bi);
1386 >
1387 >  for(j=0; j<3;j++)
1388 >  {
1389 >    if(t[j] != 1)
1390 >      continue;
1391 >    convert_dtol(b[j],boi[j]);
1392 >  }
1393 >  for(j=0; j< 3; j++)
1394 >  {
1395 >    k = (j+1)%3;
1396 >    if(t[k]==0)
1397 >      continue;
1398 >    if(t[k] == -1)
1399 >      {
1400 >        max_index(db[j],&d);
1401 >        for(i=0; i< 3; i++)
1402 >          dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1403 >        t[k] = 0;
1404 >      }
1405 >    else
1406 >      if(t[j] != 1)
1407 >        for(i=0; i< 3; i++)
1408 >          dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1409 >    else
1410 >      for(i=0; i< 3; i++)
1411 >        dbi[j][i] = boi[k][i] - boi[j][i];
1412      
1413 <    d =  DIST_SQ(p,p0);
1414 <    d1 = DIST_SQ(p,p1);
1415 <    if(d < d1)
1413 >  }
1414 > }
1415 >
1416 >
1417 > bary_dtolb(b,bi,in_tri)
1418 > double b[3][3];
1419 > BCOORD bi[3][3];
1420 > int in_tri[3];
1421 > {
1422 >  int i,j;
1423 >
1424 >  for(j=0; j<3;j++)
1425 >  {
1426 >    if(!in_tri[j])
1427 >      continue;
1428 >    for(i=0; i < 2;i++)
1429      {
1430 <      d1 = DIST_SQ(p,p2);
1431 <      i = (d1 < d)?p2id:p0id;
1430 >    if(b[j][i] <= 0.0)
1431 >    {
1432 >      bi[j][i]= 0;
1433      }
1434      else
1435 <    {
1436 <      d = DIST_SQ(p,p2);
1437 <      i = (d < d1)? p2id:p1id;
1435 >      if(b[j][i] >= 1.0)
1436 >      {
1437 >        bi[j][i]= MAXBCOORD;
1438 >      }
1439 >      else
1440 >        bi[j][i] = (BCOORD)(b[j][i]*MAXBCOORD);
1441      }
1442 <    return(i);
1442 >    bi[j][2] = MAXBCOORD - bi[j][0] - bi[j][1];
1443 >    if(bi[j][2] < 0)
1444 >      {
1445 >        bi[j][2] = 0;
1446 >        bi[j][1] = MAXBCOORD - bi[j][0];
1447 >      }
1448 >  }
1449   }
1450 +
1451 +
1452 +
1453 +
1454  
1455  
1456  

Diff Legend

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