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.3 by gwlarson, Tue Aug 25 11:03: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)
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   */
82   double
83   tri_normal(v0,v1,v2,n,norm)
84   FVECT v0,v1,v2,n;
85 < char norm;
85 > int norm;
86   {
87    double mag;
88  
# Line 63 | Line 92 | char 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 | char norm;
100  
101    if(!norm)
102       return(0);
75 –
103    
104    mag = normalize(n);
105  
# Line 80 | Line 107 | char norm;
107   }
108  
109  
110 < tri_plane_equation(v0,v1,v2,n,nd,norm)
111 <   FVECT v0,v1,v2,n;
112 <   double *nd;
113 <   char norm;
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 – int
94 – point_relative_to_plane(p,n,nd)
95 –   FVECT p,n;
96 –   double nd;
97 – {
98 –    double d;
99 –    
100 –    d = p[0]*n[0] + p[1]*n[1] + p[2]*n[2] + nd;
101 –    if(d < 0)
102 –       return(-1);
103 –    if(ZERO(d))
104 –       return(0);
105 –    else
106 –       return(1);
107 – }
120  
121 < /* From quad_edge-code */
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 < point_in_circle_thru_origin(p,p0,p1)
128 < FVECT p;
129 < FVECT p0,p1;
114 < {
115 <
116 <    double dp0,dp1;
117 <    double dp,det;
118 <    
119 <    dp0 = DOT_VEC2(p0,p0);
120 <    dp1 = DOT_VEC2(p1,p1);
121 <    dp  = DOT_VEC2(p,p);    
122 <
123 <    det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
124 <    
125 <    return (det > 0);
126 < }
127 <
128 <
129 <
130 < point_on_sphere(ps,p,c)
131 < FVECT ps,p,c;
132 < {
133 <    VSUB(ps,p,c);    
134 <    normalize(ps);
135 < }
136 <
137 <
138 < int
139 < intersect_vector_plane(v,plane_n,plane_d,tptr,r)
140 <   FVECT v,plane_n;
141 <   double plane_d;
127 > intersect_vector_plane(v,peq,tptr,r)
128 >   FVECT v;
129 >   FPEQ peq;
130     double *tptr;
131     FVECT r;
132   {
133 <  double t;
133 >  double t,d;
134    int hit;
135      /*
136        Plane is Ax + By + Cz +D = 0:
# Line 152 | 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 <    t =  plane_d/-(DOT(plane_n,v));
144 <    if(t >0 || ZERO(t))
145 <       hit = 1;
146 <    else
147 <       hit = 0;
148 <    r[0] = v[0]*t;
149 <    r[1] = v[1]*t;
150 <    r[2] = v[2]*t;
143 >  d = -(DOT(FP_N(peq),v));
144 >  if(ZERO(d))
145 >  {
146 >      t = FHUGE;
147 >      hit = 0;
148 >  }
149 >  else
150 >  {
151 >      t =  FP_D(peq)/d;
152 >      if(t < 0 )
153 >         hit = 0;
154 >      else
155 >         hit = 1;
156 >      if(r)
157 >         {
158 >             r[0] = v[0]*t;
159 >             r[1] = v[1]*t;
160 >             r[2] = v[2]*t;
161 >         }
162 >  }
163      if(tptr)
164         *tptr = t;
165    return(hit);
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;
172 <   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 184 | 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(ZERO(t) || t >0)
188 <       hit = 1;
189 <    else
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;
195  
196 <  VSUM(r,orig,dir,t);
196 >  if(r)
197 >     VSUM(r,orig,dir,t);
198  
199      if(pd)
200         *pd = t;
# Line 199 | Line 203 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
203  
204  
205   int
206 < point_in_cone(p,p0,p1,p2)
207 < FVECT p;
208 < FVECT p0,p1,p2;
206 > intersect_ray_oplane(orig,dir,n,pd,r)
207 >   FVECT orig,dir;
208 >   FVECT n;
209 >   double *pd;
210 >   FVECT r;
211   {
212 <    FVECT n;
213 <    FVECT np,x_axis,y_axis;
214 <    double d1,d2,d;
215 <    
216 <    /* Find the equation of the circle defined by the intersection
217 <       of the cone with the plane defined by p1,p2,p3- project p into
218 <       that plane and do an in-circle test in the plane
212 >  double t,d;
213 >  int hit;
214 >    /*
215 >      Plane is Ax + By + Cz +D = 0:
216 >      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
217 >    */
218 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
219 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
220 >       line is  l = p1 + (p2-p1)t
221       */
222 <    
223 <    /* find the equation of the plane defined by p1-p3 */
224 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
222 >    /* Solve for t: */
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 <    /* define a coordinate system on the plane: the x axis is in
233 <       the direction of np2-np1, and the y axis is calculated from
220 <       n cross x-axis
221 <     */
222 <    /* Project p onto the plane */
223 <    if(!intersect_vector_plane(p,n,d,NULL,np))
224 <        return(FALSE);
232 >  if(r)
233 >     VSUM(r,orig,dir,t);
234  
235 <    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
236 <    VSUB(x_axis,p1,p0);
237 <    normalize(x_axis);
238 <    /* The y axis is  */
230 <    VCROSS(y_axis,n,x_axis);
231 <    normalize(y_axis);
235 >    if(pd)
236 >       *pd = t;
237 >  return(hit);
238 > }
239  
240 <    VSUB(p1,p1,p0);
241 <    VSUB(p2,p2,p0);
242 <    VSUB(np,np,p0);
243 <    
244 <    p1[0] = VLEN(p1);
245 <    p1[1] = 0;
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 <    d1 = DOT(p2,x_axis);
251 <    d2 = DOT(p2,y_axis);
252 <    p2[0] = d1;
253 <    p2[1] = d2;
254 <    
255 <    d1 = DOT(np,x_axis);
256 <    d2 = DOT(np,y_axis);
257 <    np[0] = d1;
248 <    np[1] = d2;
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 <    /* perform the in-circle test in the new coordinate system */
260 <    return(point_in_circle_thru_origin(np,p1,p2));
259 > int
260 > intersect_edge_plane(e0,e1,peq,pd,r)
261 >   FVECT e0,e1;
262 >   FPEQ peq;
263 >   double *pd;
264 >   FVECT r;
265 > {
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 >    /* 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 >  VSUM(r,e0,d,t);
286 >
287 >    if(pd)
288 >       *pd = t;
289 >  return(hit);
290   }
291  
292   int
293 < test_point_against_spherical_tri(v0,v1,v2,p,n,nset,which,sides)
293 > point_set_in_stri(v0,v1,v2,p,n,nset,sides)
294   FVECT v0,v1,v2,p;
295   FVECT n[3];
296 < char *nset;
297 < char *which;
260 < char sides[3];
296 > int *nset;
297 > int sides[3];
298  
299   {
300 <    float d;
264 <
300 >    double d;
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 */
308      d  = DOT(n[0],p);
309  
310 <    if(ZERO(d))
311 <       sides[0] = GT_EDGE;
312 <    else
313 <       if(d > 0)
314 <      {
279 <          sides[0] =  GT_OUT;
280 <          sides[1] = sides[2] = GT_INVALID;
281 <          return(FALSE);
310 >    if(d > 0.0)
311 >     {
312 >       sides[0] =  GT_OUT;
313 >       sides[1] = sides[2] = GT_INVALID;
314 >       return(FALSE);
315        }
316      else
317         sides[0] = GT_INTERIOR;
# Line 286 | Line 319 | char 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 */
326      d  = DOT(n[1],p);
327 <    if(ZERO(d))
327 >    if(d > 0.0)
328      {
296 –        sides[1] = GT_EDGE;
297 –        /* If on plane 0-and on plane 1: lies on edge */
298 –        if(sides[0] == GT_EDGE)
299 –        {
300 –            *which = 1;
301 –            sides[2] = GT_INVALID;
302 –            return(GT_EDGE);
303 –        }
304 –    }
305 –    else if(d > 0)
306 –    {
329          sides[1] = GT_OUT;
330          sides[2] = GT_INVALID;
331          return(FALSE);
# Line 313 | Line 335 | char sides[3];
335      /* Test next edge */
336      if(!NTH_BIT(*nset,2))
337      {
338 <
317 <        VCROSS(n[2],v0,v2);
338 >        VCROSS(n[2],v2,v0);
339          SET_NTH_BIT(*nset,2);
340      }
341      /* Test the point for sidedness */
342      d  = DOT(n[2],p);
343 <    if(ZERO(d))
343 >    if(d > 0.0)
344      {
345 <        sides[2] = GT_EDGE;
346 <
326 <        /* If on plane 0 and 2: lies on edge 0*/
327 <        if(sides[0] == GT_EDGE)
328 <           {
329 <               *which = 0;
330 <               return(GT_EDGE);
331 <           }
332 <        /* If on plane 1 and 2: lies on edge  2*/
333 <        if(sides[1] == GT_EDGE)
334 <           {
335 <               *which = 2;
336 <               return(GT_EDGE);
337 <           }
338 <        /* otherwise: on face 2 */
339 <        else
340 <           {
341 <               *which = 2;
342 <               return(GT_FACE);
343 <           }
345 >      sides[2] = GT_OUT;
346 >      return(FALSE);
347      }
345 –    else if(d > 0)
346 –      {
347 –        sides[2] = GT_OUT;
348 –        return(FALSE);
349 –      }
350 –    /* If on edge */
348      else
349         sides[2] = GT_INTERIOR;
353 –    
354 –    /* If on plane 0 only: on face 0 */
355 –    if(sides[0] == GT_EDGE)
356 –    {
357 –        *which = 0;
358 –        return(GT_FACE);
359 –    }
360 –    /* If on plane 1 only: on face 1 */
361 –    if(sides[1] == GT_EDGE)
362 –    {
363 –        *which = 1;
364 –        return(GT_FACE);
365 –    }
350      /* Must be interior to the pyramid */
351      return(GT_INTERIOR);
352   }
353  
354  
355  
372 –
373 – int
374 – test_single_point_against_spherical_tri(v0,v1,v2,p,which)
375 – FVECT v0,v1,v2,p;
376 – char *which;
377 – {
378 –    float d;
379 –    FVECT n;  
380 –    char sides[3];
381 –
382 –    /* First test if point coincides with any of the vertices */
383 –    if(EQUAL_VEC3(p,v0))
384 –    {
385 –        *which = 0;
386 –        return(GT_VERTEX);
387 –    }
388 –    if(EQUAL_VEC3(p,v1))
389 –    {
390 –        *which = 1;
391 –        return(GT_VERTEX);
392 –    }
393 –    if(EQUAL_VEC3(p,v2))
394 –    {
395 –        *which = 2;
396 –        return(GT_VERTEX);
397 –    }
398 –    VCROSS(n,v1,v0);
399 –    /* Test the point for sidedness */
400 –    d  = DOT(n,p);
401 –    if(ZERO(d))
402 –       sides[0] = GT_EDGE;
403 –    else
404 –       if(d > 0)
405 –          return(FALSE);
406 –       else
407 –          sides[0] = GT_INTERIOR;
408 –    /* Test next edge */
409 –    VCROSS(n,v2,v1);
410 –    /* Test the point for sidedness */
411 –    d  = DOT(n,p);
412 –    if(ZERO(d))
413 –    {
414 –        sides[1] = GT_EDGE;
415 –        /* If on plane 0-and on plane 1: lies on edge */
416 –        if(sides[0] == GT_EDGE)
417 –        {
418 –            *which = 1;
419 –            return(GT_VERTEX);
420 –        }
421 –    }
422 –    else if(d > 0)
423 –       return(FALSE);
424 –    else
425 –       sides[1] = GT_INTERIOR;
426 –
427 –    /* Test next edge */
428 –    VCROSS(n,v0,v2);
429 –    /* Test the point for sidedness */
430 –    d  = DOT(n,p);
431 –    if(ZERO(d))
432 –    {
433 –        sides[2] = GT_EDGE;
434 –        
435 –        /* If on plane 0 and 2: lies on edge 0*/
436 –        if(sides[0] == GT_EDGE)
437 –        {
438 –            *which = 0;
439 –            return(GT_VERTEX);
440 –        }
441 –        /* If on plane 1 and 2: lies on edge  2*/
442 –        if(sides[1] == GT_EDGE)
443 –        {
444 –            *which = 2;
445 –            return(GT_VERTEX);
446 –        }
447 –        /* otherwise: on face 2 */
448 –        else
449 –       {
450 –           return(GT_FACE);
451 –       }
452 –    }
453 –    else if(d > 0)
454 –       return(FALSE);
455 –    /* Must be interior to the pyramid */
456 –    return(GT_FACE);
457 – }
458 –
459 – int
460 – test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
461 – FVECT t0,t1,t2,p0,p1,p2;
462 – char *nset;
463 – FVECT n[3];
464 – FVECT avg;
465 – char pt_sides[3][3];
466 –
467 – {
468 –    char below_plane[3],on_edge,test;
469 –    char which;
470 –
471 –    SUM_3VEC3(avg,t0,t1,t2);
472 –    on_edge = 0;
473 –    *nset = 0;
474 –    /* Test vertex v[i] against triangle j*/
475 –    /* Check if v[i] lies below plane defined by avg of 3 vectors
476 –       defining triangle
477 –       */
478 –
479 –    /* test point 0 */
480 –    if(DOT(avg,p0) < 0)
481 –      below_plane[0] = 1;
482 –    else
483 –      below_plane[0]=0;
484 –    /* Test if b[i] lies in or on triangle a */
485 –    test = test_point_against_spherical_tri(t0,t1,t2,p0,
486 –                                                 n,nset,&which,pt_sides[0]);
487 –    /* If pts[i] is interior: done */
488 –    if(!below_plane[0])
489 –      {
490 –        if(test == GT_INTERIOR)
491 –          return(TRUE);
492 –        /* Remember if b[i] fell on one of the 3 defining planes */
493 –        if(test)
494 –          on_edge++;
495 –      }
496 –    /* Now test point 1*/
497 –
498 –    if(DOT(avg,p1) < 0)
499 –      below_plane[1] = 1;
500 –    else
501 –      below_plane[1]=0;
502 –    /* Test if b[i] lies in or on triangle a */
503 –    test = test_point_against_spherical_tri(t0,t1,t2,p1,
504 –                                                 n,nset,&which,pt_sides[1]);
505 –    /* If pts[i] is interior: done */
506 –    if(!below_plane[1])
507 –    {
508 –      if(test == GT_INTERIOR)
509 –        return(TRUE);
510 –      /* Remember if b[i] fell on one of the 3 defining planes */
511 –      if(test)
512 –        on_edge++;
513 –    }
514 –    
515 –    /* Now test point 2 */
516 –    if(DOT(avg,p2) < 0)
517 –      below_plane[2] = 1;
518 –    else
519 –      below_plane[2]=0;
520 –        /* Test if b[i] lies in or on triangle a */
521 –    test = test_point_against_spherical_tri(t0,t1,t2,p2,
522 –                                                 n,nset,&which,pt_sides[2]);
523 –
524 –    /* If pts[i] is interior: done */
525 –    if(!below_plane[2])
526 –      {
527 –        if(test == GT_INTERIOR)
528 –          return(TRUE);
529 –        /* Remember if b[i] fell on one of the 3 defining planes */
530 –        if(test)
531 –          on_edge++;
532 –      }
533 –
534 –    /* If all three points below separating plane: trivial reject */
535 –    if(below_plane[0] && below_plane[1] && below_plane[2])
536 –       return(FALSE);
537 –    /* Accept if all points lie on a triangle vertex/edge edge- accept*/
538 –    if(on_edge == 3)
539 –       return(TRUE);
540 –    /* Now check vertices in a against triangle b */
541 –    return(FALSE);
542 – }
543 –
544 –
356   set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
357     FVECT t0,t1,t2,p0,p1,p2;
358 <   char test[3];
359 <   char sides[3][3];
360 <   char nset;
358 >   int test[3];
359 >   int sides[3][3];
360 >   int nset;
361     FVECT n[3];
362   {
363 <    char t;
363 >    int t;
364      double d;
365  
366      
# Line 558 | 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)
375 >      if(d >= 0.0)
376          SET_NTH_BIT(test[0],0);
377      }
378      else
# Line 571 | 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)
388 >        if(d >= 0.0)
389            SET_NTH_BIT(test[0],1);
390      }
391      else
# Line 584 | 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)
401 >      if(d >= 0.0)
402          SET_NTH_BIT(test[0],2);
403      }
404      else
# Line 600 | 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)
417 >      if(d >= 0.0)
418          SET_NTH_BIT(test[1],0);
419      }
420      else
# Line 614 | 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)
431 >      if(d >= 0.0)
432          SET_NTH_BIT(test[1],1);
433      }
434      else
# Line 628 | 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)
445 >      if(d >= 0.0)
446          SET_NTH_BIT(test[1],2);
447      }
448      else
# Line 644 | 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)
461 >      if(d >= 0.0)
462          SET_NTH_BIT(test[2],0);
463      }
464      else
# Line 657 | 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)
474 >      if(d >= 0.0)
475          SET_NTH_BIT(test[2],1);
476      }
477      else
# Line 670 | 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)
487 >      if(d >= 0.0)
488          SET_NTH_BIT(test[2],2);
489      }
490      else
# Line 681 | 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 < spherical_tri_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 <  char which,test,n_set[2];
510 <  char sides[2][3][3],i,j,inext,jnext;
691 <  char tests[2][3];
692 <  FVECT n[2][3],p,avg[2];
693 <
694 <  /* Test the vertices of triangle a against the pyramid formed by triangle
695 <     b and the origin. If any vertex of a is interior to triangle b, or
696 <     if all 3 vertices of a are ON the edges of b,return TRUE. Remember
697 <     the results of the edge normal and sidedness tests for later.
698 <   */
699 < if(test_vertices_for_tri_inclusion(a1,a2,a3,b1,b2,b3,
700 <                                    &(n_set[0]),n[0],avg[0],sides[1]))
701 <     return(TRUE);
702 <  
703 < if(test_vertices_for_tri_inclusion(b1,b2,b3,a1,a2,a3,
704 <                                    &(n_set[1]),n[1],avg[1],sides[0]))
705 <     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 */
720 <    if(tests[1][j] & tests[1][jnext])
721 <      continue;
722 <    for(i=0;i<3;i++)
723 <    {
724 <      inext = (i+1)%3;
725 <      /* IF edge a doesnt cross any great circles of b, punt */
726 <      if(tests[0][i] & tests[0][inext])
727 <        continue;
728 <      /* Now find the great circles that cross and test */
729 <      if((NTH_BIT(tests[0][i],j)^(NTH_BIT(tests[0][inext],j)))
730 <          && (NTH_BIT(tests[1][j],i)^NTH_BIT(tests[1][jnext],i)))
731 <      {
732 <        VCROSS(p,n[0][i],n[1][j]);
733 <                    
734 <        /* If zero cp= done */
735 <        if(ZERO_VEC3(p))
736 <          continue;
737 <        /* check above both planes */
738 <        if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
739 <          {
740 <            NEGATE_VEC3(p);
741 <            if(DOT(avg[0],p) < 0 || DOT(avg[1],p) < 0)
742 <              continue;
743 <          }
744 <        return(TRUE);
745 <      }
746 <    }
747 <  }
748 <  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,wptr)
537 > ray_intersect_tri(orig,dir,v0,v1,v2,pt)
538   FVECT orig,dir;
539   FVECT v0,v1,v2;
540   FVECT pt;
756 – char *wptr;
541   {
542 <  FVECT p0,p1,p2,p,n;
543 <  char type,which;
544 <  double pd;
761 <  
762 <  point_on_sphere(p0,v0,orig);
763 <  point_on_sphere(p1,v1,orig);
764 <  point_on_sphere(p2,v2,orig);
765 <  type = test_single_point_against_spherical_tri(p0,p1,p2,dir,&which);
542 >  FVECT p0,p1,p2,p;
543 >  FPEQ peq;
544 >  int type;
545  
546 <  if(type)
546 >  VSUB(p0,v0,orig);
547 >  VSUB(p1,v1,orig);
548 >  VSUB(p2,v2,orig);
549 >
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 <      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 <  if(wptr)
774 <    *wptr = which;
775 <
776 <  return(type);
556 >  return(FALSE);
557   }
558  
559  
# Line 832 | Line 612 | FVECT fnear[4],ffar[4];
612      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
613   }
614  
835 –
836 –
837 –
615   int
616 < spherical_polygon_edge_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;
844 <    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)
854 <       return(FALSE);
855 <    /* test if edge a straddles great circle of b */
856 <    VCROSS(nb,b0,b1);
857 <    d = DOT(nb,a0);
858 <    sa0 = ZERO(d)?0:(d<0)? -1: 1;
859 <    d = DOT(nb,a1);
860 <    sa1 = ZERO(d)?0:(d<0)? -1: 1;
861 <    /* edge a entirely on one side of great circle b: edges cannot intersect*/
862 <    if(sa0*sa1 > 0)
863 <       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)
877 <        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  
884 –
885 –
655   /* Find the normalized barycentric coordinates of p relative to
656   * triangle v0,v1,v2. Return result in coord
657   */
# Line 896 | 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  
903 – int
904 – bary2d_child(coord)
905 – double coord[3];
906 – {
907 –  int i;
672  
909 –  /* First check if one of the original vertices */
910 –  for(i=0;i<3;i++)
911 –    if(EQUAL(coord[i],1.0))
912 –      return(i);
673  
674 <  /* Check if one of the new vertices: for all return center child */
675 <  if(ZERO(coord[0]) && EQUAL(coord[1],0.5))
674 >
675 > bary_parent(coord,i)
676 > BCOORD coord[3];
677 > int i;
678 > {
679 >  switch(i) {
680 >  case 0:
681 >    /* update bary for child */
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) + 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) + MAXBCOORD4;
696 >    break;
697 >    
698 >  case 3:
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("bary_parent():Invalid child\n");
706 >    break;
707 > #endif
708 >  }
709 > }
710 >
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 <    coord[0] = 1.0f;
719 <    coord[1] = 0.0f;
919 <    coord[2] = 0.0f;
920 <    return(3);
718 >    eputs("bary_from_child():Invalid child\n");
719 >    return;
720    }
721 <  if(ZERO(coord[1]) && EQUAL(coord[0],0.5))
721 >  if(next <0 || next > 3)
722    {
723 <    coord[0] = 0.0f;
724 <    coord[1] = 1.0f;
926 <    coord[2] = 0.0f;
927 <    return(3);
723 >    eputs("bary_from_child():Invalid next\n");
724 >    return;
725    }
726 <  if(ZERO(coord[2]) && EQUAL(coord[0],0.5))
727 <    {
728 <      coord[0] = 0.0f;
729 <      coord[1] = 0.0f;
730 <      coord[2] = 1.0f;
731 <      return(3);
726 > #endif
727 >  if(next == child)
728 >    return;
729 >
730 >  switch(child){
731 >  case 0:
732 >      coord[0] = 0;
733 >      coord[1] = MAXBCOORD2 - coord[1];
734 >      coord[2] = MAXBCOORD2 - coord[2];
735 >      break;
736 >  case 1:
737 >      coord[0] = MAXBCOORD2 - coord[0];
738 >      coord[1] = 0;
739 >      coord[2] = MAXBCOORD2 - coord[2];
740 >      break;
741 >  case 2:
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] = MAXBCOORD2 - coord[1];
751 >      coord[2] = MAXBCOORD2 - coord[2];
752 >      break;
753 >    case 1:
754 >      coord[0] = MAXBCOORD2 - coord[0];
755 >      coord[1] = 0;
756 >      coord[2] = MAXBCOORD2 - coord[2];
757 >      break;
758 >    case 2:
759 >      coord[0] = MAXBCOORD2 - coord[0];
760 >      coord[1] = MAXBCOORD2 - coord[1];
761 >      coord[2] = 0;
762 >      break;
763      }
764 +    break;
765 +  }
766 + }
767  
768 <  /* Otherwise return child */
769 <  if(coord[0] > 0.5)
768 > int
769 > bary_child(coord)
770 > BCOORD coord[3];
771 > {
772 >  int i;
773 >
774 >  if(coord[0] > MAXBCOORD4)
775    {
776        /* update bary for child */
777 <      coord[0] = 2.0*coord[0]- 1.0;
778 <      coord[1] *= 2.0;
779 <      coord[2] *= 2.0;
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] > 0.5)
783 >    if(coord[1] > MAXBCOORD4)
784      {
785 <      coord[0] *= 2.0;
786 <      coord[1] = 2.0*coord[1]- 1.0;
787 <      coord[2] *= 2.0;
785 >      coord[0] <<= 1;
786 >      coord[1] = (coord[1] << 1) - MAXBCOORD2;
787 >      coord[2] <<= 1;
788        return(1);
789      }
790      else
791 <      if(coord[2] > 0.5)
791 >      if(coord[2] > MAXBCOORD4)
792        {
793 <        coord[0] *= 2.0;
794 <        coord[1] *= 2.0;
795 <        coord[2] = 2.0*coord[2]- 1.0;
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] = 1.0 - 2.0*coord[0];
801 <           coord[1] = 1.0 - 2.0*coord[1];
802 <           coord[2] = 1.0 - 2.0*coord[2];
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 < max_index(v)
809 < FVECT v;
808 > bary_nth_child(coord,i)
809 > BCOORD coord[3];
810 > int i;
811   {
975 –  double a,b,c;
976 –  int i;
812  
813 <  a = fabs(v[0]);
814 <  b = fabs(v[1]);
815 <  c = fabs(v[2]);
816 <  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
817 <  return(i);
818 < }
819 <
820 <
821 <
822 < /*
823 < * int
824 < * smRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
825 < *
826 < *   Intersect the ray with triangle v0v1v2, return intersection point in r
827 < *
828 < *    Assumes orig,v0,v1,v2 are in spherical coordinates, and orig is
829 < *    unit
830 < */
831 < int
832 < traceRay(orig,dir,v0,v1,v2,r)
833 <  FVECT orig,dir;
834 <  FVECT v0,v1,v2;
1000 <  FVECT r;
1001 < {
1002 <  FVECT n,p[3],d;
1003 <  double pt[3],r_eps;
1004 <  char i;
1005 <  int which;
1006 <
1007 <  /* Find the plane equation for the triangle defined by the edge v0v1 and
1008 <   the view center*/
1009 <  VCROSS(n,v0,v1);
1010 <  /* Intersect the ray with this plane */
1011 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[0]),p[0]);
1012 <  if(i)
1013 <    which = 0;
1014 <  else
1015 <    which = -1;
1016 <
1017 <  VCROSS(n,v1,v2);
1018 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[1]),p[1]);
1019 <  if(i)
1020 <    if((which==-1) || pt[1] < pt[0])
1021 <      which = 1;
1022 <
1023 <  VCROSS(n,v2,v0);
1024 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[2]),p[2]);
1025 <  if(i)
1026 <    if((which==-1) || pt[2] < pt[which])
1027 <      which = 2;
1028 <
1029 <  if(which != -1)
1030 <  {
1031 <      /* Return point that is K*eps along projection of the ray on
1032 <         the sphere to push intersection point p[which] into next cell
1033 <      */
1034 <      normalize(p[which]);
1035 <      /* Calculate the ray perpendicular to the intersection point: approx
1036 <       the direction moved along the sphere a distance of K*epsilon*/
1037 <      r_eps = -DOT(p[which],dir);
1038 <      VSUM(n,dir,p[which],r_eps);
1039 <     /* Calculate the length along ray p[which]-dir needed to move to
1040 <         cause a move along the sphere of k*epsilon
1041 <       */
1042 <       r_eps = DOT(n,dir);
1043 <      VSUM(r,p[which],dir,(20*FTINY)/r_eps);
1044 <      normalize(r);
1045 <      return(TRUE);
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    }
1047 –  else
1048 –  {
1049 –      /* Unable to find intersection: move along ray and try again */
1050 –      r_eps = -DOT(orig,dir);
1051 –      VSUM(n,dir,orig,r_eps);
1052 –      r_eps = DOT(n,dir);
1053 –      VSUM(r,orig,dir,(20*FTINY)/r_eps);
1054 –      normalize(r);
1055 – #ifdef DEBUG
1056 –      eputs("traceRay:Ray does not intersect triangle");
1057 – #endif
1058 –      return(FALSE);
1059 –  }
836   }
1061 –
1062 –
1063 –
1064 –
1065 –
1066 –
837  
838  
839  

Diff Legend

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