/*******************************************************************************
+
+  LEDA 3.5
+
+  delaunay.c
+
+  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.
+ 
*******************************************************************************/
//-----------------------------------------------------------------------------
//
// Delaunay Triangulations
//
// K. Mehlhorn and S. N"aher (1996)
//
//-----------------------------------------------------------------------------

#include <LEDA/delaunay.h>

#include <LEDA/plane_alg.h>
#include <LEDA/stream.h>
#include <math.h>


#define IN_HALFSPACE(e,p)\
((pos_target(e) - pos_source(e)) * (p - pos_source(e)) >= 0 )


void delaunay_triang::save_state(const POINT& p) const
{ ofstream out("delaunay_error.aux");
  if (hull_edge != nil)
  { out << pos_source(cur_edge) << endl;
    out << pos_target(cur_edge) << endl;
    out << pos_source(hull_edge) << endl;
    out << pos_target(hull_edge) << endl;
    out << p << endl;
    out.flush();
  }

  write("delaunay_error.graph");
}


bool delaunay_triang::IS_NON_DIAGRAM_EDGE(edge e) const
{ edge r = reversal(e);
  
  // e1,e2,e3,e4: edges of quadriliteral with diagonal e
  edge e1 = face_cycle_succ(r);
  edge e3 = face_cycle_succ(e);

    // flip test
  POINT a = pos_source(e1);
  POINT b = pos_target(e1);
  POINT c = pos_source(e3);
  POINT d = pos_target(e3);

  if (left_turn(a,b,c) && left_turn(b,c,d) && left_turn(c,d,a)
      && left_turn(d,a,c) && cocircular(a,b,c,d) )
     return true;
  return false;
}

bool delaunay_triang::check_state(const string&) const
{ return true; }




delaunay_triang::delaunay_triang()
{ nnh = 0;
  neh = 0;
  dnh = 0;
  deh = 0; 
  pre_mvh = 0;
  post_mvh = 0;
  mark_edge_handler = 0;
  cur_edge = nil; 
  hull_edge = nil;
  init_node_marks();
  check = false;
 }

delaunay_triang& delaunay_triang::operator=(const delaunay_triang& T) 
{ nnh = 0;
  neh = 0;
  dnh = 0;
  deh = 0; 
  pre_mvh = 0;
  post_mvh = 0;
  mark_edge_handler = 0;
  GRAPH<POINT,int>::operator=(T); 
  init_hull(); 
  if (number_of_edges() > 0) 
     cur_edge = reversal(hull_edge); 
  else 
     cur_edge = nil;
  init_node_marks();
  check = false;
  return *this; 
} 
  

delaunay_triang::delaunay_triang(const GRAPH<POINT,int>& G) 
                                                   :GRAPH<POINT,int>(G)
{ nnh = 0;
  neh = 0;
  dnh = 0;
  deh = 0; 
  pre_mvh = 0;
  post_mvh = 0;
  mark_edge_handler = 0;
  edge e;
  forall_edges(e,*this) mark_edge(e,DIAGRAM_EDGE);
  init_hull(); 
  make_delaunay(); 
  if (number_of_edges() > 0) 
     cur_edge = reversal(hull_edge); 
  else 
     cur_edge = nil;
  init_node_marks();
  check = false;
 }

delaunay_triang::delaunay_triang(const delaunay_triang& T):GRAPH<POINT,int>(T)
{ nnh = 0;
  neh = 0;
  dnh = 0;
  deh = 0; 
  pre_mvh = 0;
  post_mvh = 0;
  mark_edge_handler = 0;
  init_hull(); 
  if (number_of_edges() > 0) 
     cur_edge = reversal(hull_edge); 
  else 
     cur_edge = nil;
  init_node_marks();
  check = false;
}


/*
void delaunay_triang::init(const list<POINT>& L) 
{ clear();
  hull_edge = cur_edge = nil;
  if (!L.empty())
  { DELAUNAY_TRIANG(L,*this);
    init_hull();
    cur_edge = reversal(hull_edge);
   }
  init_node_marks();
}
*/


void delaunay_triang::init(const list<POINT>& L0) 
{ 
  // construct a triangulation for the points in L0

  clear();

  if (L0.empty()) return;

  list<POINT> L = L0;

  L.sort();  

  // initialize graph with a single edge starting at the first point

  POINT last_p = L.pop();          // last visited point
  node  last_v = new_node(last_p); // last inserted node

  while (!L.empty() && last_p == L.head()) L.pop();

  if (L.empty()) return;

  last_p = L.pop();
  node v = new_node(last_p);
  edge x = new_edge(last_v,v);
  edge y = new_edge(v,last_v);
  set_reversal(x,y);
  last_v = v;


  // scan remaining points

  POINT p;
  forall(p,L) 
  { if (p == last_p) continue; 

    edge e =  last_adj_edge(last_v);

    last_v = new_node(p);
    last_p = p;

    // walk up to upper tangent
    do e = face_cycle_pred(e); while (orientation(e,p) > 0);

    // walk down to lower tangent and triangulate
    do { e = face_cycle_succ(e);
         edge x = new_edge(e,last_v);
         edge y = new_edge(last_v,source(e));
         set_reversal(x,y);
       } while (orientation(e,p) > 0);
   }

  // mark edges of convex hull

  hull_edge = last_edge();
  
  cur_edge = reversal(hull_edge);
 
  edge e0 = hull_edge;
  edge e  = e0;
  do { mark_edge(e,HULL_EDGE);
       e = face_cycle_succ(e);
     } while (e != e0); 

  make_delaunay();

  init_node_marks();
}



void delaunay_triang::init_hull()
{ 
  hull_edge = nil;

  edge e;
  forall_edges(e,*this)
    if (orientation(e,pos_target(face_cycle_succ(e))) <= 0) 
    { hull_edge = e;
      break;
     }

  if (hull_edge)
  { edge e = hull_edge;
    do { mark_edge(e,HULL_EDGE);
         e = face_cycle_succ(e);
       } 
    while (e != hull_edge);
  }
}


void delaunay_triang::make_delaunay()
{ list<edge> S = all_edges();
  make_delaunay(S);
}


void delaunay_triang::make_delaunay(list<edge>& S)
{
  // Tranforms graph into a Delaunay delaunay_triang by flipping edges.
  // Diagonals of co-circular convex quadriliterals are marked as 
  // NON_DIAGRAM_EDGE
  // We maintain a stack $S$ of edges containing diagonals which might
  // have to be flipped. 

  if (number_of_nodes() <= 3) return;

  while ( !S.empty() )
  { edge e = S.pop();
    edge r = reversal(e);

    if (is_hull_edge(e) || is_hull_edge(r)) continue;

    mark_edge(e,DIAGRAM_EDGE);
    mark_edge(r,DIAGRAM_EDGE);

    // e1,e2,e3,e4: edges of quadriliteral with diagonal e
    edge e1 = face_cycle_succ(r);
    edge e3 = face_cycle_succ(e);

    // flip test
    POINT a = pos_source(e1);
    POINT b = pos_target(e1);
    POINT c = pos_source(e3);
    POINT d = pos_target(e3);

    if (left_turn(b,d,a) && right_turn(b,d,c))
    { // the quadrilateral is convex
      int soc = side_of_circle(a,b,c,d);

      if (soc == 0) // co-circular quadirliteral(a,b,c,d) 
      { mark_edge(e,NON_DIAGRAM_EDGE);
        mark_edge(r,NON_DIAGRAM_EDGE);
       }
      if (soc > 0) // flip
      { edge e2 = face_cycle_succ(e1);
        edge e4 = face_cycle_succ(e3);
  
        S.push(e1); 
        S.push(e2); 
        S.push(e3); 
        S.push(e4); 
  
        // flip diagonal
        move_edge(e,e2,source(e4));
        move_edge(r,e4,source(e2));
      }
    }
  }

}



list<POINT> delaunay_triang::points() const
{ list <POINT> L;
  node v;
  forall_nodes(v,*this) L.append(pos(v));
  return L;
}


edge delaunay_triang::d_face_cycle_succ(edge e) const
{ e = reversal(e);
  do e = cyclic_adj_pred(e); 
     while (!is_diagram_edge(e));
  return e;
 }

edge delaunay_triang::d_face_cycle_pred(edge e) const
{ do e = cyclic_adj_succ(e); 
     while (!is_diagram_edge(e));
  return reversal(e);
 }


void delaunay_triang::check_locate(edge answer, const POINT& p) const
{ 
   if (answer == nil && dim() < 1) return;

  if ( seg(answer).contains(p) && ( is_hull_edge(answer)
       || (!is_hull_edge(answer) && !is_hull_edge(reversal(answer)) ))) return;

  if (orientation(answer,p) < 0) 
  { cerr << "orientation(" << seg(answer) << "," << p << ") < 0" << endl;
    goto error_in_locate;
   }

  if ( is_hull_edge(answer) && orientation(answer,p) > 0 ) return;

  if (dim() == 1)
  { // orientation = 0 and answer does not contain p; so beyond extreme point
    edge e = face_cycle_succ(answer);
    if ( e == reversal(answer) && !IN_HALFSPACE(e ,p) )
       return;
    else 
    { cerr << "dim = 1 error" << endl;
      goto error_in_locate;
     }
  }

  // dim == 2: answer must not be a hull edge and triangle must contain p 
  if ( orientation(answer,p) > 0 &&
       orientation(face_cycle_succ(answer), p) > 0 &&
       orientation(face_cycle_pred(answer), p) > 0 )
     return;
  else
    { cerr << "dim = 2 error" << endl;
      goto error_in_locate;
     }

error_in_locate: 

  cerr << endl;
  cerr << "An error occurred in delaunay_triang::locate(point)." << endl;
  cerr << "I save the situation into the two files " << endl;
  cerr << "     delaunay_error.graph" << endl;
  cerr << "and" << endl;
  cerr << "     delaunay_error.aux" << endl;
  cerr << endl;
  cerr << "Please send these files to leda@mpi-sb.mpg.de." << endl;
  cerr << endl;

  save_state(p);
  exit(1);
}
 

/* Kurt's locate  */

edge delaunay_triang::locate(POINT p) const
{ 
  if (number_of_edges() == 0) return nil;

  if (dim() == 1)
  { edge e0 = hull_edge;
    int orient = orientation(e0,p);
    if (orient != 0)
    { if (orient < 0) e0 = reversal(e0);
      if (check) check_locate(e0 ,p);
      return e0;
     }

    // p is collinear with the points in S. We walk 
       
    if ( !IN_HALFSPACE(e0,p) ) e0 = reversal(e0);

    // in the direction of e0. We know IN_HALFSPACE(e0,p) 
     
    edge e1 = face_cycle_succ(e0);
    while ( e1 != reversal(e0) && IN_HALFSPACE(e1,p) ) 
    { e0 = e1;  
      e1 = face_cycle_succ(e0); 
     }
     
    if (check) check_locate(e0 ,p);

    return e0;
  }
    
  // dimension is two
  
  edge e0 = cur_edge;
  if ( is_hull_edge(e0) ) e0 = reversal(e0);

  // e0 is not a hull edge

  if (p == pos_source(e0) ) return reversal(e0);
  
  int orient = orientation(e0,p);

  if (orient == 0) 
  { e0 = face_cycle_pred(e0);
    orient = orientation(e0,p);
  }

  if (orient < 0) e0 = reversal(e0);
 
  SEGMENT s(pos_source(e0),p);

  while ( true )     
  { 
    if (is_hull_edge(e0)) break;

    edge e1 = face_cycle_succ(e0);
    edge e2 = face_cycle_pred(e0);
    int d = ::orientation(s,pos_target(e1));
    if (d != 0)
    { e0 = reversal( (d > 0) ? e1 : e2 );
      int orient = orientation(e0,p);
      if ( orient > 0) continue;
      if ( orient < 0) e0 = reversal(e0); 
      // if == 0 I want to keep e0 (it might be a hull edge)
      break;
     }

    // d == 0

    e0 = reversal(e1);
    if ( orientation(e0,p) > 0) continue;

    e0 = reversal(e2);
    if ( orientation(e0,p) > 0) continue;

    // p lies on the interior of the triangle or on e1 or on e2
    if ( orientation(e1,p) == 0) 
      e0 = e1;   // p lies on e1
    else 
      e0 = e2;   // p lies on e2 or in the interior of the triangle

    // e0 is never a hull edge at this point, but its reversal might be.
    // If e0 contains p we reverse it.
    if ( orientation(e0,p) == 0 ) e0 = reversal(e0);

    break;
  }

  if (check) check_locate(e0 ,p);

  ((edge&)cur_edge) = e0;
  return e0;
}



/* Stefan's locate */

/*
edge delaunay_triang::locate(POINT p) const
{ 
  if (number_of_edges() == 0) return nil;
  edge x = cur_edge;
  edge e = reversal(x);
  int orient = orientation(e,p);

  if (orient > 0) 
  { x = e;
    e = reversal(x);
    orient = -orient;
   }

  if (number_of_edges() == 2) return e;


  SEGMENT s(pos_source(e),p);

  while (orient <= 0 && p != pos_target(x) && !is_hull_edge(x))
  { edge e1 = face_cycle_succ(x);
    edge e2 = face_cycle_pred(x);
    e = (::orientation(s,pos_target(e1)) >= 0) ? e1 : e2;
    orient = orientation(e,p);
    if (orient <= 0) x = reversal(e);
   }


  if (orient > 0)
  { edge e1 = face_cycle_succ(e);
    if (orientation(e1,p) == 0) 
       { x = e1;
         orient = 0;
        }
    else
     { edge e2 = face_cycle_pred(e);
       if (orientation(e2,p) == 0) 
       { x = e2;
         orient = 0;
        }
      }
   }


  if (orient == 0)
  { 
    if (is_hull_edge(reversal(x))) x = reversal(x);

    if (is_hull_edge(x))
    { if (!IN_HALFSPACE(x,p))
        do  x = face_cycle_pred(x);
        while (orientation(x,p) == 0 && !IN_HALFSPACE(x,p));
      else
       { edge next = face_cycle_succ(x);
         while (orientation(next,p) == 0 && IN_HALFSPACE(next,p)) 
         { x = next;
           next = face_cycle_succ(x);
          }
         if (orientation(next,p) >0 ) x = next;
        }
     }
  }

  if (check) check_locate(e0 ,p);

  ((edge&)cur_edge) = x;

  return x;
}
*/




node delaunay_triang::lookup(POINT p) const
{ if (number_of_nodes() == 1)
  { node v = first_node();
    return (pos(v) == p) ? v : nil;
   }
  edge e = locate(p);
  if (pos(source(e)) == p) return source(e);
  if (pos(target(e)) == p) return target(e);
  return nil;
}




node delaunay_triang::insert(POINT p)
{ if ( check ) save_state(p);

  node v;

  
if (number_of_nodes() == 0)  { v = new_node(p); 
if ( check && !check_state("delaunay_triang::insert") )
{ cerr << "The point inserted was " << p; 
  exit(1);
}
return v;

                                                                               }

if (number_of_nodes() == 1)
{ node w = first_node();
  if (p == pos(w)) 
  { assign(w,p);      //Kurt: das hast du vergessen
    v = w; 
if ( check && !check_state("delaunay_triang::insert") )
{ cerr << "The point inserted was " << p; 
  exit(1);
}
return v;

  }
  else
  { v = new_node(p);
    edge x = new_edge(v,w);
    edge y = new_edge(w,v);
    mark_edge(x,HULL_EDGE);
    mark_edge(y,HULL_EDGE);
    set_reversal(x,y);
    hull_edge = x;
    cur_edge = x;
    
if ( check && !check_state("delaunay_triang::insert") )
{ cerr << "The point inserted was " << p; 
  exit(1);
}
return v;

  }
}


  edge e0 = locate(p);
  bool p_on_e0 = seg(e0).contains(p);
  bool outer_face = is_hull_edge(e0);

  
if (p_on_e0)
{ if (p == pos_source(e0)) { assign(source(e0),p); return source(e0); }
  if (p == pos_target(e0)) { assign(target(e0),p); return target(e0); }
} 

if (dim() == 1 && orientation(e0,p) == 0)
{ v = new_node(p);
  edge x = new_edge(v,target(e0));
  edge y = new_edge(target(e0),v);
  mark_edge(x,HULL_EDGE);
  mark_edge(y,HULL_EDGE);
  set_reversal(x,y);

  if (p_on_e0)
  { x = new_edge(v,source(e0));  
    y = new_edge(source(e0),v);
    mark_edge(x,HULL_EDGE);
    mark_edge(y,HULL_EDGE);
    set_reversal(x,y);
    hull_edge = cur_edge = x;

    del_edge(reversal(e0));
    del_edge(e0);
   }

  
if ( check && !check_state("delaunay_triang::insert") )
{ cerr << "The point inserted was " << p; 
  exit(1);
}
return v;

}



  
  v  = new_node(p);
  edge e1 = e0;
  edge e2 = e0;
  list<edge> E;

  if (outer_face) //  move e1/e2 to compute upper/lower tangents
  { do e1 = face_cycle_pred(e1); while (orientation(e1,p) > 0);
    do e2 = face_cycle_succ(e2); while (orientation(e2,p) > 0);
   }

  // insert edges between v and target(e1) ... source(e2)
  edge e = e1;
  do { e = face_cycle_succ(e);
       edge x = new_edge(e,v);
       edge y = new_edge(v,source(e));
       set_reversal(x,y);
       mark_edge(e,DIAGRAM_EDGE);
       E.append(e);
       E.append(x);
     } while (e != e2);

  if (outer_face) // mark last visited and new edges as hull edges
  { mark_edge(face_cycle_succ(e1),HULL_EDGE);
    mark_edge(face_cycle_pred(e2),HULL_EDGE);
    mark_edge(e2,HULL_EDGE);
    hull_edge = e2;
    cur_edge = reversal(e2);  // Kurt: essential when dimension jumps
   }

  make_delaunay(E); // removes degenerated triangles 
  
if ( check && !check_state("delaunay_triang::insert") )
{ cerr << "The point inserted was " << p; 
  exit(1);
}
return v;


   

}



bool delaunay_triang::simplify_face(node v, list<edge>& E)
{ 
  // flip edges adjacent to v until no flips possible anymore or degree(v)
  // == min_deg ( 3 if v is a non-hull node,  2 if v is a hull node)  
  // appends all flipped edges to E and  returns true iff v is hull node

  int min_deg = 3;

  edge e;
  forall_adj_edges(e,v) 
  { E.append(face_cycle_succ(e));
    if (is_hull_edge(e)) min_deg = 2;
   }

  // e1,e2,e3 :  three consecutive edges incident to v
  edge e1 = first_adj_edge(v);
  edge e2 = cyclic_adj_succ(e1);
  edge e3 = cyclic_adj_succ(e2);

  // count : number of traversed edges since last flip
  int count = 0; 

  while ( count++ < outdeg(v) && outdeg(v) > min_deg )
  { if ( ::orientation(pos_target(e1),pos_target(e3),pos_source(e2)) >= 0 &&
         ::orientation(pos_target(e1),pos_target(e3),pos_target(e2)) < 0)
      { move_edge(e2,reversal(e1),target(e3),before);
        move_edge(reversal(e2),reversal(e3),target(e1));
        E.append(e2);
        count = 0;
       }
    else 
       e1 = e2;

    e2 = e3;
    e3 = cyclic_adj_succ(e2);
  }

  return min_deg == 2;
}



void delaunay_triang::del(node v)
{  
  if (v == nil) 
     error_handler(1,"delaunay_triang::del: nil argument.");

  if (number_of_nodes() == 0) 
     error_handler(1,"delaunay_triang::del: graph is empty.");

  if (check) save_state(inf(v));

  if ( dim() < 2 )
  { 
    if ( outdeg(v) == 2)
    { node s = target(first_adj_edge(v));
      node t = target(last_adj_edge(v));
      edge x = new_edge(s,t);  
      edge y = new_edge(t,s);
      mark_edge(x,HULL_EDGE);
      mark_edge(y,HULL_EDGE);
      set_reversal(x,y);
     }

    del_node(v);
    cur_edge = hull_edge = first_edge();    
  
    if ( check && !check_state("delaunay_triang::del(node v)") )
    { cerr << "deleted the node with position " << pos(v);
      exit(1);
     }

    return;
  }

  list<edge> E;

  if ( simplify_face(v,E) )
  { // hull node: adjust marks of new hull edges and their reversals
    edge e,x;
    forall_adj_edges(e,v)
    { x = face_cycle_succ(e);
      mark_edge(x,HULL_EDGE);
      if ( !is_hull_edge(reversal(x)) ) mark_edge(reversal(x),DIAGRAM_EDGE);
     }
     hull_edge = x;
   }

  del_node(v);

  cur_edge = E.head();
  make_delaunay(E);
  
  if ( check && !check_state("delaunay_triang::del(node v)") )
  { cerr << "deleted the node with position " << pos(v);
    exit(1);
   }
}



void delaunay_triang::compute_voronoi(GRAPH<CIRCLE,POINT>& VD)
{ 
  VD.clear();

  if (number_of_nodes() < 2) return;

  // create Voronoi nodes

  edge_array<node> vnode(*this,nil);

  // for outer face

  edge  e = hull_edge;
  POINT a = pos_source(e);
  edge  x = e;
  do { POINT b = pos_target(x);
       vnode[x] =  VD.new_node(CIRCLE(a,center(a,b),b));
       a = b;
       x = face_cycle_succ(x);
     } while ( x != e);

  // for all other faces

  forall_edges(e,*this)
  { 
    if (vnode[e] || !is_diagram_edge(e)) continue;

    edge  x = face_cycle_succ(e);
    POINT a = pos_source(e);
    POINT b = pos_target(e);
    POINT c = pos_target(x);
    node  v = VD.new_node(CIRCLE(a,b,c));

    x = e;
    do { vnode[x] = v;
         x = d_face_cycle_succ(x);
        } while( x != e);
   }

  // construct Voronoi edges

  // for outer face

  e = hull_edge;
  x = e;
  do { edge r = reversal(x); 
       POINT p = pos_target(x);
       VD.new_edge(vnode[x],vnode[r],p);
       x = cyclic_adj_pred(r);
     } while ( x != e);

  // for all other faces

  forall_edges(e,*this)
  { node v = vnode[e]; 

    if (!is_diagram_edge(e) || VD.outdeg(v) > 0) continue;

    POINT p = pos_target(e);
    edge x = e;
    do { edge r = reversal(x); 
         VD.new_edge(v,vnode[r],p);
         x = d_face_cycle_succ(x);
       } while (x != e);
   }
   // add a check for correctness
}


/*
node delaunay_triang::nearest_neighbor(POINT p) const
{ 
  if (number_of_nodes() == 0) return nil;
  if (number_of_nodes() == 1) return first_node();

  edge e = locate(p);

  if ( is_hull_edge(e) ) 
  { while ( !IN_HALFSPACE(e,p) )          e = face_cycle_pred(e);
    while ( !IN_HALFSPACE(reverse(e),p) ) e = face_cycle_succ(e);
  }

  unmark_all_nodes();   // takes constant time !

  node min_v = source(e);
  RAT_TYPE min_d = p.sqr_dist(pos(min_v));

  list<node> L;
  L.append(source(e)); 
  L.append(target(e)); 
  mark_node(source(e));
  mark_node(target(e));

  while ( !L.empty() )
  { node v = L.pop();

    if ( p.sqr_dist(pos(v)) < min_d )
    { min_v = v;
      min_d = p.sqr_dist(pos(v));
    }

    forall_adj_edges(e,v)
    { node w = target(e);
      if ( IN_HALFSPACE(e,p) && IN_HALFSPACE(reversal(e),p) && 
           supporting_line(e).sqr_dist(p) < min_d  &&
           !is_marked(target(e)) )
      { L.append(w); 
        mark_node(w);
       }
    }
  }
  return min_v;

}
*/



node delaunay_triang::nearest_neighbor(POINT p) const
{
  if (number_of_nodes() == 0) return nil;
  if (number_of_nodes() == 1) return first_node();

  edge e = locate(p);

  if (p == pos_source(e)) return source(e);
  if (p == pos_target(e)) return target(e);

  // insert p and search neighbors of v

  node v = ((delaunay_triang*)this)->insert(p);  // cast away const

  e = first_adj_edge(v);

  node     min_v = target(e);
  RAT_TYPE min_d = p.sqr_dist(pos(min_v));

  while ((e = adj_succ(e)) != nil)
  { RAT_TYPE d = p.sqr_dist(pos_target(e));
    if (d < min_d)
    { min_d = d;
      min_v = target(e);
     }
   }

  ((delaunay_triang*)this)->del(v);

  return min_v;
}



void delaunay_triang::dfs(node s, const CIRCLE& C, list<node>& L) const
{ L.append(s);
  mark_node(s);
  node u;
  forall_adj_nodes(u,s)
      if (!is_marked(u) && ! C.outside(pos(u))) dfs(u,C,L);
}
 
 
 
list<node> delaunay_triang::range_search(const CIRCLE& C) const
{ 
  list<node> L;

  int orient = C.orientation();
  if (orient == 0)
      error_handler(1,"delaunay_triang::range_search: circle must be proper");

  if (number_of_nodes() == 0) return L;
  if ( number_of_nodes() == 1 && C.side_of(pos(first_node())) * orient >= 0 ) 
  { L.append(first_node());
    return L;
  }

  POINT p = C.center();
  edge e = locate(p);
  node v = nil;  
  bool new_v = true;

  if (p == pos_source(e)) { new_v = false; v = source(e); }
  if (p == pos_target(e)) { new_v = false; v = target(e); }

  // insert p and search neighbors of v

  if (new_v)
    v = ((delaunay_triang*)this) -> insert(p); // cast away const

  unmark_all_nodes(); // takes constant time !

  dfs(v,C,L);

  if (new_v)
  { L.pop();   
   ((delaunay_triang*)this) -> del(v);
   }
  return L;
}
 


list<node> delaunay_triang::range_search(const POINT& a, 
                                          const POINT& b,
                                          const POINT& c) const
{ int orient = ::orientation(a,b,c);
  CIRCLE C(a,b,c);
  list<node> L = range_search(C);
  list_item it = L.first_item();
  while (it != nil)
  { POINT p = pos(L[it]);
    list_item next_it = L.succ(it);
    if ( ::orientation(a,b,p) == - orient ||
         ::orientation(b,c,p) == - orient ||
         ::orientation(c,a,p) == - orient )      
       L.del_item(it);
    it = next_it;
  }
  return L;
}



list<node> delaunay_triang::range_search(const POINT& a, const POINT& c) const
{
  POINT b(c.xcoord(),a.ycoord());
  POINT d(a.xcoord(),c.ycoord());

  if (::orientation(a,b,c) < 0) 
  { POINT tmp = b;
    b = d;
    d = tmp;
   }

  CIRCLE C(a,b,c);

  list<node> L = range_search(C);
  list_item it = L.first_item();
  while (it != nil)
  { POINT p = pos(L[it]);
    list_item next_it = L.succ(it);
    if ( ::orientation(a,b,p) < 0 || ::orientation(b,c,p) < 0 ||
         ::orientation(c,d,p) < 0 || ::orientation(d,a,p) < 0 )
       L.del_item(it);
    it = next_it;
  }
  return L;
}
    
    

// from graph_alg.h
extern list<edge> MIN_SPANNING_TREE(const graph&, 
                                    int(*cmp)(const edge&,const edge&));
static const delaunay_triang* T_tmp;

static int cmp_edge_length(const edge& e1, const edge& e2)
{ RAT_TYPE l1 = T_tmp->pos_source(e1).sqr_dist(T_tmp->pos_target(e1)); 
  RAT_TYPE l2 = T_tmp->pos_source(e2).sqr_dist(T_tmp->pos_target(e2)); 
  return compare(l1,l2);
 }
  
list<edge> delaunay_triang::minimum_spanning_tree() const 
{ T_tmp = this; 
  return MIN_SPANNING_TREE(*this,cmp_edge_length);
}





void delaunay_triang::draw_nodes(void (*draw_node)(const POINT&)) 
{ node v;
  forall_nodes(v,*this) draw_node(pos(v)); 
 }


void delaunay_triang::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&))
{ 
    node v = source(e);
    node w = target(e);

    if (index(v) > index(w))
    { node tmp = v;
      v = w;
      w = tmp;
     }

    POINT p = pos(v);
    POINT q = pos(w);
    if (is_hull_edge(e) || is_hull_edge(reversal(e)))
       draw_hull_edge(p,q);
    else
       if (is_diagram_edge(e))
          draw_diagram_edge(p,q);
       else
          draw_triang_edge(p,q);
}



void delaunay_triang::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&))
{ edge e;
  forall_edges(e,*this) 
    if (index(source(e)) < index(target(e)))
       draw_edge(e,draw_diagram_edge,draw_triang_edge,draw_hull_edge);
}

void delaunay_triang::draw_edges(const list<edge>& L,
                      void (*draw_edge)(const POINT&,const POINT&) )
{ edge e;
  forall(e,L) 
    draw_edge(pos_source(e),pos_target(e));
}


void delaunay_triang::draw_hull(void (*draw_poly)(const list<POINT>&))
{ list<POINT> CH;
  edge e = get_hull_edge();
  if (e != nil)
  { edge x = e;
    do { POINT p = pos(source(x));
         CH.append(p);
         x = face_cycle_succ(x);
       } while (x != e);
    draw_poly(CH);
   }
 }


void  delaunay_triang::draw_voro_edge(const CIRCLE& c1, const CIRCLE& c2,
                                 void (*draw_edge)(const POINT&,const POINT&),
                                 void (*draw_ray) (const POINT&,const POINT&))
{
  int orient1 = c1.orientation();
  int orient2 = c2.orientation();

  if (orient1 == 0 && orient2 == 0)
    { POINT a = c1.point1();
      POINT b = c1.point3();
      POINT c = center(a,b);
      draw_ray(c,b.rotate90(c));
      draw_ray(c,a.rotate90(c));
     }
  else
    if (orient1 == 0)
       { POINT cw = c2.center();
         POINT p1 = c1.point1();
         POINT p3 = c1.point3();
         VECTOR vec = p3 - p1;
         POINT cv  = cw + vec.rotate90();
         draw_ray(cw,cv);
        }
    else
       if (orient2 == 0)
         { POINT cv = c1.center();
           POINT p1 = c2.point1();
           POINT p3 = c2.point3();
           VECTOR vec = p3 - p1;
           POINT cw  = cv + vec.rotate90();
           draw_ray(cv,cw);
          }
       else
         draw_edge(c1.center(),c2.center());
}


void  delaunay_triang::draw_voro_edges(
                       void (*draw_edge)(const POINT&,const POINT&),
                       void (*draw_ray) (const POINT&,const POINT&))
{ 
  if (number_of_edges() == 0) return;

  CIRCLE C_NULL;

  edge_array<CIRCLE> vnode(*this,C_NULL);

  // for outer face

  edge  e = hull_edge;

  POINT a = pos_source(e);
  edge  x = e;
  do { POINT b = pos_target(x);
       vnode[x] =  CIRCLE(a,center(a,b),b);
       a = b;
       x = face_cycle_succ(x);
     } while ( x != e);


  forall_edges(e,*this)
  { 
    if (!identical(vnode[e],C_NULL) || !is_diagram_edge(e)) continue;

    edge  x = face_cycle_succ(e);
    POINT a = pos_source(e);
    POINT b = pos_target(e);
    POINT c = pos_target(x);

    CIRCLE C(a,b,c);

    x = e;
    do { vnode[x] = C;
         x = d_face_cycle_succ(x);
        } while( x != e);
   }

  forall_edges(e,*this) 
  {  if (identical(vnode[e],C_NULL)) continue;
     edge r = reversal(e);
     draw_voro_edge(vnode[e],vnode[r],draw_edge,draw_ray);
     vnode[r] = C_NULL;
   }
}



void  delaunay_triang::draw_voro(const GRAPH<CIRCLE,POINT>& VD,
                                 void (*draw_site)   (const POINT&), 
                                 void (*draw_edge)(const POINT&,const POINT&),
                                 void (*draw_ray) (const POINT&,const POINT&))
{ 
  edge_array<bool> visited(VD,false); 

  node v;
  forall_nodes(v,*this) draw_site(pos(v));

  edge e;
  forall_edges(e,VD)
  { if (visited[e]) continue;
    visited[e] = visited[VD.reversal(e)] = true;
    node v = source(e);
    node w = target(e);
    draw_voro_edge(VD[v],VD[w],draw_edge,draw_ray);
  }
}
 

