ConvexArea2S.java

  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;

  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;

  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.euclidean.threed.Vector3D;
  34. import org.apache.commons.numbers.angle.Angle;
  35. import org.apache.commons.numbers.core.Precision;

  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());

  45.     /** Constant containing the area of the full spherical space. */
  46.     private static final double FULL_SIZE = 4 * Math.PI;

  48.     /** Constant containing the area of half of the spherical space. */
  49.     private static final double HALF_SIZE = Angle.TWO_PI;

  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;

  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.     }

  66.     /** {@inheritDoc} */
  67.     @Override
  68.     public Stream<GreatArc> boundaryStream() {
  69.         return getBoundaries().stream();
  70.     }

  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.         }

  82.         return paths.get(0);
  83.     }

  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();

  97.         if (numSides < 2) {
  98.             return new double[0];
  99.         }

  101.         final double[] angles = new double[numSides];

  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);

  109.             angles[i] = Math.PI - current.getCircle()
  110.                     .angle(next.getCircle(), current.getEndPoint());
  111.         }

  113.         return angles;
  114.     }

  116.     /** {@inheritDoc} */
  117.     @Override
  118.     public double getSize() {
  119.         final int numSides = getBoundaries().size();

  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();

  131.             return sum - ((angles.length - 2) * Math.PI);
  132.         }
  133.     }

  135.     /** {@inheritDoc} */
  136.     @Override
  137.     public Point2S getCentroid() {
  138.         final Vector3D weighted = getWeightedCentroidVector();
  139.         return weighted == null ? null : Point2S.from(weighted);
  140.     }

  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();

  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.             }

  170.             return computeArcPoleWeightedCentroidVector(arcs);
  171.         }
  172.     }

  174.     /** {@inheritDoc} */
  175.     @Override
  176.     public Split<ConvexArea2S> split(final Hyperplane<Point2S> splitter) {
  177.         return splitInternal(splitter, this, GreatArc.class, ConvexArea2S::new);
  178.     }

  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.     }

  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.     }

  195.     /** {@inheritDoc} */
  196.     @Override
  197.     public GreatArc trim(final HyperplaneConvexSubset<Point2S> sub) {
  198.         return (GreatArc) super.trim(sub);
  199.     }

  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.     }

  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, Precision.DoubleEquivalence)
  215.      */
  216.     public static ConvexArea2S fromVertices(final Collection<Point2S> vertices,
  217.             final Precision.DoubleEquivalence precision) {
  218.         return fromVertices(vertices, false, precision);
  219.     }

  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, Precision.DoubleEquivalence)
  228.      */
  229.     public static ConvexArea2S fromVertexLoop(final Collection<Point2S> vertices,
  230.             final Precision.DoubleEquivalence precision) {
  231.         return fromVertices(vertices, true, precision);
  232.     }

  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 Precision.DoubleEquivalence precision) {

  243.         if (vertices.isEmpty()) {
  244.             return full();
  245.         }

  247.         final List<GreatCircle> circles = new ArrayList<>();

  249.         Point2S first = null;
  250.         Point2S prev = null;
  251.         Point2S cur = null;

  253.         for (final Point2S vertex : vertices) {
  254.             cur = vertex;

  256.             if (first == null) {
  257.                 first = cur;
  258.             }

  260.             if (prev != null && !cur.eq(prev, precision)) {
  261.                 circles.add(GreatCircles.fromPoints(prev, cur, precision));
  262.             }

  264.             prev = cur;
  265.         }

  267.         if (close && cur != null && !cur.eq(first, precision)) {
  268.             circles.add(GreatCircles.fromPoints(cur, first, precision));
  269.         }

  271.         if (!vertices.isEmpty() && circles.isEmpty()) {
  272.             throw new IllegalStateException("Unable to create convex area: only a single unique vertex provided");
  273.         }

  275.         return fromBounds(circles);
  276.     }

  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());

  289.         return fromBounds(bounds);
  290.     }

  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.     }

  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.     }

  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.     }

  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();

  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();

  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()));

  353.         return centroid.withNorm(weight);
  354.     }

  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.         final Vector3D.Sum centroid = Vector3D.Sum.create();

  372.         for (final GreatArc arc : arcs) {
  373.             centroid.addScaled(arc.getSize(), arc.getCircle().getPole());
  374.         }

  376.         return centroid.get();
  377.     }

  379.     /** Compute the weighted centroid vector for the triangle or polygon formed by the given arcs
  380.      * using a triangle fan approach. This method is specifically designed for use with areas of very small size,
  381.      * where use of the standard algorithm from {@link ##computeArcPoleWeightedCentroidVector(List))} can produce
  382.      * inaccurate results. The algorithm proceeds as follows:
  383.      * <ol>
  384.      *  <li>The polygon is divided into spherical triangles using a triangle fan.</li>
  385.      *  <li>For each triangle, the vectors of the 3 spherical points are added together to approximate the direction
  386.      *      of the spherical centroid. This ensures that the computed centroid lies within the area.</li>
  387.      *  <li>The length of the weighted centroid vector is determined by computing the sum of the contributions that
  388.      *      each arc in the triangle would make to the centroid using the algorithm from
  389.      *      {@link ##computeArcPoleWeightedCentroidVector(List)}. This essentially performs part of that algorithm in
  390.      *      reverse: given a centroid direction, compute the contribution that each arc makes.</li>
  391.      *  <li>The sum of the weighted centroid vectors for each triangle is computed and returned.</li>
  392.      * </ol>
  393.      * @param arcs boundary arcs for the area; must contain at least 3 arcs
  394.      * @return the weighted centroid vector for the area
  395.      * @see ##computeArcPoleWeightedCentroidVector(List)
  396.      */
  397.     private static Vector3D computeTriangleFanWeightedCentroidVector(final List<GreatArc> arcs) {
  398.         final Iterator<GreatArc> arcIt = arcs.iterator();

  400.         final Point2S p0 = arcIt.next().getStartPoint();
  401.         final Vector3D.Unit v0 = p0.getVector();

  403.         final Vector3D.Sum areaCentroid = Vector3D.Sum.create();

  405.         GreatArc arc;
  406.         Point2S p1;
  407.         Point2S p2;
  408.         Vector3D.Unit v1;
  409.         Vector3D.Unit v2;
  410.         Vector3D.Unit triangleCentroid;
  411.         double triangleCentroidLen;
  412.         while (arcIt.hasNext()) {
  413.             arc = arcIt.next();

  415.             if (!arc.contains(p0)) {
  416.                 p1 = arc.getStartPoint();
  417.                 p2 = arc.getEndPoint();

  419.                 v1 = p1.getVector();
  420.                 v2 = p2.getVector();

  422.                 triangleCentroid = Vector3D.Sum.create()
  423.                         .add(v0)
  424.                         .add(v1)
  425.                         .add(v2)
  426.                         .get().normalize();
  427.                 triangleCentroidLen =
  428.                         computeArcCentroidContribution(v0, v1, triangleCentroid) +
  429.                         computeArcCentroidContribution(v1, v2, triangleCentroid) +
  430.                         computeArcCentroidContribution(v2, v0, triangleCentroid);

  432.                 areaCentroid.addScaled(triangleCentroidLen, triangleCentroid);
  433.             }
  434.         }

  436.         return areaCentroid.get();
  437.     }

  439.     /** Compute the contribution made by a single arc to a weighted centroid vector.
  440.      * @param a first point in the arc
  441.      * @param b second point in the arc
  442.      * @param triangleCentroid the centroid vector for the area
  443.      * @return the contribution made by the arc {@code ab} to the length of the weighted centroid vector
  444.      */
  445.     private static double computeArcCentroidContribution(final Vector3D.Unit a, final Vector3D.Unit b,
  446.             final Vector3D.Unit triangleCentroid) {
  447.         final double arcLength = a.angle(b);
  448.         final Vector3D.Unit planeNormal = a.cross(b).normalize();

  450.         return arcLength * triangleCentroid.dot(planeNormal);
  451.     }
  452. }
Created with JaCoCo 0.8.5.201910111838