/**
 * DTSweep.java
 *
 * Copyright 2014-2014 Michael Hoffer <info@michaelhoffer.de>. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without modification, are
 * permitted provided that the following conditions are met:
 *
 *    1. Redistributions of source code must retain the above copyright notice, this list of
 *       conditions and the following disclaimer.
 *
 *    2. Redistributions in binary form must reproduce the above copyright notice, this list
 *       of conditions and the following disclaimer in the documentation and/or other materials
 *       provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY Michael Hoffer <info@michaelhoffer.de> "AS IS" AND ANY EXPRESS OR IMPLIED
 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
 * FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL Michael Hoffer <info@michaelhoffer.de> OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 * The views and conclusions contained in the software and documentation are those of the
 * authors and should not be interpreted as representing official policies, either expressed
 * or implied, of Michael Hoffer <info@michaelhoffer.de>.
 */ 

package eu.mihosoft.jcsg.ext.org.poly2tri;
/* Poly2Tri
 * Copyright (c) 2009-2010, Poly2Tri Contributors
 * http://code.google.com/p/poly2tri/
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 *
 * * Redistributions of source code must retain the above copyright notice,
 *   this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of Poly2Tri nor the names of its contributors may be
 *   used to endorse or promote products derived from this software without specific
 *   prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

import static eu.mihosoft.jcsg.ext.org.poly2tri.TriangulationUtil.EPSILON;
import static eu.mihosoft.jcsg.ext.org.poly2tri.TriangulationUtil.inScanArea;
import static eu.mihosoft.jcsg.ext.org.poly2tri.TriangulationUtil.orient2d;
import static eu.mihosoft.jcsg.ext.org.poly2tri.TriangulationUtil.smartIncircle;

import java.util.List;

import eu.mihosoft.jcsg.ext.org.poly2tri.TriangulationUtil.Orientation;

import org.slf4j.Logger;
import org.slf4j.LoggerFactory;

/**
 * Sweep-line, Constrained Delauney Triangulation (CDT) See: Domiter, V. and
 * Zalik, B.(2008)'Sweep-line algorithm for constrained Delaunay triangulation',
 * International Journal of Geographical Information Science
 *
 * "FlipScan" Constrained Edge Algorithm invented by author of this code.
 *
 * Author: Thomas ???, thahlen@gmail.com
 */
class DTSweep {

    private final static Logger logger = LoggerFactory.getLogger(DTSweep.class);

    private final static double PI_div2 = Math.PI / 2;
    private final static double PI_3div4 = 3 * Math.PI / 4;

    public DTSweep() {
    }

    /**
     * Triangulate simple polygon with holes *
     */
    public static void triangulate(DTSweepContext tcx) {
        tcx.createAdvancingFront();

        sweep(tcx);

        if (tcx.getTriangulationMode() == TriangulationMode.POLYGON) {
            finalizationPolygon(tcx);
        } else {
            finalizationConvexHull(tcx);
        }

        tcx.done();
    }

    /**
     * Start sweeping the Y-sorted point set from bottom to top
     *
     * @param tcx
     */
    private static void sweep(DTSweepContext tcx) {
        List<TriangulationPoint> points;
        TriangulationPoint point;
        AdvancingFrontNode node;

        points = tcx.getPoints();

        for (int i = 1; i < points.size(); i++) {
            point = points.get(i);

            node = pointEvent(tcx, point);

            if (point.hasEdges()) {
                for (DTSweepConstraint e : point.getEdges()) {
                    if (tcx.isDebugEnabled()) {
                        tcx.getDebugContext().setActiveConstraint(e);
                    }
                    edgeEvent(tcx, e, node);
                }
            }
            tcx.update(null);
        }
    }

    /**
     * If this is a Delaunay Triangulation of a pointset we need to fill so the
     * triangle mesh gets a ConvexHull
     *
     * @param tcx
     */
    private static void finalizationConvexHull(DTSweepContext tcx) {
        AdvancingFrontNode n1, n2;
        DelaunayTriangle t1, t2;
        TriangulationPoint first, p1;

        n1 = tcx.aFront.head.next;
        n2 = n1.next;
        first = n1.point;

        turnAdvancingFrontConvex(tcx, n1, n2);

        // TODO: implement ConvexHull for lower right and left boundary
        // Lets remove triangles connected to the two "algorithm" points
        // XXX: When the first the nodes are points in a triangle we need to do a flip before 
        //      removing triangles or we will lose a valid triangle.
        //      Same for last three nodes!
        // !!! If I implement ConvexHull for lower right and left boundary this fix should not be 
        //     needed and the removed triangles will be added again by default
        n1 = tcx.aFront.tail.prev;
        if (n1.triangle.contains(n1.next.point) && n1.triangle.contains(n1.prev.point)) {
            t1 = n1.triangle.neighborAcross(n1.point);
            rotateTrianglePair(n1.triangle, n1.point, t1, t1.oppositePoint(n1.triangle, n1.point));
            tcx.mapTriangleToNodes(n1.triangle);
            tcx.mapTriangleToNodes(t1);
        }
        n1 = tcx.aFront.head.next;
        if (n1.triangle.contains(n1.prev.point) && n1.triangle.contains(n1.next.point)) {
            t1 = n1.triangle.neighborAcross(n1.point);
            rotateTrianglePair(n1.triangle, n1.point, t1, t1.oppositePoint(n1.triangle, n1.point));
            tcx.mapTriangleToNodes(n1.triangle);
            tcx.mapTriangleToNodes(t1);
        }

        // Lower right boundary 
        first = tcx.aFront.head.point;
        n2 = tcx.aFront.tail.prev;
        t1 = n2.triangle;
        p1 = n2.point;
        n2.triangle = null;
        do {
            tcx.removeFromList(t1);
            p1 = t1.pointCCW(p1);
            if (p1 == first) {
                break;
            }
            t2 = t1.neighborCCW(p1);
            t1.clear();
            t1 = t2;
        } while (true);

        // Lower left boundary
        first = tcx.aFront.head.next.point;
        p1 = t1.pointCW(tcx.aFront.head.point);
        t2 = t1.neighborCW(tcx.aFront.head.point);
        t1.clear();
        t1 = t2;
        while (p1 != first) {
            tcx.removeFromList(t1);
            p1 = t1.pointCCW(p1);
            t2 = t1.neighborCCW(p1);
            t1.clear();
            t1 = t2;
        }

        // Remove current head and tail node now that we have removed all triangles attached
        // to them. Then set new head and tail node points
        tcx.aFront.head = tcx.aFront.head.next;
        tcx.aFront.head.prev = null;
        tcx.aFront.tail = tcx.aFront.tail.prev;
        tcx.aFront.tail.next = null;

        tcx.finalizeTriangulation();
    }

    /**
     * We will traverse the entire advancing front and fill it to form a convex
     * hull.<br>
     */
    private static void turnAdvancingFrontConvex(DTSweepContext tcx,
            AdvancingFrontNode b,
            AdvancingFrontNode c) {
        AdvancingFrontNode first = b;
        while (c != tcx.aFront.tail) {
            if (tcx.isDebugEnabled()) {
                tcx.getDebugContext().setActiveNode(c);
            }

            if (orient2d(b.point, c.point, c.next.point) == Orientation.CCW) {
                // [b,c,d] Concave - fill around c
                fill(tcx, c);
                c = c.next;
            } else {
                // [b,c,d] Convex
                if (b != first && orient2d(b.prev.point, b.point, c.point) == Orientation.CCW) {
                    // [a,b,c] Concave - fill around b
                    fill(tcx, b);
                    b = b.prev;
                } else {
                    // [a,b,c] Convex - nothing to fill
                    b = c;
                    c = c.next;
                }
            }
        }
    }

    private static void finalizationPolygon(DTSweepContext tcx) {
        // Get an Internal triangle to start with
        DelaunayTriangle t = tcx.aFront.head.next.triangle;
        TriangulationPoint p = tcx.aFront.head.next.point;
        while (!t.getConstrainedEdgeCW(p)) {
            t = t.neighborCCW(p);
        }

        // Collect interior triangles constrained by edges
        tcx.meshClean(t);
    }

    /**
     * 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.
     *
     * @param tcx
     * @param point
     * @return
     */
    private static AdvancingFrontNode pointEvent(DTSweepContext tcx,
            TriangulationPoint point) {
        AdvancingFrontNode node, newNode;

        node = tcx.locateNode(point);
        if (tcx.isDebugEnabled()) {
            tcx.getDebugContext().setActiveNode(node);
        }
        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.getX() <= node.point.getX() + EPSILON) {
            fill(tcx, node);
        }
        tcx.addNode(newNode);

        fillAdvancingFront(tcx, newNode);
        return newNode;
    }

    /**
     * Creates a new front triangle and legalize it
     *
     * @param tcx
     * @param point
     * @param node
     * @return
     */
    private static AdvancingFrontNode newFrontTriangle(DTSweepContext tcx,
            TriangulationPoint point,
            AdvancingFrontNode node) {
        AdvancingFrontNode newNode;
        DelaunayTriangle triangle;

        triangle = new DelaunayTriangle(point, node.point, node.next.point);
        triangle.markNeighbor(node.triangle);
        tcx.addToList(triangle);

        newNode = new AdvancingFrontNode(point);
        newNode.next = node.next;
        newNode.prev = node;
        node.next.prev = newNode;
        node.next = newNode;

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

        if (tcx.isDebugEnabled()) {
            tcx.getDebugContext().setActiveNode(newNode);
        }

        if (!legalize(tcx, triangle)) {
            tcx.mapTriangleToNodes(triangle);
        }

        return newNode;
    }

    /**
     *
     *
     * @param tcx
     * @param edge
     * @param node
     */
    private static void edgeEvent(DTSweepContext tcx,
            DTSweepConstraint edge,
            AdvancingFrontNode node) {
        try {
            tcx.edgeEvent.constrainedEdge = edge;
            tcx.edgeEvent.right = edge.p.getX() > edge.q.getX();

            if (tcx.isDebugEnabled()) {
                tcx.getDebugContext().setPrimaryTriangle(node.triangle);
            }

            if (isEdgeSideOfTriangle(node.triangle, edge.p, edge.q)) {
                return;
            }

            // For now we will do all needed filling
            // TODO: integrate with flip process might give some better performance 
            //       but for now this avoid the issue with cases that needs both flips and fills
            fillEdgeEvent(tcx, edge, node);

            edgeEvent(tcx, edge.p, edge.q, node.triangle, edge.q);
        } catch (PointOnEdgeException e) {
            logger.warn("Skipping edge: {}", e.getMessage());
        }
    }

    private static void fillEdgeEvent(DTSweepContext tcx, DTSweepConstraint edge, AdvancingFrontNode node) {
        if (tcx.edgeEvent.right) {
            fillRightAboveEdgeEvent(tcx, edge, node);
        } else {
            fillLeftAboveEdgeEvent(tcx, edge, node);
        }
    }

    private static void fillRightConcaveEdgeEvent(DTSweepContext tcx, DTSweepConstraint edge, AdvancingFrontNode node) {
        fill(tcx, node.next);
        if (node.next.point != edge.p) {
            // Next above or below edge?
            if (orient2d(edge.q, node.next.point, edge.p) == Orientation.CCW) {
                // Below
                if (orient2d(node.point, node.next.point, node.next.next.point) == Orientation.CCW) {
                    // Next is concave
                    fillRightConcaveEdgeEvent(tcx, edge, node);
                } else {
                    // Next is convex
                }
            }
        }
    }

    private static void fillRightConvexEdgeEvent(DTSweepContext tcx, DTSweepConstraint edge, AdvancingFrontNode node) {
        // Next concave or convex?
        if (orient2d(node.next.point, node.next.next.point, node.next.next.next.point) == Orientation.CCW) {
            // Concave
            fillRightConcaveEdgeEvent(tcx, edge, node.next);
        } else {
            // Convex
            // Next above or below edge?
            if (orient2d(edge.q, node.next.next.point, edge.p) == Orientation.CCW) {
                // Below
                fillRightConvexEdgeEvent(tcx, edge, node.next);
            } else {
                // Above
            }
        }
    }

    private static void fillRightBelowEdgeEvent(DTSweepContext tcx, DTSweepConstraint edge, AdvancingFrontNode node) {
        if (tcx.isDebugEnabled()) {
            tcx.getDebugContext().setActiveNode(node);
        }
        if (node.point.getX() < edge.p.getX()) // needed?
        {
            if (orient2d(node.point, node.next.point, node.next.next.point) == Orientation.CCW) {
                // Concave 
                fillRightConcaveEdgeEvent(tcx, edge, node);
            } else {
                // Convex
                fillRightConvexEdgeEvent(tcx, edge, node);
                // Retry this one
                fillRightBelowEdgeEvent(tcx, edge, node);
            }

        }
    }

    private static void fillRightAboveEdgeEvent(DTSweepContext tcx, DTSweepConstraint edge, AdvancingFrontNode node) {
        while (node.next.point.getX() < edge.p.getX()) {
            if (tcx.isDebugEnabled()) {
                tcx.getDebugContext().setActiveNode(node);
            }
            // Check if next node is below the edge
            Orientation o1 = orient2d(edge.q, node.next.point, edge.p);
            if (o1 == Orientation.CCW) {
                fillRightBelowEdgeEvent(tcx, edge, node);
            } else {
                node = node.next;
            }
        }
    }

    private static void fillLeftConvexEdgeEvent(DTSweepContext tcx, DTSweepConstraint edge, AdvancingFrontNode node) {
        // Next concave or convex?
        if (orient2d(node.prev.point, node.prev.prev.point, node.prev.prev.prev.point) == Orientation.CW) {
            // Concave
            fillLeftConcaveEdgeEvent(tcx, edge, node.prev);
        } else {
            // Convex
            // Next above or below edge?
            if (orient2d(edge.q, node.prev.prev.point, edge.p) == Orientation.CW) {
                // Below
                fillLeftConvexEdgeEvent(tcx, edge, node.prev);
            } else {
                // Above
            }
        }
    }

    private static void fillLeftConcaveEdgeEvent(DTSweepContext tcx, DTSweepConstraint edge, AdvancingFrontNode node) {
        fill(tcx, node.prev);
        if (node.prev.point != edge.p) {
            // Next above or below edge?
            if (orient2d(edge.q, node.prev.point, edge.p) == Orientation.CW) {
                // Below
                if (orient2d(node.point, node.prev.point, node.prev.prev.point) == Orientation.CW) {
                    // Next is concave
                    fillLeftConcaveEdgeEvent(tcx, edge, node);
                } else {
                    // Next is convex
                }
            }
        }
    }

    private static void fillLeftBelowEdgeEvent(DTSweepContext tcx, DTSweepConstraint edge, AdvancingFrontNode node) {
        if (tcx.isDebugEnabled()) {
            tcx.getDebugContext().setActiveNode(node);
        }
        if (node.point.getX() > edge.p.getX()) {
            if (orient2d(node.point, node.prev.point, node.prev.prev.point) == Orientation.CW) {
                // Concave 
                fillLeftConcaveEdgeEvent(tcx, edge, node);
            } else {
                // Convex
                fillLeftConvexEdgeEvent(tcx, edge, node);
                // Retry this one
                fillLeftBelowEdgeEvent(tcx, edge, node);
            }

        }
    }

    private static void fillLeftAboveEdgeEvent(DTSweepContext tcx, DTSweepConstraint edge, AdvancingFrontNode node) {
        while (node.prev.point.getX() > edge.p.getX()) {
            if (tcx.isDebugEnabled()) {
                tcx.getDebugContext().setActiveNode(node);
            }
            // Check if next node is below the edge
            Orientation o1 = orient2d(edge.q, node.prev.point, edge.p);
            if (o1 == Orientation.CW) {
                fillLeftBelowEdgeEvent(tcx, edge, node);
            } else {
                node = node.prev;
            }
        }
    }

    private static boolean isEdgeSideOfTriangle(DelaunayTriangle triangle,
            TriangulationPoint ep,
            TriangulationPoint eq) {
        int index;
        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 edgeEvent(DTSweepContext tcx,
            TriangulationPoint ep,
            TriangulationPoint eq,
            DelaunayTriangle triangle,
            TriangulationPoint point) {
        TriangulationPoint p1, p2;

        if (tcx.isDebugEnabled()) {
            tcx.getDebugContext().setPrimaryTriangle(triangle);
        }

        if (isEdgeSideOfTriangle(triangle, ep, eq)) {
            return;
        }

        p1 = triangle.pointCCW(point);
        Orientation o1 = orient2d(eq, p1, ep);
        if (o1 == Orientation.Collinear) {
            if (triangle.contains(eq, p1)) {
                triangle.markConstrainedEdge(eq, p1);
                // We are modifying the constraint maybe it would be better to 
                // not change the given constraint and just keep a variable for the new constraint
                tcx.edgeEvent.constrainedEdge.q = p1;
                triangle = triangle.neighborAcross(point);
                edgeEvent(tcx, ep, p1, triangle, p1);
            } else {
                throw new PointOnEdgeException("EdgeEvent - Point on constrained edge not supported yet");
            }
            if (tcx.isDebugEnabled()) {
                logger.info("EdgeEvent - Point on constrained edge");
            }
            return;
        }

        p2 = triangle.pointCW(point);
        Orientation o2 = orient2d(eq, p2, ep);
        if (o2 == Orientation.Collinear) {
            if (triangle.contains(eq, p2)) {
                triangle.markConstrainedEdge(eq, p2);
                // We are modifying the constraint maybe it would be better to 
                // not change the given constraint and just keep a variable for the new constraint
                tcx.edgeEvent.constrainedEdge.q = p2;
                triangle = triangle.neighborAcross(point);
                edgeEvent(tcx, ep, p2, triangle, p2);
            } else {
                throw new PointOnEdgeException("EdgeEvent - Point on constrained edge not supported yet");
            }
            if (tcx.isDebugEnabled()) {
                logger.info("EdgeEvent - Point on constrained edge");
            }
            return;
        }

        if (o1 == o2) {
            // Need to decide if we are rotating CW or CCW to get to a triangle
            // that will cross edge
            if (o1 == Orientation.CW) {
                triangle = triangle.neighborCCW(point);
            } else {
                triangle = triangle.neighborCW(point);
            }
            edgeEvent(tcx, ep, eq, triangle, point);
        } else {
            // This triangle crosses constraint so lets flippin start!
            flipEdgeEvent(tcx, ep, eq, triangle, point);
        }
    }

    private static void flipEdgeEvent(DTSweepContext tcx,
            TriangulationPoint ep,
            TriangulationPoint eq,
            DelaunayTriangle t,
            TriangulationPoint p) {
        TriangulationPoint op, newP;
        DelaunayTriangle ot;
        boolean inScanArea;

        ot = t.neighborAcross(p);
        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 RuntimeException("[BUG:FIXME] FLIP failed due to missing triangle");
        }

        if (t.getConstrainedEdgeAcross(p)) {
            throw new RuntimeException("Intersecting Constraints");
        }

        if (tcx.isDebugEnabled()) {
            tcx.getDebugContext().setPrimaryTriangle(t);
            tcx.getDebugContext().setSecondaryTriangle(ot);
        } // TODO: remove

        inScanArea = 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) {
                    if (tcx.isDebugEnabled()) {
                        System.out.println("[FLIP] - constrained edge done");
                    } // TODO: remove                    
                    t.markConstrainedEdge(ep, eq);
                    ot.markConstrainedEdge(ep, eq);
                    legalize(tcx, t);
                    legalize(tcx, ot);
                } else {
                    if (tcx.isDebugEnabled()) {
                        System.out.println("[FLIP] - subedge done");
                    } // TODO: remove
                    // XXX: I think one of the triangles should be legalized here?                    
                }
            } else {
                if (tcx.isDebugEnabled()) {
                    System.out.println("[FLIP] - flipping and continuing with triangle still crossing edge");
                } // TODO: remove
                Orientation o = orient2d(eq, op, ep);
                t = nextFlipTriangle(tcx, o, t, ot, p, op);
                flipEdgeEvent(tcx, ep, eq, t, p);
            }
        } else {
            newP = nextFlipPoint(ep, eq, ot, op);
            flipScanEdgeEvent(tcx, ep, eq, t, ot, newP);
            edgeEvent(tcx, ep, eq, t, p);
        }
    }

    /**
     * When we need to traverse from one triangle to the next we need the point
     * in current triangle that is the opposite point to the next triangle.
     *
     * @param ep
     * @param eq
     * @param ot
     * @param op
     * @return
     */
    private static TriangulationPoint nextFlipPoint(TriangulationPoint ep,
            TriangulationPoint eq,
            DelaunayTriangle ot,
            TriangulationPoint op) {
        Orientation o2d = orient2d(eq, op, ep);
        if (o2d == Orientation.CW) {
            // Right
            return ot.pointCCW(op);
        } else if (o2d == Orientation.CCW) {
            // Left                
            return ot.pointCW(op);
        } else {
            // TODO: implement support for point on constraint edge
            throw new PointOnEdgeException("Point on constrained edge not supported yet");
        }
    }

    /**
     * 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
     *
     * @param tcx
     * @param o - should be the result of an orient2d( eq, op, ep )
     * @param t - triangle 1
     * @param ot - triangle 2
     * @param p - a point shared by both triangles
     * @param op - another point shared by both triangles
     * @return returns the triangle still intersecting the edge
     */
    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.dEdge[edgeIndex] = true;
            legalize(tcx, ot);
            ot.clearDelunayEdges();
            return t;
        }
        // t is not crossing edge after flip
        edgeIndex = t.edgeIndex(p, op);
        t.dEdge[edgeIndex] = true;
        legalize(tcx, t);
        t.clearDelunayEdges();
        return ot;
    }

    /**
     * Scan part of the FlipScan algorithm<br>
     * When a triangle pair isn't flippable we will scan for the next point that
     * is inside the flip triangle scan area. When found we generate a new
     * flipEdgeEvent
     *
     * @param tcx
     * @param ep - last point on the edge we are traversing
     * @param eq - first point on the edge we are traversing
     * @param flipTriangle - the current triangle sharing the point eq with edge
     * @param t
     * @param p
     */
    private static void flipScanEdgeEvent(DTSweepContext tcx,
            TriangulationPoint ep,
            TriangulationPoint eq,
            DelaunayTriangle flipTriangle,
            DelaunayTriangle t,
            TriangulationPoint p) {
        DelaunayTriangle ot;
        TriangulationPoint op, newP;
        boolean inScanArea;

        ot = t.neighborAcross(p);
        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 RuntimeException("[BUG:FIXME] FLIP failed due to missing triangle");
        }

        if (tcx.isDebugEnabled()) {
            System.out.println("[FLIP:SCAN] - scan next point"); // TODO: remove
            tcx.getDebugContext().setPrimaryTriangle(t);
            tcx.getDebugContext().setSecondaryTriangle(ot);
        }

        inScanArea = inScanArea(eq,
                flipTriangle.pointCCW(eq),
                flipTriangle.pointCW(eq),
                op);
        if (inScanArea) {
            // flip with new edge op->eq
            flipEdgeEvent(tcx, eq, op, ot, op);
            // TODO: Actually I just figured out that it should be possible to 
            //       improve this by getting the next ot and op before the the above 
            //       flip and continue the flipScanEdgeEvent here
            // set new ot and op here and loop back to inScanArea test
            // also need to set a new flipTriangle first
            // Turns out at first glance that this is somewhat complicated
            // so it will have to wait.
        } else {
            newP = nextFlipPoint(ep, eq, ot, op);
            flipScanEdgeEvent(tcx, ep, eq, flipTriangle, ot, newP);
        }
    }

    /**
     * Fills holes in the Advancing Front
     *
     *
     * @param tcx
     * @param n
     */
    private static void fillAdvancingFront(DTSweepContext tcx, AdvancingFrontNode n) {
        AdvancingFrontNode node;
        double angle;

        // Fill right holes
        node = n.next;
        while (node.hasNext()) {
            if (isLargeHole(node)) {
                break;
            }
            fill(tcx, node);
            node = node.next;
        }

        // Fill left holes
        node = n.prev;
        while (node.hasPrevious()) {
            if (isLargeHole(node)) {
                break;
            }
            fill(tcx, node);
            node = node.prev;
        }

        // Fill right basins
        if (n.hasNext() && n.next.hasNext()) {
            angle = basinAngle(n);
            if (angle < PI_3div4) {
                fillBasin(tcx, n);
            }
        }
    }

    /**
     * @param node
     * @return true if hole angle exceeds 90 degrees
     */
    private static boolean isLargeHole(AdvancingFrontNode node) {
        double angle = angle(node.point, node.next.point, node.prev.point);
        //XXX: don't see angle being in range [-pi/2,0] due to how advancing front works
//        return (angle > PI_div2) || (angle < -PI_div2); 
        return (angle > PI_div2) || (angle < 0);

        // ISSUE 48: http://code.google.com/p/poly2tri/issues/detail?id=48
        // TODO: Adding this fix suggested in issues 48 caused some 
        //       triangulations to fail so commented it out for now.
        //
        // Also haven't been able to produce a triangulation that gives the
        // problem described in issue 48.
//        AdvancingFrontNode nextNode = node.next;
//        AdvancingFrontNode prevNode = node.prev;
//        if( !AngleExceeds90Degrees(node.point, 
//                                   nextNode.point, 
//                                   prevNode.point))
//        {
//            return false;
//        }
//
//        // Check additional points on front.
//        AdvancingFrontNode next2Node = nextNode.next;
//        // "..Plus.." because only want angles on same side as point being added.
//        if(    (next2Node != null) 
//            && !AngleExceedsPlus90DegreesOrIsNegative(node.point, 
//                                                      next2Node.point, 
//                                                      prevNode.point))
//        {
//            return false;
//        }
//
//        AdvancingFrontNode prev2Node = prevNode.prev;
//        // "..Plus.." because only want angles on same side as point being added.
//        if(    (prev2Node != null) 
//            && !AngleExceedsPlus90DegreesOrIsNegative(node.point, 
//                                                      nextNode.point, 
//                                                      prev2Node.point))
//        {
//            return false;
//        }
//        return true;
    }

//    private static boolean AngleExceeds90Degrees
//    (
//        TriangulationPoint origin, 
//        TriangulationPoint pa, 
//        TriangulationPoint pb
//    )
//    {
//        double angle = angle(origin, pa, pb);
//        return (angle > PI_div2) || (angle < -PI_div2);
//    }
//
//
//    private static boolean AngleExceedsPlus90DegreesOrIsNegative
//    (
//        TriangulationPoint origin, 
//        TriangulationPoint pa, 
//        TriangulationPoint pb
//    )
//    {
//        double angle = angle(origin, pa, pb);
//        return (angle > PI_div2) || (angle < 0);
//    }
    /**
     * Fills a basin that has formed on the Advancing Front to the right of
     * given node.<br>
     * First we decide a left,bottom and right node that forms the boundaries of
     * the basin. Then we do a reqursive fill.
     *
     * @param tcx
     * @param node - starting node, this or next node will be left node
     */
    private static void fillBasin(DTSweepContext tcx, AdvancingFrontNode node) {
        if (orient2d(node.point, node.next.point, node.next.next.point) == Orientation.CCW) {
            tcx.basin.leftNode = node;
        } else {
            tcx.basin.leftNode = node.next;
        }

        // Find the bottom and right node
        tcx.basin.bottomNode = tcx.basin.leftNode;
        while (tcx.basin.bottomNode.hasNext()
                && tcx.basin.bottomNode.point.getY() >= tcx.basin.bottomNode.next.point.getY()) {
            tcx.basin.bottomNode = tcx.basin.bottomNode.next;
        }
        if (tcx.basin.bottomNode == tcx.basin.leftNode) {
            // No valid basin
            return;
        }

        tcx.basin.rightNode = tcx.basin.bottomNode;
        while (tcx.basin.rightNode.hasNext()
                && tcx.basin.rightNode.point.getY() < tcx.basin.rightNode.next.point.getY()) {
            tcx.basin.rightNode = tcx.basin.rightNode.next;
        }
        if (tcx.basin.rightNode == tcx.basin.bottomNode) {
            // No valid basins
            return;
        }

        tcx.basin.width = tcx.basin.rightNode.getPoint().getX() - tcx.basin.leftNode.getPoint().getX();
        tcx.basin.leftHighest = tcx.basin.leftNode.getPoint().getY() > tcx.basin.rightNode.getPoint().getY();

        fillBasinReq(tcx, tcx.basin.bottomNode);
    }

    /**
     * Recursive algorithm to fill a Basin with triangles
     *
     * @param tcx
     * @param node - bottomNode
     * @param cnt - counter used to alternate on even and odd numbers
     */
    private static void fillBasinReq(DTSweepContext tcx, AdvancingFrontNode node) {
        // if shallow stop filling
        if (isShallow(tcx, node)) {
            return;
        }

        fill(tcx, node);
        if (node.prev == tcx.basin.leftNode && node.next == tcx.basin.rightNode) {
            return;
        } else if (node.prev == tcx.basin.leftNode) {
            Orientation o = orient2d(node.point, node.next.point, node.next.next.point);
            if (o == Orientation.CW) {
                return;
            }
            node = node.next;
        } else if (node.next == tcx.basin.rightNode) {
            Orientation o = orient2d(node.point, node.prev.point, node.prev.prev.point);
            if (o == Orientation.CCW) {
                return;
            }
            node = node.prev;
        } else {
            // Continue with the neighbor node with lowest Y value
            if (node.prev.point.getY() < node.next.point.getY()) {
                node = node.prev;
            } else {
                node = node.next;
            }
        }
        fillBasinReq(tcx, node);
    }

    private static boolean isShallow(DTSweepContext tcx, AdvancingFrontNode node) {
        double height;

        if (tcx.basin.leftHighest) {
            height = tcx.basin.leftNode.getPoint().getY() - node.getPoint().getY();
        } else {
            height = tcx.basin.rightNode.getPoint().getY() - node.getPoint().getY();
        }
        if (tcx.basin.width > height) {
            return true;
        }
        return false;
    }

    /**
     *
     * @param node - middle node
     * @return the angle between p-a and p-b in range [-pi,pi]
     */
    private static double angle(TriangulationPoint p,
            TriangulationPoint a,
            TriangulationPoint b) {
        // XXX: do we really need a signed angle for holeAngle?
        //      could possible save some cycles here
        /* Complex plane
         * ab = cosA +i*sinA
         * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx)
         * atan2(y,x) computes the principal value of the argument function
         * applied to the complex number x+iy
         * Where x = ax*bx + ay*by
         *       y = ax*by - ay*bx
         */
        final double px = p.getX();
        final double py = p.getY();
        final double ax = a.getX() - px;
        final double ay = a.getY() - py;
        final double bx = b.getX() - px;
        final double by = b.getY() - py;
        return Math.atan2(ax * by - ay * bx, ax * bx + ay * by);
    }

    /**
     * The basin angle is decided against the horizontal line [1,0]
     */
    private static double basinAngle(AdvancingFrontNode node) {
        double ax = node.point.getX() - node.next.next.point.getX();
        double ay = node.point.getY() - node.next.next.point.getY();
        return Math.atan2(ay, ax);
    }

    /**
     * Adds a triangle to the advancing front to fill a hole.
     *
     * @param tcx
     * @param node - middle node, that is the bottom of the hole
     */
    private static void fill(DTSweepContext tcx, AdvancingFrontNode node) {
        DelaunayTriangle triangle = new DelaunayTriangle(node.prev.point,
                node.point,
                node.next.point);
        // TODO: should copy the cEdge value from neighbor triangles
        //       for now cEdge values are copied during the legalize 
        triangle.markNeighbor(node.prev.triangle);
        triangle.markNeighbor(node.triangle);
        tcx.addToList(triangle);

        // Update the advancing front
        node.prev.next = node.next;
        node.next.prev = node.prev;
        tcx.removeNode(node);

        // If it was legalized the triangle has already been mapped
        if (!legalize(tcx, triangle)) {
            tcx.mapTriangleToNodes(triangle);
        }
    }

    /**
     * Returns true if triangle was legalized
     */
    private static boolean legalize(DTSweepContext tcx,
            DelaunayTriangle t) {
        int oi;
        boolean inside;
        TriangulationPoint p, op;
        DelaunayTriangle ot;
        // To legalize a triangle we start by finding if any of the three edges
        // violate the Delaunay condition
        for (int i = 0; i < 3; i++) {
            // TODO: fix so that cEdge is always valid when creating new triangles then we can check it here
            //       instead of below with ot
            if (t.dEdge[i]) {
                continue;
            }
            ot = t.neighbors[i];
            if (ot != null) {
                p = t.points[i];
                op = ot.oppositePoint(t, p);
                oi = ot.index(op);
                // If this is a Constrained Edge or a Delaunay Edge(only during recursive legalization)
                // then we should not try to legalize
                if (ot.cEdge[oi] || ot.dEdge[oi]) {
                    t.cEdge[i] = ot.cEdge[oi]; // XXX: have no good way of setting this property when creating new triangles so lets set it here                     
                    continue;
                }
                inside = smartIncircle(p,
                        t.pointCCW(p),
                        t.pointCW(p),
                        op);
                if (inside) {
                    boolean notLegalized;

                    // Lets mark this shared edge as Delaunay 
                    t.dEdge[i] = true;
                    ot.dEdge[oi] = true;

                    // Lets rotate shared edge one vertex CW to legalize it
                    rotateTrianglePair(t, p, ot, op);

                    // We now got one valid Delaunay Edge shared by two triangles
                    // This gives us 4 new edges to check for Delaunay
                    // Make sure that triangle to node mapping is done only one time for a specific triangle
                    notLegalized = !legalize(tcx, t);
                    if (notLegalized) {
                        tcx.mapTriangleToNodes(t);
                    }
                    notLegalized = !legalize(tcx, ot);
                    if (notLegalized) {
                        tcx.mapTriangleToNodes(ot);
                    }

                    // Reset the Delaunay edges, since they only are valid Delaunay edges
                    // until we add a new triangle or point.
                    // XXX: need to think about this. Can these edges be tried after we 
                    //      return to previous recursive level?
                    t.dEdge[i] = false;
                    ot.dEdge[oi] = false;

                    // If triangle have been legalized no need to check the other edges since
                    // the recursive legalization will handles those so we can end here.
                    return true;
                }
            }
        }
        return false;
    }

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

        boolean ce1, ce2, ce3, ce4;
        ce1 = t.getConstrainedEdgeCCW(p);
        ce2 = t.getConstrainedEdgeCW(p);
        ce3 = ot.getConstrainedEdgeCCW(op);
        ce4 = ot.getConstrainedEdgeCW(op);

        boolean de1, de2, de3, de4;
        de1 = t.getDelunayEdgeCCW(p);
        de2 = t.getDelunayEdgeCW(p);
        de3 = ot.getDelunayEdgeCCW(op);
        de4 = ot.getDelunayEdgeCW(op);

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

        // Remap dEdge
        ot.setDelunayEdgeCCW(p, de1);
        t.setDelunayEdgeCW(p, de2);
        t.setDelunayEdgeCCW(op, de3);
        ot.setDelunayEdgeCW(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.clearNeighbors();
        ot.clearNeighbors();
        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);
    }
}
