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.1 by gwlarson, Wed Aug 19 17:45:24 1998 UTC vs.
Revision 3.5 by gwlarson, Mon Sep 14 10:33:46 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,pd,r)
128 > intersect_vector_plane(v,plane_n,plane_d,tptr,r)
129     FVECT v,plane_n;
130     double plane_d;
131 <   double *pd;
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,pd,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;
152 <    if(pd)
153 <       *pd = 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);
167   }
168  
# Line 179 | Line 180 | intersect_ray_plane(orig,dir,plane_n,plane_d,pd,r)
180        Plane is Ax + By + Cz +D = 0:
181        plane[0] = A,plane[1] = B,plane[2] = C,plane[3] = D
182      */
183 <     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0 */
184 <    /* t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
185 <    /* line is  l = p1 + (p2-p1)t */
183 >     /*  A(orig[0] + dxt) + B(orig[1] + dyt) + C(orig[2] + dzt) + pd = 0
184 >         t = -(DOT(plane_n,orig)+ plane_d)/(DOT(plane_n,d))
185 >       line is  l = p1 + (p2-p1)t
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 219 | 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 251 | 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;
259 < char sides[3];
296 > int *nset;
297 > int sides[3];
298  
299   {
300 <    float d;
263 <
300 >    double d;
301      /* Find the normal to the triangle ORIGIN:v0:v1 */
302      if(!NTH_BIT(*nset,0))
303      {
# Line 270 | 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 <      {
278 <          sides[0] =  GT_OUT;
279 <          sides[1] = sides[2] = GT_INVALID;
280 <          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 290 | 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      {
295 –        sides[1] = GT_EDGE;
296 –        /* If on plane 0-and on plane 1: lies on edge */
297 –        if(sides[0] == GT_EDGE)
298 –        {
299 –            *which = 1;
300 –            sides[2] = GT_INVALID;
301 –            return(GT_EDGE);
302 –        }
303 –    }
304 –    else if(d > 0)
305 –    {
329          sides[1] = GT_OUT;
330          sides[2] = GT_INVALID;
331          return(FALSE);
# Line 312 | Line 335 | char sides[3];
335      /* Test next edge */
336      if(!NTH_BIT(*nset,2))
337      {
315 –
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 <
325 <        /* If on plane 0 and 2: lies on edge 0*/
326 <        if(sides[0] == GT_EDGE)
327 <           {
328 <               *which = 0;
329 <               return(GT_EDGE);
330 <           }
331 <        /* If on plane 1 and 2: lies on edge  2*/
332 <        if(sides[1] == GT_EDGE)
333 <           {
334 <               *which = 2;
335 <               return(GT_EDGE);
336 <           }
337 <        /* otherwise: on face 2 */
338 <        else
339 <           {
340 <               *which = 2;
341 <               return(GT_FACE);
342 <           }
345 >      sides[2] = GT_OUT;
346 >      return(FALSE);
347      }
344 –    else if(d > 0)
345 –      {
346 –        sides[2] = GT_OUT;
347 –        return(FALSE);
348 –      }
349 –    /* If on edge */
348      else
349         sides[2] = GT_INTERIOR;
352 –    
353 –    /* If on plane 0 only: on face 0 */
354 –    if(sides[0] == GT_EDGE)
355 –    {
356 –        *which = 0;
357 –        return(GT_FACE);
358 –    }
359 –    /* If on plane 1 only: on face 1 */
360 –    if(sides[1] == GT_EDGE)
361 –    {
362 –        *which = 1;
363 –        return(GT_FACE);
364 –    }
350      /* Must be interior to the pyramid */
351      return(GT_INTERIOR);
352   }
# Line 370 | 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;
375 – char *which;
360   {
361 <    float d;
361 >    double d;
362      FVECT n;  
379 –    char sides[3];
363  
381 –    /* First test if point coincides with any of the vertices */
382 –    if(EQUAL_VEC3(p,v0))
383 –    {
384 –        *which = 0;
385 –        return(GT_VERTEX);
386 –    }
387 –    if(EQUAL_VEC3(p,v1))
388 –    {
389 –        *which = 1;
390 –        return(GT_VERTEX);
391 –    }
392 –    if(EQUAL_VEC3(p,v2))
393 –    {
394 –        *which = 2;
395 –        return(GT_VERTEX);
396 –    }
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
403 <       if(d > 0)
404 <          return(FALSE);
405 <       else
406 <          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))
412 <    {
413 <        sides[1] = GT_EDGE;
414 <        /* If on plane 0-and on plane 1: lies on edge */
415 <        if(sides[0] == GT_EDGE)
416 <        {
417 <            *which = 1;
418 <            return(GT_VERTEX);
419 <        }
420 <    }
421 <    else if(d > 0)
374 >    if(d > 0.0)
375         return(FALSE);
423 –    else
424 –       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))
431 <    {
432 <        sides[2] = GT_EDGE;
433 <        
434 <        /* If on plane 0 and 2: lies on edge 0*/
435 <        if(sides[0] == GT_EDGE)
436 <        {
437 <            *which = 0;
438 <            return(GT_VERTEX);
439 <        }
440 <        /* If on plane 1 and 2: lies on edge  2*/
441 <        if(sides[1] == GT_EDGE)
442 <        {
443 <            *which = 2;
444 <            return(GT_VERTEX);
445 <        }
446 <        /* otherwise: on face 2 */
447 <        else
448 <       {
449 <           return(GT_FACE);
450 <       }
451 <    }
452 <    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;
468 <    char which;
396 >    int below_plane[3],test;
397  
398      SUM_3VEC3(avg,t0,t1,t2);
471 –    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 476 | 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,
485 <                                                 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);
491 –        /* Remember if b[i] fell on one of the 3 defining planes */
492 –        if(test)
493 –          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,
503 <                                                 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);
509 –      /* Remember if b[i] fell on one of the 3 defining planes */
510 –      if(test)
511 –        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,
521 <                                                 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);
528 –        /* Remember if b[i] fell on one of the 3 defining planes */
529 –        if(test)
530 –          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);
536 –    /* Accept if all points lie on a triangle vertex/edge edge- accept*/
537 –    if(on_edge == 3)
538 –       return(TRUE);
451      /* Now check vertices in a against triangle b */
452      return(FALSE);
453   }
# Line 543 | 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 560 | 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 573 | 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 586 | 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 602 | 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 616 | 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 630 | 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 646 | 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 659 | 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 672 | 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 682 | 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 695 | 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,
699 <                                    &(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,
703 <                                    &(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 748 | 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;
755 – char *wptr;
665   {
666    FVECT p0,p1,p2,p,n;
758 –  char type,which;
667    double pd;
668 <  
761 <  point_on_sphere(p0,v0,orig);
762 <  point_on_sphere(p1,v1,orig);
763 <  point_on_sphere(p2,v2,orig);
764 <  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)
773 <    *wptr = which;
774 <
775 <  return(type);
680 >  return(FALSE);
681   }
682  
683  
# Line 835 | 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 879 | Line 784 | FVECT a0,a1,b0,b1;
784        return(FALSE);
785      return(TRUE);
786   }
787 +
788 +
789 +
790 + /* Find the normalized barycentric coordinates of p relative to
791 + * triangle v0,v1,v2. Return result in coord
792 + */
793 + bary2d(x1,y1,x2,y2,x3,y3,px,py,coord)
794 + double x1,y1,x2,y2,x3,y3;
795 + double px,py;
796 + double coord[3];
797 + {
798 +  double a;
799 +
800 +  a =  (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1);
801 +  coord[0] = ((x2 - px) * (y3 - py) - (x3 - px) * (y2 - py)) / a;
802 +  coord[1] = ((x3 - px) * (y1 - py) - (x1 - px) * (y3 - py)) / a;
803 +  coord[2]  = 1.0 - coord[0] - coord[1];
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 + bary_child(coord)
845 + double coord[3];
846 + {
847 +  int i;
848 +
849 +  if(coord[0] > 0.5)
850 +  {
851 +      /* update bary for child */
852 +      coord[0] = 2.0*coord[0]- 1.0;
853 +      coord[1] *= 2.0;
854 +      coord[2] *= 2.0;
855 +      return(0);
856 +  }
857 +  else
858 +    if(coord[1] > 0.5)
859 +    {
860 +      coord[0] *= 2.0;
861 +      coord[1] = 2.0*coord[1]- 1.0;
862 +      coord[2] *= 2.0;
863 +      return(1);
864 +    }
865 +    else
866 +      if(coord[2] > 0.5)
867 +      {
868 +        coord[0] *= 2.0;
869 +        coord[1] *= 2.0;
870 +        coord[2] = 2.0*coord[2]- 1.0;
871 +        return(2);
872 +      }
873 +      else
874 +         {
875 +           coord[0] = 1.0 - 2.0*coord[0];
876 +           coord[1] = 1.0 - 2.0*coord[1];
877 +           coord[2] = 1.0 - 2.0*coord[2];
878 +           return(3);
879 +         }
880 + }
881 +
882 + /* Coord was the ith child of the parent: set the coordinate
883 +   relative to the parent
884 + */
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;
1019 + {
1020 +  double a,b,c;
1021 +  int i;
1022 +
1023 +  a = fabs(v[0]);
1024 +  b = fabs(v[1]);
1025 +  c = fabs(v[2]);
1026 +  i = (a>=b)?((a>=c)?0:2):((b>=c)?1:2);  
1027 +  return(i);
1028 + }
1029 +
1030 + int
1031 + closest_point_in_tri(p0,p1,p2,p,p0id,p1id,p2id)
1032 + FVECT p0,p1,p2,p;
1033 + int p0id,p1id,p2id;
1034 + {
1035 +    double d,d1;
1036 +    int i;
1037 +    
1038 +    d =  DIST_SQ(p,p0);
1039 +    d1 = DIST_SQ(p,p1);
1040 +    if(d < d1)
1041 +    {
1042 +      d1 = DIST_SQ(p,p2);
1043 +      i = (d1 < d)?p2id:p0id;
1044 +    }
1045 +    else
1046 +    {
1047 +      d = DIST_SQ(p,p2);
1048 +      i = (d < d1)? p2id:p1id;
1049 +    }
1050 +    return(i);
1051 + }
1052 +
1053 +
1054 +
1055 +
1056 +
1057 +
1058  

Diff Legend

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