home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch2_2 / ptinply3.c next >
C/C++ Source or Header  |  1995-04-04  |  4KB  |  167 lines
  1. #include <math.h>
  2. #include <stdlib.h>
  3. #include <stdio.h>
  4.  
  5. #ifndef max
  6. #define max(a,b) ((a)>(b)?(a):(b))
  7. #define min(a,b) ((a)<(b)?(a):(b))
  8. #endif
  9. #define PI 3.141592653589793324
  10. #define GeoZeroVec(v) ((v).x = (v).y = (v).z = 0.0)
  11. #define GeoMultVec(a,b,c) do {(c).x = a*(b).x; (c).y = a*(b).y;    (c).z = a*(b).z; } while (0)
  12. #define Geo_Vet(a,b,c) do {(c).x = (b).x-(a).x; (c).y = (b).y-(a).y; (c).z = (b).z-(a).z;} while (0)
  13.  
  14. typedef double Rdouble;
  15. typedef float  Rfloat;
  16. typedef struct _GeoPoint { Rfloat x, y, z; } GeoPoint;
  17.  
  18. /*=========================  Geometrical Procedures  ======================= */
  19.  
  20. Rdouble GeoDotProd ( GeoPoint *vec0, GeoPoint *vec1 )
  21. {
  22.  return ( vec0->x * vec1->x + vec0->y * vec1->y + vec0->z * vec1->z );
  23. }
  24.  
  25. void GeoCrossProd ( GeoPoint *in0, GeoPoint *in1, GeoPoint *out )
  26. {
  27.  out->x = (in0->y * in1->z) - (in0->z * in1->y);
  28.  out->y = (in0->z * in1->x) - (in0->x * in1->z);
  29.  out->z = (in0->x * in1->y) - (in0->y * in1->x);
  30. }
  31.  
  32. Rdouble GeoTripleProd ( GeoPoint *vec0, GeoPoint *vec1, GeoPoint *vec2 )
  33. {
  34.  GeoPoint tmp;
  35.  
  36.  GeoCrossProd ( vec0, vec1, &tmp );
  37.  return ( GeoDotProd( &tmp, vec2 ) );
  38. }
  39.  
  40. Rdouble GeoVecLen ( GeoPoint *vec )
  41. {
  42.  return sqrt ( GeoDotProd ( vec, vec ) );
  43. }
  44.  
  45. int GeoPolyNormal ( int    n_verts, GeoPoint *verts, GeoPoint *n ) 
  46. {
  47.  int      i;
  48.  Rfloat      n_size;        
  49.  GeoPoint v0, v1, p; 
  50.  
  51.  GeoZeroVec ( *n );
  52.  Geo_Vet ( verts[0], verts[1], v0 );
  53.  for ( i = 2; i < n_verts; i++ )
  54.      {
  55.       Geo_Vet ( verts[0], verts[i], v1 );
  56.       GeoCrossProd ( &v0, &v1, &p );
  57.       n->x += p.x; n->y += p.y; n->z += p.z;
  58.       v0 = v1;
  59.      }
  60.  
  61.  n_size = GeoVecLen ( n );
  62.  if ( n_size > 0.0 )
  63.     {
  64.      GeoMultVec ( 1/n_size, *n, *n );
  65.      return 1;
  66.     }
  67.  else
  68.      return 0;
  69. }
  70.  
  71. /*=========================  geo_solid_angle  =========================*/
  72. /* 
  73.   Calculates the solid angle given by the spherical projection of 
  74.   a 3D plane polygon
  75. */
  76.  
  77. Rdouble geo_solid_angle ( 
  78.         int      n_vert,  /* number of vertices */
  79.         GeoPoint *verts,  /* vertex coordinates list */
  80.         GeoPoint *p )     /* point to be tested */
  81. {
  82.  int      i;
  83.  Rdouble  area = 0.0, ang, s, l1, l2;
  84.  GeoPoint p1, p2, r1, a, b, n1, n2;
  85.  GeoPoint plane;
  86.  
  87.  if ( n_vert < 3 ) return 0.0;
  88.  
  89.  GeoPolyNormal ( n_vert, verts, &plane );
  90.  
  91.  /* 
  92.     WARNING: at this point, a practical implementation should check
  93.     whether p is too close to the polygon plane. If it is, then
  94.     there are two possibilities: 
  95.       a) if the projection of p onto the plane is outside the 
  96.          polygon, then area zero should be returned;
  97.       b) otherwise, p is on the polyhedron boundary.
  98.  */ 
  99.  
  100.  p2 = verts[n_vert-1];  /* last vertex */
  101.  p1 = verts[0];         /* first vertex */
  102.  Geo_Vet ( p1, p2, a ); /* a = p2 - p1 */  
  103.  
  104.  for ( i = 0; i < n_vert; i++ )
  105.      {
  106.       Geo_Vet(*p, p1, r1); 
  107.       p2 = verts[(i+1)%n_vert];
  108.       Geo_Vet ( p1, p2, b );
  109.       GeoCrossProd ( &a, &r1, &n1 );
  110.       GeoCrossProd ( &r1, &b, &n2 );
  111.     
  112.       l1 = GeoVecLen ( &n1 );
  113.       l2 = GeoVecLen ( &n2 );
  114.       s  = GeoDotProd ( &n1, &n2 ) / ( l1 * l2 );
  115.       ang = acos ( max(-1.0,min(1.0,s)) );
  116.       s = GeoTripleProd( &b, &a, &plane );
  117.       area += s > 0.0 ? PI - ang : PI + ang;
  118.      
  119.       GeoMultVec ( -1.0, b, a );
  120.       p1 = p2;
  121.      }
  122.  
  123.  area -= PI*(n_vert-2);
  124.  
  125.  return ( GeoDotProd ( &plane, &r1 ) > 0.0 ) ? -area : area; 
  126. }
  127.  
  128. /* ======================  main ========================== */
  129.  
  130. int main ( void ) 
  131. {
  132.  FILE     *f;
  133.  char     s[32];
  134.  int      nv, j;
  135.  GeoPoint verts[100], p;
  136.  Rdouble  Area =0.0;
  137.      
  138.  fprintf ( stdout, "\nFile Name: " );
  139.  gets ( s );
  140.  if ( (f = fopen ( s, "r" )) == NULL )
  141.     {
  142.      fprintf ( stdout, "Can not open the Polyhedron file \n" );
  143.      exit ( 1 );
  144.     }
  145.  fprintf ( stdout, "\nPoint to be tested: " );
  146.  fscanf( stdin, "%f %f %f", &p.x, &p.y, &p.z );
  147.  
  148.  while ( fscanf ( f, "%d", &nv ) == 1 )
  149.        {    
  150.         for ( j = 0; j < nv; j++ )
  151.             if ( fscanf ( f, "%f %f %f", 
  152.                  &verts[j].x, &verts[j].y, &verts[j].z ) != 3 )
  153.                {
  154.                 fprintf ( stdout, "Invalid Polyhedron file \n" );
  155.                 exit ( 2 );
  156.                }
  157.  
  158.         Area += geo_solid_angle ( nv, verts, &p ); 
  159.        }
  160.         
  161.  fprintf ( stdout, "\n  Area = %12.4lf spherical radians.\n", Area);
  162.  fprintf ( stdout, "\n  The point is %s", 
  163.          ( (Area > 2*PI) || (Area < -2*PI) )? "inside" : "outside" );
  164.  fprintf ( stdout, "the given polyhedron \n" );
  165.  return 1;
  166. }
  167.