View Javadoc

    1   /*
    2    * Licensed to the Apache Software Foundation (ASF) under one or more
    3    * contributor license agreements.  See the NOTICE file distributed with
    4    * this work for additional information regarding copyright ownership.
    5    * The ASF licenses this file to You under the Apache License, Version 2.0
    6    * (the "License"); you may not use this file except in compliance with
    7    * the License.  You may obtain a copy of the License at
    8    *
    9    *      http://www.apache.org/licenses/LICENSE-2.0
    10   *
    11   * Unless required by applicable law or agreed to in writing, software
    12   * distributed under the License is distributed on an "AS IS" BASIS,
    13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    14   * See the License for the specific language governing permissions and
    15   * limitations under the License.
    16   */
    17  package org.apache.commons.geometry.euclidean.threed;
    18
    19  import java.text.MessageFormat;
    20  import java.util.ArrayList;
    21  import java.util.Arrays;
    22  import java.util.Collection;
    23  import java.util.Iterator;
    24  import java.util.List;
    25  import java.util.function.BiFunction;
    26
    27  import org.apache.commons.geometry.core.partitioning.HyperplaneBoundedRegion;
    28  import org.apache.commons.geometry.core.partitioning.Split;
    29  import org.apache.commons.geometry.core.partitioning.SplitLocation;
    30  import org.apache.commons.geometry.core.precision.DoublePrecisionContext;
    31  import org.apache.commons.geometry.euclidean.threed.line.Line3D;
    32  import org.apache.commons.geometry.euclidean.threed.line.LineConvexSubset3D;
    33  import org.apache.commons.geometry.euclidean.twod.ConvexArea;
    34  import org.apache.commons.geometry.euclidean.twod.Line;
    35  import org.apache.commons.geometry.euclidean.twod.LineConvexSubset;
    36  import org.apache.commons.geometry.euclidean.twod.Lines;
    37  import org.apache.commons.geometry.euclidean.twod.RegionBSPTree2D;
    38  import org.apache.commons.geometry.euclidean.twod.Vector2D;
    39  import org.apache.commons.geometry.euclidean.twod.path.LinePath;
    40
    41  /** Class containing factory methods for constructing {@link Plane} and {@link PlaneSubset} instances.
    42   */
    43  public final class Planes {
    44
    45      /** Utility class; no instantiation. */
    46      private Planes() {
    47      }
    48
    49      /** Build a plane from a point and two (on plane) vectors.
    50       * @param p the provided point (on plane)
    51       * @param u u vector (on plane)
    52       * @param v v vector (on plane)
    53       * @param precision precision context used to compare floating point values
    54       * @return a new plane
    55       * @throws IllegalArgumentException if the norm of the given values is zero, NaN, or infinite.
    56       */
    57      public static EmbeddingPlane fromPointAndPlaneVectors(final Vector3D p, final Vector3D u, final Vector3D v,
    58              final DoublePrecisionContext precision) {
    59          final Vector3D.Unit uNorm = u.normalize();
    60          final Vector3D.Unit vNorm = uNorm.orthogonal(v);
    61          final Vector3D.Unit wNorm = uNorm.cross(vNorm).normalize();
    62          final double originOffset = -p.dot(wNorm);
    63
    64          return new EmbeddingPlane(uNorm, vNorm, wNorm, originOffset, precision);
    65      }
    66
    67      /** Build a plane from a normal.
    68       * Chooses origin as point on plane.
    69       * @param normal normal direction to the plane
    70       * @param precision precision context used to compare floating point values
    71       * @return a new plane
    72       * @throws IllegalArgumentException if the norm of the given values is zero, NaN, or infinite.
    73       */
    74      public static Plane fromNormal(final Vector3D normal, final DoublePrecisionContext precision) {
    75          return fromPointAndNormal(Vector3D.ZERO, normal, precision);
    76      }
    77
    78      /** Build a plane from a point and a normal.
    79       *
    80       * @param p point belonging to the plane
    81       * @param normal normal direction to the plane
    82       * @param precision precision context used to compare floating point values
    83       * @return a new plane
    84       * @throws IllegalArgumentException if the norm of the given values is zero, NaN, or infinite.
    85       */
    86      public static Plane fromPointAndNormal(final Vector3D p, final Vector3D normal,
    87              final DoublePrecisionContext precision) {
    88          final Vector3D.Unit unitNormal = normal.normalize();
    89          final double originOffset = -p.dot(unitNormal);
    90
    91          return new Plane(unitNormal, originOffset, precision);
    92      }
    93
    94      /** Build a plane from three points.
    95       * <p>
    96       * The plane is oriented in the direction of {@code (p2-p1) ^ (p3-p1)}
    97       * </p>
    98       *
    99       * @param p1 first point belonging to the plane
    100      * @param p2 second point belonging to the plane
    101      * @param p3 third point belonging to the plane
    102      * @param precision precision context used to compare floating point values
    103      * @return a new plane
    104      * @throws IllegalArgumentException if the points do not define a unique plane
    105      */
    106     public static Plane fromPoints(final Vector3D p1, final Vector3D p2, final Vector3D p3,
    107             final DoublePrecisionContext precision) {
    108         return fromPoints(Arrays.asList(p1, p2, p3), precision);
    109     }
    110
    111     /** Construct a plane from a collection of points lying on the plane. The plane orientation is
    112      * determined by the overall orientation of the point sequence. For example, if the points wind
    113      * around the z-axis in a counter-clockwise direction, then the plane normal will point up the
    114      * +z axis. If the points wind in the opposite direction, then the plane normal will point down
    115      * the -z axis. The {@code u} vector for the plane is set to the first non-zero vector between
    116      * points in the sequence (ie, the first direction in the path).
    117      *
    118      * @param pts collection of sequenced points lying on the plane
    119      * @param precision precision context used to compare floating point values
    120      * @return a new plane containing the given points
    121      * @throws IllegalArgumentException if the given collection does not contain at least 3 points or the
    122      *      points do not define a unique plane
    123      */
    124     public static Plane fromPoints(final Collection<Vector3D> pts, final DoublePrecisionContext precision) {
    125         return new PlaneBuilder(pts, precision).build();
    126     }
    127
    128     /** Create a new plane subset from a plane and an embedded convex subspace area.
    129      * @param plane embedding plane for the area
    130      * @param area area embedded in the plane
    131      * @return a new convex sub plane instance
    132      */
    133     public static PlaneConvexSubset subsetFromConvexArea(final EmbeddingPlane plane, final ConvexArea area) {
    134         if (area.isFinite()) {
    135             // prefer a vertex-based representation for finite areas
    136             final List<Vector3D> vertices = plane.toSpace(area.getVertices());
    137             return fromConvexPlanarVertices(plane, vertices);
    138         }
    139
    140         return new EmbeddedAreaPlaneConvexSubset(plane, area);
    141     }
    142
    143     /** Create a new convex polygon from the given sequence of vertices. The vertices must define a unique
    144      * plane, meaning that at least 3 unique vertices must be given. The given sequence is assumed to be closed,
    145      * ie that an edge exists between the last vertex and the first.
    146      * @param pts collection of points defining the convex polygon
    147      * @param precision precision context used to compare floating point values
    148      * @return a new convex polygon defined by the given sequence of vertices
    149      * @throws IllegalArgumentException if fewer than 3 vertices are given or the vertices do not define a
    150      *       unique plane
    151      * @see #fromPoints(Collection, DoublePrecisionContext)
    152      */
    153     public static ConvexPolygon3D convexPolygonFromVertices(final Collection<Vector3D> pts,
    154             final DoublePrecisionContext precision) {
    155         final List<Vector3D> vertices = new ArrayList<>(pts.size());
    156         final Plane plane = new PlaneBuilder(pts, precision).buildForConvexPolygon(vertices);
    157
    158         // make sure that the first point is not repeated at the end
    159         final Vector3D firstPt = vertices.get(0);
    160         final Vector3D lastPt = vertices.get(vertices.size() - 1);
    161         if (firstPt.eq(lastPt, precision)) {
    162             vertices.remove(vertices.size() - 1);
    163         }
    164
    165         if (vertices.size() == 3) {
    166             return new SimpleTriangle3D(plane, vertices.get(0), vertices.get(1), vertices.get(2));
    167         }
    168         return new VertexListConvexPolygon3D(plane, vertices);
    169     }
    170
    171     /** Construct a triangle from three vertices. The triangle plane is oriented such that the points
    172      * are arranged in a counter-clockwise order when looking down the plane normal.
    173      * @param p1 first vertex
    174      * @param p2 second vertex
    175      * @param p3 third vertex
    176      * @param precision precision context used for floating point comparisons
    177      * @return a triangle constructed from the three vertices
    178      * @throws IllegalArgumentException if the points do not define a unique plane
    179      */
    180     public static Triangle3D triangleFromVertices(final Vector3D p1, final Vector3D p2, final Vector3D p3,
    181             final DoublePrecisionContext precision) {
    182         final Plane plane = fromPoints(p1, p2, p3, precision);
    183         return new SimpleTriangle3D(plane, p1, p2, p3);
    184     }
    185
    186     /** Construct a list of {@link Triangle3D} instances from a set of vertices and arrays of face indices.
    187      * For example, the following code constructs a list of triangles forming a square pyramid.
    188      * <pre>
    189      * DoublePrecisionContext precision = new EpsilonDoublePrecisionContext(1e-10);
    190      *
    191      * Vector3D[] vertices = {
    192      *      Vector3D.ZERO,
    193      *      Vector3D.of(1, 0, 0),
    194      *      Vector3D.of(1, 1, 0),
    195      *      Vector3D.of(0, 1, 0),
    196      *      Vector3D.of(0.5, 0.5, 4)
    197      * };
    198      *
    199      * int[][] faceIndices = {
    200      *      {0, 2, 1},
    201      *      {0, 3, 2},
    202      *      {0, 1, 4},
    203      *      {1, 2, 4},
    204      *      {2, 3, 4},
    205      *      {3, 0, 4}
    206      * };
    207      *
    208      * List&lt;Triangle3D&gt; triangles = Planes.indexedTriangles(vertices, faceIndices, TEST_PRECISION);
    209      * </pre>
    210      * @param vertices vertices available for use in triangle construction
    211      * @param faceIndices array of indices for each triangular face; each entry in the array is an array of
    212      *      3 index values into {@code vertices}, defining the 3 vertices that will be used to construct the
    213      *      triangle
    214      * @param precision precision context used for floating point comparisons
    215      * @return a list of triangles constructed from the set of vertices and face indices
    216      * @throws IllegalArgumentException if any face index array does not contain exactly 3 elements or a set
    217      *      of 3 vertices do not define a plane
    218      * @throws IndexOutOfBoundsException if any index into {@code vertices} is out of bounds
    219      */
    220     public static List<Triangle3D> indexedTriangles(final Vector3D[] vertices, final int[][] faceIndices,
    221             final DoublePrecisionContext precision) {
    222         return indexedTriangles(Arrays.asList(vertices), faceIndices, precision);
    223     }
    224
    225     /** Construct a list of {@link Triangle3D} instances from a set of vertices and arrays of face indices.
    226      * @param vertices vertices available for use in triangle construction
    227      * @param faceIndices array of indices for each triangular face; each entry in the array is an array of
    228      *      3 index values into {@code vertices}, defining the 3 vertices that will be used to construct the
    229      *      triangle
    230      * @param precision precision context used for floating point comparisons
    231      * @return a list of triangles constructed from the set of vertices and face indices
    232      * @throws IllegalArgumentException if any face index array does not contain exactly 3 elements or a set
    233      *      of 3 vertices do not define a plane
    234      * @throws IndexOutOfBoundsException if any index into {@code vertices} is out of bounds
    235      * @see #indexedTriangles(Vector3D[], int[][], DoublePrecisionContext)
    236      */
    237     public static List<Triangle3D> indexedTriangles(final List<Vector3D> vertices, final int[][] faceIndices,
    238             final DoublePrecisionContext precision) {
    239
    240         final int numFaces = faceIndices.length;
    241         final List<Triangle3D> triangles = new ArrayList<>(numFaces);
    242
    243         int[] face;
    244         for (int i = 0; i < numFaces; ++i) {
    245             face = faceIndices[i];
    246             if (face.length != 3) {
    247                 throw new IllegalArgumentException(MessageFormat.format(
    248                         "Invalid number of vertex indices for face at index {0}: expected 3 but found {1}",
    249                         i, face.length));
    250             }
    251
    252             triangles.add(triangleFromVertices(
    253                         vertices.get(face[0]),
    254                         vertices.get(face[1]),
    255                         vertices.get(face[2]),
    256                         precision
    257                     ));
    258         }
    259
    260         return triangles;
    261     }
    262
    263     /** Construct a list of {@link ConvexPolygon3D} instances from a set of vertices and arrays of face indices. Each
    264      * face must contain at least 3 vertices but the number of vertices per face does not need to be constant.
    265      * For example, the following code constructs a list of convex polygons forming a square pyramid.
    266      * Note that the first face (the pyramid base) uses a different number of vertices than the other faces.
    267      * <pre>
    268      * DoublePrecisionContext precision = new EpsilonDoublePrecisionContext(1e-10);
    269      *
    270      * Vector3D[] vertices = {
    271      *      Vector3D.ZERO,
    272      *      Vector3D.of(1, 0, 0),
    273      *      Vector3D.of(1, 1, 0),
    274      *      Vector3D.of(0, 1, 0),
    275      *      Vector3D.of(0.5, 0.5, 4)
    276      * };
    277      *
    278      * int[][] faceIndices = {
    279      *      {0, 3, 2, 1}, // square base
    280      *      {0, 1, 4},
    281      *      {1, 2, 4},
    282      *      {2, 3, 4},
    283      *      {3, 0, 4}
    284      * };
    285      *
    286      * List&lt;ConvexPolygon3D&gt; polygons = Planes.indexedConvexPolygons(vertices, faceIndices, precision);
    287      * </pre>
    288      * @param vertices vertices available for use in convex polygon construction
    289      * @param faceIndices array of indices for each triangular face; each entry in the array is an array of
    290      *      at least 3 index values into {@code vertices}, defining the vertices that will be used to construct the
    291      *      convex polygon
    292      * @param precision precision context used for floating point comparisons
    293      * @return a list of convex polygons constructed from the set of vertices and face indices
    294      * @throws IllegalArgumentException if any face index array does not contain at least 3 elements or a set
    295      *      of vertices do not define a planar convex polygon
    296      * @throws IndexOutOfBoundsException if any index into {@code vertices} is out of bounds
    297      */
    298     public static List<ConvexPolygon3D> indexedConvexPolygons(final Vector3D[] vertices, final int[][] faceIndices,
    299             final DoublePrecisionContext precision) {
    300         return indexedConvexPolygons(Arrays.asList(vertices), faceIndices, precision);
    301     }
    302
    303     /** Construct a list of {@link ConvexPolygon3D} instances from a set of vertices and arrays of face indices. Each
    304      * face must contain at least 3 vertices but the number of vertices per face does not need to be constant.
    305      * @param vertices vertices available for use in convex polygon construction
    306      * @param faceIndices array of indices for each triangular face; each entry in the array is an array of
    307      *      at least 3 index values into {@code vertices}, defining the vertices that will be used to construct the
    308      *      convex polygon
    309      * @param precision precision context used for floating point comparisons
    310      * @return a list of convex polygons constructed from the set of vertices and face indices
    311      * @throws IllegalArgumentException if any face index array does not contain at least 3 elements or a set
    312      *      of vertices do not define a planar convex polygon
    313      * @throws IndexOutOfBoundsException if any index into {@code vertices} is out of bounds
    314      * @see #indexedConvexPolygons(Vector3D[], int[][], DoublePrecisionContext)
    315      */
    316     public static List<ConvexPolygon3D> indexedConvexPolygons(final List<Vector3D> vertices, final int[][] faceIndices,
    317             final DoublePrecisionContext precision) {
    318         final int numFaces = faceIndices.length;
    319         final List<ConvexPolygon3D> polygons = new ArrayList<>(numFaces);
    320         final List<Vector3D> faceVertices = new ArrayList<>();
    321
    322         int[] face;
    323         for (int i = 0; i < numFaces; ++i) {
    324             face = faceIndices[i];
    325             if (face.length < 3) {
    326                 throw new IllegalArgumentException(MessageFormat.format(
    327                         "Invalid number of vertex indices for face at index {0}: required at least 3 but found {1}",
    328                         i, face.length));
    329             }
    330
    331             for (final int vertexIndex : face) {
    332                 faceVertices.add(vertices.get(vertexIndex));
    333             }
    334
    335             polygons.add(convexPolygonFromVertices(
    336                         faceVertices,
    337                         precision
    338                     ));
    339
    340             faceVertices.clear();
    341         }
    342
    343         return polygons;
    344     }
    345
    346     /** Get the boundaries of a 3D region created by extruding a polygon defined by a list of vertices. The ends
    347      * ("top" and "bottom") of the extruded 3D region are flat while the sides follow the boundaries of the original
    348      * 2D region.
    349      * @param vertices vertices forming the 2D polygon to extrude
    350      * @param plane plane to extrude the 2D polygon from
    351      * @param extrusionVector vector to extrude the polygon vertices through
    352      * @param precision precision context used to construct the 3D region boundaries
    353      * @return the boundaries of the extruded 3D region
    354      * @throws IllegalStateException if {@code vertices} contains only a single unique vertex
    355      * @throws IllegalArgumentException if regions of non-zero size cannot be produced with the
    356      *      given plane and extrusion vector. This occurs when the extrusion vector has zero length
    357      *      or is orthogonal to the plane normal
    358      * @see LinePath#fromVertexLoop(Collection, DoublePrecisionContext)
    359      * @see #extrude(LinePath, EmbeddingPlane, Vector3D, DoublePrecisionContext)
    360      */
    361     public static List<PlaneConvexSubset> extrudeVertexLoop(final List<Vector2D> vertices,
    362             final EmbeddingPlane plane, final Vector3D extrusionVector, final DoublePrecisionContext precision) {
    363         final LinePath path = LinePath.fromVertexLoop(vertices, precision);
    364         return extrude(path, plane, extrusionVector, precision);
    365     }
    366
    367     /** Get the boundaries of the 3D region created by extruding a 2D line path. The ends ("top" and "bottom") of
    368      * the extruded 3D region are flat while the sides follow the boundaries of the original 2D region. The path is
    369      * converted to a BSP tree before extrusion.
    370      * @param path path to extrude
    371      * @param plane plane to extrude the path from
    372      * @param extrusionVector vector to extrude the polygon points through
    373      * @param precision precision precision context used to construct the 3D region boundaries
    374      * @return the boundaries of the extruded 3D region
    375      * @throws IllegalArgumentException if regions of non-zero size cannot be produced with the
    376      *      given plane and extrusion vector. This occurs when the extrusion vector has zero length
    377      *      or is orthogonal to the plane normal
    378      * @see #extrude(RegionBSPTree2D, EmbeddingPlane, Vector3D, DoublePrecisionContext)
    379      */
    380     public static List<PlaneConvexSubset> extrude(final LinePath path, final EmbeddingPlane plane,
    381             final Vector3D extrusionVector, final DoublePrecisionContext precision) {
    382         return extrude(path.toTree(), plane, extrusionVector, precision);
    383     }
    384
    385     /** Get the boundaries of the 3D region created by extruding a 2D region. The ends ("top" and "bottom") of
    386      * the extruded 3D region are flat while the sides follow the boundaries of the original 2D region.
    387      * @param region region to extrude
    388      * @param plane plane to extrude the region from
    389      * @param extrusionVector vector to extrude the region points through
    390      * @param precision precision precision context used to construct the 3D region boundaries
    391      * @return the boundaries of the extruded 3D region
    392      * @throws IllegalArgumentException if regions of non-zero size cannot be produced with the
    393      *      given plane and extrusion vector. This occurs when the extrusion vector has zero length
    394      *      or is orthogonal to the plane normal
    395      */
    396     public static List<PlaneConvexSubset> extrude(final RegionBSPTree2D region, final EmbeddingPlane plane,
    397             final Vector3D extrusionVector, final DoublePrecisionContext precision) {
    398         return new PlaneRegionExtruder(plane, extrusionVector, precision).extrude(region);
    399     }
    400
    401     /** Get the unique intersection of the plane subset with the given line. Null is
    402      * returned if no unique intersection point exists (ie, the line and plane are
    403      * parallel or coincident) or the line does not intersect the plane subset.
    404      * @param planeSubset plane subset to intersect with
    405      * @param line line to intersect with this plane subset
    406      * @return the unique intersection point between the line and this plane subset
    407      *      or null if no such point exists.
    408      */
    409     static Vector3D intersection(final PlaneSubset planeSubset, final Line3D line) {
    410         final Vector3D pt = planeSubset.getPlane().intersection(line);
    411         return (pt != null && planeSubset.contains(pt)) ? pt : null;
    412     }
    413
    414     /** Get the unique intersection of the plane subset with the given line subset. Null
    415      * is returned if the underlying line and plane do not have a unique intersection
    416      * point (ie, they are parallel or coincident) or the intersection point is unique
    417      * but is not contained in both the line subset and plane subset.
    418      * @param planeSubset plane subset to intersect with
    419      * @param lineSubset line subset to intersect with
    420      * @return the unique intersection point between this plane subset and the argument or
    421      *      null if no such point exists.
    422      */
    423     static Vector3D intersection(final PlaneSubset planeSubset, final LineConvexSubset3D lineSubset) {
    424         final Vector3D pt = intersection(planeSubset, lineSubset.getLine());
    425         return (pt != null && lineSubset.contains(pt)) ? pt : null;
    426     }
    427
    428     /** Validate that the actual plane contains the same points as the expected plane, throwing an exception if not.
    429      * The subspace orientations of embedding planes are not considered.
    430      * @param expected the expected plane
    431      * @param actual the actual plane
    432      * @throws IllegalArgumentException if the actual plane is not equivalent to the expected plane
    433      */
    434     static void validatePlanesEquivalent(final Plane expected, final Plane actual) {
    435         if (!expected.eq(actual, expected.getPrecision())) {
    436             throw new IllegalArgumentException("Arguments do not represent the same plane. Expected " +
    437                     expected + " but was " + actual + ".");
    438         }
    439     }
    440
    441     /** Generic split method that uses performs the split using the subspace region of the plane subset.
    442      * @param splitter splitting hyperplane
    443      * @param subset the plane subset being split
    444      * @param factory function used to create new plane subset instances
    445      * @param <T> Plane subset implementation type
    446      * @return the result of the split operation
    447      */
    448     static <T extends PlaneSubset> Split<T> subspaceSplit(final Plane splitter, final T subset,
    449             final BiFunction<EmbeddingPlane, HyperplaneBoundedRegion<Vector2D>, T> factory) {
    450
    451         final EmbeddingPlane thisPlane = subset.getPlane().getEmbedding();
    452
    453         final Line3D intersection = thisPlane.intersection(splitter);
    454         if (intersection == null) {
    455             return getNonIntersectingSplitResult(splitter, subset);
    456         } else {
    457             final EmbeddingPlane embeddingPlane = subset.getPlane().getEmbedding();
    458
    459             // the lines intersect; split the subregion
    460             final Vector3D intersectionOrigin = intersection.getOrigin();
    461             final Vector2D subspaceP1 = embeddingPlane.toSubspace(intersectionOrigin);
    462             final Vector2D subspaceP2 = embeddingPlane.toSubspace(intersectionOrigin.add(intersection.getDirection()));
    463
    464             final Line subspaceSplitter = Lines.fromPoints(subspaceP1, subspaceP2, thisPlane.getPrecision());
    465
    466             final Split<? extends HyperplaneBoundedRegion<Vector2D>> split =
    467                     subset.getEmbedded().getSubspaceRegion().split(subspaceSplitter);
    468             final SplitLocation subspaceSplitLoc = split.getLocation();
    469
    470             if (SplitLocation.MINUS == subspaceSplitLoc) {
    471                 return new Split<>(subset, null);
    472             } else if (SplitLocation.PLUS == subspaceSplitLoc) {
    473                 return new Split<>(null, subset);
    474             }
    475
    476             final T minus = (split.getMinus() != null) ? factory.apply(thisPlane, split.getMinus()) : null;
    477             final T plus = (split.getPlus() != null) ? factory.apply(thisPlane, split.getPlus()) : null;
    478
    479             return new Split<>(minus, plus);
    480         }
    481     }
    482
    483     /** Get a split result for cases where the splitting plane and the plane containing the subset being split
    484      * do not intersect. Callers are responsible for ensuring that the planes involved do not actually intersect.
    485      * @param <T> Plane subset implementation type
    486      * @param splitter plane performing the splitting
    487      * @param subset subset being split
    488      * @return the split result for the non-intersecting split
    489      */
    490     private static <T extends PlaneSubset> Split<T> getNonIntersectingSplitResult(
    491             final Plane splitter, final T subset) {
    492         final Plane plane = subset.getPlane();
    493
    494         final double offset = splitter.offset(plane);
    495         final int comp = plane.getPrecision().compare(offset, 0.0);
    496
    497         if (comp < 0) {
    498             return new Split<>(subset, null);
    499         } else if (comp > 0) {
    500             return new Split<>(null, subset);
    501         } else {
    502             return new Split<>(null, null);
    503         }
    504     }
    505
    506     /** Construct a convex polygon 3D from a plane and a list of vertices lying in the plane. Callers are
    507      * responsible for ensuring that the vertices lie in the plane and define a convex polygon.
    508      * @param plane the plane containing the convex polygon
    509      * @param vertices vertices defining the closed, convex polygon. The must must contain at least 3 unique
    510      *      vertices and should not include the start vertex at the end of the list.
    511      * @return a new convex polygon instance
    512      * @throws IllegalArgumentException if the size of {@code vertices} if less than 3
    513      */
    514     static ConvexPolygon3D fromConvexPlanarVertices(final Plane plane, final List<Vector3D> vertices) {
    515         final int size = vertices.size();
    516
    517         if (size == 3) {
    518             return new SimpleTriangle3D(plane, vertices.get(0), vertices.get(1), vertices.get(2));
    519         }
    520
    521         return new VertexListConvexPolygon3D(plane, vertices);
    522     }
    523
    524     /** Convert a convex polygon defined by a plane and list of points into a triangle fan.
    525      * @param plane plane containing the convex polygon
    526      * @param vertices vertices defining the convex polygon
    527      * @return a triangle fan representing the same area as the convex polygon
    528      * @throws IllegalArgumentException if fewer than 3 vertices are given
    529      */
    530     static List<Triangle3D> convexPolygonToTriangleFan(final Plane plane, final List<Vector3D> vertices) {
    531         final int size = vertices.size();
    532         if (size < 3) {
    533             throw new IllegalArgumentException("Cannot create triangle fan: 3 or more vertices are required " +
    534                     "but found only " + vertices.size());
    535         }
    536
    537         final List<Triangle3D> triangles = new ArrayList<>(size - 2);
    538
    539         final int fanIdx = findBestTriangleFanIndex(vertices);
    540         int vertexIdx = (fanIdx + 1) % size;
    541
    542         final Vector3D fanBase = vertices.get(fanIdx);
    543         Vector3D vertexA = vertices.get(vertexIdx);
    544         Vector3D vertexB;
    545
    546         vertexIdx = (vertexIdx + 1) % size;
    547         while (vertexIdx != fanIdx) {
    548             vertexB = vertices.get(vertexIdx);
    549
    550             // add directly as a triangle instance to avoid computation of the plane again
    551             triangles.add(new SimpleTriangle3D(plane, fanBase, vertexA, vertexB));
    552
    553             vertexA = vertexB;
    554             vertexIdx = (vertexIdx + 1) % size;
    555         }
    556
    557         return triangles;
    558     }
    559
    560     /** Find the index of the best vertex to use as the base for a triangle fan split of the convex polygon
    561      * defined by the given vertices. The best vertex is the one that forms the largest interior angle in the
    562      * polygon since a split at that point will help prevent the creation of very thin triangles.
    563      * @param vertices vertices defining the convex polygon; must not be empty
    564      * @return the index of the best vertex to use as the base for a triangle fan split of the convex polygon
    565      */
    566     private static int findBestTriangleFanIndex(final List<Vector3D> vertices) {
    567         final Iterator<Vector3D> it = vertices.iterator();
    568
    569         Vector3D curPt = it.next();
    570         Vector3D nextPt;
    571
    572         final Vector3D lastVec = vertices.get(vertices.size() - 1).directionTo(curPt);
    573         Vector3D incomingVec = lastVec;
    574         Vector3D outgoingVec;
    575
    576         int bestIdx = 0;
    577         double bestDot = -1.0;
    578
    579         int idx = 0;
    580         double dot;
    581         while (it.hasNext()) {
    582             nextPt = it.next();
    583             outgoingVec = curPt.directionTo(nextPt);
    584
    585             dot = incomingVec.dot(outgoingVec);
    586             if (dot > bestDot) {
    587                 bestIdx = idx;
    588                 bestDot = dot;
    589             }
    590
    591             curPt = nextPt;
    592             incomingVec = outgoingVec;
    593
    594             ++idx;
    595         }
    596
    597         // handle the last vertex on its own
    598         dot = incomingVec.dot(lastVec);
    599         if (dot > bestDot) {
    600             bestIdx = idx;
    601         }
    602
    603         return bestIdx;
    604     }
    605
    606     /** Internal helper class used to construct planes from sequences of points. Instances can be also be
    607      * configured to collect lists of unique points found during plane construction and validate that the
    608      * defined region is convex.
    609      */
    610     private static final class PlaneBuilder {
    611
    612         /** The point sequence to build a plane for. */
    613         private final Collection<Vector3D> pts;
    614
    615         /** Precision context used for floating point comparisons. */
    616         private final DoublePrecisionContext precision;
    617
    618         /** The start point from the point sequence. */
    619         private Vector3D startPt;
    620
    621         /** The previous point from the point sequence. */
    622         private Vector3D prevPt;
    623
    624         /** The previous vector from the point sequence, preceding from the {@code startPt} to {@code prevPt}. */
    625         private Vector3D prevVector;
    626
    627         /** The computed {@code normal} vector for the plane. */
    628         private Vector3D.Unit normal;
    629
    630         /** The x component of the sum of all cross products from adjacent vectors in the point sequence. */
    631         private double crossSumX;
    632
    633         /** The y component of the sum of all cross products from adjacent vectors in the point sequence. */
    634         private double crossSumY;
    635
    636         /** The z component of the sum of all cross products from adjacent vectors in the point sequence. */
    637         private double crossSumZ;
    638
    639         /** If true, an exception will be thrown if the point sequence is discovered to be non-convex. */
    640         private boolean requireConvex = false;
    641
    642         /** List that unique vertices discovered in the input sequence will be added to. */
    643         private List<Vector3D> uniqueVertexOutput;
    644
    645         /** Construct a new build instance for the given point sequence and precision context.
    646          * @param pts point sequence
    647          * @param precision precision context used to perform floating point comparisons
    648          */
    649         PlaneBuilder(final Collection<Vector3D> pts, final DoublePrecisionContext precision) {
    650             this.pts = pts;
    651             this.precision = precision;
    652         }
    653
    654         /** Build a plane from the configured point sequence.
    655          * @return a plane built from the configured point sequence
    656          * @throws IllegalArgumentException if the points do not define a plane
    657          */
    658         Plane build() {
    659             if (pts.size() < 3) {
    660                 throw nonPlanar();
    661             }
    662
    663             pts.forEach(this::processPoint);
    664
    665             return createPlane();
    666         }
    667
    668         /** Build a plane from the configured point sequence, validating that the points form a convex region
    669          * and adding all discovered unique points to the given list.
    670          * @param vertexOutput list that unique points discovered in the point sequence will be added to
    671          * @return a plane created from the configured point sequence
    672          * @throws IllegalArgumentException if the points do not define a plane or the {@code requireConvex}
    673          *      flag is true and the points do not define a convex area
    674          */
    675         Plane buildForConvexPolygon(final List<Vector3D> vertexOutput) {
    676             this.requireConvex = true;
    677             this.uniqueVertexOutput = vertexOutput;
    678
    679             return build();
    680         }
    681
    682         /** Process a point from the point sequence.
    683          * @param pt
    684          * @throws IllegalArgumentException if the points do not define a plane or the {@code requireConvex}
    685          *      flag is true and the points do not define a convex area
    686          */
    687         private void processPoint(final Vector3D pt) {
    688             if (prevPt == null) {
    689                 startPt = pt;
    690                 prevPt = pt;
    691
    692                 if (uniqueVertexOutput != null) {
    693                     uniqueVertexOutput.add(pt);
    694                 }
    695
    696             } else if (!prevPt.eq(pt, precision)) { // skip duplicate points
    697                 final Vector3D vec = startPt.vectorTo(pt);
    698
    699                 if (prevVector != null) {
    700                     processCrossProduct(prevVector.cross(vec));
    701                 }
    702
    703                 if (uniqueVertexOutput != null) {
    704                     uniqueVertexOutput.add(pt);
    705                 }
    706
    707                 prevPt = pt;
    708                 prevVector = vec;
    709             }
    710         }
    711
    712         /** Process the computed cross product of two vectors from the input point sequence. The vectors
    713          * start at the first point in the sequence and point to adjacent points later in the sequence.
    714          * @param cross the cross product of two vectors from the input point sequence
    715          * @throws IllegalArgumentException if the points do not define a plane or the {@code requireConvex}
    716          *      flag is true and the points do not define a convex area
    717          */
    718         private void processCrossProduct(final Vector3D cross) {
    719             crossSumX += cross.getX();
    720             crossSumY += cross.getY();
    721             crossSumZ += cross.getZ();
    722
    723             final double crossNorm = cross.norm();
    724
    725             if (!precision.eqZero(crossNorm)) {
    726                 // the cross product has non-zero magnitude
    727                 if (normal == null) {
    728                     // save the first non-zero cross product as our normal
    729                     normal = cross.normalize();
    730                 } else {
    731                     final double crossDot = normal.dot(cross) / crossNorm;
    732
    733                     // check non-planar before non-convex since the former is a more general type
    734                     // of issue
    735                     if (!precision.eq(1.0, Math.abs(crossDot))) {
    736                         throw nonPlanar();
    737                     } else if (requireConvex && crossDot < 0) {
    738                         throw nonConvex();
    739                     }
    740                 }
    741             }
    742         }
    743
    744         /** Construct the plane instance using the value gathered during point processing.
    745          * @return the created plane instance
    746          * @throws IllegalArgumentException if the point do not define a plane
    747          */
    748         private Plane createPlane() {
    749             if (normal == null) {
    750                 throw nonPlanar();
    751             }
    752
    753             // flip the normal if needed to match the overall orientation of the points
    754             if (normal.dot(Vector3D.of(crossSumX, crossSumY, crossSumZ)) < 0) {
    755                 normal = normal.negate();
    756             }
    757
    758             // construct the plane
    759             final double originOffset = -startPt.dot(normal);
    760
    761             return new Plane(normal, originOffset, precision);
    762         }
    763
    764         /** Return an exception with a message stating that the points given to this builder do not
    765          * define a plane.
    766          * @return an exception stating that the points do not define a plane
    767          */
    768         private IllegalArgumentException nonPlanar() {
    769             return new IllegalArgumentException("Points do not define a plane: " + pts);
    770         }
    771
    772         /** Return an exception with a message stating that the points given to this builder do not
    773          * define a convex region.
    774          * @return an exception stating that the points do not define a plane
    775          */
    776         private IllegalArgumentException nonConvex() {
    777             return new IllegalArgumentException("Points do not define a convex region: " + pts);
    778         }
    779     }
    780
    781     /** Class designed to create 3D regions by taking a 2D region and extruding from a base plane
    782      * through an extrusion vector. The ends ("top" and "bottom") of the extruded 3D region are flat
    783      * while the sides follow the boundaries of the original 2D region.
    784      */
    785     private static final class PlaneRegionExtruder {
    786         /** Base plane to extrude from. */
    787         private final EmbeddingPlane basePlane;
    788
    789         /** Vector to extrude along; the extruded plane is translated from the base plane by this amount. */
    790         private final Vector3D extrusionVector;
    791
    792         /** True if the extrusion vector points to the plus side of the base plane. */
    793         private final boolean extrudingOnPlusSide;
    794
    795         /** Precision context used to create boundaries. */
    796         private final DoublePrecisionContext precision;
    797
    798         /** Construct a new instance that performs extrusions from {@code basePlane} along {@code extrusionVector}.
    799          * @param basePlane base plane to extrude from
    800          * @param extrusionVector vector to extrude along
    801          * @param precision precision context used to construct boundaries
    802          * @throws IllegalArgumentException if the given extrusion vector and plane produce regions
    803          *      of zero size
    804          */
    805         PlaneRegionExtruder(final EmbeddingPlane basePlane, final Vector3D extrusionVector,
    806                 final DoublePrecisionContext precision) {
    807
    808             this.basePlane = basePlane;
    809
    810             // Extruded plane; this forms the end of the 3D region opposite the base plane.
    811             EmbeddingPlane extrudedPlane = basePlane.translate(extrusionVector);
    812
    813             if (basePlane.contains(extrudedPlane)) {
    814                 throw new IllegalArgumentException(
    815                         "Extrusion vector produces regions of zero size: extrusionVector= " +
    816                                 extrusionVector + ",  plane= " + basePlane);
    817             }
    818
    819             this.extrusionVector = extrusionVector;
    820             this.extrudingOnPlusSide = basePlane.getNormal().dot(extrusionVector) > 0;
    821
    822             this.precision = precision;
    823         }
    824
    825         /** Extrude the given 2D BSP tree using the configured base plane and extrusion vector.
    826          * @param subspaceRegion region to extrude
    827          * @return the boundaries of the extruded region
    828          */
    829         public List<PlaneConvexSubset> extrude(final RegionBSPTree2D subspaceRegion) {
    830             final List<PlaneConvexSubset> extrudedBoundaries = new ArrayList<>();
    831
    832             // add the boundaries
    833             addEnds(subspaceRegion, extrudedBoundaries);
    834             addSides(subspaceRegion, extrudedBoundaries);
    835
    836             return extrudedBoundaries;
    837         }
    838
    839         /** Add the end ("top" and "bottom") of the extruded subspace region to the result list.
    840          * @param subspaceRegion subspace region being extruded.
    841          * @param result list to add the boundary results to
    842          */
    843         private void addEnds(final RegionBSPTree2D subspaceRegion, final List<PlaneConvexSubset> result) {
    844             // add the base boundaries
    845             final List<ConvexArea> baseAreas = subspaceRegion.toConvex();
    846
    847             final List<PlaneConvexSubset> baseList = new ArrayList<>(baseAreas.size());
    848             final List<PlaneConvexSubset> extrudedList = new ArrayList<>(baseAreas.size());
    849
    850             final AffineTransformMatrix3D extrudeTransform = AffineTransformMatrix3D.createTranslation(extrusionVector);
    851
    852             PlaneConvexSubset base;
    853             for (final ConvexArea area : baseAreas) {
    854                 base = subsetFromConvexArea(basePlane, area);
    855                 if (extrudingOnPlusSide) {
    856                     base = base.reverse();
    857                 }
    858
    859                 baseList.add(base);
    860                 extrudedList.add(base.transform(extrudeTransform).reverse());
    861             }
    862
    863             result.addAll(baseList);
    864             result.addAll(extrudedList);
    865         }
    866
    867         /** Add the side boundaries of the extruded region to the result list.
    868          * @param subspaceRegion subspace region being extruded.
    869          * @param result list to add the boundary results to
    870          */
    871         private void addSides(final RegionBSPTree2D subspaceRegion, final List<PlaneConvexSubset> result) {
    872             Vector2D subStartPt;
    873             Vector2D subEndPt;
    874
    875             PlaneConvexSubset boundary;
    876             for (final LinePath path : subspaceRegion.getBoundaryPaths()) {
    877                 for (final LineConvexSubset lineSubset : path.getElements()) {
    878                     subStartPt = lineSubset.getStartPoint();
    879                     subEndPt = lineSubset.getEndPoint();
    880
    881                     boundary = (subStartPt != null && subEndPt != null) ?
    882                             extrudeSideFinite(basePlane.toSpace(subStartPt), basePlane.toSpace(subEndPt)) :
    883                             extrudeSideInfinite(lineSubset);
    884
    885                     result.add(boundary);
    886                 }
    887             }
    888         }
    889
    890         /** Extrude a single, finite boundary forming one of the sides of the extruded region.
    891          * @param startPt start point of the boundary
    892          * @param endPt end point of the boundary
    893          * @return the extruded region side boundary
    894          */
    895         private ConvexPolygon3D extrudeSideFinite(final Vector3D startPt, final Vector3D endPt) {
    896             final Vector3D extrudedStartPt = startPt.add(extrusionVector);
    897             final Vector3D extrudedEndPt = endPt.add(extrusionVector);
    898
    899             final List<Vector3D> vertices = extrudingOnPlusSide ?
    900                     Arrays.asList(startPt, endPt, extrudedEndPt, extrudedStartPt) :
    901                     Arrays.asList(startPt, extrudedStartPt, extrudedEndPt, endPt);
    902
    903             return convexPolygonFromVertices(vertices, precision);
    904         }
    905
    906         /** Extrude a single, infinite boundary forming one of the sides of the extruded region.
    907          * @param lineSubset line subset to extrude
    908          * @return the extruded region side boundary
    909          */
    910         private PlaneConvexSubset extrudeSideInfinite(final LineConvexSubset lineSubset) {
    911             final Vector2D subLinePt = lineSubset.getLine().getOrigin();
    912             final Vector2D subLineDir = lineSubset.getLine().getDirection();
    913
    914             final Vector3D linePt = basePlane.toSpace(subLinePt);
    915             final Vector3D lineDir = linePt.vectorTo(basePlane.toSpace(subLinePt.add(subLineDir)));
    916
    917             final EmbeddingPlane sidePlane;
    918             if (extrudingOnPlusSide) {
    919                 sidePlane = fromPointAndPlaneVectors(linePt, lineDir, extrusionVector, precision);
    920             } else {
    921                 sidePlane = fromPointAndPlaneVectors(linePt, extrusionVector, lineDir, precision);
    922             }
    923
    924             final Vector2D sideLineOrigin = sidePlane.toSubspace(linePt);
    925             final Vector2D sideLineDir = sideLineOrigin.vectorTo(sidePlane.toSubspace(linePt.add(lineDir)));
    926
    927             final Vector2D extrudedSideLineOrigin = sidePlane.toSubspace(linePt.add(extrusionVector));
    928
    929             final Vector2D sideExtrusionDir = sidePlane.toSubspace(sidePlane.getOrigin().add(extrusionVector))
    930                     .normalize();
    931
    932             // construct a list of lines forming the bounds of the extruded subspace region
    933             final List<Line> lines = new ArrayList<>();
    934
    935             // add the top and bottom lines (original and extruded)
    936             if (extrudingOnPlusSide) {
    937                 lines.add(Lines.fromPointAndDirection(sideLineOrigin, sideLineDir, precision));
    938                 lines.add(Lines.fromPointAndDirection(extrudedSideLineOrigin, sideLineDir.negate(), precision));
    939             } else {
    940                 lines.add(Lines.fromPointAndDirection(sideLineOrigin, sideLineDir.negate(), precision));
    941                 lines.add(Lines.fromPointAndDirection(extrudedSideLineOrigin, sideLineDir, precision));
    942             }
    943
    944             // if we have a point on the original line, then connect the two
    945             final Vector2D startPt = lineSubset.getStartPoint();
    946             final Vector2D endPt = lineSubset.getEndPoint();
    947             if (startPt != null) {
    948                 lines.add(Lines.fromPointAndDirection(
    949                         sidePlane.toSubspace(basePlane.toSpace(startPt)),
    950                         extrudingOnPlusSide ? sideExtrusionDir.negate() : sideExtrusionDir,
    951                         precision));
    952             } else if (endPt != null) {
    953                 lines.add(Lines.fromPointAndDirection(
    954                         sidePlane.toSubspace(basePlane.toSpace(endPt)),
    955                         extrudingOnPlusSide ? sideExtrusionDir : sideExtrusionDir.negate(),
    956                         precision));
    957             }
    958
    959             return subsetFromConvexArea(sidePlane, ConvexArea.fromBounds(lines));
    960         }
    961     }
    962 }