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

namespace GisSharpBlog.NetTopologySuite.Algorithm
{
    /// <summary>
    /// Specifies and implements various fundamental Computational Geometric algorithms.
    /// The algorithms supplied in this class are robust for double-precision floating point.
    /// </summary>
    public static class CGAlgorithms
    {
        /// <summary> 
        /// A value that indicates an orientation of clockwise, or a right turn.
        /// </summary>
        public const int Clockwise = -1;

        /// <summary> 
        /// A value that indicates an orientation of clockwise, or a right turn.
        /// </summary>
        public const int Right = Clockwise;

        /// <summary>
        /// A value that indicates an orientation of counterclockwise, or a left turn.
        /// </summary>
        public const int CounterClockwise = 1;

        /// <summary>
        /// A value that indicates an orientation of counterclockwise, or a left turn.
        /// </summary>
        public const int Left = CounterClockwise;

        /// <summary>
        /// A value that indicates an orientation of collinear, or no turn (straight).
        /// </summary>
        public const int Collinear = 0;

        /// <summary>
        /// A value that indicates an orientation of collinear, or no turn (straight).
        /// </summary>
        public const int Straight = Collinear;

        /// <summary> 
        /// Returns the index of the direction of the point <c>q</c>
        /// relative to a vector specified by <c>p1-p2</c>.
        /// </summary>
        /// <param name="p1">The origin point of the vector.</param>
        /// <param name="p2">The final point of the vector.</param>
        /// <param name="q">The point to compute the direction to.</param>
        /// <returns> 
        /// 1 if q is counter-clockwise (left) from p1-p2,
        /// -1 if q is clockwise (right) from p1-p2,
        /// 0 if q is collinear with p1-p2.
        /// </returns>
        public static int OrientationIndex(ICoordinate p1, ICoordinate p2, ICoordinate q)
        {
            // travelling along p1->p2, turn counter clockwise to get to q return 1,
            // travelling along p1->p2, turn clockwise to get to q return -1,
            // p1, p2 and q are colinear return 0.
            double dx1 = p2.X - p1.X;
            double dy1 = p2.Y - p1.Y;
            double dx2 = q.X - p2.X;
            double dy2 = q.Y - p2.Y;
            return RobustDeterminant.SignOfDet2x2(dx1, dy1, dx2, dy2);
        }

        /// <summary> 
        /// Test whether a point lies inside a ring.
        /// The ring may be oriented in either direction.
        /// If the point lies on the ring boundary the result of this method is unspecified.
        /// This algorithm does not attempt to first check the point against the envelope
        /// of the ring.
        /// </summary>
        /// <param name="p">Point to check for ring inclusion.</param>
        /// <param name="ring">Assumed to have first point identical to last point.</param>
        /// <returns><c>true</c> if p is inside ring.</returns>
        public static bool IsPointInRing(ICoordinate p, ICoordinate[] ring)
        {
            if (p == null)
            {
                throw new ArgumentNullException("p", "coordinate 'p' is null");
            }
            if (ring == null)
            {
                throw new ArgumentNullException("ring", "ring 'ring' is null");
            }

            int crossings = 0;
            /*
            *  For each segment l = (i-1, i), see if it crosses ray from test point in positive x direction.
            */
            for (int i = 1; i < ring.Length; i++)
            {
                int i1 = i - 1;
                ICoordinate p1 = ring[i];
                ICoordinate p2 = ring[i1];

                if (!IsValidCoordinate(p1) || !IsValidCoordinate(p2))
                {
                    continue;
                }

                if (((p1.Y > p.Y) && (p2.Y <= p.Y)) || ((p2.Y > p.Y) && (p1.Y <= p.Y)))
                {
                    double x1 = p1.X - p.X;
                    double y1 = p1.Y - p.Y;
                    double x2 = p2.X - p.X;
                    double y2 = p2.Y - p.Y;
                    /*
                    *  segment straddles x axis, so compute intersection.
                    */
                    double xInt = RobustDeterminant.SignOfDet2x2(x1, y1, x2, y2)/(y2 - y1);

                    /*
                    *  crosses ray if strictly positive intersection.
                    */
                    if (0.0 < xInt)
                    {
                        crossings++;
                    }
                }
            }

            /*
            *  p is inside if number of crossings is odd.
            */
            if ((crossings%2) == 1)
            {
                return true;
            }
            else
            {
                return false;
            }
        }

        /// <summary> 
        /// Test whether a point lies on the line segments defined by a
        /// list of coordinates.
        /// </summary>
        /// <param name="p"></param>
        /// <param name="pt"></param>
        /// <returns> 
        /// <c>true</c> true if
        /// the point is a vertex of the line or lies in the interior of a line
        /// segment in the linestring.
        /// </returns>
        public static bool IsOnLine(ICoordinate p, ICoordinate[] pt)
        {
            LineIntersector lineIntersector = new RobustLineIntersector();
            for (int i = 1; i < pt.Length; i++)
            {
                ICoordinate p0 = pt[i - 1];
                ICoordinate p1 = pt[i];

                if (!IsValidCoordinate(p0) || !IsValidCoordinate(p1))
                {
                    continue;
                }

                lineIntersector.ComputeIntersection((Coordinate) p, (Coordinate) p0, (Coordinate) p1);
                if (lineIntersector.HasIntersection)
                {
                    return true;
                }
            }
            return false;
        }

        /// <summary>
        /// Computes whether a ring defined by an array of <see cref="Coordinate" />s is oriented counter-clockwise.
        /// The list of points is assumed to have the first and last points equal.
        /// This will handle coordinate lists which contain repeated points.
        /// This algorithm is only guaranteed to work with valid rings.
        /// If the ring is invalid (e.g. self-crosses or touches),
        /// the computed result may not be correct.
        /// </summary>>
        /// <param name="ring"></param>
        /// <returns></returns>
        public static bool IsCCW(ICoordinate[] ring)
        {
            // # of points without closing endpoint
            int nPts = ring.Length - 1;

            // find highest point
            ICoordinate hiPt = ring[0];
            int hiIndex = 0;
            for (int i = 1; i <= nPts; i++)
            {
                ICoordinate p = ring[i];
                if (p.Y > hiPt.Y)
                {
                    hiPt = p;
                    hiIndex = i;
                }
            }

            if (!IsValidCoordinate(hiPt))
            {
                return false;
            }

            // find distinct point before highest point
            int iPrev = hiIndex;
            do
            {
                iPrev = iPrev - 1;
                if (iPrev < 0)
                {
                    iPrev = nPts;
                }
            } while (ring[iPrev].Equals2D(hiPt) && iPrev != hiIndex);

            // find distinct point after highest point
            int iNext = hiIndex;
            do
                iNext = (iNext + 1)%nPts; while (ring[iNext].Equals2D(hiPt) && iNext != hiIndex);

            ICoordinate prev = ring[iPrev];
            ICoordinate next = ring[iNext];

            /*
             * This check catches cases where the ring contains an A-B-A configuration of points.
             * This can happen if the ring does not contain 3 distinct points
             * (including the case where the input array has fewer than 4 elements),
             * or it contains coincident line segments.
             */
            if (prev.Equals2D(hiPt) || next.Equals2D(hiPt) || prev.Equals2D(next))
            {
                return false;
            }

            int disc = ComputeOrientation(prev, hiPt, next);

            /*
             *  If disc is exactly 0, lines are collinear.  There are two possible cases:
             *  (1) the lines lie along the x axis in opposite directions
             *  (2) the lines lie on top of one another
             *
             *  (1) is handled by checking if next is left of prev ==> CCW
             *  (2) will never happen if the ring is valid, so don't check for it
             *  (Might want to assert this)
             */
            bool isCCW = false;
            if (disc == 0)
            {
                // poly is CCW if prev x is right of next x
                isCCW = (prev.X > next.X);
            }
            else
            {
                // if area is positive, points are ordered CCW
                isCCW = (disc > 0);
            }
            return isCCW;
        }

        public static bool IsValidCoordinate(ICoordinate coordinate)
        {
            return !(double.IsInfinity(coordinate.X) || double.IsInfinity(coordinate.Y) || double.IsNaN(coordinate.X) || double.IsNaN(coordinate.Y));
        }

        /// <summary>
        /// Computes the orientation of a point q to the directed line segment p1-p2.
        /// The orientation of a point relative to a directed line segment indicates
        /// which way you turn to get to q after travelling from p1 to p2.
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="q"></param>
        /// <returns> 
        /// 1 if q is counter-clockwise from p1-p2,
        /// -1 if q is clockwise from p1-p2,
        /// 0 if q is collinear with p1-p2-
        /// </returns>
        public static int ComputeOrientation(ICoordinate p1, ICoordinate p2, ICoordinate q)
        {
            return OrientationIndex(p1, p2, q);
        }

        /// <summary> 
        /// Computes the distance from a point p to a line segment AB.
        /// Note: NON-ROBUST!
        /// </summary>
        /// <param name="p">The point to compute the distance for.</param>
        /// <param name="A">One point of the line.</param>
        /// <param name="B">Another point of the line (must be different to A).</param>
        /// <returns> The distance from p to line segment AB.</returns>
        public static double DistancePointLine(ICoordinate p, ICoordinate A, ICoordinate B)
        {
            // if start == end, then use pt distance
            if (A.Equals(B))
            {
                return p.Distance(A);
            }

            // otherwise use comp.graphics.algorithms Frequently Asked Questions method
            /*(1)     	      AC dot AB
                        r =   ---------
                              ||AB||^2
             
		                r has the following meaning:
		                r=0 Point = A
		                r=1 Point = B
		                r<0 Point is on the backward extension of AB
		                r>1 Point is on the forward extension of AB
		                0<r<1 Point is interior to AB
	        */

            double r = ((p.X - A.X)*(B.X - A.X) + (p.Y - A.Y)*(B.Y - A.Y))
                       /
                       ((B.X - A.X)*(B.X - A.X) + (B.Y - A.Y)*(B.Y - A.Y));

            if (r <= 0.0)
            {
                return p.Distance(A);
            }
            if (r >= 1.0)
            {
                return p.Distance(B);
            }

            /*(2)
		                    (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
		                s = -----------------------------
		             	                Curve^2

		                Then the distance from C to Point = |s|*Curve.
	        */

            double s = ((A.Y - p.Y)*(B.X - A.X) - (A.X - p.X)*(B.Y - A.Y))
                       /
                       ((B.X - A.X)*(B.X - A.X) + (B.Y - A.Y)*(B.Y - A.Y));

            return Math.Abs(s)*Math.Sqrt(((B.X - A.X)*(B.X - A.X) + (B.Y - A.Y)*(B.Y - A.Y)));
        }

        /// <summary> 
        /// Computes the perpendicular distance from a point p
        /// to the (infinite) line containing the points AB
        /// </summary>
        /// <param name="p">The point to compute the distance for.</param>
        /// <param name="A">One point of the line.</param>
        /// <param name="B">Another point of the line (must be different to A).</param>
        /// <returns>The perpendicular distance from p to line AB.</returns>
        public static double DistancePointLinePerpendicular(ICoordinate p, ICoordinate A, ICoordinate B)
        {
            // use comp.graphics.algorithms Frequently Asked Questions method
            /*(2)
                            (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
                        s = -----------------------------
                                         Curve^2

                        Then the distance from C to Point = |s|*Curve.
            */

            double s = ((A.Y - p.Y)*(B.X - A.X) - (A.X - p.X)*(B.Y - A.Y))
                       /
                       ((B.X - A.X)*(B.X - A.X) + (B.Y - A.Y)*(B.Y - A.Y));

            return Math.Abs(s)*Math.Sqrt(((B.X - A.X)*(B.X - A.X) + (B.Y - A.Y)*(B.Y - A.Y)));
        }

        /// <summary> 
        /// Computes the distance from a line segment AB to a line segment CD.
        /// Note: NON-ROBUST!
        /// </summary>
        /// <param name="A">A point of one line.</param>
        /// <param name="B">The second point of the line (must be different to A).</param>
        /// <param name="C">One point of the line.</param>
        /// <param name="D">Another point of the line (must be different to A).</param>
        /// <returns>The distance from line segment AB to line segment CD.</returns>
        public static double DistanceLineLine(ICoordinate A, ICoordinate B, ICoordinate C, ICoordinate D)
        {
            // check for zero-length segments
            if (A.Equals(B))
            {
                return DistancePointLine(A, C, D);
            }
            if (C.Equals(D))
            {
                return DistancePointLine(D, A, B);
            }

            // AB and CD are line segments
            /* from comp.graphics.algo

	            Solving the above for r and s yields
				            (Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)
	                    r = ----------------------------- (eqn 1)
				            (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)

		 	                (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
		                s = ----------------------------- (eqn 2)
			                (Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
	            Let Point be the position vector of the intersection point, then
		            Point=A+r(B-A) or
		            Px=Ax+r(Bx-Ax)
		            Py=Ay+r(By-Ay)
	            By examining the values of r & s, you can also determine some other
                limiting conditions:
		            If 0<=r<=1 & 0<=s<=1, intersection exists
		            r<0 or r>1 or s<0 or s>1 line segments do not intersect
		            If the denominator in eqn 1 is zero, AB & CD are parallel
		            If the numerator in eqn 1 is also zero, AB & CD are collinear.

	        */
            double r_top = (A.Y - C.Y)*(D.X - C.X) - (A.X - C.X)*(D.Y - C.Y);
            double r_bot = (B.X - A.X)*(D.Y - C.Y) - (B.Y - A.Y)*(D.X - C.X);

            double s_top = (A.Y - C.Y)*(B.X - A.X) - (A.X - C.X)*(B.Y - A.Y);
            double s_bot = (B.X - A.X)*(D.Y - C.Y) - (B.Y - A.Y)*(D.X - C.X);

            if ((r_bot == 0) || (s_bot == 0))
            {
                return Math.Min(DistancePointLine(A, C, D),
                                Math.Min(DistancePointLine(B, C, D),
                                         Math.Min(DistancePointLine(C, A, B),
                                                  DistancePointLine(D, A, B))));
            }

            double s = s_top/s_bot;
            double r = r_top/r_bot;

            if ((r < 0) || (r > 1) || (s < 0) || (s > 1))
            {
                //no intersection
                return Math.Min(DistancePointLine(A, C, D),
                                Math.Min(DistancePointLine(B, C, D),
                                         Math.Min(DistancePointLine(C, A, B),
                                                  DistancePointLine(D, A, B))));
            }

            return 0.0; //intersection exists
        }

        /// <summary>
        /// Returns the signed area for a ring.  The area is positive ifthe ring is oriented CW.
        /// </summary>
        /// <param name="ring"></param>
        /// <returns></returns>
        public static double SignedArea(ICoordinate[] ring)
        {
            if (ring.Length < 3)
            {
                return 0.0;
            }

            double sum = 0.0;
            for (int i = 0; i < ring.Length - 1; i++)
            {
                double bx = ring[i].X;
                double by = ring[i].Y;
                double cx = ring[i + 1].X;
                double cy = ring[i + 1].Y;
                sum += (bx + cx)*(cy - by);
            }
            return -sum/2.0;
        }

        /// <summary> 
        /// Computes the length of a linestring specified by a sequence of points.
        /// </summary>
        /// <param name="pts">The points specifying the linestring.</param>
        /// <returns>The length of the linestring.</returns>
        public static double Length(ICoordinateSequence pts)
        {
            if (pts.Count < 1)
            {
                return 0.0;
            }

            double sum = 0.0;
            for (int i = 1; i < pts.Count; i++)
            {
                sum += pts.GetCoordinate(i).Distance(pts.GetCoordinate(i - 1));
            }

            return sum;
        }
    }
}