commons-geometry-spherical/src/main/java/org/apache/commons/geometry/spherical/twod/ConvexArea2S.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.geometry.spherical.twod;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.Collections;
import java.util.Iterator;
import java.util.List;
import java.util.stream.Collectors;
import java.util.stream.Stream;

import org.apache.commons.geometry.core.Transform;
import org.apache.commons.geometry.core.partitioning.AbstractConvexHyperplaneBoundedRegion;
import org.apache.commons.geometry.core.partitioning.Hyperplane;
import org.apache.commons.geometry.core.partitioning.HyperplaneConvexSubset;
import org.apache.commons.geometry.core.partitioning.Split;
import org.apache.commons.geometry.core.precision.DoublePrecisionContext;
import org.apache.commons.geometry.euclidean.threed.Vector3D;
import org.apache.commons.numbers.angle.PlaneAngleRadians;

/** Class representing a convex area in 2D spherical space. The boundaries of this
 * area, if any, are composed of convex great circle arcs.
 */
public final class ConvexArea2S extends AbstractConvexHyperplaneBoundedRegion<Point2S, GreatArc>
    implements BoundarySource2S {
    /** Instance representing the full spherical area. */
    private static final ConvexArea2S FULL = new ConvexArea2S(Collections.emptyList());

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

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

    /** Empirically determined threshold for computing the weighted centroid vector using the
     * triangle fan approach. Areas with boundary sizes under this value use the triangle fan
     * method to increase centroid accuracy.
     */
    private static final double TRIANGLE_FAN_CENTROID_COMPUTE_THRESHOLD = 1e-2;

    /** Construct an instance from its boundaries. Callers are responsible for ensuring
     * that the given path represents the boundary of a convex area. No validation is
     * performed.
     * @param boundaries the boundaries of the convex area
     */
    private ConvexArea2S(final List<GreatArc> boundaries) {
        super(boundaries);
    }

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

    /** Get a path instance representing the boundary of the area. The path is oriented
     * so that the minus sides of the arcs lie on the inside of the area.
     * @return the boundary path of the area
     */
    public GreatArcPath getBoundaryPath() {
        final List<GreatArcPath> paths = InteriorAngleGreatArcConnector.connectMinimized(getBoundaries());
        if (paths.isEmpty()) {
            return GreatArcPath.empty();
        }

        return paths.get(0);
    }

    /** Get an array of interior angles for the area. An empty array is returned if there
     * are no boundary intersections (ie, it has only one boundary or no boundaries at all).
     *
     * <p>The order of the angles corresponds with the order of the boundaries returned
     * by {@link #getBoundaries()}: if {@code i} is an index into the boundaries list,
     * then {@code angles[i]} is the angle between boundaries {@code i} and {@code (i+1) % boundariesSize}.</p>
     * @return an array of interior angles for the area
     */
    public double[] getInteriorAngles() {
        final List<GreatArc> arcs = getBoundaryPath().getArcs();
        final int numSides = arcs.size();

        if (numSides < 2) {
            return new double[0];
        }

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

        GreatArc current;
        GreatArc next;
        for (int i = 0; i < numSides; ++i) {
            current = arcs.get(i);
            next = arcs.get((i + 1) % numSides);

            angles[i] = PlaneAngleRadians.PI - current.getCircle()
                    .angle(next.getCircle(), current.getEndPoint());
        }

        return angles;
    }

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

        if (numSides == 0) {
            return FULL_SIZE;
        } else if (numSides == 1) {
            return HALF_SIZE;
        } else {
            // use the extended version of Girard's theorem
            // https://en.wikipedia.org/wiki/Spherical_trigonometry#Girard's_theorem
            final double[] angles = getInteriorAngles();
            final double sum = Arrays.stream(angles).sum();

            return sum - ((angles.length - 2) * PlaneAngleRadians.PI);
        }
    }

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

    /** Return the weighted centroid vector of the area. The returned vector points in the direction of the
     * centroid point on the surface of the unit sphere with the length of the vector proportional to the
     * effective mass of the area at the centroid. By adding the weighted centroid vectors of multiple
     * convex areas, a single centroid can be computed for the combined area.
     * @return weighted centroid vector.
     * @see <a href="https://archive.org/details/centroidinertiat00broc">
     *  <em>The Centroid and Inertia Tensor for a Spherical Triangle</em> - John E. Brock</a>
     */
    Vector3D getWeightedCentroidVector() {
        final List<GreatArc> arcs = getBoundaries();
        final int numBoundaries = arcs.size();

        switch (numBoundaries) {
        case 0:
            // full space; no centroid
            return null;
        case 1:
            // hemisphere
            return computeHemisphereWeightedCentroidVector(arcs.get(0));
        case 2:
            // lune
            return computeLuneWeightedCentroidVector(arcs.get(0), arcs.get(1));
        default:
            // triangle or other convex polygon
            if (getBoundarySize() < TRIANGLE_FAN_CENTROID_COMPUTE_THRESHOLD) {
                return computeTriangleFanWeightedCentroidVector(arcs);
            }

            return computeArcPoleWeightedCentroidVector(arcs);
        }
    }

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

    /** Return a BSP tree representing the same region as this instance.
     */
    @Override
    public RegionBSPTree2S toTree() {
        return RegionBSPTree2S.from(getBoundaries(), true);
    }

    /** Return a new instance transformed by the argument.
     * @param transform transform to apply
     * @return a new instance transformed by the argument
     */
    public ConvexArea2S transform(final Transform<Point2S> transform) {
        return transformInternal(transform, this, GreatArc.class, ConvexArea2S::new);
    }

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

    /** Return an instance representing the full spherical 2D space.
     * @return an instance representing the full spherical 2D space.
     */
    public static ConvexArea2S full() {
        return FULL;
    }

    /** Construct a convex area by creating great circles between adjacent vertices. The vertices must be given
     * in a counter-clockwise around order the interior of the shape. If the area is intended to be closed, the
     * beginning point must be repeated at the end of the path.
     * @param vertices vertices to use to construct the area
     * @param precision precision context used to create new great circle instances
     * @return a convex area constructed using great circles between adjacent vertices
     * @see #fromVertexLoop(Collection, DoublePrecisionContext)
     */
    public static ConvexArea2S fromVertices(final Collection<Point2S> vertices,
            final DoublePrecisionContext precision) {
        return fromVertices(vertices, false, precision);
    }

    /** Construct a convex area by creating great circles between adjacent vertices. An implicit great circle is
     * created between the last vertex given and the first one, if needed. The vertices must be given in a
     * counter-clockwise around order the interior of the shape.
     * @param vertices vertices to use to construct the area
     * @param precision precision context used to create new great circles instances
     * @return a convex area constructed using great circles between adjacent vertices
     * @see #fromVertices(Collection, DoublePrecisionContext)
     */
    public static ConvexArea2S fromVertexLoop(final Collection<Point2S> vertices,
            final DoublePrecisionContext precision) {
        return fromVertices(vertices, true, precision);
    }

    /** Construct a convex area from great circles between adjacent vertices.
     * @param vertices vertices to use to construct the area
     * @param close if true, an additional great circle will be created between the last and first vertex
     * @param precision precision context used to create new great circle instances
     * @return a convex area constructed using great circles between adjacent vertices
     */
    public static ConvexArea2S fromVertices(final Collection<Point2S> vertices, final boolean close,
            final DoublePrecisionContext precision) {

        if (vertices.isEmpty()) {
            return full();
        }

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

        Point2S first = null;
        Point2S prev = null;
        Point2S cur = null;

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

            if (first == null) {
                first = cur;
            }

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

            prev = cur;
        }

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

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

        return fromBounds(circles);
    }

    /** Construct a convex area from an arc path. The area represents the intersection of all of the negative
     * half-spaces of the great circles in the path. The boundaries of the returned area may therefore not match
     * the arcs in the path.
     * @param path path to construct the area from
     * @return a convex area constructed from the great circles in the given path
     */
    public static ConvexArea2S fromPath(final GreatArcPath path) {
        final List<GreatCircle> bounds = path.getArcs().stream()
            .map(GreatArc::getCircle)
            .collect(Collectors.toList());

        return fromBounds(bounds);
    }

    /** Create a convex area formed by the intersection of the negative half-spaces of the
     * given bounding great circles. The returned instance represents the area that is on the
     * minus side of all of the given circles. Note that this method does not support areas
     * of zero size (ie, infinitely thin areas or points.)
     * @param bounds great circles used to define the convex area
     * @return a new convex area instance representing the area on the minus side of all
     *      of the bounding great circles or an instance representing the full area if no
     *      circles are given
     * @throws IllegalArgumentException if the given set of bounding great circles do not form a convex area,
     *      meaning that there is no region that is on the minus side of all of the bounding circles.
     */
    public static ConvexArea2S fromBounds(final GreatCircle... bounds) {
        return fromBounds(Arrays.asList(bounds));
    }

    /** Create a convex area formed by the intersection of the negative half-spaces of the
     * given bounding great circles. The returned instance represents the area that is on the
     * minus side of all of the given circles. Note that this method does not support areas
     * of zero size (ie, infinitely thin areas or points.)
     * @param bounds great circles used to define the convex area
     * @return a new convex area instance representing the area on the minus side of all
     *      of the bounding great circles or an instance representing the full area if no
     *      circles are given
     * @throws IllegalArgumentException if the given set of bounding great circles do not form a convex area,
     *      meaning that there is no region that is on the minus side of all of the bounding circles.
     */
    public static ConvexArea2S fromBounds(final Iterable<GreatCircle> bounds) {
        final List<GreatArc> arcs = new ConvexRegionBoundaryBuilder<>(GreatArc.class).build(bounds);
        return arcs.isEmpty() ?
                full() :
                new ConvexArea2S(arcs);
    }

    /** Compute the weighted centroid vector for the hemisphere formed by the given arc.
     * @param arc arc defining the hemisphere
     * @return the weighted centroid vector for the hemisphere
     * @see #getWeightedCentroidVector()
     */
    private static Vector3D computeHemisphereWeightedCentroidVector(final GreatArc arc) {
        return arc.getCircle().getPole().withNorm(HALF_SIZE);
    }

    /** Compute the weighted centroid vector for the lune formed by the given arcs.
     * @param a first arc for the lune
     * @param b second arc for the lune
     * @return the weighted centroid vector for the lune
     * @see #getWeightedCentroidVector()
     */
    private static Vector3D computeLuneWeightedCentroidVector(final GreatArc a, final GreatArc b) {
        final Point2S aMid = a.getCentroid();
        final Point2S bMid = b.getCentroid();

        // compute the centroid vector as the exact center of the lune to avoid inaccurate
        // results with very small regions
        final Vector3D.Unit centroid = aMid.slerp(bMid, 0.5).getVector();

        // compute the weight using the reverse of the algorithm from computeArcPoleWeightedCentroidVector()
        final double weight =
            (a.getSize() * centroid.dot(a.getCircle().getPole())) +
            (b.getSize() * centroid.dot(b.getCircle().getPole()));

        return centroid.withNorm(weight);
    }

    /** Compute the weighted centroid vector for the triangle or polygon formed by the given arcs
     * by adding together the arc pole vectors multiplied by their respective arc lengths. This
     * algorithm is described in the paper <a href="https://archive.org/details/centroidinertiat00broc">
     * <em>The Centroid and Inertia Tensor for a Spherical Triangle</em></a> by John E Brock.
     *
     * <p>Note: This algorithm works well in general but is susceptible to floating point errors
     * on very small areas. In these cases, the computed centroid may not be in the expected location
     * and may even be outside of the area. The {@link #computeTriangleFanWeightedCentroidVector(List)}
     * method can produce more accurate results in these cases.</p>
     * @param arcs boundary arcs for the area
     * @return the weighted centroid vector for the area
     * @see #computeTriangleFanWeightedCentroidVector(List)
     */
    private static Vector3D computeArcPoleWeightedCentroidVector(final List<GreatArc> arcs) {
        Vector3D centroid = Vector3D.ZERO;

        Vector3D arcContribution;
        for (final GreatArc arc : arcs) {
            arcContribution = arc.getCircle().getPole().withNorm(arc.getSize());

            centroid = centroid.add(arcContribution);
        }

        return centroid;
    }

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

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

        Vector3D areaCentroid = Vector3D.ZERO;

        GreatArc arc;
        Point2S p1;
        Point2S p2;
        Vector3D.Unit v1;
        Vector3D.Unit v2;
        Vector3D.Unit triangleCentroid;
        double triangleCentroidLen;
        while (arcIt.hasNext()) {
            arc = arcIt.next();

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

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

                triangleCentroid = v0.add(v1).add(v2).normalize();
                triangleCentroidLen =
                        computeArcCentroidContribution(v0, v1, triangleCentroid) +
                        computeArcCentroidContribution(v1, v2, triangleCentroid) +
                        computeArcCentroidContribution(v2, v0, triangleCentroid);

                areaCentroid = areaCentroid.add(triangleCentroid.withNorm(triangleCentroidLen));
            }
        }

        return areaCentroid;
    }

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

        return arcLength * triangleCentroid.dot(planeNormal);
    }
}