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.spherical.twod;
    18
    19  import java.util.ArrayList;
    20  import java.util.Arrays;
    21  import java.util.Collection;
    22  import java.util.Collections;
    23  import java.util.Iterator;
    24  import java.util.List;
    25  import java.util.stream.Collectors;
    26  import java.util.stream.Stream;
    27
    28  import org.apache.commons.geometry.core.Transform;
    29  import org.apache.commons.geometry.core.partitioning.AbstractConvexHyperplaneBoundedRegion;
    30  import org.apache.commons.geometry.core.partitioning.Hyperplane;
    31  import org.apache.commons.geometry.core.partitioning.HyperplaneConvexSubset;
    32  import org.apache.commons.geometry.core.partitioning.Split;
    33  import org.apache.commons.geometry.core.precision.DoublePrecisionContext;
    34  import org.apache.commons.geometry.euclidean.threed.Vector3D;
    35  import org.apache.commons.numbers.angle.PlaneAngleRadians;
    36
    37  /** Class representing a convex area in 2D spherical space. The boundaries of this
    38   * area, if any, are composed of convex great circle arcs.
    39   */
    40  public final class ConvexArea2S extends AbstractConvexHyperplaneBoundedRegion<Point2S, GreatArc>
    41      implements BoundarySource2S {
    42      /** Instance representing the full spherical area. */
    43      private static final ConvexArea2S FULL = new ConvexArea2S(Collections.emptyList());
    44
    45      /** Constant containing the area of the full spherical space. */
    46      private static final double FULL_SIZE = 4 * PlaneAngleRadians.PI;
    47
    48      /** Constant containing the area of half of the spherical space. */
    49      private static final double HALF_SIZE = PlaneAngleRadians.TWO_PI;
    50
    51      /** Empirically determined threshold for computing the weighted centroid vector using the
    52       * triangle fan approach. Areas with boundary sizes under this value use the triangle fan
    53       * method to increase centroid accuracy.
    54       */
    55      private static final double TRIANGLE_FAN_CENTROID_COMPUTE_THRESHOLD = 1e-2;
    56
    57      /** Construct an instance from its boundaries. Callers are responsible for ensuring
    58       * that the given path represents the boundary of a convex area. No validation is
    59       * performed.
    60       * @param boundaries the boundaries of the convex area
    61       */
    62      private ConvexArea2S(final List<GreatArc> boundaries) {
    63          super(boundaries);
    64      }
    65
    66      /** {@inheritDoc} */
    67      @Override
    68      public Stream<GreatArc> boundaryStream() {
    69          return getBoundaries().stream();
    70      }
    71
    72      /** Get a path instance representing the boundary of the area. The path is oriented
    73       * so that the minus sides of the arcs lie on the inside of the area.
    74       * @return the boundary path of the area
    75       */
    76      public GreatArcPath getBoundaryPath() {
    77          final List<GreatArcPath> paths = InteriorAngleGreatArcConnector.connectMinimized(getBoundaries());
    78          if (paths.isEmpty()) {
    79              return GreatArcPath.empty();
    80          }
    81
    82          return paths.get(0);
    83      }
    84
    85      /** Get an array of interior angles for the area. An empty array is returned if there
    86       * are no boundary intersections (ie, it has only one boundary or no boundaries at all).
    87       *
    88       * <p>The order of the angles corresponds with the order of the boundaries returned
    89       * by {@link #getBoundaries()}: if {@code i} is an index into the boundaries list,
    90       * then {@code angles[i]} is the angle between boundaries {@code i} and {@code (i+1) % boundariesSize}.</p>
    91       * @return an array of interior angles for the area
    92       */
    93      public double[] getInteriorAngles() {
    94          final List<GreatArc> arcs = getBoundaryPath().getArcs();
    95          final int numSides = arcs.size();
    96
    97          if (numSides < 2) {
    98              return new double[0];
    99          }
    100
    101         final double[] angles = new double[numSides];
    102
    103         GreatArc current;
    104         GreatArc next;
    105         for (int i = 0; i < numSides; ++i) {
    106             current = arcs.get(i);
    107             next = arcs.get((i + 1) % numSides);
    108
    109             angles[i] = PlaneAngleRadians.PI - current.getCircle()
    110                     .angle(next.getCircle(), current.getEndPoint());
    111         }
    112
    113         return angles;
    114     }
    115
    116     /** {@inheritDoc} */
    117     @Override
    118     public double getSize() {
    119         final int numSides = getBoundaries().size();
    120
    121         if (numSides == 0) {
    122             return FULL_SIZE;
    123         } else if (numSides == 1) {
    124             return HALF_SIZE;
    125         } else {
    126             // use the extended version of Girard's theorem
    127             // https://en.wikipedia.org/wiki/Spherical_trigonometry#Girard's_theorem
    128             final double[] angles = getInteriorAngles();
    129             final double sum = Arrays.stream(angles).sum();
    130
    131             return sum - ((angles.length - 2) * PlaneAngleRadians.PI);
    132         }
    133     }
    134
    135     /** {@inheritDoc} */
    136     @Override
    137     public Point2S getCentroid() {
    138         final Vector3D weighted = getWeightedCentroidVector();
    139         return weighted == null ? null : Point2S.from(weighted);
    140     }
    141
    142     /** Return the weighted centroid vector of the area. The returned vector points in the direction of the
    143      * centroid point on the surface of the unit sphere with the length of the vector proportional to the
    144      * effective mass of the area at the centroid. By adding the weighted centroid vectors of multiple
    145      * convex areas, a single centroid can be computed for the combined area.
    146      * @return weighted centroid vector.
    147      * @see <a href="https://archive.org/details/centroidinertiat00broc">
    148      *  <em>The Centroid and Inertia Tensor for a Spherical Triangle</em> - John E. Brock</a>
    149      */
    150     Vector3D getWeightedCentroidVector() {
    151         final List<GreatArc> arcs = getBoundaries();
    152         final int numBoundaries = arcs.size();
    153
    154         switch (numBoundaries) {
    155         case 0:
    156             // full space; no centroid
    157             return null;
    158         case 1:
    159             // hemisphere
    160             return computeHemisphereWeightedCentroidVector(arcs.get(0));
    161         case 2:
    162             // lune
    163             return computeLuneWeightedCentroidVector(arcs.get(0), arcs.get(1));
    164         default:
    165             // triangle or other convex polygon
    166             if (getBoundarySize() < TRIANGLE_FAN_CENTROID_COMPUTE_THRESHOLD) {
    167                 return computeTriangleFanWeightedCentroidVector(arcs);
    168             }
    169
    170             return computeArcPoleWeightedCentroidVector(arcs);
    171         }
    172     }
    173
    174     /** {@inheritDoc} */
    175     @Override
    176     public Split<ConvexArea2S> split(final Hyperplane<Point2S> splitter) {
    177         return splitInternal(splitter, this, GreatArc.class, ConvexArea2S::new);
    178     }
    179
    180     /** Return a BSP tree representing the same region as this instance.
    181      */
    182     @Override
    183     public RegionBSPTree2S toTree() {
    184         return RegionBSPTree2S.from(getBoundaries(), true);
    185     }
    186
    187     /** Return a new instance transformed by the argument.
    188      * @param transform transform to apply
    189      * @return a new instance transformed by the argument
    190      */
    191     public ConvexArea2S transform(final Transform<Point2S> transform) {
    192         return transformInternal(transform, this, GreatArc.class, ConvexArea2S::new);
    193     }
    194
    195     /** {@inheritDoc} */
    196     @Override
    197     public GreatArc trim(final HyperplaneConvexSubset<Point2S> sub) {
    198         return (GreatArc) super.trim(sub);
    199     }
    200
    201     /** Return an instance representing the full spherical 2D space.
    202      * @return an instance representing the full spherical 2D space.
    203      */
    204     public static ConvexArea2S full() {
    205         return FULL;
    206     }
    207
    208     /** Construct a convex area by creating great circles between adjacent vertices. The vertices must be given
    209      * in a counter-clockwise around order the interior of the shape. If the area is intended to be closed, the
    210      * beginning point must be repeated at the end of the path.
    211      * @param vertices vertices to use to construct the area
    212      * @param precision precision context used to create new great circle instances
    213      * @return a convex area constructed using great circles between adjacent vertices
    214      * @see #fromVertexLoop(Collection, DoublePrecisionContext)
    215      */
    216     public static ConvexArea2S fromVertices(final Collection<Point2S> vertices,
    217             final DoublePrecisionContext precision) {
    218         return fromVertices(vertices, false, precision);
    219     }
    220
    221     /** Construct a convex area by creating great circles between adjacent vertices. An implicit great circle is
    222      * created between the last vertex given and the first one, if needed. The vertices must be given in a
    223      * counter-clockwise around order the interior of the shape.
    224      * @param vertices vertices to use to construct the area
    225      * @param precision precision context used to create new great circles instances
    226      * @return a convex area constructed using great circles between adjacent vertices
    227      * @see #fromVertices(Collection, DoublePrecisionContext)
    228      */
    229     public static ConvexArea2S fromVertexLoop(final Collection<Point2S> vertices,
    230             final DoublePrecisionContext precision) {
    231         return fromVertices(vertices, true, precision);
    232     }
    233
    234     /** Construct a convex area from great circles between adjacent vertices.
    235      * @param vertices vertices to use to construct the area
    236      * @param close if true, an additional great circle will be created between the last and first vertex
    237      * @param precision precision context used to create new great circle instances
    238      * @return a convex area constructed using great circles between adjacent vertices
    239      */
    240     public static ConvexArea2S fromVertices(final Collection<Point2S> vertices, final boolean close,
    241             final DoublePrecisionContext precision) {
    242
    243         if (vertices.isEmpty()) {
    244             return full();
    245         }
    246
    247         final List<GreatCircle> circles = new ArrayList<>();
    248
    249         Point2S first = null;
    250         Point2S prev = null;
    251         Point2S cur = null;
    252
    253         for (final Point2S vertex : vertices) {
    254             cur = vertex;
    255
    256             if (first == null) {
    257                 first = cur;
    258             }
    259
    260             if (prev != null && !cur.eq(prev, precision)) {
    261                 circles.add(GreatCircles.fromPoints(prev, cur, precision));
    262             }
    263
    264             prev = cur;
    265         }
    266
    267         if (close && cur != null && !cur.eq(first, precision)) {
    268             circles.add(GreatCircles.fromPoints(cur, first, precision));
    269         }
    270
    271         if (!vertices.isEmpty() && circles.isEmpty()) {
    272             throw new IllegalStateException("Unable to create convex area: only a single unique vertex provided");
    273         }
    274
    275         return fromBounds(circles);
    276     }
    277
    278     /** Construct a convex area from an arc path. The area represents the intersection of all of the negative
    279      * half-spaces of the great circles in the path. The boundaries of the returned area may therefore not match
    280      * the arcs in the path.
    281      * @param path path to construct the area from
    282      * @return a convex area constructed from the great circles in the given path
    283      */
    284     public static ConvexArea2S fromPath(final GreatArcPath path) {
    285         final List<GreatCircle> bounds = path.getArcs().stream()
    286             .map(GreatArc::getCircle)
    287             .collect(Collectors.toList());
    288
    289         return fromBounds(bounds);
    290     }
    291
    292     /** Create a convex area formed by the intersection of the negative half-spaces of the
    293      * given bounding great circles. The returned instance represents the area that is on the
    294      * minus side of all of the given circles. Note that this method does not support areas
    295      * of zero size (ie, infinitely thin areas or points.)
    296      * @param bounds great circles used to define the convex area
    297      * @return a new convex area instance representing the area on the minus side of all
    298      *      of the bounding great circles or an instance representing the full area if no
    299      *      circles are given
    300      * @throws IllegalArgumentException if the given set of bounding great circles do not form a convex area,
    301      *      meaning that there is no region that is on the minus side of all of the bounding circles.
    302      */
    303     public static ConvexArea2S fromBounds(final GreatCircle... bounds) {
    304         return fromBounds(Arrays.asList(bounds));
    305     }
    306
    307     /** Create a convex area formed by the intersection of the negative half-spaces of the
    308      * given bounding great circles. The returned instance represents the area that is on the
    309      * minus side of all of the given circles. Note that this method does not support areas
    310      * of zero size (ie, infinitely thin areas or points.)
    311      * @param bounds great circles used to define the convex area
    312      * @return a new convex area instance representing the area on the minus side of all
    313      *      of the bounding great circles or an instance representing the full area if no
    314      *      circles are given
    315      * @throws IllegalArgumentException if the given set of bounding great circles do not form a convex area,
    316      *      meaning that there is no region that is on the minus side of all of the bounding circles.
    317      */
    318     public static ConvexArea2S fromBounds(final Iterable<GreatCircle> bounds) {
    319         final List<GreatArc> arcs = new ConvexRegionBoundaryBuilder<>(GreatArc.class).build(bounds);
    320         return arcs.isEmpty() ?
    321                 full() :
    322                 new ConvexArea2S(arcs);
    323     }
    324
    325     /** Compute the weighted centroid vector for the hemisphere formed by the given arc.
    326      * @param arc arc defining the hemisphere
    327      * @return the weighted centroid vector for the hemisphere
    328      * @see #getWeightedCentroidVector()
    329      */
    330     private static Vector3D computeHemisphereWeightedCentroidVector(final GreatArc arc) {
    331         return arc.getCircle().getPole().withNorm(HALF_SIZE);
    332     }
    333
    334     /** Compute the weighted centroid vector for the lune formed by the given arcs.
    335      * @param a first arc for the lune
    336      * @param b second arc for the lune
    337      * @return the weighted centroid vector for the lune
    338      * @see #getWeightedCentroidVector()
    339      */
    340     private static Vector3D computeLuneWeightedCentroidVector(final GreatArc a, final GreatArc b) {
    341         final Point2S aMid = a.getCentroid();
    342         final Point2S bMid = b.getCentroid();
    343
    344         // compute the centroid vector as the exact center of the lune to avoid inaccurate
    345         // results with very small regions
    346         final Vector3D.Unit centroid = aMid.slerp(bMid, 0.5).getVector();
    347
    348         // compute the weight using the reverse of the algorithm from computeArcPoleWeightedCentroidVector()
    349         final double weight =
    350             (a.getSize() * centroid.dot(a.getCircle().getPole())) +
    351             (b.getSize() * centroid.dot(b.getCircle().getPole()));
    352
    353         return centroid.withNorm(weight);
    354     }
    355
    356     /** Compute the weighted centroid vector for the triangle or polygon formed by the given arcs
    357      * by adding together the arc pole vectors multiplied by their respective arc lengths. This
    358      * algorithm is described in the paper <a href="https://archive.org/details/centroidinertiat00broc">
    359      * <em>The Centroid and Inertia Tensor for a Spherical Triangle</em></a> by John E Brock.
    360      *
    361      * <p>Note: This algorithm works well in general but is susceptible to floating point errors
    362      * on very small areas. In these cases, the computed centroid may not be in the expected location
    363      * and may even be outside of the area. The {@link #computeTriangleFanWeightedCentroidVector(List)}
    364      * method can produce more accurate results in these cases.</p>
    365      * @param arcs boundary arcs for the area
    366      * @return the weighted centroid vector for the area
    367      * @see #computeTriangleFanWeightedCentroidVector(List)
    368      */
    369     private static Vector3D computeArcPoleWeightedCentroidVector(final List<GreatArc> arcs) {
    370         Vector3D centroid = Vector3D.ZERO;
    371
    372         Vector3D arcContribution;
    373         for (final GreatArc arc : arcs) {
    374             arcContribution = arc.getCircle().getPole().withNorm(arc.getSize());
    375
    376             centroid = centroid.add(arcContribution);
    377         }
    378
    379         return centroid;
    380     }
    381
    382     /** Compute the weighted centroid vector for the triangle or polygon formed by the given arcs
    383      * using a triangle fan approach. This method is specifically designed for use with areas of very small size,
    384      * where use of the standard algorithm from {@link ##computeArcPoleWeightedCentroidVector(List))} can produce
    385      * inaccurate results. The algorithm proceeds as follows:
    386      * <ol>
    387      *  <li>The polygon is divided into spherical triangles using a triangle fan.</li>
    388      *  <li>For each triangle, the vectors of the 3 spherical points are added together to approximate the direction
    389      *      of the spherical centroid. This ensures that the computed centroid lies within the area.</li>
    390      *  <li>The length of the weighted centroid vector is determined by computing the sum of the contributions that
    391      *      each arc in the triangle would make to the centroid using the algorithm from
    392      *      {@link ##computeArcPoleWeightedCentroidVector(List)}. This essentially performs part of that algorithm in
    393      *      reverse: given a centroid direction, compute the contribution that each arc makes.</li>
    394      *  <li>The sum of the weighted centroid vectors for each triangle is computed and returned.</li>
    395      * </ol>
    396      * @param arcs boundary arcs for the area; must contain at least 3 arcs
    397      * @return the weighted centroid vector for the area
    398      * @see ##computeArcPoleWeightedCentroidVector(List)
    399      */
    400     private static Vector3D computeTriangleFanWeightedCentroidVector(final List<GreatArc> arcs) {
    401         final Iterator<GreatArc> arcIt = arcs.iterator();
    402
    403         final Point2S p0 = arcIt.next().getStartPoint();
    404         final Vector3D.Unit v0 = p0.getVector();
    405
    406         Vector3D areaCentroid = Vector3D.ZERO;
    407
    408         GreatArc arc;
    409         Point2S p1;
    410         Point2S p2;
    411         Vector3D.Unit v1;
    412         Vector3D.Unit v2;
    413         Vector3D.Unit triangleCentroid;
    414         double triangleCentroidLen;
    415         while (arcIt.hasNext()) {
    416             arc = arcIt.next();
    417
    418             if (!arc.contains(p0)) {
    419                 p1 = arc.getStartPoint();
    420                 p2 = arc.getEndPoint();
    421
    422                 v1 = p1.getVector();
    423                 v2 = p2.getVector();
    424
    425                 triangleCentroid = v0.add(v1).add(v2).normalize();
    426                 triangleCentroidLen =
    427                         computeArcCentroidContribution(v0, v1, triangleCentroid) +
    428                         computeArcCentroidContribution(v1, v2, triangleCentroid) +
    429                         computeArcCentroidContribution(v2, v0, triangleCentroid);
    430
    431                 areaCentroid = areaCentroid.add(triangleCentroid.withNorm(triangleCentroidLen));
    432             }
    433         }
    434
    435         return areaCentroid;
    436     }
    437
    438     /** Compute the contribution made by a single arc to a weighted centroid vector.
    439      * @param a first point in the arc
    440      * @param b second point in the arc
    441      * @param triangleCentroid the centroid vector for the area
    442      * @return the contribution made by the arc {@code ab} to the length of the weighted centroid vector
    443      */
    444     private static double computeArcCentroidContribution(final Vector3D.Unit a, final Vector3D.Unit b,
    445             final Vector3D.Unit triangleCentroid) {
    446         final double arcLength = a.angle(b);
    447         final Vector3D.Unit planeNormal = a.cross(b).normalize();
    448
    449         return arcLength * triangleCentroid.dot(planeNormal);
    450     }
    451 }