From 3aa338492de1780e8730f206eb9957f4ac442ac4 Mon Sep 17 00:00:00 2001
From: Thomas <thomas.musset@pasteur.fr>
Date: Mon, 3 Jul 2023 16:14:56 +0200
Subject: [PATCH] updated pom, optimized code for java 11
---
pom.xml | 2 +-
.../java/plugins/adufour/quickhull/Face.java | 999 +++---
.../plugins/adufour/quickhull/FaceList.java | 85 +-
.../plugins/adufour/quickhull/HalfEdge.java | 461 ++-
.../plugins/adufour/quickhull/QuickHull.java | 28 +-
.../adufour/quickhull/QuickHull2D.java | 127 +-
.../adufour/quickhull/QuickHull3D.java | 2670 ++++++++---------
.../plugins/adufour/quickhull/Vertex.java | 102 +-
.../plugins/adufour/quickhull/VertexList.java | 256 +-
9 files changed, 2236 insertions(+), 2494 deletions(-)
diff --git a/pom.xml b/pom.xml
index 8e9aeec..b5454c5 100644
--- a/pom.xml
+++ b/pom.xml
@@ -11,7 +11,7 @@
</parent>
<artifactId>quickhull</artifactId>
- <version>1.0.3</version>
+ <version>2.0.0</version>
<packaging>jar</packaging>
diff --git a/src/main/java/plugins/adufour/quickhull/Face.java b/src/main/java/plugins/adufour/quickhull/Face.java
index e38ee98..042eea7 100644
--- a/src/main/java/plugins/adufour/quickhull/Face.java
+++ b/src/main/java/plugins/adufour/quickhull/Face.java
@@ -1,539 +1,486 @@
+/*
+ * Copyright (c) 2010-2023. Institut Pasteur.
+ *
+ * This file is part of Icy.
+ * Icy is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Icy 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 Icy. If not, see <https://www.gnu.org/licenses/>.
+ */
+
package plugins.adufour.quickhull;
import javax.vecmath.Vector3d;
/**
* Basic triangular face used to form the hull.
- *
+ *
* <p>
* The information stored for each face consists of a planar normal, a planar offset, and a
* doubly-linked list of three <a href=HalfEdge>HalfEdges</a> which surround the face in a
* counter-clockwise direction.
- *
+ *
* @author John E. Lloyd, Fall 2004
*/
-class Face
-{
- HalfEdge he0;
- private Vector3d normal;
- double area;
- private Vector3d centroid;
- double planeOffset;
- int index;
- int numVerts;
-
- Face next;
-
- static final int VISIBLE = 1;
- static final int NON_CONVEX = 2;
- static final int DELETED = 3;
-
- int mark = VISIBLE;
-
- Vertex outside;
-
- public void computeCentroid(Vector3d centroid)
- {
- centroid.set(0, 0, 0);
- HalfEdge he = he0;
- do
- {
- centroid.add(he.head().pnt);
- he = he.next;
- } while (he != he0);
- centroid.scale(1 / (double) numVerts);
- }
-
- public void computeNormal(Vector3d normal, double minArea)
- {
- computeNormal(normal);
-
- if (area < minArea)
- {
- System.out.println("area=" + area);
- // make the normal more robust by removing
- // components parallel to the longest edge
-
- HalfEdge hedgeMax = null;
- double lenSqrMax = 0;
- HalfEdge hedge = he0;
- do
- {
- double lenSqr = hedge.lengthSquared();
- if (lenSqr > lenSqrMax)
- {
- hedgeMax = hedge;
- lenSqrMax = lenSqr;
- }
- hedge = hedge.next;
- } while (hedge != he0);
-
- Vector3d p2 = hedgeMax.head().pnt;
- Vector3d p1 = hedgeMax.tail().pnt;
- double lenMax = Math.sqrt(lenSqrMax);
- double ux = (p2.x - p1.x) / lenMax;
- double uy = (p2.y - p1.y) / lenMax;
- double uz = (p2.z - p1.z) / lenMax;
- double dot = normal.x * ux + normal.y * uy + normal.z * uz;
- normal.x -= dot * ux;
- normal.y -= dot * uy;
- normal.z -= dot * uz;
-
- normal.normalize();
- }
- }
-
- public void computeNormal(Vector3d normal)
- {
- HalfEdge he1 = he0.next;
- HalfEdge he2 = he1.next;
-
- Vector3d p0 = he0.head().pnt;
- Vector3d p2 = he1.head().pnt;
-
- double d2x = p2.x - p0.x;
- double d2y = p2.y - p0.y;
- double d2z = p2.z - p0.z;
-
- normal.set(0, 0, 0);
-
- numVerts = 2;
-
- while (he2 != he0)
- {
- double d1x = d2x;
- double d1y = d2y;
- double d1z = d2z;
-
- p2 = he2.head().pnt;
- d2x = p2.x - p0.x;
- d2y = p2.y - p0.y;
- d2z = p2.z - p0.z;
-
- normal.x += d1y * d2z - d1z * d2y;
- normal.y += d1z * d2x - d1x * d2z;
- normal.z += d1x * d2y - d1y * d2x;
-
- he1 = he2;
- he2 = he2.next;
- numVerts++;
- }
- area = normal.length();
- normal.scale(1 / area);
- }
-
- private void computeNormalAndCentroid()
- {
- computeNormal(normal);
- computeCentroid(centroid);
- planeOffset = normal.dot(centroid);
- int numv = 0;
- HalfEdge he = he0;
- do
- {
- numv++;
- he = he.next;
- } while (he != he0);
- if (numv != numVerts)
- {
- throw new RuntimeException("face " + getVertexString() + " numVerts=" + numVerts + " should be " + numv);
- }
- }
-
- private void computeNormalAndCentroid(double minArea)
- {
- computeNormal(normal, minArea);
- computeCentroid(centroid);
- planeOffset = normal.dot(centroid);
- }
-
- public static Face createTriangle(Vertex v0, Vertex v1, Vertex v2)
- {
- return createTriangle(v0, v1, v2, 0);
- }
-
- /**
- * Constructs a triangule Face from vertices v0, v1, and v2.
- *
- * @param v0
- * first vertex
- * @param v1
- * second vertex
- * @param v2
- * third vertex
- */
- public static Face createTriangle(Vertex v0, Vertex v1, Vertex v2, double minArea)
- {
- Face face = new Face();
- HalfEdge he0 = new HalfEdge(v0, face);
- HalfEdge he1 = new HalfEdge(v1, face);
- HalfEdge he2 = new HalfEdge(v2, face);
-
- he0.prev = he2;
- he0.next = he1;
- he1.prev = he0;
- he1.next = he2;
- he2.prev = he1;
- he2.next = he0;
-
- face.he0 = he0;
-
- // compute the normal and offset
- face.computeNormalAndCentroid(minArea);
- return face;
- }
-
- public static Face create(Vertex[] vtxArray, int[] indices)
- {
- Face face = new Face();
- HalfEdge hePrev = null;
- for (int i = 0; i < indices.length; i++)
- {
- HalfEdge he = new HalfEdge(vtxArray[indices[i]], face);
- if (hePrev != null)
- {
- he.setPrev(hePrev);
- hePrev.setNext(he);
- }
- else
- {
- face.he0 = he;
- }
- hePrev = he;
- }
- face.he0.setPrev(hePrev);
- hePrev.setNext(face.he0);
-
- // compute the normal and offset
- face.computeNormalAndCentroid();
- return face;
- }
-
- public Face()
- {
- normal = new Vector3d();
- centroid = new Vector3d();
- mark = VISIBLE;
- }
-
- /**
- * Gets the i-th half-edge associated with the face.
- *
- * @param i
- * the half-edge index, in the range 0-2.
- * @return the half-edge
- */
- public HalfEdge getEdge(int i)
- {
- HalfEdge he = he0;
- while (i > 0)
- {
- he = he.next;
- i--;
- }
- while (i < 0)
- {
- he = he.prev;
- i++;
- }
- return he;
- }
-
- public HalfEdge getFirstEdge()
- {
- return he0;
- }
-
- /**
- * Finds the half-edge within this face which has tail <code>vt</code> and head
- * <code>vh</code>.
- *
- * @param vt
- * tail point
- * @param vh
- * head point
- * @return the half-edge, or null if none is found.
- */
- public HalfEdge findEdge(Vertex vt, Vertex vh)
- {
- HalfEdge he = he0;
- do
- {
- if (he.head() == vh && he.tail() == vt)
- {
- return he;
- }
- he = he.next;
- } while (he != he0);
- return null;
- }
-
- /**
- * Computes the distance from a point p to the plane of this face.
- *
- * @param p
- * the point
- * @return distance from the point to the plane
- */
- public double distanceToPlane(Vector3d p)
- {
- return normal.x * p.x + normal.y * p.y + normal.z * p.z - planeOffset;
- }
-
- /**
- * Returns the normal of the plane associated with this face.
- *
- * @return the planar normal
- */
- public Vector3d getNormal()
- {
- return normal;
- }
-
- public Vector3d getCentroid()
- {
- return centroid;
- }
-
- public int numVertices()
- {
- return numVerts;
- }
-
- public String getVertexString()
- {
- String s = null;
- HalfEdge he = he0;
- do
- {
- if (s == null)
- {
- s = "" + he.head().index;
- }
- else
- {
- s += " " + he.head().index;
- }
- he = he.next;
- } while (he != he0);
- return s;
- }
-
- public void getVertexIndices(int[] idxs)
- {
- HalfEdge he = he0;
- int i = 0;
- do
- {
- idxs[i++] = he.head().index;
- he = he.next;
- } while (he != he0);
- }
-
- private Face connectHalfEdges(HalfEdge hedgePrev, HalfEdge hedge)
- {
- Face discardedFace = null;
-
- if (hedgePrev.oppositeFace() == hedge.oppositeFace())
- { // then there is a redundant edge that we can get rid off
-
- Face oppFace = hedge.oppositeFace();
- HalfEdge hedgeOpp;
-
- if (hedgePrev == he0)
- {
- he0 = hedge;
- }
- if (oppFace.numVertices() == 3)
- { // then we can get rid of the opposite face altogether
- hedgeOpp = hedge.getOpposite().prev.getOpposite();
-
- oppFace.mark = DELETED;
- discardedFace = oppFace;
- }
- else
- {
- hedgeOpp = hedge.getOpposite().next;
-
- if (oppFace.he0 == hedgeOpp.prev)
- {
- oppFace.he0 = hedgeOpp;
- }
- hedgeOpp.prev = hedgeOpp.prev.prev;
- hedgeOpp.prev.next = hedgeOpp;
- }
- hedge.prev = hedgePrev.prev;
- hedge.prev.next = hedge;
-
- hedge.opposite = hedgeOpp;
- hedgeOpp.opposite = hedge;
-
- // oppFace was modified, so need to recompute
- oppFace.computeNormalAndCentroid();
- }
- else
- {
- hedgePrev.next = hedge;
- hedge.prev = hedgePrev;
- }
- return discardedFace;
- }
-
- void checkConsistency()
- {
- // do a sanity check on the face
- HalfEdge hedge = he0;
- double maxd = 0;
- int numv = 0;
-
- if (numVerts < 3)
- {
- throw new RuntimeException("degenerate face: " + getVertexString());
- }
- do
- {
- HalfEdge hedgeOpp = hedge.getOpposite();
- if (hedgeOpp == null)
- {
- throw new RuntimeException("face " + getVertexString() + ": " + "unreflected half edge " + hedge.getVertexString());
- }
- else if (hedgeOpp.getOpposite() != hedge)
- {
- throw new RuntimeException("face " + getVertexString() + ": " + "opposite half edge " + hedgeOpp.getVertexString() + " has opposite " + hedgeOpp.getOpposite().getVertexString());
- }
- if (hedgeOpp.head() != hedge.tail() || hedge.head() != hedgeOpp.tail())
- {
- throw new RuntimeException("face " + getVertexString() + ": " + "half edge " + hedge.getVertexString() + " reflected by " + hedgeOpp.getVertexString());
- }
- Face oppFace = hedgeOpp.face;
- if (oppFace == null)
- {
- throw new RuntimeException("face " + getVertexString() + ": " + "no face on half edge " + hedgeOpp.getVertexString());
- }
- else if (oppFace.mark == DELETED)
- {
- throw new RuntimeException("face " + getVertexString() + ": " + "opposite face " + oppFace.getVertexString() + " not on hull");
- }
- double d = Math.abs(distanceToPlane(hedge.head().pnt));
- if (d > maxd)
- {
- maxd = d;
- }
- numv++;
- hedge = hedge.next;
- } while (hedge != he0);
-
- if (numv != numVerts)
- {
- throw new RuntimeException("face " + getVertexString() + " numVerts=" + numVerts + " should be " + numv);
- }
-
- }
-
- public int mergeAdjacentFace(HalfEdge hedgeAdj, Face[] discarded)
- {
- Face oppFace = hedgeAdj.oppositeFace();
- int numDiscarded = 0;
-
- discarded[numDiscarded++] = oppFace;
- oppFace.mark = DELETED;
-
- HalfEdge hedgeOpp = hedgeAdj.getOpposite();
-
- HalfEdge hedgeAdjPrev = hedgeAdj.prev;
- HalfEdge hedgeAdjNext = hedgeAdj.next;
- HalfEdge hedgeOppPrev = hedgeOpp.prev;
- HalfEdge hedgeOppNext = hedgeOpp.next;
-
- while (hedgeAdjPrev.oppositeFace() == oppFace)
- {
- hedgeAdjPrev = hedgeAdjPrev.prev;
- hedgeOppNext = hedgeOppNext.next;
- }
-
- while (hedgeAdjNext.oppositeFace() == oppFace)
- {
- hedgeOppPrev = hedgeOppPrev.prev;
- hedgeAdjNext = hedgeAdjNext.next;
- }
-
- HalfEdge hedge;
-
- for (hedge = hedgeOppNext; hedge != hedgeOppPrev.next; hedge = hedge.next)
- {
- hedge.face = this;
- }
-
- if (hedgeAdj == he0)
- {
- he0 = hedgeAdjNext;
- }
-
- // handle the half edges at the head
- Face discardedFace;
-
- discardedFace = connectHalfEdges(hedgeOppPrev, hedgeAdjNext);
- if (discardedFace != null)
- {
- discarded[numDiscarded++] = discardedFace;
- }
-
- // handle the half edges at the tail
- discardedFace = connectHalfEdges(hedgeAdjPrev, hedgeOppNext);
- if (discardedFace != null)
- {
- discarded[numDiscarded++] = discardedFace;
- }
-
- computeNormalAndCentroid();
- checkConsistency();
-
- return numDiscarded;
- }
-
- public void triangulate(FaceList newFaces, double minArea)
- {
- HalfEdge hedge;
-
- if (numVertices() < 4)
- {
- return;
- }
-
- Vertex v0 = he0.head();
-
- hedge = he0.next;
- HalfEdge oppPrev = hedge.opposite;
- Face face0 = null;
-
- for (hedge = hedge.next; hedge != he0.prev; hedge = hedge.next)
- {
- Face face = createTriangle(v0, hedge.prev.head(), hedge.head(), minArea);
- face.he0.next.setOpposite(oppPrev);
- face.he0.prev.setOpposite(hedge.opposite);
- oppPrev = face.he0;
- newFaces.add(face);
- if (face0 == null)
- {
- face0 = face;
- }
- }
- hedge = new HalfEdge(he0.prev.prev.head(), this);
- hedge.setOpposite(oppPrev);
-
- hedge.prev = he0;
- hedge.prev.next = hedge;
-
- hedge.next = he0.prev;
- hedge.next.prev = hedge;
-
- computeNormalAndCentroid(minArea);
- checkConsistency();
-
- for (Face face = face0; face != null; face = face.next)
- {
- face.checkConsistency();
- }
-
- }
+class Face {
+ HalfEdge he0;
+ private final Vector3d normal;
+ double area;
+ private final Vector3d centroid;
+ double planeOffset;
+ int index;
+ int numVerts;
+
+ Face next;
+
+ static final int VISIBLE = 1;
+ static final int NON_CONVEX = 2;
+ static final int DELETED = 3;
+
+ int mark = VISIBLE;
+
+ Vertex outside;
+
+ public void computeCentroid(final Vector3d centroid) {
+ centroid.set(0, 0, 0);
+ HalfEdge he = he0;
+ do {
+ centroid.add(he.head().pnt);
+ he = he.next;
+ } while (he != he0);
+ centroid.scale(1 / (double) numVerts);
+ }
+
+ public void computeNormal(final Vector3d normal, final double minArea) {
+ computeNormal(normal);
+
+ if (area < minArea) {
+ System.out.println("area=" + area);
+ // make the normal more robust by removing
+ // components parallel to the longest edge
+
+ HalfEdge hedgeMax = null;
+ double lenSqrMax = 0;
+ HalfEdge hedge = he0;
+ do {
+ final double lenSqr = hedge.lengthSquared();
+ if (lenSqr > lenSqrMax) {
+ hedgeMax = hedge;
+ lenSqrMax = lenSqr;
+ }
+ hedge = hedge.next;
+ } while (hedge != he0);
+
+ assert hedgeMax != null;
+ final Vector3d p2 = hedgeMax.head().pnt;
+ final Vector3d p1 = hedgeMax.tail().pnt;
+ final double lenMax = Math.sqrt(lenSqrMax);
+ final double ux = (p2.x - p1.x) / lenMax;
+ final double uy = (p2.y - p1.y) / lenMax;
+ final double uz = (p2.z - p1.z) / lenMax;
+ final double dot = normal.x * ux + normal.y * uy + normal.z * uz;
+ normal.x -= dot * ux;
+ normal.y -= dot * uy;
+ normal.z -= dot * uz;
+
+ normal.normalize();
+ }
+ }
+
+ public void computeNormal(final Vector3d normal) {
+ HalfEdge he1 = he0.next;
+ HalfEdge he2 = he1.next;
+
+ final Vector3d p0 = he0.head().pnt;
+ Vector3d p2 = he1.head().pnt;
+
+ double d2x = p2.x - p0.x;
+ double d2y = p2.y - p0.y;
+ double d2z = p2.z - p0.z;
+
+ normal.set(0, 0, 0);
+
+ numVerts = 2;
+
+ while (he2 != he0) {
+ final double d1x = d2x;
+ final double d1y = d2y;
+ final double d1z = d2z;
+
+ p2 = he2.head().pnt;
+ d2x = p2.x - p0.x;
+ d2y = p2.y - p0.y;
+ d2z = p2.z - p0.z;
+
+ normal.x += d1y * d2z - d1z * d2y;
+ normal.y += d1z * d2x - d1x * d2z;
+ normal.z += d1x * d2y - d1y * d2x;
+
+ he1 = he2;
+ he2 = he2.next;
+ numVerts++;
+ }
+ area = normal.length();
+ normal.scale(1 / area);
+ }
+
+ private void computeNormalAndCentroid() {
+ computeNormal(normal);
+ computeCentroid(centroid);
+ planeOffset = normal.dot(centroid);
+ int numv = 0;
+ HalfEdge he = he0;
+ do {
+ numv++;
+ he = he.next;
+ } while (he != he0);
+ if (numv != numVerts) {
+ throw new RuntimeException("face " + getVertexString() + " numVerts=" + numVerts + " should be " + numv);
+ }
+ }
+
+ private void computeNormalAndCentroid(final double minArea) {
+ computeNormal(normal, minArea);
+ computeCentroid(centroid);
+ planeOffset = normal.dot(centroid);
+ }
+
+ public static Face createTriangle(final Vertex v0, final Vertex v1, final Vertex v2) {
+ return createTriangle(v0, v1, v2, 0);
+ }
+
+ /**
+ * Constructs a triangule Face from vertices v0, v1, and v2.
+ *
+ * @param v0 first vertex
+ * @param v1 second vertex
+ * @param v2 third vertex
+ */
+ public static Face createTriangle(final Vertex v0, final Vertex v1, final Vertex v2, final double minArea) {
+ final Face face = new Face();
+ final HalfEdge he0 = new HalfEdge(v0, face);
+ final HalfEdge he1 = new HalfEdge(v1, face);
+ final HalfEdge he2 = new HalfEdge(v2, face);
+
+ he0.prev = he2;
+ he0.next = he1;
+ he1.prev = he0;
+ he1.next = he2;
+ he2.prev = he1;
+ he2.next = he0;
+
+ face.he0 = he0;
+
+ // compute the normal and offset
+ face.computeNormalAndCentroid(minArea);
+ return face;
+ }
+
+ public static Face create(final Vertex[] vtxArray, final int[] indices) {
+ final Face face = new Face();
+ HalfEdge hePrev = null;
+ for (final int j : indices) {
+ final HalfEdge he = new HalfEdge(vtxArray[j], face);
+ if (hePrev != null) {
+ he.setPrev(hePrev);
+ hePrev.setNext(he);
+ }
+ else {
+ face.he0 = he;
+ }
+ hePrev = he;
+ }
+ assert hePrev != null;
+ face.he0.setPrev(hePrev);
+ hePrev.setNext(face.he0);
+
+ // compute the normal and offset
+ face.computeNormalAndCentroid();
+ return face;
+ }
+
+ public Face() {
+ normal = new Vector3d();
+ centroid = new Vector3d();
+ mark = VISIBLE;
+ }
+
+ /**
+ * Gets the i-th half-edge associated with the face.
+ *
+ * @param i the half-edge index, in the range 0-2.
+ * @return the half-edge
+ */
+ public HalfEdge getEdge(int i) {
+ HalfEdge he = he0;
+ while (i > 0) {
+ he = he.next;
+ i--;
+ }
+ while (i < 0) {
+ he = he.prev;
+ i++;
+ }
+ return he;
+ }
+
+ public HalfEdge getFirstEdge() {
+ return he0;
+ }
+
+ /**
+ * Finds the half-edge within this face which has tail <code>vt</code> and head
+ * <code>vh</code>.
+ *
+ * @param vt tail point
+ * @param vh head point
+ * @return the half-edge, or null if none is found.
+ */
+ public HalfEdge findEdge(final Vertex vt, final Vertex vh) {
+ HalfEdge he = he0;
+ do {
+ if (he.head() == vh && he.tail() == vt) {
+ return he;
+ }
+ he = he.next;
+ } while (he != he0);
+ return null;
+ }
+
+ /**
+ * Computes the distance from a point p to the plane of this face.
+ *
+ * @param p the point
+ * @return distance from the point to the plane
+ */
+ public double distanceToPlane(final Vector3d p) {
+ return normal.x * p.x + normal.y * p.y + normal.z * p.z - planeOffset;
+ }
+
+ /**
+ * Returns the normal of the plane associated with this face.
+ *
+ * @return the planar normal
+ */
+ public Vector3d getNormal() {
+ return normal;
+ }
+
+ public Vector3d getCentroid() {
+ return centroid;
+ }
+
+ public int numVertices() {
+ return numVerts;
+ }
+
+ public String getVertexString() {
+ final StringBuilder s = new StringBuilder();
+ HalfEdge he = he0;
+ do {
+ if (s.length() == 0) {
+ s.append(he.head().index);
+ }
+ else {
+ s.append(" ").append(he.head().index);
+ }
+ he = he.next;
+ } while (he != he0);
+ return s.toString();
+ }
+
+ public void getVertexIndices(final int[] idxs) {
+ HalfEdge he = he0;
+ int i = 0;
+ do {
+ idxs[i++] = he.head().index;
+ he = he.next;
+ } while (he != he0);
+ }
+
+ private Face connectHalfEdges(final HalfEdge hedgePrev, final HalfEdge hedge) {
+ Face discardedFace = null;
+
+ if (hedgePrev.oppositeFace() == hedge.oppositeFace()) { // then there is a redundant edge that we can get rid off
+
+ final Face oppFace = hedge.oppositeFace();
+ final HalfEdge hedgeOpp;
+
+ if (hedgePrev == he0) {
+ he0 = hedge;
+ }
+ if (oppFace.numVertices() == 3) { // then we can get rid of the opposite face altogether
+ hedgeOpp = hedge.getOpposite().prev.getOpposite();
+
+ oppFace.mark = DELETED;
+ discardedFace = oppFace;
+ }
+ else {
+ hedgeOpp = hedge.getOpposite().next;
+
+ if (oppFace.he0 == hedgeOpp.prev) {
+ oppFace.he0 = hedgeOpp;
+ }
+ hedgeOpp.prev = hedgeOpp.prev.prev;
+ hedgeOpp.prev.next = hedgeOpp;
+ }
+ hedge.prev = hedgePrev.prev;
+ hedge.prev.next = hedge;
+
+ hedge.opposite = hedgeOpp;
+ hedgeOpp.opposite = hedge;
+
+ // oppFace was modified, so need to recompute
+ oppFace.computeNormalAndCentroid();
+ }
+ else {
+ hedgePrev.next = hedge;
+ hedge.prev = hedgePrev;
+ }
+ return discardedFace;
+ }
+
+ void checkConsistency() {
+ // do a sanity check on the face
+ HalfEdge hedge = he0;
+ double maxd = 0;
+ int numv = 0;
+
+ if (numVerts < 3) {
+ throw new RuntimeException("degenerate face: " + getVertexString());
+ }
+ do {
+ final HalfEdge hedgeOpp = hedge.getOpposite();
+ if (hedgeOpp == null) {
+ throw new RuntimeException("face " + getVertexString() + ": " + "unreflected half edge " + hedge.getVertexString());
+ }
+ else if (hedgeOpp.getOpposite() != hedge) {
+ throw new RuntimeException("face " + getVertexString() + ": " + "opposite half edge " + hedgeOpp.getVertexString() + " has opposite " + hedgeOpp.getOpposite().getVertexString());
+ }
+ if (hedgeOpp.head() != hedge.tail() || hedge.head() != hedgeOpp.tail()) {
+ throw new RuntimeException("face " + getVertexString() + ": " + "half edge " + hedge.getVertexString() + " reflected by " + hedgeOpp.getVertexString());
+ }
+ final Face oppFace = hedgeOpp.face;
+ if (oppFace == null) {
+ throw new RuntimeException("face " + getVertexString() + ": " + "no face on half edge " + hedgeOpp.getVertexString());
+ }
+ else if (oppFace.mark == DELETED) {
+ throw new RuntimeException("face " + getVertexString() + ": " + "opposite face " + oppFace.getVertexString() + " not on hull");
+ }
+ final double d = Math.abs(distanceToPlane(hedge.head().pnt));
+ if (d > maxd) {
+ maxd = d;
+ }
+ numv++;
+ hedge = hedge.next;
+ } while (hedge != he0);
+
+ if (numv != numVerts) {
+ throw new RuntimeException("face " + getVertexString() + " numVerts=" + numVerts + " should be " + numv);
+ }
+
+ }
+
+ public int mergeAdjacentFace(final HalfEdge hedgeAdj, final Face[] discarded) {
+ final Face oppFace = hedgeAdj.oppositeFace();
+ int numDiscarded = 0;
+
+ discarded[numDiscarded++] = oppFace;
+ oppFace.mark = DELETED;
+
+ final HalfEdge hedgeOpp = hedgeAdj.getOpposite();
+
+ HalfEdge hedgeAdjPrev = hedgeAdj.prev;
+ HalfEdge hedgeAdjNext = hedgeAdj.next;
+ HalfEdge hedgeOppPrev = hedgeOpp.prev;
+ HalfEdge hedgeOppNext = hedgeOpp.next;
+
+ while (hedgeAdjPrev.oppositeFace() == oppFace) {
+ hedgeAdjPrev = hedgeAdjPrev.prev;
+ hedgeOppNext = hedgeOppNext.next;
+ }
+
+ while (hedgeAdjNext.oppositeFace() == oppFace) {
+ hedgeOppPrev = hedgeOppPrev.prev;
+ hedgeAdjNext = hedgeAdjNext.next;
+ }
+
+ HalfEdge hedge;
+
+ for (hedge = hedgeOppNext; hedge != hedgeOppPrev.next; hedge = hedge.next) {
+ hedge.face = this;
+ }
+
+ if (hedgeAdj == he0) {
+ he0 = hedgeAdjNext;
+ }
+
+ // handle the half edges at the head
+ Face discardedFace;
+
+ discardedFace = connectHalfEdges(hedgeOppPrev, hedgeAdjNext);
+ if (discardedFace != null) {
+ discarded[numDiscarded++] = discardedFace;
+ }
+
+ // handle the half edges at the tail
+ discardedFace = connectHalfEdges(hedgeAdjPrev, hedgeOppNext);
+ if (discardedFace != null) {
+ discarded[numDiscarded++] = discardedFace;
+ }
+
+ computeNormalAndCentroid();
+ checkConsistency();
+
+ return numDiscarded;
+ }
+
+ public void triangulate(final FaceList newFaces, final double minArea) {
+ HalfEdge hedge;
+
+ if (numVertices() < 4) {
+ return;
+ }
+
+ final Vertex v0 = he0.head();
+
+ hedge = he0.next;
+ HalfEdge oppPrev = hedge.opposite;
+ Face face0 = null;
+
+ for (hedge = hedge.next; hedge != he0.prev; hedge = hedge.next) {
+ final Face face = createTriangle(v0, hedge.prev.head(), hedge.head(), minArea);
+ face.he0.next.setOpposite(oppPrev);
+ face.he0.prev.setOpposite(hedge.opposite);
+ oppPrev = face.he0;
+ newFaces.add(face);
+ if (face0 == null) {
+ face0 = face;
+ }
+ }
+ hedge = new HalfEdge(he0.prev.prev.head(), this);
+ hedge.setOpposite(oppPrev);
+
+ hedge.prev = he0;
+ hedge.prev.next = hedge;
+
+ hedge.next = he0.prev;
+ hedge.next.prev = hedge;
+
+ computeNormalAndCentroid(minArea);
+ checkConsistency();
+
+ for (Face face = face0; face != null; face = face.next) {
+ face.checkConsistency();
+ }
+
+ }
}
diff --git a/src/main/java/plugins/adufour/quickhull/FaceList.java b/src/main/java/plugins/adufour/quickhull/FaceList.java
index 3bcc85b..1713a12 100644
--- a/src/main/java/plugins/adufour/quickhull/FaceList.java
+++ b/src/main/java/plugins/adufour/quickhull/FaceList.java
@@ -1,46 +1,59 @@
+/*
+ * Copyright (c) 2010-2023. Institut Pasteur.
+ *
+ * This file is part of Icy.
+ * Icy is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Icy 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 Icy. If not, see <https://www.gnu.org/licenses/>.
+ */
+
package plugins.adufour.quickhull;
/**
* Maintains a single-linked list of faces for use by QuickHull3D
*/
-class FaceList
-{
- private Face head;
- private Face tail;
+class FaceList {
+ private Face head;
+ private Face tail;
- /**
- * Clears this list.
- */
- public void clear()
- {
- head = tail = null;
- }
+ /**
+ * Clears this list.
+ */
+ public void clear() {
+ head = tail = null;
+ }
- /**
- * Adds a vertex to the end of this list.
- */
- public void add (Face vtx)
- {
- if (head == null)
- { head = vtx;
- }
- else
- { tail.next = vtx;
- }
- vtx.next = null;
- tail = vtx;
- }
+ /**
+ * Adds a vertex to the end of this list.
+ */
+ public void add(final Face vtx) {
+ if (head == null) {
+ head = vtx;
+ }
+ else {
+ tail.next = vtx;
+ }
+ vtx.next = null;
+ tail = vtx;
+ }
- public Face first()
- {
- return head;
- }
+ public Face first() {
+ return head;
+ }
- /**
- * Returns true if this list is empty.
- */
- public boolean isEmpty()
- {
- return head == null;
- }
+ /**
+ * Returns true if this list is empty.
+ */
+ public boolean isEmpty() {
+ return head == null;
+ }
}
diff --git a/src/main/java/plugins/adufour/quickhull/HalfEdge.java b/src/main/java/plugins/adufour/quickhull/HalfEdge.java
index 22c6af2..1290f57 100644
--- a/src/main/java/plugins/adufour/quickhull/HalfEdge.java
+++ b/src/main/java/plugins/adufour/quickhull/HalfEdge.java
@@ -1,15 +1,19 @@
/*
- * Copyright John E. Lloyd, 2003. All rights reserved. Permission
- * to use, copy, and modify, without fee, is granted for non-commercial
- * and research purposes, provided that this copyright notice appears
- * in all copies.
+ * Copyright (c) 2010-2023. Institut Pasteur.
*
- * This software is distributed "as is", without any warranty, including
- * any implied warranty of merchantability or fitness for a particular
- * use. The authors assume 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.
+ * This file is part of Icy.
+ * Icy is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Icy 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 Icy. If not, see <https://www.gnu.org/licenses/>.
*/
package plugins.adufour.quickhull;
@@ -18,238 +22,211 @@ import javax.vecmath.Vector3d;
/**
* Represents the half-edges that surround each face in a counter-clockwise direction.
- *
+ *
* @author John E. Lloyd, Fall 2004
*/
-class HalfEdge
-{
- /**
- * The vertex associated with the head of this half-edge.
- */
- Vertex vertex;
-
- /**
- * Triangular face associated with this half-edge.
- */
- Face face;
-
- /**
- * Next half-edge in the triangle.
- */
- HalfEdge next;
-
- /**
- * Previous half-edge in the triangle.
- */
- HalfEdge prev;
-
- /**
- * Half-edge associated with the opposite triangle adjacent to this edge.
- */
- HalfEdge opposite;
-
- /**
- * Constructs a HalfEdge with head vertex <code>v</code> and left-hand triangular face
- * <code>f</code>.
- *
- * @param v
- * head vertex
- * @param f
- * left-hand triangular face
- */
- public HalfEdge(Vertex v, Face f)
- {
- vertex = v;
- face = f;
- }
-
- public HalfEdge()
- {
- }
-
- /**
- * Sets the value of the next edge adjacent (counter-clockwise) to this one within the triangle.
- *
- * @param edge
- * next adjacent edge
- */
- public void setNext(HalfEdge edge)
- {
- next = edge;
- }
-
- /**
- * Gets the value of the next edge adjacent (counter-clockwise) to this one within the triangle.
- *
- * @return next adjacent edge
- */
- public HalfEdge getNext()
- {
- return next;
- }
-
- /**
- * Sets the value of the previous edge adjacent (clockwise) to this one within the triangle.
- *
- * @param edge
- * previous adjacent edge
- */
- public void setPrev(HalfEdge edge)
- {
- prev = edge;
- }
-
- /**
- * Gets the value of the previous edge adjacent (clockwise) to this one within the triangle.
- *
- * @return previous adjacent edge
- */
- public HalfEdge getPrev()
- {
- return prev;
- }
-
- /**
- * Returns the triangular face located to the left of this half-edge.
- *
- * @return left-hand triangular face
- */
- public Face getFace()
- {
- return face;
- }
-
- /**
- * Returns the half-edge opposite to this half-edge.
- *
- * @return opposite half-edge
- */
- public HalfEdge getOpposite()
- {
- return opposite;
- }
-
- /**
- * Sets the half-edge opposite to this half-edge.
- *
- * @param edge
- * opposite half-edge
- */
- public void setOpposite(HalfEdge edge)
- {
- opposite = edge;
- edge.opposite = this;
- }
-
- /**
- * Returns the head vertex associated with this half-edge.
- *
- * @return head vertex
- */
- public Vertex head()
- {
- return vertex;
- }
-
- /**
- * Returns the tail vertex associated with this half-edge.
- *
- * @return tail vertex
- */
- public Vertex tail()
- {
- return prev != null ? prev.vertex : null;
- }
-
- /**
- * Returns the opposite triangular face associated with this half-edge.
- *
- * @return opposite triangular face
- */
- public Face oppositeFace()
- {
- return opposite != null ? opposite.face : null;
- }
-
- /**
- * Produces a string identifying this half-edge by the point index values of its tail and head
- * vertices.
- *
- * @return identifying string
- */
- public String getVertexString()
- {
- if (tail() != null)
- {
- return "" + tail().index + "-" + head().index;
- }
- else
- {
- return "?-" + head().index;
- }
- }
-
- /**
- * Returns the length of this half-edge.
- *
- * @return half-edge length
- */
- public double length()
- {
- if (tail() != null)
- {
- Vector3d v = new Vector3d();
- v.sub(head().pnt, tail().pnt);
- return v.length();
- }
- else
- {
- return -1;
- }
- }
-
- /**
- * Returns the length squared of this half-edge.
- *
- * @return half-edge length squared
- */
- public double lengthSquared()
- {
- if (tail() != null)
- {
- Vector3d v = new Vector3d();
- v.sub(head().pnt, tail().pnt);
- return v.lengthSquared();
- }
- else
- {
- return -1;
- }
- }
-
- // /**
- // * Computes nrml . (del0 X del1), where del0 and del1
- // * are the direction vectors along this halfEdge, and the
- // * halfEdge he1.
- // *
- // * A product > 0 indicates a left turn WRT the normal
- // */
- // public double turnProduct (HalfEdge he1, Vector3d nrml)
- // {
- // Point3d pnt0 = tail().pnt;
- // Point3d pnt1 = head().pnt;
- // Point3d pnt2 = he1.head().pnt;
-
- // double del0x = pnt1.x - pnt0.x;
- // double del0y = pnt1.y - pnt0.y;
- // double del0z = pnt1.z - pnt0.z;
-
- // double del1x = pnt2.x - pnt1.x;
- // double del1y = pnt2.y - pnt1.y;
- // double del1z = pnt2.z - pnt1.z;
-
- // return (nrml.x*(del0y*del1z - del0z*del1y) +
- // nrml.y*(del0z*del1x - del0x*del1z) +
- // nrml.z*(del0x*del1y - del0y*del1x));
- // }
+class HalfEdge {
+ /**
+ * The vertex associated with the head of this half-edge.
+ */
+ Vertex vertex;
+
+ /**
+ * Triangular face associated with this half-edge.
+ */
+ Face face;
+
+ /**
+ * Next half-edge in the triangle.
+ */
+ HalfEdge next;
+
+ /**
+ * Previous half-edge in the triangle.
+ */
+ HalfEdge prev;
+
+ /**
+ * Half-edge associated with the opposite triangle adjacent to this edge.
+ */
+ HalfEdge opposite;
+
+ /**
+ * Constructs a HalfEdge with head vertex <code>v</code> and left-hand triangular face
+ * <code>f</code>.
+ *
+ * @param v head vertex
+ * @param f left-hand triangular face
+ */
+ public HalfEdge(final Vertex v, final Face f) {
+ vertex = v;
+ face = f;
+ }
+
+ public HalfEdge() {
+ }
+
+ /**
+ * Sets the value of the next edge adjacent (counter-clockwise) to this one within the triangle.
+ *
+ * @param edge next adjacent edge
+ */
+ public void setNext(final HalfEdge edge) {
+ next = edge;
+ }
+
+ /**
+ * Gets the value of the next edge adjacent (counter-clockwise) to this one within the triangle.
+ *
+ * @return next adjacent edge
+ */
+ public HalfEdge getNext() {
+ return next;
+ }
+
+ /**
+ * Sets the value of the previous edge adjacent (clockwise) to this one within the triangle.
+ *
+ * @param edge previous adjacent edge
+ */
+ public void setPrev(final HalfEdge edge) {
+ prev = edge;
+ }
+
+ /**
+ * Gets the value of the previous edge adjacent (clockwise) to this one within the triangle.
+ *
+ * @return previous adjacent edge
+ */
+ public HalfEdge getPrev() {
+ return prev;
+ }
+
+ /**
+ * Returns the triangular face located to the left of this half-edge.
+ *
+ * @return left-hand triangular face
+ */
+ public Face getFace() {
+ return face;
+ }
+
+ /**
+ * Returns the half-edge opposite to this half-edge.
+ *
+ * @return opposite half-edge
+ */
+ public HalfEdge getOpposite() {
+ return opposite;
+ }
+
+ /**
+ * Sets the half-edge opposite to this half-edge.
+ *
+ * @param edge opposite half-edge
+ */
+ public void setOpposite(final HalfEdge edge) {
+ opposite = edge;
+ edge.opposite = this;
+ }
+
+ /**
+ * Returns the head vertex associated with this half-edge.
+ *
+ * @return head vertex
+ */
+ public Vertex head() {
+ return vertex;
+ }
+
+ /**
+ * Returns the tail vertex associated with this half-edge.
+ *
+ * @return tail vertex
+ */
+ public Vertex tail() {
+ return prev != null ? prev.vertex : null;
+ }
+
+ /**
+ * Returns the opposite triangular face associated with this half-edge.
+ *
+ * @return opposite triangular face
+ */
+ public Face oppositeFace() {
+ return opposite != null ? opposite.face : null;
+ }
+
+ /**
+ * Produces a string identifying this half-edge by the point index values of its tail and head
+ * vertices.
+ *
+ * @return identifying string
+ */
+ public String getVertexString() {
+ if (tail() != null) {
+ return tail().index + "-" + head().index;
+ }
+ else {
+ return "?-" + head().index;
+ }
+ }
+
+ /**
+ * Returns the length of this half-edge.
+ *
+ * @return half-edge length
+ */
+ public double length() {
+ if (tail() != null) {
+ final Vector3d v = new Vector3d();
+ v.sub(head().pnt, tail().pnt);
+ return v.length();
+ }
+ else {
+ return -1;
+ }
+ }
+
+ /**
+ * Returns the length squared of this half-edge.
+ *
+ * @return half-edge length squared
+ */
+ public double lengthSquared() {
+ if (tail() != null) {
+ final Vector3d v = new Vector3d();
+ v.sub(head().pnt, tail().pnt);
+ return v.lengthSquared();
+ }
+ else {
+ return -1;
+ }
+ }
+
+ // /**
+ // * Computes nrml . (del0 X del1), where del0 and del1
+ // * are the direction vectors along this halfEdge, and the
+ // * halfEdge he1.
+ // *
+ // * A product > 0 indicates a left turn WRT the normal
+ // */
+ // public double turnProduct (HalfEdge he1, Vector3d nrml)
+ // {
+ // Point3d pnt0 = tail().pnt;
+ // Point3d pnt1 = head().pnt;
+ // Point3d pnt2 = he1.head().pnt;
+
+ // double del0x = pnt1.x - pnt0.x;
+ // double del0y = pnt1.y - pnt0.y;
+ // double del0z = pnt1.z - pnt0.z;
+
+ // double del1x = pnt2.x - pnt1.x;
+ // double del1y = pnt2.y - pnt1.y;
+ // double del1z = pnt2.z - pnt1.z;
+
+ // return (nrml.x*(del0y*del1z - del0z*del1y) +
+ // nrml.y*(del0z*del1x - del0x*del1z) +
+ // nrml.z*(del0x*del1y - del0y*del1x));
+ // }
}
diff --git a/src/main/java/plugins/adufour/quickhull/QuickHull.java b/src/main/java/plugins/adufour/quickhull/QuickHull.java
index c87cf06..1be8ced 100644
--- a/src/main/java/plugins/adufour/quickhull/QuickHull.java
+++ b/src/main/java/plugins/adufour/quickhull/QuickHull.java
@@ -1,3 +1,21 @@
+/*
+ * Copyright (c) 2010-2023. Institut Pasteur.
+ *
+ * This file is part of Icy.
+ * Icy is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Icy 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 Icy. If not, see <https://www.gnu.org/licenses/>.
+ */
+
package plugins.adufour.quickhull;
import icy.plugin.abstract_.Plugin;
@@ -5,14 +23,12 @@ import icy.plugin.interface_.PluginLibrary;
/**
* Main class of the QuickHull library for Icy
- *
+ * <p>
* {@link QuickHull2D}
* {@link QuickHull3D}
- *
- * @author Alexandre Dufour
*
+ * @author Alexandre Dufour
*/
-public class QuickHull extends Plugin implements PluginLibrary
-{
-
+public class QuickHull extends Plugin implements PluginLibrary {
+
}
diff --git a/src/main/java/plugins/adufour/quickhull/QuickHull2D.java b/src/main/java/plugins/adufour/quickhull/QuickHull2D.java
index 8d601d1..2f2542b 100644
--- a/src/main/java/plugins/adufour/quickhull/QuickHull2D.java
+++ b/src/main/java/plugins/adufour/quickhull/QuickHull2D.java
@@ -1,3 +1,21 @@
+/*
+ * Copyright (c) 2010-2023. Institut Pasteur.
+ *
+ * This file is part of Icy.
+ * Icy is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Icy 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 Icy. If not, see <https://www.gnu.org/licenses/>.
+ */
+
package plugins.adufour.quickhull;
import icy.plugin.abstract_.Plugin;
@@ -7,104 +25,91 @@ import java.awt.geom.Point2D;
import java.util.ArrayList;
import java.util.List;
-public class QuickHull2D extends Plugin
-{
- public static List<Point2D> computeConvexEnvelope(List<Point2D> points)
- {
- ArrayList<Point2D> envelope = new ArrayList<Point2D>();
-
+public class QuickHull2D extends Plugin {
+ public static List<Point2D> computeConvexEnvelope(final List<Point2D> points) {
+ final ArrayList<Point2D> envelope = new ArrayList<>();
+
// find two points: right (bottom) and left (top)
-
+
Point2D l = points.get(0);
Point2D r = points.get(0);
-
- for (int i = 1; i < points.size(); i++)
- {
- Point2D p = points.get(i);
+
+ for (int i = 1; i < points.size(); i++) {
+ final Point2D p = points.get(i);
if (r.getX() > p.getX() || (r.getX() == p.getX() && r.getY() > p.getY())) r = p;
if (l.getX() < p.getX() || (l.getX() == p.getX() && l.getY() < p.getY())) l = p;
}
-
+
// create neighbor lists
-
- ArrayList<Point2D> neighbors1 = new ArrayList<Point2D>();
- ArrayList<Point2D> neighbors2 = new ArrayList<Point2D>();
-
- for (Point2D p : points)
- {
+
+ final ArrayList<Point2D> neighbors1 = new ArrayList<>();
+ final ArrayList<Point2D> neighbors2 = new ArrayList<>();
+
+ for (final Point2D p : points) {
if (p == l || p == r) continue;
-
- int upper = Line2D.relativeCCW(r.getX(), r.getY(), l.getX(), l.getY(), p.getX(), p.getY());
-
- if (upper > 0)
- {
+
+ final int upper = Line2D.relativeCCW(r.getX(), r.getY(), l.getX(), l.getY(), p.getX(), p.getY());
+
+ if (upper > 0) {
neighbors1.add(p);
}
- else if (upper < 0)
- {
+ else if (upper < 0) {
neighbors2.add(p);
}
}
-
+
envelope.add(r);
-
+
quickhull(envelope, r, l, neighbors1);
-
+
envelope.add(l);
-
+
quickhull(envelope, l, r, neighbors2);
-
+
return envelope;
}
-
- private static void quickhull(ArrayList<Point2D> envelope, Point2D a, Point2D b, ArrayList<Point2D> neighbors)
- {
+
+ private static void quickhull(final ArrayList<Point2D> envelope, final Point2D a, final Point2D b, final ArrayList<Point2D> neighbors) {
if (neighbors.size() == 0) return;
-
- Point2D c = farthestpoint(a, b, neighbors);
-
- ArrayList<Point2D> al1 = new ArrayList<Point2D>();
- ArrayList<Point2D> al2 = new ArrayList<Point2D>();
-
- for (Point2D pt : neighbors)
- {
+
+ final Point2D c = farthestpoint(a, b, neighbors);
+
+ final ArrayList<Point2D> al1 = new ArrayList<>();
+ final ArrayList<Point2D> al2 = new ArrayList<>();
+
+ for (final Point2D pt : neighbors) {
if (pt == a || pt == b) continue;
-
- if (Line2D.relativeCCW(a.getX(), a.getY(), c.getX(), c.getY(), pt.getX(), pt.getY()) > 0)
- {
+
+ if (Line2D.relativeCCW(a.getX(), a.getY(), c.getX(), c.getY(), pt.getX(), pt.getY()) > 0) {
al1.add(pt);
}
- else if (Line2D.relativeCCW(c.getX(), c.getY(), b.getX(), b.getY(), pt.getX(), pt.getY()) > 0)
- {
+ else if (Line2D.relativeCCW(c.getX(), c.getY(), b.getX(), b.getY(), pt.getX(), pt.getY()) > 0) {
al2.add(pt);
}
}
-
+
quickhull(envelope, a, c, al1);
-
+
envelope.add(c);
-
+
quickhull(envelope, c, b, al2);
}
-
- private static Point2D farthestpoint(Point2D a, Point2D b, ArrayList<Point2D> points)
- {
+
+ private static Point2D farthestpoint(final Point2D a, final Point2D b, final ArrayList<Point2D> points) {
double maxD = -1;
Point2D maxP = null;
-
- for (Point2D p : points)
- {
+
+ for (final Point2D p : points) {
if (p == a || p == b) continue;
-
- double dist = Line2D.ptLineDistSq(a.getX(), a.getY(), b.getX(), b.getY(), p.getX(), p.getY());
-
- if (dist > maxD)
- {
+
+ final double dist = Line2D.ptLineDistSq(a.getX(), a.getY(), b.getX(), b.getY(), p.getX(), p.getY());
+
+ if (dist > maxD) {
maxD = dist;
maxP = p;
}
}
-
+
return maxP;
}
}
diff --git a/src/main/java/plugins/adufour/quickhull/QuickHull3D.java b/src/main/java/plugins/adufour/quickhull/QuickHull3D.java
index 7821d00..dac0a2f 100644
--- a/src/main/java/plugins/adufour/quickhull/QuickHull3D.java
+++ b/src/main/java/plugins/adufour/quickhull/QuickHull3D.java
@@ -1,40 +1,46 @@
-/**
- * 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.
+/*
+ * Copyright (c) 2010-2023. Institut Pasteur.
+ *
+ * This file is part of Icy.
+ * Icy is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
*
- * 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.
+ * Icy 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 Icy. If not, see <https://www.gnu.org/licenses/>.
*/
package plugins.adufour.quickhull;
-import java.util.*;
-import java.io.*;
-
import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;
+import java.io.*;
+import java.util.Arrays;
+import java.util.Iterator;
+import java.util.Vector;
/**
- * Computes the convex hull of a set of three dimensional points.
- *
+ * Computes the convex hull of a three-dimensional points set.
+ *
* <p>
- * The algorithm is a three dimensional implementation of Quickhull, as described in Barber, Dobkin,
+ * 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),
@@ -53,17 +59,17 @@ import javax.vecmath.Vector3d;
* 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>
- *
+ * <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.
- *
+ *
* <h3><a href=distTol>Robustness </a></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
@@ -74,1451 +80,1219 @@ import javax.vecmath.Vector3d;
* 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
+ * 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
+ * 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 Vector<Face> faces = new Vector<Face>(16);
- protected Vector<HalfEdge> horizon = new Vector<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 (Iterator<Face> it = faces.iterator(); it.hasNext();)
- {
- HalfEdge he = it.next().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();
- Vector<Integer> indexList = new Vector<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, Integer.valueOf((int) stok.nval));
- }
- faceIndices[i] = new int[indexList.size()];
- int k = 0;
- for (Integer integer : indexList)
- {
- faceIndices[i][k++] = integer;
- }
- }
- setHull(coords, nump, faceIndices, numf);
- }
- catch (Exception e)
- {
- e.printStackTrace();
- System.exit(1);
- }
- }
-
- /**
- * 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++)
- {
- Vector3d 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;
-
- {
- double diff;
-
- diff = maxVtxs[0].pnt.x - minVtxs[0].pnt.x;
- if (diff > max)
- {
- max = diff;
- imax = 0;
- }
-
- diff = maxVtxs[1].pnt.y - minVtxs[1].pnt.y;
- if (diff > max)
- {
- max = diff;
- imax = 0;
- }
-
- diff = maxVtxs[1].pnt.z - minVtxs[1].pnt.z;
- if (diff > max)
- {
- max = diff;
- imax = 0;
- }
- }
-
- 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.lengthSquared();
- 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 coplanar");
- }
- 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] = new Point3d(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 = new Point3d(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++)
- {
- Vector3d 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;
-
- 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 (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(Vector3d eyePnt, HalfEdge edge0, Face face, Vector<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, Vector<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<Face> 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++)
- {
- Vector3d 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;
- }
+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 final Face[] discardedFaces = new Face[3];
+
+ private final Vertex[] maxVtxs = new Vertex[3];
+ private final Vertex[] minVtxs = new Vertex[3];
+
+ protected Vector<Face> faces = new Vector<>(16);
+ protected Vector<HalfEdge> horizon = new Vector<>(16);
+
+ private final FaceList newFaces = new FaceList();
+ private final VertexList unclaimed = new VertexList();
+ private final 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 if 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(final 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(final double tol) {
+ explicitTolerance = tol;
+ }
+
+ /**
+ * Returns the explicit distance tolerance.
+ *
+ * @return explicit tolerance
+ * @see #setExplicitDistanceTolerance
+ */
+ public double getExplicitDistanceTolerance() {
+ return explicitTolerance;
+ }
+
+ private void addPointToFace(final Vertex vtx, final Face face) {
+ vtx.face = face;
+
+ if (face.outside == null) {
+ claimed.add(vtx);
+ }
+ else {
+ claimed.insertBefore(vtx, face.outside);
+ }
+ face.outside = vtx;
+ }
+
+ private void removePointFromFace(final Vertex vtx, final 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(final 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 point set 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 QuickHull3D(final double[] coords) throws IllegalArgumentException {
+ build(coords, coords.length / 3);
+ }
+
+ /**
+ * Creates a convex hull object and initializes it to the convex hull of a point set.
+ *
+ * @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(final Point3d[] points) throws IllegalArgumentException {
+ build(points, points.length);
+ }
+
+ private HalfEdge findHalfEdge(final Vertex tail, final Vertex head) {
+ // brute force ... OK, since setHull is not used much
+ for (final Face face : faces) {
+ final HalfEdge he = face.findEdge(tail, head);
+ if (he != null) {
+ return he;
+ }
+ }
+ return null;
+ }
+
+ protected void setHull(final double[] coords, final int nump, final int[][] faceIndices, final int numf) {
+ initBuffers(nump);
+ setPoints(coords, nump);
+ computeMaxAndMin();
+ for (int i = 0; i < numf; i++) {
+ final Face face = Face.create(pointBuffer, faceIndices[i]);
+ HalfEdge he = face.he0;
+ do {
+ final 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(final Process proc) throws IOException {
+ boolean wrote = false;
+ final InputStream es = proc.getErrorStream();
+ while (es.available() > 0) {
+ System.out.write(es.read());
+ wrote = true;
+ }
+ if (wrote) {
+ System.out.println();
+ }
+ }
+
+ protected void setFromQhull(final double[] coords, final int nump, final boolean triangulate) {
+ String commandStr = "./qhull i";
+ if (triangulate) {
+ commandStr += " -Qt";
+ }
+ try {
+ final Process proc = Runtime.getRuntime().exec(commandStr);
+ final PrintStream ps = new PrintStream(proc.getOutputStream());
+ final StreamTokenizer stok = new StreamTokenizer(new InputStreamReader(proc.getInputStream()));
+
+ ps.println("3 " + nump);
+ for (int i = 0; i < nump; i++) {
+ ps.println(coords[i * 3] + " " + coords[i * 3 + 1] + " " + coords[i * 3 + 2]);
+ }
+ ps.flush();
+ ps.close();
+ final Vector<Integer> indexList = new Vector<>(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);
+ }
+ final int numf = (int) stok.nval;
+ stok.nextToken(); // clear EOL
+ final 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, (int) stok.nval);
+ }
+ faceIndices[i] = new int[indexList.size()];
+ int k = 0;
+ for (final Integer integer : indexList) {
+ faceIndices[i][k++] = integer;
+ }
+ }
+ setHull(coords, nump, faceIndices, numf);
+ }
+ catch (final Exception e) {
+ e.printStackTrace();
+ System.exit(1);
+ }
+ }
+
+ /**
+ * Constructs the convex hull of a point set 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(final double[] coords) throws IllegalArgumentException {
+ build(coords, coords.length / 3);
+ }
+
+ /**
+ * Constructs the convex hull of a point set 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(final double[] coords, final 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 point set.
+ *
+ * @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(final Point3d[] points) throws IllegalArgumentException {
+ build(points, points.length);
+ }
+
+ /**
+ * Constructs the convex hull of a point set.
+ *
+ * @param points input points
+ * @param nump number of input points
+ * @throws IllegalArgumentException the number of input points is less than four or greater than the length of
+ * <code>points</code>, or the points appear to be coincident, colinear, or
+ * coplanar.
+ */
+ public void build(final Point3d[] points, final 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() {
+ final double minArea = 1000 * charLength * DOUBLE_PREC;
+ newFaces.clear();
+ for (final 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(final int nump) {
+ if (pointBuffer.length < nump) {
+ final 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;
+ }
+
+ protected void setPoints(final double[] coords, final int nump) {
+ for (int i = 0; i < nump; i++) {
+ final Vertex vtx = pointBuffer[i];
+ vtx.pnt.set(coords[i * 3], coords[i * 3 + 1], coords[i * 3 + 2]);
+ vtx.index = i;
+ }
+ }
+
+ protected void setPoints(final Point3d[] pnts, final int nump) {
+ for (int i = 0; i < nump; i++) {
+ final Vertex vtx = pointBuffer[i];
+ vtx.pnt.set(pnts[i]);
+ vtx.index = i;
+ }
+ }
+
+ protected void computeMaxAndMin() {
+ final Vector3d max = new Vector3d();
+ final 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++) {
+ final Vector3d 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;
+ final int imax = 0;
+
+ {
+ double diff;
+
+ diff = maxVtxs[0].pnt.x - minVtxs[0].pnt.x;
+ if (diff > max) {
+ max = diff;
+ //imax = 0;
+ }
+
+ diff = maxVtxs[1].pnt.y - minVtxs[1].pnt.y;
+ if (diff > max) {
+ max = diff;
+ //imax = 0;
+ }
+
+ diff = maxVtxs[1].pnt.z - minVtxs[1].pnt.z;
+ if (diff > max) {
+ max = diff;
+ //imax = 0;
+ }
+ }
+
+ if (max <= tolerance) {
+ throw new IllegalArgumentException("Input points appear to be coincident");
+ }
+ final 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
+ final Vector3d u01 = new Vector3d();
+ final Vector3d diff02 = new Vector3d();
+ final Vector3d nrml = new Vector3d();
+ final 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);
+ final double lenSqr = xprod.lengthSquared();
+ 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 coplanar");
+ }
+ nrml.normalize();
+
+ double maxDist = 0;
+ final double d0 = vtx[2].pnt.dot(nrml);
+ for (int i = 0; i < numPoints; i++) {
+ final 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);
+ }
+
+ final 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++) {
+ final 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++) {
+ final 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++) {
+ final 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++) {
+ final 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() {
+ final Point3d[] vtxs = new Point3d[numVertices];
+ for (int i = 0; i < numVertices; i++) {
+ vtxs[i] = new Point3d(pointBuffer[vertexPointIndices[i]].pnt);
+ }
+ return vtxs;
+ }
+
+ /**
+ * Returns the coordinates of this hull vertex points.
+ *
+ * @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(final double[] coords) {
+ for (int i = 0; i < numVertices; i++) {
+ final Point3d pnt = new Point3d(pointBuffer[vertexPointIndices[i]].pnt);
+ coords[i * 3] = 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() {
+ final 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(final int indexFlags) {
+ final int[][] allFaces = new int[faces.size()][];
+ int k = 0;
+ for (final 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 the 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(final 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 the 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(final PrintStream ps, int indexFlags) {
+ if ((indexFlags & INDEXED_FROM_ZERO) == 0) {
+ indexFlags |= INDEXED_FROM_ONE;
+ }
+ for (int i = 0; i < numVertices; i++) {
+ final Vector3d pnt = pointBuffer[vertexPointIndices[i]].pnt;
+ ps.println("v " + pnt.x + " " + pnt.y + " " + pnt.z);
+ }
+ for (final Face face : faces) {
+ final int[] indices = new int[face.numVertices()];
+ getFaceIndices(indices, face, indexFlags);
+
+ ps.print("f");
+ for (final int index : indices) {
+ ps.print(" " + index);
+ }
+ ps.println();
+ }
+ }
+
+ private void getFaceIndices(final int[] indices, final Face face, final int flags) {
+ final boolean ccw = ((flags & CLOCKWISE) == 0);
+ final boolean indexedFromOne = ((flags & INDEXED_FROM_ONE) != 0);
+ final 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(final 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) {
+ final 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(final Face face, final Face absorbingFace) {
+ final 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;
+ final 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(final HalfEdge he) {
+ return he.face.distanceToPlane(he.opposite.face.getCentroid());
+ }
+
+ private boolean doAdjacentMerge(final Face face, final int mergeType) {
+ HalfEdge hedge = face.he0;
+
+ boolean convex = true;
+ do {
+ final 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 (debug) {
+ System.out.println(" merging " + face.getVertexString() + " and " + oppFace.getVertexString());
+ }
+
+ final 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(final Vector3d eyePnt, HalfEdge edge0, final Face face, final Vector<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 {
+ final 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(final Vertex eyeVtx, final HalfEdge he) {
+ final 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(final FaceList newFaces, final Vertex eyeVtx, final Vector<HalfEdge> horizon) {
+ newFaces.clear();
+
+ HalfEdge hedgeSidePrev = null;
+ HalfEdge hedgeSideBegin = null;
+
+ for (final HalfEdge horizonHe : horizon) {
+ final 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;
+ }
+ assert hedgeSideBegin != null;
+ hedgeSideBegin.next.setOpposite(hedgeSidePrev);
+ }
+
+ protected Vertex nextPointToAdd() {
+ if (!claimed.isEmpty()) {
+ final Face eyeFace = claimed.first().face;
+ Vertex eyeVtx = null;
+ double maxDist = 0;
+ for (Vertex vtx = eyeFace.outside; vtx != null && vtx.face == eyeFace; vtx = vtx.next) {
+ final double dist = eyeFace.distanceToPlane(vtx.pnt);
+ if (dist > maxDist) {
+ maxDist = dist;
+ eyeVtx = vtx;
+ }
+ }
+ return eyeVtx;
+ }
+ else {
+ return null;
+ }
+ }
+
+ protected void addPointToHull(final 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(final Face face, final int mark) {
+ final 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 (final Iterator<Face> it = faces.iterator(); it.hasNext(); ) {
+ final 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++) {
+ final Vertex vtx = pointBuffer[i];
+ if (vtx.index == 0) {
+ vertexPointIndices[numVertices] = i;
+ vtx.index = numVertices++;
+ }
+ }
+ }
+
+ protected boolean checkFaceConvexity(final Face face, final double tol, final 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(final double tol, final PrintStream ps) {
+ // check edge convexity
+ boolean convex = true;
+ for (final 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(final 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(final PrintStream ps, final double tol) {
+ // check to make sure all edges are fully connected
+ // and that the edges are convex
+ double dist;
+ final double pointTol = 10 * tol;
+
+ if (!checkFaces(tolerance, ps)) {
+ return false;
+ }
+
+ // check point inclusion
+
+ for (int i = 0; i < numPoints; i++) {
+ final Vector3d pnt = pointBuffer[i].pnt;
+ for (final 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;
+ }
}
diff --git a/src/main/java/plugins/adufour/quickhull/Vertex.java b/src/main/java/plugins/adufour/quickhull/Vertex.java
index e70f2ce..b652810 100644
--- a/src/main/java/plugins/adufour/quickhull/Vertex.java
+++ b/src/main/java/plugins/adufour/quickhull/Vertex.java
@@ -1,3 +1,21 @@
+/*
+ * Copyright (c) 2010-2023. Institut Pasteur.
+ *
+ * This file is part of Icy.
+ * Icy is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Icy 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 Icy. If not, see <https://www.gnu.org/licenses/>.
+ */
+
package plugins.adufour.quickhull;
import javax.vecmath.Vector3d;
@@ -8,48 +26,46 @@ import javax.vecmath.Vector3d;
*
* @author John E. Lloyd, Fall 2004
*/
-class Vertex
-{
- /**
- * Spatial point associated with this vertex.
- */
- Vector3d pnt;
-
- /**
- * Back index into an array.
- */
- int index;
-
- /**
- * List forward link.
- */
- Vertex prev;
-
- /**
- * List backward link.
- */
- Vertex next;
-
- /**
- * Current face that this vertex is outside of.
- */
- Face face;
-
- /**
- * Constructs a vertex and sets its coordinates to 0.
- */
- public Vertex()
- { pnt = new Vector3d();
- }
-
- /**
- * Constructs a vertex with the specified coordinates
- * and index.
- */
- public Vertex (double x, double y, double z, int idx)
- {
- pnt = new Vector3d(x, y, z);
- index = idx;
- }
+class Vertex {
+ /**
+ * Spatial point associated with this vertex.
+ */
+ Vector3d pnt;
+
+ /**
+ * Back index into an array.
+ */
+ int index;
+
+ /**
+ * List forward link.
+ */
+ Vertex prev;
+
+ /**
+ * List backward link.
+ */
+ Vertex next;
+
+ /**
+ * Current face that this vertex is outside.
+ */
+ Face face;
+
+ /**
+ * Constructs a vertex and sets its coordinates to 0.
+ */
+ public Vertex() {
+ pnt = new Vector3d();
+ }
+
+ /**
+ * Constructs a vertex with the specified coordinates
+ * and index.
+ */
+ public Vertex(final double x, final double y, final double z, final int idx) {
+ pnt = new Vector3d(x, y, z);
+ index = idx;
+ }
}
diff --git a/src/main/java/plugins/adufour/quickhull/VertexList.java b/src/main/java/plugins/adufour/quickhull/VertexList.java
index fc6c29f..f6d6247 100644
--- a/src/main/java/plugins/adufour/quickhull/VertexList.java
+++ b/src/main/java/plugins/adufour/quickhull/VertexList.java
@@ -1,137 +1,131 @@
+/*
+ * Copyright (c) 2010-2023. Institut Pasteur.
+ *
+ * This file is part of Icy.
+ * Icy is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Icy 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 Icy. If not, see <https://www.gnu.org/licenses/>.
+ */
+
package plugins.adufour.quickhull;
/**
* Maintains a double-linked list of vertices for use by QuickHull3D
*/
-class VertexList
-{
- private Vertex head;
- private Vertex tail;
-
- /**
- * Clears this list.
- */
- public void clear()
- {
- head = tail = null;
- }
-
- /**
- * Adds a vertex to the end of this list.
- */
- public void add(Vertex vtx)
- {
- if (head == null)
- {
- head = vtx;
- }
- else
- {
- tail.next = vtx;
- }
- vtx.prev = tail;
- vtx.next = null;
- tail = vtx;
- }
-
- /**
- * Adds a chain of vertices to the end of this list.
- */
- public void addAll(Vertex vtx)
- {
- if (head == null)
- {
- head = vtx;
- }
- else
- {
- tail.next = vtx;
- }
- vtx.prev = tail;
- while (vtx.next != null)
- {
- vtx = vtx.next;
- }
- tail = vtx;
- }
-
- /**
- * Deletes a vertex from this list.
- */
- public void delete(Vertex vtx)
- {
- if (vtx.prev == null)
- {
- head = vtx.next;
- }
- else
- {
- vtx.prev.next = vtx.next;
- }
- if (vtx.next == null)
- {
- tail = vtx.prev;
- }
- else
- {
- vtx.next.prev = vtx.prev;
- }
- }
-
- /**
- * Deletes a chain of vertices from this list.
- */
- public void delete(Vertex vtx1, Vertex vtx2)
- {
- if (vtx1.prev == null)
- {
- head = vtx2.next;
- }
- else
- {
- vtx1.prev.next = vtx2.next;
- }
- if (vtx2.next == null)
- {
- tail = vtx1.prev;
- }
- else
- {
- vtx2.next.prev = vtx1.prev;
- }
- }
-
- /**
- * Inserts a vertex into this list before another specificed vertex.
- */
- public void insertBefore(Vertex vtx, Vertex next)
- {
- vtx.prev = next.prev;
- if (next.prev == null)
- {
- head = vtx;
- }
- else
- {
- next.prev.next = vtx;
- }
- vtx.next = next;
- next.prev = vtx;
- }
-
- /**
- * Returns the first element in this list.
- */
- public Vertex first()
- {
- return head;
- }
-
- /**
- * Returns true if this list is empty.
- */
- public boolean isEmpty()
- {
- return head == null;
- }
+class VertexList {
+ private Vertex head;
+ private Vertex tail;
+
+ /**
+ * Clears this list.
+ */
+ public void clear() {
+ head = tail = null;
+ }
+
+ /**
+ * Adds a vertex to the end of this list.
+ */
+ public void add(final Vertex vtx) {
+ if (head == null) {
+ head = vtx;
+ }
+ else {
+ tail.next = vtx;
+ }
+ vtx.prev = tail;
+ vtx.next = null;
+ tail = vtx;
+ }
+
+ /**
+ * Adds a chain of vertices to the end of this list.
+ */
+ public void addAll(Vertex vtx) {
+ if (head == null) {
+ head = vtx;
+ }
+ else {
+ tail.next = vtx;
+ }
+ vtx.prev = tail;
+ while (vtx.next != null) {
+ vtx = vtx.next;
+ }
+ tail = vtx;
+ }
+
+ /**
+ * Deletes a vertex from this list.
+ */
+ public void delete(final Vertex vtx) {
+ if (vtx.prev == null) {
+ head = vtx.next;
+ }
+ else {
+ vtx.prev.next = vtx.next;
+ }
+ if (vtx.next == null) {
+ tail = vtx.prev;
+ }
+ else {
+ vtx.next.prev = vtx.prev;
+ }
+ }
+
+ /**
+ * Deletes a chain of vertices from this list.
+ */
+ public void delete(final Vertex vtx1, final Vertex vtx2) {
+ if (vtx1.prev == null) {
+ head = vtx2.next;
+ }
+ else {
+ vtx1.prev.next = vtx2.next;
+ }
+ if (vtx2.next == null) {
+ tail = vtx1.prev;
+ }
+ else {
+ vtx2.next.prev = vtx1.prev;
+ }
+ }
+
+ /**
+ * Inserts a vertex into this list before another specificed vertex.
+ */
+ public void insertBefore(final Vertex vtx, final Vertex next) {
+ vtx.prev = next.prev;
+ if (next.prev == null) {
+ head = vtx;
+ }
+ else {
+ next.prev.next = vtx;
+ }
+ vtx.next = next;
+ next.prev = vtx;
+ }
+
+ /**
+ * Returns the first element in this list.
+ */
+ public Vertex first() {
+ return head;
+ }
+
+ /**
+ * Returns true if this list is empty.
+ */
+ public boolean isEmpty() {
+ return head == null;
+ }
}
--
GitLab