chull.c

凸壳算法

/*
This code is described in "Computational Geometry in C" (Second Edition),
Chapter 4.  It is not written to be comprehensible without the 
explanation in that book.

Input: 3n integer coordinates for the points.
Output: the 3D convex hull, in postscript with embedded comments
        showing the vertices and faces.

Compile: gcc -o chull chull.c (or simply: make)

Written by Joseph O'Rourke, with contributions by 
  Kristy Anderson, John Kutcher, Catherine Schevon, Susan Weller.
Last modified: May 2000
Questions to [e-mail].

--------------------------------------------------------------------
This code is Copyright 2000 by Joseph O'Rourke.  It may be freely 
redistributed in its entirety provided that this copyright notice is 
not removed.
--------------------------------------------------------------------
*/
#include <stdio.h>
#include <math.h>

/*Define Boolean type */
typedef	enum { FALSE, TRUE }	bool;

/* Define vertex indices. */
#define X   0
#define Y   1
#define Z   2

/* Define structures for vertices, edges and faces */
typedef struct tVertexStructure tsVertex;
typedef tsVertex *tVertex;

typedef struct tEdgeStructure tsEdge;
typedef tsEdge *tEdge;

typedef struct tFaceStructure tsFace;
typedef tsFace *tFace;

struct tVertexStructure {
   int      v[3];
   int	    vnum;
   tEdge    duplicate;	        /* pointer to incident cone edge (or NULL) */
   bool     onhull;		/* T iff point on hull. */
   bool	    mark;		/* T iff point already processed. */
   tVertex  next, prev;
};

struct tEdgeStructure {
   tFace    adjface[2];
   tVertex  endpts[2];
   tFace    newface;            /* pointer to incident cone face. */
   bool     delete;		/* T iff edge should be delete. */
   tEdge    next, prev;
};

struct tFaceStructure {
   tEdge    edge[3];
   tVertex  vertex[3];
   bool	    visible;	        /* T iff face visible from new point. */
   tFace    next, prev;
};

/* Define flags */
#define ONHULL   	TRUE
#define REMOVED  	TRUE
#define VISIBLE  	TRUE
#define PROCESSED	TRUE
#define SAFE		1000000		/* Range of safe coord values. */

/* Global variable definitions */
tVertex vertices = NULL;
tEdge edges    	 = NULL;
tFace faces    	 = NULL;
bool debug = FALSE;
bool check = FALSE;

/* Function declarations */
tVertex MakeNullVertex( void );
void    ReadVertices( void );
void    Print( void );
void    SubVec( int a[3], int b[3], int c[3]);
void    DoubleTriangle( void );
void    ConstructHull( void );
bool	AddOne( tVertex p );
int     VolumeSign(tFace f, tVertex p);
int 	Volumei( tFace f, tVertex p );
tFace	MakeConeFace( tEdge e, tVertex p );
void    MakeCcw( tFace f, tEdge e, tVertex p );
tEdge   MakeNullEdge( void );
tFace   MakeNullFace( void );
tFace   MakeFace( tVertex v0, tVertex v1, tVertex v2, tFace f );
void    CleanUp( tVertex *pvnext );
void    CleanEdges( void );
void    CleanFaces( void );
void    CleanVertices( tVertex *pvnext );
bool	Collinear( tVertex a, tVertex b, tVertex c );
void    CheckEuler(int V, int E, int F );
void	PrintPoint( tVertex p );
void    Checks( void );
void	Consistency( void );
void	Convexity( void );
void	PrintOut( tVertex v );
void	PrintVertices( void );
void	PrintEdges( void );
void	PrintFaces( void );
void	CheckEndpts ( void );
void	EdgeOrderOnFaces ( void );

#include "macros.h"

/*-------------------------------------------------------------------*/
main( int argc, char *argv[] )
{

  if ( argc > 1 && argv[1][0] == '-' ) {
    if( argv[1][1] ==  'd' ) {
      debug = TRUE;
      check = TRUE;
      fprintf( stderr, "Debug and check mode\n");
    }
    if( argv[1][1] == 'c' ) {
      check = TRUE;
      fprintf( stderr, "Check mode\n");
    }
  }
  else if ( argc > 1 && argv[1][0] != '-' ) {
    printf ("Usage:  %s -d[ebug] c[heck]\n", *argv );
    printf ("x y z coords of vertices from stdin\n");
    exit(1);
  }
  
   ReadVertices();
   DoubleTriangle();
   ConstructHull();
   EdgeOrderOnFaces();
   Print();
}

/*---------------------------------------------------------------------
MakeNullVertex: Makes a vertex, nulls out fields.
---------------------------------------------------------------------*/
tVertex	MakeNullVertex( void )
{
   tVertex  v;
   
   NEW( v, tsVertex );
   v->duplicate = NULL;
   v->onhull = !ONHULL;
   v->mark = !PROCESSED;
   ADD( vertices, v );

   return v;
}

/*---------------------------------------------------------------------
ReadVertices: Reads in the vertices, and links them into a circular
list with MakeNullVertex.  There is no need for the # of vertices to be
the first line: the function looks for EOF instead.  Sets the global
variable vertices via the ADD macro.
---------------------------------------------------------------------*/
void	ReadVertices( void )
{
   tVertex  v;
   int      x, y, z;
   int	    vnum = 0;

   while ( scanf ("%d %d %d", &x, &y, &z ) != EOF )  {
      v = MakeNullVertex();
      v->v[X] = x;
      v->v[Y] = y;
      v->v[Z] = z;
      v->vnum = vnum++;
      if ( ( abs(x) > SAFE ) || ( abs(y) > SAFE ) || ( abs(z) > SAFE ) ) {
         printf("Coordinate of vertex below might be too large: run with -d flag\n");
         PrintPoint(v);
      }
   }
}

/*---------------------------------------------------------------------
Print: Prints out the vertices and the faces.  Uses the vnum indices 
corresponding to the order in which the vertices were input.
Output is in PostScript format.
---------------------------------------------------------------------*/
void	Print( void )
{
   /* Pointers to vertices, edges, faces. */
   tVertex  v;
   tEdge    e;
   tFace    f;
   int xmin, ymin, xmax, ymax;
   int a[3], b[3];  /* used to compute normal vector */
   /* Counters for Euler's formula. */
   int 	V = 0, E = 0 , F = 0;
   /* Note: lowercase==pointer, uppercase==counter. */

   /*-- find X min & max --*/
   v = vertices;
   xmin = xmax = v->v[X];
   do {
      if( v->v[X] > xmax ) xmax = v->v[X];
      else
	 if( v->v[X] < xmin ) xmin = v->v[X];
      v = v->next;
   } while ( v != vertices );
	
   /*-- find Y min & max --*/
   v = vertices;
   ymin = ymax = v->v[Y];
   do {
      if( v->v[Y] > ymax ) ymax = v->v[Y];
      else
	 if( v->v[Y] < ymin ) ymin = v->v[Y];
      v = v->next;
   } while ( v != vertices );
	
   /* PostScript header */
   printf("%%!PS\n");
   printf("%%%%BoundingBox: %d %d %d %d\n", 
	  xmin, ymin, xmax, ymax);
   printf(".00 .00 setlinewidth\n");
   printf("%d %d translate\n", -xmin+72, -ymin+72 );
   /* The +72 shifts the figure one inch from the lower left corner */

   /* Vertices. */
   v = vertices;
   do {                                 
      if( v->mark ) V++;           
      v = v->next;
   } while ( v != vertices );
   printf("\n%%%% Vertices:\tV = %d\n", V);
   printf("%%%% index:\tx\ty\tz\n");
   do {                                 
      printf( "%%%% %5d:\t%d\t%d\t%d\n", 
	     v->vnum, v->v[X], v->v[Y], v->v[Z] );
      v = v->next;
   } while ( v != vertices );
	
   /* Faces. */
   /* visible faces are printed as PS output */
   f = faces;
   do {
      ++F;                              
      f  = f ->next;
   } while ( f  != faces );
   printf("\n%%%% Faces:\tF = %d\n", F );
   printf("%%%% Visible faces only: \n");
   do {           
      /* Print face only if it is visible: if normal vector >= 0 */
      SubVec( f->vertex[1]->v, f->vertex[0]->v, a );
      SubVec( f->vertex[2]->v, f->vertex[1]->v, b );	  
      if(( a[0] * b[1] - a[1] * b[0] ) >= 0 )
      {
	 printf("%%%% vnums:  %d  %d  %d\n", 
		f->vertex[0]->vnum, 
		f->vertex[1]->vnum, 
		f->vertex[2]->vnum);
	 printf("newpath\n");
	 printf("%d\t%d\tmoveto\n", 
		f->vertex[0]->v[X], f->vertex[0]->v[Y] );
	 printf("%d\t%d\tlineto\n", 
		f->vertex[1]->v[X], f->vertex[1]->v[Y] );
	 printf("%d\t%d\tlineto\n", 
		f->vertex[2]->v[X], f->vertex[2]->v[Y] );
	 printf("closepath stroke\n\n");
      }
      f = f->next;
   } while ( f != faces );

   /* prints a list of all faces */
   printf("%%%% List of all faces: \n");
   printf("%%%%\tv0\tv1\tv2\t(vertex indices)\n");
   do {
      printf("%%%%\t%d\t%d\t%d\n",
	     f->vertex[0]->vnum,
	     f->vertex[1]->vnum,
	     f->vertex[2]->vnum );
      f = f->next;
   } while ( f != faces );
	
   /* Edges. */	
   e = edges;
   do {
      E++;
      e = e->next;
   } while ( e != edges );
   printf("\n%%%% Edges:\tE = %d\n", E );
   /* Edges not printed out (but easily added). */

   printf("\nshowpage\n\n");

   check = TRUE;
   CheckEuler( V, E, F );

}

/*---------------------------------------------------------------------
SubVec:  Computes a - b and puts it into c.
---------------------------------------------------------------------*/
void    SubVec( int a[3], int b[3], int c[3])
{
   int  i;

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

}

/*---------------------------------------------------------------------
 DoubleTriangle builds the initial double triangle.  It first finds 3 
 noncollinear points and makes two faces out of them, in opposite order.
 It then finds a fourth point that is not coplanar with that face.  The  
 vertices are stored in the face structure in counterclockwise order so 
 that the volume between the face and the point is negative. Lastly, the
 3 newfaces to the fourth point are constructed and the data structures
 are cleaned up. 
---------------------------------------------------------------------*/
void    DoubleTriangle( void )
{
   tVertex  v0, v1, v2, v3, t;
   tFace    f0, f1 = NULL;
   tEdge    e0, e1, e2, s;
   int      vol;
	
   /* Find 3 noncollinear points. */
   v0 = vertices;
   while ( Collinear( v0, v0->next, v0->next->next ) )
      if ( ( v0 = v0->next ) == vertices )
         printf("DoubleTriangle:  All points are Collinear!\n"), exit(0);
   v1 = v0->next;
   v2 = v1->next;
	
   /* Mark the vertices as processed. */
   v0->mark = PROCESSED;
   v1->mark = PROCESSED;
   v2->mark = PROCESSED;
   
   /* Create the two "twin" faces. */
   f0 = MakeFace( v0, v1, v2, f1 );
   f1 = MakeFace( v2, v1, v0, f0 );

   /* Link adjacent face fields. */
   f0->edge[0]->adjface[1] = f1;
   f0->edge[1]->adjface[1] = f1;
   f0->edge[2]->adjface[1] = f1;
   f1->edge[0]->adjface[1] = f0;
   f1->edge[1]->adjface[1] = f0;
   f1->edge[2]->adjface[1] = f0;
	
   /* Find a fourth, noncoplanar point to form tetrahedron. */
   v3 = v2->next;
   vol = VolumeSign( f0, v3 );
   while ( !vol )   {
      if ( ( v3 = v3->next ) == v0 ) 
         printf("DoubleTriangle:  All points are coplanar!\n"), exit(0);
      vol = VolumeSign( f0, v3 );
   }
	
   /* Insure that v3 will be the first added. */
   vertices = v3;
   if ( debug ) {
      fprintf(stderr, "DoubleTriangle: finished. Head repositioned at v3.\n");
      PrintOut( vertices );
   }

	
}