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.5 by gwlarson, Mon Sep 14 10:33:46 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 736 | Line 612 | FVECT fnear[4],ffar[4];
612      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
613   }
614  
739 –
740 –
741 –
615   int
616 < sedge_intersect(a0,a1,b0,b1)
617 < FVECT a0,a1,b0,b1;
616 > max_index(v,r)
617 > FVECT v;
618 > double *r;
619   {
620 <    FVECT na,nb,avga,avgb,p;
621 <    double d;
748 <    int sb0,sb1,sa0,sa1;
620 >  double p[3];
621 >  int i;
622  
623 <    /* First test if edge b straddles great circle of a */
624 <    VCROSS(na,a0,a1);
625 <    d = DOT(na,b0);
626 <    sb0 = ZERO(d)?0:(d<0)? -1: 1;
627 <    d = DOT(na,b1);
628 <    sb1 = ZERO(d)?0:(d<0)? -1: 1;
629 <    /* edge b entirely on one side of great circle a: edges cannot intersect*/
630 <    if(sb0*sb1 > 0)
758 <       return(FALSE);
759 <    /* test if edge a straddles great circle of b */
760 <    VCROSS(nb,b0,b1);
761 <    d = DOT(nb,a0);
762 <    sa0 = ZERO(d)?0:(d<0)? -1: 1;
763 <    d = DOT(nb,a1);
764 <    sa1 = ZERO(d)?0:(d<0)? -1: 1;
765 <    /* edge a entirely on one side of great circle b: edges cannot intersect*/
766 <    if(sa0*sa1 > 0)
767 <       return(FALSE);
623 >  p[0] = fabs(v[0]);
624 >  p[1] = fabs(v[1]);
625 >  p[2] = fabs(v[2]);
626 >  i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);  
627 >  if(r)
628 >    *r = p[i];
629 >  return(i);
630 > }
631  
632 <    /* Find one of intersection points of the great circles */
633 <    VCROSS(p,na,nb);
634 <    /* If they lie on same great circle: call an intersection */
635 <    if(ZERO_VEC3(p))
636 <       return(TRUE);
637 <
638 <    VADD(avga,a0,a1);
639 <    VADD(avgb,b0,b1);
640 <    if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
632 > int
633 > closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
634 > FVECT p0,p1,p2,p;
635 > int p0id,p1id,p2id;
636 > {
637 >    double d,d1;
638 >    int i;
639 >    
640 >    d =  DIST_SQ(p,p0);
641 >    d1 = DIST_SQ(p,p1);
642 >    if(d < d1)
643      {
644 <      NEGATE_VEC3(p);
645 <      if(DOT(avga,p) < 0 || DOT(avgb,p) < 0)
781 <        return(FALSE);
644 >      d1 = DIST_SQ(p,p2);
645 >      i = (d1 < d)?p2id:p0id;
646      }
647 <    if((!sb0 || !sb1) && (!sa0 || !sa1))
648 <      return(FALSE);
649 <    return(TRUE);
647 >    else
648 >    {
649 >      d = DIST_SQ(p,p2);
650 >      i = (d < d1)? p2id:p1id;
651 >    }
652 >    return(i);
653   }
654  
788 –
789 –
655   /* Find the normalized barycentric coordinates of p relative to
656   * triangle v0,v1,v2. Return result in coord
657   */
# Line 800 | Line 665 | double coord[3];
665    a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
666    coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
667    coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
668 <  coord[2]  = 1.0 - coord[0] - coord[1];
668 >  coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
669  
670   }
671  
807 – bary_ith_child(coord,i)
808 – double coord[3];
809 – int i;
810 – {
811 –
812 –  switch(i){
813 –  case 0:
814 –      /* update bary for child */
815 –      coord[0] = 2.0*coord[0]- 1.0;
816 –      coord[1] *= 2.0;
817 –      coord[2] *= 2.0;
818 –      break;
819 –  case 1:
820 –    coord[0] *= 2.0;
821 –    coord[1] = 2.0*coord[1]- 1.0;
822 –    coord[2] *= 2.0;
823 –    break;
824 –  case 2:
825 –    coord[0] *= 2.0;
826 –    coord[1] *= 2.0;
827 –    coord[2] = 2.0*coord[2]- 1.0;
828 –    break;
829 –  case 3:
830 –    coord[0] = 1.0 - 2.0*coord[0];
831 –    coord[1] = 1.0 - 2.0*coord[1];
832 –    coord[2] = 1.0 - 2.0*coord[2];
833 –    break;
834 – #ifdef DEBUG
835 –  default:
836 –    eputs("bary_ith_child():Invalid child\n");
837 –    break;
838 – #endif
839 –  }
840 – }
672  
673  
843 – int
844 – bary_child(coord)
845 – double coord[3];
846 – {
847 –  int i;
674  
849 –  if(coord[0] > 0.5)
850 –  {
851 –      /* update bary for child */
852 –      coord[0] = 2.0*coord[0]- 1.0;
853 –      coord[1] *= 2.0;
854 –      coord[2] *= 2.0;
855 –      return(0);
856 –  }
857 –  else
858 –    if(coord[1] > 0.5)
859 –    {
860 –      coord[0] *= 2.0;
861 –      coord[1] = 2.0*coord[1]- 1.0;
862 –      coord[2] *= 2.0;
863 –      return(1);
864 –    }
865 –    else
866 –      if(coord[2] > 0.5)
867 –      {
868 –        coord[0] *= 2.0;
869 –        coord[1] *= 2.0;
870 –        coord[2] = 2.0*coord[2]- 1.0;
871 –        return(2);
872 –      }
873 –      else
874 –         {
875 –           coord[0] = 1.0 - 2.0*coord[0];
876 –           coord[1] = 1.0 - 2.0*coord[1];
877 –           coord[2] = 1.0 - 2.0*coord[2];
878 –           return(3);
879 –         }
880 – }
881 –
882 – /* Coord was the ith child of the parent: set the coordinate
883 –   relative to the parent
884 – */
675   bary_parent(coord,i)
676 < double coord[3];
676 > BCOORD coord[3];
677   int i;
678   {
889 –
679    switch(i) {
680    case 0:
681      /* update bary for child */
682 <    coord[0] = coord[0]*0.5 + 0.5;
683 <    coord[1] *= 0.5;
684 <    coord[2] *= 0.5;
682 >    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
683 >    coord[1] >>= 1;
684 >    coord[2] >>= 1;
685      break;
686    case 1:
687 <    coord[0] *= 0.5;
688 <    coord[1]  = coord[1]*0.5 + 0.5;
689 <    coord[2] *= 0.5;
687 >    coord[0] >>= 1;
688 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
689 >    coord[2] >>= 1;
690      break;
691      
692    case 2:
693 <    coord[0] *= 0.5;
694 <    coord[1] *= 0.5;
695 <    coord[2] = coord[2]*0.5 + 0.5;
693 >    coord[0] >>= 1;
694 >    coord[1] >>= 1;
695 >    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
696      break;
697      
698    case 3:
699 <    coord[0] = 0.5 - 0.5*coord[0];
700 <    coord[1] = 0.5 - 0.5*coord[1];
701 <    coord[2] = 0.5 - 0.5*coord[2];
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:
# Line 920 | Line 709 | int i;
709   }
710  
711   bary_from_child(coord,child,next)
712 < double coord[3];
712 > BCOORD coord[3];
713   int child,next;
714   {
715   #ifdef DEBUG
# Line 940 | Line 729 | int child,next;
729  
730    switch(child){
731    case 0:
732 <    switch(next){
733 <    case 1:
734 <      coord[0] += 1.0;
946 <      coord[1] -= 1.0;
732 >      coord[0] = 0;
733 >      coord[1] = MAXBCOORD2 - coord[1];
734 >      coord[2] = MAXBCOORD2 - coord[2];
735        break;
948 –    case 2:
949 –      coord[0] += 1.0;
950 –      coord[2] -= 1.0;
951 –      break;
952 –    case 3:
953 –      coord[0] *= -1.0;
954 –      coord[1] = 1 - coord[1];
955 –      coord[2] = 1 - coord[2];
956 –      break;
957 –
958 –    }
959 –    break;
736    case 1:
737 <    switch(next){
738 <    case 0:
739 <      coord[0] -= 1.0;
964 <      coord[1] += 1.0;
737 >      coord[0] = MAXBCOORD2 - coord[0];
738 >      coord[1] = 0;
739 >      coord[2] = MAXBCOORD2 - coord[2];
740        break;
966 –    case 2:
967 –      coord[1] += 1.0;
968 –      coord[2] -= 1.0;
969 –      break;
970 –    case 3:
971 –      coord[0] = 1 - coord[0];
972 –      coord[1] *= -1.0;
973 –      coord[2] = 1 - coord[2];
974 –      break;
975 –    }
976 –    break;
741    case 2:
742 <    switch(next){
743 <    case 0:
744 <      coord[0] -= 1.0;
981 <      coord[2] += 1.0;
982 <      break;
983 <    case 1:
984 <      coord[1] -= 1.0;
985 <      coord[2] += 1.0;
986 <      break;
987 <    case 3:
988 <      coord[0] = 1 - coord[0];
989 <      coord[1] = 1 - coord[1];
990 <      coord[2] *= -1.0;
991 <      break;
992 <    }
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] *= -1.0;
750 <      coord[1] = 1 - coord[1];
751 <      coord[2] = 1 - coord[2];
749 >      coord[0] = 0;
750 >      coord[1] = MAXBCOORD2 - coord[1];
751 >      coord[2] = MAXBCOORD2 - coord[2];
752        break;
753      case 1:
754 <      coord[0] = 1 - coord[0];
755 <      coord[1] *= -1.0;
756 <      coord[2] = 1 - coord[2];
754 >      coord[0] = MAXBCOORD2 - coord[0];
755 >      coord[1] = 0;
756 >      coord[2] = MAXBCOORD2 - coord[2];
757        break;
758      case 2:
759 <      coord[0] = 1 - coord[0];
760 <      coord[1] = 1 - coord[1];
761 <      coord[2] *= -1.0;
759 >      coord[0] = MAXBCOORD2 - coord[0];
760 >      coord[1] = MAXBCOORD2 - coord[1];
761 >      coord[2] = 0;
762        break;
763      }
764      break;
# Line 1014 | Line 766 | int child,next;
766   }
767  
768   int
769 < max_index(v)
770 < FVECT v;
769 > bary_child(coord)
770 > BCOORD coord[3];
771   {
1020 –  double a,b,c;
772    int i;
773  
774 <  a = fabs(v[0]);
775 <  b = fabs(v[1]);
776 <  c = fabs(v[2]);
777 <  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
778 <  return(i);
779 < }
780 <
781 < int
782 < closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
783 < FVECT p0,p1,p2,p;
1033 < int p0id,p1id,p2id;
1034 < {
1035 <    double d,d1;
1036 <    int i;
1037 <    
1038 <    d =  DIST_SQ(p,p0);
1039 <    d1 = DIST_SQ(p,p1);
1040 <    if(d < d1)
774 >  if(coord[0] > MAXBCOORD4)
775 >  {
776 >      /* update bary for child */
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] > MAXBCOORD4)
784      {
785 <      d1 = DIST_SQ(p,p2);
786 <      i = (d1 < d)?p2id:p0id;
785 >      coord[0] <<= 1;
786 >      coord[1] = (coord[1] << 1) - MAXBCOORD2;
787 >      coord[2] <<= 1;
788 >      return(1);
789      }
790      else
791 <    {
792 <      d = DIST_SQ(p,p2);
793 <      i = (d < d1)? p2id:p1id;
794 <    }
795 <    return(i);
791 >      if(coord[2] > MAXBCOORD4)
792 >      {
793 >        coord[0] <<= 1;
794 >        coord[1] <<= 1;
795 >        coord[2] = (coord[2] << 1) - MAXBCOORD2;
796 >        return(2);
797 >      }
798 >      else
799 >         {
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 + bary_nth_child(coord,i)
809 + BCOORD coord[3];
810 + int i;
811 + {
812  
813 <
814 <
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 >  }
836 > }
837  
838  
839  

Diff Legend

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