/// <summary>
 /// Give two points in any order. Will always be ordered so
 /// that q.y > p.y and q.x > p.x if same y value 
 /// </summary>
 public DTSweepConstraint(TriangulationPoint p1, TriangulationPoint p2)
 {
     P = p1;
     Q = p2;
     if (p1.y > p2.y)
     {
         Q = p1;
         P = p2;
     }
     else if (p1.y == p2.y)
     {
         if (p1.x > p2.x)
         {
             Q = p1;
             P = p2;
         }
         else if (p1.x == p2.x)
         {
             //                logger.info( "Failed to create constraint {}={}", p1, p2 );
             //                throw new DuplicatePointException( p1 + "=" + p2 );
             //                return;
         }
     }
     Q.AddEdge(this);
 }

		public static bool InScanArea(
			TriangulationPoint pa, TriangulationPoint pb, TriangulationPoint pc, TriangulationPoint pd)
		{
			double pdx = pd.X;
			double pdy = pd.Y;
			double adx = pa.X - pdx;
			double ady = pa.Y - pdy;
			double bdx = pb.X - pdx;
			double bdy = pb.Y - pdy;

			double adxbdy = adx * bdy;
			double bdxady = bdx * ady;
			double oabd = adxbdy - bdxady;
			//        oabd = orient2d(pa,pb,pd);
			if (oabd <= 0)
			{
				return false;
			}

			double cdx = pc.X - pdx;
			double cdy = pc.Y - pdy;

			double cdxady = cdx * ady;
			double adxcdy = adx * cdy;
			double ocad = cdxady - adxcdy;
			//      ocad = orient2d(pc,pa,pd);
			if (ocad <= 0)
			{
				return false;
			}
			return true;
		}

		/// <summary>
		/// Give two points in any order. Will always be ordered so
		/// that q.y > p.y and q.x > p.x if same y value 
		/// </summary>
		public DTSweepConstraint(TriangulationPoint p1, TriangulationPoint p2)
		{
			this.P = p1;
			this.Q = p2;
			if (p1.Y > p2.Y)
			{
				this.Q = p1;
				this.P = p2;
			}
			else if (p1.Y == p2.Y)
			{
				if (p1.X > p2.X)
				{
					this.Q = p1;
					this.P = p2;
				}
				else if (p1.X == p2.X)
				{
					//                logger.info( "Failed to create constraint {}={}", p1, p2 );
					//                throw new DuplicatePointException( p1 + "=" + p2 );
					//                return;
				}
			}
			this.Q.AddEdge(this);
		}

 public PointOnEdgeException( string message, TriangulationPoint a, TriangulationPoint b, TriangulationPoint c )
     : base(message)
 {
     A=a;
     B=b;
     C=c;
 }

		/// <summary>
		///   Requirements:
		/// 1. a,b and c form a triangle.
		/// 2. a and d is know to be on opposite side of bc
		/// <code>
		///                a
		///                +
		///               / \
		///              /   \
		///            b/     \c
		///            +-------+ 
		///           /    B    \  
		///          /           \ 
		/// </code>
		///    Facts:
		///  d has to be in area B to have a chance to be inside the circle formed by a,b and c
		///  d is outside B if orient2d(a,b,d) or orient2d(c,a,d) is CW
		///  This preknowledge gives us a way to optimize the incircle test
		/// </summary>
		/// <param name="pa">triangle point, opposite d</param>
		/// <param name="pb">triangle point</param>
		/// <param name="pc">triangle point</param>
		/// <param name="pd">point opposite a</param>
		/// <returns>true if d is inside circle, false if on circle edge</returns>
		public static bool SmartIncircle( TriangulationPoint pa, TriangulationPoint pb, TriangulationPoint pc, TriangulationPoint pd ) {
			double pdx = pd.X;
			double pdy = pd.Y;
			double adx = pa.X - pdx;
			double ady = pa.Y - pdy;
			double bdx = pb.X - pdx;
			double bdy = pb.Y - pdy;

			double adxbdy = adx * bdy;
			double bdxady = bdx * ady;
			double oabd = adxbdy - bdxady;
			//        oabd = orient2d(pa,pb,pd);
			if (oabd <= 0) return false;

			double cdx = pc.X - pdx;
			double cdy = pc.Y - pdy;

			double cdxady = cdx * ady;
			double adxcdy = adx * cdy;
			double ocad = cdxady - adxcdy;
			//      ocad = orient2d(pc,pa,pd);
			if (ocad <= 0) return false;

			double bdxcdy = bdx * cdy;
			double cdxbdy = cdx * bdy;

			double alift = adx * adx + ady * ady;
			double blift = bdx * bdx + bdy * bdy;
			double clift = cdx * cdx + cdy * cdy;

			double det = alift * (bdxcdy - cdxbdy) + blift * ocad + clift * oabd;

			return det > 0;
		}

		public PointOnEdgeException(string message, TriangulationPoint a, TriangulationPoint b, TriangulationPoint c)
			: base(message)
		{
			this.A = a;
			this.B = b;
			this.C = c;
		}

 /// Forumla to calculate signed area
 /// Positive if CCW
 /// Negative if CW
 /// 0 if collinear
 /// A[P1,P2,P3]  =  (x1*y2 - y1*x2) + (x2*y3 - y2*x3) + (x3*y1 - y3*x1)
 ///              =  (x1-x3)*(y2-y3) - (y1-y3)*(x2-x3)
 public static Orientation Orient2d( TriangulationPoint pa, TriangulationPoint pb, TriangulationPoint pc )
 {
     double detleft = (pa.X - pc.X) * (pb.Y - pc.Y);
     double detright = (pa.Y - pc.Y) * (pb.X - pc.X);
     double val = detleft - detright;
     if (val > -EPSILON && val < EPSILON) {
         return Orientation.Collinear;
     } else if (val > 0) {
         return Orientation.CCW;
     }
     return Orientation.CW;
 }

 //TODO: Port note - different implementation
 public int IndexCW(TriangulationPoint p)
 {
     int index = IndexOf(p);
     switch (index)
     {
         case 0:
             return 2;
         case 1:
             return 0;
         default:
             return 1;
     }
 }

        /*
        public static bool InScanArea(TriangulationPoint pa, TriangulationPoint pb, TriangulationPoint pc,
                                      TriangulationPoint pd)
        {
            double pdx = pd.X;
            double pdy = pd.Y;
            double adx = pa.X - pdx;
            double ady = pa.Y - pdy;
            double bdx = pb.X - pdx;
            double bdy = pb.Y - pdy;

            double adxbdy = adx*bdy;
            double bdxady = bdx*ady;
            double oabd = adxbdy - bdxady;
            //        oabd = orient2d(pa,pb,pd);
            if (oabd <= 0)
            {
                return false;
            }

            double cdx = pc.X - pdx;
            double cdy = pc.Y - pdy;

            double cdxady = cdx*ady;
            double adxcdy = adx*cdy;
            double ocad = cdxady - adxcdy;
            //      ocad = orient2d(pc,pa,pd);
            if (ocad <= 0)
            {
                return false;
            }
            return true;
        }
        */
        public static bool InScanArea(TriangulationPoint pa, TriangulationPoint pb, TriangulationPoint pc, TriangulationPoint pd)
        {
            double oadb = (pa.X - pb.X) * (pd.Y - pb.Y) - (pd.X - pb.X) * (pa.Y - pb.Y);
            if (oadb >= -EPSILON)
            {
                return false;
            }

            double oadc = (pa.X - pc.X) * (pd.Y - pc.Y) - (pd.X - pc.X) * (pa.Y - pc.Y);
            if (oadc <= EPSILON)
            {
                return false;
            }
            return true;
        }

 private static bool AngleExceeds90Degrees(TriangulationPoint origin, TriangulationPoint pa, TriangulationPoint pb)
 {
     double angle = Angle(origin, pa, pb);
     bool exceeds90Degrees = ((angle > PI_div2) || (angle < -PI_div2));
     return exceeds90Degrees;
 }

        /// <summary>
        /// Rotates a triangle pair one vertex CW
        ///       n2                    n2
        ///  P +-----+             P +-----+
        ///    | t  /|               |\  t |  
        ///    |   / |               | \   |
        ///  n1|  /  |n3           n1|  \  |n3
        ///    | /   |    after CW   |   \ |
        ///    |/ oT |               | oT \|
        ///    +-----+ oP            +-----+
        ///       n4                    n4
        /// </summary>
        private static void RotateTrianglePair(DelaunayTriangle t, TriangulationPoint p, DelaunayTriangle ot, TriangulationPoint op)
        {
            DelaunayTriangle n1 = t.NeighborCCW(p);
            DelaunayTriangle n2 = t.NeighborCW(p);
            DelaunayTriangle n3 = ot.NeighborCCW(op);
            DelaunayTriangle n4 = ot.NeighborCW(op);

            bool ce1 = t.GetConstrainedEdgeCCW(p);
            bool ce2 = t.GetConstrainedEdgeCW(p);
            bool ce3 = ot.GetConstrainedEdgeCCW(op);
            bool ce4 = ot.GetConstrainedEdgeCW(op);

            bool de1 = t.GetDelaunayEdgeCCW(p);
            bool de2 = t.GetDelaunayEdgeCW(p);
            bool de3 = ot.GetDelaunayEdgeCCW(op);
            bool de4 = ot.GetDelaunayEdgeCW(op);

            t.Legalize(p, op);
            ot.Legalize(op, p);

            // Remap dEdge
            ot.SetDelaunayEdgeCCW(p, de1);
            t.SetDelaunayEdgeCW(p, de2);
            t.SetDelaunayEdgeCCW(op, de3);
            ot.SetDelaunayEdgeCW(op, de4);

            // Remap cEdge
            ot.SetConstrainedEdgeCCW(p, ce1);
            t.SetConstrainedEdgeCW(p, ce2);
            t.SetConstrainedEdgeCCW(op, ce3);
            ot.SetConstrainedEdgeCW(op, ce4);

            // Remap neighbors
            // XXX: might optimize the markNeighbor by keeping track of
            //      what side should be assigned to what neighbor after the
            //      rotation. Now mark neighbor does lots of testing to find
            //      the right side.
            t.Neighbors.Clear();
            ot.Neighbors.Clear();
            if (n1 != null) ot.MarkNeighbor(n1);
            if (n2 != null) t.MarkNeighbor(n2);
            if (n3 != null) t.MarkNeighbor(n3);
            if (n4 != null) ot.MarkNeighbor(n4);
            t.MarkNeighbor(ot);
        }

        /// <summary>
        /// Find closes node to the left of the new point and
        /// create a new triangle. If needed new holes and basins
        /// will be filled to.
        /// </summary>
        private static AdvancingFrontNode PointEvent(DTSweepContext tcx, TriangulationPoint point)
        {
            AdvancingFrontNode node = tcx.LocateNode(point);
            AdvancingFrontNode newNode = NewFrontTriangle(tcx, point, node);

            // Only need to check +epsilon since point never have smaller
            // x value than node due to how we fetch nodes from the front
            if (point.X <= node.Point.X + TriangulationUtil.EPSILON)
            {
                Fill(tcx, node);
            }

            tcx.AddNode(newNode);

            FillAdvancingFront(tcx, newNode);
            return newNode;
        }

		public DelaunayTriangle( TriangulationPoint p1, TriangulationPoint p2, TriangulationPoint p3 )
		{
			points[0] = p1;
			points[1] = p2;
			points[2] = p3;
		}

 /// <summary>
 /// After a flip we have two triangles and know that only one will still be
 /// intersecting the edge. So decide which to contiune with and legalize the other
 /// </summary>
 /// <param name="tcx"></param>
 /// <param name="o">should be the result of an TriangulationUtil.orient2d( eq, op, ep )</param>
 /// <param name="t">triangle 1</param>
 /// <param name="ot">triangle 2</param>
 /// <param name="p">a point shared by both triangles</param>
 /// <param name="op">another point shared by both triangles</param>
 /// <returns>returns the triangle still intersecting the edge</returns>
 private static DelaunayTriangle NextFlipTriangle(DTSweepContext tcx, Orientation o, DelaunayTriangle t, DelaunayTriangle ot, TriangulationPoint p, TriangulationPoint op)
 {
     int edgeIndex;
     if (o == Orientation.CCW)
     {
         // ot is not crossing edge after flip
         edgeIndex = ot.EdgeIndex(p, op);
         ot.EdgeIsDelaunay[edgeIndex] = true;
         Legalize(tcx, ot);
         ot.EdgeIsDelaunay.Clear();
         return t;
     }
     // t is not crossing edge after flip
     edgeIndex = t.EdgeIndex(p, op);
     t.EdgeIsDelaunay[edgeIndex] = true;
     Legalize(tcx, t);
     t.EdgeIsDelaunay.Clear();
     return ot;
 }

        /// <summary>
        /// Creates a new front triangle and legalize it
        /// </summary>
        private static AdvancingFrontNode NewFrontTriangle(DTSweepContext tcx, TriangulationPoint point, AdvancingFrontNode node)
        {
            DelaunayTriangle triangle = new DelaunayTriangle(point, node.Point, node.Next.Point);
            triangle.MarkNeighbor(node.Triangle);
            tcx.Triangles.Add(triangle);

            AdvancingFrontNode newNode = new AdvancingFrontNode(point);
            newNode.Next = node.Next;
            newNode.Prev = node;
            node.Next.Prev = newNode;
            node.Next = newNode;

            tcx.AddNode(newNode); // XXX: BST

            if (!Legalize(tcx, triangle))
            {
                tcx.MapTriangleToNodes(triangle);
            }

            return newNode;
        }

 private static bool IsEdgeSideOfTriangle(DelaunayTriangle triangle, TriangulationPoint ep, TriangulationPoint eq)
 {
     int index = triangle.EdgeIndex(ep, eq);
     if (index != -1)
     {
         triangle.MarkConstrainedEdge(index);
         triangle = triangle.Neighbors[index];
         if (triangle != null)
         {
             triangle.MarkConstrainedEdge(ep, eq);
         }
         return true;
     }
     return false;
 }

        private static void FlipEdgeEvent(DTSweepContext tcx, TriangulationPoint ep, TriangulationPoint eq, DelaunayTriangle t, TriangulationPoint p)
        {
            DelaunayTriangle ot = t.NeighborAcross(p);
            TriangulationPoint op = ot.OppositePoint(t, p);

            if (ot == null)
            {
                // If we want to integrate the fillEdgeEvent do it here
                // With current implementation we should never get here
                throw new InvalidOperationException("[BUG:FIXME] FLIP failed due to missing triangle");
            }

            if (t.GetConstrainedEdgeAcross(p))
            {
                throw new Exception("Intersecting Constraints");
            }

            bool inScanArea = TriangulationUtil.InScanArea(p, t.PointCCW(p), t.PointCW(p), op);
            if (inScanArea)
            {
                // Lets rotate shared edge one vertex CW
                RotateTrianglePair(t, p, ot, op);
                tcx.MapTriangleToNodes(t);
                tcx.MapTriangleToNodes(ot);

                if (p == eq && op == ep)
                {
                    if (eq == tcx.EdgeEvent.ConstrainedEdge.Q
                        && ep == tcx.EdgeEvent.ConstrainedEdge.P)
                    {
                        t.MarkConstrainedEdge(ep, eq);
                        ot.MarkConstrainedEdge(ep, eq);
                        Legalize(tcx, t);
                        Legalize(tcx, ot);
                    }
                    else
                    {
                        // XXX: I think one of the triangles should be legalized here?
                    }
                }
                else
                {
                    // TODO: remove
                    Orientation o = TriangulationUtil.Orient2d(eq, op, ep);
                    t = NextFlipTriangle(tcx, o, t, ot, p, op);
                    FlipEdgeEvent(tcx, ep, eq, t, p);
                }
            }
            else
            {
                TriangulationPoint newP = NextFlipPoint(ep, eq, ot, op);
                FlipScanEdgeEvent(tcx, ep, eq, t, ot, newP);
                EdgeEvent(tcx, ep, eq, t, p);
            }
        }

 private static bool AngleExceedsPlus90DegreesOrIsNegative(TriangulationPoint origin, TriangulationPoint pa, TriangulationPoint pb)
 {
     double angle = Angle(origin, pa, pb);
     bool exceedsPlus90DegreesOrIsNegative = (angle > PI_div2) || (angle < 0);
     return exceedsPlus90DegreesOrIsNegative;
 }

 public void AddSteinerPoint(TriangulationPoint point)
 {
     if (_steinerPoints == null)
     {
         _steinerPoints = new List<TriangulationPoint>();
     }
     _steinerPoints.Add(point);
 }

		public int IndexOf( TriangulationPoint p )
		{
			int i = points.IndexOf( p );
			if( i == -1 ) throw new Exception( "Calling index with a point that doesn't exist in triangle" );
			return i;
		}