#ifndef TDelaunay2D_hxx_seen
#define TDelaunay2D_hxx_seen

#include <set>
#include <queue>

#include "IMesh.hxx"

namespace COMET {
    class IDelaunay2D;
    class IDelaunayTriangle;
};

///  The Delaunay tree is a randomized geometric data structure computing the
///  Delaunay triangulation. This structure holds insertion and
///  queries. Deletions are not supported in this version.
///
/// \verbatim
///  Olivier.Devillers@sophia.inria.fr
///  2004 route des Lucioles,
///  BP 93 06902 Sophia Antipolis Cedex 
///  (33) 93 65 77 63
///  Fax (33) 93 65 76 43
///  http://www.inria.fr:/prisme/personnel/devillers/devillers.html
///
///  Reference:
///                                                              
///    J.-D. Boissonnat and M. Teillaud,
///    On the randomized construction of the Delaunay tree,
///    Theoret. Comput. Sci. 112:339--354, 1993.
///
///    O. Devillers, S. Meiser and M. Teillaud,
///    Fully dynamic Delaunay triangulation in logarithmic expected time 
///       per operation, 
///    Comput. Geom. Theory Appl. 2(2):55-80, 1992.
///
///    O. Devillers,
///    Robust and efficient implementation of the Delaunay tree,
///   INRIA Research Report 1619, 1992.
/// \endverbatim
class COMET::IDelaunay2D : public COMET::IMesh {
public:
    class ITriangleIterator 
        : public std::iterator<std::input_iterator_tag,IDelaunayTriangle> {
    public:
        ITriangleIterator();
        ITriangleIterator(const ITriangleIterator&);
        ITriangleIterator(IDelaunayTriangle*);
        ~ITriangleIterator();
        IDelaunayTriangle& operator *() const;
        IDelaunayTriangle* operator ->() const;
        ITriangleIterator& operator ++();
        bool operator == (ITriangleIterator& rhs);
        bool operator != (ITriangleIterator& rhs);
    private:
        friend class IDelaunay2D;
        void PackQueue(void);
        std::queue<IDelaunayTriangle*> fQueue;
        std::set<IDelaunayTriangle*> fSeen;
    };
    
    /// Create an empty tree.
    IDelaunay2D(void);
    virtual ~IDelaunay2D(void);

    /// Insert a new point into the tree.
    virtual IMeshPoint* AddPoint(IMeshPoint*);

    /// Close up the triangulation after the last point is added.  This will
    /// correctly connect the edges to the IMeshPoint objects so that the mesh
    /// is well defined.  If this is called a second time, it will reconnect
    /// any edges that have been removed from the mesh.  
    virtual void Close();

    /// Provide an ITriangleIterator through the triangles.
    ITriangleIterator BeginTriangles();

    /// Provide an end of the triangle tree marker.
    ITriangleIterator& EndTriangles();

    /// Called by the IDelaunayTriangle constructors to add themselves into
    /// the owned list
    void AddTriangle(IDelaunayTriangle* t) {
        fChildren.push_back(t);
    }

    /// Initialized the triangulation
    void Init();

    void Paint(Option_t *option="");

private :
    /// number of the current operation
    int fNumber;

    /// the root of the delaunay_tree
    IDelaunayTriangle* fRoot;

    /// Keep track of all of the triangles owned by this tree so they can be
    /// deleted at the end.  It's not really possible to walk the graph and
    /// delete triangles since it has loops.  This works around the loops.
    std::vector<IDelaunayTriangle*> fChildren;

    static ITriangleIterator fEnd;
};

/// A class to represent a single triangle in the Delaunay triangulation.  
class COMET::IDelaunayTriangle {
    friend class COMET::IDelaunay2D::ITriangleIterator;
public:
    class IDT_flag {
    public :
        IDT_flag (void) {fFlag = (int) 0;}
        void infinite(int i) {fFlag |= (int) i;}
        void last_finite(void) {fFlag |= 8;}
        /// Mark that a triangle has been split.  After a split, a triangle
        /// has daughters, and it's vertices are no longer part of the
        /// triangulation.  
        void Split(void) {fFlag |= 16;}
        /// Check if a triangle has been split.
        int HasSplit(void) {return fFlag & 16;}
        int is_infinite(void) {return fFlag & 7;}
        int is_last_finite(void) {return fFlag & 8;}
        int is_degenerate(void) { return fFlag & 32;}
        /// Mark a particular triangle as bogus.  This works around the
        /// problem that the first points can't all be colinear.
        void Bogus(void) {fFlag |= 64;}
        /// Check if a triangle is bogus.
        int IsBogus(void) {return fFlag & 64;}
        int GetValue(void) {return fFlag;}
    private :
        int fFlag;
    };

    class IDT_list {
    public :
        IDT_list(IDT_list* l, IDelaunayTriangle* k): fNext(l), fKey(k) { }
        ~IDT_list(void);
        
        /// Return the next element in the single linked list.
        IDT_list* GetNext(void) {return fNext;}
        
        /// Get the node associated with the list element.
        IDelaunayTriangle* GetNode(void) {return fKey;}
        
    private :
        IDT_list *fNext;
        IDelaunayTriangle *fKey;
    };
    
public:
    virtual ~IDelaunayTriangle();

    /// Create an empty triangle which doesn't have any points, and isn't
    /// connected to any other triangles.  This is used to create the root
    /// triangle at the stat of the process.
    IDelaunayTriangle(IDelaunay2D* owner);

    /// Create from the parent triangle by taking it's points.
    IDelaunayTriangle(IDelaunay2D* owner, IDelaunayTriangle* parent, int); 

    /// Create a neighbor triangle from the parent and the new point.
    IDelaunayTriangle(IDelaunay2D* owner,
                      IDelaunayTriangle* parent, IMeshPoint* point, int);

    /// Check if a point is inside the circumcircle of the current node and
    /// return true if it is.  If this is true, the triangle will need to be
    /// split.
    bool Conflict(IMeshPoint *);

    /// Return the node number.
    unsigned int GetNumber() {return fNumber;}

    /// Set the node number.
    void SetNumber(unsigned int n) {fNumber = n;}

    /// Get the neighbor nodes
    IDelaunayTriangle* GetNeighbor(int i) {
        if (i<0||i>2) throw;
        return fNeighbors[i];
    }

    /// Set the neighbor nodes
    void SetNeighbor(int i, IDelaunayTriangle* n) {
        if (i<0||i>2) throw;
        fNeighbors[i] = n;
    }

    /// Get the corner points of the triangle.
    IMeshPoint* GetVertex(int i) {
        if (i<0||i>2) throw;
        return fVertices[i];
    }

    /// Find the center and radius of the circle defined by the points of this
    /// triangle.  The Delaunay tree guarantees that there won't be any other
    /// points inside of this circle.
    void Circle(double& xc, double& yc, double& r);

    IDT_flag& GetFlag(void) {return fFlag;}

    IDelaunayTriangle* FindConflict(IMeshPoint* p);

    int ClockwiseNeighbor(IMeshPoint *p) {
        return ((p==fVertices[0]) ? 2 : ((p==fVertices[1]) ? 0 : 1)); 
    }

    int NeighborIndex(IDelaunayTriangle *n) { 
        return ((fNeighbors[0]==n) ? 0 : ((fNeighbors[1]==n) ? 1 : 2));
    }

private :
    /// the first vertex is the creator, that is finite
    /// except for the root and its neighbors
    IDT_flag fFlag;

    /// The index number of this node.
    int fNumber;

    /// The points for this node.
    IMeshPoint* fVertices[3];

    /// The three neighbors of this node.
    IDelaunayTriangle* fNeighbors[3];

    /// The daughers of this node.
    IDT_list* fDaughters;
};

std::ostream& operator << (std::ostream& s, COMET::IDelaunayTriangle& t);
#endif