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.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 +    double dp,dp1;
40 +    int test,test1;
41 +    FVECT nv0,nv1,nv2;
42      FVECT cp01,cp12,cp;
43 <    
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 <    if(DOT(cp,v1) < 0.0)
57 <       return(FALSE);
58 <    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 63 | 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]);
67 <
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 72 | Line 100 | int norm;
100  
101    if(!norm)
102       return(0);
75 –
103    
104    mag = normalize(n);
105  
# Line 80 | Line 107 | int norm;
107   }
108  
109  
110 < tri_plane_equation(v0,v1,v2,n,nd,norm)
111 <   FVECT v0,v1,v2,n;
112 <   double *nd;
110 > tri_plane_equation(v0,v1,v2,peqptr,norm)
111 >   FVECT v0,v1,v2;
112 >   FPEQ *peqptr;
113     int norm;
114   {
115 <    tri_normal(v0,v1,v2,n,norm);
115 >    tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
116  
117 <    *nd = -(DOT(n,v0));
117 >    FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
118   }
119  
93 – /* From quad_edge-code */
94 – int
95 – point_in_circle_thru_origin(p,p0,p1)
96 – FVECT p;
97 – FVECT p0,p1;
98 – {
120  
100 –    double dp0,dp1;
101 –    double dp,det;
102 –    
103 –    dp0 = DOT_VEC2(p0,p0);
104 –    dp1 = DOT_VEC2(p1,p1);
105 –    dp  = DOT_VEC2(p,p);    
106 –
107 –    det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
108 –    
109 –    return (det > 0);
110 – }
111 –
112 –
113 –
114 – point_on_sphere(ps,p,c)
115 – FVECT ps,p,c;
116 – {
117 –    VSUB(ps,p,c);    
118 –    normalize(ps);
119 – }
120 –
121 –
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
124    * t, if parallel, returns t=FHUGE
125    */
126   int
127 < intersect_vector_plane(v,plane_n,plane_d,tptr,r)
128 <   FVECT v,plane_n;
129 <   double plane_d;
127 > intersect_vector_plane(v,peq,tptr,r)
128 >   FVECT v;
129 >   FPEQ peq;
130     double *tptr;
131     FVECT r;
132   {
# Line 141 | Line 140 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
140      /* line is  l = p1 + (p2-p1)t, p1=origin */
141  
142      /* Solve for t: */
143 <  d = -(DOT(plane_n,v));
143 >  d = -(DOT(FP_N(peq),v));
144    if(ZERO(d))
145    {
146        t = FHUGE;
# Line 149 | Line 148 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
148    }
149    else
150    {
151 <      t =  plane_d/d;
151 >      t =  FP_D(peq)/d;
152        if(t < 0 )
153           hit = 0;
154        else
# Line 167 | Line 166 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
166   }
167  
168   int
169 < intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
169 > intersect_ray_plane(orig,dir,peq,pd,r)
170     FVECT orig,dir;
171 <   FVECT plane_n;
173 <   double plane_d;
171 >   FPEQ peq;
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 185 | Line 183 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
183         line is  l = p1 + (p2-p1)t
184       */
185      /* Solve for t: */
186 <    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,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 201 | Line 203 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
203  
204  
205   int
206 < intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
207 <   FVECT e0,e1;
208 <   FVECT plane_n;
207 <   double plane_d;
206 > intersect_ray_oplane(orig,dir,n,pd,r)
207 >   FVECT orig,dir;
208 >   FVECT n;
209     double *pd;
210     FVECT r;
211   {
212 <  double t;
212 >  double t,d;
213    int hit;
214 <  FVECT d;
214 <  /*
214 >    /*
215        Plane is Ax + By + Cz +D = 0:
216        plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
217      */
# Line 220 | Line 220 | intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
220         line is  l = p1 + (p2-p1)t
221       */
222      /* Solve for t: */
223 <  VSUB(d,e1,e0);
224 <  t =  -(DOT(plane_n,e0) + plane_d)/(DOT(plane_n,d));
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
230         hit = 1;
231  
232 <  VSUM(r,e0,d,t);
232 >  if(r)
233 >     VSUM(r,orig,dir,t);
234  
235      if(pd)
236         *pd = t;
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 < point_in_cone(p,p0,p1,p2)
261 < FVECT p;
262 < FVECT p0,p1,p2;
260 > intersect_edge_plane(e0,e1,peq,pd,r)
261 >   FVECT e0,e1;
262 >   FPEQ peq;
263 >   double *pd;
264 >   FVECT r;
265   {
266 <    FVECT n;
267 <    FVECT np,x_axis,y_axis;
268 <    double d1,d2,d;
269 <    
270 <    /* Find the equation of the circle defined by the intersection
271 <       of the cone with the plane defined by p1,p2,p3- project p into
272 <       that plane and do an in-circle test in the plane
266 >  double t;
267 >  int hit;
268 >  FVECT d;
269 >  /*
270 >      Plane is Ax + By + Cz +D = 0:
271 >      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
272 >    */
273 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
274 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
275 >       line is  l = p1 + (p2-p1)t
276       */
277 <    
278 <    /* find the equation of the plane defined by p1-p3 */
279 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
277 >    /* Solve for t: */
278 >  VSUB(d,e1,e0);
279 >  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
280 >    if(t < 0)
281 >       hit = 0;
282 >    else
283 >       hit = 1;
284  
285 <    /* define a coordinate system on the plane: the x axis is in
256 <       the direction of np2-np1, and the y axis is calculated from
257 <       n cross x-axis
258 <     */
259 <    /* Project p onto the plane */
260 <    /* NOTE: check this: does sideness check?*/
261 <    if(!intersect_vector_plane(p,n,d,NULL,np))
262 <        return(FALSE);
285 >  VSUM(r,e0,d,t);
286  
287 <    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
288 <    VSUB(x_axis,p1,p0);
289 <    normalize(x_axis);
267 <    /* The y axis is  */
268 <    VCROSS(y_axis,n,x_axis);
269 <    normalize(y_axis);
270 <
271 <    VSUB(p1,p1,p0);
272 <    VSUB(p2,p2,p0);
273 <    VSUB(np,np,p0);
274 <    
275 <    p1[0] = VLEN(p1);
276 <    p1[1] = 0;
277 <
278 <    d1 = DOT(p2,x_axis);
279 <    d2 = DOT(p2,y_axis);
280 <    p2[0] = d1;
281 <    p2[1] = d2;
282 <    
283 <    d1 = DOT(np,x_axis);
284 <    d2 = DOT(np,y_axis);
285 <    np[0] = d1;
286 <    np[1] = d2;
287 <
288 <    /* perform the in-circle test in the new coordinate system */
289 <    return(point_in_circle_thru_origin(np,p1,p2));
287 >    if(pd)
288 >       *pd = t;
289 >  return(hit);
290   }
291  
292   int
# Line 301 | 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 319 | 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 335 | 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 353 | Line 353 | int sides[3];
353  
354  
355  
356 –
357 – int
358 – point_in_stri(v0,v1,v2,p)
359 – FVECT v0,v1,v2,p;
360 – {
361 –    double d;
362 –    FVECT n;  
363 –
364 –    VCROSS(n,v1,v0);
365 –    /* Test the point for sidedness */
366 –    d  = DOT(n,p);
367 –    if(d > 0.0)
368 –      return(FALSE);
369 –
370 –    /* Test next edge */
371 –    VCROSS(n,v2,v1);
372 –    /* Test the point for sidedness */
373 –    d  = DOT(n,p);
374 –    if(d > 0.0)
375 –       return(FALSE);
376 –
377 –    /* Test next edge */
378 –    VCROSS(n,v0,v2);
379 –    /* Test the point for sidedness */
380 –    d  = DOT(n,p);
381 –    if(d > 0.0)
382 –       return(FALSE);
383 –    /* Must be interior to the pyramid */
384 –    return(GT_INTERIOR);
385 – }
386 –
387 – int
388 – vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
389 – FVECT t0,t1,t2,p0,p1,p2;
390 – int *nset;
391 – FVECT n[3];
392 – FVECT avg;
393 – int pt_sides[3][3];
394 –
395 – {
396 –    int below_plane[3],test;
397 –
398 –    SUM_3VEC3(avg,t0,t1,t2);
399 –    *nset = 0;
400 –    /* Test vertex v[i] against triangle j*/
401 –    /* Check if v[i] lies below plane defined by avg of 3 vectors
402 –       defining triangle
403 –       */
404 –
405 –    /* test point 0 */
406 –    if(DOT(avg,p0) < 0.0)
407 –      below_plane[0] = 1;
408 –    else
409 –      below_plane[0] = 0;
410 –    /* Test if b[i] lies in or on triangle a */
411 –    test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
412 –    /* If pts[i] is interior: done */
413 –    if(!below_plane[0])
414 –      {
415 –        if(test == GT_INTERIOR)
416 –          return(TRUE);
417 –      }
418 –    /* Now test point 1*/
419 –
420 –    if(DOT(avg,p1) < 0.0)
421 –      below_plane[1] = 1;
422 –    else
423 –      below_plane[1]=0;
424 –    /* Test if b[i] lies in or on triangle a */
425 –    test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
426 –    /* If pts[i] is interior: done */
427 –    if(!below_plane[1])
428 –    {
429 –      if(test == GT_INTERIOR)
430 –        return(TRUE);
431 –    }
432 –    
433 –    /* Now test point 2 */
434 –    if(DOT(avg,p2) < 0.0)
435 –      below_plane[2] = 1;
436 –    else
437 –      below_plane[2] = 0;
438 –        /* Test if b[i] lies in or on triangle a */
439 –    test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
440 –
441 –    /* If pts[i] is interior: done */
442 –    if(!below_plane[2])
443 –      {
444 –        if(test == GT_INTERIOR)
445 –          return(TRUE);
446 –      }
447 –
448 –    /* If all three points below separating plane: trivial reject */
449 –    if(below_plane[0] && below_plane[1] && below_plane[2])
450 –       return(FALSE);
451 –    /* Now check vertices in a against triangle b */
452 –    return(FALSE);
453 – }
454 –
455 –
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 469 | 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 482 | 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 495 | 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 511 | 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 525 | 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 539 | 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 555 | 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 568 | 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 581 | 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 592 | 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;
602 <  int tests[2][3];
603 <  FVECT n[2][3],p,avg[2];
604 <
605 <  /* Test the vertices of triangle a against the pyramid formed by triangle
606 <     b and the origin. If any vertex of a is interior to triangle b, or
607 <     if all 3 vertices of a are ON the edges of b,return TRUE. Remember
608 <     the results of the edge normal and sidedness tests for later.
609 <   */
610 < if(vertices_in_stri(a1,a2,a3,b1,b2,b3,&(n_set[0]),n[0],avg[0],sides[1]))
611 <     return(TRUE);
612 <  
613 < if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
614 <     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 */
629 <    if(tests[1][j] & tests[1][jnext])
630 <      continue;
631 <    for(i=0;i<3;i++)
632 <    {
633 <      inext = (i+1)%3;
634 <      /* IF edge a doesnt cross any great circles of b, punt */
635 <      if(tests[0][i] & tests[0][inext])
636 <        continue;
637 <      /* Now find the great circles that cross and test */
638 <      if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
639 <          && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
640 <      {
641 <        VCROSS(p,n[0][i],n[1][j]);
642 <                    
643 <        /* If zero cp= done */
644 <        if(ZERO_VEC3(p))
645 <          continue;
646 <        /* check above both planes */
647 <        if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
648 <          {
649 <            NEGATE_VEC3(p);
650 <            if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
651 <              continue;
652 <          }
653 <        return(TRUE);
654 <      }
655 <    }
656 <  }
657 <  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;
539   FVECT v0,v1,v2;
540   FVECT pt;
541   {
542 <  FVECT p0,p1,p2,p,n;
543 <  double pd;
542 >  FVECT p0,p1,p2,p;
543 >  FPEQ peq;
544    int type;
545  
546    VSUB(p0,v0,orig);
# Line 674 | Line 550 | FVECT pt;
550    if(point_in_stri(p0,p1,p2,dir))
551    {
552        /* Intersect the ray with the triangle plane */
553 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
554 <      return(intersect_ray_plane(orig,dir,n,pd,NULL,pt));
553 >      tri_plane_equation(v0,v1,v2,&peq,FALSE);
554 >      return(intersect_ray_plane(orig,dir,peq,NULL,pt));
555    }
556    return(FALSE);
557   }
# Line 776 | Line 652 | int p0id,p1id,p2id;
652      return(i);
653   }
654  
779 –
780 – int
781 – sedge_intersect(a0,a1,b0,b1)
782 – FVECT a0,a1,b0,b1;
783 – {
784 –    FVECT na,nb,avga,avgb,p;
785 –    double d;
786 –    int sb0,sb1,sa0,sa1;
787 –
788 –    /* First test if edge b straddles great circle of a */
789 –    VCROSS(na,a0,a1);
790 –    d = DOT(na,b0);
791 –    sb0 = ZERO(d)?0:(d<0)? -1: 1;
792 –    d = DOT(na,b1);
793 –    sb1 = ZERO(d)?0:(d<0)? -1: 1;
794 –    /* edge b entirely on one side of great circle a: edges cannot intersect*/
795 –    if(sb0*sb1 > 0)
796 –       return(FALSE);
797 –    /* test if edge a straddles great circle of b */
798 –    VCROSS(nb,b0,b1);
799 –    d = DOT(nb,a0);
800 –    sa0 = ZERO(d)?0:(d<0)? -1: 1;
801 –    d = DOT(nb,a1);
802 –    sa1 = ZERO(d)?0:(d<0)? -1: 1;
803 –    /* edge a entirely on one side of great circle b: edges cannot intersect*/
804 –    if(sa0*sa1 > 0)
805 –       return(FALSE);
806 –
807 –    /* Find one of intersection points of the great circles */
808 –    VCROSS(p,na,nb);
809 –    /* If they lie on same great circle: call an intersection */
810 –    if(ZERO_VEC3(p))
811 –       return(TRUE);
812 –
813 –    VADD(avga,a0,a1);
814 –    VADD(avgb,b0,b1);
815 –    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
816 –    {
817 –      NEGATE_VEC3(p);
818 –      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
819 –        return(FALSE);
820 –    }
821 –    if((!sb0 || !sb1) && (!sa0 || !sa1))
822 –      return(FALSE);
823 –    return(TRUE);
824 – }
825 –
826 –
827 –
655   /* Find the normalized barycentric coordinates of p relative to
656   * triangle v0,v1,v2. Return result in coord
657   */
# Line 842 | Line 669 | double coord[3];
669  
670   }
671  
845 – bary_ith_child(coord,i)
846 – double coord[3];
847 – int i;
848 – {
849 –
850 –  switch(i){
851 –  case 0:
852 –      /* update bary for child */
853 –      coord[0] = 2.0*coord[0]- 1.0;
854 –      coord[1] *= 2.0;
855 –      coord[2] *= 2.0;
856 –      break;
857 –  case 1:
858 –    coord[0] *= 2.0;
859 –    coord[1] = 2.0*coord[1]- 1.0;
860 –    coord[2] *= 2.0;
861 –    break;
862 –  case 2:
863 –    coord[0] *= 2.0;
864 –    coord[1] *= 2.0;
865 –    coord[2] = 2.0*coord[2]- 1.0;
866 –    break;
867 –  case 3:
868 –    coord[0] = 1.0 - 2.0*coord[0];
869 –    coord[1] = 1.0 - 2.0*coord[1];
870 –    coord[2] = 1.0 - 2.0*coord[2];
871 –    break;
872 – #ifdef DEBUG
873 –  default:
874 –    eputs("bary_ith_child():Invalid child\n");
875 –    break;
876 – #endif
877 –  }
878 – }
672  
673  
881 – int
882 – bary_child(coord)
883 – double coord[3];
884 – {
885 –  int i;
674  
887 –  if(coord[0] > 0.5)
888 –  {
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 –      return(0);
894 –  }
895 –  else
896 –    if(coord[1] > 0.5)
897 –    {
898 –      coord[0] *= 2.0;
899 –      coord[1] = 2.0*coord[1]- 1.0;
900 –      coord[2] *= 2.0;
901 –      return(1);
902 –    }
903 –    else
904 –      if(coord[2] > 0.5)
905 –      {
906 –        coord[0] *= 2.0;
907 –        coord[1] *= 2.0;
908 –        coord[2] = 2.0*coord[2]- 1.0;
909 –        return(2);
910 –      }
911 –      else
912 –         {
913 –           coord[0] = 1.0 - 2.0*coord[0];
914 –           coord[1] = 1.0 - 2.0*coord[1];
915 –           coord[2] = 1.0 - 2.0*coord[2];
916 –           return(3);
917 –         }
918 – }
919 –
920 – /* Coord was the ith child of the parent: set the coordinate
921 –   relative to the parent
922 – */
675   bary_parent(coord,i)
924 – double coord[3];
925 – int i;
926 – {
927 –
928 –  switch(i) {
929 –  case 0:
930 –    /* update bary for child */
931 –    coord[0] = coord[0]*0.5 + 0.5;
932 –    coord[1] *= 0.5;
933 –    coord[2] *= 0.5;
934 –    break;
935 –  case 1:
936 –    coord[0] *= 0.5;
937 –    coord[1]  = coord[1]*0.5 + 0.5;
938 –    coord[2] *= 0.5;
939 –    break;
940 –    
941 –  case 2:
942 –    coord[0] *= 0.5;
943 –    coord[1] *= 0.5;
944 –    coord[2] = coord[2]*0.5 + 0.5;
945 –    break;
946 –    
947 –  case 3:
948 –    coord[0] = 0.5 - 0.5*coord[0];
949 –    coord[1] = 0.5 - 0.5*coord[1];
950 –    coord[2] = 0.5 - 0.5*coord[2];
951 –    break;
952 – #ifdef DEBUG
953 –  default:
954 –    eputs("bary_parent():Invalid child\n");
955 –    break;
956 – #endif
957 –  }
958 – }
959 –
960 – bary_from_child(coord,child,next)
961 – double coord[3];
962 – int child,next;
963 – {
964 – #ifdef DEBUG
965 –  if(child <0 || child > 3)
966 –  {
967 –    eputs("bary_from_child():Invalid child\n");
968 –    return;
969 –  }
970 –  if(next <0 || next > 3)
971 –  {
972 –    eputs("bary_from_child():Invalid next\n");
973 –    return;
974 –  }
975 – #endif
976 –  if(next == child)
977 –    return;
978 –
979 –  switch(child){
980 –  case 0:
981 –    switch(next){
982 –    case 1:
983 –      coord[0] += 1.0;
984 –      coord[1] -= 1.0;
985 –      break;
986 –    case 2:
987 –      coord[0] += 1.0;
988 –      coord[2] -= 1.0;
989 –      break;
990 –    case 3:
991 –      coord[0] *= -1.0;
992 –      coord[1] = 1 - coord[1];
993 –      coord[2] = 1 - coord[2];
994 –      break;
995 –
996 –    }
997 –    break;
998 –  case 1:
999 –    switch(next){
1000 –    case 0:
1001 –      coord[0] -= 1.0;
1002 –      coord[1] += 1.0;
1003 –      break;
1004 –    case 2:
1005 –      coord[1] += 1.0;
1006 –      coord[2] -= 1.0;
1007 –      break;
1008 –    case 3:
1009 –      coord[0] = 1 - coord[0];
1010 –      coord[1] *= -1.0;
1011 –      coord[2] = 1 - coord[2];
1012 –      break;
1013 –    }
1014 –    break;
1015 –  case 2:
1016 –    switch(next){
1017 –    case 0:
1018 –      coord[0] -= 1.0;
1019 –      coord[2] += 1.0;
1020 –      break;
1021 –    case 1:
1022 –      coord[1] -= 1.0;
1023 –      coord[2] += 1.0;
1024 –      break;
1025 –    case 3:
1026 –      coord[0] = 1 - coord[0];
1027 –      coord[1] = 1 - coord[1];
1028 –      coord[2] *= -1.0;
1029 –      break;
1030 –    }
1031 –    break;
1032 –  case 3:
1033 –    switch(next){
1034 –    case 0:
1035 –      coord[0] *= -1.0;
1036 –      coord[1] = 1 - coord[1];
1037 –      coord[2] = 1 - coord[2];
1038 –      break;
1039 –    case 1:
1040 –      coord[0] = 1 - coord[0];
1041 –      coord[1] *= -1.0;
1042 –      coord[2] = 1 - coord[2];
1043 –      break;
1044 –    case 2:
1045 –      coord[0] = 1 - coord[0];
1046 –      coord[1] = 1 - coord[1];
1047 –      coord[2] *= -1.0;
1048 –      break;
1049 –    }
1050 –    break;
1051 –  }
1052 – }
1053 –
1054 –
1055 – baryi_parent(coord,i)
676   BCOORD coord[3];
677   int i;
678   {
1059 –
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 1111 | 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 1147 | 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 < /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
808 <   dir unbounded to [-MAXLONG,MAXLONG]
809 < */
810 < bary_dtol(b,db,bi,dbi,t)
1192 < double b[3],db[3][3];
1193 < BCOORD bi[3];
1194 < BDIR dbi[3][3];
1195 < TINT t[3];
807 > int
808 > bary_nth_child(coord,i)
809 > BCOORD coord[3];
810 > int i;
811   {
1197 –  int i,id,j;
1198 –  double d;
812  
813 <  for(i=0; i < 2;i++)
814 <  {
815 <    if(b[i] <= 0.0)
816 <    {
817 <      bi[i]= 0;
818 <    }
819 <    else
820 <      if(b[i] >= 1.0)
821 <      {
822 <        bi[i]= MAXBCOORD;
823 <      }
824 <      else
825 <        bi[i] = (BCOORD)(b[i]*MAXBCOORD);
813 >  switch(i){
814 >  case 0:
815 >    /* update bary for child */
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) - MAXBCOORD2;
823 >    coord[2] <<= 1;
824 >    break;
825 >  case 2:
826 >    coord[0] <<= 1;
827 >    coord[1] <<= 1;
828 >    coord[2] = (coord[2] << 1) - MAXBCOORD2;
829 >    break;
830 >  case 3:
831 >    coord[0] = MAXBCOORD2 - (coord[0] << 1);
832 >    coord[1] = MAXBCOORD2 - (coord[1] << 1);
833 >    coord[2] = MAXBCOORD2 - (coord[2] << 1);
834 >    break;
835    }
1214 –  bi[2] = MAXBCOORD - bi[0] - bi[1];
1215 –
1216 –  if(bi[2] < 0)
1217 –  {
1218 –      bi[2] = 0;
1219 –      bi[1] = MAXBCOORD - bi[0];
1220 –  }
1221 –  for(j=0; j< 3; j++)
1222 –  {
1223 –    if(t[j]==0)
1224 –      continue;
1225 –    if(t[j] == HUGET)
1226 –       max_index(db[j],&d);
1227 –    for(i=0; i< 3; i++)
1228 –       if(t[j] != HUGET)
1229 –          dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1230 –       else
1231 –          dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1232 –  }
836   }
1234 –
1235 –
1236 –
1237 –
1238 –
1239 –
837  
838  
839  

Diff Legend

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