inhedron.c

判断点是否在多面体之中

] );
   }
   printf("n = %3d vertices read\n",n);
   putchar('\n');

   return n;
}

/*---------------------------------------------------------------------
a - b ==> c.
---------------------------------------------------------------------*/
void    SubVec( tPointi a, tPointi b, tPointi c )
{
   int i;

   for( i = 0; i < DIM; i++ )
      c[i] = a[i] - b[i];
}

/*---------------------------------------------------------------------
Returns the dot product of the two input vectors.
---------------------------------------------------------------------*/
double	Dot( tPointi a, tPointd b )
{
    int i;
    double sum = 0.0;

    for( i = 0; i < DIM; i++ )
       sum += a[i] * b[i];

    return  sum;
}

/*---------------------------------------------------------------------
Compute the cross product of (b-a)x(c-a) and place into N.
---------------------------------------------------------------------*/
void	NormalVec( tPointi a, tPointi b, tPointi c, tPointd N )
{
    N[X] = ( c[Z] - a[Z] ) * ( b[Y] - a[Y] ) -
           ( b[Z] - a[Z] ) * ( c[Y] - a[Y] );
    N[Y] = ( b[Z] - a[Z] ) * ( c[X] - a[X] ) -
           ( b[X] - a[X] ) * ( c[Z] - a[Z] );
    N[Z] = ( b[X] - a[X] ) * ( c[Y] - a[Y] ) -
           ( b[Y] - a[Y] ) * ( c[X] - a[X] );
}

/* Reads in the number of faces of the polyhedron and their indices from stdin,
    and returns the number n. */
int ReadFaces( void )
{
  int	i,j,k, n;
  int   w; /* temp storage for coordinate. */
  
  do {
    scanf( "%d", &n );
    if ( n <= PMAX )
      break;
    printf("Error in read_vertex:  too many points; max is %d\n", PMAX);
  }
  while ( 1 );

  printf( "Faces:\n" );
  printf( "  i   i0  i1  i2\n");
  for ( i = 0; i < n; i++ ) {
    scanf( "%d %d %d", &Faces[i][0], &Faces[i][1], &Faces[i][2] );
    printf( "%3d:%3d%4d%4d\n", i, Faces[i][0], Faces[i][1], Faces[i][2] );
    /* Compute bounding box. */
    /* Initialize to first vertex. */
    for ( j=0; j < 3; j++ ) {
       Box[i][0][j] = Vertices[ Faces[i][0] ][j];
       Box[i][1][j] = Vertices[ Faces[i][0] ][j];
    }
    /* Check k=1,2 vertices of face. */
    for ( k=1; k < 3; k++ )
    for ( j=0; j < 3; j++ ) {
       w = Vertices[ Faces[i][k] ][j];
       if ( w < Box[i][0][j] ) Box[i][0][j] = w;
       if ( w > Box[i][1][j] ) Box[i][1][j] = w;
    }
    /* printf("Bounding box: (%d,%d,%d);(%d,%d,%d)\n",
       Box[i][0][0],
       Box[i][0][1],
       Box[i][0][2],
       Box[i][1][0],
       Box[i][1][1],
       Box[i][1][2] );
    */
  }
  printf("n = %3d faces read\n",n);
  putchar('\n');

  return n;
}

/* Assumption: p lies in the plane containing T.
    Returns a char:
     'V': the query point p coincides with a Vertex of triangle T.
     'E': the query point p is in the relative interior of an Edge of triangle T.
     'F': the query point p is in the relative interior of a Face of triangle T.
     '0': the query point p does not intersect (misses) triangle T.
*/

char 	InTri3D( tPointi T, int m, tPointi p )
{
   int i;           /* Index for X,Y,Z           */
   int j;           /* Index for X,Y             */
   int k;           /* Index for triangle vertex */
   tPointi pp;      /* projected p */
   tPointi Tp[3];   /* projected T: three new vertices */

   /* Project out coordinate m in both p and the triangular face */
   j = 0;
   for ( i = 0; i < DIM; i++ ) {
     if ( i != m ) {    /* skip largest coordinate */
       pp[j] = p[i];
       for ( k = 0; k < 3; k++ )
	Tp[k][j] = Vertices[T[k]][i];
       j++;
     }
   }
   return( InTri2D( Tp, pp ) );
}

char 	InTri2D( tPointi Tp[3], tPointi pp )
{
   int area0, area1, area2;

   /* compute three AreaSign() values for pp w.r.t. each edge of the face in 2D */
   area0 = AreaSign( pp, Tp[0], Tp[1] );
   area1 = AreaSign( pp, Tp[1], Tp[2] );
   area2 = AreaSign( pp, Tp[2], Tp[0] );
   printf("area0=%d  area1=%d  area2=%d\n",area0,area1,area2);

   if ( ( area0 == 0 ) && ( area1 > 0 ) && ( area2 > 0 ) ||
        ( area1 == 0 ) && ( area0 > 0 ) && ( area2 > 0 ) ||
        ( area2 == 0 ) && ( area0 > 0 ) && ( area1 > 0 ) ) 
     return 'E';

   if ( ( area0 == 0 ) && ( area1 < 0 ) && ( area2 < 0 ) ||
        ( area1 == 0 ) && ( area0 < 0 ) && ( area2 < 0 ) ||
        ( area2 == 0 ) && ( area0 < 0 ) && ( area1 < 0 ) )
     return 'E';                 
   
   if ( ( area0 >  0 ) && ( area1 > 0 ) && ( area2 > 0 ) ||
        ( area0 <  0 ) && ( area1 < 0 ) && ( area2 < 0 ) )
     return 'F';

   if ( ( area0 == 0 ) && ( area1 == 0 ) && ( area2 == 0 ) )
     fprintf( stderr, "Error in InTriD\n" ), exit(EXIT_FAILURE);

   if ( ( area0 == 0 ) && ( area1 == 0 ) ||
        ( area0 == 0 ) && ( area2 == 0 ) ||
        ( area1 == 0 ) && ( area2 == 0 ) )
     return 'V';

   else  
     return '0';  
}

int     AreaSign( tPointi a, tPointi b, tPointi c )  
{
    double area2;

    area2 = ( b[0] - a[0] ) * (double)( c[1] - a[1] ) -
            ( c[0] - a[0] ) * (double)( b[1] - a[1] );

    /* The area should be an integer. */
    if      ( area2 >  0.5 ) return  1;
    else if ( area2 < -0.5 ) return -1;
    else                     return  0;
}                            

char    SegTriInt( tPointi T, tPointi q, tPointi r, tPointd p )
{
    int code = '?';
    int m = -1;

    code = SegPlaneInt( T, q, r, p, &m );
    printf("SegPlaneInt code=%c, m=%d; p=(%lf,%lf,%lf)\n", code,m,p[X],p[Y],p[Z]
);

    if      ( code == '0')
       return '0';
    else if ( code == 'q')
       return InTri3D( T, m, q );
    else if ( code == 'r')
       return InTri3D( T, m, r );
    else if ( code == 'p' )
       return InPlane( T, m, q, r, p );
    else if ( code == '1' )
       return SegTriCross( T, q, r );
    else /* Error */
       return code;
}

char	InPlane( tPointi T, int m, tPointi q, tPointi r, tPointd p)
{
    /* NOT IMPLEMENTED */
    return 'p';
}

/*---------------------------------------------------------------------
The signed volumes of three tetrahedra are computed, determined
by the segment qr, and each edge of the triangle.  
Returns a char:
   'v': the open segment includes a vertex of T.
   'e': the open segment includes a point in the relative interior of an edge
   of T.
   'f': the open segment includes a point in the relative interior of a face
   of T.
   '0': the open segment does not intersect triangle T.
---------------------------------------------------------------------*/

char SegTriCross( tPointi T, tPointi q, tPointi r )
{
   int vol0, vol1, vol2;
   
   vol0 = VolumeSign( q, Vertices[ T[0] ], Vertices[ T[1] ], r ); 
   vol1 = VolumeSign( q, Vertices[ T[1] ], Vertices[ T[2] ], r ); 
   vol2 = VolumeSign( q, Vertices[ T[2] ], Vertices[ T[0] ], r );
 
   printf( "SegTriCross:  vol0 = %d; vol1 = %d; vol2 = %d\n", 
      vol0, vol1, vol2 ); 
     
   /* Same sign: segment intersects interior of triangle. */
   if ( ( ( vol0 > 0 ) && ( vol1 > 0 ) && ( vol2 > 0 ) ) || 
        ( ( vol0 < 0 ) && ( vol1 < 0 ) && ( vol2 < 0 ) ) )
      return 'f';
   
   /* Opposite sign: no intersection between segment and triangle */
   if ( ( ( vol0 > 0 ) || ( vol1 > 0 ) || ( vol2 > 0 ) ) &&
        ( ( vol0 < 0 ) || ( vol1 < 0 ) || ( vol2 < 0 ) ) )
      return '0';

   else if ( ( vol0 == 0 ) && ( vol1 == 0 ) && ( vol2 == 0 ) )
     fprintf( stderr, "Error 1 in SegTriCross\n" ), exit(EXIT_FAILURE);
   
   /* Two zeros: segment intersects vertex. */
   else if ( ( ( vol0 == 0 ) && ( vol1 == 0 ) ) || 
             ( ( vol0 == 0 ) && ( vol2 == 0 ) ) || 
             ( ( vol1 == 0 ) && ( vol2 == 0 ) ) )
      return 'v';

   /* One zero: segment intersects edge. */
   else if ( ( vol0 == 0 ) || ( vol1 == 0 ) || ( vol2 == 0 ) )
      return 'e';
   
   else
     fprintf( stderr, "Error 2 in SegTriCross\n" ), exit(EXIT_FAILURE);
}

int 	VolumeSign( tPointi a, tPointi b, tPointi c, tPointi d )
{ 
   double vol;
   double ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz;
   double bxdx, bydy, bzdz, cxdx, cydy, czdz;

   ax = a[X];
   ay = a[Y];
   az = a[Z];
   bx = b[X];
   by = b[Y];
   bz = b[Z];
   cx = c[X]; 
   cy = c[Y];
   cz = c[Z];
   dx = d[X];
   dy = d[Y];
   dz = d[Z];

   bxdx=bx-dx;
   bydy=by-dy;
   bzdz=bz-dz;
   cxdx=cx-dx;
   cydy=cy-dy;
   czdz=cz-dz;
   vol =   (az-dz) * (bxdx*cydy - bydy*cxdx)
         + (ay-dy) * (bzdz*cxdx - bxdx*czdz)
         + (ax-dx) * (bydy*czdz - bzdz*cydy);


   /* The volume should be an integer. */
   if      ( vol > 0.5 )   return  1;
   else if ( vol < -0.5 )  return -1;
   else                    return  0;
}

/*
  This function returns a char:
    '0': the segment [ab] does not intersect (completely misses) the 
         bounding box surrounding the n-th triangle T.  It lies
         strictly to one side of one of the six supporting planes.
    '?': status unknown: the segment may or may not intersect T.
*/
char BoxTest ( int n, tPointi a, tPointi b )
{
   int i; /* Coordinate index */
   int w;

   for ( i=0; i < DIM; i++ ) {
       w = Box[ n ][0][i]; /* min: lower left */
       if ( (a[i] < w) && (b[i] < w) ) return '0';
       w = Box[ n ][1][i]; /* max: upper right */
       if ( (a[i] > w) && (b[i] > w) ) return '0';
   }
   return '?';
}