/*******************************************************************************
+
+  LEDA 3.5
+
+  delaunay.h
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
/*{\Manpage {delaunay_triang} {} {Delaunay Triangulations} }*/

#include <LEDA/graph.h>

#ifdef delaunay_triang
#undef delaunay_triang
#endif

#if defined(USE_FLOAT_KERNEL)

#define delaunay_triang float_delaunay_triang

#include <LEDA/point.h>
#include <LEDA/segment.h>
#include <LEDA/line.h>
#include <LEDA/circle.h>

typedef double  RAT_TYPE;
typedef vector  VECTOR;
typedef point   POINT;
typedef segment SEGMENT;
typedef line    LINE;
typedef circle  CIRCLE;

#else


#include <LEDA/rat_point.h>
#include <LEDA/rat_segment.h>
#include <LEDA/rat_line.h>
#include <LEDA/rat_circle.h>

typedef rational    RAT_TYPE;
typedef rat_vector  VECTOR;
typedef rat_point   POINT;
typedef rat_segment SEGMENT;
typedef rat_line    LINE;
typedef rat_circle  CIRCLE;

#endif

/*{/Msubst 
POINT   rat_point
SEGMENT rat_segment
LINE    rat_line
CIRCLE  rat_circle
}*/



/*{\Mdefinition

An instance $T$ of data type |\Mname| is a planar embedded bidirected
graph (map) representing the {\em Delaunay Triangulation} of its vertex set.
The position of a vertex $v$ is given by |T.pos(v)| and we use 
$S = \{ |T.pos(v)| \mid v \in T \}$ to denote the underlying point set.
Each face of $T$ (except for the outer face) is  a triangle whose 
circumscribing circle does not contain any point of $S$ in its interior.
For every edge $e$, the sequence 
\[e, |T.face_cycle_succ(e)|, |T.face_cycle_succ(T.face_cycle_succ(e))|,\ldots\]
traces the boundary of the face to the left of $e$.
The edges of the outer face of $T$ form the convex hull of $S$; the trace of 
the convex hull is clockwise.
The subgraph obtained from $T$ by removing all diagonals of
co-circular quadrilaterals  is called the {\em Delaunay Diagram}
of $S$.

|\Mname| provides all {\em constant} graph operations, 
e.g., |T.reversal(e)| returns the reversal of edge $e$, |T.all_edges()|
returns the list of all edges of $T$, and |forall_edges(e,T)| iterates over
all edges of $T$. In addition, |\Mname| provides
operations for inserting and deleting points,
point location, nearest neighbor searches, and navigation in both
the triangulation and the diagram.

|\Mname|s are essentially objects of type |GRAPH<POINT,int>|, where the 
node information is the position of the node and the edge information is 
irrelevant. For a graph $G$ of type |GRAPH<POINT,int>| the function
|Is_Delaunay(G)| tests whether $G$ is a Delaunay triangulation of its vertices.

The data type |\Mname| is illustrated by |delaunay_demo| in the LEDA demo directory.

}*/


class delaunay_triang : public  GRAPH<POINT,int>
{

   enum edge_kind { DIAGRAM_EDGE = 0,  NON_DIAGRAM_EDGE = 1, HULL_EDGE = 2 };


   edge cur_edge;
   edge hull_edge;

   bool check;

   // for temporarily marking nodes

   int cur_mark;
   node_map<int> mark;


   void (*nnh)(node);
   void (*neh)(edge);
   void (*dnh)(node);
   void (*deh)(edge);

   void (*pre_mvh)(edge,node,node);
   void (*post_mvh)(edge,node,node);
   void (*mark_edge_handler)(edge);

   void post_new_node_handler(node v) { if (nnh) nnh(v); }
   void pre_del_node_handler(node v)  { if (dnh) dnh(v); }

   void post_new_edge_handler(edge e) { if (neh) neh(e); }
   void pre_del_edge_handler(edge e)  { if (deh) deh(e); }

   void pre_move_edge_handler(edge e, node v, node w) 
   { if (pre_mvh) pre_mvh(e,v,w); }

   void post_move_edge_handler(edge e, node v, node w) 
   { if (post_mvh) post_mvh(e,v,w); }


   void init_node_marks()
   { mark.init(*this,-1);
     cur_mark = 0;
    }

   void mark_node(node v) const 
   { ((node_map<int>&)mark)[v] = cur_mark; }

   void unmark_node(node v) const 
   { ((node_map<int>&)mark)[v] = cur_mark-1; }

   bool is_marked(node v)  const { return mark[v] == cur_mark; }

   void unmark_all_nodes() const 
   { ((int&)cur_mark)++; 
     if ( cur_mark == MAXINT)
        ((delaunay_triang*)this) -> init_node_marks();
   }


   void mark_edge(edge e, edge_kind k) 
   { assign(e,k); 
     if (mark_edge_handler) mark_edge_handler(e);
    }


   void init_hull();
   void make_delaunay(list<edge>&);
   void make_delaunay();

   bool simplify_face(node v, list<edge>& E);

   //void del_hull_node(node);

   void check_locate(edge answer,const POINT& p) const;

   node  new_node(POINT p) 
   { node v = GRAPH<POINT,int>::new_node(p); 
     mark[v] = -1;
     return v;
    }

   edge  new_edge(node v, node w) 
   { return GRAPH<POINT,int>::new_edge(v,w,0); }

   edge  new_edge(edge e, node w) 
   { return GRAPH<POINT,int>::new_edge(e,w,0,0); }

   void  del_node(node v)  { GRAPH<POINT,int>::del_node(v); }
   void  del_edge(edge e)  { GRAPH<POINT,int>::del_edge(e); }

  
   void  draw_voro_edge(const CIRCLE& c1, const CIRCLE& c2,
                        void (*draw_edge)(const POINT&,const POINT&),
                        void (*draw_ray) (const POINT&,const POINT&));


   void dfs(node s, const CIRCLE& C, list<node>& L) const;
   /* a procedure used in range_search(const CIRCLE& C) */
 

public:

/*{\Mcreation T }*/

  delaunay_triang();
/*{\Mcreate  creates an empty Delaunay Triangulation |\Mvar|.}*/ 

  delaunay_triang(const list<POINT>& S) { init(S); }
/*{\Mcreate  creates a Delaunay Triangulation |\Mvar| of the points in $S$. 
             If $S$ contains multiple occurrences of points only the last 
             occurence of each point is retained.}*/


  delaunay_triang(const GRAPH<POINT,int>& G);
/*{\Mcreate  initializes |\Mvar| with a copy of $G$.\\
             \precond |Is_Delaunay(G)| is true.}*/


  delaunay_triang(const delaunay_triang&);

  delaunay_triang& operator=(const delaunay_triang&);
 ~delaunay_triang() {}

// for testing
  edge get_cur_edge()  const { return cur_edge; }
  void set_cur_edge(edge e)  { cur_edge = e;    }
  void set_hull_edge(edge e) { hull_edge = e;   }


/*{\Moperations 2.5 4.5 }*/

  void  set_new_node_handler(void (*f)(node)=0) { nnh = f; }
  void  set_new_edge_handler(void (*f)(edge)=0) { neh = f; }
  void  set_del_node_handler(void (*f)(node)=0) { dnh = f; }
  void  set_del_edge_handler(void (*f)(edge)=0) { deh = f; }

  void set_pre_move_edge_handler(void (*f)(edge,node,node)=0) { pre_mvh = f; }
  void set_post_move_edge_handler(void (*f)(edge,node,node)=0) { post_mvh = f; }

  void set_mark_edge_handler(void (*f)(edge)=0) { mark_edge_handler = f; }


  void  init(const list<POINT>& L);
/*{\Mop makes |\Mvar|  a Delaunay Triangulation of the points in $S$.}*/


  POINT  pos(node v) const;
/*{\Mop  returns the position of node $v$. }*/

  POINT  pos_source(edge e) const;
/*{\Mop  returns the position of $source(e)$. }*/

  POINT  pos_target(edge e) const;
/*{\Mop  returns the position of $target(e)$. }*/

  SEGMENT seg(edge e) const;
/*{\Mop   returns the line segment corresponding to edge $e$ 
          (|SEGMENT(T.pos_source(e),T.pos_target(e))|). }*/

  LINE    supporting_line(edge e) const;
/*{\Mop   returns the supporting line of edge $e$ 
          (|LINE(T.pos_source(e),T.pos_target(e))|). }*/


  int     orientation(edge e, POINT p) const;
/*{\Mop   returns $orientation(T.seg(e),p)$.}*/

  int     dim() const;
/*{\Mop   returns $-1$ if $S$ is empty, returns $0$ if $S$ consists of only one point,  
          returns $1$ if $S$ consists of at least two points and all points in $S$ 
          are collinear, and returns $2$ otherwise.
}*/

  list<POINT>  points() const;
/*{\Mop  returns $S$. }*/

  node   insert(POINT p);
/*{\Mop  inserts point $p$ into |\Mvar| and returns the corresponding node. 
         More precisely, if there is already a node |v| in |\Mvar| positioned 
         at $p$ (i.e., |pos(v)| is equal to |p|) then |pos(v)| is changed to 
         |p| (i.e., |pos(v)| is made identical to |p|) and if there is no 
         such node then a new node $v$ with |pos(v) = p| is added to |\Mvar|. 
         In either case, $v$ is returned.}*/

  void   del(node v);
/*{\Mop  removes node $v$, i.e., makes |\Mvar| a Delaunay triangulation 
         of $S-\{|pos(v)|\}$. }*/

  edge   locate(POINT p) const;
/*{\Mop  returns an edge |e| of |T| that contains |p| or that borders the face
         that contains $p$. In the former case, the edge is either a hull edge 
         or an edge in the interior of the triangulation (i.e., neither the edge
         nor its reversal is a hull edge). In the latter case we have 
         |orientation(e,p) > 0| except if all points of |T| are collinear and 
         |p| lies on the induced line. In this case |target(e)| is visible from 
         |p|. The function returns |nil| if |T| has no edge.}*/

  node   lookup(POINT p) const;
/*{\Mop  if |\Mvar| contains a node $v$ with $pos(v) = p$ the result is $v$
         otherwise the result is |nil|. }*/

  node   nearest_neighbor(POINT p) const;
/*{\Mop  computes a node $v$ of |\Mvar| that is closest to $p$,
         i.e., $|dist(p,pos(v))| = \min\{\ |dist(p,pos(u))| \mid u\in T\ \}$.}*/ 
  list<node> range_search(const CIRCLE& C) const;
/*{\Mop returns the list of all nodes contained in the closure of disk $C$.\\
        \precond $C$ must be a proper circle (not a straight line).}*/


  list<node> range_search(const POINT& a, const POINT& b, const POINT& c) const;
/*{\Mop returns the list of all nodes contained in the closure of the triangle 
        $(a,b,c)$.\\
        \precond $a$, $b$, and $c$ must not be collinear.}*/

  list<node> range_search(const POINT& a, const POINT& b) const;
/*{\Mop returns the list of all nodes contained in the closure of the rectangle 
        with diagonal $(a,b)$. }*/



  edge   get_hull_edge() const;
/*{\Mop  returns an edge of the outer face  of |\Mvar| (i.e., an edge of the 
         convex hull). }*/

  bool   is_hull_edge(edge e) const;
/*{\Mop  decides whether $e$ belongs to the convex hull of |\Mvar|. }*/

  bool   is_diagram_edge(edge e) const;
/*{\Mop  decides whether $e$ is an edge of the Delaunay {\em diagram} 
         of |\Mvar|. }*/

  bool   IS_NON_DIAGRAM_EDGE(edge e) const;
/* computes wether $e$ is an edge of the Delaunay triangulation which is not an edge of the Delaunay diagram.*/


  edge   d_face_cycle_succ(edge e) const;
/*{\Mop  returns face cycle successor of $e$ in the Delaunay diagram 
         of |\Mvar|. \precond  $e$ belongs to the Delaunay diagram.}*/

  edge   d_face_cycle_pred(edge e) const;
/*{\Mop  returns the face cycle predecessor of $e$ in the Delaunay diagram 
         of |\Mvar|. \precond  $e$ belongs to the Delaunay diagram.}*/

  void   compute_voronoi(GRAPH<CIRCLE,POINT>& V);
/*{\Mop  computes the corresponding Voronoi diagram $V$. 
         Each node of $VD$ is labeled with its defining circle.
         Each edge is labeled with the site lying in the face to
         its left. }*/

list<edge> minimum_spanning_tree() const; 
/*{\Mop returns the list of edges of |\Mvar| that comprise a minimum spanning tree of |S|. }*/


  bool   empty() { return number_of_nodes() == 0; }
/*{\Mop  decides whether |\Mvar| is empty. }*/

  void   clear() { GRAPH<POINT,int>::clear(); cur_edge = hull_edge = nil; }
/*{\Mop  makes |\Mvar| empty. }*/


/*{\Mtext \headerline{Drawing Routines}
\setlength{\typewidth}{0cm}
\setlength{\callwidth}{3cm}
\computewidths
}*/

  void   draw_nodes(void (*draw_node)(const POINT&)); 
/*{\Mopl draws all nodes of |\Mvar| by calling 
         |draw_node(pos(v))| for every node $v$. }*/


  void   draw_edge(edge e,
                   void (*draw_diagram_edge)(const POINT&,const POINT&),
                   void (*draw_triang_edge) (const POINT&,const POINT&),
                   void (*draw_hull_edge)   (const POINT&,const POINT&));

/*{\Mopl draws edges |e| of |\Mvar| by calling 
         |draw_diagram_edge(pos_source(e),pos_target(e)| if |e| is a diagram edge,
         |draw_hull_edge(pos_source(e),pos_target(e)| if |e| is a hull edge
         and |draw_triang_edge(pos_source(e),pos_target(e)| if |e| is a
         non-diagram edge. }*/


  void   draw_edges( void (*draw_diagram_edge)(const POINT&,const POINT&),
                     void (*draw_triang_edge) (const POINT&,const POINT&),
                     void (*draw_hull_edge)   (const POINT&,const POINT&));
/*{\Mopl draws all edges of |\Mvar| by calling the corresponding draw 
         function. }*/


  void   draw_edges(const list<edge>& L,
                    void (*draw_edge)(const POINT&, const POINT&));
/*{\Mopl draws all edges in $L$ by calling 
         |draw_edge(pos_source(e),pos_target(e)| for every edge $e \in L$. }*/


  void   draw_voro_edges(void (*draw_edge)(const POINT&,const POINT&),
                         void (*draw_ray) (const POINT&,const POINT&));

/*{\Mopl draws the edges of the Vorononoi Diagram of |\Mvar| by calling ... }*/


  void   draw_hull(void (*draw_poly)(const list<POINT>&));
/*{\Mopl draws the convex hull of |\Mvar| by calling ... }*/


  void   draw_voro(const GRAPH<CIRCLE,POINT>&,
                   void (*draw_node)(const POINT&), 
                   void (*draw_edge)(const POINT&,const POINT&),
                   void (*draw_ray) (const POINT&,const POINT&));



void save_state(const POINT& p) const;
/*{\Xop saves the state of the data structure on 
files [[delaunay_error.graph]] and [[delaunay_error.aux]]. 
The graph part is written on the former file and
hull_edge and cur_edge are written on the latter file.}*/


bool check_state(const string& location) const;
/*{\Xop checks the current state of the data structure. 
If it finds an error it writes a diagnostic messages to |cerr|.}*/

void set_verify()
{ check = true; }
/*{\Xop enables checking mode.}*/

void unset_verify()
{ check = false; }
/*{\Xop disables checking mode.}*/

};





inline POINT delaunay_triang::pos(node v) const
{ return inf(v); }

inline POINT delaunay_triang::pos_source(edge e) const
{ return inf(source(e)); }

inline POINT delaunay_triang::pos_target(edge e) const
{ return inf(target(e)); }


inline SEGMENT delaunay_triang::seg(edge e) const           
{ return SEGMENT(inf(source(e)),inf(target(e))); }

inline LINE delaunay_triang::supporting_line(edge e) const           
{ return LINE(inf(source(e)),inf(target(e))); }

inline edge  delaunay_triang::get_hull_edge() const
{ return hull_edge; }


inline bool  delaunay_triang::is_diagram_edge(edge e) const
{ return !(inf(e) & NON_DIAGRAM_EDGE); }

inline bool  delaunay_triang::is_hull_edge(edge e) const
{ return inf(e) & HULL_EDGE; }


inline int delaunay_triang::orientation(edge e, POINT p) const
{ return ::orientation(pos(source(e)),pos(target(e)),p); }

inline int  delaunay_triang::dim() const
{ int n = number_of_nodes();
  if (n <= 1) 
     return n - 1;
  else
     return (is_hull_edge(reversal(hull_edge))) ? 1 : 2;
}

/*{\Mimplementation
The main ingredients for the implementation are Delaunay flipping, segment walking, and plane sweep.

The constructor |delaunay_triang(list<POINT> S)| first constructs a triangulation of |S| by sweeping and then makes the triangulation Delaunay by
a sequence of Delaunay flips. 

|Locate| walks through the triangulation along the segment from some fixed point of $T$ to the query point. |Insert| first locates the point, then updates the triangulation locally, and finally performs flips to reestablish the Delaunay property. |Delete| deletes the node, retriangulates the resulting face, and then performs flips. Nearest neighbor searching, circular range queries, and triangular range queries insert the query point into the triangulation, then perform an appropriate graph search on the triangulation, and finally remove
the query point.

All algorithms show good expected behaviour.

For details we refer the reader to the LEDA implementation report "Delaunay triangulations".

\begin{figure}

\caption{A Delaunay triangulation and a Voronoi diagram}
\end{figure}
}*/

