/**
* 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 quickhull3d;
import java.util.*;
import java.io.*;
/**
* Computes the convex hull of a set of three dimensional points.
*
* <p>The algorithm is a three dimensional implementation of Quickhull, as
* described in Barber, Dobkin, and Huhdanpaa, <a
* href=http://citeseer.ist.psu.edu/barber96quickhull.html> ``The Quickhull
* Algorithm for Convex Hulls''</a> (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 <a
* href=http://www.qhull.org>qhull</a>.
*
* <p>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:
* <pre>
* // 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 < vertices.length; i++)
* { for (int k = 0; k < faceIndices[i].length; k++)
* { System.out.print (faceIndices[i][k] + " ");
* }
* System.out.println ("");
* }
* </pre>
* 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.
*
* <h3><a name=distTol>Robustness</h3> 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 <a
* href=http://www.qhull.org>qhull</a>, 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 <i>below</i> 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}.
*
* <p>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
* <code>true</code>. These tests have been successful for a large number of
* trials and so we are confident that QuickHull3D is reasonably robust.
*
* <h3>Merged Faces</h3> 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 <a
* href=http://www.qhull.org>qhull</a>).
*
* <h3>Degenerate Input</h3>It 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 */
public 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;
protected int findIndex = -1;
// estimated size of the point set
protected double charLength;
protected boolean debug = false;
protected Vertex[] pointBuffer = new Vertex[0];
protected int[] vertexPointIndices = new int[0];
private Face[] discardedFaces = new Face[3];
private Vertex[] maxVtxs = new Vertex[3];
private Vertex[] minVtxs = new Vertex[3];
protected List<Face> faces = new ArrayList<Face>(16);
protected List<HalfEdge> horizon = new ArrayList<HalfEdge>(16);
private FaceList newFaces = new FaceList();
private VertexList unclaimed = new VertexList();
private VertexList claimed = new VertexList();
protected int numVertices;
protected int numFaces;
protected int numPoints;
protected double explicitTolerance = AUTOMATIC_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 <a href=#distTol>numerical imprecision</a>. 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;
}
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;
}
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);
}
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);
}
private HalfEdge findHalfEdge (Vertex tail, Vertex head)
{
// brute force ... OK, since setHull is not used much
for (Face face:faces)
{ HalfEdge he = face.findEdge (tail, head);
if (he != null)
{ return he;
}
}
return null;
}
protected void setHull (double[] coords, int nump,
int[][] faceIndices, int numf)
{
initBuffers (nump);
setPoints (coords, nump);
computeMaxAndMin ();
for (int i=0; i<numf; i++)
{ Face face = Face.create (pointBuffer, faceIndices[i]);
HalfEdge he = face.he0;
do
{ HalfEdge heOpp = findHalfEdge (he.head(), he.tail());
if (heOpp != null)
{ he.setOpposite (heOpp);
}
he = he.next;
}
while (he != face.he0);
faces.add (face);
}
}
private void printQhullErrors (Process proc)
throws IOException
{
boolean wrote = false;
InputStream es = proc.getErrorStream();
while (es.available() > 0)
{ System.out.write (es.read());
wrote = true;
}
if (wrote)
{ System.out.println("");
}
}
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; i<nump; i++)
{ ps.println (
coords[i*3+0] + " " +
coords[i*3+1] + " " +
coords[i*3+2]);
}
ps.flush();
ps.close();
List<Integer> indexList = new ArrayList<Integer>(3);
stok.eolIsSignificant(true);
printQhullErrors (proc);
do
{ stok.nextToken();
}
while (stok.sval == null ||
!stok.sval.startsWith ("MERGEexact"));
for (int i=0; i<4; i++)
{ stok.nextToken();
}
if (stok.ttype != StreamTokenizer.TT_NUMBER)
{ System.out.println ("Expecting number of faces");
System.exit(1);
}
int numf = (int)stok.nval;
stok.nextToken(); // clear EOL
int[][] faceIndices = new int[numf][];
for (int i=0; i<numf; i++)
{ indexList.clear();
while (stok.nextToken() != StreamTokenizer.TT_EOL)
{ if (stok.ttype != StreamTokenizer.TT_NUMBER)
{ System.out.println ("Expecting face index");
System.exit(1);
}
indexList.add (0, new Integer((int)stok.nval));
}
faceIndices[i] = new int[indexList.size()];
int k = 0;
for (Integer it: indexList )
{ faceIndices[i][k++] = it.intValue();
}
}
setHull (coords, nump, faceIndices, numf);
}
catch (Exception e)
{ e.printStackTrace();
System.exit(1);
}
}
private void printPoints (PrintStream ps)
{
for (int i=0; i<numPoints; i++)
{ Point3d pnt = pointBuffer[i].pnt;
ps.println (pnt.x + ", " + pnt.y + ", " + pnt.z + ",");
}
}
/**
* Constructs 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 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 void build (double[] coords)
throws IllegalArgumentException
{
build (coords, coords.length/3);
}
/**
* Constructs 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 must be at least three times
* <code>nump</code>.
* @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 <code>coords</code>,
* 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 <code>points</code>, 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");
}
initBuffers (nump);
setPoints (points, nump);
buildHull ();
}
/**
* Triangulates any non-triangular hull faces. In some cases, due to
* precision issues, the resulting triangles may be very thin or small,
* and hence appear to be non-convex (this same limitation is present
* in <a href=http://www.qhull.org>qhull</a>).
*/
public void triangulate ()
{
double minArea = 1000*charLength*DOUBLE_PREC;
newFaces.clear();
for (Face face: faces)
{
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);
// }
// }
protected void initBuffers (int nump)
{
if (pointBuffer.length < nump)
{ Vertex[] newBuffer = new Vertex[nump];
vertexPointIndices = new int[nump];
for (int i=0; i<pointBuffer.length; i++)
{ newBuffer[i] = pointBuffer[i];
}
for (int i=pointBuffer.length; i<nump; i++)
{ newBuffer[i] = new Vertex();
}
pointBuffer = newBuffer;
}
faces.clear();
claimed.clear();
numFaces = 0;
numPoints = nump;
}
protected void setPoints (double[] coords, int nump)
{
for (int i=0; i<nump; i++)
{
Vertex vtx = pointBuffer[i];
vtx.pnt.set (coords[i*3+0], coords[i*3+1], coords[i*3+2]);
vtx.index = i;
}
}
protected void setPoints (Point3d[] pnts, int nump)
{
for (int i=0; i<nump; i++)
{
Vertex vtx = pointBuffer[i];
vtx.pnt.set (pnts[i]);
vtx.index = i;
}
}
protected void computeMaxAndMin ()
{
Vector3d max = new Vector3d();
Vector3d min = new Vector3d();
for (int i=0; i<3; i++)
{ maxVtxs[i] = minVtxs[i] = pointBuffer[0];
}
max.set (pointBuffer[0].pnt);
min.set (pointBuffer[0].pnt);
for (int i=1; i<numPoints; i++)
{ Point3d pnt = pointBuffer[i].pnt;
if (pnt.x > 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;
maxVtxs[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.
*/
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<numPoints; i++)
{ diff02.sub (pointBuffer[i].pnt, vtx[0].pnt);
xprod.cross (u01, diff02);
double lenSqr = xprod.normSquared();
if (lenSqr > 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<numPoints; i++)
{ double dist = Math.abs (pointBuffer[i].pnt.dot(nrml) - d0);
if (dist > 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<numPoints; i++)
{ Vertex v = pointBuffer[i];
if (v == vtx[0] || v == vtx[1] || v == vtx[2] || v == vtx[3])
{ continue;
}
maxDist = tolerance;
Face maxFace = null;
for (int k=0; k<4; k++)
{ double dist = tris[k].distanceToPlane (v.pnt);
if (dist > 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; i<numVertices; i++)
{ vtxs[i] = pointBuffer[vertexPointIndices[i]].pnt;
}
return vtxs;
}
/**
* Returns the coordinates of the vertex points of this hull.
*
* @param coords returns the x, y, z coordinates of each vertex.
* This length of this array must be at least three times
* the number of vertices.
* @return the number of vertices
* @see QuickHull3D#getVertices()
* @see QuickHull3D#getFaces()
*/
public int getVertices(double[] coords)
{
for (int i=0; i<numVertices; i++)
{ Point3d pnt = pointBuffer[vertexPointIndices[i]].pnt;
coords[i*3+0] = pnt.x;
coords[i*3+1] = pnt.y;
coords[i*3+2] = pnt.z;
}
return numVertices;
}
/**
* Returns an array specifing the index of each hull vertex
* with respect to the original input points.
*
* @return vertex indices with respect to the original points
*/
public int[] getVertexPointIndices()
{
int[] indices = new int[numVertices];
for (int i=0; i<numVertices; i++)
{ indices[i] = vertexPointIndices[i];
}
return indices;
}
/**
* Returns the number of faces in this hull.
*
* @return number of faces
*/
public int getNumFaces()
{
return faces.size();
}
/**
* Returns the faces associated with this hull.
*
* <p>Each face is represented by an integer array which gives the
* indices of the vertices. These indices are numbered
* relative to the
* hull vertices, are zero-based,
* and are arranged counter-clockwise. More control
* over the index format can be obtained using
* {@link #getFaces(int) getFaces(indexFlags)}.
*
* @return array of integer arrays, giving the vertex
* indices for each face.
* @see QuickHull3D#getVertices()
* @see QuickHull3D#getFaces(int)
*/
public int[][] getFaces ()
{
return getFaces(0);
}
/**
* Returns the faces associated with this hull.
*
* <p>Each face is represented by an integer array which gives the
* indices of the vertices. By default, these indices are numbered with
* respect to the hull vertices (as opposed to the input points), are
* zero-based, and are arranged counter-clockwise. However, this
* can be changed by setting {@link #POINT_RELATIVE
* POINT_RELATIVE}, {@link #INDEXED_FROM_ONE INDEXED_FROM_ONE}, or
* {@link #CLOCKWISE CLOCKWISE} in the indexFlags parameter.
*
* @param indexFlags specifies index characteristics (0 results
* in the default)
* @return array of integer arrays, giving the vertex
* indices for each face.
* @see QuickHull3D#getVertices()
*/
public int[][] getFaces (int indexFlags)
{
int[][] allFaces = new int[faces.size()][];
int k = 0;
for (Face face: faces )
{
allFaces[k] = new int[face.numVertices()];
getFaceIndices (allFaces[k], face, indexFlags);
k++;
}
return allFaces;
}
/**
* Prints the vertices and faces of this hull to the stream ps.
*
* <p>
* This is done using the Alias Wavefront .obj file
* format, with the vertices printed first (each preceding by
* the letter <code>v</code>), followed by the vertex indices
* for each face (each
* preceded by the letter <code>f</code>).
*
* <p>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.
*
* <p> This is done using the Alias Wavefront .obj file format, with
* the vertices printed first (each preceding by the letter
* <code>v</code>), followed by the vertex indices for each face (each
* preceded by the letter <code>f</code>).
*
* <p>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<numVertices; i++)
{ Point3d pnt = pointBuffer[vertexPointIndices[i]].pnt;
ps.println ("v " + pnt.x + " " + pnt.y + " " + pnt.z);
}
for (Face face: faces )
{
int[] indices = new int[face.numVertices()];
getFaceIndices (indices, face, indexFlags);
ps.print ("f");
for (int k=0; k<indices.length; k++)
{ ps.print (" " + indices[k]);
}
ps.println ("");
}
}
private void getFaceIndices (int[] indices, Face face, int flags)
{
boolean ccw = ((flags & CLOCKWISE) == 0);
boolean indexedFromOne = ((flags & INDEXED_FROM_ONE) != 0);
boolean pointRelative = ((flags & POINT_RELATIVE) != 0);
HalfEdge hedge = face.he0;
int k = 0;
do
{ int idx = hedge.head().index;
if (pointRelative)
{ idx = vertexPointIndices[idx];
}
if (indexedFromOne)
{ idx++;
}
indices[k++] = idx;
hedge = (ccw ? hedge.next : hedge.prev);
}
while (hedge != face.he0);
}
protected void resolveUnclaimedPoints (FaceList newFaces)
{
Vertex vtxNext = unclaimed.first();
for (Vertex vtx=vtxNext; vtx!=null; vtx=vtxNext)
{ vtxNext = vtx.next;
double maxDist = tolerance;
Face maxFace = null;
for (Face newFace=newFaces.first(); newFace != null;
newFace=newFace.next)
{
if (newFace.mark == Face.VISIBLE)
{ double dist = newFace.distanceToPlane(vtx.pnt);
if (dist > 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");
}
}
}
}
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);
}
}
}
}
}
private static final int NONCONVEX_WRT_LARGER_FACE = 1;
private static final int NONCONVEX = 2;
protected double oppFaceDistance (HalfEdge he)
{
return he.face.distanceToPlane (he.opposite.face.getCentroid());
}
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<numd; i++)
{ deleteFacePoints (discardedFaces[i], face);
}
if (debug)
{ System.out.println (
" result: " + face.getVertexString());
}
return true;
}
hedge = hedge.next;
}
while (hedge != face.he0);
if (!convex)
{ face.mark = Face.NON_CONVEX;
}
return false;
}
protected void calculateHorizon (Point3d eyePnt, HalfEdge edge0, Face face, List<HalfEdge> horizon)
{
// oldFaces.add (face);
deleteFacePoints (face, null);
face.mark = Face.DELETED;
if (debug)
{ System.out.println (" visiting face " + face.getVertexString());
}
HalfEdge edge;
if (edge0 == null)
{ edge0 = face.getEdge(0);
edge = edge0;
}
else
{ edge = edge0.getNext();
}
do
{ Face oppFace = edge.oppositeFace();
if (oppFace.mark == Face.VISIBLE)
{ if (oppFace.distanceToPlane (eyePnt) > 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);
}
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);
}
protected void addNewFaces (FaceList newFaces, Vertex eyeVtx, List<HalfEdge> horizon)
{
newFaces.clear();
HalfEdge hedgeSidePrev = null;
HalfEdge hedgeSideBegin = null;
for (HalfEdge horizonHe : horizon )
{
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;
}
hedgeSideBegin.next.setOpposite (hedgeSidePrev);
}
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;
}
}
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);
}
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");
}
}
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);
}
protected void reindexFacesAndVertices()
{
for (int i=0; i<numPoints; i++)
{ pointBuffer[i].index = -1;
}
// remove inactive faces and mark active vertices
numFaces = 0;
for (Iterator it=faces.iterator(); it.hasNext(); )
{ Face face = (Face)it.next();
if (face.mark != Face.VISIBLE)
{ it.remove();
}
else
{ markFaceVertices (face, 0);
numFaces++;
}
}
// reindex vertices
numVertices = 0;
for (int i=0; i<numPoints; i++)
{ Vertex vtx = pointBuffer[i];
if (vtx.index == 0)
{ vertexPointIndices[numVertices] = i;
vtx.index = numVertices++;
}
}
}
protected boolean checkFaceConvexity (
Face face, double tol, PrintStream ps)
{
double dist;
HalfEdge he = face.he0;
do
{ face.checkConsistency();
// make sure edge is convex
dist = oppFaceDistance (he);
if (dist > 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;
}
protected boolean checkFaces(double tol, PrintStream ps)
{
// check edge convexity
boolean convex = true;
for (Face face : faces)
{
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 <code>null</code> 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 <i>tol</i>.
* 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 <i>tol</i> 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 <i>tol</i>.
*
* <p>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 <code>null</code> 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<numPoints; i++)
{ Point3d pnt = pointBuffer[i].pnt;
for (Face face : faces)
{
if (face.mark == Face.VISIBLE)
{
dist = face.distanceToPlane (pnt);
if (dist > pointTol)
{ if (ps != null)
{ ps.println (
"Point " + i + " " + dist + " above face " +
face.getVertexString());
}
return false;
}
}
}
}
return true;
}
}