chull.cpp

求空间点集的三维凸包。chull.cpp, chull.h。来自计算几何英文原著作者个人网站。

#include <stdio.h>
#include <assert.h>
#include <stdlib.h>

#include "chull.h"

#define SWAP(t,x,y) { t = x; x = y; y = t; }
#define ADD( head, p )  if ( head )  { \
				p->next = head; \
				p->prev = head->prev; \
				head->prev = p; \
				p->prev->next = p; \
			} \
			else { \
				head = p; \
				head->next = head->prev = p; \
			}

#define DELETE( head, p ) if ( head )  { \
				if ( head == head->next ) \
					head = NULL;  \
				else if ( p == head ) \
					head = head->next; \
				p->next->prev = p->prev;  \
				p->prev->next = p->next; \
				delete p; \
			}

Chull3D::Chull3D (float *v, int n)
{
  int i;
  assert (v);

  /* init vertices */
  vertices = NULL;
  for (i=0; i<n; i++)
    add_vertex (new Chull3D_vertex (v[3*i], v[3*i+1], v[3*i+2]));
  edges = NULL;
  faces = NULL;
}

Chull3D::~Chull3D ()
{
  Chull3D_vertex *v, *vt;
  Chull3D_edge *e, *et;
  Chull3D_face *f, *ft;

  // clean vertices
  v = vertices;
  do {
	  vt = v;
	  v = v->next;
	  delete_vertex (vt);
  } while (vertices->next != vertices);
  delete_vertex (vertices);

  // clean edges
  e = edges;
  do {
	  et = e;
	  e = e->next;
	  delete_edge (et);
  } while (edges->next != edges);
  delete_edge (edges);

  // clean faces
  f = faces;
  do {
	  ft = f;
	  f = f->next;
	  delete_face (ft);
  } while (faces->next != faces);
  delete_face (faces);
}

void Chull3D::add_vertex    (Chull3D_vertex *v) { ADD (vertices, v); }
void Chull3D::add_edge      (Chull3D_edge *e)   { ADD (edges, e); }
void Chull3D::add_face      (Chull3D_face *f)   { ADD (faces, f); }
void Chull3D::delete_vertex (Chull3D_vertex *v) { DELETE (vertices, v); }
void Chull3D::delete_edge   (Chull3D_edge *e)   { DELETE (edges, e); }
void Chull3D::delete_face   (Chull3D_face *f)   { DELETE (faces, f); }

/**********/
/* output */
/**********/
int
Chull3D::get_n_vertices (void)
{
  Chull3D_vertex *v = vertices;
  int n=0;
  if (!vertices)
    return 0;
  else
    do { v = v->next; n++; } while (v != vertices);
  return n;
}

int
Chull3D::get_n_faces (void)
{
  Chull3D_face *f = faces;
  int n=0;
  if (!faces)
    return 0;
  else
    do { f = f->next; n++; } while (f != faces);
  return n;
}

int
Chull3D::get_vertex_index (Chull3D_vertex *v)
{
  Chull3D_vertex *v_walk = vertices;
  int index = 0;
  do {
    if (v_walk == v)
      return index;
    index++;
    v_walk = v_walk->next;
  } while (v_walk != vertices);
  return -1;
}

int
Chull3D::get_convex_hull (float **v, int *nv, int **f, int *nf)
{
  *nv = get_n_vertices ();
  *nf = get_n_faces ();
  
  /* memory allocation */
  float *cv_vertices = (float*)malloc(3*(*nv)*sizeof(float));
  int   *cv_faces    = (int*)malloc(3*(*nf)*sizeof(int));
  if (!cv_vertices || !cv_faces)
    {
      v = NULL; f = NULL;
      *nv = *nf = 0;
      return 0;
    }
  
  /* vertices */
  Chull3D_vertex *v_walk=vertices;
  int i=0;
  do {
    cv_vertices[3*i]   = v_walk->pt[0];
    cv_vertices[3*i+1] = v_walk->pt[1];
    cv_vertices[3*i+2] = v_walk->pt[2];
    v_walk = v_walk->next;
    i++;
  } while (v_walk != vertices);

  /* faces */
  i = 0;
  Chull3D_face *f_walk = faces;
  do {
    cv_faces[3*i]   = get_vertex_index (f_walk->vertices[0]);
    cv_faces[3*i+1] = get_vertex_index (f_walk->vertices[1]);
    cv_faces[3*i+2] = get_vertex_index (f_walk->vertices[2]);
    f_walk = f_walk->next;
    i++;
  } while (f_walk != faces);
  
  *v = cv_vertices;
  *f = cv_faces;

  return 1;
}

void
Chull3D::export_obj (char *filename)
{
  FILE *ptr;
  ptr = fopen (filename, "w");
  if (!ptr) return;

  /* vertices */
  Chull3D_vertex *v_walk = vertices;
  do {
    fprintf (ptr, "v %f %f %f\n", v_walk->pt[0], v_walk->pt[1], v_walk->pt[2]);
    v_walk = v_walk->next;
  } while (v_walk != vertices);

  /* faces */
  Chull3D_face *f_walk = faces;
  do {
    fprintf (ptr, "f %d %d %d\n",
	     get_vertex_index (f_walk->vertices[0])+1,
	     get_vertex_index (f_walk->vertices[1])+1,
	     get_vertex_index (f_walk->vertices[2])+1);
    f_walk = f_walk->next;
  } while (f_walk != faces);

  fclose (ptr);
}

/*****************/
/*** Algorithm ***/
/*****************/
void
Chull3D::compute (void)
{
  double_triangle ();
  construct_hull ();
}

/* builds the initial double triangle */
int
Chull3D::double_triangle (void)
{
  Chull3D_vertex *v0, *v1, *v2, *v3;
  int vol;

  /* find 3 non collinear points */
  v0 = vertices;
  while (are_collinear (v0, v0->next, v0->next->next))
    if ( (v0=v0->next) == vertices)
      {
	printf ("All the vertices are collinear\n");
	return 1;
      }
  v1 = v0->next;
  v2 = v1->next;

  /* mark the vertices as processed */
  v0->processed = 1;
  v1->processed = 1;
  v2->processed = 1;

  /* create the two "twins" faces */
  Chull3D_face *f0, *f1 = NULL;
  f0 = new Chull3D_face (this, v0, v1, v2, f1);
  f1 = new Chull3D_face (this, v2, v1, v0, f0);

  /* link adjacent face fields */
  f0->edges[0]->adj_faces[1] = f1;
  f0->edges[1]->adj_faces[1] = f1;
  f0->edges[2]->adj_faces[1] = f1;
  f1->edges[0]->adj_faces[1] = f0;
  f1->edges[1]->adj_faces[1] = f0;
  f1->edges[2]->adj_faces[1] = f0;

  /* find a fourth, non coplanar point to form tetrahedron */
  v3 = v2->next;
  vol = volume_sign (f0, v3);
  while (!vol)
    {
      if ( (v3=v3->next) == v0)
	{
	  printf ("All the vertices are coplanar\n");
	  return 1;
	}
      vol = volume_sign (f0, v3);
    }

  /* insure that v3 will be the first added */
  vertices = v3;
  
  return 0;
}

/*
 * construct_hull adds the vertices to the hull one at a time.
 */
int
Chull3D::construct_hull (void)
{
  Chull3D_vertex *v, *vnext;

  v = vertices;
  do {
    vnext = v->next;
    if (!v->processed)
      {
	v->processed = 1;
	add_one (v);
	clean_up (vnext);
      }
    v = vnext;
  } while (v != vertices);
  
  return 0;
}

/*
 * are_collinear checks to see if the three points given are
 * collinear by checking to see if each element of the cross
 * product is zero.
 */
int
Chull3D::are_collinear (Chull3D_vertex *v1, Chull3D_vertex *v2, Chull3D_vertex *v3)
{
  return
    (( v3->pt[2] - v1->pt[2] ) * ( v2->pt[1] - v1->pt[1] ) -
     ( v2->pt[2] - v1->pt[2] ) * ( v3->pt[1] - v1->pt[1] ) == 0
     && ( v2->pt[2] - v1->pt[2] ) * ( v3->pt[0] - v1->pt[0] ) -
     ( v2->pt[0] - v1->pt[0] ) * ( v3->pt[2] - v1->pt[2] ) == 0
     && ( v2->pt[0] - v1->pt[0] ) * ( v3->pt[1] - v1->pt[1] ) -
     ( v2->pt[1] - v1->pt[1] ) * ( v3->pt[0] - v1->pt[0] ) == 0);
}

/*
 * volume_sign returns the sign of the volume of the tetrahedron
 * determined by f and p.
 * Volume_sign is +1 iff p is on the negative side of f, where the
 * positive side is determined by the rh-rule. So the volume is
 * positive if the ccw normal to f points outside the tetrahedron.
 * The final fewer-multiplications form is due to Bob Williamson.
 */
int
Chull3D::volume_sign (Chull3D_face *f, Chull3D_vertex *v)
{
  float ax, ay, az, bx, by, bz, cx, cy, cz;
  float vol;

   ax = f->vertices[0]->pt[0] - v->pt[0];
   ay = f->vertices[0]->pt[1] - v->pt[1];
   az = f->vertices[0]->pt[2] - v->pt[2];
   bx = f->vertices[1]->pt[0] - v->pt[0];
   by = f->vertices[1]->pt[1] - v->pt[1];
   bz = f->vertices[1]->pt[2] - v->pt[2];
   cx = f->vertices[2]->pt[0] - v->pt[0];
   cy = f->vertices[2]->pt[1] - v->pt[1];
   cz = f->vertices[2]->pt[2] - v->pt[2];

   vol =   ax * (by*cz - bz*cy)
         + ay * (bz*cx - bx*cz)
         + az * (bx*cy - by*cx);

   if      ( vol >  0.0 ) return  1;
   else if ( vol <  0.0 ) return -1;
   else                   return  0;
}