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.7 by gwlarson, Tue Oct 6 18:16:53 1998 UTC vs.
Revision 3.12 by gwlarson, Thu Jun 10 15:22:22 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 cp,cp01,cp12,v10,v02;
40 <    double dp;
41 <    /*
42 <      VSUB(v10,v1,v0);
43 <      VSUB(v02,v0,v2);
44 <      VCROSS(cp,v10,v02);
45 <   */
46 <      /* test sign of (v0Xv1)X(v1Xv2). v1 */
47 <    VCROSS(cp01,v0,v1);
48 <    VCROSS(cp12,v1,v2);
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 >    /* 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(ZERO(dp) || dp < 0.0)
59 <       return(FALSE);
60 <    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)
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 70 | 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 78 | Line 100 | int norm;
100  
101    if(!norm)
102       return(0);
81 –
103    
104    mag = normalize(n);
105  
# Line 87 | 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 96 | Line 117 | tri_plane_equation(v0,v1,v2,peqptr,norm)
117      FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
118   }
119  
99 – /* From quad_edge-code */
100 – int
101 – point_in_circle_thru_origin(p,p0,p1)
102 – FVECT p;
103 – FVECT p0,p1;
104 – {
120  
106 –    double dp0,dp1;
107 –    double dp,det;
108 –    
109 –    dp0 = DOT_VEC2(p0,p0);
110 –    dp1 = DOT_VEC2(p1,p1);
111 –    dp  = DOT_VEC2(p,p);    
112 –
113 –    det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
114 –    
115 –    return (det > 0);
116 – }
117 –
118 –
119 –
120 – point_on_sphere(ps,p,c)
121 – FVECT ps,p,c;
122 – {
123 –    VSUB(ps,p,c);    
124 –    normalize(ps);
125 – }
126 –
127 –
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 179 | Line 172 | intersect_ray_plane(orig,dir,peq,pd,r)
172     double *pd;
173     FVECT r;
174   {
175 <  double t;
175 >  double t,d;
176    int hit;
177      /*
178        Plane is Ax + By + Cz +D = 0:
# Line 190 | Line 183 | intersect_ray_plane(orig,dir,peq,pd,r)
183         line is  l = p1 + (p2-p1)t
184       */
185      /* Solve for t: */
186 <    t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/(DOT(FP_N(peq),dir));
187 <    if(t < 0)
186 >  d = DOT(FP_N(peq),dir);
187 >  if(ZERO(d))
188 >     return(0);
189 >  t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
190 >
191 >  if(t < 0)
192         hit = 0;
193      else
194         hit = 1;
# Line 212 | Line 209 | intersect_ray_oplane(orig,dir,n,pd,r)
209     double *pd;
210     FVECT r;
211   {
212 <  double t;
212 >  double t,d;
213    int hit;
214      /*
215        Plane is Ax + By + Cz +D = 0:
# Line 223 | Line 220 | intersect_ray_oplane(orig,dir,n,pd,r)
220         line is  l = p1 + (p2-p1)t
221       */
222      /* Solve for t: */
223 <    t =  -(DOT(n,orig))/(DOT(n,dir));
223 >    d= DOT(n,dir);
224 >    if(ZERO(d))
225 >       return(0);
226 >    t =  -(DOT(n,orig))/d;
227      if(t < 0)
228         hit = 0;
229      else
# Line 237 | 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 271 | Line 289 | intersect_edge_plane(e0,e1,peq,pd,r)
289    return(hit);
290   }
291  
274 –
292   int
276 – point_in_cone(p,p0,p1,p2)
277 – FVECT p;
278 – FVECT p0,p1,p2;
279 – {
280 –    FVECT np,x_axis,y_axis;
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
286 –       that plane and do an in-circle test in the plane
287 –     */
288 –    
289 –    /* find the equation of the plane defined by p1-p3 */
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
294 –       n cross x-axis
295 –     */
296 –    /* Project p onto the plane */
297 –    /* NOTE: check this: does sideness check?*/
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,FP_N(peq),x_axis);
306 –    normalize(y_axis);
307 –
308 –    VSUB(p1,p1,p0);
309 –    VSUB(p2,p2,p0);
310 –    VSUB(np,np,p0);
311 –    
312 –    p1[0] = VLEN(p1);
313 –    p1[1] = 0;
314 –
315 –    d1 = DOT(p2,x_axis);
316 –    d2 = DOT(p2,y_axis);
317 –    p2[0] = d1;
318 –    p2[1] = d2;
319 –    
320 –    d1 = DOT(np,x_axis);
321 –    d2 = DOT(np,y_axis);
322 –    np[0] = d1;
323 –    np[1] = d2;
324 –
325 –    /* perform the in-circle test in the new coordinate system */
326 –    return(point_in_circle_thru_origin(np,p1,p2));
327 – }
328 –
329 – int
293   point_set_in_stri(v0,v1,v2,p,n,nset,sides)
294   FVECT v0,v1,v2,p;
295   FVECT n[3];
# Line 338 | Line 301 | int sides[3];
301      /* Find the normal to the triangle ORIGIN:v0:v1 */
302      if(!NTH_BIT(*nset,0))
303      {
304 <        VCROSS(n[0],v1,v0);
304 >        VCROSS(n[0],v0,v1);
305          SET_NTH_BIT(*nset,0);
306      }
307      /* Test the point for sidedness */
# Line 356 | Line 319 | int sides[3];
319      /* Test next edge */
320      if(!NTH_BIT(*nset,1))
321      {
322 <        VCROSS(n[1],v2,v1);
322 >        VCROSS(n[1],v1,v2);
323          SET_NTH_BIT(*nset,1);
324      }
325      /* Test the point for sidedness */
# Line 372 | Line 335 | int sides[3];
335      /* Test next edge */
336      if(!NTH_BIT(*nset,2))
337      {
338 <        VCROSS(n[2],v0,v2);
338 >        VCROSS(n[2],v2,v0);
339          SET_NTH_BIT(*nset,2);
340      }
341      /* Test the point for sidedness */
# Line 390 | Line 353 | int sides[3];
353  
354  
355  
393 –
394 – int
395 – point_in_stri(v0,v1,v2,p)
396 – FVECT v0,v1,v2,p;
397 – {
398 –    double d;
399 –    FVECT n;  
400 –
401 –    VCROSS(n,v1,v0);
402 –    /* Test the point for sidedness */
403 –    d  = DOT(n,p);
404 –    if(d > 0.0)
405 –      return(FALSE);
406 –
407 –    /* Test next edge */
408 –    VCROSS(n,v2,v1);
409 –    /* Test the point for sidedness */
410 –    d  = DOT(n,p);
411 –    if(d > 0.0)
412 –       return(FALSE);
413 –
414 –    /* Test next edge */
415 –    VCROSS(n,v0,v2);
416 –    /* Test the point for sidedness */
417 –    d  = DOT(n,p);
418 –    if(d > 0.0)
419 –       return(FALSE);
420 –    /* Must be interior to the pyramid */
421 –    return(GT_INTERIOR);
422 – }
423 –
424 – int
425 – vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
426 – FVECT t0,t1,t2,p0,p1,p2;
427 – int *nset;
428 – FVECT n[3];
429 – FVECT avg;
430 – int pt_sides[3][3];
431 –
432 – {
433 –    int below_plane[3],test;
434 –
435 –    SUM_3VEC3(avg,t0,t1,t2);
436 –    *nset = 0;
437 –    /* Test vertex v[i] against triangle j*/
438 –    /* Check if v[i] lies below plane defined by avg of 3 vectors
439 –       defining triangle
440 –       */
441 –
442 –    /* test point 0 */
443 –    if(DOT(avg,p0) < 0.0)
444 –      below_plane[0] = 1;
445 –    else
446 –      below_plane[0] = 0;
447 –    /* Test if b[i] lies in or on triangle a */
448 –    test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
449 –    /* If pts[i] is interior: done */
450 –    if(!below_plane[0])
451 –      {
452 –        if(test == GT_INTERIOR)
453 –          return(TRUE);
454 –      }
455 –    /* Now test point 1*/
456 –
457 –    if(DOT(avg,p1) < 0.0)
458 –      below_plane[1] = 1;
459 –    else
460 –      below_plane[1]=0;
461 –    /* Test if b[i] lies in or on triangle a */
462 –    test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
463 –    /* If pts[i] is interior: done */
464 –    if(!below_plane[1])
465 –    {
466 –      if(test == GT_INTERIOR)
467 –        return(TRUE);
468 –    }
469 –    
470 –    /* Now test point 2 */
471 –    if(DOT(avg,p2) < 0.0)
472 –      below_plane[2] = 1;
473 –    else
474 –      below_plane[2] = 0;
475 –        /* Test if b[i] lies in or on triangle a */
476 –    test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
477 –
478 –    /* If pts[i] is interior: done */
479 –    if(!below_plane[2])
480 –      {
481 –        if(test == GT_INTERIOR)
482 –          return(TRUE);
483 –      }
484 –
485 –    /* If all three points below separating plane: trivial reject */
486 –    if(below_plane[0] && below_plane[1] && below_plane[2])
487 –       return(FALSE);
488 –    /* Now check vertices in a against triangle b */
489 –    return(FALSE);
490 – }
491 –
492 –
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 506 | Line 369 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
369      if(sides[0][0] == GT_INVALID)
370      {
371        if(!NTH_BIT(nset,0))
372 <        VCROSS(n[0],t1,t0);
372 >        VCROSS(n[0],t0,t1);
373        /* Test the point for sidedness */
374        d  = DOT(n[0],p0);
375        if(d >= 0.0)
# Line 519 | Line 382 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
382      if(sides[0][1] == GT_INVALID)
383      {
384        if(!NTH_BIT(nset,1))
385 <        VCROSS(n[1],t2,t1);
385 >        VCROSS(n[1],t1,t2);
386          /* Test the point for sidedness */
387          d  = DOT(n[1],p0);
388          if(d >= 0.0)
# Line 532 | Line 395 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
395      if(sides[0][2] == GT_INVALID)
396      {
397        if(!NTH_BIT(nset,2))
398 <        VCROSS(n[2],t0,t2);
398 >        VCROSS(n[2],t2,t0);
399        /* Test the point for sidedness */
400        d  = DOT(n[2],p0);
401        if(d >= 0.0)
# Line 548 | Line 411 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
411      if(sides[1][0] == GT_INVALID)
412      {
413        if(!NTH_BIT(nset,0))
414 <        VCROSS(n[0],t1,t0);
414 >        VCROSS(n[0],t0,t1);
415        /* Test the point for sidedness */
416        d  = DOT(n[0],p1);
417        if(d >= 0.0)
# Line 562 | Line 425 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
425      if(sides[1][1] == GT_INVALID)
426      {
427        if(!NTH_BIT(nset,1))
428 <        VCROSS(n[1],t2,t1);
428 >        VCROSS(n[1],t1,t2);
429        /* Test the point for sidedness */
430        d  = DOT(n[1],p1);
431        if(d >= 0.0)
# Line 576 | Line 439 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
439      if(sides[1][2] == GT_INVALID)
440      {
441        if(!NTH_BIT(nset,2))
442 <        VCROSS(n[2],t0,t2);
442 >        VCROSS(n[2],t2,t0);
443        /* Test the point for sidedness */
444        d  = DOT(n[2],p1);
445        if(d >= 0.0)
# Line 592 | Line 455 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
455      if(sides[2][0] == GT_INVALID)
456      {
457        if(!NTH_BIT(nset,0))
458 <        VCROSS(n[0],t1,t0);
458 >        VCROSS(n[0],t0,t1);
459        /* Test the point for sidedness */
460        d  = DOT(n[0],p2);
461        if(d >= 0.0)
# Line 605 | Line 468 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
468      if(sides[2][1] == GT_INVALID)
469      {
470        if(!NTH_BIT(nset,1))
471 <        VCROSS(n[1],t2,t1);
471 >        VCROSS(n[1],t1,t2);
472        /* Test the point for sidedness */
473        d  = DOT(n[1],p2);
474        if(d >= 0.0)
# Line 618 | Line 481 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
481      if(sides[2][2] == GT_INVALID)
482      {
483        if(!NTH_BIT(nset,2))
484 <        VCROSS(n[2],t0,t2);
484 >        VCROSS(n[2],t2,t0);
485        /* Test the point for sidedness */
486        d  = DOT(n[2],p2);
487        if(d >= 0.0)
# Line 629 | 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;
639 <  int tests[2][3];
640 <  FVECT n[2][3],p,avg[2];
641 <
642 <  /* Test the vertices of triangle a against the pyramid formed by triangle
643 <     b and the origin. If any vertex of a is interior to triangle b, or
644 <     if all 3 vertices of a are ON the edges of b,return TRUE. Remember
645 <     the results of the edge normal and sidedness tests for later.
646 <   */
647 < if(vertices_in_stri(a1,a2,a3,b1,b2,b3,&(n_set[0]),n[0],avg[0],sides[1]))
648 <     return(TRUE);
649 <  
650 < if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
651 <     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 */
666 <    if(tests[1][j] & tests[1][jnext])
667 <      continue;
668 <    for(i=0;i<3;i++)
669 <    {
670 <      inext = (i+1)%3;
671 <      /* IF edge a doesnt cross any great circles of b, punt */
672 <      if(tests[0][i] & tests[0][inext])
673 <        continue;
674 <      /* Now find the great circles that cross and test */
675 <      if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
676 <          && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
677 <      {
678 <        VCROSS(p,n[0][i],n[1][j]);
679 <                    
680 <        /* If zero cp= done */
681 <        if(ZERO_VEC3(p))
682 <          continue;
683 <        /* check above both planes */
684 <        if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
685 <          {
686 <            NEGATE_VEC3(p);
687 <            if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
688 <              continue;
689 <          }
690 <        return(TRUE);
691 <      }
692 <    }
693 <  }
694 <  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 813 | Line 652 | int p0id,p1id,p2id;
652      return(i);
653   }
654  
816 –
817 – int
818 – sedge_intersect(a0,a1,b0,b1)
819 – FVECT a0,a1,b0,b1;
820 – {
821 –    FVECT na,nb,avga,avgb,p;
822 –    double d;
823 –    int sb0,sb1,sa0,sa1;
824 –
825 –    /* First test if edge b straddles great circle of a */
826 –    VCROSS(na,a0,a1);
827 –    d = DOT(na,b0);
828 –    sb0 = ZERO(d)?0:(d<0)? -1: 1;
829 –    d = DOT(na,b1);
830 –    sb1 = ZERO(d)?0:(d<0)? -1: 1;
831 –    /* edge b entirely on one side of great circle a: edges cannot intersect*/
832 –    if(sb0*sb1 > 0)
833 –       return(FALSE);
834 –    /* test if edge a straddles great circle of b */
835 –    VCROSS(nb,b0,b1);
836 –    d = DOT(nb,a0);
837 –    sa0 = ZERO(d)?0:(d<0)? -1: 1;
838 –    d = DOT(nb,a1);
839 –    sa1 = ZERO(d)?0:(d<0)? -1: 1;
840 –    /* edge a entirely on one side of great circle b: edges cannot intersect*/
841 –    if(sa0*sa1 > 0)
842 –       return(FALSE);
843 –
844 –    /* Find one of intersection points of the great circles */
845 –    VCROSS(p,na,nb);
846 –    /* If they lie on same great circle: call an intersection */
847 –    if(ZERO_VEC3(p))
848 –       return(TRUE);
849 –
850 –    VADD(avga,a0,a1);
851 –    VADD(avgb,b0,b1);
852 –    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
853 –    {
854 –      NEGATE_VEC3(p);
855 –      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
856 –        return(FALSE);
857 –    }
858 –    if((!sb0 || !sb1) && (!sa0 || !sa1))
859 –      return(FALSE);
860 –    return(TRUE);
861 – }
862 –
863 –
864 –
655   /* Find the normalized barycentric coordinates of p relative to
656   * triangle v0,v1,v2. Return result in coord
657   */
# Line 879 | Line 669 | double coord[3];
669  
670   }
671  
882 – bary_ith_child(coord,i)
883 – double coord[3];
884 – int i;
885 – {
886 –
887 –  switch(i){
888 –  case 0:
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 –      break;
894 –  case 1:
895 –    coord[0] *= 2.0;
896 –    coord[1] = 2.0*coord[1]- 1.0;
897 –    coord[2] *= 2.0;
898 –    break;
899 –  case 2:
900 –    coord[0] *= 2.0;
901 –    coord[1] *= 2.0;
902 –    coord[2] = 2.0*coord[2]- 1.0;
903 –    break;
904 –  case 3:
905 –    coord[0] = 1.0 - 2.0*coord[0];
906 –    coord[1] = 1.0 - 2.0*coord[1];
907 –    coord[2] = 1.0 - 2.0*coord[2];
908 –    break;
909 – #ifdef DEBUG
910 –  default:
911 –    eputs("bary_ith_child():Invalid child\n");
912 –    break;
913 – #endif
914 –  }
915 – }
672  
673  
918 – int
919 – bary_child(coord)
920 – double coord[3];
921 – {
922 –  int i;
674  
924 –  if(coord[0] > 0.5)
925 –  {
926 –      /* update bary for child */
927 –      coord[0] = 2.0*coord[0]- 1.0;
928 –      coord[1] *= 2.0;
929 –      coord[2] *= 2.0;
930 –      return(0);
931 –  }
932 –  else
933 –    if(coord[1] > 0.5)
934 –    {
935 –      coord[0] *= 2.0;
936 –      coord[1] = 2.0*coord[1]- 1.0;
937 –      coord[2] *= 2.0;
938 –      return(1);
939 –    }
940 –    else
941 –      if(coord[2] > 0.5)
942 –      {
943 –        coord[0] *= 2.0;
944 –        coord[1] *= 2.0;
945 –        coord[2] = 2.0*coord[2]- 1.0;
946 –        return(2);
947 –      }
948 –      else
949 –         {
950 –           coord[0] = 1.0 - 2.0*coord[0];
951 –           coord[1] = 1.0 - 2.0*coord[1];
952 –           coord[2] = 1.0 - 2.0*coord[2];
953 –           return(3);
954 –         }
955 – }
956 –
957 – /* Coord was the ith child of the parent: set the coordinate
958 –   relative to the parent
959 – */
675   bary_parent(coord,i)
961 – double coord[3];
962 – int i;
963 – {
964 –
965 –  switch(i) {
966 –  case 0:
967 –    /* update bary for child */
968 –    coord[0] = coord[0]*0.5 + 0.5;
969 –    coord[1] *= 0.5;
970 –    coord[2] *= 0.5;
971 –    break;
972 –  case 1:
973 –    coord[0] *= 0.5;
974 –    coord[1]  = coord[1]*0.5 + 0.5;
975 –    coord[2] *= 0.5;
976 –    break;
977 –    
978 –  case 2:
979 –    coord[0] *= 0.5;
980 –    coord[1] *= 0.5;
981 –    coord[2] = coord[2]*0.5 + 0.5;
982 –    break;
983 –    
984 –  case 3:
985 –    coord[0] = 0.5 - 0.5*coord[0];
986 –    coord[1] = 0.5 - 0.5*coord[1];
987 –    coord[2] = 0.5 - 0.5*coord[2];
988 –    break;
989 – #ifdef DEBUG
990 –  default:
991 –    eputs("bary_parent():Invalid child\n");
992 –    break;
993 – #endif
994 –  }
995 – }
996 –
997 – bary_from_child(coord,child,next)
998 – double coord[3];
999 – int child,next;
1000 – {
1001 – #ifdef DEBUG
1002 –  if(child <0 || child > 3)
1003 –  {
1004 –    eputs("bary_from_child():Invalid child\n");
1005 –    return;
1006 –  }
1007 –  if(next <0 || next > 3)
1008 –  {
1009 –    eputs("bary_from_child():Invalid next\n");
1010 –    return;
1011 –  }
1012 – #endif
1013 –  if(next == child)
1014 –    return;
1015 –
1016 –  switch(child){
1017 –  case 0:
1018 –    switch(next){
1019 –    case 1:
1020 –      coord[0] += 1.0;
1021 –      coord[1] -= 1.0;
1022 –      break;
1023 –    case 2:
1024 –      coord[0] += 1.0;
1025 –      coord[2] -= 1.0;
1026 –      break;
1027 –    case 3:
1028 –      coord[0] *= -1.0;
1029 –      coord[1] = 1 - coord[1];
1030 –      coord[2] = 1 - coord[2];
1031 –      break;
1032 –
1033 –    }
1034 –    break;
1035 –  case 1:
1036 –    switch(next){
1037 –    case 0:
1038 –      coord[0] -= 1.0;
1039 –      coord[1] += 1.0;
1040 –      break;
1041 –    case 2:
1042 –      coord[1] += 1.0;
1043 –      coord[2] -= 1.0;
1044 –      break;
1045 –    case 3:
1046 –      coord[0] = 1 - coord[0];
1047 –      coord[1] *= -1.0;
1048 –      coord[2] = 1 - coord[2];
1049 –      break;
1050 –    }
1051 –    break;
1052 –  case 2:
1053 –    switch(next){
1054 –    case 0:
1055 –      coord[0] -= 1.0;
1056 –      coord[2] += 1.0;
1057 –      break;
1058 –    case 1:
1059 –      coord[1] -= 1.0;
1060 –      coord[2] += 1.0;
1061 –      break;
1062 –    case 3:
1063 –      coord[0] = 1 - coord[0];
1064 –      coord[1] = 1 - coord[1];
1065 –      coord[2] *= -1.0;
1066 –      break;
1067 –    }
1068 –    break;
1069 –  case 3:
1070 –    switch(next){
1071 –    case 0:
1072 –      coord[0] *= -1.0;
1073 –      coord[1] = 1 - coord[1];
1074 –      coord[2] = 1 - coord[2];
1075 –      break;
1076 –    case 1:
1077 –      coord[0] = 1 - coord[0];
1078 –      coord[1] *= -1.0;
1079 –      coord[2] = 1 - coord[2];
1080 –      break;
1081 –    case 2:
1082 –      coord[0] = 1 - coord[0];
1083 –      coord[1] = 1 - coord[1];
1084 –      coord[2] *= -1.0;
1085 –      break;
1086 –    }
1087 –    break;
1088 –  }
1089 – }
1090 –
1091 –
1092 – baryi_parent(coord,i)
676   BCOORD coord[3];
677   int i;
678   {
1096 –
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 1148 | 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 1184 | 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 1231 | 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   }
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];
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 –    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 –    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 –  }
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 –
1455 –
1456 –
1457 –
837  
838  
839  

Diff Legend

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