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.4 by gwlarson, Fri Sep 11 11:52:25 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 <  
668 >  int type;
669 >
670    point_on_sphere(p0,v0,orig);
671    point_on_sphere(p1,v1,orig);
672    point_on_sphere(p2,v2,orig);
673 <  type = test_single_point_against_spherical_tri(p0,p1,p2,dir,&which);
674 <
767 <  if(type)
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 836 | Line 740 | FVECT fnear[4],ffar[4];
740  
741  
742   int
743 < spherical_polygon_edge_intersect(a0,a1,b0,b1)
743 > sedge_intersect(a0,a1,b0,b1)
744   FVECT a0,a1,b0,b1;
745   {
746      FVECT na,nb,avga,avgb,p;
# Line 900 | Line 804 | double coord[3];
804  
805   }
806  
807 + bary_ith_child(coord,i)
808 + double coord[3];
809 + int i;
810 + {
811 +
812 +  switch(i){
813 +  case 0:
814 +      /* update bary for child */
815 +      coord[0] = 2.0*coord[0]- 1.0;
816 +      coord[1] *= 2.0;
817 +      coord[2] *= 2.0;
818 +      break;
819 +  case 1:
820 +    coord[0] *= 2.0;
821 +    coord[1] = 2.0*coord[1]- 1.0;
822 +    coord[2] *= 2.0;
823 +    break;
824 +  case 2:
825 +    coord[0] *= 2.0;
826 +    coord[1] *= 2.0;
827 +    coord[2] = 2.0*coord[2]- 1.0;
828 +    break;
829 +  case 3:
830 +    coord[0] = 1.0 - 2.0*coord[0];
831 +    coord[1] = 1.0 - 2.0*coord[1];
832 +    coord[2] = 1.0 - 2.0*coord[2];
833 +    break;
834 + #ifdef DEBUG
835 +  default:
836 +    eputs("bary_ith_child():Invalid child\n");
837 +    break;
838 + #endif
839 +  }
840 + }
841 +
842 +
843   int
844 < bary2d_child(coord)
844 > bary_child(coord)
845   double coord[3];
846   {
847    int i;
848  
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 */
849    if(coord[0] > 0.5)
850    {
851        /* update bary for child */
# Line 968 | Line 879 | double coord[3];
879           }
880   }
881  
882 + /* Coord was the ith child of the parent: set the coordinate
883 +   relative to the parent
884 + */
885 + bary_parent(coord,i)
886 + double coord[3];
887 + int i;
888 + {
889 +
890 +  switch(i) {
891 +  case 0:
892 +    /* update bary for child */
893 +    coord[0] = coord[0]*0.5 + 0.5;
894 +    coord[1] *= 0.5;
895 +    coord[2] *= 0.5;
896 +    break;
897 +  case 1:
898 +    coord[0] *= 0.5;
899 +    coord[1]  = coord[1]*0.5 + 0.5;
900 +    coord[2] *= 0.5;
901 +    break;
902 +    
903 +  case 2:
904 +    coord[0] *= 0.5;
905 +    coord[1] *= 0.5;
906 +    coord[2] = coord[2]*0.5 + 0.5;
907 +    break;
908 +    
909 +  case 3:
910 +    coord[0] = 0.5 - 0.5*coord[0];
911 +    coord[1] = 0.5 - 0.5*coord[1];
912 +    coord[2] = 0.5 - 0.5*coord[2];
913 +    break;
914 + #ifdef DEBUG
915 +  default:
916 +    eputs("bary_parent():Invalid child\n");
917 +    break;
918 + #endif
919 +  }
920 + }
921 +
922 + bary_from_child(coord,child,next)
923 + double coord[3];
924 + int child,next;
925 + {
926 + #ifdef DEBUG
927 +  if(child <0 || child > 3)
928 +  {
929 +    eputs("bary_from_child():Invalid child\n");
930 +    return;
931 +  }
932 +  if(next <0 || next > 3)
933 +  {
934 +    eputs("bary_from_child():Invalid next\n");
935 +    return;
936 +  }
937 + #endif
938 +  if(next == child)
939 +    return;
940 +
941 +  switch(child){
942 +  case 0:
943 +    switch(next){
944 +    case 1:
945 +      coord[0] += 1.0;
946 +      coord[1] -= 1.0;
947 +      break;
948 +    case 2:
949 +      coord[0] += 1.0;
950 +      coord[2] -= 1.0;
951 +      break;
952 +    case 3:
953 +      coord[0] *= -1.0;
954 +      coord[1] = 1 - coord[1];
955 +      coord[2] = 1 - coord[2];
956 +      break;
957 +
958 +    }
959 +    break;
960 +  case 1:
961 +    switch(next){
962 +    case 0:
963 +      coord[0] -= 1.0;
964 +      coord[1] += 1.0;
965 +      break;
966 +    case 2:
967 +      coord[1] += 1.0;
968 +      coord[2] -= 1.0;
969 +      break;
970 +    case 3:
971 +      coord[0] = 1 - coord[0];
972 +      coord[1] *= -1.0;
973 +      coord[2] = 1 - coord[2];
974 +      break;
975 +    }
976 +    break;
977 +  case 2:
978 +    switch(next){
979 +    case 0:
980 +      coord[0] -= 1.0;
981 +      coord[2] += 1.0;
982 +      break;
983 +    case 1:
984 +      coord[1] -= 1.0;
985 +      coord[2] += 1.0;
986 +      break;
987 +    case 3:
988 +      coord[0] = 1 - coord[0];
989 +      coord[1] = 1 - coord[1];
990 +      coord[2] *= -1.0;
991 +      break;
992 +    }
993 +    break;
994 +  case 3:
995 +    switch(next){
996 +    case 0:
997 +      coord[0] *= -1.0;
998 +      coord[1] = 1 - coord[1];
999 +      coord[2] = 1 - coord[2];
1000 +      break;
1001 +    case 1:
1002 +      coord[0] = 1 - coord[0];
1003 +      coord[1] *= -1.0;
1004 +      coord[2] = 1 - coord[2];
1005 +      break;
1006 +    case 2:
1007 +      coord[0] = 1 - coord[0];
1008 +      coord[1] = 1 - coord[1];
1009 +      coord[2] *= -1.0;
1010 +      break;
1011 +    }
1012 +    break;
1013 +  }
1014 + }
1015 +
1016   int
1017   max_index(v)
1018   FVECT v;
# Line 986 | Line 1031 | FVECT v;
1031  
1032   /*
1033   * int
1034 < * smRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
1034 > * traceRay(FVECT orig, FVECT dir,FVECT v0,FVECT v1,FVECT v2,FVECT r)
1035   *
1036   *   Intersect the ray with triangle v0v1v2, return intersection point in r
1037   *
# Line 1001 | Line 1046 | traceRay(orig,dir,v0,v1,v2,r)
1046   {
1047    FVECT n,p[3],d;
1048    double pt[3],r_eps;
1049 <  char i;
1049 >  int i;
1050    int which;
1051  
1052    /* Find the plane equation for the triangle defined by the edge v0v1 and
# Line 1053 | Line 1098 | traceRay(orig,dir,v0,v1,v2,r)
1098        VSUM(r,orig,dir,(20*FTINY)/r_eps);
1099        normalize(r);
1100   #ifdef DEBUG
1101 <      eputs("traceRay:Ray does not intersect triangle");
1101 >      eputs("traceRay:Ray does not intersect triangle\n");
1102   #endif
1103        return(FALSE);
1104    }
1105   }
1106  
1107  
1108 + int
1109 + closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
1110 + FVECT p0,p1,p2,p;
1111 + int p0id,p1id,p2id;
1112 + {
1113 +    double d,d1;
1114 +    int i;
1115 +    
1116 +    d =  DIST_SQ(p,p0);
1117 +    d1 = DIST_SQ(p,p1);
1118 +    if(d < d1)
1119 +    {
1120 +      d1 = DIST_SQ(p,p2);
1121 +      i = (d1 < d)?p2id:p0id;
1122 +    }
1123 +    else
1124 +    {
1125 +      d = DIST_SQ(p,p2);
1126 +      i = (d < d1)? p2id:p1id;
1127 +    }
1128 +    return(i);
1129 + }
1130  
1131  
1132  

Diff Legend

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