eu.mihosoft.vrl.v3d.ext.quickhull3d.QuickHull3D Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of JavaCad Show documentation
A Java based CSG Cad library
There is a newer version: 2.0.0
/**
  * Copyright John E. Lloyd, 2004. All rights reserved. Permission to use,
  * copy, modify and redistribute is granted, provided that this copyright
  * notice is retained and the author is given credit whenever appropriate.
  *
  * This  software is distributed "as is", without any warranty, including 
  * any implied warranty of merchantability or fitness for a particular
  * use. The author assumes no responsibility for, and shall not be liable
  * for, any special, indirect, or consequential damages, or any damages
  * whatsoever, arising out of or in connection with the use of this
  * software.
  */
package eu.mihosoft.vrl.v3d.ext.quickhull3d;

import java.util.*;
import java.io.*;

// TODO: Auto-generated Javadoc
/**
 * Computes the convex hull of a set of three dimensional points.
 *
 *  The algorithm is a three dimensional implementation of Quickhull, as
 * described in Barber, Dobkin, and Huhdanpaa,  ``The Quickhull
 * Algorithm for Convex Hulls''  (ACM Transactions on Mathematical Software,
 * Vol. 22, No. 4, December 1996), and has a complexity of O(n log(n)) with
 * respect to the number of points. A well-known C implementation of Quickhull
 * that works for arbitrary dimensions is provided by qhull .
 *
 *  A hull is constructed by providing a set of points
 * to either a constructor or a
 * {@link #build(Point3d[]) build} method. After
 * the hull is built, its vertices and faces can be retrieved
 * using {@link #getVertices()
 * getVertices} and {@link #getFaces() getFaces}.
 * A typical usage might look like this:
 *  
 *   // x y z coordinates of 6 points
 *   Point3d[] points = new Point3d[] 
 *    { new Point3d (0.0,  0.0,  0.0),
 *      new Point3d (1.0,  0.5,  0.0),
 *      new Point3d (2.0,  0.0,  0.0),
 *      new Point3d (0.5,  0.5,  0.5),
 *      new Point3d (0.0,  0.0,  2.0),
 *      new Point3d (0.1,  0.2,  0.3),
 *      new Point3d (0.0,  2.0,  0.0),
 *    };
 *
 *   QuickHull3D hull = new QuickHull3D();
 *   hull.build (points);
 *
 *   System.out.println ("Vertices:");
 *   Point3d[] vertices = hull.getVertices();
 *   for (int i = 0; i < vertices.length; i++)
 *    { Point3d pnt = vertices[i];
 *      System.out.println (pnt.x + " " + pnt.y + " " + pnt.z);
 *    }
 *
 *   System.out.println ("Faces:");
 *   int[][] faceIndices = hull.getFaces();
 *   for (int i = 0; i < faceIndices.length; i++)
 *    { for (int k = 0; k < faceIndices[i].length; k++)
 *       { System.out.print (faceIndices[i][k] + " ");
 *       }
 *      System.out.println ("");
 *    }
 *  
 * As a convenience, there are also {@link #build(double[]) build}
 * and {@link #getVertices(double[]) getVertex} methods which
 * pass point information using an array of doubles.
 *
 * Robustness
 Because this algorithm uses floating
 * point arithmetic, it is potentially vulnerable to errors arising from
 * numerical imprecision.  We address this problem in the same way as qhull , by merging faces whose edges are not
 * clearly convex. A face is convex if its edges are convex, and an edge is
 * convex if the centroid of each adjacent plane is clearly  below  the
 * plane of the other face. The centroid is considered below a plane if its
 * distance to the plane is less than the negative of a {@link
 * #getDistanceTolerance() distance tolerance}.  This tolerance represents the
 * smallest distance that can be reliably computed within the available numeric
 * precision. It is normally computed automatically from the point data,
 * although an application may {@link #setExplicitDistanceTolerance set this
 * tolerance explicitly}.
 *
 *  Numerical problems are more likely to arise in situations where data
 * points lie on or within the faces or edges of the convex hull. We have
 * tested QuickHull3D for such situations by computing the convex hull of a
 * random point set, then adding additional randomly chosen points which lie
 * very close to the hull vertices and edges, and computing the convex
 * hull again. The hull is deemed correct if {@link #check check} returns
 * true.  These tests have been successful for a large number of
 * trials and so we are confident that QuickHull3D is reasonably robust.
 *
 * Merged Faces The merging of faces means that the faces returned by
 * QuickHull3D may be convex polygons instead of triangles. If triangles are
 * desired, the application may {@link #triangulate triangulate} the faces, but
 * it should be noted that this may result in triangles which are very small or
 * thin and hence difficult to perform reliable convexity tests on. In other
 * words, triangulating a merged face is likely to restore the numerical
 * problems which the merging process removed. Hence is it
 * possible that, after triangulation, {@link #check check} will fail (the same
 * behavior is observed with triangulated output from qhull ).
 *
 * Degenerate InputIt is assumed that the input points
 * are non-degenerate in that they are not coincident, colinear, or
 * colplanar, and thus the convex hull has a non-zero volume.
 * If the input points are detected to be degenerate within
 * the {@link #getDistanceTolerance() distance tolerance}, an
 * IllegalArgumentException will be thrown.
 *
 * @author John E. Lloyd, Fall 2004 */
class QuickHull3D
{
	/**
	 * Specifies that (on output) vertex indices for a face should be
	 * listed in clockwise order.
	 */
	public static final int CLOCKWISE = 0x1;

	/**
	 * Specifies that (on output) the vertex indices for a face should be
	 * numbered starting from 1.
	 */
	public static final int INDEXED_FROM_ONE = 0x2;

	/**
	 * Specifies that (on output) the vertex indices for a face should be
	 * numbered starting from 0.
	 */
	public static final int INDEXED_FROM_ZERO = 0x4;

	/**
	 * Specifies that (on output) the vertex indices for a face should be
	 * numbered with respect to the original input points.
	 */
	public static final int POINT_RELATIVE = 0x8;

	/**
	 * Specifies that the distance tolerance should be
	 * computed automatically from the input point data.
	 */
	public static final double AUTOMATIC_TOLERANCE = -1;

	/** The find index. */
	protected int findIndex = -1;

	/** The char length. */
	// estimated size of the point set
	protected double charLength;

	/** The debug. */
	protected boolean debug = false;

	/** The point buffer. */
	protected Vertex[] pointBuffer = new Vertex[0];
	
	/** The vertex point indices. */
	protected int[] vertexPointIndices = new int[0];
	
	/** The discarded faces. */
	private Face[] discardedFaces = new Face[3];

	/** The max vtxs. */
	private Vertex[] maxVtxs = new Vertex[3];
	
	/** The min vtxs. */
	private Vertex[] minVtxs = new Vertex[3];

	/** The faces. */
	protected Vector faces = new Vector(16);
	
	/** The horizon. */
	protected Vector horizon = new Vector(16);

	/** The new faces. */
	private FaceList newFaces = new FaceList();
	
	/** The unclaimed. */
	private VertexList unclaimed = new VertexList();
	
	/** The claimed. */
	private VertexList claimed = new VertexList();

	/** The num vertices. */
	protected int numVertices;
	
	/** The num faces. */
	protected int numFaces;
	
	/** The num points. */
	protected int numPoints;

	/** The explicit tolerance. */
	protected double explicitTolerance = AUTOMATIC_TOLERANCE;
	
	/** The tolerance. */
	protected double tolerance;

	/**
	 * Returns true if debugging is enabled.
	 *
	 * @return true is debugging is enabled
	 * @see QuickHull3D#setDebug
	 */
	public boolean getDebug()
	 {
	   return debug;
	 }

	/**
	 * Enables the printing of debugging diagnostics.
	 *
	 * @param enable if true, enables debugging
	 */
	public void setDebug (boolean enable)
	 { 
	   debug = enable;
	 }

	/**
	 * Precision of a double.
	 */
	static private final double DOUBLE_PREC = 2.2204460492503131e-16;


	/**
	 * Returns the distance tolerance that was used for the most recently
	 * computed hull. The distance tolerance is used to determine when
	 * faces are unambiguously convex with respect to each other, and when
	 * points are unambiguously above or below a face plane, in the
	 * presence of  =#distTol>numerical imprecision . Normally,
	 * this tolerance is computed automatically for each set of input
	 * points, but it can be set explicitly by the application.
	 *
	 * @return distance tolerance
	 * @see QuickHull3D#setExplicitDistanceTolerance
	 */
	public double getDistanceTolerance()
	 {
	   return tolerance;
	 }

	/**
	 * Sets an explicit distance tolerance for convexity tests.
	 * If {@link #AUTOMATIC_TOLERANCE AUTOMATIC_TOLERANCE}
	 * is specified (the default), then the tolerance will be computed
	 * automatically from the point data.
	 *
	 * @param tol explicit tolerance
	 * @see #getDistanceTolerance
	 */
	public void setExplicitDistanceTolerance(double tol)
	 { 
	   explicitTolerance = tol;
	 }

	/**
	 * Returns the explicit distance tolerance.
	 *
	 * @return explicit tolerance
	 * @see #setExplicitDistanceTolerance
	 */
	public double getExplicitDistanceTolerance()
	 {
	   return explicitTolerance;
	 }

	/**
	 * Adds the point to face.
	 *
	 * @param vtx the vtx
	 * @param face the face
	 */
	private void addPointToFace (Vertex vtx, Face face)
	 {
	   vtx.face = face;

	   if (face.outside == null)
	    { claimed.add (vtx);
	    }
	   else
	    { claimed.insertBefore (vtx, face.outside); 
	    }
	   face.outside = vtx;
	 }

	/**
	 * Removes the point from face.
	 *
	 * @param vtx the vtx
	 * @param face the face
	 */
	private void removePointFromFace (Vertex vtx, Face face)
	 {
	   if (vtx == face.outside)
	    { if (vtx.next != null && vtx.next.face == face)
	       { face.outside = vtx.next;
	       }
	      else
	       { face.outside = null; 
	       }
	    }
	   claimed.delete (vtx);
	 }

	/**
	 * Removes the all points from face.
	 *
	 * @param face the face
	 * @return the vertex
	 */
	private Vertex removeAllPointsFromFace (Face face)
	 {
	   if (face.outside != null)
	    { 
	      Vertex end = face.outside;
	      while (end.next != null && end.next.face == face)
	       { end = end.next;
	       }
	      claimed.delete (face.outside, end);
	      end.next = null;
	      return face.outside;
	    }
	   else
	    { return null; 
	    }
	 }

	/**
	 * Creates an empty convex hull object.
	 */
	public QuickHull3D ()
	 { 
	 }

	/**
	 * Creates a convex hull object and initializes it to the convex hull
	 * of a set of points whose coordinates are given by an
	 * array of doubles.
	 *
	 * @param coords x, y, and z coordinates of each input
	 * point. The length of this array will be three times
	 * the the number of input points.
	 * @throws IllegalArgumentException the number of input points is less
	 * than four, or the points appear to be coincident, colinear, or
	 * coplanar.
	 */
	public QuickHull3D (double[] coords)
	   throws IllegalArgumentException
	 {
	   build (coords, coords.length/3);
	 }

	/**
	 * Creates a convex hull object and initializes it to the convex hull
	 * of a set of points.
	 *
	 * @param points input points.
	 * @throws IllegalArgumentException the number of input points is less
	 * than four, or the points appear to be coincident, colinear, or
	 * coplanar.
	 */
	public QuickHull3D (Point3d[] points)
	   throws IllegalArgumentException
	 {
	   build (points, points.length);
	 }

	/**
	 * Find half edge.
	 *
	 * @param tail the tail
	 * @param head the head
	 * @return the half edge
	 */
	private HalfEdge findHalfEdge (Vertex tail, Vertex head)
	 { 
	   // brute force ... OK, since setHull is not used much
	   for (Iterator it=faces.iterator(); it.hasNext(); ) 
	    { HalfEdge he = ((Face)it.next()).findEdge (tail, head);
	      if (he != null)
	       { return he; 
	       }
	    }
	   return null;
	 }

 	/**
	  * Sets the hull.
	  *
	  * @param coords the coords
	  * @param nump the nump
	  * @param faceIndices the face indices
	  * @param numf the numf
	  */
	 protected void setHull (double[] coords, int nump,
				int[][] faceIndices, int numf)
 	 {
 	   initBuffers (nump);
	   setPoints (coords, nump);
	   computeMaxAndMin ();
	   for (int i=0; i 0)
	    { System.out.write (es.read());
	      wrote = true;
	    }
	   if (wrote)
	    { System.out.println("");
	    }
	 }

	/**
	 * Sets the from qhull.
	 *
	 * @param coords the coords
	 * @param nump the nump
	 * @param triangulate the triangulate
	 */
	protected void setFromQhull (double[] coords, int nump,
				     boolean triangulate)
	 {
	   String commandStr = "./qhull i";
	   if (triangulate)
	    { commandStr += " -Qt"; 
	    }
	   try
	    { 
	      Process proc = Runtime.getRuntime().exec (commandStr);
	      PrintStream ps = new PrintStream (proc.getOutputStream());
	      StreamTokenizer stok =
		 new StreamTokenizer (
		    new InputStreamReader (proc.getInputStream()));

	      ps.println ("3 " + nump);
	      for (int i=0; inump.
	 * @param nump number of input points
	 * @throws IllegalArgumentException the number of input points is less
	 * than four or greater than 1/3 the length of coords,
	 * or the points appear to be coincident, colinear, or
	 * coplanar.
	 */
	public void build (double[] coords, int nump)
	   throws IllegalArgumentException
	 {
	   if (nump < 4)
	    { throw new IllegalArgumentException (
		 "Less than four input points specified");
	    }
	   if (coords.length/3 < nump)
	    { throw new IllegalArgumentException (
		 "Coordinate array too small for specified number of points"); 
	    }
	   initBuffers (nump);
	   setPoints (coords, nump);
	   buildHull ();
	 }

	/**
	 * Constructs the convex hull of a set of points.
	 *
	 * @param points input points
	 * @throws IllegalArgumentException the number of input points is less
	 * than four, or the points appear to be coincident, colinear, or
	 * coplanar.
	 */
	public void build (Point3d[] points)
	   throws IllegalArgumentException
	 {
	   build (points, points.length);
	 }

	/**
	 * Constructs the convex hull of a set of points.
	 *
	 * @param points input points
	 * @param nump number of input points
	 * @throws IllegalArgumentException the number of input points is less
	 * than four or greater then the length of points, or the
	 * points appear to be coincident, colinear, or coplanar.
	 */
	public void build(Point3d[] points, int nump) throws IllegalArgumentException {
		if (nump < 4) {
			throw new IllegalArgumentException("Less than four input points specified");
		}
		if (points.length < nump) {
			throw new IllegalArgumentException("Point array too small for specified number of points");
		}
		HashSet set = new HashSet();
		for(int i=0;iqhull ).
	 */
	public void triangulate ()
	 {
	   double minArea = 1000*charLength*DOUBLE_PREC;
	   newFaces.clear();
	   for (Iterator it=faces.iterator(); it.hasNext(); ) 
	    { Face face = (Face)it.next();
	      if (face.mark == Face.VISIBLE)
	       { 
		 face.triangulate (newFaces, minArea);
		 // splitFace (face);
	       }
	    }
	   for (Face face=newFaces.first(); face!=null; face=face.next)
	    { faces.add (face);
	    }
	 }

// 	private void splitFace (Face face)
// 	 {
//  	   Face newFace = face.split();
//  	   if (newFace != null)
//  	    { newFaces.add (newFace);
//  	      splitFace (newFace);
//  	      splitFace (face);
//  	    }
// 	 }

	/**
 * Inits the buffers.
 *
 * @param nump the nump
 */
protected void initBuffers (int nump)
	 {
	   if (pointBuffer.length < nump)
	    { Vertex[] newBuffer = new Vertex[nump];
	      vertexPointIndices = new int[nump];
	      for (int i=0; i max.x)
	       { max.x = pnt.x;
		 maxVtxs[0] = pointBuffer[i];
	       }
	      else if (pnt.x < min.x)
	       { min.x = pnt.x;
		 minVtxs[0] = pointBuffer[i];
	       }
	      if (pnt.y > max.y)
	       { max.y = pnt.y;
		 maxVtxs[1] = pointBuffer[i];
	       }
	      else if (pnt.y < min.y)
	       { min.y = pnt.y;
		 minVtxs[1] = pointBuffer[i];
	       }
	      if (pnt.z > max.z)
	       { max.z = pnt.z;
		 maxVtxs[2] = pointBuffer[i];
	       }
	      else if (pnt.z < min.z)
	       { min.z = pnt.z;
		 minVtxs[2] = pointBuffer[i];
	       }
	    }

	   // this epsilon formula comes from QuickHull, and I'm
	   // not about to quibble.
	   charLength = Math.max(max.x-min.x, max.y-min.y);
	   charLength = Math.max(max.z-min.z, charLength);
	   if (explicitTolerance == AUTOMATIC_TOLERANCE)
	    { tolerance =
		 3*DOUBLE_PREC*(Math.max(Math.abs(max.x),Math.abs(min.x))+
				Math.max(Math.abs(max.y),Math.abs(min.y))+
				Math.max(Math.abs(max.z),Math.abs(min.z)));
	    }
	   else
	    { tolerance = explicitTolerance; 
	    }
	 }

	/**
	 * Creates the initial simplex from which the hull will be built.
	 *
	 * @throws IllegalArgumentException the illegal argument exception
	 */
	protected void createInitialSimplex ()
	   throws IllegalArgumentException
	 {
	   double max = 0;
	   int imax = 0;

	   for (int i=0; i<3; i++)
	    { double diff = maxVtxs[i].pnt.get(i)-minVtxs[i].pnt.get(i);
	      if (diff > max)
	       { max = diff;
		 imax = i;
	       }
 	    }

	   if (max <= tolerance)
	    { throw new IllegalArgumentException (
"Input points appear to be coincident");
	    }
	   Vertex[] vtx = new Vertex[4];
	   // set first two vertices to be those with the greatest
	   // one dimensional separation

	   vtx[0] = maxVtxs[imax];
	   vtx[1] = minVtxs[imax];

	   // set third vertex to be the vertex farthest from
	   // the line between vtx0 and vtx1
	   Vector3d u01 = new Vector3d();
	   Vector3d diff02 = new Vector3d();
	   Vector3d nrml = new Vector3d();
	   Vector3d xprod = new Vector3d();
	   double maxSqr = 0;
	   u01.sub (vtx[1].pnt, vtx[0].pnt);
	   u01.normalize();
	   for (int i=0; i maxSqr &&
		  pointBuffer[i] != vtx[0] &&  // paranoid
		  pointBuffer[i] != vtx[1])
	       { maxSqr = lenSqr; 
		 vtx[2] = pointBuffer[i];
		 nrml.set (xprod);
	       }
	    }
	   if (Math.sqrt(maxSqr) <= 100*tolerance)
	    { throw new IllegalArgumentException (
"Input points appear to be colinear");
	    }
	   nrml.normalize();


	   double maxDist = 0;
	   double d0 = vtx[2].pnt.dot (nrml);
	   for (int i=0; i maxDist &&
		  pointBuffer[i] != vtx[0] &&  // paranoid
		  pointBuffer[i] != vtx[1] &&
		  pointBuffer[i] != vtx[2])
	       { maxDist = dist;
		 vtx[3] = pointBuffer[i];
	       }
	    }
	   if (Math.abs(maxDist) <= 100*tolerance)
	    { throw new IllegalArgumentException (
"Input points appear to be coplanar"); 
	    }

	   if (debug)
	    { System.out.println ("initial vertices:");
	      System.out.println (vtx[0].index + ": " + vtx[0].pnt);
	      System.out.println (vtx[1].index + ": " + vtx[1].pnt);
	      System.out.println (vtx[2].index + ": " + vtx[2].pnt);
	      System.out.println (vtx[3].index + ": " + vtx[3].pnt);
	    }

	   Face[] tris = new Face[4];

	   if (vtx[3].pnt.dot (nrml) - d0 < 0)
	    { tris[0] = Face.createTriangle (vtx[0], vtx[1], vtx[2]);
	      tris[1] = Face.createTriangle (vtx[3], vtx[1], vtx[0]);
	      tris[2] = Face.createTriangle (vtx[3], vtx[2], vtx[1]);
	      tris[3] = Face.createTriangle (vtx[3], vtx[0], vtx[2]);

	      for (int i=0; i<3; i++)
	       { int k = (i+1)%3;
		 tris[i+1].getEdge(1).setOpposite (tris[k+1].getEdge(0));
		 tris[i+1].getEdge(2).setOpposite (tris[0].getEdge(k));
	       }
	    }
	   else
	    { tris[0] = Face.createTriangle (vtx[0], vtx[2], vtx[1]);
	      tris[1] = Face.createTriangle (vtx[3], vtx[0], vtx[1]);
	      tris[2] = Face.createTriangle (vtx[3], vtx[1], vtx[2]);
	      tris[3] = Face.createTriangle (vtx[3], vtx[2], vtx[0]);

	      for (int i=0; i<3; i++)
	       { int k = (i+1)%3;
		 tris[i+1].getEdge(0).setOpposite (tris[k+1].getEdge(1));
		 tris[i+1].getEdge(2).setOpposite (tris[0].getEdge((3-i)%3));
	       }
	    }


 	   for (int i=0; i<4; i++)
 	    { faces.add (tris[i]); 
 	    }

	   for (int i=0; i maxDist)
		  { maxFace = tris[k];
		    maxDist = dist;
		  }
	       }
	      if (maxFace != null)
	       { addPointToFace (v, maxFace);
	       }	      
	    }
	 }

	/**
	 * Returns the number of vertices in this hull.
	 *
	 * @return number of vertices
	 */
	public int getNumVertices()
	 {
	   return numVertices;
	 }

	/**
	 * Returns the vertex points in this hull.
	 *
	 * @return array of vertex points
	 * @see QuickHull3D#getVertices(double[])
	 * @see QuickHull3D#getFaces()
	 */
 	public Point3d[] getVertices()
 	 {
 	   Point3d[] vtxs = new Point3d[numVertices];
 	   for (int i=0; iv), followed by the vertex indices
	 * for each face (each
	 * preceded by the letter f).
	 *
	 *  The face indices are numbered with respect to the hull vertices
	 * (as opposed to the input points), with a lowest index of 1, and are
	 * arranged counter-clockwise. More control over the index format can
	 * be obtained using
	 * {@link #print(PrintStream,int) print(ps,indexFlags)}.
	 *
	 * @param ps stream used for printing
	 * @see QuickHull3D#print(PrintStream,int)
	 * @see QuickHull3D#getVertices()
	 * @see QuickHull3D#getFaces()
	 */
	public void print (PrintStream ps)
	 {
	   print (ps, 0);
	 }

	/**
	 * Prints the vertices and faces of this hull to the stream ps.
	 *
	 *   This is done using the Alias Wavefront .obj file format, with
	 * the vertices printed first (each preceding by the letter
	 * v), followed by the vertex indices for each face (each
	 * preceded by the letter f).
	 *
	 *  By default, the face indices are numbered with respect to the
	 * hull vertices (as opposed to the input points), with a lowest index
	 * of 1, and are arranged counter-clockwise. However, this
	 * can be changed by setting {@link #POINT_RELATIVE POINT_RELATIVE},
	 * {@link #INDEXED_FROM_ONE INDEXED_FROM_ZERO}, or {@link #CLOCKWISE
	 * CLOCKWISE} in the indexFlags parameter.
	 *
	 * @param ps stream used for printing
	 * @param indexFlags specifies index characteristics
	 * (0 results in the default).
	 * @see QuickHull3D#getVertices()
	 * @see QuickHull3D#getFaces()
	 */
	public void print (PrintStream ps, int indexFlags)
	 {
	   if ((indexFlags & INDEXED_FROM_ZERO) == 0)
	    { indexFlags |= INDEXED_FROM_ONE;
	    }
	   for (int i=0; i maxDist)
		     { maxDist = dist;
		       maxFace = newFace;
		     }
		    if (maxDist > 1000*tolerance)
		     { break;
		     }
		  }
	       }
	      if (maxFace != null)
	       { 
		 addPointToFace (vtx, maxFace);
 		 if (debug && vtx.index == findIndex)
 		  { System.out.println (findIndex + " CLAIMED BY " +
 		     maxFace.getVertexString()); 
 		  }
	       }
	      else
	       { if (debug && vtx.index == findIndex)
		  { System.out.println (findIndex + " DISCARDED"); 
		  } 
	       }
	    }
	 }

	/**
	 * Delete face points.
	 *
	 * @param face the face
	 * @param absorbingFace the absorbing face
	 */
	protected void deleteFacePoints (Face face, Face absorbingFace)
	 {
	   Vertex faceVtxs = removeAllPointsFromFace (face);
	   if (faceVtxs != null)
	    { 
	      if (absorbingFace == null)
	       { unclaimed.addAll (faceVtxs);
	       }
	      else
	       { Vertex vtxNext = faceVtxs;
		 for (Vertex vtx=vtxNext; vtx!=null; vtx=vtxNext)
		  { vtxNext = vtx.next;
		    double dist = absorbingFace.distanceToPlane (vtx.pnt);
		    if (dist > tolerance)
		     { 
		       addPointToFace (vtx, absorbingFace);
		     }
		    else
		     { 
		       unclaimed.add (vtx);
		     }
		  }
	       }
	    }
	 }

	/** The Constant NONCONVEX_WRT_LARGER_FACE. */
	private static final int NONCONVEX_WRT_LARGER_FACE = 1;
	
	/** The Constant NONCONVEX. */
	private static final int NONCONVEX = 2;

	/**
	 * Opp face distance.
	 *
	 * @param he the he
	 * @return the double
	 */
	protected double oppFaceDistance (HalfEdge he)
	 {
	   return he.face.distanceToPlane (he.opposite.face.getCentroid());
	 }

	/**
	 * Do adjacent merge.
	 *
	 * @param face the face
	 * @param mergeType the merge type
	 * @return true, if successful
	 */
	private boolean doAdjacentMerge (Face face, int mergeType)
	 {
	   HalfEdge hedge = face.he0;

	   boolean convex = true;
	   do
	    { Face oppFace = hedge.oppositeFace();
	      boolean merge = false;
	      double dist1, dist2;

	      if (mergeType == NONCONVEX)
	       { // then merge faces if they are definitively non-convex
		 if (oppFaceDistance (hedge) > -tolerance ||
		     oppFaceDistance (hedge.opposite) > -tolerance)
		  { merge = true;
		  }
	       }
	      else // mergeType == NONCONVEX_WRT_LARGER_FACE
	       { // merge faces if they are parallel or non-convex
		 // wrt to the larger face; otherwise, just mark
		 // the face non-convex for the second pass.
		 if (face.area > oppFace.area)
		  { if ((dist1 = oppFaceDistance (hedge)) > -tolerance) 
		     { merge = true;
		     }
		    else if (oppFaceDistance (hedge.opposite) > -tolerance)
		     { convex = false;
		     }
		  }
		 else
		  { if (oppFaceDistance (hedge.opposite) > -tolerance)
		     { merge = true;
		     }
		    else if (oppFaceDistance (hedge) > -tolerance) 
		     { convex = false;
		     }
		  }
	       }

	      if (merge)
	       { if (debug)
		  { System.out.println (
		    "  merging " + face.getVertexString() + "  and  " +
		    oppFace.getVertexString());
		  }

		 int numd = face.mergeAdjacentFace (hedge, discardedFaces);
		 for (int i=0; i tolerance)
		  { calculateHorizon (eyePnt, edge.getOpposite(),
				      oppFace, horizon);
		  }
		 else
		  { horizon.add (edge);
		    if (debug)
		     { System.out.println ("  adding horizon edge " +
					   edge.getVertexString());
		     }
		  }
	       }
	      edge = edge.getNext();
	    }
	   while (edge != edge0);
	 }

	/**
	 * Adds the adjoining face.
	 *
	 * @param eyeVtx the eye vtx
	 * @param he the he
	 * @return the half edge
	 */
	private HalfEdge addAdjoiningFace (
	   Vertex eyeVtx, HalfEdge he)
	 { 
	   Face face = Face.createTriangle (
	      eyeVtx, he.tail(), he.head());
	   faces.add (face);
	   face.getEdge(-1).setOpposite(he.getOpposite());
	   return face.getEdge(0);
	 }

	/**
	 * Adds the new faces.
	 *
	 * @param newFaces the new faces
	 * @param eyeVtx the eye vtx
	 * @param horizon the horizon
	 */
	protected void addNewFaces (
	   FaceList newFaces, Vertex eyeVtx, Vector horizon)
	 { 
	   newFaces.clear();

	   HalfEdge hedgeSidePrev = null;
	   HalfEdge hedgeSideBegin = null;

		for (Iterator it = horizon.iterator(); it.hasNext();) {
			try {
				HalfEdge horizonHe = (HalfEdge) it.next();
				HalfEdge hedgeSide = addAdjoiningFace(eyeVtx, horizonHe);
				if (debug) {
					System.out.println("new face: " + hedgeSide.face.getVertexString());
				}
				if (hedgeSidePrev != null) {
					hedgeSide.next.setOpposite(hedgeSidePrev);
				} else {
					hedgeSideBegin = hedgeSide;
				}
				newFaces.add(hedgeSide.getFace());
				hedgeSidePrev = hedgeSide;
			}catch(Throwable t) {
				t.printStackTrace();
			}
		}
		hedgeSideBegin.next.setOpposite(hedgeSidePrev);
	}

	/**
	 * Next point to add.
	 *
	 * @return the vertex
	 */
	protected Vertex nextPointToAdd()
	 {
	   if (!claimed.isEmpty())
	    { Face eyeFace = claimed.first().face;
	      Vertex eyeVtx = null;
	      double maxDist = 0;
	      for (Vertex vtx=eyeFace.outside;
		   vtx != null && vtx.face==eyeFace;
		   vtx = vtx.next)
	       { double dist = eyeFace.distanceToPlane(vtx.pnt);
		 if (dist > maxDist)
		  { maxDist = dist;
		    eyeVtx = vtx;
		  }
	       }
	      return eyeVtx;
	    }
	   else
	    { return null;
	    }
	 }
	
	/**
	 * Adds the point to hull.
	 *
	 * @param eyeVtx the eye vtx
	 */
	protected void addPointToHull(Vertex eyeVtx)
	 {
	     horizon.clear();
	     unclaimed.clear();
	      
	     if (debug)
	      { System.out.println ("Adding point: " + eyeVtx.index);
		System.out.println (
		   " which is " + eyeVtx.face.distanceToPlane(eyeVtx.pnt) +
		   " above face " + eyeVtx.face.getVertexString());
	      }
	     removePointFromFace (eyeVtx, eyeVtx.face);
	     calculateHorizon (eyeVtx.pnt, null, eyeVtx.face, horizon);
	     newFaces.clear();
	     addNewFaces (newFaces, eyeVtx, horizon);
	     
	     // first merge pass ... merge faces which are non-convex
	     // as determined by the larger face
	     
	     for (Face face = newFaces.first(); face!=null; face=face.next)
	      { 
		if (face.mark == Face.VISIBLE)
		 { while (doAdjacentMerge(face, NONCONVEX_WRT_LARGER_FACE))
		      ;
		 }
	      }		 
	     // second merge pass ... merge faces which are non-convex
	     // wrt either face	     
	     for (Face face = newFaces.first(); face!=null; face=face.next)
	      { 
 		if (face.mark == Face.NON_CONVEX)
		 { face.mark = Face.VISIBLE;
		   while (doAdjacentMerge(face, NONCONVEX))
		      ;
 		 }
 	      }	
	     resolveUnclaimedPoints(newFaces);
	 }

	/**
	 * Builds the hull.
	 */
	protected void buildHull ()
	 {
	   int cnt = 0;
	   Vertex eyeVtx;

	   computeMaxAndMin ();
	   createInitialSimplex ();
	   while ((eyeVtx = nextPointToAdd()) != null)
	    { addPointToHull (eyeVtx);
	      cnt++;
	      if (debug)
	       { System.out.println ("iteration " + cnt + " done"); 
	       }
	    }
	   reindexFacesAndVertices();
	   if (debug)
	    { System.out.println ("hull done");
	    }
	 }

	/**
	 * Mark face vertices.
	 *
	 * @param face the face
	 * @param mark the mark
	 */
	private void markFaceVertices (Face face, int mark)
	 {
	   HalfEdge he0 = face.getFirstEdge();
	   HalfEdge he = he0;
	   do
	    { he.head().index = mark;
	      he = he.next;
	    }
	   while (he != he0);
	 }

	/**
	 * Reindex faces and vertices.
	 */
	protected void reindexFacesAndVertices()
	 { 
	   for (int i=0; i tol)
	       { if (ps != null)
		  { ps.println ("Edge " + he.getVertexString() +
				" non-convex by " + dist);
		  }
		 return false;
	       }
	      dist = oppFaceDistance (he.opposite);
	      if (dist > tol)
	       { if (ps != null)
		  { ps.println ("Opposite edge " +
				he.opposite.getVertexString() +
				" non-convex by " + dist);
		  }
		 return false;
	       }
	      if (he.next.oppositeFace() == he.oppositeFace())
	       { if (ps != null)
		  { ps.println ("Redundant vertex " + he.head().index +
				" in face " + face.getVertexString());
		  }
		 return false;
	       }
	      he = he.next;
	    }
	   while (he != face.he0);	   
	   return true;
	 }

	/**
	 * Check faces.
	 *
	 * @param tol the tol
	 * @param ps the ps
	 * @return true, if successful
	 */
	protected boolean checkFaces(double tol, PrintStream ps)
	 { 
	   // check edge convexity
	   boolean convex = true;
	   for (Iterator it=faces.iterator(); it.hasNext(); ) 
	    { Face face = (Face)it.next();
	      if (face.mark == Face.VISIBLE)
	       { if (!checkFaceConvexity (face, tol, ps))
		  { convex = false;
		  }
	       }
	    }
	   return convex;
	 }

	/**
	 * Checks the correctness of the hull using the distance tolerance
	 * returned by {@link QuickHull3D#getDistanceTolerance
	 * getDistanceTolerance}; see
	 * {@link QuickHull3D#check(PrintStream,double)
	 * check(PrintStream,double)} for details.
	 *
	 * @param ps print stream for diagnostic messages; may be
	 * set to null if no messages are desired.
	 * @return true if the hull is valid
	 * @see QuickHull3D#check(PrintStream,double)
	 */
	public boolean check (PrintStream ps)
	 {
	   return check (ps, getDistanceTolerance());
	 }

	/**
	 * Checks the correctness of the hull. This is done by making sure that
	 * no faces are non-convex and that no points are outside any face.
	 * These tests are performed using the distance tolerance  tol .
	 * Faces are considered non-convex if any edge is non-convex, and an
	 * edge is non-convex if the centroid of either adjoining face is more
	 * than  tol  above the plane of the other face. Similarly,
	 * a point is considered outside a face if its distance to that face's
	 * plane is more than 10 times  tol .
	 *
	 *  If the hull has been {@link #triangulate triangulated},
	 * then this routine may fail if some of the resulting
	 * triangles are very small or thin.
	 *
	 * @param ps print stream for diagnostic messages; may be
	 * set to null if no messages are desired.
	 * @param tol distance tolerance
	 * @return true if the hull is valid
	 * @see QuickHull3D#check(PrintStream)
	 */
	public boolean check (PrintStream ps, double tol)

	 {
	   // check to make sure all edges are fully connected
	   // and that the edges are convex
	   double dist;
	   double pointTol = 10*tol;

	   if (!checkFaces(tolerance, ps))
	    { return false; 
	    }

	   // check point inclusion

	   for (int i=0; i pointTol)
		     { if (ps != null)
			{ ps.println (
			     "Point " + i + " " + dist + " above face " +
			     face.getVertexString());
			}
		       return false;
		     }
		  }
	       }
	    }
	   return true;
	 }
}
Related Artifacts