/*
    * 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.euclidean.threed.Vector3D;
  import org.apache.commons.numbers.angle.Angle;
  import org.apache.commons.numbers.core.Precision;
  
  /** 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 * Math.PI;
  
      /** Constant containing the area of half of the spherical space. */
      private static final double HALF_SIZE = Angle.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] = Math.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) * Math.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, Precision.DoubleEquivalence)
      */
     public static ConvexArea2S fromVertices(final Collection<Point2S> vertices,
             final Precision.DoubleEquivalence 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, Precision.DoubleEquivalence)
      */
     public static ConvexArea2S fromVertexLoop(final Collection<Point2S> vertices,
             final Precision.DoubleEquivalence 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 Precision.DoubleEquivalence 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) {
         final Vector3D.Sum centroid = Vector3D.Sum.create();
 
         for (final GreatArc arc : arcs) {
             centroid.addScaled(arc.getSize(), arc.getCircle().getPole());
         }
 
         return centroid.get();
     }
 
     /** 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();
 
         final Vector3D.Sum areaCentroid = Vector3D.Sum.create();
 
         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 = Vector3D.Sum.create()
                         .add(v0)
                         .add(v1)
                         .add(v2)
                         .get().normalize();
                 triangleCentroidLen =
                         computeArcCentroidContribution(v0, v1, triangleCentroid) +
                         computeArcCentroidContribution(v1, v2, triangleCentroid) +
                         computeArcCentroidContribution(v2, v0, triangleCentroid);
 
                 areaCentroid.addScaled(triangleCentroidLen, triangleCentroid);
             }
         }
 
         return areaCentroid.get();
     }
 
     /** 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);
     }
 }