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.6 by gwlarson, Wed Sep 16 18:16:28 1998 UTC

# Line 39 | Line 39 | FVECT v0,v1,v2;
39      VCROSS(cp01,v0,v1);
40      VCROSS(cp12,v1,v2);
41      VCROSS(cp,cp01,cp12);
42 <    if(DOT(cp,v1) < 0)
42 >    if(DOT(cp,v1) < 0.0)
43         return(FALSE);
44      return(TRUE);
45   }
# Line 53 | Line 53 | FVECT v0,v1,v2;
53   double
54   tri_normal(v0,v1,v2,n,norm)
55   FVECT v0,v1,v2,n;
56 < char norm;
56 > int norm;
57   {
58    double mag;
59  
# Line 83 | Line 83 | char norm;
83   tri_plane_equation(v0,v1,v2,n,nd,norm)
84     FVECT v0,v1,v2,n;
85     double *nd;
86 <   char norm;
86 >   int norm;
87   {
88      tri_normal(v0,v1,v2,n,norm);
89  
90      *nd = -(DOT(n,v0));
91   }
92  
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 –
93   /* From quad_edge-code */
94   int
95   point_in_circle_thru_origin(p,p0,p1)
# Line 135 | Line 119 | FVECT ps,p,c;
119   }
120  
121  
122 + /* returns TRUE if ray from origin in direction v intersects plane defined
123 +  * by normal plane_n, and plane_d. If plane is not parallel- returns
124 +  * intersection point if r != NULL. If tptr!= NULL returns value of
125 +  * t, if parallel, returns t=FHUGE
126 +  */
127   int
128   intersect_vector_plane(v,plane_n,plane_d,tptr,r)
129     FVECT v,plane_n;
# Line 142 | Line 131 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
131     double *tptr;
132     FVECT r;
133   {
134 <  double t;
134 >  double t,d;
135    int hit;
136      /*
137        Plane is Ax + By + Cz +D = 0:
# Line 152 | Line 141 | intersect_vector_plane(v,plane_n,plane_d,tptr,r)
141      /* line is  l = p1 + (p2-p1)t, p1=origin */
142  
143      /* Solve for t: */
144 <    t =  plane_d/-(DOT(plane_n,v));
145 <    if(t >0 || ZERO(t))
146 <       hit = 1;
147 <    else
148 <       hit = 0;
149 <    r[0] = v[0]*t;
150 <    r[1] = v[1]*t;
151 <    r[2] = v[2]*t;
144 >  d = -(DOT(plane_n,v));
145 >  if(ZERO(d))
146 >  {
147 >      t = FHUGE;
148 >      hit = 0;
149 >  }
150 >  else
151 >  {
152 >      t =  plane_d/d;
153 >      if(t < 0 )
154 >         hit = 0;
155 >      else
156 >         hit = 1;
157 >      if(r)
158 >         {
159 >             r[0] = v[0]*t;
160 >             r[1] = v[1]*t;
161 >             r[2] = v[2]*t;
162 >         }
163 >  }
164      if(tptr)
165         *tptr = t;
166    return(hit);
# Line 185 | Line 186 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
186       */
187      /* Solve for t: */
188      t =  -(DOT(plane_n,orig) + plane_d)/(DOT(plane_n,dir));
189 <    if(ZERO(t) || t >0)
190 <       hit = 1;
189 >    if(t < 0)
190 >       hit = 0;
191      else
192 +       hit = 1;
193 +
194 +  if(r)
195 +     VSUM(r,orig,dir,t);
196 +
197 +    if(pd)
198 +       *pd = t;
199 +  return(hit);
200 + }
201 +
202 +
203 + int
204 + intersect_edge_plane(e0,e1,plane_n,plane_d,pd,r)
205 +   FVECT e0,e1;
206 +   FVECT plane_n;
207 +   double plane_d;
208 +   double *pd;
209 +   FVECT r;
210 + {
211 +  double t;
212 +  int hit;
213 +  FVECT d;
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 +    /* Solve for t: */
223 +  VSUB(d,e1,e0);
224 +  t =  -(DOT(plane_n,e0) + plane_d)/(DOT(plane_n,d));
225 +    if(t < 0)
226         hit = 0;
227 +    else
228 +       hit = 1;
229  
230 <  VSUM(r,orig,dir,t);
230 >  VSUM(r,e0,d,t);
231  
232      if(pd)
233         *pd = t;
# Line 220 | Line 257 | FVECT p0,p1,p2;
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);
263  
# Line 252 | Line 290 | FVECT p0,p1,p2;
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      {
# Line 271 | Line 307 | char sides[3];
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 291 | Line 324 | char sides[3];
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      {
316 –
338          VCROSS(n[2],v0,v2);
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   }
# Line 371 | Line 355 | char sides[3];
355  
356  
357   int
358 < test_single_point_against_spherical_tri(v0,v1,v2,p,which)
358 > point_in_stri(v0,v1,v2,p)
359   FVECT v0,v1,v2,p;
376 – char *which;
360   {
361 <    float d;
361 >    double d;
362      FVECT n;  
380 –    char sides[3];
363  
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 –    }
364      VCROSS(n,v1,v0);
365      /* Test the point for sidedness */
366      d  = DOT(n,p);
367 <    if(ZERO(d))
368 <       sides[0] = GT_EDGE;
369 <    else
404 <       if(d > 0)
405 <          return(FALSE);
406 <       else
407 <          sides[0] = GT_INTERIOR;
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(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)
374 >    if(d > 0.0)
375         return(FALSE);
424 –    else
425 –       sides[1] = GT_INTERIOR;
376  
377      /* Test next edge */
378      VCROSS(n,v0,v2);
379      /* Test the point for sidedness */
380      d  = DOT(n,p);
381 <    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)
381 >    if(d > 0.0)
382         return(FALSE);
383      /* Must be interior to the pyramid */
384 <    return(GT_FACE);
384 >    return(GT_INTERIOR);
385   }
386  
387   int
388 < test_vertices_for_tri_inclusion(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
388 > vertices_in_stri(t0,t1,t2,p0,p1,p2,nset,n,avg,pt_sides)
389   FVECT t0,t1,t2,p0,p1,p2;
390 < char *nset;
390 > int *nset;
391   FVECT n[3];
392   FVECT avg;
393 < char pt_sides[3][3];
393 > int pt_sides[3][3];
394  
395   {
396 <    char below_plane[3],on_edge,test;
469 <    char which;
396 >    int below_plane[3],test;
397  
398      SUM_3VEC3(avg,t0,t1,t2);
472 –    on_edge = 0;
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
# Line 477 | Line 403 | char pt_sides[3][3];
403         */
404  
405      /* test point 0 */
406 <    if(DOT(avg,p0) < 0)
406 >    if(DOT(avg,p0) < 0.0)
407        below_plane[0] = 1;
408      else
409 <      below_plane[0]=0;
409 >      below_plane[0] = 0;
410      /* Test if b[i] lies in or on triangle a */
411 <    test = test_point_against_spherical_tri(t0,t1,t2,p0,
486 <                                                 n,nset,&which,pt_sides[0]);
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);
492 –        /* Remember if b[i] fell on one of the 3 defining planes */
493 –        if(test)
494 –          on_edge++;
417        }
418      /* Now test point 1*/
419  
420 <    if(DOT(avg,p1) < 0)
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 = test_point_against_spherical_tri(t0,t1,t2,p1,
504 <                                                 n,nset,&which,pt_sides[1]);
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);
510 –      /* Remember if b[i] fell on one of the 3 defining planes */
511 –      if(test)
512 –        on_edge++;
431      }
432      
433      /* Now test point 2 */
434 <    if(DOT(avg,p2) < 0)
434 >    if(DOT(avg,p2) < 0.0)
435        below_plane[2] = 1;
436      else
437 <      below_plane[2]=0;
437 >      below_plane[2] = 0;
438          /* Test if b[i] lies in or on triangle a */
439 <    test = test_point_against_spherical_tri(t0,t1,t2,p2,
522 <                                                 n,nset,&which,pt_sides[2]);
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);
529 –        /* Remember if b[i] fell on one of the 3 defining planes */
530 –        if(test)
531 –          on_edge++;
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);
537 –    /* Accept if all points lie on a triangle vertex/edge edge- accept*/
538 –    if(on_edge == 3)
539 –       return(TRUE);
451      /* Now check vertices in a against triangle b */
452      return(FALSE);
453   }
# Line 544 | Line 455 | char pt_sides[3][3];
455  
456   set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,n)
457     FVECT t0,t1,t2,p0,p1,p2;
458 <   char test[3];
459 <   char sides[3][3];
460 <   char nset;
458 >   int test[3];
459 >   int sides[3][3];
460 >   int nset;
461     FVECT n[3];
462   {
463 <    char t;
463 >    int t;
464      double d;
465  
466      
# Line 561 | Line 472 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
472          VCROSS(n[0],t1,t0);
473        /* Test the point for sidedness */
474        d  = DOT(n[0],p0);
475 <      if(d >= 0)
475 >      if(d >= 0.0)
476          SET_NTH_BIT(test[0],0);
477      }
478      else
# Line 574 | Line 485 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
485          VCROSS(n[1],t2,t1);
486          /* Test the point for sidedness */
487          d  = DOT(n[1],p0);
488 <        if(d >= 0)
488 >        if(d >= 0.0)
489            SET_NTH_BIT(test[0],1);
490      }
491      else
# Line 587 | Line 498 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
498          VCROSS(n[2],t0,t2);
499        /* Test the point for sidedness */
500        d  = DOT(n[2],p0);
501 <      if(d >= 0)
501 >      if(d >= 0.0)
502          SET_NTH_BIT(test[0],2);
503      }
504      else
# Line 603 | Line 514 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
514          VCROSS(n[0],t1,t0);
515        /* Test the point for sidedness */
516        d  = DOT(n[0],p1);
517 <      if(d >= 0)
517 >      if(d >= 0.0)
518          SET_NTH_BIT(test[1],0);
519      }
520      else
# Line 617 | Line 528 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
528          VCROSS(n[1],t2,t1);
529        /* Test the point for sidedness */
530        d  = DOT(n[1],p1);
531 <      if(d >= 0)
531 >      if(d >= 0.0)
532          SET_NTH_BIT(test[1],1);
533      }
534      else
# Line 631 | Line 542 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
542          VCROSS(n[2],t0,t2);
543        /* Test the point for sidedness */
544        d  = DOT(n[2],p1);
545 <      if(d >= 0)
545 >      if(d >= 0.0)
546          SET_NTH_BIT(test[1],2);
547      }
548      else
# Line 647 | Line 558 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
558          VCROSS(n[0],t1,t0);
559        /* Test the point for sidedness */
560        d  = DOT(n[0],p2);
561 <      if(d >= 0)
561 >      if(d >= 0.0)
562          SET_NTH_BIT(test[2],0);
563      }
564      else
# Line 660 | Line 571 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
571          VCROSS(n[1],t2,t1);
572        /* Test the point for sidedness */
573        d  = DOT(n[1],p2);
574 <      if(d >= 0)
574 >      if(d >= 0.0)
575          SET_NTH_BIT(test[2],1);
576      }
577      else
# Line 673 | Line 584 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
584          VCROSS(n[2],t0,t2);
585        /* Test the point for sidedness */
586        d  = DOT(n[2],p2);
587 <      if(d >= 0)
587 >      if(d >= 0.0)
588          SET_NTH_BIT(test[2],2);
589      }
590      else
# Line 683 | Line 594 | set_sidedness_tests(t0,t1,t2,p0,p1,p2,test,sides,nset,
594  
595  
596   int
597 < spherical_tri_intersect(a1,a2,a3,b1,b2,b3)
597 > stri_intersect(a1,a2,a3,b1,b2,b3)
598   FVECT a1,a2,a3,b1,b2,b3;
599   {
600 <  char which,test,n_set[2];
601 <  char sides[2][3][3],i,j,inext,jnext;
602 <  char tests[2][3];
600 >  int which,test,n_set[2];
601 >  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
# Line 696 | Line 607 | FVECT a1,a2,a3,b1,b2,b3;
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(test_vertices_for_tri_inclusion(a1,a2,a3,b1,b2,b3,
700 <                                    &(n_set[0]),n[0],avg[0],sides[1]))
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(test_vertices_for_tri_inclusion(b1,b2,b3,a1,a2,a3,
704 <                                    &(n_set[1]),n[1],avg[1],sides[0]))
613 > if(vertices_in_stri(b1,b2,b3,a1,a2,a3,&(n_set[1]),n[1],avg[1],sides[0]))
614       return(TRUE);
615  
616  
# Line 749 | Line 658 | FVECT a1,a2,a3,b1,b2,b3;
658   }
659  
660   int
661 < ray_intersect_tri(orig,dir,v0,v1,v2,pt,wptr)
661 > ray_intersect_tri(orig,dir,v0,v1,v2,pt)
662   FVECT orig,dir;
663   FVECT v0,v1,v2;
664   FVECT pt;
756 – char *wptr;
665   {
666    FVECT p0,p1,p2,p,n;
759 –  char type,which;
667    double pd;
668 <  
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);
668 >  int type;
669  
670 <  if(type)
670 >  VSUB(p0,v0,orig);
671 >  VSUB(p1,v1,orig);
672 >  VSUB(p2,v2,orig);
673 >
674 >  if(point_in_stri(p0,p1,p2,dir))
675    {
676        /* Intersect the ray with the triangle plane */
677        tri_plane_equation(v0,v1,v2,n,&pd,FALSE);
678 <      intersect_ray_plane(orig,dir,n,pd,NULL,pt);        
678 >      return(intersect_ray_plane(orig,dir,n,pd,NULL,pt));
679    }
680 <  if(wptr)
774 <    *wptr = which;
775 <
776 <  return(type);
680 >  return(FALSE);
681   }
682  
683  
# Line 832 | Line 736 | FVECT fnear[4],ffar[4];
736      ffar[3][2] =  width*nhv[2] - height*nvv[2] + far*ndv[2] + vp[2] ;
737   }
738  
739 + int
740 + max_index(v,r)
741 + FVECT v;
742 + double *r;
743 + {
744 +  double p[3];
745 +  int i;
746  
747 +  p[0] = fabs(v[0]);
748 +  p[1] = fabs(v[1]);
749 +  p[2] = fabs(v[2]);
750 +  i = (p[0]>=p[1])?((p[0]>=p[2])?0:2):((p[1]>=p[2])?1:2);  
751 +  if(r)
752 +    *r = p[i];
753 +  return(i);
754 + }
755  
756 + int
757 + closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
758 + FVECT p0,p1,p2,p;
759 + int p0id,p1id,p2id;
760 + {
761 +    double d,d1;
762 +    int i;
763 +    
764 +    d =  DIST_SQ(p,p0);
765 +    d1 = DIST_SQ(p,p1);
766 +    if(d < d1)
767 +    {
768 +      d1 = DIST_SQ(p,p2);
769 +      i = (d1 < d)?p2id:p0id;
770 +    }
771 +    else
772 +    {
773 +      d = DIST_SQ(p,p2);
774 +      i = (d < d1)? p2id:p1id;
775 +    }
776 +    return(i);
777 + }
778  
779 +
780   int
781 < spherical_polygon_edge_intersect(a0,a1,b0,b1)
781 > sedge_intersect(a0,a1,b0,b1)
782   FVECT a0,a1,b0,b1;
783   {
784      FVECT na,nb,avga,avgb,p;
# Line 896 | Line 838 | double coord[3];
838    a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
839    coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
840    coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
841 <  coord[2]  = 1.0 - coord[0] - coord[1];
841 >  coord[2] = ((x1 - px) * (y2 - py) - (x2 - px) * (y1 - py)) / a;
842  
843   }
844  
845 + bary_ith_child(coord,i)
846 + double coord[3];
847 + int i;
848 + {
849 +
850 +  switch(i){
851 +  case 0:
852 +      /* update bary for child */
853 +      coord[0] = 2.0*coord[0]- 1.0;
854 +      coord[1] *= 2.0;
855 +      coord[2] *= 2.0;
856 +      break;
857 +  case 1:
858 +    coord[0] *= 2.0;
859 +    coord[1] = 2.0*coord[1]- 1.0;
860 +    coord[2] *= 2.0;
861 +    break;
862 +  case 2:
863 +    coord[0] *= 2.0;
864 +    coord[1] *= 2.0;
865 +    coord[2] = 2.0*coord[2]- 1.0;
866 +    break;
867 +  case 3:
868 +    coord[0] = 1.0 - 2.0*coord[0];
869 +    coord[1] = 1.0 - 2.0*coord[1];
870 +    coord[2] = 1.0 - 2.0*coord[2];
871 +    break;
872 + #ifdef DEBUG
873 +  default:
874 +    eputs("bary_ith_child():Invalid child\n");
875 +    break;
876 + #endif
877 +  }
878 + }
879 +
880 +
881   int
882 < bary2d_child(coord)
882 > bary_child(coord)
883   double coord[3];
884   {
885    int i;
886  
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);
913 –
914 –  /* Check if one of the new vertices: for all return center child */
915 –  if(ZERO(coord[0]) && EQUAL(coord[1],0.5))
916 –  {
917 –    coord[0] = 1.0f;
918 –    coord[1] = 0.0f;
919 –    coord[2] = 0.0f;
920 –    return(3);
921 –  }
922 –  if(ZERO(coord[1]) && EQUAL(coord[0],0.5))
923 –  {
924 –    coord[0] = 0.0f;
925 –    coord[1] = 1.0f;
926 –    coord[2] = 0.0f;
927 –    return(3);
928 –  }
929 –  if(ZERO(coord[2]) && EQUAL(coord[0],0.5))
930 –    {
931 –      coord[0] = 0.0f;
932 –      coord[1] = 0.0f;
933 –      coord[2] = 1.0f;
934 –      return(3);
935 –    }
936 –
937 –  /* Otherwise return child */
887    if(coord[0] > 0.5)
888    {
889        /* update bary for child */
# Line 968 | Line 917 | double coord[3];
917           }
918   }
919  
920 < int
921 < max_index(v)
922 < FVECT v;
920 > /* Coord was the ith child of the parent: set the coordinate
921 >   relative to the parent
922 > */
923 > bary_parent(coord,i)
924 > double coord[3];
925 > int i;
926   {
975 –  double a,b,c;
976 –  int i;
927  
928 <  a = fabs(v[0]);
929 <  b = fabs(v[1]);
930 <  c = fabs(v[2]);
931 <  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
932 <  return(i);
928 >  switch(i) {
929 >  case 0:
930 >    /* update bary for child */
931 >    coord[0] = coord[0]*0.5 + 0.5;
932 >    coord[1] *= 0.5;
933 >    coord[2] *= 0.5;
934 >    break;
935 >  case 1:
936 >    coord[0] *= 0.5;
937 >    coord[1]  = coord[1]*0.5 + 0.5;
938 >    coord[2] *= 0.5;
939 >    break;
940 >    
941 >  case 2:
942 >    coord[0] *= 0.5;
943 >    coord[1] *= 0.5;
944 >    coord[2] = coord[2]*0.5 + 0.5;
945 >    break;
946 >    
947 >  case 3:
948 >    coord[0] = 0.5 - 0.5*coord[0];
949 >    coord[1] = 0.5 - 0.5*coord[1];
950 >    coord[2] = 0.5 - 0.5*coord[2];
951 >    break;
952 > #ifdef DEBUG
953 >  default:
954 >    eputs("bary_parent():Invalid child\n");
955 >    break;
956 > #endif
957 >  }
958   }
959  
960 + bary_from_child(coord,child,next)
961 + double coord[3];
962 + int child,next;
963 + {
964 + #ifdef DEBUG
965 +  if(child <0 || child > 3)
966 +  {
967 +    eputs("bary_from_child():Invalid child\n");
968 +    return;
969 +  }
970 +  if(next <0 || next > 3)
971 +  {
972 +    eputs("bary_from_child():Invalid next\n");
973 +    return;
974 +  }
975 + #endif
976 +  if(next == child)
977 +    return;
978  
979 +  switch(child){
980 +  case 0:
981 +    switch(next){
982 +    case 1:
983 +      coord[0] += 1.0;
984 +      coord[1] -= 1.0;
985 +      break;
986 +    case 2:
987 +      coord[0] += 1.0;
988 +      coord[2] -= 1.0;
989 +      break;
990 +    case 3:
991 +      coord[0] *= -1.0;
992 +      coord[1] = 1 - coord[1];
993 +      coord[2] = 1 - coord[2];
994 +      break;
995  
996 < /*
997 < * int
998 < * smRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
999 < *
1000 < *   Intersect the ray with triangle v0v1v2, return intersection point in r
1001 < *
1002 < *    Assumes orig,v0,v1,v2 are in spherical coordinates, and orig is
1003 < *    unit
1004 < */
996 >    }
997 >    break;
998 >  case 1:
999 >    switch(next){
1000 >    case 0:
1001 >      coord[0] -= 1.0;
1002 >      coord[1] += 1.0;
1003 >      break;
1004 >    case 2:
1005 >      coord[1] += 1.0;
1006 >      coord[2] -= 1.0;
1007 >      break;
1008 >    case 3:
1009 >      coord[0] = 1 - coord[0];
1010 >      coord[1] *= -1.0;
1011 >      coord[2] = 1 - coord[2];
1012 >      break;
1013 >    }
1014 >    break;
1015 >  case 2:
1016 >    switch(next){
1017 >    case 0:
1018 >      coord[0] -= 1.0;
1019 >      coord[2] += 1.0;
1020 >      break;
1021 >    case 1:
1022 >      coord[1] -= 1.0;
1023 >      coord[2] += 1.0;
1024 >      break;
1025 >    case 3:
1026 >      coord[0] = 1 - coord[0];
1027 >      coord[1] = 1 - coord[1];
1028 >      coord[2] *= -1.0;
1029 >      break;
1030 >    }
1031 >    break;
1032 >  case 3:
1033 >    switch(next){
1034 >    case 0:
1035 >      coord[0] *= -1.0;
1036 >      coord[1] = 1 - coord[1];
1037 >      coord[2] = 1 - coord[2];
1038 >      break;
1039 >    case 1:
1040 >      coord[0] = 1 - coord[0];
1041 >      coord[1] *= -1.0;
1042 >      coord[2] = 1 - coord[2];
1043 >      break;
1044 >    case 2:
1045 >      coord[0] = 1 - coord[0];
1046 >      coord[1] = 1 - coord[1];
1047 >      coord[2] *= -1.0;
1048 >      break;
1049 >    }
1050 >    break;
1051 >  }
1052 > }
1053 >
1054 >
1055 > baryi_parent(coord,i)
1056 > BCOORD coord[3];
1057 > int i;
1058 > {
1059 >
1060 >  switch(i) {
1061 >  case 0:
1062 >    /* update bary for child */
1063 >    coord[0] = (coord[0] >> 1) + MAXBCOORD2;
1064 >    coord[1] >>= 1;
1065 >    coord[2] >>= 1;
1066 >    break;
1067 >  case 1:
1068 >    coord[0] >>= 1;
1069 >    coord[1]  = (coord[1] >> 1) + MAXBCOORD2;
1070 >    coord[2] >>= 1;
1071 >    break;
1072 >    
1073 >  case 2:
1074 >    coord[0] >>= 1;
1075 >    coord[1] >>= 1;
1076 >    coord[2] = (coord[2] >> 1) + MAXBCOORD2;
1077 >    break;
1078 >    
1079 >  case 3:
1080 >    coord[0] = MAXBCOORD2 - (coord[0] >> 1);
1081 >    coord[1] = MAXBCOORD2 - (coord[1] >> 1);
1082 >    coord[2] = MAXBCOORD2 - (coord[2] >> 1);
1083 >    break;
1084 > #ifdef DEBUG
1085 >  default:
1086 >    eputs("baryi_parent():Invalid child\n");
1087 >    break;
1088 > #endif
1089 >  }
1090 > }
1091 >
1092 > baryi_from_child(coord,child,next)
1093 > BCOORD coord[3];
1094 > int child,next;
1095 > {
1096 > #ifdef DEBUG
1097 >  if(child <0 || child > 3)
1098 >  {
1099 >    eputs("baryi_from_child():Invalid child\n");
1100 >    return;
1101 >  }
1102 >  if(next <0 || next > 3)
1103 >  {
1104 >    eputs("baryi_from_child():Invalid next\n");
1105 >    return;
1106 >  }
1107 > #endif
1108 >  if(next == child)
1109 >    return;
1110 >
1111 >  switch(child){
1112 >  case 0:
1113 >      coord[0] = 0;
1114 >      coord[1] = MAXBCOORD - coord[1];
1115 >      coord[2] = MAXBCOORD - coord[2];
1116 >      break;
1117 >  case 1:
1118 >      coord[0] = MAXBCOORD - coord[0];
1119 >      coord[1] = 0;
1120 >      coord[2] = MAXBCOORD - coord[2];
1121 >      break;
1122 >  case 2:
1123 >      coord[0] = MAXBCOORD - coord[0];
1124 >      coord[1] = MAXBCOORD - coord[1];
1125 >      coord[2] = 0;
1126 >    break;
1127 >  case 3:
1128 >    switch(next){
1129 >    case 0:
1130 >      coord[0] = 0;
1131 >      coord[1] = MAXBCOORD - coord[1];
1132 >      coord[2] = MAXBCOORD - coord[2];
1133 >      break;
1134 >    case 1:
1135 >      coord[0] = MAXBCOORD - coord[0];
1136 >      coord[1] = 0;
1137 >      coord[2] = MAXBCOORD - coord[2];
1138 >      break;
1139 >    case 2:
1140 >      coord[0] = MAXBCOORD - coord[0];
1141 >      coord[1] = MAXBCOORD - coord[1];
1142 >      coord[2] = 0;
1143 >      break;
1144 >    }
1145 >    break;
1146 >  }
1147 > }
1148 >
1149   int
1150 < traceRay(orig,dir,v0,v1,v2,r)
1151 <  FVECT orig,dir;
999 <  FVECT v0,v1,v2;
1000 <  FVECT r;
1150 > baryi_child(coord)
1151 > BCOORD coord[3];
1152   {
1153 <  FVECT n,p[3],d;
1003 <  double pt[3],r_eps;
1004 <  char i;
1005 <  int which;
1153 >  int i;
1154  
1155 <  /* Find the plane equation for the triangle defined by the edge v0v1 and
1156 <   the view center*/
1157 <  VCROSS(n,v0,v1);
1158 <  /* Intersect the ray with this plane */
1159 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[0]),p[0]);
1160 <  if(i)
1161 <    which = 0;
1155 >  if(coord[0] > MAXBCOORD2)
1156 >  {
1157 >      /* update bary for child */
1158 >      coord[0] = (coord[0]<< 1) - MAXBCOORD;
1159 >      coord[1] <<= 1;
1160 >      coord[2] <<= 1;
1161 >      return(0);
1162 >  }
1163    else
1164 <    which = -1;
1164 >    if(coord[1] > MAXBCOORD2)
1165 >    {
1166 >      coord[0] <<= 1;
1167 >      coord[1] = (coord[1] << 1) - MAXBCOORD;
1168 >      coord[2] <<= 1;
1169 >      return(1);
1170 >    }
1171 >    else
1172 >      if(coord[2] > MAXBCOORD2)
1173 >      {
1174 >        coord[0] <<= 1;
1175 >        coord[1] <<= 1;
1176 >        coord[2] = (coord[2] << 1) - MAXBCOORD;
1177 >        return(2);
1178 >      }
1179 >      else
1180 >         {
1181 >           coord[0] = MAXBCOORD - (coord[0] << 1);
1182 >           coord[1] = MAXBCOORD - (coord[1] << 1);
1183 >           coord[2] = MAXBCOORD - (coord[2] << 1);
1184 >           return(3);
1185 >         }
1186 > }
1187  
1188 <  VCROSS(n,v1,v2);
1189 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[1]),p[1]);
1190 <  if(i)
1191 <    if((which==-1) || pt[1] < pt[0])
1192 <      which = 1;
1188 > /* convert barycentric coordinate b in [-eps,1+eps] to [0,MAXLONG],
1189 >   dir unbounded to [-MAXLONG,MAXLONG]
1190 > */
1191 > bary_dtol(b,db,bi,dbi,t)
1192 > double b[3],db[3][3];
1193 > BCOORD bi[3];
1194 > BDIR dbi[3][3];
1195 > TINT t[3];
1196 > {
1197 >  int i,id,j;
1198 >  double d;
1199  
1200 <  VCROSS(n,v2,v0);
1201 <  i = intersect_ray_plane(orig,dir,n,0.0,&(pt[2]),p[2]);
1202 <  if(i)
1203 <    if((which==-1) || pt[2] < pt[which])
1204 <      which = 2;
1200 >  for(i=0; i < 2;i++)
1201 >  {
1202 >    if(b[i] <= 0.0)
1203 >    {
1204 >      bi[i]= 0;
1205 >    }
1206 >    else
1207 >      if(b[i] >= 1.0)
1208 >      {
1209 >        bi[i]= MAXBCOORD;
1210 >      }
1211 >      else
1212 >        bi[i] = (BCOORD)(b[i]*MAXBCOORD);
1213 >  }
1214 >  bi[2] = MAXBCOORD - bi[0] - bi[1];
1215  
1216 <  if(which != -1)
1216 >  if(bi[2] < 0)
1217    {
1218 <      /* Return point that is K*eps along projection of the ray on
1219 <         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);
1218 >      bi[2] = 0;
1219 >      bi[1] = MAXBCOORD - bi[0];
1220    }
1221 <  else
1221 >  for(j=0; j< 3; j++)
1222    {
1223 <      /* Unable to find intersection: move along ray and try again */
1224 <      r_eps = -DOT(orig,dir);
1225 <      VSUM(n,dir,orig,r_eps);
1226 <      r_eps = DOT(n,dir);
1227 <      VSUM(r,orig,dir,(20*FTINY)/r_eps);
1228 <      normalize(r);
1229 < #ifdef DEBUG
1230 <      eputs("traceRay:Ray does not intersect triangle");
1231 < #endif
1058 <      return(FALSE);
1223 >    if(t[j]==0)
1224 >      continue;
1225 >    if(t[j] == HUGET)
1226 >       max_index(db[j],&d);
1227 >    for(i=0; i< 3; i++)
1228 >       if(t[j] != HUGET)
1229 >          dbi[j][i] = (BDIR)(db[j][i]*MAXBDIR);
1230 >       else
1231 >          dbi[j][i] = (BDIR)(db[j][i]/d*MAXBDIR);
1232    }
1233   }
1234  

Diff Legend

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