using GeoAPI.Geometries;
using GisSharpBlog.NetTopologySuite.Geometries;

namespace GisSharpBlog.NetTopologySuite.Algorithm
{
    /// <summary> 
    /// Computes the centroid of an area point.
    /// Algorithm:
    /// Based on the usual algorithm for calculating
    /// the centroid as a weighted sum of the centroids
    /// of a decomposition of the area into (possibly overlapping) triangles.
    /// The algorithm has been extended to handle holes and multi-polygons.
    /// See <see href="http://www.faqs.org/faqs/graphics/algorithms-faq"/>
    /// for further details of the basic approach.
    /// </summary>
    public class CentroidArea
    {
        private ICoordinate basePt = null; // the point all triangles are based at
        private ICoordinate triangleCent3 = new Coordinate(); // temporary variable to hold centroid of triangle
        private double areasum2 = 0; // Partial area sum
        private readonly ICoordinate cg3 = new Coordinate(); // partial centroid sum

        /// <summary>
        /// 
        /// </summary>
        public CentroidArea() {}

        /// <summary>
        /// 
        /// </summary>
        public ICoordinate Centroid
        {
            get
            {
                ICoordinate cent = new Coordinate();
                cent.X = cg3.X/3/areasum2;
                cent.Y = cg3.Y/3/areasum2;
                return cent;
            }
        }

        /// <summary> 
        /// Adds the area defined by a Geometry to the centroid total.
        /// If the point has no area it does not contribute to the centroid.
        /// </summary>
        /// <param name="geom">The point to add.</param>
        public void Add(IGeometry geom)
        {
            if (geom is IPolygon)
            {
                IPolygon poly = (IPolygon) geom;
                BasePoint = poly.ExteriorRing.GetCoordinateN(0);
                Add(poly);
            }
            else if (geom is IGeometryCollection)
            {
                IGeometryCollection gc = (IGeometryCollection) geom;
                foreach (IGeometry geometry in gc.Geometries)
                {
                    Add(geometry);
                }
            }
        }

        /// <summary> 
        /// Adds the area defined by an array of
        /// coordinates.  The array must be a ring;
        /// i.e. end with the same coordinate as it starts with.
        /// </summary>
        /// <param name="ring">An array of Coordinates.</param>
        public void Add(ICoordinate[] ring)
        {
            BasePoint = ring[0];
            AddShell(ring);
        }

        /// <summary>
        /// 
        /// </summary>
        private ICoordinate BasePoint
        {
            get
            {
                return basePt;
            }
            set
            {
                if (basePt == null)
                {
                    basePt = value;
                }
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="poly"></param>
        private void Add(IPolygon poly)
        {
            AddShell(poly.ExteriorRing.Coordinates);
            foreach (ILineString ls in poly.InteriorRings)
            {
                AddHole(ls.Coordinates);
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="pts"></param>
        private void AddShell(ICoordinate[] pts)
        {
            bool isPositiveArea = !CGAlgorithms.IsCCW(pts);
            for (int i = 0; i < pts.Length - 1; i++)
            {
                AddTriangle(basePt, pts[i], pts[i + 1], isPositiveArea);
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="pts"></param>
        private void AddHole(ICoordinate[] pts)
        {
            bool isPositiveArea = CGAlgorithms.IsCCW(pts);
            for (int i = 0; i < pts.Length - 1; i++)
            {
                AddTriangle(basePt, pts[i], pts[i + 1], isPositiveArea);
            }
        }

        /// <summary>
        /// 
        /// </summary>
        /// <param name="p0"></param>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="isPositiveArea"></param>
        private void AddTriangle(ICoordinate p0, ICoordinate p1, ICoordinate p2, bool isPositiveArea)
        {
            double sign = (isPositiveArea) ? 1.0 : -1.0;
            Centroid3(p0, p1, p2, ref triangleCent3);
            double area2 = Area2(p0, p1, p2);
            cg3.X += sign*area2*triangleCent3.X;
            cg3.Y += sign*area2*triangleCent3.Y;
            areasum2 += sign*area2;
        }

        /// <summary> 
        /// Returns three times the centroid of the triangle p1-p2-p3.
        /// The factor of 3 is
        /// left in to permit division to be avoided until later.
        /// </summary>
        private static void Centroid3(ICoordinate p1, ICoordinate p2, ICoordinate p3, ref ICoordinate c)
        {
            c.X = p1.X + p2.X + p3.X;
            c.Y = p1.Y + p2.Y + p3.Y;
            return;
        }

        /// <summary>
        /// Returns twice the signed area of the triangle p1-p2-p3,
        /// positive if a,b,c are oriented ccw, and negative if cw.
        /// </summary>
        private static double Area2(ICoordinate p1, ICoordinate p2, ICoordinate p3)
        {
            return (p2.X - p1.X)*(p3.Y - p1.Y) - (p3.X - p1.X)*(p2.Y - p1.Y);
        }
    }
}