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, determinedby 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 '?';}/* irint not available in some libraries, so... */int	irint( double x ){	return (int) rint( x );}