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.11 by gwlarson, Sun Jan 10 10:27:47 1999 UTC

# Line 27 | Line 27 | FVECT v1,v2;
27   {
28     return(EQUAL(v1[0],v2[0]) && EQUAL(v1[1],v2[1])&& EQUAL(v1[2],v2[2]));
29   }
30 + #if 0
31 + extern FVECT Norm[500];
32 + extern int Ncnt;
33 + #endif
34  
31 –
35   int
36   convex_angle(v0,v1,v2)
37   FVECT v0,v1,v2;
38   {
39 <    FVECT cp01,cp12,cp;
40 <    
39 >    FVECT cp,cp01,cp12,v10,v02;
40 >    double dp;
41 >
42      /* test sign of (v0Xv1)X(v1Xv2). v1 */
43      VCROSS(cp01,v0,v1);
44      VCROSS(cp12,v1,v2);
45      VCROSS(cp,cp01,cp12);
46 <    if(DOT(cp,v1) < 0)
47 <       return(FALSE);
46 >        
47 >    dp = DOT(cp,v1);
48 > #if 0
49 >    VCOPY(Norm[Ncnt++],cp01);
50 >    VCOPY(Norm[Ncnt++],cp12);
51 >    VCOPY(Norm[Ncnt++],cp);
52 > #endif
53 >    if(ZERO(dp) || dp < 0.0)
54 >      return(FALSE);
55      return(TRUE);
56   }
57  
58   /* calculates the normal of a face contour using Newell's formula. e
59  
60 <               a =  SUMi (yi - yi+1)(zi + zi+1)
60 >               a =  SUMi (yi - yi+1)(zi + zi+1);
61                 b =  SUMi (zi - zi+1)(xi + xi+1)
62                 c =  SUMi (xi - xi+1)(yi + yi+1)
63   */
64   double
65   tri_normal(v0,v1,v2,n,norm)
66   FVECT v0,v1,v2,n;
67 < char norm;
67 > int norm;
68   {
69    double mag;
70  
# Line 63 | Line 74 | char norm;
74    
75    n[1] = (v0[2] - v1[2]) * (v0[0] + v1[0]) +
76             (v1[2] - v2[2]) * (v1[0] + v2[0]) +
77 <           (v2[2] - v0[2]) * (v2[0] + v0[0]);
67 <
77 >            (v2[2] - v0[2]) * (v2[0] + v0[0]);
78    
79    n[2] = (v0[1] + v1[1]) * (v0[0] - v1[0]) +
80           (v1[1] + v2[1]) * (v1[0] - v2[0]) +
# Line 72 | Line 82 | char norm;
82  
83    if(!norm)
84       return(0);
75 –
85    
86    mag = normalize(n);
87  
# Line 80 | Line 89 | char norm;
89   }
90  
91  
92 < tri_plane_equation(v0,v1,v2,n,nd,norm)
93 <   FVECT v0,v1,v2,n;
94 <   double *nd;
95 <   char norm;
92 > tri_plane_equation(v0,v1,v2,peqptr,norm)
93 >   FVECT v0,v1,v2;
94 >   FPEQ *peqptr;
95 >   int norm;
96   {
97 <    tri_normal(v0,v1,v2,n,norm);
97 >    tri_normal(v0,v1,v2,FP_N(*peqptr),norm);
98  
99 <    *nd = -(DOT(n,v0));
99 >    FP_D(*peqptr) = -(DOT(FP_N(*peqptr),v0));
100   }
101  
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 – }
108 –
102   /* From quad_edge-code */
103   int
104   point_in_circle_thru_origin(p,p0,p1)
# Line 122 | Line 115 | FVECT p0,p1;
115  
116      det = -dp0*CROSS_VEC2(p1,p) + dp1*CROSS_VEC2(p0,p) - dp*CROSS_VEC2(p0,p1);
117      
118 <    return (det > 0);
118 >    return (det > (1e-10));
119   }
120  
121  
122 <
122 > double
123   point_on_sphere(ps,p,c)
124   FVECT ps,p,c;
125   {
126 +  double d;
127      VSUB(ps,p,c);    
128 <    normalize(ps);
128 >    d= normalize(ps);
129 >    return(d);
130   }
131  
132  
133 + /* returns TRUE if ray from origin in direction v intersects plane defined
134 +  * by normal plane_n, and plane_d. If plane is not parallel- returns
135 +  * intersection point if r != NULL. If tptr!= NULL returns value of
136 +  * t, if parallel, returns t=FHUGE
137 +  */
138   int
139 < intersect_vector_plane(v,plane_n,plane_d,tptr,r)
140 <   FVECT v,plane_n;
141 <   double plane_d;
139 > intersect_vector_plane(v,peq,tptr,r)
140 >   FVECT v;
141 >   FPEQ peq;
142     double *tptr;
143     FVECT r;
144   {
145 <  double t;
145 >  double t,d;
146    int hit;
147      /*
148        Plane is Ax + By + Cz +D = 0:
# Line 152 | Line 152 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
152      /* line is  l = p1 + (p2-p1)t, p1=origin */
153  
154      /* Solve for t: */
155 <    t =  plane_d/-(DOT(plane_n,v));
156 <    if(t >0 || ZERO(t))
157 <       hit = 1;
158 <    else
159 <       hit = 0;
160 <    r[0] = v[0]*t;
161 <    r[1] = v[1]*t;
162 <    r[2] = v[2]*t;
155 >  d = -(DOT(FP_N(peq),v));
156 >  if(ZERO(d))
157 >  {
158 >      t = FHUGE;
159 >      hit = 0;
160 >  }
161 >  else
162 >  {
163 >      t =  FP_D(peq)/d;
164 >      if(t < 0 )
165 >         hit = 0;
166 >      else
167 >         hit = 1;
168 >      if(r)
169 >         {
170 >             r[0] = v[0]*t;
171 >             r[1] = v[1]*t;
172 >             r[2] = v[2]*t;
173 >         }
174 >  }
175      if(tptr)
176         *tptr = t;
177    return(hit);
178   }
179  
180   int
181 < intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
181 > intersect_ray_plane(orig,dir,peq,pd,r)
182     FVECT orig,dir;
183 <   FVECT plane_n;
172 <   double plane_d;
183 >   FPEQ peq;
184     double *pd;
185     FVECT r;
186   {
187 <  double t;
187 >  double t,d;
188    int hit;
189      /*
190        Plane is Ax + By + Cz +D = 0:
# Line 184 | Line 195 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
195         line is  l = p1 + (p2-p1)t
196       */
197      /* Solve for t: */
198 <    t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
199 <    if(ZERO(t) || t >0)
198 >  d = DOT(FP_N(peq),dir);
199 >  if(ZERO(d))
200 >     return(0);
201 >  t =  -(DOT(FP_N(peq),orig) + FP_D(peq))/d;
202 >
203 >  if(t < 0)
204 >       hit = 0;
205 >    else
206         hit = 1;
207 +
208 +  if(r)
209 +     VSUM(r,orig,dir,t);
210 +
211 +    if(pd)
212 +       *pd = t;
213 +  return(hit);
214 + }
215 +
216 +
217 + int
218 + intersect_ray_oplane(orig,dir,n,pd,r)
219 +   FVECT orig,dir;
220 +   FVECT n;
221 +   double *pd;
222 +   FVECT r;
223 + {
224 +  double t,d;
225 +  int hit;
226 +    /*
227 +      Plane is Ax + By + Cz +D = 0:
228 +      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
229 +    */
230 +     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
231 +         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
232 +       line is  l = p1 + (p2-p1)t
233 +     */
234 +    /* Solve for t: */
235 +    d= DOT(n,dir);
236 +    if(ZERO(d))
237 +       return(0);
238 +    t =  -(DOT(n,orig))/d;
239 +    if(t < 0)
240 +       hit = 0;
241      else
242 +       hit = 1;
243 +
244 +  if(r)
245 +     VSUM(r,orig,dir,t);
246 +
247 +    if(pd)
248 +       *pd = t;
249 +  return(hit);
250 + }
251 +
252 + /* Assumption: know crosses plane:dont need to check for 'on' case */
253 + intersect_edge_coord_plane(v0,v1,w,r)
254 + FVECT v0,v1;
255 + int w;
256 + FVECT r;
257 + {
258 +  FVECT dv;
259 +  int wnext;
260 +  double t;
261 +
262 +  VSUB(dv,v1,v0);
263 +  t = -v0[w]/dv[w];
264 +  r[w] = 0.0;
265 +  wnext = (w+1)%3;
266 +  r[wnext] = v0[wnext] + dv[wnext]*t;
267 +  wnext = (w+2)%3;
268 +  r[wnext] = v0[wnext] + dv[wnext]*t;
269 + }
270 +
271 + int
272 + intersect_edge_plane(e0,e1,peq,pd,r)
273 +   FVECT e0,e1;
274 +   FPEQ peq;
275 +   double *pd;
276 +   FVECT r;
277 + {
278 +  double t;
279 +  int hit;
280 +  FVECT d;
281 +  /*
282 +      Plane is Ax + By + Cz +D = 0:
283 +      plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
284 +    */
285 +     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
286 +         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
287 +       line is  l = p1 + (p2-p1)t
288 +     */
289 +    /* Solve for t: */
290 +  VSUB(d,e1,e0);
291 +  t =  -(DOT(FP_N(peq),e0) + FP_D(peq))/(DOT(FP_N(peq),d));
292 +    if(t < 0)
293         hit = 0;
294 +    else
295 +       hit = 1;
296  
297 <  VSUM(r,orig,dir,t);
297 >  VSUM(r,e0,d,t);
298  
299      if(pd)
300         *pd = t;
# Line 203 | Line 307 | point_in_cone(p,p0,p1,p2)
307   FVECT p;
308   FVECT p0,p1,p2;
309   {
206 –    FVECT n;
310      FVECT np,x_axis,y_axis;
311 <    double d1,d2,d;
311 >    double d1,d2;
312 >    FPEQ peq;
313      
314      /* Find the equation of the circle defined by the intersection
315         of the cone with the plane defined by p1,p2,p3- project p into
316         that plane and do an in-circle test in the plane
317       */
318      
319 <    /* find the equation of the plane defined by p1-p3 */
320 <    tri_plane_equation(p0,p1,p2,n,&d,FALSE);
319 >    /* find the equation of the plane defined by p0-p2 */
320 >    tri_plane_equation(p0,p1,p2,&peq,FALSE);
321  
322      /* define a coordinate system on the plane: the x axis is in
323         the direction of np2-np1, and the y axis is calculated from
324         n cross x-axis
325       */
326      /* Project p onto the plane */
327 <    if(!intersect_vector_plane(p,n,d,NULL,np))
327 >    /* NOTE: check this: does sideness check?*/
328 >    if(!intersect_vector_plane(p,peq,NULL,np))
329          return(FALSE);
330  
331 <    /* create coordinate system on  plane: p2-p1 defines the x_axis*/
331 >    /* create coordinate system on  plane: p1-p0 defines the x_axis*/
332      VSUB(x_axis,p1,p0);
333      normalize(x_axis);
334      /* The y axis is  */
335 <    VCROSS(y_axis,n,x_axis);
335 >    VCROSS(y_axis,FP_N(peq),x_axis);
336      normalize(y_axis);
337  
338      VSUB(p1,p1,p0);
339      VSUB(p2,p2,p0);
340      VSUB(np,np,p0);
341      
342 <    p1[0] = VLEN(p1);
342 >    p1[0] = DOT(p1,x_axis);
343      p1[1] = 0;
344  
345      d1 = DOT(p2,x_axis);
# Line 252 | Line 357 | FVECT p0,p1,p2;
357   }
358  
359   int
360 < test_point_against_spherical_tri(v0,v1,v2,p,n,nset,which,sides)
360 > point_set_in_stri(v0,v1,v2,p,n,nset,sides)
361   FVECT v0,v1,v2,p;
362   FVECT n[3];
363 < char *nset;
364 < char *which;
260 < char sides[3];
363 > int *nset;
364 > int sides[3];
365  
366   {
367 <    float d;
264 <
367 >    double d;
368      /* Find the normal to the triangle ORIGIN:v0:v1 */
369      if(!NTH_BIT(*nset,0))
370      {
371 <        VCROSS(n[0],v1,v0);
371 >        VCROSS(n[0],v0,v1);
372          SET_NTH_BIT(*nset,0);
373      }
374      /* Test the point for sidedness */
375      d  = DOT(n[0],p);
376  
377 <    if(ZERO(d))
378 <       sides[0] = GT_EDGE;
379 <    else
380 <       if(d > 0)
381 <      {
279 <          sides[0] =  GT_OUT;
280 <          sides[1] = sides[2] = GT_INVALID;
281 <          return(FALSE);
377 >    if(d > 0.0)
378 >     {
379 >       sides[0] =  GT_OUT;
380 >       sides[1] = sides[2] = GT_INVALID;
381 >       return(FALSE);
382        }
383      else
384         sides[0] = GT_INTERIOR;
# Line 286 | Line 386 | char sides[3];
386      /* Test next edge */
387      if(!NTH_BIT(*nset,1))
388      {
389 <        VCROSS(n[1],v2,v1);
389 >        VCROSS(n[1],v1,v2);
390          SET_NTH_BIT(*nset,1);
391      }
392      /* Test the point for sidedness */
393      d  = DOT(n[1],p);
394 <    if(ZERO(d))
394 >    if(d > 0.0)
395      {
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 –    {
396          sides[1] = GT_OUT;
397          sides[2] = GT_INVALID;
398          return(FALSE);
# Line 313 | Line 402 | char sides[3];
402      /* Test next edge */
403      if(!NTH_BIT(*nset,2))
404      {
405 <
317 <        VCROSS(n[2],v0,v2);
405 >        VCROSS(n[2],v2,v0);
406          SET_NTH_BIT(*nset,2);
407      }
408      /* Test the point for sidedness */
409      d  = DOT(n[2],p);
410 <    if(ZERO(d))
410 >    if(d > 0.0)
411      {
412 <        sides[2] = GT_EDGE;
413 <
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 <           }
412 >      sides[2] = GT_OUT;
413 >      return(FALSE);
414      }
345 –    else if(d > 0)
346 –      {
347 –        sides[2] = GT_OUT;
348 –        return(FALSE);
349 –      }
350 –    /* If on edge */
415      else
416         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 –    }
417      /* Must be interior to the pyramid */
418      return(GT_INTERIOR);
419   }
420  
421  
422  
423 <
423 >
424   int
425 < test_single_point_against_spherical_tri(v0,v1,v2,p,which)
425 > point_in_stri(v0,v1,v2,p)
426   FVECT v0,v1,v2,p;
376 – char *which;
427   {
428 <    float d;
428 >    double d;
429      FVECT n;  
380 –    char sides[3];
430  
431 <    /* 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);
431 >    VCROSS(n,v0,v1);
432      /* Test the point for sidedness */
433      d  = DOT(n,p);
434 <    if(ZERO(d))
435 <       sides[0] = GT_EDGE;
436 <    else
404 <       if(d > 0)
405 <          return(FALSE);
406 <       else
407 <          sides[0] = GT_INTERIOR;
434 >    if(d > 0.0)
435 >      return(FALSE);
436 >
437      /* Test next edge */
438 <    VCROSS(n,v2,v1);
438 >    VCROSS(n,v1,v2);
439      /* Test the point for sidedness */
440      d  = DOT(n,p);
441 <    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)
441 >    if(d > 0.0)
442         return(FALSE);
424 –    else
425 –       sides[1] = GT_INTERIOR;
443  
444      /* Test next edge */
445 <    VCROSS(n,v0,v2);
445 >    VCROSS(n,v2,v0);
446      /* Test the point for sidedness */
447      d  = DOT(n,p);
448 <    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)
448 >    if(d > 0.0)
449         return(FALSE);
450      /* Must be interior to the pyramid */
451 <    return(GT_FACE);
451 >    return(GT_INTERIOR);
452   }
453  
454   int
455 < test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
455 > vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
456   FVECT t0,t1,t2,p0,p1,p2;
457 < char *nset;
457 > int *nset;
458   FVECT n[3];
459   FVECT avg;
460 < char pt_sides[3][3];
460 > int pt_sides[3][3];
461  
462   {
463 <    char below_plane[3],on_edge,test;
469 <    char which;
463 >    int below_plane[3],test;
464  
465      SUM_3VEC3(avg,t0,t1,t2);
472 –    on_edge = 0;
466      *nset = 0;
467      /* Test vertex v[i] against triangle j*/
468      /* Check if v[i] lies below plane defined by avg of 3 vectors
# Line 477 | Line 470 | char pt_sides[3][3];
470         */
471  
472      /* test point 0 */
473 <    if(DOT(avg,p0) < 0)
473 >    if(DOT(avg,p0) < 0.0)
474        below_plane[0] = 1;
475      else
476 <      below_plane[0]=0;
476 >      below_plane[0] = 0;
477      /* Test if b[i] lies in or on triangle a */
478 <    test = test_point_against_spherical_tri(t0,t1,t2,p0,
486 <                                                 n,nset,&which,pt_sides[0]);
478 >    test = point_set_in_stri(t0,t1,t2,p0,n,nset,pt_sides[0]);
479      /* If pts[i] is interior: done */
480      if(!below_plane[0])
481        {
482          if(test == GT_INTERIOR)
483            return(TRUE);
492 –        /* Remember if b[i] fell on one of the 3 defining planes */
493 –        if(test)
494 –          on_edge++;
484        }
485      /* Now test point 1*/
486  
487 <    if(DOT(avg,p1) < 0)
487 >    if(DOT(avg,p1) < 0.0)
488        below_plane[1] = 1;
489      else
490        below_plane[1]=0;
491      /* Test if b[i] lies in or on triangle a */
492 <    test = test_point_against_spherical_tri(t0,t1,t2,p1,
504 <                                                 n,nset,&which,pt_sides[1]);
492 >    test = point_set_in_stri(t0,t1,t2,p1,n,nset,pt_sides[1]);
493      /* If pts[i] is interior: done */
494      if(!below_plane[1])
495      {
496        if(test == GT_INTERIOR)
497          return(TRUE);
510 –      /* Remember if b[i] fell on one of the 3 defining planes */
511 –      if(test)
512 –        on_edge++;
498      }
499      
500      /* Now test point 2 */
501 <    if(DOT(avg,p2) < 0)
501 >    if(DOT(avg,p2) < 0.0)
502        below_plane[2] = 1;
503      else
504 <      below_plane[2]=0;
504 >      below_plane[2] = 0;
505          /* Test if b[i] lies in or on triangle a */
506 <    test = test_point_against_spherical_tri(t0,t1,t2,p2,
522 <                                                 n,nset,&which,pt_sides[2]);
506 >    test = point_set_in_stri(t0,t1,t2,p2,n,nset,pt_sides[2]);
507  
508      /* If pts[i] is interior: done */
509      if(!below_plane[2])
510        {
511          if(test == GT_INTERIOR)
512            return(TRUE);
529 –        /* Remember if b[i] fell on one of the 3 defining planes */
530 –        if(test)
531 –          on_edge++;
513        }
514  
515      /* If all three points below separating plane: trivial reject */
516      if(below_plane[0] && below_plane[1] && below_plane[2])
517         return(FALSE);
537 –    /* Accept if all points lie on a triangle vertex/edge edge- accept*/
538 –    if(on_edge == 3)
539 –       return(TRUE);
518      /* Now check vertices in a against triangle b */
519      return(FALSE);
520   }
# Line 544 | Line 522 | char pt_sides[3][3];
522  
523   set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
524     FVECT t0,t1,t2,p0,p1,p2;
525 <   char test[3];
526 <   char sides[3][3];
527 <   char nset;
525 >   int test[3];
526 >   int sides[3][3];
527 >   int nset;
528     FVECT n[3];
529   {
530 <    char t;
530 >    int t;
531      double d;
532  
533      
# Line 558 | Line 536 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
536      if(sides[0][0] == GT_INVALID)
537      {
538        if(!NTH_BIT(nset,0))
539 <        VCROSS(n[0],t1,t0);
539 >        VCROSS(n[0],t0,t1);
540        /* Test the point for sidedness */
541        d  = DOT(n[0],p0);
542 <      if(d >= 0)
542 >      if(d >= 0.0)
543          SET_NTH_BIT(test[0],0);
544      }
545      else
# Line 571 | Line 549 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
549      if(sides[0][1] == GT_INVALID)
550      {
551        if(!NTH_BIT(nset,1))
552 <        VCROSS(n[1],t2,t1);
552 >        VCROSS(n[1],t1,t2);
553          /* Test the point for sidedness */
554          d  = DOT(n[1],p0);
555 <        if(d >= 0)
555 >        if(d >= 0.0)
556            SET_NTH_BIT(test[0],1);
557      }
558      else
# Line 584 | Line 562 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
562      if(sides[0][2] == GT_INVALID)
563      {
564        if(!NTH_BIT(nset,2))
565 <        VCROSS(n[2],t0,t2);
565 >        VCROSS(n[2],t2,t0);
566        /* Test the point for sidedness */
567        d  = DOT(n[2],p0);
568 <      if(d >= 0)
568 >      if(d >= 0.0)
569          SET_NTH_BIT(test[0],2);
570      }
571      else
# Line 600 | Line 578 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
578      if(sides[1][0] == GT_INVALID)
579      {
580        if(!NTH_BIT(nset,0))
581 <        VCROSS(n[0],t1,t0);
581 >        VCROSS(n[0],t0,t1);
582        /* Test the point for sidedness */
583        d  = DOT(n[0],p1);
584 <      if(d >= 0)
584 >      if(d >= 0.0)
585          SET_NTH_BIT(test[1],0);
586      }
587      else
# Line 614 | Line 592 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
592      if(sides[1][1] == GT_INVALID)
593      {
594        if(!NTH_BIT(nset,1))
595 <        VCROSS(n[1],t2,t1);
595 >        VCROSS(n[1],t1,t2);
596        /* Test the point for sidedness */
597        d  = DOT(n[1],p1);
598 <      if(d >= 0)
598 >      if(d >= 0.0)
599          SET_NTH_BIT(test[1],1);
600      }
601      else
# Line 628 | Line 606 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
606      if(sides[1][2] == GT_INVALID)
607      {
608        if(!NTH_BIT(nset,2))
609 <        VCROSS(n[2],t0,t2);
609 >        VCROSS(n[2],t2,t0);
610        /* Test the point for sidedness */
611        d  = DOT(n[2],p1);
612 <      if(d >= 0)
612 >      if(d >= 0.0)
613          SET_NTH_BIT(test[1],2);
614      }
615      else
# Line 644 | Line 622 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
622      if(sides[2][0] == GT_INVALID)
623      {
624        if(!NTH_BIT(nset,0))
625 <        VCROSS(n[0],t1,t0);
625 >        VCROSS(n[0],t0,t1);
626        /* Test the point for sidedness */
627        d  = DOT(n[0],p2);
628 <      if(d >= 0)
628 >      if(d >= 0.0)
629          SET_NTH_BIT(test[2],0);
630      }
631      else
# Line 657 | Line 635 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
635      if(sides[2][1] == GT_INVALID)
636      {
637        if(!NTH_BIT(nset,1))
638 <        VCROSS(n[1],t2,t1);
638 >        VCROSS(n[1],t1,t2);
639        /* Test the point for sidedness */
640        d  = DOT(n[1],p2);
641 <      if(d >= 0)
641 >      if(d >= 0.0)
642          SET_NTH_BIT(test[2],1);
643      }
644      else
# Line 670 | Line 648 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
648      if(sides[2][2] == GT_INVALID)
649      {
650        if(!NTH_BIT(nset,2))
651 <        VCROSS(n[2],t0,t2);
651 >        VCROSS(n[2],t2,t0);
652        /* Test the point for sidedness */
653        d  = DOT(n[2],p2);
654 <      if(d >= 0)
654 >      if(d >= 0.0)
655          SET_NTH_BIT(test[2],2);
656      }
657      else
# Line 683 | Line 661 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
661  
662  
663   int
664 < spherical_tri_intersect(a1,a2,a3,b1,b2,b3)
664 > stri_intersect(a1,a2,a3,b1,b2,b3)
665   FVECT a1,a2,a3,b1,b2,b3;
666   {
667 <  char which,test,n_set[2];
668 <  char sides[2][3][3],i,j,inext,jnext;
669 <  char tests[2][3];
670 <  FVECT n[2][3],p,avg[2];
667 >  int which,test,n_set[2];
668 >  int sides[2][3][3],i,j,inext,jnext;
669 >  int tests[2][3];
670 >  FVECT n[2][3],p,avg[2],t1,t2,t3;
671  
672 + #ifdef DEBUG
673 +    tri_normal(b1,b2,b3,p,FALSE);
674 +    if(DOT(p,b1) > 0)
675 +      {
676 +       VCOPY(t3,b1);
677 +       VCOPY(t2,b2);
678 +       VCOPY(t1,b3);
679 +      }
680 +    else
681 + #endif
682 +      {
683 +       VCOPY(t1,b1);
684 +       VCOPY(t2,b2);
685 +       VCOPY(t3,b3);
686 +      }
687 +
688    /* Test the vertices of triangle a against the pyramid formed by triangle
689       b and the origin. If any vertex of a is interior to triangle b, or
690       if all 3 vertices of a are ON the edges of b,return TRUE. Remember
691       the results of the edge normal and sidedness tests for later.
692     */
693 < if(test_vertices_for_tri_inclusion(a1,a2,a3,b1,b2,b3,
700 <                                    &(n_set[0]),n[0],avg[0],sides[1]))
693 > if(vertices_in_stri(a1,a2,a3,t1,t2,t3,&(n_set[0]),n[0],avg[0],sides[1]))
694       return(TRUE);
695    
696 < if(test_vertices_for_tri_inclusion(b1,b2,b3,a1,a2,a3,
704 <                                    &(n_set[1]),n[1],avg[1],sides[0]))
696 > if(vertices_in_stri(t1,t2,t3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
697       return(TRUE);
698  
699  
700 <  set_sidedness_tests(b1,b2,b3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
700 >  set_sidedness_tests(t1,t2,t3,a1,a2,a3,tests[0],sides[0],n_set[1],n[1]);
701    if(tests[0][0]&tests[0][1]&tests[0][2])
702      return(FALSE);
703  
704 <  set_sidedness_tests(a1,a2,a3,b1,b2,b3,tests[1],sides[1],n_set[0],n[0]);
704 >  set_sidedness_tests(a1,a2,a3,t1,t2,t3,tests[1],sides[1],n_set[0],n[0]);
705    if(tests[1][0]&tests[1][1]&tests[1][2])
706      return(FALSE);
707  
# Line 749 | Line 741 | FVECT a1,a2,a3,b1,b2,b3;
741   }
742  
743   int
744 < ray_intersect_tri(orig,dir,v0,v1,v2,pt,wptr)
744 > ray_intersect_tri(orig,dir,v0,v1,v2,pt)
745   FVECT orig,dir;
746   FVECT v0,v1,v2;
747   FVECT pt;
756 – char *wptr;
748   {
749 <  FVECT p0,p1,p2,p,n;
750 <  char type,which;
751 <  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);
749 >  FVECT p0,p1,p2,p;
750 >  FPEQ peq;
751 >  int type;
752  
753 <  if(type)
753 >  VSUB(p0,v0,orig);
754 >  VSUB(p1,v1,orig);
755 >  VSUB(p2,v2,orig);
756 >
757 >  if(point_in_stri(p0,p1,p2,dir))
758    {
759        /* Intersect the ray with the triangle plane */
760 <      tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
761 <      intersect_ray_plane(orig,dir,n,pd,NULL,pt);        
760 >      tri_plane_equation(v0,v1,v2,&peq,FALSE);
761 >      return(intersect_ray_plane(orig,dir,peq,NULL,pt));
762    }
763 <  if(wptr)
774 <    *wptr = which;
775 <
776 <  return(type);
763 >  return(FALSE);
764   }
765  
766  
# Line 832 | Line 819 | FVECT fnear[4],ffar[4];
819      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
820   }
821  
822 + int
823 + max_index(v,r)
824 + FVECT v;
825 + double *r;
826 + {
827 +  double p[3];
828 +  int i;
829  
830 +  p[0] = fabs(v[0]);
831 +  p[1] = fabs(v[1]);
832 +  p[2] = fabs(v[2]);
833 +  i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);  
834 +  if(r)
835 +    *r = p[i];
836 +  return(i);
837 + }
838  
839 + int
840 + closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
841 + FVECT p0,p1,p2,p;
842 + int p0id,p1id,p2id;
843 + {
844 +    double d,d1;
845 +    int i;
846 +    
847 +    d =  DIST_SQ(p,p0);
848 +    d1 = DIST_SQ(p,p1);
849 +    if(d < d1)
850 +    {
851 +      d1 = DIST_SQ(p,p2);
852 +      i = (d1 < d)?p2id:p0id;
853 +    }
854 +    else
855 +    {
856 +      d = DIST_SQ(p,p2);
857 +      i = (d < d1)? p2id:p1id;
858 +    }
859 +    return(i);
860 + }
861  
862 +
863   int
864 < spherical_polygon_edge_intersect(a0,a1,b0,b1)
864 > sedge_intersect(a0,a1,b0,b1)
865   FVECT a0,a1,b0,b1;
866   {
867      FVECT na,nb,avga,avgb,p;
# Line 882 | Line 907 | FVECT a0,a1,b0,b1;
907   }
908  
909  
885 –
910   /* Find the normalized barycentric coordinates of p relative to
911   * triangle v0,v1,v2. Return result in coord
912   */
# Line 896 | Line 920 | double coord[3];
920    a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
921    coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
922    coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
923 <  coord[2]  = 1.0 - coord[0] - coord[1];
923 >  coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
924  
925   }
926  
903 – int
904 – bary2d_child(coord)
905 – double coord[3];
906 – {
907 –  int i;
927  
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);
928  
929 <  /* Check if one of the new vertices: for all return center child */
930 <  if(ZERO(coord[0]) && EQUAL(coord[1],0.5))
929 >
930 > bary_parent(coord,i)
931 > BCOORD coord[3];
932 > int i;
933 > {
934 >  switch(i) {
935 >  case 0:
936 >    /* update bary for child */
937 >    coord[0] = (coord[0] >> 1) + MAXBCOORD4;
938 >    coord[1] >>= 1;
939 >    coord[2] >>= 1;
940 >    break;
941 >  case 1:
942 >    coord[0] >>= 1;
943 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD4;
944 >    coord[2] >>= 1;
945 >    break;
946 >    
947 >  case 2:
948 >    coord[0] >>= 1;
949 >    coord[1] >>= 1;
950 >    coord[2] = (coord[2] >> 1) + MAXBCOORD4;
951 >    break;
952 >    
953 >  case 3:
954 >    coord[0] = MAXBCOORD4 - (coord[0] >> 1);
955 >    coord[1] = MAXBCOORD4 - (coord[1] >> 1);
956 >    coord[2] = MAXBCOORD4 - (coord[2] >> 1);
957 >    break;
958 > #ifdef DEBUG
959 >  default:
960 >    eputs("bary_parent():Invalid child\n");
961 >    break;
962 > #endif
963 >  }
964 > }
965 >
966 > bary_from_child(coord,child,next)
967 > BCOORD coord[3];
968 > int child,next;
969 > {
970 > #ifdef DEBUG
971 >  if(child <0 || child > 3)
972    {
973 <    coord[0] = 1.0f;
974 <    coord[1] = 0.0f;
919 <    coord[2] = 0.0f;
920 <    return(3);
973 >    eputs("bary_from_child():Invalid child\n");
974 >    return;
975    }
976 <  if(ZERO(coord[1]) && EQUAL(coord[0],0.5))
976 >  if(next <0 || next > 3)
977    {
978 <    coord[0] = 0.0f;
979 <    coord[1] = 1.0f;
926 <    coord[2] = 0.0f;
927 <    return(3);
978 >    eputs("bary_from_child():Invalid next\n");
979 >    return;
980    }
981 <  if(ZERO(coord[2]) && EQUAL(coord[0],0.5))
982 <    {
983 <      coord[0] = 0.0f;
984 <      coord[1] = 0.0f;
985 <      coord[2] = 1.0f;
986 <      return(3);
981 > #endif
982 >  if(next == child)
983 >    return;
984 >
985 >  switch(child){
986 >  case 0:
987 >      coord[0] = 0;
988 >      coord[1] = MAXBCOORD2 - coord[1];
989 >      coord[2] = MAXBCOORD2 - coord[2];
990 >      break;
991 >  case 1:
992 >      coord[0] = MAXBCOORD2 - coord[0];
993 >      coord[1] = 0;
994 >      coord[2] = MAXBCOORD2 - coord[2];
995 >      break;
996 >  case 2:
997 >      coord[0] = MAXBCOORD2 - coord[0];
998 >      coord[1] = MAXBCOORD2 - coord[1];
999 >      coord[2] = 0;
1000 >    break;
1001 >  case 3:
1002 >    switch(next){
1003 >    case 0:
1004 >      coord[0] = 0;
1005 >      coord[1] = MAXBCOORD2 - coord[1];
1006 >      coord[2] = MAXBCOORD2 - coord[2];
1007 >      break;
1008 >    case 1:
1009 >      coord[0] = MAXBCOORD2 - coord[0];
1010 >      coord[1] = 0;
1011 >      coord[2] = MAXBCOORD2 - coord[2];
1012 >      break;
1013 >    case 2:
1014 >      coord[0] = MAXBCOORD2 - coord[0];
1015 >      coord[1] = MAXBCOORD2 - coord[1];
1016 >      coord[2] = 0;
1017 >      break;
1018      }
1019 +    break;
1020 +  }
1021 + }
1022  
1023 <  /* Otherwise return child */
1024 <  if(coord[0] > 0.5)
1023 > int
1024 > bary_child(coord)
1025 > BCOORD coord[3];
1026 > {
1027 >  int i;
1028 >
1029 >  if(coord[0] > MAXBCOORD4)
1030    {
1031        /* update bary for child */
1032 <      coord[0] = 2.0*coord[0]- 1.0;
1033 <      coord[1] *= 2.0;
1034 <      coord[2] *= 2.0;
1032 >      coord[0] = (coord[0]<< 1) - MAXBCOORD2;
1033 >      coord[1] <<= 1;
1034 >      coord[2] <<= 1;
1035        return(0);
1036    }
1037    else
1038 <    if(coord[1] > 0.5)
1038 >    if(coord[1] > MAXBCOORD4)
1039      {
1040 <      coord[0] *= 2.0;
1041 <      coord[1] = 2.0*coord[1]- 1.0;
1042 <      coord[2] *= 2.0;
1040 >      coord[0] <<= 1;
1041 >      coord[1] = (coord[1] << 1) - MAXBCOORD2;
1042 >      coord[2] <<= 1;
1043        return(1);
1044      }
1045      else
1046 <      if(coord[2] > 0.5)
1046 >      if(coord[2] > MAXBCOORD4)
1047        {
1048 <        coord[0] *= 2.0;
1049 <        coord[1] *= 2.0;
1050 <        coord[2] = 2.0*coord[2]- 1.0;
1048 >        coord[0] <<= 1;
1049 >        coord[1] <<= 1;
1050 >        coord[2] = (coord[2] << 1) - MAXBCOORD2;
1051          return(2);
1052        }
1053        else
1054           {
1055 <           coord[0] = 1.0 - 2.0*coord[0];
1056 <           coord[1] = 1.0 - 2.0*coord[1];
1057 <           coord[2] = 1.0 - 2.0*coord[2];
1055 >           coord[0] = MAXBCOORD2 - (coord[0] << 1);
1056 >           coord[1] = MAXBCOORD2 - (coord[1] << 1);
1057 >           coord[2] = MAXBCOORD2 - (coord[2] << 1);
1058             return(3);
1059           }
1060   }
1061  
1062   int
1063 < max_index(v)
1064 < FVECT v;
1063 > bary_nth_child(coord,i)
1064 > BCOORD coord[3];
1065 > int i;
1066   {
975 –  double a,b,c;
976 –  int i;
1067  
1068 <  a = fabs(v[0]);
1069 <  b = fabs(v[1]);
1070 <  c = fabs(v[2]);
1071 <  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
1072 <  return(i);
1073 < }
1074 <
1075 <
1076 <
1077 < /*
1078 < * int
1079 < * smRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
1080 < *
1081 < *   Intersect the ray with triangle v0v1v2, return intersection point in r
1082 < *
1083 < *    Assumes orig,v0,v1,v2 are in spherical coordinates, and orig is
1084 < *    unit
1085 < */
1086 < int
1087 < traceRay(orig,dir,v0,v1,v2,r)
1088 <  FVECT orig,dir;
1089 <  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);
1068 >  switch(i){
1069 >  case 0:
1070 >    /* update bary for child */
1071 >    coord[0] = (coord[0]<< 1) - MAXBCOORD2;
1072 >    coord[1] <<= 1;
1073 >    coord[2] <<= 1;
1074 >    break;
1075 >  case 1:
1076 >    coord[0] <<= 1;
1077 >    coord[1] = (coord[1] << 1) - MAXBCOORD2;
1078 >    coord[2] <<= 1;
1079 >    break;
1080 >  case 2:
1081 >    coord[0] <<= 1;
1082 >    coord[1] <<= 1;
1083 >    coord[2] = (coord[2] << 1) - MAXBCOORD2;
1084 >    break;
1085 >  case 3:
1086 >    coord[0] = MAXBCOORD2 - (coord[0] << 1);
1087 >    coord[1] = MAXBCOORD2 - (coord[1] << 1);
1088 >    coord[2] = MAXBCOORD2 - (coord[2] << 1);
1089 >    break;
1090    }
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 –  }
1091   }
1061 –
1062 –
1063 –
1064 –
1065 –
1066 –
1092  
1093  
1094  

Diff Legend

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