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.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    VSUB(p0,v0,orig);
# Line 674 | Line 711 | FVECT pt;
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 1185 | Line 1222 | BCOORD coord[3];
1222           }
1223   }
1224  
1225 < /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1226 <   dir unbounded to [-MAXLONG,MAXLONG]
1227 < */
1228 < bary_dtol(b,db,bi,dbi,t)
1229 < double b[3],db[3][3];
1225 > int
1226 > baryi_nth_child(coord,i)
1227 > BCOORD coord[3];
1228 > int i;
1229 > {
1230 >
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 >
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 >    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];
1194 – BDIR dbi[3][3];
1195 – TINT t[3];
1300   {
1301 <  int i,id,j;
1198 <  double d;
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] = MAXBCOORD - bi[0] - bi[1];
1327 <
1216 <  if(bi[2] < 0)
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 +    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 + /* 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 +  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 <    if(t[j]==0)
1395 >    k = (j+1)%3;
1396 >    if(t[k]==0)
1397        continue;
1398 <    if(t[j] == HUGET)
1399 <       max_index(db[j],&d);
1400 <    for(i=0; i< 3; i++)
1401 <       if(t[j] != HUGET)
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 <          dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1409 >    else
1410 >      for(i=0; i< 3; i++)
1411 >        dbi[j][i] = boi[k][i] - boi[j][i];
1412 >    
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 +    if(b[j][i] <= 0.0)
1431 +    {
1432 +      bi[j][i]= 0;
1433 +    }
1434 +    else
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 +    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  

Diff Legend

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