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.util.ArrayList;
    20  import java.util.Iterator;
    21  import java.util.List;
    22
    23  import org.apache.commons.geometry.core.RegionLocation;
    24  import org.apache.commons.geometry.core.partitioning.Hyperplane;
    25  import org.apache.commons.geometry.core.partitioning.HyperplaneLocation;
    26  import org.apache.commons.geometry.core.partitioning.Split;
    27  import org.apache.commons.geometry.euclidean.internal.Vectors;
    28  import org.apache.commons.geometry.euclidean.threed.line.Lines3D;
    29  import org.apache.commons.geometry.euclidean.twod.ConvexArea;
    30  import org.apache.commons.geometry.euclidean.twod.Vector2D;
    31  import org.apache.commons.numbers.core.Precision;
    32
    33  /** Abstract base class for {@link ConvexPolygon3D} implementations.
    34   */
    35  abstract class AbstractConvexPolygon3D extends AbstractPlaneSubset implements ConvexPolygon3D {
    36
    37      /** Plane containing the convex polygon. */
    38      private final Plane plane;
    39
    40      /** Simple constructor.
    41       * @param plane the plane containing the convex polygon
    42       */
    43      AbstractConvexPolygon3D(final Plane plane) {
    44          this.plane = plane;
    45      }
    46
    47      /** {@inheritDoc} */
    48      @Override
    49      public Plane getPlane() {
    50          return plane;
    51      }
    52
    53      /** {@inheritDoc}
    54       *
    55       *  <p>This method always returns {@code false}.</p>
    56       */
    57      @Override
    58      public boolean isFull() {
    59          return false;
    60      }
    61
    62      /** {@inheritDoc}
    63       *
    64       *  <p>This method always returns {@code false}.</p>
    65       */
    66      @Override
    67      public boolean isEmpty() {
    68          return false;
    69      }
    70
    71      /** {@inheritDoc} */
    72      @Override
    73      public double getSize() {
    74          // see http://geomalgorithms.com/a01-_area.html#3D-Planar-Polygons
    75          final List<Vector3D> vertices = getVertices();
    76
    77          double crossSumX = 0.0;
    78          double crossSumY = 0.0;
    79          double crossSumZ = 0.0;
    80
    81          Vector3D prevPt = vertices.get(vertices.size() - 1);
    82          Vector3D cross;
    83          for (final Vector3D curPt : vertices) {
    84              cross = prevPt.cross(curPt);
    85
    86              crossSumX += cross.getX();
    87              crossSumY += cross.getY();
    88              crossSumZ += cross.getZ();
    89
    90              prevPt = curPt;
    91          }
    92
    93          return 0.5 * plane.getNormal().dot(Vector3D.of(crossSumX, crossSumY, crossSumZ));
    94      }
    95
    96      /** {@inheritDoc} */
    97      @Override
    98      public Vector3D getCentroid() {
    99          final List<Vector3D> vertices = getVertices();
    100
    101         double areaSum = 0.0;
    102         double scaledCentroidSumX = 0.0;
    103         double scaledCentroidSumY = 0.0;
    104         double scaledCentroidSumZ = 0.0;
    105
    106         final Iterator<Vector3D> it = vertices.iterator();
    107
    108         final Vector3D startPt = it.next();
    109
    110         Vector3D prevPt = it.next();
    111         Vector3D curPt;
    112
    113         Vector3D prevVec = startPt.vectorTo(prevPt);
    114         Vector3D curVec;
    115
    116         double triArea;
    117         Vector3D triCentroid;
    118         while (it.hasNext()) {
    119             curPt = it.next();
    120             curVec = startPt.vectorTo(curPt);
    121
    122             triArea = 0.5 * prevVec.cross(curVec).norm();
    123             triCentroid = Vector3D.centroid(startPt, prevPt, curPt);
    124
    125             areaSum += triArea;
    126
    127             scaledCentroidSumX += triArea * triCentroid.getX();
    128             scaledCentroidSumY += triArea * triCentroid.getY();
    129             scaledCentroidSumZ += triArea * triCentroid.getZ();
    130
    131             prevPt = curPt;
    132             prevVec = curVec;
    133         }
    134
    135         if (areaSum > 0) {
    136             final double scale = 1 / areaSum;
    137             return Vector3D.of(
    138                         scale * scaledCentroidSumX,
    139                         scale * scaledCentroidSumY,
    140                         scale * scaledCentroidSumZ
    141                     );
    142         }
    143
    144         // zero area, which means that the points are all linear; return the point midway between the
    145         // min and max points
    146         final Vector3D min = Vector3D.min(vertices);
    147         final Vector3D max = Vector3D.max(vertices);
    148
    149         return min.lerp(max, 0.5);
    150     }
    151
    152     /** {@inheritDoc} */
    153     @Override
    154     public Bounds3D getBounds() {
    155         return Bounds3D.from(getVertices());
    156     }
    157
    158     /** {@inheritDoc} */
    159     @Override
    160     public RegionLocation classify(final Vector3D pt) {
    161         if (plane.contains(pt)) {
    162             final List<Vector3D> vertices = getVertices();
    163             final Precision.DoubleEquivalence precision = plane.getPrecision();
    164
    165             final Vector3D normal = plane.getNormal();
    166             Vector3D edgeVec;
    167             Vector3D edgePlusVec;
    168             Vector3D testVec;
    169
    170             Vector3D offsetVec;
    171             double offsetSign;
    172             double offset;
    173             int cmp;
    174
    175             boolean onBoundary = false;
    176
    177             Vector3D startVertex = vertices.get(vertices.size() - 1);
    178             for (final Vector3D nextVertex : vertices) {
    179
    180                 edgeVec = startVertex.vectorTo(nextVertex);
    181                 edgePlusVec = edgeVec.cross(normal);
    182
    183                 testVec = startVertex.vectorTo(pt);
    184
    185                 offsetVec = testVec.reject(edgeVec);
    186                 offsetSign = Math.signum(offsetVec.dot(edgePlusVec));
    187                 offset = offsetSign * offsetVec.norm();
    188
    189                 cmp = precision.compare(offset, 0.0);
    190                 if (cmp > 0) {
    191                     // the point is on the plus side (outside) of a boundary
    192                     return RegionLocation.OUTSIDE;
    193                 } else if (cmp == 0) {
    194                     onBoundary = true;
    195                 }
    196
    197                 startVertex = nextVertex;
    198             }
    199
    200             if (onBoundary) {
    201                 // the point is not on the outside of any boundaries and is directly on at least one
    202                 return RegionLocation.BOUNDARY;
    203             }
    204
    205             // the point is on the inside of all boundaries
    206             return RegionLocation.INSIDE;
    207         }
    208
    209         // the point is not on the plane
    210         return RegionLocation.OUTSIDE;
    211     }
    212
    213     /** {@inheritDoc} */
    214     @Override
    215     public Vector3D closest(final Vector3D pt) {
    216         final Vector3D normal = plane.getNormal();
    217         final Precision.DoubleEquivalence precision = plane.getPrecision();
    218
    219         final List<Vector3D> vertices = getVertices();
    220
    221         final Vector3D projPt = plane.project(pt);
    222
    223         Vector3D edgeVec;
    224         Vector3D edgePlusVec;
    225         Vector3D testVec;
    226
    227         Vector3D offsetVec;
    228         double offsetSign;
    229         double offset;
    230         int cmp;
    231
    232         Vector3D boundaryVec;
    233         double boundaryPointT;
    234         Vector3D boundaryPoint;
    235         double boundaryPointDistSq;
    236
    237         double closestBoundaryPointDistSq = Double.POSITIVE_INFINITY;
    238         Vector3D closestBoundaryPoint = null;
    239
    240         Vector3D startVertex = vertices.get(vertices.size() - 1);
    241         for (final Vector3D nextVertex : vertices) {
    242
    243             edgeVec = startVertex.vectorTo(nextVertex);
    244             edgePlusVec = edgeVec.cross(normal);
    245
    246             testVec = startVertex.vectorTo(projPt);
    247
    248             offsetVec = testVec.reject(edgeVec);
    249             offsetSign = Math.signum(offsetVec.dot(edgePlusVec));
    250             offset = offsetSign * offsetVec.norm();
    251
    252             cmp = precision.compare(offset, 0.0);
    253             if (cmp >= 0) {
    254                 // the point is on directly on the boundary or on its plus side; project the point onto the
    255                 // boundary, taking care to restrict the point to the actual extent of the boundary,
    256                 // and select the point with the shortest distance
    257                 boundaryVec = testVec.subtract(offsetVec);
    258                 boundaryPointT =
    259                         Math.signum(boundaryVec.dot(edgeVec)) * (boundaryVec.norm() / Vectors.checkedNorm(edgeVec));
    260                 boundaryPointT = Math.max(0, Math.min(1, boundaryPointT));
    261
    262                 boundaryPoint = startVertex.lerp(nextVertex, boundaryPointT);
    263
    264                 boundaryPointDistSq = boundaryPoint.distanceSq(projPt);
    265                 if (boundaryPointDistSq < closestBoundaryPointDistSq) {
    266                     closestBoundaryPointDistSq = boundaryPointDistSq;
    267                     closestBoundaryPoint = boundaryPoint;
    268                 }
    269             }
    270
    271             startVertex = nextVertex;
    272         }
    273
    274         if (closestBoundaryPoint != null) {
    275             // the point is on the outside of the polygon; return the closest point on the boundary
    276             return closestBoundaryPoint;
    277         }
    278
    279         // the projected point is on the inside of all boundaries and therefore on the inside of the subset
    280         return projPt;
    281     }
    282
    283     /** {@inheritDoc} */
    284     @Override
    285     public PlaneConvexSubset.Embedded getEmbedded() {
    286         final EmbeddingPlane embeddingPlane = plane.getEmbedding();
    287         final List<Vector2D> subspaceVertices = embeddingPlane.toSubspace(getVertices());
    288         final ConvexArea area = ConvexArea.convexPolygonFromVertices(subspaceVertices,
    289                 embeddingPlane.getPrecision());
    290
    291         return new EmbeddedAreaPlaneConvexSubset(embeddingPlane, area);
    292     }
    293
    294     /** {@inheritDoc} */
    295     @Override
    296     public Split<PlaneConvexSubset> split(final Hyperplane<Vector3D> splitter) {
    297         final Plane splitterPlane = (Plane) splitter;
    298         final List<Vector3D> vertices = getVertices();
    299
    300         final int size = vertices.size();
    301
    302         int minusPlusTransitionIdx = -1;
    303         Vector3D minusPlusInsertVertex = null;
    304
    305         int plusMinusTransitionIdx = -1;
    306         Vector3D plusMinusInsertVertex = null;
    307
    308         int transitionCount = 0;
    309
    310         Vector3D curVertex;
    311         HyperplaneLocation curLoc;
    312
    313         int lastSideIdx = -1;
    314         Vector3D lastSideVertex = null;
    315         HyperplaneLocation lastSideLoc = null;
    316
    317         int lastBoundaryIdx = -1;
    318
    319         for (int i = 0; i <= size || transitionCount == 1; ++i) {
    320
    321             curVertex = vertices.get(i % size);
    322             curLoc = splitter.classify(curVertex);
    323
    324             if (lastSideLoc == HyperplaneLocation.MINUS && curLoc == HyperplaneLocation.PLUS) {
    325                 // transitioned from minus side to plus side
    326                 minusPlusTransitionIdx = Math.max(lastSideIdx, lastBoundaryIdx);
    327                 ++transitionCount;
    328
    329                 if (lastBoundaryIdx < 0) {
    330                     // no shared boundary point; compute a new vertex
    331                     minusPlusInsertVertex = splitterPlane.intersection(
    332                             Lines3D.fromPoints(lastSideVertex, curVertex, splitterPlane.getPrecision()));
    333                 }
    334             } else if (lastSideLoc == HyperplaneLocation.PLUS && curLoc == HyperplaneLocation.MINUS) {
    335                 // transitioned from plus side to minus side
    336                 plusMinusTransitionIdx = Math.max(lastSideIdx, lastBoundaryIdx);
    337                 ++transitionCount;
    338
    339                 if (lastBoundaryIdx < 0) {
    340                     // no shared boundary point; compute a new vertex
    341                     plusMinusInsertVertex = splitterPlane.intersection(
    342                             Lines3D.fromPoints(lastSideVertex, curVertex, splitterPlane.getPrecision()));
    343                 }
    344             }
    345
    346             if (curLoc == HyperplaneLocation.ON) {
    347                 lastBoundaryIdx = i;
    348             } else {
    349                 lastBoundaryIdx = -1;
    350
    351                 lastSideIdx = i;
    352                 lastSideVertex = curVertex;
    353                 lastSideLoc = curLoc;
    354             }
    355         }
    356
    357         if (minusPlusTransitionIdx > -1 && plusMinusTransitionIdx > -1) {
    358             // we've split; compute the vertex list for each side
    359             final List<Vector3D> minusVertices =  buildPolygonSplitVertexList(
    360                     plusMinusTransitionIdx, plusMinusInsertVertex,
    361                     minusPlusTransitionIdx, minusPlusInsertVertex, vertices);
    362             final List<Vector3D> plusVertices = buildPolygonSplitVertexList(
    363                     minusPlusTransitionIdx, minusPlusInsertVertex,
    364                     plusMinusTransitionIdx, plusMinusInsertVertex, vertices);
    365
    366             // delegate back to the Planes factory methods to determine the concrete types
    367             // for each side of the split
    368             return new Split<>(
    369                     Planes.fromConvexPlanarVertices(plane, minusVertices),
    370                     Planes.fromConvexPlanarVertices(plane, plusVertices));
    371
    372         } else if (lastSideLoc == HyperplaneLocation.PLUS) {
    373             // we lie entirely on the plus side of the splitter
    374             return new Split<>(null, this);
    375         } else if (lastSideLoc == HyperplaneLocation.MINUS) {
    376             // we lie entirely on the minus side of the splitter
    377             return new Split<>(this, null);
    378         }
    379
    380         // we lie entirely on the splitter
    381         return new Split<>(null, null);
    382     }
    383
    384     /** Internal method for building a vertex list for one side of a split result. The method is
    385      * designed to make the fewest allocations possible.
    386      * @param enterIdx the index of the vertex from {@code vertices} immediately before the polygon transitioned
    387      *      to being fully entered into this side of the split result. If no point from {@code vertices} lay
    388      *      directly on the splitting plane while entering this side and a new vertex had to be computed for the
    389      *      split result, then this index will be the last vertex on the opposite side of the split. If a vertex
    390      *      did lie directly on the splitting plane, then this index will point to that vertex.
    391      * @param newEnterPt the newly-computed point to be added as the first vertex in the split result; may
    392      *      be null if no such point exists
    393      * @param exitIdx the index of the vertex from {@code vertices} immediately before the polygon transitioned
    394      *      to being fully exited from this side of the split result. If no point from {@code vertices} lay
    395      *      directly on the splitting plane while exiting this side and a new vertex had to be computed for the
    396      *      split result, then this index will be the last vertex on the this side of the split. If a vertex did
    397      *      lie directly on the splitting plane, then this index will point to that vertex.
    398      * @param newExitPt the newly-computed point to be added as the last vertex in the split result; may
    399      *      be null if no such point exists
    400      * @param vertices the original list of vertices that this split result originated from; this list is
    401      *      not modified by this operation
    402      * @return the list of vertices for the split result
    403      */
    404     private List<Vector3D> buildPolygonSplitVertexList(final int enterIdx, final Vector3D newEnterPt,
    405             final int exitIdx, final Vector3D newExitPt, final List<? extends Vector3D> vertices) {
    406
    407         final int size = vertices.size();
    408
    409         final boolean hasNewEnterPt = newEnterPt != null;
    410         final boolean hasNewExitPt = newExitPt != null;
    411
    412         final int startIdx = (hasNewEnterPt ? enterIdx + 1 : enterIdx) % size;
    413         final int endIdx = exitIdx % size;
    414
    415         final boolean hasWrappedIndices = endIdx < startIdx;
    416
    417         final int resultSize = (hasWrappedIndices ? endIdx + size : endIdx) - startIdx + 1;
    418         final List<Vector3D> result = new ArrayList<>(resultSize);
    419
    420         if (hasNewEnterPt) {
    421             result.add(newEnterPt);
    422         }
    423
    424         if (hasWrappedIndices) {
    425             result.addAll(vertices.subList(startIdx, size));
    426             result.addAll(vertices.subList(0, endIdx + 1));
    427         } else {
    428             result.addAll(vertices.subList(startIdx, endIdx + 1));
    429         }
    430
    431         if (hasNewExitPt) {
    432             result.add(newExitPt);
    433         }
    434
    435         return result;
    436     }
    437
    438     /** {@inheritDoc} */
    439     @Override
    440     public String toString() {
    441         final StringBuilder sb = new StringBuilder();
    442         sb.append(getClass().getSimpleName())
    443             .append("[normal= ")
    444             .append(getPlane().getNormal())
    445             .append(", vertices= ")
    446             .append(getVertices())
    447             .append(']');
    448
    449         return sb.toString();
    450     }
    451 }