// ******************************************************************************
   //
   // Title:       Force Field X.
   // Description: Force Field X - Software for Molecular Biophysics.
   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
   //
   // This file is part of Force Field X.
   //
   // Force Field X is free software; you can redistribute it and/or modify it
  // under the terms of the GNU General Public License version 3 as published by
  // the Free Software Foundation.
  //
  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  // details.
  //
  // You should have received a copy of the GNU General Public License along with
  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
  // Place, Suite 330, Boston, MA 02111-1307 USA
  //
  // Linking this library statically or dynamically with other modules is making a
  // combined work based on this library. Thus, the terms and conditions of the
  // GNU General Public License cover the whole combination.
  //
  // As a special exception, the copyright holders of this library give you
  // permission to link this library with independent modules to produce an
  // executable, regardless of the license terms of these independent modules, and
  // to copy and distribute the resulting executable under terms of your choice,
  // provided that you also meet, for each linked independent module, the terms
  // and conditions of the license of that module. An independent module is a
  // module which is not derived from or based on this library. If you modify this
  // library, you may extend this exception to your version of the library, but
  // you are not obligated to do so. If you do not wish to do so, delete this
  // exception statement from your version.
  //
  // ******************************************************************************
  package ffx.numerics.quickhull;
  
  import java.io.PrintStream;
  import java.util.Arrays;
  import java.util.Iterator;
  import java.util.Vector;
  import java.util.logging.Level;
  import java.util.logging.Logger;
  
  /**
   * 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(&quot;Vertices:&quot;);
   * Point3d[] vertices = hull.getVertices();
   * for (int i = 0; i &lt; vertices.length; i++) {
   *     Point3d pnt = vertices[i];
   *     System.out.println(pnt.x + &quot; &quot; + pnt.y + &quot; &quot; + pnt.z);
   * }
   *
   * System.out.println(&quot;Faces:&quot;);
   * int[][] faceIndices = hull.getFaces();
   * for (int i = 0; i &lt; faceIndices.length; i++) {
   *     for (int k = 0; k &lt; faceIndices[i].length; k++) {
   *         System.out.print(faceIndices[i][k] + &quot; &quot;);
   *     }
   *     System.out.println(&quot;&quot;);
   * }
   * </pre>
   * <p>
   * 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.
   * <p>
  * <b>Robustness</b>
  * <p>
  * 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.
  * <p>
  * <b>Merged Faces</b>
  * <p>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>).
  * <p>
  * <b>Degenerate Input</b>
  * <p>
  * 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
  * @author Michael J. Schnieders
  * @since 1.0
  */
 public class QuickHull3D {
 
   /**
    * Logger to log to.
    */
   private static final Logger logger = Logger.getLogger(QuickHull3D.class.getName());
 
   /**
    * 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 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 Vector<Face> faces = new Vector<>(16);
 
   protected Vector<HalfEdge> horizon = new Vector<>(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;
 
   /**
    * Precision of a double.
    */
   private static 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);
   }
 
   /**
    * Finds the half-edge within the current hull whose tail/head match the given vertices.
    * Used during setHull to establish opposite half-edge links between faces.
    *
    * @param tail the tail vertex of the desired half-edge
    * @param head the head vertex of the desired half-edge
    * @return the matching HalfEdge, or null if none is found
    */
   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;
   }
 
   /**
    * Initializes the hull from precomputed face indices and point coordinates.
    * This is primarily used for testing or reconstructing a hull state.
    *
    * @param coords      flat array of x,y,z coordinates (length >= 3*nump)
    * @param nump        number of points in coords
    * @param faceIndices face vertex indices (each row lists vertices of a face in CCW order)
    * @param numf        number of faces to add
    */
   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);
     }
   }
 
   /**
    * print all points to the print stream (very point a line)
    *
    * @param ps the print stream to write to
    */
   public 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);
     }
   }
 
   /**
    * (Re)initializes internal buffers and lists for a hull build of nump points.
    * Ensures pointBuffer capacity, clears face/claim lists, and sets counters.
    *
    * @param nump number of input points to prepare for
    */
   protected void initBuffers(int nump) {
     if (pointBuffer.length < nump) {
       Vertex[] newBuffer = new Vertex[nump];
       vertexPointIndices = new int[nump];
       System.arraycopy(pointBuffer, 0, newBuffer, 0, pointBuffer.length);
       for (int i = pointBuffer.length; i < nump; i++) {
         newBuffer[i] = new Vertex();
       }
       pointBuffer = newBuffer;
     }
     faces.clear();
     claimed.clear();
     numFaces = 0;
     numPoints = nump;
   }
 
   /**
    * Populates the internal vertex buffer from a flat coordinate array.
    *
    * @param coords flat array of x,y,z coordinates (length >= 3*nump)
    * @param nump   number of points to load
    */
   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;
     }
   }
 
   /**
    * Populates the internal vertex buffer from an array of Point3d.
    *
    * @param pnts array of Point3d
    * @param nump number of points to load
    */
   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;
     }
   }
 
   /**
    * Computes per-axis min/max vertices and sets characteristic length and tolerance.
    * Used to seed the initial simplex and numerical thresholds.
    */
   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;
         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 tetrahedral simplex from extremal points, throwing if points are
    * coincident/colinear/coplanar within the current tolerance.
    *
    * @throws IllegalArgumentException if a non-degenerate simplex cannot be formed
    */
   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 (logger.isLoggable(Level.FINE)) {
       logger.fine("initial vertices:");
       logger.fine(vtx[0].index + ": " + vtx[0].pnt);
       logger.fine(vtx[1].index + ": " + vtx[1].pnt);
       logger.fine(vtx[2].index + ": " + vtx[2].pnt);
       logger.fine(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));
       }
     }
 
     faces.addAll(Arrays.asList(tris).subList(0, 4));
 
     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];
     System.arraycopy(vertexPointIndices, 0, indices, 0, numVertices);
     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 index : indices) {
         ps.print(" " + index);
       }
       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);
   }
 
   /**
    * Assigns previously unclaimed vertices to the newly created faces most distant above them.
    *
    * @param newFaces list of newly created faces during the current expansion
    */
   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 (logger.isLoggable(Level.FINE) && vtx.index == findIndex) {
           logger.fine(findIndex + " CLAIMED BY " + maxFace.getVertexString());
         }
       } else {
         if (logger.isLoggable(Level.FINE) && vtx.index == findIndex) {
           logger.fine(findIndex + " DISCARDED");
         }
       }
     }
   }
 
   /**
    * Removes all outside vertices from a face and either discards them as unclaimed
    * or reassigns them to an absorbing face if they lie above it.
    *
    * @param face          face whose outside vertices are being removed
    * @param absorbingFace face to receive reassigned vertices; may be null
    */
   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;
 
   /**
    * Computes the signed distance from the centroid of the opposite face across this edge
    * to the plane of this edge's face. Positive indicates non-convexity.
    *
    * @param he the half-edge whose adjacent faces are considered
    * @return signed distance of opposite face centroid to this face plane
    */
   protected double oppFaceDistance(HalfEdge he) {
     return he.face.distanceToPlane(he.opposite.face.getCentroid());
   }
 
   /**
    * Attempts to merge a face with one of its adjacent faces if they are sufficiently
    * coplanar or non-convex per the specified mergeType.
    *
    * @param face      the face under consideration
    * @param mergeType NONCONVEX or NONCONVEX_WRT_LARGER_FACE
    * @return true if a merge occurred; false otherwise
    */
   private boolean doAdjacentMerge(Face face, int mergeType) {
     HalfEdge hedge = face.he0;
 
     boolean convex = true;
     do {
       Face oppFace = hedge.oppositeFace();
       boolean merge = false;
 
       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 (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 (logger.isLoggable(Level.FINE)) {
           logger.fine("  merging " + face.getVertexString() + "  and  " + oppFace.getVertexString());
         }
 
         int numd = face.mergeAdjacentFace(hedge, discardedFaces);
         for (int i = 0; i < numd; i++) {
           deleteFacePoints(discardedFaces[i], face);
         }
         if (logger.isLoggable(Level.FINE)) {
           logger.fine("  result: " + face.getVertexString());
         }
         return true;
       }
       hedge = hedge.next;
     } while (hedge != face.he0);
     if (!convex) {
       face.mark = Face.NON_CONVEX;
     }
     return false;
   }
 
   /**
    * Recursively computes the horizon (boundary) edges visible from a given eye point,
    * marking visited faces deleted and collecting the bordering edges.
    *
    * @param eyePnt  the point from which visibility is tested
    * @param edge0   starting edge within the face (may be null to start at edge 0)
    * @param face    the current face being visited
    * @param horizon output collection of horizon half-edges
    */
   protected void calculateHorizon(Point3d eyePnt, HalfEdge edge0, Face face, Vector<HalfEdge> horizon) {
     // oldFaces.add (face);
     deleteFacePoints(face, null);
     face.mark = Face.DELETED;
     if (logger.isLoggable(Level.FINE)) {
       logger.fine("  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 (logger.isLoggable(Level.FINE)) {
             logger.fine("  adding horizon edge " + edge.getVertexString());
           }
         }
       }
       edge = edge.getNext();
     } while (edge != edge0);
   }
 
   /**
    * Creates a new triangular face adjoining an existing horizon edge using the eye vertex,
    * links it into the hull, and returns the half-edge from eyeVtx to he.tail().
    *
    * @param eyeVtx the vertex being added to the hull
    * @param he     a horizon half-edge to adjoin
    * @return the half-edge associated with the new face for chaining
    */
   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);
   }
 
   /**
    * Builds and links the ring of new faces around the horizon using the eye vertex,
    * recording them in the provided newFaces list.
    *
    * @param newFaces list that will collect the created faces
    * @param eyeVtx   the vertex being added to the hull
    * @param horizon  ordered list of horizon half-edges
    */
   protected void addNewFaces(FaceList newFaces, Vertex eyeVtx, Vector<HalfEdge> horizon) {
     newFaces.clear();
 
     HalfEdge hedgeSidePrev = null;
     HalfEdge hedgeSideBegin = null;
 
     for (HalfEdge horizonHe : horizon) {
       HalfEdge hedgeSide = addAdjoiningFace(eyeVtx, horizonHe);
       if (logger.isLoggable(Level.FINE)) {
         logger.fine("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);
   }
 
   /**
    * Selects the next vertex to add: the farthest point above some face that currently
    * has outside points claimed.
    *
    * @return the next Vertex to add to the hull, or null if none remain
    */
   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;
     }
   }
 
   /**
    * Incorporates the specified eye vertex into the hull: computes the horizon, creates
    * new faces, merges non-convex adjacencies, and resolves/reassigns outside points.
    *
    * @param eyeVtx the vertex to add to the hull
    */
   protected void addPointToHull(Vertex eyeVtx) {
     horizon.clear();
     unclaimed.clear();
 
     if (logger.isLoggable(Level.FINE)) {
       logger.fine("Adding point: " + eyeVtx.index);
       logger.fine(" 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)) {
           // Do nothing.
         }
       }
     }
     // 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)) {
           // Do nothing.
         }
       }
     }
     resolveUnclaimedPoints(newFaces);
   }
 
   /**
    * Builds the full convex hull by repeatedly selecting and adding extreme points
    * until no outside points remain; then reindexes faces and vertices.
    */
   protected void buildHull() {
     int cnt = 0;
     Vertex eyeVtx;
 
     computeMaxAndMin();
     createInitialSimplex();
     while ((eyeVtx = nextPointToAdd()) != null) {
       addPointToHull(eyeVtx);
       cnt++;
       logger.fine("iteration " + cnt + " done");
     }
     reindexFacesAndVertices();
     logger.fine("hull done");
   }
 
   /**
    * Marks all vertices around a face's boundary with the provided integer mark value.
    *
    * @param face the face whose boundary vertices will be marked
    * @param mark the mark value to assign
    */
   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);
   }
 
   /**
    * Removes deleted faces, marks active vertices, and assigns new contiguous indices
    * to both faces and vertices for compact output.
    */
   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<Face> it = faces.iterator(); it.hasNext(); ) {
       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++;
       }
     }
   }
 
   /**
    * Verifies that a face is locally convex by checking distances between opposite face
    * centroids and face planes, and ensuring no redundant vertices exist.
    *
    * @param face the face to validate
    * @param tol  distance tolerance for convexity checks
    * @param ps   optional PrintStream for diagnostics (null to suppress)
    * @return true if the face passes local convexity checks
    */
   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;
   }
 
   /**
    * Performs convexity validation across all visible faces using a specified tolerance.
    *
    * @param tol distance tolerance for convexity checks
    * @param ps  optional PrintStream for diagnostics (null to suppress)
    * @return true if all visible faces pass convexity checks
    */
   protected boolean checkFaces(double tol, PrintStream ps) {
     // check edge convexity
     boolean convex = true;
     for (Face face : faces) {
       if (face.mark == Face.VISIBLE && !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;
   }
 }