chull.c

凸壳算法

   while ( faces && faces->visible ) { 
      f = faces;
      DELETE( faces, f );
   }
   f = faces->next;
   do {
      if ( f->visible ) {
	 t = f;
	 f = f->next;
	 DELETE( faces, t );
      }
      else f = f->next;
   } while ( f != faces );
}

/*---------------------------------------------------------------------
CleanVertices runs through the vertex list and deletes the 
vertices that are marked as processed but are not incident to any 
undeleted edges. 
The pointer to vnext, pvnext, is used to alter vnext in
ConstructHull() if we are about to delete vnext.
---------------------------------------------------------------------*/
void	CleanVertices( tVertex *pvnext )
{
   tEdge    e;
   tVertex  v, t;
	
   /* Mark all vertices incident to some undeleted edge as on the hull. */
   e = edges;
   do {
      e->endpts[0]->onhull = e->endpts[1]->onhull = ONHULL;
      e = e->next;
   } while (e != edges);
	
   /* Delete all vertices that have been processed but
      are not on the hull. */
   while ( vertices && vertices->mark && !vertices->onhull ) { 
      /* If about to delete vnext, advance it first. */
      if ( v == *pvnext )
         *pvnext = v->next;
      v = vertices;
      DELETE( vertices, v );
   }
   v = vertices->next;
   do {
      if ( v->mark && !v->onhull ) {    
	 t = v; 
	 v = v->next;
	 DELETE( vertices, t )
      }
      else v = v->next;
   } while ( v != vertices );
	
   /* Reset flags. */
   v = vertices;
   do {
      v->duplicate = NULL; 
      v->onhull = !ONHULL; 
      v = v->next;
   } while ( v != vertices );
}

/*---------------------------------------------------------------------
Collinear checks to see if the three points given are collinear,
by checking to see if each element of the cross product is zero.
---------------------------------------------------------------------*/
bool	Collinear( tVertex a, tVertex b, tVertex c )
{
   return 
         ( c->v[Z] - a->v[Z] ) * ( b->v[Y] - a->v[Y] ) -
         ( b->v[Z] - a->v[Z] ) * ( c->v[Y] - a->v[Y] ) == 0
      && ( b->v[Z] - a->v[Z] ) * ( c->v[X] - a->v[X] ) -
         ( b->v[X] - a->v[X] ) * ( c->v[Z] - a->v[Z] ) == 0
      && ( b->v[X] - a->v[X] ) * ( c->v[Y] - a->v[Y] ) -
         ( b->v[Y] - a->v[Y] ) * ( c->v[X] - a->v[X] ) == 0  ;
}

/*---------------------------------------------------------------------
Consistency runs through the edge list and checks that all
adjacent faces have their endpoints in opposite order.  This verifies
that the vertices are in counterclockwise order.
---------------------------------------------------------------------*/
void	Consistency( void )
{
   register tEdge  e;
   register int    i, j;

   e = edges;

   do {
      /* find index of endpoint[0] in adjacent face[0] */
      for ( i = 0; e->adjface[0]->vertex[i] != e->endpts[0]; ++i )
	 ;
   
      /* find index of endpoint[0] in adjacent face[1] */
      for ( j = 0; e->adjface[1]->vertex[j] != e->endpts[0]; ++j )
	 ;

      /* check if the endpoints occur in opposite order */
      if ( !( e->adjface[0]->vertex[ (i+1) % 3 ] ==
	      e->adjface[1]->vertex[ (j+2) % 3 ] ||
	      e->adjface[0]->vertex[ (i+2) % 3 ] ==
	      e->adjface[1]->vertex[ (j+1) % 3 ] )  )
	 break;
      e = e->next;

   } while ( e != edges );

   if ( e != edges )
      fprintf( stderr, "Checks: edges are NOT consistent.\n");
   else
      fprintf( stderr, "Checks: edges consistent.\n");

}

/*---------------------------------------------------------------------
Convexity checks that the volume between every face and every
point is negative.  This shows that each point is inside every face
and therefore the hull is convex.
---------------------------------------------------------------------*/
void	Convexity( void )
{
   register tFace    f;
   register tVertex  v;
   int               vol;

   f = faces;
   
   do {
      v = vertices;
      do {
	 if ( v->mark ) {
	    vol = VolumeSign( f, v );
	    if ( vol < 0 )
	       break;
	 }
	 v = v->next;
      } while ( v != vertices );

      f = f->next;

   } while ( f != faces );

   if ( f != faces )
      fprintf( stderr, "Checks: NOT convex.\n");
   else if ( check ) 
      fprintf( stderr, "Checks: convex.\n");
}

/*---------------------------------------------------------------------
CheckEuler checks Euler's relation, as well as its implications when
all faces are known to be triangles.  Only prints positive information
when debug is true, but always prints negative information.
---------------------------------------------------------------------*/
void	CheckEuler( int V, int E, int F )
{
   if ( check )
      fprintf( stderr, "Checks: V, E, F = %d %d %d:\t", V, E, F);

   if ( (V - E + F) != 2 )
      fprintf( stderr, "Checks: V-E+F != 2\n");
   else if ( check )
      fprintf( stderr, "V-E+F = 2\t");


   if ( F != (2 * V - 4) )
      fprintf( stderr, "Checks: F=%d != 2V-4=%d; V=%d\n",
	      F, 2*V-4, V);
   else if ( check ) 
      fprintf( stderr, "F = 2V-4\t");
   
   if ( (2 * E) != (3 * F) )
      fprintf( stderr, "Checks: 2E=%d != 3F=%d; E=%d, F=%d\n",
	      2*E, 3*F, E, F );
   else if ( check ) 
      fprintf( stderr, "2E = 3F\n");
}

/*-------------------------------------------------------------------*/
void	Checks( void )
{
   tVertex  v;
   tEdge    e;
   tFace    f;
   int 	   V = 0, E = 0 , F = 0;

   Consistency();
   Convexity();
   if ( v = vertices )
      do {
         if (v->mark) V++;
	 v = v->next;
      } while ( v != vertices );
   if ( e = edges )
      do {
         E++;
	 e = e->next;
      } while ( e != edges );
   if ( f = faces )
      do {
         F++;
	 f  = f ->next;
      } while ( f  != faces );
   CheckEuler( V, E, F );
   CheckEndpts();
}


/*===================================================================
These functions are used whenever the debug flag is set.
They print out the entire contents of each data structure.  
Printing is to standard error.  To grab the output in a file in the csh, 
use this:
	chull < i.file >&! o.file
=====================================================================*/
/*-------------------------------------------------------------------*/
void	PrintOut( tVertex v )
{
   fprintf( stderr, "\nHead vertex %d = %6x :\n", v->vnum, v );
   PrintVertices();
   PrintEdges();
   PrintFaces();
}

/*-------------------------------------------------------------------*/
void	PrintVertices( void )
{
   tVertex  temp;

   temp = vertices;
   fprintf (stderr, "Vertex List\n");
   if (vertices) do {
      fprintf(stderr,"  addr %6x\t", vertices );
      fprintf(stderr,"  vnum %4d", vertices->vnum );
      fprintf(stderr,"   (%6d,%6d,%6d)",vertices->v[X],
	      vertices->v[Y], vertices->v[Z] );
      fprintf(stderr,"   active:%3d", vertices->onhull );
      fprintf(stderr,"   dup:%5x", vertices->duplicate );
      fprintf(stderr,"   mark:%2d\n", vertices->mark );
      vertices = vertices->next;
   } while ( vertices != temp );

}

/*-------------------------------------------------------------------*/
void	PrintEdges( void )
{
   tEdge  temp;
   int 	  i;
	
   temp = edges;
   fprintf (stderr, "Edge List\n");
   if (edges) do {
      fprintf( stderr, "  addr: %6x\t", edges );
      fprintf( stderr, "adj: ");
      for (i=0; i<2; ++i) 
	 fprintf( stderr, "%6x", edges->adjface[i] );
      fprintf( stderr, "  endpts:");
      for (i=0; i<2; ++i) 
	 fprintf( stderr, "%4d", edges->endpts[i]->vnum);
      fprintf( stderr, "  del:%3d\n", edges->delete );
      edges = edges->next; 
   } while (edges != temp );

}

/*-------------------------------------------------------------------*/
void	PrintFaces( void )
{
   int 	  i;
   tFace  temp;

   temp = faces;
   fprintf (stderr, "Face List\n");
   if (faces) do {
      fprintf(stderr, "  addr: %10x  ", faces );
      fprintf(stderr, "  edges:");
      for( i=0; i<3; ++i )
	 fprintf(stderr, "%10x ", faces->edge[i] );
      fprintf(stderr, "  vert:");
      for ( i=0; i<3; ++i)
	 fprintf(stderr, "%4d", faces->vertex[i]->vnum );
      fprintf(stderr, "  vis: %d\n", faces->visible );
      faces= faces->next;
   } while ( faces != temp );

}

/*-------------------------------------------------------------------
Checks that, for each face, for each i={0,1,2}, the [i]th vertex of
that face is either the [0]th or [1]st endpoint of the [ith] edge of
the face.
-------------------------------------------------------------------*/
void	CheckEndpts ( void )
{
   int 	   i;
   tFace   fstart;
   tEdge   e;
   tVertex v;
   bool error = FALSE;

   fstart = faces;
   if (faces) do {
      for( i=0; i<3; ++i ) {
         v = faces->vertex[i];
         e = faces->edge[i];
         if ( v != e->endpts[0] && v != e->endpts[1] ) {
            error = TRUE;
            fprintf(stderr,"CheckEndpts: Error!\n");
            fprintf(stderr,"  addr: %8x;", faces );
            fprintf(stderr,"  edges:");
            fprintf(stderr,"(%3d,%3d)", 
               e->endpts[0]->vnum,
               e->endpts[1]->vnum);
            fprintf(stderr,"\n");
         }
      }
      faces= faces->next;
   } while ( faces != fstart );

   if ( error )
     fprintf(stderr,"Checks: ERROR found and reported above.\n");
   else
     fprintf(stderr,"Checks: All endpts of all edges of all faces check.\n");

}

/*------------------------------------------------------------------
  EdgeOrderOnFaces: puts e0 between v0 and v1, e1 between v1 and v2,
  e2 between v2 and v0 on each face.  This should be unnecessary, alas.
  Not used in code, but useful for other purposes.
------------------------------------------------------------------*/
void    EdgeOrderOnFaces ( void ) {
  tFace f = faces;
  tEdge new;
  int i,j;

  do {
    for (i = 0; i < 3; i++) {
      if (!(((f->edge[i]->endpts[0] == f->vertex[i]) &&
             (f->edge[i]->endpts[1] == f->vertex[(i+1)%3])) ||
            ((f->edge[i]->endpts[1] == f->vertex[i]) &&
             (f->edge[i]->endpts[0] == f->vertex[(i+1)%3])))) {
        /* Change the order of the edges on the face: */
        for (j = 0; j < 3; j ++) {
          /* find the edge that should be there */
          if (((f->edge[j]->endpts[0] == f->vertex[i]) &&
               (f->edge[j]->endpts[1] == f->vertex[(i+1)%3])) ||
              ((f->edge[j]->endpts[1] == f->vertex[i]) &&
               (f->edge[j]->endpts[0] == f->vertex[(i+1)%3]))) {
            /* Swap it with the one erroneously put into its place: */
            if ( debug )
            fprintf(stderr,
              "Making a swap in EdgeOrderOnFaces: F(%d,%d,%d): e[%d] and e[%d]\n",
                    f->vertex[0]->vnum,
                    f->vertex[1]->vnum,
                    f->vertex[2]->vnum,
                    i, j);
            new = f->edge[i];
            f->edge[i] = f->edge[j];
            f->edge[j] = new;
          }
        }
      }
    }
    f = f->next;
  } while (f != faces);

}