Ticket #363: lemon-sikgraf_ea56aa8243b1.patch

lemon/Makefile.am

@@ -110 +110 @@
         lemon/network_simplex.h \
         lemon/pairing_heap.h \
         lemon/path.h \
+        lemon/planar_graph.h \
         lemon/planarity.h \
+        lemon/plane_graph.h \
         lemon/preflow.h \
         lemon/quad_heap.h \
         lemon/radix_heap.h \

lemon/bits/bezier.h

@@ -19 +19 @@
 #ifndef LEMON_BEZIER_H
 #define LEMON_BEZIER_H
 
+#include <cmath>
+
 //\ingroup misc
 //\file
 //\brief Classes to compute with Bezier curves.
@@ -66 +68 @@
   Point norm() const { return rot90(p2-p1); }
   Point grad(double) const { return grad(); }
   Point norm(double t) const { return rot90(grad(t)); }
+  Box<double> boundingBox() const {
+    double minx, miny;
+    double maxx, maxy;
+    if (p1.x < p2.x) {
+      minx = p1.x; maxx = p2.x;
+    } else {
+      minx = p2.x; maxx = p1.x;
+    }
+    if (p1.y < p2.y) {
+      miny = p1.y; maxy = p2.y;
+    } else {
+      miny = p2.y; maxy = p1.y;
+    }
+    return Box<double>(minx, miny, maxx, maxy);
+  }
 };
 
 class Bezier2 : public BezierBase
@@ -100 +117 @@
   Bezier1 norm() const { return Bezier1(2.0*rot90(p2-p1),2.0*rot90(p3-p2)); }
   Point grad(double t) const { return grad()(t); }
   Point norm(double t) const { return rot90(grad(t)); }
+  Box<double> boundingBox() const {
+    double minx, miny;
+    double maxx, maxy;
+    if (p1.x < p3.x) {
+      minx = p1.x; maxx = p3.x;
+    } else {
+      minx = p3.x; maxx = p1.x;
+    }
+    if (p1.y < p3.y) {
+      miny = p1.y; maxy = p3.y;
+    } else {
+      miny = p3.y; maxy = p1.y;
+    }
+    if (p1.x + p3.x - 2 * p2.x != 0.0) {
+      double t = (p1.x - p2.x) / (p1.x + p3.x - 2 * p2.x);
+      if (t >= 0.0 && t <= 1.0) {
+        double x = ((1-t)*(1-t))*p1.x+(2*(1-t)*t)*p2.x+(t*t)*p3.x;
+        if (x < minx) minx = x;
+        if (x > maxx) maxx = x;
+      }
+    }
+    if (p1.y + p3.y - 2 * p2.y != 0.0) {
+      double t = (p1.y - p2.y) / (p1.y + p3.y - 2 * p2.y);
+      if (t >= 0.0 && t <= 1.0) {
+        double y = ((1-t)*(1-t))*p1.y+(2*(1-t)*t)*p2.y+(t*t)*p3.y;
+        if (y < miny) miny = y;
+        if (y > maxy) maxy = y;
+      }
+    }
+    return Box<double>(minx, miny, maxx, maxy);
+  }
 };
 
 class Bezier3 : public BezierBase
@@ -150 +198 @@
                                   3.0*rot90(p4-p3)); }
   Point grad(double t) const { return grad()(t); }
   Point norm(double t) const { return rot90(grad(t)); }
+  Box<double> boundingBox() const {
+    double minx, miny;
+    double maxx, maxy;
+    if (p1.x < p4.x) {
+      minx = p1.x; maxx = p4.x;
+    } else {
+      minx = p4.x; maxx = p1.x;
+    }
+    if (p1.y < p4.y) {
+      miny = p1.y; maxy = p4.y;
+    } else {
+      miny = p4.y; maxy = p1.y;
+    }
+    if (p1.x - 3 * p2.x + 3 * p3.x - p4.x != 0.0) {
+      double disc = p2.x * p2.x + p3.x * p3.x - p1.x * p3.x 
+        + p1.x * p4.x - p2.x * p3.x - p2.x * p4.x;
+      if (disc >= 0.0) {
+        double t = (p1.x + p3.x - 2 * p2.x + std::sqrt(disc)) / 
+          (p1.x - 3 * p2.x + 3 * p3.x - p4.x);
+        if (t >= 0.0 && t <= 1.0) {
+          double x = ((1-t)*(1-t)*(1-t))*p1.x+(3*t*(1-t)*(1-t))*p2.x+
+            (3*t*t*(1-t))*p3.x+(t*t*t)*p4.x;
+          if (x < minx) minx = x;
+          if (x > maxx) maxx = x;
+        }
+        t = (p1.x + p3.x - 2 * p2.x - std::sqrt(disc)) / 
+          (p1.x - 3 * p2.x + 3 * p3.x - p4.x);
+        if (t >= 0.0 && t <= 1.0) {
+          double x = ((1-t)*(1-t)*(1-t))*p1.x+(3*t*(1-t)*(1-t))*p2.x+
+            (3*t*t*(1-t))*p3.x+(t*t*t)*p4.x;
+          if (x < minx) minx = x;
+          if (x > maxx) maxx = x;
+        }
+      }
+    } else if (p1.x + p3.x - 2 * p2.x != 0.0) {
+      double t = (p1.x - p2.x) / (2 * p1.x + 2 * p3.x - 4 * p2.x);
+      if (t >= 0.0 && t <= 1.0) {
+        double x = ((1-t)*(1-t)*(1-t))*p1.x+(3*t*(1-t)*(1-t))*p2.x+
+          (3*t*t*(1-t))*p3.x+(t*t*t)*p4.x;
+        if (x < minx) minx = x;
+        if (x > maxx) maxx = x;
+      }      
+    }
+    if (p1.y - 3 * p2.y + 3 * p3.y - p4.y != 0.0) {
+      double disc = p2.y * p2.y + p3.y * p3.y - p1.y * p3.y 
+        + p1.y * p4.y - p2.y * p3.y - p2.y * p4.y;
+      if (disc >= 0.0) {
+        double t = (p1.y + p3.y - 2 * p2.y + std::sqrt(disc)) / 
+          (p1.y - 3 * p2.y + 3 * p3.y - p4.y);
+        if (t >= 0.0 && t <= 1.0) {
+          double y = ((1-t)*(1-t)*(1-t))*p1.y+(3*t*(1-t)*(1-t))*p2.y+
+            (3*t*t*(1-t))*p3.y+(t*t*t)*p4.y;
+          if (y < miny) miny = y;
+          if (y > maxy) maxy = y;
+        }
+        t = (p1.y + p3.y - 2 * p2.y - std::sqrt(disc)) / 
+          (p1.y - 3 * p2.y + 3 * p3.y - p4.y);
+        if (t >= 0.0 && t <= 1.0) {
+          double y = ((1-t)*(1-t)*(1-t))*p1.y+(3*t*(1-t)*(1-t))*p2.y+
+            (3*t*t*(1-t))*p3.y+(t*t*t)*p4.y;
+          if (y < miny) miny = y;
+          if (y > maxy) maxy = y;
+        }
+      }
+    } else if (p1.y + p3.y - 2 * p2.y != 0.0) {
+      double t = (p1.y - p2.y) / (2 * p1.y + 2 * p3.y - 4 * p2.y);
+      if (t >= 0.0 && t <= 1.0) {
+        double y = ((1-t)*(1-t)*(1-t))*p1.y+(3*t*(1-t)*(1-t))*p2.y+
+          (3*t*t*(1-t))*p3.y+(t*t*t)*p4.y;
+        if (y < miny) miny = y;
+        if (y > maxy) maxy = y;
+      }      
+    }
+    return Box<double>(minx, miny, maxx, maxy);
+  }
 
   template<class R,class F,class S,class D>
   R recSplit(F &_f,const S &_s,D _d) const

lemon/graph_to_eps.h

@@ -103 +103 @@
   double _arrowLength, _arrowWidth;
 
   bool _showNodes, _showArcs;
+  bool _customEdgeShapes;
 
   bool _enableParallel;
   double _parArcDist;
@@ -154 +155 @@
     _nodeScale(.01), _xBorder(10), _yBorder(10), _scale(1.0),
     _nodeBorderQuotient(.1),
     _drawArrows(false), _arrowLength(1), _arrowWidth(0.3),
-    _showNodes(true), _showArcs(true),
+    _showNodes(true), _showArcs(true), _customEdgeShapes(false),
     _enableParallel(false), _parArcDist(1),
     _showNodeText(false), _nodeTexts(false), _nodeTextSize(1),
     _showNodePsText(false), _nodePsTexts(false), _nodePsTextsPreamble(0),
@@ -201 +202 @@
 
   using T::_showNodes;
   using T::_showArcs;
+  using T::_customEdgeShapes;
 
   using T::_enableParallel;
   using T::_parArcDist;
@@ -579 +581 @@
   GraphToEps<T> &hideArcs(bool b=true) {_showArcs=!b;return *this;}
   ///Hides the nodes
   GraphToEps<T> &hideNodes(bool b=true) {_showNodes=!b;return *this;}
+  ///Uses default edge shapes
+  GraphToEps<T> &customEdgeShapes(bool b=false) {
+    _customEdgeShapes=b;return *this;
+  }
+
 
   ///Sets the size of the node texts
   GraphToEps<T> &nodeTextSize(double d) {_nodeTextSize=d;return *this;}
@@ -655 +662 @@
 public:
   ~GraphToEps() { }
 
-  ///Draws the graph.
+  template<typename GR, typename D>
+  void edgeToEps(const GR &g, std::ostream &os, const typename GR::Arc e,
+    const Color &col, double width, D t) {
+  }
 
-  ///Like other functions using
-  ///\ref named-templ-func-param "named template parameters",
-  ///this function calls the algorithm itself, i.e. in this case
-  ///it draws the graph.
-  void run() {
+  template<typename GR>
+  void edgeToEps(const GR &g, std::ostream &os, const typename GR::Arc e,
+    const Color &col, double width, bool b) {
+    g.edgeToEps(os,e,_arcColors[e],_arcWidths[e]*_arcWidthScale);
+  }
+
+  template<typename D>
+  void inner_run(D t) {
     const double EPSILON=1e-9;
     if(dontPrint) return;
 
@@ -868 +881 @@
 
     if(_showArcs) {
       os << "%Arcs:\ngsave\n";
-      if(_enableParallel) {
+      if (_customEdgeShapes) {
+        for(ArcIt e(g);e!=INVALID;++e) {
+          if((!_undirected||g.source(e) == g.target(e)||g.source(e)<g.target(e))
+              &&_arcWidths[e]>0) {
+            edgeToEps(g,os,e,_arcColors[e],_arcWidths[e]*_arcWidthScale,t);
+          }
+        }
+      } else if (_enableParallel) {
         std::vector<Arc> el;
         for(ArcIt e(g);e!=INVALID;++e)
           if((!_undirected||g.source(e)<g.target(e))&&_arcWidths[e]>0
@@ -1056 +1076 @@
     if(_pleaseRemoveOsStream) {delete &os;}
   }
 
+  ///Draws the graph.
+
+  ///Like other functions using
+  ///\ref named-templ-func-param "named template parameters",
+  ///this function calls the algorithm itself, i.e. in this case
+  ///it draws the graph.
+  void run(int i = 0) {
+    this->inner_run(i); //for normal graphs
+  }
+
+  void run(bool b) {
+    this->inner_run(b); //for plane graphs with custom edge shape support
+  }
+
   ///\name Aliases
   ///These are just some aliases to other parameter setting functions.
 
@@ -1150 +1184 @@
 ///\sa graphToEps(GR &g, std::ostream& os)
 template<class GR>
 GraphToEps<DefaultGraphToEpsTraits<GR> >
-graphToEps(GR &g,const char *file_name)
+graphToEps(const GR &g,const char *file_name)
 {
   std::ostream* os = new std::ofstream(file_name);
   if (!(*os)) {

new file lemon/planar_graph.h

@@ +1 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2010
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_PLANAR_GRAPH_H
+#define LEMON_PLANAR_GRAPH_H
+
+///\ingroup graphs
+///\file
+///\brief PlanarGraph classes. UNDER CONSTRUCTION.
+
+#include <lemon/core.h>
+#include <lemon/error.h>
+#include <lemon/bits/graph_extender.h>
+#include <lemon/planarity.h>
+
+#include <vector>
+#include <list>
+#include <set>
+
+#ifdef REMOVE_BEFORE_RELEASE
+#include <iostream>
+#endif
+
+namespace lemon {
+
+  class PlanarGraphBase {
+
+  protected:
+
+    struct NodeT {
+      int first_out, last_out;
+      int prev, next;
+      int outer_face;
+      int component;
+    };
+
+    struct ArcT {
+      int target;
+      int left_face;
+      int prev_out, next_out;
+    };
+
+    struct FaceT {
+      int first_arc;
+      int prev, next;
+    };
+
+    struct ComponentT {
+      int num_nodes;
+      int prev, next;
+    };
+
+    std::vector<NodeT> nodes;
+    int first_node;
+    int first_free_node;
+
+    std::vector<ArcT> arcs;
+    int first_free_arc;
+
+    std::vector<FaceT> faces;
+    int first_face;
+    int first_free_face;
+
+    std::vector<ComponentT> components;
+    int first_component;
+    int first_free_component;
+
+  public:
+
+    typedef PlanarGraphBase Graph;
+
+    class Node {
+      friend class PlanarGraphBase;
+    protected:
+
+      int id;
+      explicit Node(int pid) {
+        id = pid;
+      }
+
+    public:
+      Node() {}
+      Node (Invalid) {
+        id = -1;
+      }
+      bool operator==(const Node& node) const {
+        return id == node.id;
+      }
+      bool operator!=(const Node& node) const {
+        return id != node.id;
+      }
+      bool operator<(const Node& node) const {
+        return id < node.id;
+      }
+    };
+
+    class Edge {
+      friend class PlanarGraphBase;
+    protected:
+
+      int id;
+      explicit Edge(int pid) {
+        id = pid;
+      }
+
+    public:
+      Edge() {}
+      Edge (Invalid) {
+        id = -1;
+      }
+      bool operator==(const Edge& edge) const {
+        return id == edge.id;
+      }
+      bool operator!=(const Edge& edge) const {
+        return id != edge.id;
+      }
+      bool operator<(const Edge& edge) const {
+        return id < edge.id;
+      }
+    };
+
+    class Arc {
+      friend class PlanarGraphBase;
+    protected:
+
+      int id;
+      explicit Arc(int pid) {
+        id = pid;
+      }
+
+    public:
+      operator Edge() const {
+        return id != -1 ? edgeFromId(id / 2) : INVALID;
+      }
+
+      Arc() {}
+      Arc (Invalid) {
+        id = -1;
+      }
+      bool operator==(const Arc& arc) const {
+        return id == arc.id;
+      }
+      bool operator!=(const Arc& arc) const {
+        return id != arc.id;
+      }
+      bool operator<(const Arc& arc) const {
+        return id < arc.id;
+      }
+    };
+
+    class Face {
+      friend class PlanarGraphBase;
+    protected:
+
+      int id;
+      explicit Face(int pid) {
+        id = pid;
+      }
+
+    public:
+      Face() {}
+      Face (Invalid) {
+        id = -1;
+      }
+      bool operator==(const Face& face) const {
+        return id == face.id;
+      }
+      bool operator!=(const Face& face) const {
+        return id != face.id;
+      }
+      bool operator<(const Face& face) const {
+        return id < face.id;
+      }
+    };
+
+    PlanarGraphBase()
+        : nodes(), first_node(-1), first_free_node(-1),
+        arcs(), first_free_arc(-1),
+        faces(), first_face(-1), first_free_face(-1),
+        components(), first_component(-1), first_free_component(-1) {
+    }
+
+
+    int maxNodeId() const {
+      return nodes.size()-1;
+    }
+    int maxEdgeId() const {
+      return arcs.size() / 2 - 1;
+    }
+    int maxArcId() const {
+      return arcs.size()-1;
+    }
+    int maxFaceId() const {
+      return faces.size()-1;
+    }
+
+    int numComponents() const {
+      int result = 0;
+      int c = first_component;
+      while (c != -1) {
+        ++result;
+        c = components[c].next;
+      }
+      return result;
+    }
+
+    Node source(Arc e) const {
+      return Node(arcs[e.id ^ 1].target);
+    }
+    Node target(Arc e) const {
+      return Node(arcs[e.id].target);
+    }
+    Face leftFace(Arc e) const {
+      return Face(arcs[e.id].left_face);
+    }
+    Face rightFace(Arc e) const {
+      return Face(arcs[e.id ^ 1].left_face);
+    }
+    Face outerFace(Node n) const {
+      return Face(nodes[n.id].outer_face);
+    }
+
+    Node u(Edge e) const {
+      return Node(arcs[2 * e.id].target);
+    }
+    Node v(Edge e) const {
+      return Node(arcs[2 * e.id + 1].target);
+    }
+    Face w1(Edge e) const {
+      return Face(arcs[2 * e.id].left_face);
+    }
+    Face w2(Edge e) const {
+      return Face(arcs[2 * e.id + 1].left_face);
+    }
+
+    static bool direction(Arc e) {
+      return (e.id & 1) == 1;
+    }
+
+    static Arc direct(Edge e, bool d) {
+      return Arc(e.id * 2 + (d ? 1 : 0));
+    }
+
+    //Primitives to use in iterators
+    void first(Node& node) const {
+      node.id = first_node;
+    }
+
+    void next(Node& node) const {
+      node.id = nodes[node.id].next;
+    }
+
+    void first(Arc& e) const {
+      int n = first_node;
+      while (n != -1 && nodes[n].first_out == -1) {
+        n = nodes[n].next;
+      }
+      e.id = (n == -1) ? -1 : nodes[n].first_out;
+    }
+
+    void next(Arc& e) const {
+      if (arcs[e.id].next_out != -1) {
+        e.id = arcs[e.id].next_out;
+      } else {
+        int n = nodes[arcs[e.id ^ 1].target].next;
+        while (n != -1 && nodes[n].first_out == -1) {
+          n = nodes[n].next;
+        }
+        e.id = (n == -1) ? -1 : nodes[n].first_out;
+      }
+    }
+
+    void first(Edge& e) const {
+      int n = first_node;
+      while (n != -1) {
+        e.id = nodes[n].first_out;
+        while ((e.id & 1) != 1) {
+          e.id = arcs[e.id].next_out;
+        }
+        if (e.id != -1) {
+          e.id /= 2;
+          return;
+        }
+        n = nodes[n].next;
+      }
+      e.id = -1;
+    }
+
+    void next(Edge& e) const {
+      int n = arcs[e.id * 2].target;
+      e.id = arcs[(e.id * 2) | 1].next_out;
+      while ((e.id & 1) != 1) {
+        e.id = arcs[e.id].next_out;
+      }
+      if (e.id != -1) {
+        e.id /= 2;
+        return;
+      }
+      n = nodes[n].next;
+      while (n != -1) {
+        e.id = nodes[n].first_out;
+        while ((e.id & 1) != 1) {
+          e.id = arcs[e.id].next_out;
+        }
+        if (e.id != -1) {
+          e.id /= 2;
+          return;
+        }
+        n = nodes[n].next;
+      }
+      e.id = -1;
+    }
+
+    void firstOut(Arc &e, const Node& v) const {
+      e.id = nodes[v.id].first_out;
+    }
+    void nextOut(Arc &e) const {
+      e.id = arcs[e.id].next_out;
+    }
+
+    void lastOut(Arc &e, const Node& v) const {
+      e.id = nodes[v.id].last_out;
+    }
+    void prevOut(Arc &e) const {
+      e.id = arcs[e.id].prev_out;
+    }
+
+    void firstIn(Arc &e, const Node& v) const {
+      e.id = ((nodes[v.id].first_out) ^ 1);
+      if (e.id == -2) e.id = -1;
+    }
+    void nextIn(Arc &e) const {
+      e.id = ((arcs[e.id ^ 1].next_out) ^ 1);
+      if (e.id == -2) e.id = -1;
+    }
+    void lastIn(Arc &e, const Node& v) const {
+      e.id = ((nodes[v.id].last_out) ^ 1);
+      if (e.id == -2) e.id = -1;
+    }
+    void prevIn(Arc &e) const {
+      e.id = ((arcs[e.id ^ 1].prev_out) ^ 1);
+      if (e.id == -2) e.id = -1;
+    }
+
+    void firstCcwF(Arc &e, const Face &f) const {
+      e.id = faces[f.id].first_arc;
+      turnRight(e);
+      setToOpposite(e);
+    }
+    void nextCcwF(Arc &e) const {
+      if (e.id == faces[arcs[e.id].left_face].first_arc) e = INVALID;
+      else {
+        turnRight(e);
+        setToOpposite(e);
+      }
+    }
+    void firstCwF(Arc &e, const Face &f) const {
+      e.id = faces[f.id].first_arc;
+    }
+    void nextCwF(Arc &e) const {
+      turnLeft(e);
+      setToOpposite(e);
+      if (e.id == faces[arcs[e.id].left_face].first_arc) e = INVALID;
+    }
+
+    void firstInF(Arc &e, const Face &f) const {
+      e.id = faces[f.id].first_arc;
+    }
+    void nextInF(Arc &e) const {
+      setToOpposite(e);
+      turnRight(e);
+      if (e.id == faces[arcs[e.id].left_face].first_arc) e = INVALID;
+    }
+    void lastInF(Arc &e, const Face &f) const {
+      e.id = faces[f.id].first_arc;
+      setToOpposite(e);
+      turnRight(e);
+    }
+    void prevInF(Arc &e) const {
+      if (e.id == faces[arcs[e.id].left_face].first_arc) {
+        e = INVALID;
+        return;
+      }
+      setToOpposite(e);
+      turnRight(e);
+    }
+
+    void firstOutF(Arc &e, const Face &f) const {
+      e.id = faces[e.id].first_arc ^ 1;
+    }
+    void nextOutF(Arc &e) const {
+      turnRight(e);
+      setToOpposite(e);
+      if (e.id == faces[arcs[e.id ^ 1].left_face].first_arc) e = INVALID;
+    }
+
+    void first(Arc &arc, const Face& face) const {
+      arc.id = faces[face.id].first_arc;
+    }
+
+    void turnLeftF(Arc &e) const {
+      setToOpposite(e);
+      turnRight(e);
+    }
+
+    void turnRightF(Arc &e) const {
+      turnLeft(e);
+      setToOpposite(e);
+    }
+
+    void turnLeft(Arc &e) const {
+      if (arcs[e.id].next_out > -1) {
+        e.id = arcs[e.id].next_out;
+      } else {
+        e.id = nodes[arcs[e.id ^ 1].target].first_out;
+      }
+    }
+    void turnRight(Arc &e) const {
+      if (arcs[e.id].prev_out > -1) {
+        e.id = arcs[e.id].prev_out;
+      } else {
+        e.id = nodes[arcs[e.id ^ 1].target].last_out;
+      }
+    }
+    void setToOpposite(Arc &a) const {
+      if (a.id != -1) a.id ^= 1;
+    }
+
+    void firstInc(Edge &e, bool& d, const Node& v) const {
+      int a = nodes[v.id].first_out;
+      if (a != -1 ) {
+        e.id = a / 2;
+        d = ((a & 1) == 1);
+      } else {
+        e.id = -1;
+        d = true;
+      }
+    }
+    void nextInc(Edge &e, bool& d) const {
+      int a = (arcs[(e.id * 2) | (d ? 1 : 0)].next_out);
+      if (a != -1 ) {
+        e.id = a / 2;
+        d = ((a & 1) == 1);
+      } else {
+        e.id = -1;
+        d = true;
+      }
+    }
+
+    void first(Face& face) const {
+      face.id = first_face;
+    }
+
+    void next(Face& face) const {
+      face.id = faces[face.id].next;
+    }
+
+    Arc arcAt(Node n, Edge e) {
+      if (e.id == -1) return INVALID;
+      return Arc((e.id*2) | (arcs[e.id*2].target == n.id?1:0));
+    }
+
+    static int id(Node v) {
+      return v.id;
+    }
+    static int id(Arc e) {
+      return e.id;
+    }
+    static int id(Edge e) {
+      return e.id;
+    }
+    static int id(Face f) {
+      return f.id;
+    }
+
+    static Node nodeFromId(int id) {
+      return Node(id);
+    }
+    static Arc arcFromId(int id) {
+      return Arc(id);
+    }
+    static Edge edgeFromId(int id) {
+      return Edge(id);
+    }
+    static Face faceFromId(int id) {
+      return Face(id);
+    }
+
+    bool valid(Node n) const {
+      return n.id >= 0 && n.id < static_cast<int>(nodes.size()) &&
+             nodes[n.id].prev != -2;
+    }
+
+    bool valid(Arc a) const {
+      return a.id >= 0 && a.id < static_cast<int>(arcs.size()) &&
+             arcs[a.id].prev_out != -2;
+    }
+
+    bool valid(Edge e) const {
+      return e.id >= 0 && 2 * e.id < static_cast<int>(arcs.size()) &&
+             arcs[2 * e.id].prev_out != -2;
+    }
+
+    bool valid(Face f) const {
+      return f.id >= 0 && f.id < static_cast<int>(faces.size()) &&
+             faces[f.id].prev != -2;
+    }
+
+    Node addNode() {
+      int n = link_node();
+
+      nodes[n].component = addComponent();
+      components[nodes[n].component].num_nodes = 1;
+      nodes[n].outer_face = addFace().id;
+
+      return Node(n);
+    }
+
+    Edge addEdge(Node u, Node v, Edge e_u, Edge e_v) {
+
+      Arc p_u = arcAt(u,e_u);
+      Arc p_v = arcAt(v,e_v);
+
+      if (p_u.id > -1 && p_v.id > -1 && arcs[p_u.id].left_face != arcs[p_v.id].
+          left_face && nodes[u.id].component == nodes[v.id].component) return
+          INVALID;
+
+      int n = link_edge();
+
+      arcs[n].target = u.id;
+      arcs[n | 1].target = v.id;
+
+      arcs[n].prev_out = p_v.id;
+      if (p_v.id > -1) {
+        arcs[n].next_out = arcs[p_v.id].next_out;
+        arcs[p_v.id].next_out = n;
+      } else {
+        arcs[n].next_out = nodes[v.id].first_out;
+        nodes[v.id].first_out = n;
+      }
+      if (arcs[n].next_out > -1) {
+        arcs[arcs[n].next_out].prev_out = n;
+      } else {
+        nodes[v.id].last_out = n;
+      }
+
+      arcs[n | 1].prev_out = p_u.id;
+      if (p_u.id > -1) {
+        arcs[n | 1].next_out = arcs[p_u.id].next_out;
+        arcs[p_u.id].next_out = n | 1;
+      } else {
+        arcs[n | 1].next_out = nodes[u.id].first_out;
+        nodes[u.id].first_out = n | 1;
+      }
+      if (arcs[n | 1].next_out > -1) {
+        arcs[arcs[n | 1].next_out].prev_out = n | 1;
+      } else {
+        nodes[u.id].last_out = n | 1;
+      }
+
+      int ca = nodes[u.id].component;
+      int cb = nodes[v.id].component;
+
+      //Add the extra face, if needed
+      if (p_u.id > -1 & p_v.id > -1) {
+        int oldf = arcs[p_u.id].left_face;
+        int oldfb = arcs[p_v.id].left_face;
+        arcs[n].left_face = arcs[n | 1].left_face = oldf;
+        Face f = addFace();
+        faces[f.id].first_arc = n | 1;
+        faces[oldf].first_arc = n;
+        Arc arc(n | 1);
+        wall_paint(arc,f.id,arc);
+        if (nodes[v.id].component != nodes[u.id].component) {
+          erase(Face(oldf));
+          erase(Face(oldfb));
+          if (components[ca].num_nodes > components[cb].num_nodes) {
+            component_relabel(v,ca);
+            eraseComponent(cb);
+          } else {
+            component_relabel(u,cb);
+            eraseComponent(ca);
+          }
+        }
+      } else if (p_u.id > -1) {
+        arcs[n].left_face = arcs[n | 1].left_face = arcs[p_u.id].left_face;
+        faces[arcs[n].left_face].first_arc = n | 1;
+        nodes[v.id].component = nodes[u.id].component;
+        ++components[ca].num_nodes;
+        eraseComponent(cb);
+        erase(Face(nodes[v.id].outer_face));
+      } else if (p_v.id > -1) {
+        arcs[n].left_face = arcs[n | 1].left_face = arcs[p_v.id].left_face;
+        faces[arcs[n].left_face].first_arc = n | 1;
+        ++components[cb].num_nodes;
+        eraseComponent(ca);
+        nodes[u.id].component = nodes[v.id].component;
+        erase(Face(nodes[u.id].outer_face));
+      } else if (u.id != v.id) {    //both prevs are INVALID
+        Face f(nodes[u.id].outer_face);
+        arcs[n].left_face = arcs[n | 1].left_face = f.id;
+        faces[f.id].first_arc = n | 1;
+        ++components[ca].num_nodes;
+        eraseComponent(cb);
+        nodes[v.id].component = nodes[u.id].component;
+        erase(Face(nodes[v.id].outer_face));
+      } else {  //insert a loop to a new node
+        Face f(nodes[u.id].outer_face);
+        Face f2 = addFace();
+        arcs[n].left_face = f2.id;
+        arcs[n | 1].left_face = f.id;
+        faces[f.id].first_arc = n | 1;
+        faces[f2.id].first_arc = n;
+      }
+
+      nodes[u.id].outer_face = nodes[v.id].outer_face = -1;
+      return Edge(n / 2);
+    }
+
+    void erase(const Node& node) {
+      int n = node.id;
+
+      eraseComponent(nodes[n].component);
+      erase(Face(nodes[n].outer_face));
+
+      unlink_node(n);
+    }
+
+    void erase(const Edge& edge) {
+      int n = edge.id * 2;
+
+      //"retreat" the incident faces' first arcs
+      int fl = arcs[n].left_face;
+      if ((faces[fl].first_arc | 1) == (n | 1)) {
+        Arc e(faces[fl].first_arc);
+        turnRightF(e);
+        if ((e.id | 1) == (n | 1)) turnRightF(e);
+        faces[fl].first_arc = e.id;
+      }
+
+      int fr = arcs[n | 1].left_face;
+
+      bool comp_split = false;
+      if (fr != fl) {
+        Arc arc(faces[fr].first_arc);
+        wall_paint(arc,fl,arc);
+        if ((faces[fl].first_arc | 1) == (n | 1))
+          faces[fl].first_arc = faces[fr].first_arc;
+        erase(Face(fr));
+      } else {
+        comp_split = true;
+        if ((arcs[n].next_out > -1 || arcs[n].prev_out > -1) &&
+                 (arcs[n | 1].next_out > -1 || arcs[n | 1].prev_out > -1)) {
+          Arc arc(n);
+          Arc ed = arc;
+          ed.id ^= 1;
+          turnRightF(arc);
+          Face f = addFace();
+          wall_paint(arc,f.id,ed);
+          faces[f.id].first_arc = arc.id;
+          turnRightF(ed);
+          faces[fr].first_arc = ed.id;
+        }
+      }
+
+      if (arcs[n].next_out != -1) {
+        arcs[arcs[n].next_out].prev_out = arcs[n].prev_out;
+      } else {
+        nodes[arcs[n].target].last_out = arcs[n].prev_out;
+      }
+
+      if (arcs[n].prev_out != -1) {
+        arcs[arcs[n].prev_out].next_out = arcs[n].next_out;
+      } else {
+        nodes[arcs[n | 1].target].first_out = arcs[n].next_out;
+      }
+
+      if (arcs[n | 1].next_out != -1) {
+        arcs[arcs[n | 1].next_out].prev_out = arcs[n | 1].prev_out;
+      } else {
+        nodes[arcs[n | 1].target].last_out = arcs[n | 1].prev_out;
+      }
+
+      if (arcs[n | 1].prev_out != -1) {
+        arcs[arcs[n | 1].prev_out].next_out = arcs[n | 1].next_out;
+      } else {
+        nodes[arcs[n].target].first_out = arcs[n | 1].next_out;
+      }
+
+      unlink_edge(edge);
+
+      if (comp_split) component_relabel(Node(arcs[n | 1].target),
+        addComponent());
+      if (nodes[arcs[n].target].first_out == -1 && nodes[arcs[n | 1].target].
+              first_out == -1) {
+        int t1 = arcs[n].target;
+        int t2 = arcs[n | 1].target;
+        if (t1 != t2) {
+          nodes[t1].outer_face = addFace().id;
+          nodes[t2].outer_face = fr;
+        } else {
+          faces[fl].first_arc = -1;
+          nodes[t1].outer_face = fl;
+        }
+      } else if (nodes[arcs[n].target].first_out == -1)
+        nodes[arcs[n].target].outer_face = addFace().id;
+      else if (nodes[arcs[n | 1].target].first_out == -1)
+        nodes[arcs[n | 1].target].outer_face = addFace().id;
+
+    }
+
+    void clear() {
+      arcs.clear();
+      nodes.clear();
+      faces.clear();
+      components.clear();
+      first_node = first_free_node = first_free_arc = first_face =
+        first_free_face = first_component = first_free_component = -1;
+    }
+
+    Node split(Edge e) {
+      Node v;
+      Edge e2;
+      inner_split1(v,e2);
+      return inner_split2(e,v,e2);
+    }
+
+    Node split(Node n1, Edge e1, Edge e2, bool b) {
+      Node n2;
+      inner_split1(n1,n2);
+      Edge a3;
+      inner_split2(n1,n2,e1,e2,b,a3);
+      if (!b && a3 != INVALID) erase(a3);
+      return n2;
+    }
+
+    void contract(Edge e) {
+      Node n1 = u(e);
+      Node n2 = v(e);
+      Node nd;
+      bool simple;
+      inner_contract1(n1,n2,nd,simple);
+      inner_contract2(e,n1,n2,nd,simple);
+    }
+
+  protected:
+
+    //Adds a node to the linked list and sets default properties. Used by
+    //addNode and other functions that take care of the planar properties by
+    //themselves.
+    int link_node() {
+      int n;
+      if (first_free_node==-1) {
+        n = nodes.size();
+        nodes.push_back(NodeT());
+      } else {
+        n = first_free_node;
+        first_free_node = nodes[n].next;
+      }
+
+      nodes[n].next = first_node;
+      if (first_node != -1) nodes[first_node].prev = n;
+      first_node = n;
+      nodes[n].prev = -1;
+      nodes[n].first_out = nodes[n].last_out = -1;
+      return n;
+    }
+
+    //Removes a node from the linked list. Used by
+    //erase(Node) and other functions that take care of the planar properties by
+    //themselves.
+    void unlink_node(int n) {
+      if (nodes[n].next != -1) {
+        nodes[nodes[n].next].prev = nodes[n].prev;
+      }
+
+      if (nodes[n].prev != -1) {
+        nodes[nodes[n].prev].next = nodes[n].next;
+      } else {
+        first_node = nodes[n].next;
+      }
+
+      nodes[n].next = first_free_node;
+      first_free_node = n;
+      nodes[n].prev = -2;
+    }
+
+    //Adds an edge to the linked list. Used by
+    //addEdge and other functions that take care of the planar properties by
+    //themselves.
+    int link_edge() {
+      int n;
+      if (first_free_arc == -1) {
+        n = arcs.size();
+        arcs.push_back(ArcT());
+        arcs.push_back(ArcT());
+      } else {
+        n = first_free_arc;
+        first_free_arc = arcs[n].next_out;
+      }
+      return n;
+    }
+
+    //Removes an edge from the linked list. Used by
+    //erase(Edge) and other functions that take care of the planar properties by
+    //themselves.
+    void unlink_edge(Edge e) {
+      int n = e.id*2;
+      arcs[n].next_out = first_free_arc;
+      first_free_arc = n;
+      arcs[n].prev_out = -2;
+      arcs[n | 1].prev_out = -2;
+    }
+
+    //Primitives of split(Edge)
+    void inner_split1(Node &v, Edge &e2) {
+      v = Node(link_node());
+      e2 = Arc(link_edge());
+    }
+
+    Node inner_split2(Edge e, Node v, Edge e2) {
+      Arc a(e.id*2);
+      int b = e2.id*2;
+      nodes[v.id].component = nodes[arcs[a.id].target].component;
+      ++components[nodes[v.id].component].num_nodes;
+      nodes[v.id].first_out = a.id;
+      nodes[v.id].last_out = b | 1;
+
+      arcs[b] = arcs[a.id];
+      arcs[b].target = v.id;
+      if (arcs[a.id].next_out > -1)
+        arcs[arcs[a.id].next_out].prev_out = b;
+      else
+        nodes[arcs[a.id | 1].target].last_out = b;
+      if (arcs[a.id].prev_out > -1)
+        arcs[arcs[a.id].prev_out].next_out = b;
+      else
+        nodes[arcs[a.id | 1].target].first_out = b;
+
+      arcs[b | 1] = arcs[a.id | 1];
+      arcs[b | 1].next_out = -1;
+      arcs[b | 1].prev_out = a.id;
+
+      arcs[a.id].next_out = b | 1;
+      arcs[a.id].prev_out = -1;
+      arcs[a.id | 1].target = v.id;
+
+      return v;
+    }
+
+    //Primitives of contract(Edge)
+    void inner_contract1(Node n1, Node n2, Node &to_delete, bool &simple) {
+      if (nodes[n1.id].first_out == nodes[n1.id].last_out) {
+        simple = true;
+        to_delete = n1;
+      } else if (nodes[n2.id].first_out == nodes[n2.id].last_out) {
+        simple = true;
+        to_delete = n2;
+      } else {
+        simple = false;
+        to_delete = n2;
+      }
+    }
+
+    void inner_contract2(Edge e, Node n1, Node n2, Node &to_delete, bool
+      &simple) {
+      if (simple) {
+        erase(e);
+        erase(to_delete);
+        return;
+      }
+      Arc a(e.id*2);
+      int fl = arcs[a.id].left_face;
+      if ((faces[fl].first_arc | 1) == a.id | 1) {
+        Arc b = a;
+        turnRightF(b);
+        faces[fl].first_arc = b.id;
+      }
+      int fr = arcs[a.id | 1].left_face;
+      if ((faces[fr].first_arc | 1) == a.id | 1) {
+        Arc b(a.id | 1);
+        turnRightF(b);
+        faces[fr].first_arc = b.id;
+      }
+      Arc b = a;
+      turnLeft(b);
+      while (b != a) {
+        arcs[b.id ^ 1].target = n1.id;
+        if (arcs[b.id].next_out > -1)
+          b.id = arcs[b.id].next_out;
+        else
+          b.id = nodes[n2.id].first_out;
+      }
+      arcs[nodes[n2.id].last_out].next_out = nodes[n2.id].first_out;
+      arcs[nodes[n2.id].first_out].prev_out = nodes[n2.id].last_out;
+      if (arcs[a.id | 1].prev_out > -1) {
+        arcs[arcs[a.id | 1].prev_out].next_out = arcs[a.id].next_out;
+      } else {
+        nodes[n1.id].first_out = arcs[a.id].next_out;
+      }
+      arcs[arcs[a.id].next_out].prev_out = arcs[a.id | 1].prev_out;
+      if (arcs[a.id | 1].next_out > -1) {
+        arcs[arcs[a.id | 1].next_out].prev_out = arcs[a.id].prev_out;
+      } else {
+        nodes[n1.id].last_out = arcs[a.id].prev_out;
+      }
+      arcs[arcs[a.id].prev_out].next_out = arcs[a.id | 1].next_out;
+
+      unlink_edge(e);
+      --components[nodes[n2.id].component].num_nodes;
+      unlink_node(n2.id);
+    }
+
+    //Primitives of split(Node)
+    void inner_split1(Node n1, Node &n2) {
+      n2 = Node(link_node());
+    }
+
+    void inner_split2(Node n1, Node n2, Edge e1, Edge e2, bool b, Edge &e3) {
+      e3 = INVALID;
+      if (nodes[n1.id].first_out == -1) {
+        if (b) e3 = addEdge(n1,n2,INVALID,INVALID);
+        return;
+      }
+      arcs[nodes[n1.id].first_out].prev_out = nodes[n1.id].last_out;
+      arcs[nodes[n1.id].last_out].next_out = nodes[n1.id].first_out;
+      nodes[n2.id].component = nodes[n1.id].component;
+      ++components[nodes[n2.id].component].num_nodes;
+      Arc a1(e1.id*2);
+      if (arcs[a1.id].target != n1.id) a1.id |= 1;
+      Arc a2(e2.id*2);
+      if (arcs[a2.id].target != n1.id) a2.id |= 1;
+      Arc a3(link_edge());
+      e3 = a3;
+      arcs[a3.id].target = n1.id;
+      arcs[a3.id | 1].target = n2.id;
+      arcs[a3.id].left_face = arcs[a1.id].left_face;
+      arcs[a3.id | 1].left_face = arcs[a2.id].left_face;
+
+      arcs[a3.id].prev_out = arcs[a2.id ^ 1].prev_out;
+      arcs[arcs[a2.id ^ 1].prev_out].next_out = a3.id;
+      arcs[a3.id | 1].prev_out = arcs[a1.id ^ 1].prev_out;
+      arcs[arcs[a1.id ^ 1].prev_out].next_out = a3.id | 1;
+      nodes[n1.id].first_out = a2.id ^ 1;
+      arcs[a1.id ^ 1].prev_out = -1;
+      nodes[n1.id].last_out = a3.id | 1;
+      arcs[a3.id | 1].next_out = -1;
+      nodes[n2.id].first_out = a1.id ^ 1;
+      arcs[a2.id ^ 1].prev_out = -1;
+      nodes[n2.id].last_out = a3.id;
+      arcs[a3.id].next_out = -1;
+
+      Arc a4(a1.id);
+      while (a4.id != (a3.id | 1)) {
+        arcs[a4.id].target = n2.id;
+        a4.id = arcs[a4.id ^ 1].next_out ^ 1;
+      }
+    }
+
+    //Relabels the "wall" of a face with a new face ID.
+    void wall_paint(Arc arc, int f_id, Arc ed) {
+      do {
+        arcs[arc.id].left_face = f_id;
+        turnRightF(arc);
+      } while (arc != ed);
+    }
+
+    //Relabels a component with a new component ID.
+    void component_relabel(Node node, int comp_id) {
+      std::vector<int> ns(nodes.size());
+      std::list<int> q;
+      q.push_back(node.id);
+      ns[node.id] = 1;
+      int cnt = 0;
+      while (!q.empty()) {
+        int n = q.front();
+        ns[n] = 2;
+        --components[nodes[n].component].num_nodes;
+        nodes[n].component = comp_id;
+        ++cnt;
+        q.pop_front();
+        Arc arc;
+        firstOut(arc,Node(n));
+        while (arc.id > -1) {
+          int m = arcs[arc.id].target;
+          if (ns[m] == 0 && nodes[m].component != comp_id) {
+            ns[m] = 1;
+            q.push_back(m);
+          }
+          nextOut(arc);
+        }
+      }
+      components[comp_id].num_nodes += cnt;
+    }
+
+    Face addFace() {
+      int n;
+
+      if (first_free_face==-1) {
+        n = faces.size();
+        faces.push_back(FaceT());
+      } else {
+        n = first_free_face;
+        first_free_face = faces[n].next;
+      }
+
+      faces[n].next = first_face;
+      if (first_face != -1) faces[first_face].prev = n;
+      first_face = n;
+      faces[n].prev = -1;
+
+      faces[n].first_arc = -1;
+
+      return Face(n);
+    }
+
+    void erase(const Face& face) {
+      int n = face.id;
+
+      if (faces[n].next != -1) {
+        faces[faces[n].next].prev = faces[n].prev;
+      }
+
+      if (faces[n].prev != -1) {
+        faces[faces[n].prev].next = faces[n].next;
+      } else {
+        first_face = faces[n].next;
+      }
+
+      faces[n].next = first_free_face;
+      first_free_face = n;
+      faces[n].prev = -2;
+    }
+
+    int addComponent() {
+      int n;
+
+      if (first_free_component==-1) {
+        n = components.size();
+        components.push_back(ComponentT());
+      } else {
+        n = first_free_component;
+        first_free_component = components[n].next;
+      }
+
+      components[n].next = first_component;
+      if (first_component != -1) components[first_component].prev = n;
+      first_component = n;
+      components[n].prev = -1;
+
+      components[n].num_nodes = 0;
+
+      return n;
+    }
+
+    void eraseComponent(const int n) {
+
+      if (components[n].next != -1) {
+        components[components[n].next].prev = components[n].prev;
+      }
+
+      if (components[n].prev != -1) {
+        components[components[n].prev].next = components[n].next;
+      } else {
+        first_component = components[n].next;
+      }
+
+      components[n].next = first_free_component;
+      first_free_component = n;
+      components[n].prev = -2;
+    }
+
+    int idComponent(Node n)
+    {
+      return nodes[n.id].component;
+    }
+
+#ifdef REMOVE_BEFORE_RELEASE
+  public:
+    void print() {
+      std::cout << "Nodes: " << std::endl;
+      for (int i=0; i<nodes.size(); ++i)
+        std::cout << i << ":"
+        << " fo=" << nodes[i].first_out
+        << " lo=" << nodes[i].last_out
+        << " pr=" << nodes[i].prev
+        << " nx=" << nodes[i].next
+        << " co=" << nodes[i].component
+        << " of=" << nodes[i].outer_face
+        <<std::endl;
+      std::cout << "Arcs: " << std::endl;
+      for (int i=0; i<arcs.size(); ++i) {
+        if (arcs[i].next_out > -2) {
+          std::cout << i << ":"
+          << " tg=" << arcs[i].target
+          << " po=" << arcs[i].prev_out
+          << " no=" << arcs[i].next_out
+          << " lf=" << arcs[i].left_face
+          <<std::endl;
+        } else std::cout << i << ": (deleted)" << std::endl;
+      }
+      std::cout << "Faces: " << std::endl;
+      for (int i=0; i<faces.size(); ++i) {
+        if (faces[i].prev > -2) {
+          std::cout << i
+            << " pr=" << faces[i].prev
+            << " nx=" << faces[i].next
+            << " fa=" << faces[i].first_arc
+            <<std::endl;
+        } else std::cout << i << ": (deleted)" << std::endl;
+      }
+      std::cout << "Components: " << std::endl;
+      for (int i=0; i<components.size(); ++i) {
+        std::cout << i << ':';
+        if (components[i].prev > -2) {
+          std::cout
+            << " pr=" << components[i].prev
+            << " nx=" << components[i].next
+            << " nn=" << components[i].num_nodes
+            <<std::endl;
+        } else std::cout << " (deleted)" << std::endl;
+      }
+    }
+#endif
+
+  };
+
+  // Face counting:
+
+  namespace _core_bits {
+
+    template <typename Graph, typename Enable = void>
+    struct CountFacesSelector {
+      static int count(const Graph &g) {
+        return countItems<Graph, typename Graph::Face>(g);
+      }
+    };
+
+    template <typename Graph>
+    struct CountFacesSelector<
+      Graph, typename
+      enable_if<typename Graph::FaceNumTag, void>::type>
+    {
+      static int count(const Graph &g) {
+        return g.faceNum();
+      }
+    };
+  }
+
+  /// \brief Function to count the faces in the graph.
+  ///
+  /// This function counts the faces in the graph.
+  /// The complexity of the function is <em>O</em>(<em>n</em>), but for some
+  /// graph structures it is specialized to run in <em>O</em>(1).
+  ///
+  /// \note If the graph contains a \c faceNum() member function and a
+  /// \c FaceNumTag tag then this function calls directly the member
+  /// function to query the cardinality of the face set.
+  template <typename Graph>
+  inline int countFaces(const Graph& g) {
+    return _core_bits::CountFacesSelector<Graph>::count(g);
+  }
+
+  template <typename Graph, typename DegIt>
+  inline int countFaceDegree(const Graph& _g, const typename Graph::Face& _n) {
+    int num = 0;
+    for (DegIt it(_g, _n); it != INVALID; ++it) {
+      ++num;
+    }
+    return num;
+  }
+  /// \brief Function to count the number of the boundary arcs of face \c f.
+  ///
+  /// This function counts the number of the boundary arcs of face \c f
+  /// in the undirected graph \c g.
+  template <typename Graph>
+  inline int countBoundaryArcs(const Graph& g,  const typename Graph::Face& f) {
+    return countFaceDegree<Graph, typename Graph::CwBoundaryArcIt>(g, f);
+  }
+
+
+  template <typename GR, typename Enable = void>
+  struct FaceNotifierIndicator {
+    typedef InvalidType Type;
+  };
+  template <typename GR>
+  struct FaceNotifierIndicator<
+    GR,
+    typename enable_if<typename GR::FaceNotifier::Notifier, void>::type
+  > {
+    typedef typename GR::FaceNotifier Type;
+  };
+
+  template <typename GR>
+  class ItemSetTraits<GR, typename GR::Face> {
+  public:
+
+    typedef GR Graph;
+    typedef GR Digraph;
+
+    typedef typename GR::Face Item;
+    typedef typename GR::FaceIt ItemIt;
+
+    typedef typename FaceNotifierIndicator<GR>::Type ItemNotifier;
+
+    template <typename V>
+    class Map : public GR::template FaceMap<V> {
+      typedef typename GR::template FaceMap<V> Parent;
+
+    public:
+      typedef typename GR::template FaceMap<V> Type;
+      typedef typename Parent::Value Value;
+
+      Map(const GR& _digraph) : Parent(_digraph) {}
+      Map(const GR& _digraph, const Value& _value)
+        : Parent(_digraph, _value) {}
+
+     };
+
+  };
+
+  template<typename Base>
+  class PlanarGraphExtender : public Base{
+
+    typedef Base Parent;
+
+  public:
+    typedef PlanarGraphExtender Graph;
+
+    typedef typename Parent::Node Node;
+    typedef typename Parent::Arc Arc;
+    typedef typename Parent::Edge Edge;
+    typedef typename Parent::Face Face;
+
+    typedef typename Parent::NodeNotifier NodeNotifier;
+    typedef typename Parent::EdgeNotifier EdgeNotifier;
+    typedef typename Parent::ArcNotifier ArcNotifier;
+    typedef AlterationNotifier<PlanarGraphExtender, Face> FaceNotifier;
+
+  protected:
+    mutable FaceNotifier face_notifier;
+
+  public:
+
+    int maxId(Face) const {
+      return Parent::maxFaceId();
+    }
+
+    NodeNotifier& notifier(Node) const {
+      return Parent::node_notifier;
+    }
+
+    ArcNotifier& notifier(Arc) const {
+      return Parent::arc_notifier;
+    }
+
+    EdgeNotifier& notifier(Edge) const {
+      return Parent::edge_notifier;
+    }
+
+    FaceNotifier& notifier(Face) const {
+      return face_notifier;
+    }
+
+    template <typename _Value>
+    class FaceMap
+    : public MapExtender<DefaultMap<Graph, Face, _Value> > {
+      typedef MapExtender<DefaultMap<Graph, Face, _Value> > Parent;
+
+    public:
+      explicit FaceMap(const Graph& graph)
+      : Parent(graph) {}
+      FaceMap(const Graph& graph, const _Value& value)
+      : Parent(graph, value) {}
+
+    private:
+      FaceMap& operator=(const FaceMap& cmap) {
+        return operator=<FaceMap>(cmap);
+      }
+
+      template <typename CMap>
+      FaceMap& operator=(const CMap& cmap) {
+        Parent::operator=(cmap);
+        return *this;
+      }
+
+    };
+
+
+  /// Iterator class for the faces.
+
+  /// This iterator goes through the faces of a planar graph.
+  class FaceIt : public Face {
+      const Graph* _graph;
+    public:
+
+      FaceIt() {}
+
+      FaceIt(Invalid i) : Face(i) { }
+
+      explicit FaceIt(const Graph& graph) : _graph(&graph) {
+        _graph->first(static_cast<Face&>(*this));
+      }
+
+      FaceIt(const Graph& graph, const Face& face)
+          : Face(face), _graph(&graph) {}
+
+      FaceIt& operator++() {
+        _graph->next(*this);
+        return *this;
+      }
+
+    };
+
+
+  /// Iterator class for the common faces of two nodes
+
+  /// This iterator goes through the common faces of two nodes
+  class CommonFacesIt : public Face {
+      const Graph* _graph;
+      const Node _n1, _n2;
+      std::set<Face> _f1, _f2;
+      typename std::set<Face>::const_iterator _fi2;
+    public:
+
+      CommonFacesIt() {}
+
+      CommonFacesIt(Invalid i) : Face(i) { }
+
+      explicit CommonFacesIt(const Graph& graph, Node n1, Node n2) : _graph(
+        &graph), _n1(n1), _n2(n2), _f1(), _f2(){
+        Arc _i1;
+        _graph->firstOut(_i1,_n1);
+        while (_i1 != INVALID) {
+          _f1.insert(_graph->leftFace(_i1));
+          _graph->nextOut(_i1);
+        }
+        Arc _i2;
+        _graph->firstOut(_i2,_n2);
+        while (_i2 != INVALID) {
+          _f2.insert(_graph->leftFace(_i2));
+          _graph->nextOut(_i2);
+        }
+
+        _fi2 = _f2.begin();
+        while (_fi2 != _f2.end() && _f1.count(*_fi2) == 0)
+          ++_fi2;
+        static_cast<Face&>(*this) = (_fi2 != _f2.end())?*_fi2:INVALID;
+      }
+
+      CommonFacesIt& operator++() {
+        ++_fi2;
+        while (_fi2 != _f2.end() && _f1.count(*_fi2) == 0)
+          ++_fi2;
+        static_cast<Face&>(*this) = (_fi2 != _f2.end())?*_fi2:INVALID;
+        return *this;
+      }
+
+    };
+
+  /// Iterator class for the common nodes of two facess
+
+  /// This iterator goes through the common nodes of two faces
+  class CommonNodesIt : public Node {
+      const Graph* _graph;
+      const Face _f1, _f2;
+      std::set<Node> _ns1,_ns2;
+      typename std::set<Node>::const_iterator _ni2;
+    public:
+
+      CommonNodesIt() {}
+
+      CommonNodesIt(Invalid i) : Face(i) { }
+
+      explicit CommonNodesIt(const Graph& graph, Face f1, Face f2) : _graph(
+        &graph), _f1(f1), _f2(f2), _ns1(), _ns2(){
+        Arc _i1;
+        _graph->firstCwF(_i1,_f1);
+        while (_i1 != INVALID) {
+          _ns1.insert(_graph->source(_i1));
+          _graph->nextCwF(_i1);
+        }
+        Arc _i2;
+        _graph->firstCwF(_i2,_f2);
+        while (_i2 != INVALID) {
+          _ns2.insert(_graph->source(_i2));
+          _graph->nextCwF(_i2);
+        }
+
+        _ni2 = _ns2.begin();
+        while (_ni2 != _ns2.end() && _ns1.count(*_ni2) == 0)
+          ++_ni2;
+        static_cast<Node&>(*this) = (_ni2 != _ns2.end())?*_ni2:INVALID;
+      }
+
+      CommonNodesIt& operator++() {
+        ++_ni2;
+        while (_ni2 != _ns2.end() && _ns1.count(*_ni2) == 0)
+          ++_ni2;
+        static_cast<Node&>(*this) = (_ni2 != _ns2.end())?*_ni2:INVALID;
+        return *this;
+      }
+
+    };
+
+  /// Iterator class for the arcs on the boundary of a face.
+
+  /// This iterator goes through the arcs on the boundary of a face clockwise.
+  class CwBoundaryArcIt : public Arc {
+      const Graph* _graph;
+      Face _face;
+      Arc f_arc;
+    public:
+
+      CwBoundaryArcIt() { }
+
+      CwBoundaryArcIt(Invalid i) : Arc(i) { }
+
+      CwBoundaryArcIt(const Graph& graph, const Face& face)
+          : _graph(&graph), _face(face) {
+        _graph->firstCwF(static_cast<Arc&>(*this),face);
+        f_arc = *this;
+      }
+
+      CwBoundaryArcIt(const Graph& graph, const Arc& arc)
+          : Arc(arc), _graph(&graph) {}
+
+      CwBoundaryArcIt& operator++() {
+        _graph->nextCwF(*this);
+        return *this;
+      }
+
+    };
+
+  /// Iterator class for the arcs around a node.
+
+  /// This iterator goes through the arcs around a node anticlockwise.
+  class CcwArcIt : public Arc {
+      const Graph* _graph;
+      const Node _node;
+    public:
+
+      CcwArcIt() { }
+
+      CcwArcIt(Invalid i) : Arc(i) { }
+
+      CcwArcIt(const Graph& graph, const Node& node)
+          : _graph(&graph), _node(node) {
+        _graph->firstCcw(*this, node);
+      }
+
+      CcwArcIt(const Graph& graph, const Arc& arc)
+          : Arc(arc), _graph(&graph) {}
+
+      CcwArcIt& operator++() {
+        _graph->nextCcw(*this, _node);
+        return *this;
+      }
+
+    };
+
+    void clear() {
+      notifier(Face()).clear();
+      Parent::clear();
+    }
+
+    template <typename Graph, typename NodeRefMap, typename EdgeRefMap>
+    void build(const Graph& graph, NodeRefMap& nodeRef,
+            EdgeRefMap& edgeRef) {
+      Parent::build(graph, nodeRef, edgeRef);
+      notifier(Face()).build();
+    }
+
+    PlanarGraphExtender() {
+      face_notifier.setContainer(*this);
+    }
+
+    ~PlanarGraphExtender() {
+      face_notifier.clear();
+    }
+
+  };
+
+  typedef PlanarGraphExtender<GraphExtender<PlanarGraphBase> >
+    ExtendedPlanarGraphBase;
+
+
+  /// \addtogroup graphs
+  /// @{
+
+  ///An undirected planar graph structure.
+
+  ///\ref PlanarGraph is a versatile and fast undirected planar graph
+  ///implementation based on linked lists that are stored in
+  ///\c std::vector structures. It maintains a topology of nodes, edges
+  ///and faces.
+  ///
+  ///This type fully conforms to the \ref concepts::Graph "Graph concept"
+  ///and it also provides several useful additional functionalities, including
+  ///those specific to planar graphs.
+  ///Most of its member functions and nested classes are documented
+  ///only in the concept class.
+  ///
+  ///This class provides only linear time counting for nodes, edges, arcs and
+  ///faces.
+  ///
+  ///A disconnected planar graph has have an outer face for each of its
+  ///components, effectively turning them into separate graphs. Each component
+  ///has a corresponding component in the dual.
+  ///\sa concepts::Graph
+  ///\sa PlaneGraph
+  class PlanarGraph : public ExtendedPlanarGraphBase {
+    typedef ExtendedPlanarGraphBase Parent;
+
+  public:
+    /// Graphs are \e not copy constructible. Use GraphCopy instead.
+    PlanarGraph(const PlanarGraph &) = delete;
+    /// \brief Assignment of a graph to another one is \e not allowed.
+    /// Use GraphCopy instead.
+    void operator=(const PlanarGraph &) = delete;
+
+    /// \brief Constructor
+
+    /// Constructor.
+    ///
+    PlanarGraph() {}
+
+
+    ///\brief Constructor that copies a planar embedding
+
+    ///Constructor that copies a planar embedding.
+    template<typename Graph>
+    PlanarGraph(const Graph &g, const lemon::PlanarEmbedding<Graph> &em) {
+      typename Graph::template NodeMap<Node> nm(g);
+      typename Graph::template ArcMap<Arc> am(g);
+      this->copyEmbedding(g,em,nm,am);
+    }
+
+    template<typename Graph>
+    PlanarGraph(const Graph &g, const lemon::PlanarEmbedding<Graph> &em,
+        typename Graph::template NodeMap<Node> &nm,
+        typename Graph::template ArcMap<Arc> &am) {
+      this->copyEmbedding(g,em,nm,am);
+    }
+
+    typedef Parent::OutArcIt IncEdgeIt;
+
+  protected:
+
+    template<typename Graph>
+    void copyEmbedding(const Graph &g, const lemon::PlanarEmbedding<Graph> &em,
+        typename Graph::template NodeMap<Node> &nm,
+        typename Graph::template ArcMap<Arc> &am) {
+      for (typename Graph::NodeIt it(g); it != INVALID; ++it) {
+        nm[it] = addNode();
+      }
+      for (typename Graph::ArcIt it(g); it != INVALID; ++it) {
+        am[it] = INVALID;
+      }
+      for (typename Graph::NodeIt it(g); it != INVALID; ++it) {
+        for (typename Graph::OutArcIt it2(g,it); it2 != INVALID; ++it2) {
+          if (am[it2] == INVALID) {
+            typename Graph::Arc p_u = em.next(it2);
+            while (am[p_u] == INVALID && it2 != p_u) {
+              p_u = em.next(p_u);
+            }
+            typename Graph::Arc ra = g.oppositeArc(it2);
+            typename Graph::Arc p_v = em.next(ra);
+            while (am[p_v] == INVALID && ra != p_v) {
+              p_v = em.next(p_v);
+            }
+            am[it2] = direct(addEdge(nm[g.source(it2)],nm[g.target(it2)],
+                    am[p_u],am[p_v]),nm[g.source(it2)]);
+            am[g.oppositeArc(it2)] = oppositeArc(am[it2]);
+          }
+        }
+      }
+    }
+
+    void edge_add_notify(Edge edge) {
+      notifier(Edge()).add(edge);
+      std::vector<Arc> ev;
+      ev.push_back(Parent::direct(edge, true));
+      ev.push_back(Parent::direct(edge, false));
+      notifier(Arc()).add(ev);
+    }
+
+    void edge_erase_notify(Edge edge) {
+      std::vector<Arc> av;
+      av.push_back(Parent::direct(edge, true));
+      av.push_back(Parent::direct(edge, false));
+      notifier(Arc()).erase(av);
+      notifier(Edge()).erase(edge);
+    }
+
+  public:
+    /// \brief Add a new node to the graph.
+    ///
+    /// This function adds a new node to the graph. A new outer face is created
+    /// for the node.
+    /// \return The new node.
+    Node addNode() {
+      Node n = PlanarGraphBase::addNode();
+      notifier(Node()).add(n);
+      notifier(Face()).add(outerFace(n));
+      return n;
+    }
+
+    /// \brief Add a new edge to the graph.
+    ///
+    /// This function adds a new edge to the graph between nodes
+    /// \c u and \c v with inherent orientation from node \c u to
+    /// node \c v. \c p_u and \c p_v are the edges that directly precede the new
+    /// edge in anticlockwise order at the nodes \c u and \c v, respectively.
+    /// INVALID should be passed as \c p_u or \c p_v if and only if the
+    /// respective node is isolated.
+    ///
+    /// If \c u and \c v are in the same component, \c p_u and \c p_v must share
+    /// the same left face (when looking from \c u or \c v). This face will be
+    /// split in two. If \c u and \c v are in different components, one of the
+    /// outer faces will be deleted.
+    /// \note It can be costly to join components unless one of them is an
+    /// isolated node.
+    /// \note Even if a face is not deleted in the process, BoundaryArcIt
+    /// iterators might run an incomplete lap or revisit previously passed edges
+    /// if edges are added to or removed from the boundary of the face.
+    /// \return The new edge, or INVALID if the parameters don't meet the
+    /// preconditions.
+    Edge addEdge(Node u, Node v, Edge p_u, Edge p_v) {
+      Face f_u = INVALID, f_v = INVALID;
+      if (p_u != INVALID) {
+        Arc a_u = direct(p_u,u);
+        f_u = leftFace(a_u);
+      }
+      if (p_v != INVALID) {
+        Arc a_v = direct(p_v,v);
+        f_v = leftFace(a_v);
+      }
+      Face o_u = outerFace(u);
+      Face o_v = outerFace(v);
+      Edge edge = PlanarGraphBase::addEdge(u, v, p_u, p_v);
+      if (edge == INVALID) return edge;
+      if (!valid(f_u)) notifier(Face()).erase(f_u);
+      if (!valid(f_v) && f_u != f_v) notifier(Face()).erase(f_v);
+      if (!valid(o_u)) notifier(Face()).erase(o_u);
+      if (!valid(o_v)) notifier(Face()).erase(o_v);
+      edge_add_notify(edge);
+      Face e_l = leftFace(direct(edge,u));
+      Face e_r = rightFace(direct(edge,u));
+      if (e_l != f_u && e_l != f_v && e_l != o_u && e_l != o_v)
+        notifier(Face()).add(e_l);
+      if (e_r != f_u && e_r != f_v && e_r != o_u && e_r != o_v && e_r != e_l)
+        notifier(Face()).add(e_r);
+      return edge;
+    }
+
+    ///\brief Erase a node from the graph.
+    ///
+    /// This function erases the given node along with its incident arcs
+    /// from the graph.
+    ///
+    /// \note All iterators referencing the removed node or the incident
+    /// edges are invalidated, of course.
+    void erase(Node n) {
+      notifier(Face()).erase(outerFace(n));
+      Parent::erase(n);
+    }
+
+    ///\brief Erase an edge from the graph.
+    ///
+    /// This function erases the given edge from the graph. The faces on the two
+    /// sides are merged into one, unless the edge holds two components
+    /// together. In the latter case, a new outer face is added and the
+    /// component is split in two.
+    ///
+    /// \note It can be costly to split a component, unless one of the resulting
+    /// components is an isolated node.
+    /// \note All iterators referencing the removed edge are invalidated,
+    /// of course.
+    /// \note Even if a face is not deleted in the process, BoundaryArcIt
+    /// iterators might run an incomplete lap or revisit previously passed edges
+    /// if edges are added to or removed from the boundary of the face.
+    void erase(Edge e) {
+      Node n1 = u(e);
+      Node n2 = v(e);
+      Face f_l = leftFace(direct(e,n1));
+      Face f_r = rightFace(direct(e,n1));
+      Parent::erase(e);
+      Face o_u = outerFace(n1);
+      Face o_v = outerFace(n2);
+      if (!valid(f_l)) notifier(Face()).erase(f_l);
+      if (!valid(f_r)) notifier(Face()).erase(f_r);
+      if (valid(o_u) && o_u != f_l && o_u != f_r)
+        notifier(Face()).add(o_u);
+      if (valid(o_v) && o_v != f_l && o_v != f_r)
+        notifier(Face()).add(o_v);
+    }
+    /// Node validity check
+
+    /// This function gives back \c true if the given node is valid,
+    /// i.e. it is a real node of the graph.
+    ///
+    /// \warning A removed node could become valid again if new nodes are
+    /// added to the graph.
+    bool valid(Node n) const {
+      return Parent::valid(n);
+    }
+    /// Edge validity check
+
+    /// This function gives back \c true if the given edge is valid,
+    /// i.e. it is a real edge of the graph.
+    ///
+    /// \warning A removed edge could become valid again if new edges are
+    /// added to the graph.
+    bool valid(Edge e) const {
+      return Parent::valid(e);
+    }
+    /// Arc validity check
+
+    /// This function gives back \c true if the given arc is valid,
+    /// i.e. it is a real arc of the graph.
+    ///
+    /// \warning A removed arc could become valid again if new edges are
+    /// added to the graph.
+    bool valid(Arc a) const {
+      return Parent::valid(a);
+    }
+
+    /// Face validity check
+
+    /// This function gives back \c true if the given face is valid,
+    /// i.e. it is a real face of the graph.
+    ///
+    /// \warning A removed face could become valid again if new faces are
+    /// added to the graph.
+    bool valid(Face f) const {
+      return Parent::valid(f);
+    }
+
+    /// \brief Change the first node of an edge.
+    ///
+    /// Planar graphs don't support changing the endpoints of edges.
+    void changeU(Edge e, Node n) = delete;
+    /// \brief Change the second node of an edge.
+    ///
+    /// Planar graphs don't support changing the endpoints of edges.
+    void changeV(Edge e, Node n) = delete;
+
+    /// \brief Contract two nodes.
+    ///
+    /// This function contracts the two ends of the given edge.
+    /// One of the nodes on \c e is removed, but instead of deleting
+    /// its incident edges, they are joined to the other node.
+    /// \c e will be deleted.
+    ///
+    /// The deleted node is v(e) unless the degree of v(e) is higher than 1 and
+    /// the degree of u(e) is 1, in which case u(e) is deleted.
+    ///
+    /// \note All edge and arc iterators whose base node is
+    /// the deleted one are invalidated.
+    /// Moreover all iterators referencing the removed node or
+    /// edge are also invalidated. Other iterators remain valid. However, they
+    /// might run an incomplete lap or revisit previously passed edges if
+    /// incremented.
+    ///
+    ///\warning This functionality cannot be used together with the
+    ///Snapshot feature.
+    void contract(Edge e) {
+      Node n1 = u(e);
+      Node n2 = v(e);
+      if (n1 == n2) return;
+      Node nd;
+      bool simple;
+      inner_contract1(n1,n2,nd,simple);
+      notifier(Node()).erase(nd);
+      edge_erase_notify(e);
+      inner_contract2(e,n1,n2,nd,simple);
+    }
+
+    /// \brief Split an edge.
+    ///
+    ///This function splits the given edge. First, a new node \c v is
+    ///added to the graph, then the target node of the original edge
+    ///is set to \c v. Finally, an edge from \c v to the original target
+    ///is added.
+    ///\return The newly created node.
+    ///
+    ///\note Iterators referencing the original edge are
+    ///invalidated. Other iterators remain valid.
+    ///
+    ///\warning This functionality cannot be used together with the
+    ///Snapshot feature.
+    Node split(Edge e) {
+      Node v;
+      Edge e2;
+      inner_split1(v,e2);
+      notifier(Node()).add(v);
+      edge_add_notify(e2);
+      return inner_split2(e,v,e2);
+    }
+
+    ///Split a node.
+
+    ///This function splits the given node. First, a new node is added
+    ///to the graph, then all edges in anticlockwise order from \c e1 until but
+    ///not including \c e2
+    ///are re-anchored from \c n to the new node.
+    ///If the second parameter \c connect is \c true (this is the default
+    ///value), then a new edge from node \c n to the newly created node
+    ///is also added. \c e1 and \c e2 must be INVALID if and only if \c n is
+    ///isolated. If \c connect is false, the component of \c n might be split in
+    ///two. It can be costly to split components unless one of the resulting
+    ///components is an isolated node.
+    ///\return The newly created node.
+    ///
+    ///\note All iterators remain valid but some might run an incomplete lap or
+    ///revisit previously passed nodes when incremented.
+    ///
+    ///\warning This functionality cannot be used together with the
+    ///Snapshot feature.
+    Node split(Node n1, Edge e1, Edge e2, bool connect) {
+      Node n2;
+      inner_split1(n1,n2);
+      notifier(Node()).add(n2);
+      Edge a3;
+      inner_split2(n1,n2,e1,e2,connect,a3);
+      if (!connect) {
+        if (a3 != INVALID)
+          erase(a3);
+      } else {
+        edge_add_notify(a3);
+      }
+      return n2;
+    }
+
+
+    ///Clear the graph.
+
+    ///This function erases all nodes and arcs from the graph.
+    ///
+    ///\note All iterators of the graph are invalidated, of course.
+    void clear() {
+      Parent::clear();
+    }
+
+    /// Reserve memory for nodes.
+
+    /// Using this function, it is possible to avoid superfluous memory
+    /// allocation: if you know that the graph you want to build will
+    /// be large (e.g. it will contain millions of nodes and/or edges),
+    /// then it is worth reserving space for this amount before starting
+    /// to build the graph.
+    /// \sa reserveEdge()
+    void reserveNode(int n) {
+      nodes.reserve(n);
+    };
+
+    /// Reserve memory for edges.
+
+    /// Using this function, it is possible to avoid superfluous memory
+    /// allocation: if you know that the graph you want to build will
+    /// be large (e.g. it will contain millions of nodes and/or edges),
+    /// then it is worth reserving space for this amount before starting
+    /// to build the graph.
+    /// \sa reserveNode()
+    void reserveEdge(int m) {
+      arcs.reserve(2 * m);
+    };
+
+    /// Reserve memory for faces.
+
+    /// Using this function, it is possible to avoid superfluous memory
+    /// allocation: if you know that the graph you want to build will
+    /// be large (e.g. it will contain millions of nodes and/or edges),
+    /// then it is worth reserving space for this amount before starting
+    /// to build the graph.
+    /// \sa reserveFace()
+    void reserveFace(int n) {
+      faces.reserve(n);
+    };
+
+    class DualBase {
+      const Graph *_graph;
+    protected:
+      void initialize(const Graph *graph) {
+        _graph = graph;
+      }
+    public:
+
+      typedef PlanarGraph::Face Node;
+      typedef PlanarGraph::Arc Arc;
+      typedef PlanarGraph::Edge Edge;
+      typedef PlanarGraph::Node Face;
+
+      int maxNodeId() const {
+        return _graph->maxFaceId();
+      }
+      int maxArcId() const {
+        return _graph->maxArcId();
+      }
+      int maxFaceId() const {
+        return _graph->maxNodeId();
+      }
+
+      Node source(Arc e) const {
+        return _graph->leftFace(e);
+      }
+      Node target(Arc e) const {
+        return _graph->rightFace(e);
+      }
+      Face leftFace(Arc e) const {
+        return _graph->target(e);
+      }
+      Face rightFace(Arc e) const {
+        return _graph->source(e);
+      }
+      Arc direct(const Edge &edge, const Node &node) const {
+        return _graph->direct(edge, _graph->w1(edge) == node);
+      }
+
+      void first(Node &i) const {
+        _graph->first(i);
+      }
+      void next(Node &i) const {
+        _graph->next(i);
+      }
+      void first(Arc &i) const {
+        _graph->first(i);
+      }
+      void next(Arc &i) const {
+        _graph->next(i);
+      }
+      void firstCcw(Arc& i, const Node& n) const {
+        _graph->lastInF(i, n);
+      }
+      void nextCcw(Arc& i, const Node &n) const {
+        _graph->prevInF(i);
+      }
+      void firstIn(Arc& i, const Node& n) const {
+        _graph->firstInF(i, n);
+      }
+      void nextIn(Arc& i) const {
+        _graph->nextInF(i);
+      }
+      void firstCwF(Arc& i, const Face& n) const {
+        _graph->lastIn(i, n);
+      }
+      void nextCwF(Arc& i) const {
+        _graph->prevIn(i);
+      }
+      void firstOut(Arc& i, const Node& n ) const {
+        _graph->firstOutF(i, n);
+      }
+      void nextOut(Arc& i) const {
+        _graph->nextOutF(i);
+      }
+      void first(Face &i) const {
+        _graph->first(i);
+      }
+      void next(Face &i) const {
+        _graph->next(i);
+      }
+
+      static int id(Node v) {
+        return PlanarGraph::id(v);
+      }
+      static int id(Arc e) {
+        return PlanarGraph::id(e);
+      }
+      static int id(Face f) {
+        return PlanarGraph::id(f);
+      }
+      static Node nodeFromId(int id) {
+        return PlanarGraph::faceFromId(id);
+      }
+      static Arc arcFromId(int id) {
+        return PlanarGraph::arcFromId(id);
+      }
+      static Face faceFromId(int id) {
+        return PlanarGraph::nodeFromId(id);
+      }
+
+      bool valid(Node n) const {
+        return _graph->valid(n);
+      }
+      bool valid(Arc n) const {
+        return _graph->valid(n);
+      }
+      bool valid(Face n) const {
+        return _graph->valid(n);
+      }
+
+    };
+
+    typedef PlanarGraphExtender<GraphExtender<DualBase> > ExtendedDualBase;
+
+  /// Adaptor class for the dual of a planar graph.
+
+  /// This is an adaptor class for the dual of a planar graph.
+  class Dual : public ExtendedDualBase {
+    public:
+      Dual(const PlanarGraph &graph) {
+        initialize(&graph);
+      }
+
+    };
+    /// \brief Class to make a snapshot of the graph and restore
+    /// it later.
+    ///
+    /// Class to make a snapshot of the graph and restore it later.
+    ///
+    /// The newly added nodes and edges can be removed
+    /// using the restore() function.
+    ///
+    /// \note After a state is restored, you cannot restore a later state,
+    /// i.e. you cannot add the removed nodes and edges again using
+    /// another Snapshot instance.
+    ///
+    /// \warning Node and edge deletions and other modifications
+    /// (e.g. changing the end-nodes of edges or contracting nodes)
+    /// cannot be restored. These events invalidate the snapshot.
+    /// However, the edges and nodes that were added to the graph after
+    /// making the current snapshot can be removed without invalidating it.
+    class Snapshot {
+    protected:
+
+      typedef Parent::NodeNotifier NodeNotifier;
+
+    class NodeObserverProxy : public NodeNotifier::ObserverBase {
+      public:
+
+        NodeObserverProxy(Snapshot& _snapshot)
+            : snapshot(_snapshot) {}
+
+        using NodeNotifier::ObserverBase::attach;
+        using NodeNotifier::ObserverBase::detach;
+        using NodeNotifier::ObserverBase::attached;
+
+      protected:
+
+        virtual void add(const Node& node) {
+          snapshot.addNode(node);
+        }
+        virtual void add(const std::vector<Node>& nodes) {
+          for (int i = nodes.size() - 1; i >= 0; ++i) {
+            snapshot.addNode(nodes[i]);
+          }
+        }
+        virtual void erase(const Node& node) {
+          snapshot.eraseNode(node);
+        }
+        virtual void erase(const std::vector<Node>& nodes) {
+          for (int i = 0; i < int(nodes.size()); ++i) {
+            snapshot.eraseNode(nodes[i]);
+          }
+        }
+        virtual void build() {
+          Node node;
+          std::vector<Node> nodes;
+          for (notifier()->first(node); node != INVALID;
+               notifier()->next(node)) {
+            nodes.push_back(node);
+          }
+          for (int i = nodes.size() - 1; i >= 0; --i) {
+            snapshot.addNode(nodes[i]);
+          }
+        }
+        virtual void clear() {
+          Node node;
+          for (notifier()->first(node); node != INVALID;
+               notifier()->next(node)) {
+            snapshot.eraseNode(node);
+          }
+        }
+
+        Snapshot& snapshot;
+      };
+
+    class EdgeObserverProxy : public EdgeNotifier::ObserverBase {
+      public:
+
+        EdgeObserverProxy(Snapshot& _snapshot)
+            : snapshot(_snapshot) {}
+
+        using EdgeNotifier::ObserverBase::attach;
+        using EdgeNotifier::ObserverBase::detach;
+        using EdgeNotifier::ObserverBase::attached;
+
+      protected:
+
+        virtual void add(const Edge& edge) {
+          snapshot.addEdge(edge);
+        }
+        virtual void add(const std::vector<Edge>& edges) {
+          for (int i = edges.size() - 1; i >= 0; ++i) {
+            snapshot.addEdge(edges[i]);
+          }
+        }
+        virtual void erase(const Edge& edge) {
+          snapshot.eraseEdge(edge);
+        }
+        virtual void erase(const std::vector<Edge>& edges) {
+          for (int i = 0; i < int(edges.size()); ++i) {
+            snapshot.eraseEdge(edges[i]);
+          }
+        }
+        virtual void build() {
+          Edge edge;
+          std::vector<Edge> edges;
+          for (notifier()->first(edge); edge != INVALID;
+               notifier()->next(edge)) {
+            edges.push_back(edge);
+          }
+          for (int i = edges.size() - 1; i >= 0; --i) {
+            snapshot.addEdge(edges[i]);
+          }
+        }
+        virtual void clear() {
+          Edge edge;
+          for (notifier()->first(edge); edge != INVALID;
+               notifier()->next(edge)) {
+            snapshot.eraseEdge(edge);
+          }
+        }
+
+        Snapshot& snapshot;
+      };
+
+      PlanarGraph *graph;
+
+      NodeObserverProxy node_observer_proxy;
+      EdgeObserverProxy edge_observer_proxy;
+
+      std::list<Node> added_nodes;
+      std::list<Edge> added_edges;
+
+
+      void addNode(const Node& node) {
+        added_nodes.push_front(node);
+      }
+      void eraseNode(const Node& node) {
+        std::list<Node>::iterator it =
+          std::find(added_nodes.begin(), added_nodes.end(), node);
+        if (it == added_nodes.end()) {
+          clear();
+          edge_observer_proxy.detach();
+          throw NodeNotifier::ImmediateDetach();
+        } else {
+          added_nodes.erase(it);
+        }
+      }
+
+      void addEdge(const Edge& edge) {
+        added_edges.push_front(edge);
+      }
+      void eraseEdge(const Edge& edge) {
+        std::list<Edge>::iterator it =
+          std::find(added_edges.begin(), added_edges.end(), edge);
+        if (it == added_edges.end()) {
+          clear();
+          node_observer_proxy.detach();
+          throw EdgeNotifier::ImmediateDetach();
+        } else {
+          added_edges.erase(it);
+        }
+      }
+
+      void attach(PlanarGraph &_graph) {
+        graph = &_graph;
+        node_observer_proxy.attach(graph->notifier(Node()));
+        edge_observer_proxy.attach(graph->notifier(Edge()));
+      }
+
+      void detach() {
+        node_observer_proxy.detach();
+        edge_observer_proxy.detach();
+      }
+
+      bool attached() const {
+        return node_observer_proxy.attached();
+      }
+
+      void clear() {
+        added_nodes.clear();
+        added_edges.clear();
+      }
+
+    public:
+
+      /// \brief Default constructor.
+      ///
+      /// Default constructor.
+      /// You have to call save() to actually make a snapshot.
+      Snapshot()
+          : graph(0), node_observer_proxy(*this),
+          edge_observer_proxy(*this) {}
+
+      /// \brief Constructor that immediately makes a snapshot.
+      ///
+      /// This constructor immediately makes a snapshot of the given graph.
+      Snapshot(PlanarGraph &gr)
+          : node_observer_proxy(*this),
+          edge_observer_proxy(*this) {
+        attach(gr);
+      }
+
+      /// \brief Make a snapshot.
+      ///
+      /// This function makes a snapshot of the given graph.
+      /// It can be called more than once. In case of a repeated
+      /// call, the previous snapshot gets lost.
+      void save(PlanarGraph &gr) {
+        if (attached()) {
+          detach();
+          clear();
+        }
+        attach(gr);
+      }
+
+      /// \brief Undo the changes until the last snapshot.
+      ///
+      /// This function undos the changes until the last snapshot
+      /// created by save() or Snapshot(PlanarGraph&).
+      ///
+      /// \warning This method invalidates the snapshot, i.e. repeated
+      /// restoring is not supported unless you call save() again.
+      void restore() {
+        detach();
+        for (std::list<Edge>::iterator it = added_edges.begin();
+             it != added_edges.end(); ++it) {
+          graph->erase(*it);
+        }
+        for (std::list<Node>::iterator it = added_nodes.begin();
+             it != added_nodes.end(); ++it) {
+          graph->erase(*it);
+        }
+        clear();
+      }
+
+      /// \brief Returns \c true if the snapshot is valid.
+      ///
+      /// This function returns \c true if the snapshot is valid.
+      bool valid() const {
+        return attached();
+      }
+    };
+  };
+
+  /// @}
+} //namespace lemon
+
+
+#endif

new file lemon/plane_graph.h

@@ +1 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2010
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_PLANE_GRAPH_H
+#define LEMON_PLANE_GRAPH_H
+
+///\ingroup graphs
+///\file
+///\brief PlaneGraph classes. UNDER CONSTRUCTION.
+
+#include <lemon/planarity.h>
+#include <lemon/planar_graph.h>
+#include <lemon/dim2.h>
+#include <cmath>
+#include <lemon/math.h>
+#include <map>
+
+#include "bits\graph_extender.h"
+#include "bits\bezier.h"
+
+
+namespace lemon {
+
+  template<typename Point>
+  struct LexicographicPointOrdering {
+    bool operator()(const Point &p1, const Point &p2) {
+      return p1.x < p2.x || (p1.x == p2.x && p1.y < p2.y);
+    }
+  };
+
+  ///Traits class for straight edges in a plane graph.
+
+  ///This class represents straight edges. It can be used as a template
+  ///parameter to \ref PlaneGraph.
+  class PlaneGraphStraightEdge {
+
+  private:
+    struct Data {
+      dim2::Point<double> p1, p2;
+      double angle;
+    };
+
+    static const double epsilon = 0.00000000000001;
+
+    static void getEquation(const Data &d, double &a, double &b, double &c) {
+      double vx = d.p2.x-d.p1.x, vy = d.p2.y-d.p1.y;
+      a = vy;
+      b = -vx;
+      c = vy*d.p1.x-vx*d.p1.y;
+    }
+
+    static bool hasPoint(const Data &d, dim2::Point<double> p, bool os, bool oe
+      ) {
+      if (!os && d.p1.x-epsilon < p.x && p.x < d.p1.x+epsilon &&
+        d.p1.y-epsilon < p.y && p.y < d.p1.y+epsilon) return true;
+      if (!oe && d.p2.x-epsilon < p.x && p.x < d.p2.x+epsilon &&
+        d.p2.y-epsilon < p.y && p.y < d.p2.y+epsilon) return true;
+      double x1,x2,y1,y2;
+      d.p1.x <= d.p2.x ? (x1 = d.p1.x, x2 = d.p2.x) : (x2 = d.p1.x, x1 =
+        d.p2.x);
+      d.p1.y <= d.p2.y ? (y1 = d.p1.y, y2 = d.p2.y) : (y2 = d.p1.y, y1 =
+        d.p2.y);
+      return ((x1+epsilon < p.x && p.x < x2-epsilon) ||
+             (x1-epsilon <= p.x && p.x <= x1+epsilon &&
+              x2-epsilon <= p.x && p.x <= x2+epsilon)) &&
+             ((y1+epsilon < p.y && p.y < y2-epsilon) ||
+             (y1-epsilon <= p.y && p.y <= y1+epsilon &&
+              y2-epsilon <= p.y && p.y <= y2+epsilon));
+    }
+
+  public:
+
+    typedef Data Shape;
+    typedef dim2::Point<double> Point;
+
+    static Data create() { return Data(); }
+    static Data create(dim2::Point<double> p1_, dim2::Point<double> p2_) {
+      Data d = {p1_, p2_, atan2(p2_.y-p1_.y,p2_.x-p1_.x)};
+      return d;
+    }
+
+    static Point source(const Data &d) {
+      return d.p1;
+    }
+
+    static Point target(const Data &d) {
+      return d.p2;
+    }
+
+    //Determine if c comes between a and b or not, in positive order
+    static bool between(const Data &c, const Data &a, const Data &b)
+    {
+      double aa = a.angle, ba = b.angle, ca = c.angle;
+      if (ba < aa) ba += 2*PI;
+      if (ca < aa) ca += 2*PI;
+      return ba > ca;
+    }
+
+    static Data reverse(const Data &d) {
+      return create(d.p2, d.p1);
+    }
+
+    static bool match(const Data &a, const Data &b) {
+      return (a.p1 == b.p1 && a.p2 == b.p2) || (a.p1 == b.p2 && a.p2 == b.p1);
+    }
+
+    static bool intersect(const Data &q, const Data &e, bool qos = true,
+      bool qoe = true, bool eos = true, bool eoe = true) {
+      if (match(q,e)) return true;
+      double a1,b1,c1,a2,b2,c2;
+      getEquation(q,a1,b1,c1);
+      getEquation(e,a2,b2,c2);
+      double d = b1*a2-b2*a1;
+      if (d == 0) {
+        if (a2*q.p1.x+b2*q.p1.y == c2) {
+          return hasPoint(q,e.p1,qos,qoe) || hasPoint(q,e.p2,qos,qoe) ||
+            hasPoint(e,q.p1,eos,eoe) || hasPoint(e,q.p2,eos,eoe);
+        } else {
+          return false;
+        }
+      }
+      double x = (c2*b1-c1*b2)/(d);
+      double y = (c1*a2-c2*a1)/(d);
+      dim2::Point<double> pt(x,y);
+      return hasPoint(q,pt,qos,qoe) && hasPoint(e,pt,eos,eoe);
+    }
+
+    static bool degenerate(const Data &d) {
+      return d.p1 == d.p2;
+    }
+
+    static void toEps(std::ostream &os, const Data &d, const Color &col,
+      double width) {
+      os << d.p1.x << ' '
+         << d.p1.y << ' '
+         << d.p2.x << ' '
+         << d.p2.y << ' '
+         << col.red() << ' '
+         << col.green() << ' '
+         << col.blue() << ' '
+         << width << " l\n";
+    }
+  };
+
+  ///Class for a polyline edge in a plane graph.
+
+  ///This class represents a polyline edge. It can be used as a template
+  ///parameter to \ref PlaneGraph.
+  class PlaneGraphLineStripEdge {
+  private:
+    struct Data {
+      std::vector<dim2::Point<double> > ps;
+      double angle;
+    };
+
+  public:
+
+    typedef Data Shape;
+    typedef dim2::Point<double> Point;
+
+    static Data create() { return Data(); }
+    static Data create(const std::vector<Point> &ps_) {
+      Data d;
+      d.ps = ps_;
+      if (ps_.size() > 1) {
+        const Point &p1 = ps_[0], &p2 = ps_[1];
+        d.angle = atan2(p2.y-p1.y,p2.x-p1.x);
+      }
+      return d;
+    }
+
+    static Point source(const Data &d) {
+      return d.ps.front();
+    }
+
+    static Point target(const Data &d) {
+      return d.ps.back();
+    }
+
+    //Determine if c comes between a and b or not, in positive order
+    static bool between(const Data &c, const Data &a, const Data &b)
+    {
+      double aa = a.angle, ba = b.angle, ca = c.angle;
+      if (ba < aa) ba += 2*PI;
+      if (ca < aa) ca += 2*PI;
+      return ba > ca;
+    }
+
+    static Data reverse(const Data &d) {
+      std::vector<Point> e(d.ps.size());
+      reverse_copy(d.ps.begin(), d.ps.end(), e.begin());
+      return create(e);
+    }
+
+    static bool intersect(const Data &q, const Data &e) {
+      for (int i=0; i<q.ps.size()-1; ++i) {
+        PlaneGraphStraightEdge::Shape s1 { q.ps[i], q.ps[i+1] };
+        for (int j=0; j<e.ps.size()-1; ++j) {
+          PlaneGraphStraightEdge::Shape s2 { e.ps[j], e.ps[j+1] };
+          if (PlaneGraphStraightEdge::intersect(s1,s2,i==0,i==q.ps.size()-2,
+            j==0,j==e.ps.size()-2)) return true;
+        }
+      }
+      return false;
+    }
+
+    static bool degenerate(const Data &d) {
+      if (d.ps.size() < 2) return true;
+      for (int i=0; i<d.ps.size(); ++i)
+        for (int j=i+1; j<d.ps.size(); ++j) {
+          if (d.ps[i] == d.ps[j] && (i > 0 || j<d.ps.size()-1)) return true;
+        }
+
+      for (int i=0; i<d.ps.size()-1; ++i) {
+        PlaneGraphStraightEdge::Shape s1 { d.ps[i], d.ps[i+1] };
+        for (int j=i+1; j<d.ps.size()-1; ++j) {
+          PlaneGraphStraightEdge::Shape s2 { d.ps[j], d.ps[j+1] };
+          if (PlaneGraphStraightEdge::intersect(s1,s2,true,true,true,true))
+            return true;
+        }
+      }
+
+      return false;
+    }
+
+    static void toEps(std::ostream &os, const Data &d, const Color &col,
+      double width) {
+      for (int i=0; i<d.ps.size()-1; ++i) {
+        os << d.ps[i].x << ' '
+           << d.ps[i].y << ' '
+           << d.ps[i+1].x << ' '
+           << d.ps[i+1].y << ' '
+           << col.red() << ' '
+           << col.green() << ' '
+           << col.blue() << ' '
+           << width << " l\n";
+      }
+  };
+
+  };
+
+  ///Class for a B&eacute;zier curve edge in a plane graph.
+
+  ///This class represents a B&eacute;zier curve edge of degree 2. It can be
+  ///used as a
+  ///template parameter to \ref PlaneGraph.
+  class PlaneGraphBezier2Edge {
+
+  private:
+    struct Data {
+      dim2::Bezier2 ps;
+      double angle2, angle1;
+    };
+
+  public:
+
+    typedef Data Shape;
+    typedef dim2::Point<double> Point;
+
+    static Data create() { return Data(); }
+    static Data create(const Point &p1, const Point &p2, const Point &p3) {
+      Data d;
+      d.ps = dim2::Bezier2(p1,p2,p3);
+      d.angle2 = atan2(p2.y-p1.y,p2.x-p1.x);
+      dim2::Bezier1 pps = d.ps.grad();
+      d.angle1 = atan2(pps.p2.y-pps.p1.y,pps.p2.x-pps.p1.x);
+      return d;
+    }
+
+    static Point source(const Data &d) {
+      return d.ps.p1;
+    }
+
+    static Point target(const Data &d) {
+      return d.ps.p3;
+    }
+
+    //Determine if c comes between a and b or not, in positive order
+    static bool between(const Data &c, const Data &a, const Data &b)
+    {
+      double aa = a.angle2, ba = b.angle2, ca = c.angle2;
+      if (ba < aa) ba += 2*PI;
+      if (ca < aa) ca += 2*PI;
+      if (ba > ca) return true;
+      else if (ba < ca) return false;
+      else {
+        double aa = a.angle1, ba = b.angle1, ca = c.angle1;
+        if (ba < aa) ba += 2*PI;
+        if (ca < aa) ca += 2*PI;
+        if (ba > ca) return true;
+        return ba > ca;
+      }
+    }
+
+    static Data reverse(const Data &d) {
+      return create(d.ps.p3,d.ps.p2,d.ps.p1);
+    }
+
+    static bool intersect(const Data &q, const Data &e) {
+      dim2::Box<double> box1 = q.ps.boundingBox();
+      dim2::Box<double> box2 = q.ps.boundingBox();
+      if ((box1 & box2).empty()) return false;
+      std::vector<Point> v1, v2;
+      for (int i=0; i<=100; ++i) {
+        v1.push_back(q.ps(1.0*i/100));
+        v2.push_back(e.ps(1.0*i/100));
+      }
+      return PlaneGraphLineStripEdge::intersect(
+              PlaneGraphLineStripEdge::create(v1),
+              PlaneGraphLineStripEdge::create(v2));
+    }
+
+    static bool degenerate(const Data &d) {
+      return d.ps.p1 == d.ps.p2 || d.ps.p2 == d.ps.p3 || d.ps.p1 == d.ps.p3;
+    }
+
+    static void toEps(std::ostream &os, const Data &d, const Color &col,
+      double width) {
+      dim2::Bezier3 b = d.ps;
+      os << width << " setlinewidth "
+         << col.red() << ' '
+         << col.green() << ' '
+         << col.blue() << " setrgbcolor newpath\n"
+         << b.p1.x << ' ' << b.p1.y << " moveto\n"
+         << b.p2.x << ' ' << b.p2.y << ' '
+         << b.p3.x << ' ' << b.p3.y << ' '
+         << b.p4.x << ' ' << b.p4.y << " curveto stroke\n";
+    }
+  };
+
+  ///Class for a B&eacute;zier curve edge in a plane graph.
+
+  ///This class represents a B&eacute;zier curve edge of degree 3. It can be
+  ///used as a
+  ///template parameter to \ref PlaneGraph.
+  class PlaneGraphBezier3Edge {
+
+  private:
+    struct Data {
+      dim2::Bezier3 ps;
+      double angle3, angle2, angle1;
+    };
+
+  public:
+
+    typedef Data Shape;
+    typedef dim2::Point<double> Point;
+
+    static Data create() { return Data(); }
+    static Data create(const Point &p1, const Point &p2, const Point &p3,
+        const Point &p4) {
+      Data d;
+      d.ps = dim2::Bezier3(p1,p2,p3,p4);
+      d.angle3 = atan2(p2.y-p1.y,p2.x-p1.x);
+      dim2::Bezier2 pps = d.ps.grad();
+      d.angle2 = atan2(pps.p2.y-pps.p1.y,pps.p2.x-pps.p1.x);
+      dim2::Bezier1 ppps = pps.grad();
+      d.angle3 = atan2(ppps.p2.y-ppps.p1.y,ppps.p2.x-ppps.p1.x);
+      return d;
+    }
+
+    static Point source(const Data &d) {
+      return d.ps.p1;
+    }
+
+    static Point target(const Data &d) {
+      return d.ps.p4;
+    }
+
+    //Determine if c comes between a and b or not, in positive order
+    static bool between(const Data &c, const Data &a, const Data &b)
+    {
+      double aa = a.angle3, ba = b.angle3, ca = c.angle3;
+      if (ba < aa) ba += 2*PI;
+      if (ca < aa) ca += 2*PI;
+      if (ba > ca) return true;
+      else if (ba < ca) return false;
+      else {
+        double aa = a.angle2, ba = b.angle2, ca = c.angle2;
+        if (ba < aa) ba += 2*PI;
+        if (ca < aa) ca += 2*PI;
+        if (ba > ca) return true;
+        else if (ba < ca) return false;
+        else {
+          double aa = a.angle1, ba = b.angle1, ca = c.angle1;
+          if (ba < aa) ba += 2*PI;
+          if (ca < aa) ca += 2*PI;
+          return (ba > ca);
+        }
+      }
+    }
+
+    static Data reverse(const Data &d) {
+      return create(d.ps.p4,d.ps.p3,d.ps.p2,d.ps.p1);
+    }
+
+    static bool intersect(const Data &q, const Data &e) {
+      dim2::Box<double> box1 = q.ps.boundingBox();
+      dim2::Box<double> box2 = q.ps.boundingBox();
+      if ((box1 & box2).empty()) return false;
+      std::vector<Point> v1, v2;
+      for (int i=0; i<=100; ++i) {
+        v1.push_back(q.ps(1.0*i/100));
+        v2.push_back(e.ps(1.0*i/100));
+      }
+      return PlaneGraphLineStripEdge::intersect(
+              PlaneGraphLineStripEdge::create(v1),
+              PlaneGraphLineStripEdge::create(v2));
+    }
+
+    static bool degenerate(const Data &d) {
+      if (d.ps.p1 == d.ps.p2 || d.ps.p3 == d.ps.p4)
+        return true;
+      return PlaneGraphStraightEdge::intersect(
+              PlaneGraphStraightEdge::create(d.ps.p1,d.ps.p2),
+              PlaneGraphStraightEdge::create(d.ps.p3,d.ps.p4),
+              true,true,true,true);
+    }
+
+    static void toEps(std::ostream &os, const Data &d, const Color &col,
+      double width) {
+      Point p1 = d.ps.p1;
+/*      if (d.ps.p1 == d.ps.p2) {
+        p1 = d.ps(0.5);
+      }*/
+      os << width << " setlinewidth "
+         << col.red() << ' '
+         << col.green() << ' '
+         << col.blue() << " setrgbcolor newpath\n"
+         << p1.x << ' ' << p1.y << " moveto\n"
+         << d.ps.p2.x << ' ' << d.ps.p2.y << ' '
+         << d.ps.p3.x << ' ' << d.ps.p3.y << ' '
+         << d.ps.p4.x << ' ' << d.ps.p4.y << " curveto stroke\n";
+    }
+  };
+
+  /// \addtogroup graphs
+  /// @{
+
+  ///A data structure for a graph embedded in the plane.
+
+  ///\ref PlaneGraph adds geometry features to PlanarGraph. While PlanarGraph
+  ///must be built from the topology, PlaneGraph requires coordinates for nodes
+  ///and determines the topology that fits the coordinates.
+  ///\sa PlanarGraph
+  template <typename EdgeType = PlaneGraphStraightEdge>
+  class PlaneGraph : public PlanarGraph {
+
+  private:
+
+    typedef typename EdgeType::Shape EdgeShape;
+    typedef typename EdgeType::Point EdgePoint;
+
+    NodeMap<EdgePoint> nodepos;
+    ArcMap<EdgeShape> edgeshapes;
+    
+    typedef std::map<EdgePoint, Node, LexicographicPointOrdering<EdgePoint> >
+      PNMap;
+    PNMap posnode;
+
+    static double angle(double x, double y) {
+      return atan2(y,x);
+    }
+
+    bool checkIntersections(const EdgeShape &es, Arc b) {
+      Arc a;
+      firstCwF(a, leftFace(b));
+      while (a != INVALID) {
+        if (EdgeType::intersect(es, edgeshapes[a])) return true;
+        nextCwF(a);
+      }
+      return false;
+    }
+
+  protected:
+
+    //Intersection between two faces. O(n*m)
+    bool polygonIntersect(Face f1, Face f2) {
+      Arc e1, e2;
+      for(firstCwF(e1,f1); e1 != INVALID; nextCwF(e1)) {
+        for(firstCwF(e2,f2); e2 != INVALID; nextCwF(e2)) {
+          if (EdgeType::intersect(edgeshapes[e1], edgeshapes[e2])) return true;
+        }
+      }
+      return false;
+    }
+
+  public:
+
+    ///\brief Constructor.
+
+    ///Constructor.
+    PlaneGraph() : nodepos(*this), edgeshapes(*this) {}
+
+    ///\brief Constructor that copies a planar drawing
+
+    ///Constructor that copies a planar drawing.
+    template<typename Graph>
+    PlaneGraph(const Graph &graph, const lemon::PlanarDrawing<Graph> &drawing) :
+      nodepos(*this), edgeshapes(*this) {
+        typename Graph::template NodeMap<Node> nm(graph);
+        for (typename Graph::NodeIt it(graph); it != INVALID; ++it) {
+          nm[it] = addNode(drawing[it]);
+        }
+        for (typename Graph::EdgeIt it(graph); it != INVALID; ++it) {
+          addEdge(nm[graph.u(it)], nm[graph.v(it)]);
+        }
+    }
+
+    template<typename Graph>
+    PlaneGraph(const Graph &graph, const lemon::PlanarDrawing<Graph> &drawing,
+        typename Graph::template NodeMap<Node> &nm,
+        typename Graph::template ArcMap<Arc> &em) :
+      nodepos(*this), edgeshapes(*this) {
+        for (typename Graph::NodeIt it(graph); it != INVALID; ++it) {
+          nm[it] = addNode(drawing[it]);
+        }
+        for (typename Graph::EdgeIt it(graph); it != INVALID; ++it) {
+          em[graph.direct(it,graph.u(it))] = direct(addEdge(nm[graph.u(it)
+            ], nm[graph.v(it)]),nm[graph.u(it)]);
+          em[graph.direct(it,graph.v(it))] = oppositeArc(em[graph.direct(
+            it,graph.u(it))]);
+        }
+    }
+
+    ///Add a new node to the graph at the given point.
+    ///\return The new node, or INVALID if the given point already contains a
+    ///node.
+    Node addNode(const EdgePoint &point) {
+      typename PNMap::iterator it = posnode.find(point);
+      if (it != posnode.end()) return INVALID;
+      Node n = PlanarGraph::addNode();
+      nodepos[n] = point;
+      posnode[point] = n;
+      return n;
+    }
+
+    ///Add a new node to the graph at the given point.
+    ///\return The new node, or INVALID if the given point already contains a
+    ///node.
+    Node addNode(double x, double y) {
+      return addNode(EdgePoint(x,y));
+    }
+
+    ///Add a new edge to the graph.
+
+    ///Add a new edge to the graph between the nodes \c n1 and \c n2.
+    ///This function should only be used if the edge type doesn't
+    ///need parameters other than the locations of the two nodes.
+    ///\return The new edge, or INVALID if the topology doesn't allow the
+    ///addition.
+    Edge addEdge(Node n1, Node n2) {
+      EdgeShape es = EdgeType::create(nodepos[n1],nodepos[n2]);
+      return addEdge(n1, n2, es);
+    }
+
+    ///Add a new edge to the graph.
+
+    ///Add a new edge to the graph with the given edge shape \c es.
+    ///\return The new edge, or INVALID if \c es doesn't connect the locations
+    ///of two existing nodes or the topology doesn't allow the
+    ///addition.
+    Edge addEdge(const EdgeShape &es) {
+      typename PNMap::iterator it = posnode.find(EdgeType::source(es));
+      if (it == posnode.end()) return INVALID;
+      Node n1 = it->second;
+      it = posnode.find(EdgeType::target(es));
+      if (it == posnode.end()) return INVALID;
+      Node n2 = it->second;
+      return addEdge(n1, n2, es);
+    }
+
+    ///Add a new edge to the graph.
+
+    ///Add a new edge to the graph between the nodes \c n1 and \c n2, using the
+    ///edge shape \c es. The coordinates of \c n1, \c n2 and \c es will be used
+    ///to determine the correct topology.
+    ///\note For performance reasons, this function doesn't check if the source
+    ///and target coordinates of \c es match those of \c n1 and \c n2.
+    ///\return The new edge, or INVALID if the topology doesn't allow the
+    ///addition.
+    Edge addEdge(Node n1, Node n2, const EdgeShape &es) {
+      if (EdgeType::degenerate(es)) return INVALID;
+      EdgeShape esr = EdgeType::reverse(es);
+      Arc p_u;
+      firstOut(p_u,n1);
+      if (p_u != INVALID) {
+        Arc p_u0 = p_u;
+        nextOut(p_u0);
+        while (p_u0 != INVALID && !EdgeType::between(es, edgeshapes[p_u],
+            edgeshapes[p_u0])) {
+          p_u = p_u0;
+          nextOut(p_u0);
+        }
+      }
+      Arc p_v;
+      firstOut(p_v,n2);
+      if (p_v != INVALID) {
+        Arc p_v0 = p_v;
+        nextOut(p_v0);
+        while (p_v0 != INVALID && !EdgeType::between(esr, edgeshapes[p_v],
+            edgeshapes[p_v0])) {
+          p_v = p_v0;
+          nextOut(p_v0);
+        }
+      }
+      if (p_u != INVALID && p_v != INVALID) {
+        if (idComponent(n1) == idComponent(n2)) {
+          if (leftFace(p_u) != leftFace(p_v))
+            return INVALID;
+        } else {
+          if (polygonIntersect(leftFace(p_u),leftFace(p_v)))
+            return INVALID;
+        }
+      }
+
+      if (p_u != INVALID && checkIntersections(es,p_u)) return INVALID;
+      if (p_v != INVALID && (p_u == INVALID || leftFace(p_u) != leftFace(p_v))
+        && checkIntersections(es,p_v)) return INVALID;
+
+#ifdef REMOVE_BEFORE_RELEASE
+            std::cout << "AddArc: " << id(n1) << "->" << id(n2) << ", p_u: "
+             << id(p_u) << ", p_v: " << id(p_v) << std::endl;
+#endif
+      Edge e = PlanarGraph::addEdge(n1,n2,p_u,p_v);
+      if (e != INVALID) {
+        if (n1 != n2 || EdgeType::between(esr,edgeshapes[p_u],es)) {
+          edgeshapes[direct(e,true)] = es;
+          edgeshapes[direct(e,false)] = esr;
+        } else {
+          edgeshapes[direct(e,true)] = esr;
+          edgeshapes[direct(e,false)] = es;
+        }
+      }
+      return e;
+
+    }
+
+    void erase(Node n) {
+      EdgePoint p = nodepos[n];
+      posnode.remove(p);
+      PlanarGraph::erase(n);
+    }
+
+    void erase(Edge e) {
+      PlanarGraph::erase(e);
+    }
+
+    void clear() {
+      posnode.clear();
+      PlanarGraph::clear();
+    }
+
+    ///Get a map of the node locations.
+    ///Get a map of the node locations.
+    const NodeMap<EdgePoint> &nodePosMap() const {
+      return nodepos;
+    }
+
+    ///Tell the location of a node.
+    ///Tell the location of a node.
+    const EdgePoint &locate(const Node &n) const {
+      return nodepos[n];
+    }
+
+    ///Tell the location of an arc.
+    ///Tell the location of an arc.
+    const EdgeShape &locate(const Arc &a) const {
+      return edgeshapes[a];
+    }
+
+
+    void edgeToEps(std::ostream &os, const Arc &e, const Color &col,
+      double width) const {
+      EdgeType::toEps(os,edgeshapes[e],col,width);
+    }
+
+
+    /// Iterator class for the edge shapes around a node.
+
+    /// This iterator goes through the arcs around a node anticlockwise.
+    class CcwIncIt : public Arc {
+      const PlaneGraph* _graph;
+      const Node _node;
+    public:
+
+      CcwIncIt() { }
+
+      CcwIncIt(Invalid i) : Arc(i) { }
+
+      CcwIncIt(const PlaneGraph& graph, const Node& node)
+          : _graph(&graph), _node(node) {
+        _graph->firstOut(*this, node);
+      }
+
+      CcwIncIt(const Graph& graph, const Arc& arc)
+          : Arc(arc), _graph(&graph) {}
+
+      CcwIncIt& operator++() {
+        _graph->nextOut(*this, _node);
+        return *this;
+      }
+
+      const EdgeShape &operator*() {
+        return _graph->edgeshapes[*this];
+      }
+
+    };
+
+    /// Iterator class for the edge shapes around a node.
+
+    /// This iterator goes through the arcs around a node clockwise.
+    class CwIncIt : public Arc {
+      const PlaneGraph* _graph;
+      const Node _node;
+    public:
+
+      CwIncIt() { }
+
+      CwIncIt(Invalid i) : Arc(i) { }
+
+      CwIncIt(const PlaneGraph& graph, const Node& node)
+          : _graph(&graph), _node(node) {
+        _graph->lastOut(*this, node);
+      }
+
+      CwIncIt(const Graph& graph, const Arc& arc)
+          : Arc(arc), _graph(&graph) {}
+
+      CwIncIt& operator++() {
+        _graph->prevOut(*this, _node);
+        return *this;
+      }
+
+      const EdgeShape &operator*() {
+        return _graph->edgeshapes[*this];
+      }
+
+    };
+
+    /// Iterator class for the edge shapes on the boundary of a face.
+
+    /// This iterator goes through the arcs on the boundary of a face clockwise.
+    class CwBoundaryIt {
+      const PlaneGraph* _graph;
+      Face _face;
+      Arc f_arc;
+    public:
+
+      CwBoundaryIt() { }
+
+      CwBoundaryIt(const PlaneGraph& graph, const Face& face)
+          : _graph(&graph), _face(face) {
+        _graph->firstCwF(f_arc,face);
+      }
+
+      CwBoundaryIt(const Graph& graph, const Arc& arc)
+          : f_arc(arc), _graph(&graph) {}
+
+      CwBoundaryIt& operator++() {
+        _graph->nextCwF(f_arc);
+        return *this;
+      }
+
+      const EdgeShape &operator*() {
+        return _graph->edgeshapes[f_arc];
+      }
+
+    };
+
+    /// Iterator class for the edge shapes on the boundary of a face.
+
+    /// This iterator goes through the arcs on the boundary of a face
+    /// anticlockwise.
+    class CcwBoundaryIt {
+      const PlaneGraph* _graph;
+      Face _face;
+      Arc f_arc;
+    public:
+
+      CcwBoundaryIt() { }
+
+      CcwBoundaryIt(const PlaneGraph& graph, const Face& face)
+          : _graph(&graph), _face(face) {
+        _graph->firstCcwF(f_arc,face);
+      }
+
+      CcwBoundaryIt(const Graph& graph, const Arc& arc)
+          : f_arc(arc), _graph(&graph) {}
+
+      CcwBoundaryIt& operator++() {
+        _graph->nextCcwF(f_arc);
+        return *this;
+      }
+
+      const EdgeShape &operator*() {
+        return _graph->edgeshapes[f_arc];
+      }
+
+    };
+
+  };
+  /// @}
+
+} //namespace lemon
+
+#endif

test/CMakeLists.txt

@@ -35 +35 @@
   min_cost_flow_test
   min_mean_cycle_test
   path_test
+  planar_graph_test
   planarity_test
   preflow_test
   radix_sort_test

test/Makefile.am

@@ -37 +37 @@
         test/min_cost_flow_test \
         test/min_mean_cycle_test \
         test/path_test \
+        test/planar_graph_test \
         test/planarity_test \
         test/preflow_test \
         test/radix_sort_test \
@@ -88 +89 @@
 test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
 test_min_mean_cycle_test_SOURCES = test/min_mean_cycle_test.cc
 test_path_test_SOURCES = test/path_test.cc
+test_planar_graph_test_SOURCES = test/planar_graph_test.cc
 test_planarity_test_SOURCES = test/planarity_test.cc
 test_preflow_test_SOURCES = test/preflow_test.cc
 test_radix_sort_test_SOURCES = test/radix_sort_test.cc

new file test/planar_graph_test.cc

@@ +1 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2010
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#include <lemon/concepts/graph.h>
+#include <lemon/planar_graph.h>
+#include <lemon/plane_graph.h>
+
+#include "test_tools.h"
+#include "planar_graph_test.h"
+
+using namespace lemon;
+using namespace lemon::concepts;
+
+template <class Graph>
+void checkGraphBuild() {
+  TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+  Graph G;
+  checkGraphNodeList(G, 0);
+  checkGraphEdgeList(G, 0);
+  checkGraphArcList(G, 0);
+  checkGraphFaceList(G, 0);
+
+  G.reserveNode(3);
+  G.reserveEdge(3);
+
+  Node
+    n1 = G.addNode(),
+    n2 = G.addNode(),
+    n3 = G.addNode();
+  checkGraphNodeList(G, 3);
+  checkGraphEdgeList(G, 0);
+  checkGraphArcList(G, 0);
+  checkGraphFaceList(G, 3);
+
+  Edge e1 = G.addEdge(n1, n2, INVALID, INVALID);
+  check((G.u(e1) == n1 && G.v(e1) == n2) || (G.u(e1) == n2 && G.v(e1) == n1),
+        "Wrong edge");
+
+  checkGraphNodeList(G, 3);
+  checkGraphEdgeList(G, 1);
+  checkGraphArcList(G, 2);
+  checkGraphFaceList(G, 2);
+
+  checkGraphIncEdgeArcLists(G, n1, 1);
+  checkGraphIncEdgeArcLists(G, n2, 1);
+  checkGraphIncEdgeArcLists(G, n3, 0);
+
+  checkGraphConEdgeList(G, 1);
+  checkGraphConArcList(G, 2);
+
+  Edge e2 = G.addEdge(n2, n1, e1, e1),
+       e3 = G.addEdge(n2, n3, e1, INVALID);
+
+  checkGraphNodeList(G, 3);
+  checkGraphEdgeList(G, 3);
+  checkGraphArcList(G, 6);
+  checkGraphFaceList(G, 2);
+
+  checkGraphIncEdgeArcLists(G, n1, 2);
+  checkGraphIncEdgeArcLists(G, n2, 3);
+  checkGraphIncEdgeArcLists(G, n3, 1);
+
+  checkGraphConEdgeList(G, 3);
+  checkGraphConArcList(G, 6);
+
+  checkArcDirections(G);
+
+  checkNodeIds(G);
+  checkArcIds(G);
+  checkEdgeIds(G);
+  checkFaceIds(G);
+  checkGraphNodeMap(G);
+  checkGraphArcMap(G);
+  checkGraphFaceMap(G);
+}
+
+template <class Graph>
+void checkPlaneGraphBuild() {
+  TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+  Graph G;
+  checkGraphNodeList(G, 0);
+  checkGraphEdgeList(G, 0);
+  checkGraphArcList(G, 0);
+  checkGraphFaceList(G, 0);
+
+  G.reserveNode(3);
+  G.reserveEdge(3);
+
+  Node
+    n1 = G.addNode(dim2::Point<double>(0.0,0.0)),
+    n2 = G.addNode(dim2::Point<double>(1.0,0.0)),
+    n3 = G.addNode(dim2::Point<double>(0.5,1.0));
+  checkGraphNodeList(G, 3);
+  checkGraphEdgeList(G, 0);
+  checkGraphArcList(G, 0);
+  checkGraphFaceList(G, 3);
+
+  Edge e1 = G.addEdge(n1, n2);
+  check((G.u(e1) == n1 && G.v(e1) == n2) || (G.u(e1) == n2 && G.v(e1) == n1),
+        "Wrong edge");
+
+  checkGraphNodeList(G, 3);
+  checkGraphEdgeList(G, 1);
+  checkGraphArcList(G, 2);
+  checkGraphFaceList(G, 2);
+
+  checkGraphIncEdgeArcLists(G, n1, 1);
+  checkGraphIncEdgeArcLists(G, n2, 1);
+  checkGraphIncEdgeArcLists(G, n3, 0);
+
+  checkGraphConEdgeList(G, 1);
+  checkGraphConArcList(G, 2);
+
+  Edge e2 = G.addEdge(n2, n3),
+       e3 = G.addEdge(n3, n1);
+
+  checkGraphNodeList(G, 3);
+  checkGraphEdgeList(G, 3);
+  checkGraphArcList(G, 6);
+  checkGraphFaceList(G, 2);
+
+  checkGraphIncEdgeArcLists(G, n1, 2);
+  checkGraphIncEdgeArcLists(G, n2, 2);
+  checkGraphIncEdgeArcLists(G, n3, 2);
+
+  checkGraphConEdgeList(G, 3);
+  checkGraphConArcList(G, 6);
+
+  checkArcDirections(G);
+
+  checkNodeIds(G);
+  checkArcIds(G);
+  checkEdgeIds(G);
+  checkFaceIds(G);
+  checkGraphNodeMap(G);
+  checkGraphArcMap(G);
+  checkGraphFaceMap(G);
+}
+
+template <class Graph>
+void checkGraphArcSplit() {
+  TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+  Graph G;
+  Node n0 = G.addNode(),
+       n1 = G.addNode(),
+       n2 = G.addNode(),
+       n3 = G.addNode();
+  Edge a0 = G.addEdge(n0,n1,INVALID,INVALID),
+       a1 = G.addEdge(n0,n3,a0,INVALID),
+       a2 = G.addEdge(n0,n2,a0,INVALID),
+       a3 = G.addEdge(n1,n3,a0,a1),
+       a4 = G.addEdge(n3,n2,a1,a2);
+
+  checkGraphNodeList(G, 4);
+  checkGraphEdgeList(G, 5);
+  checkGraphArcList(G, 10);
+  checkGraphFaceList(G, 3);
+
+  G.split(a1);
+
+  checkGraphNodeList(G, 5);
+  checkGraphEdgeList(G, 6);
+  checkGraphArcList(G, 12);
+  checkGraphFaceList(G, 3);
+
+  checkGraphIncEdgeArcLists(G, n0, 3);
+  checkGraphIncEdgeArcLists(G, n1, 2);
+  checkGraphIncEdgeArcLists(G, n2, 2);
+  checkGraphIncEdgeArcLists(G, n3, 3);
+
+  checkGraphBoundaryArcList(G, G.leftFace(G.direct(a0,true)), 4);
+  checkGraphBoundaryArcList(G, G.leftFace(G.direct(a0,false)), 4);
+  checkGraphBoundaryArcList(G, G.leftFace(G.direct(a1,false)), 4);
+
+  checkGraphConEdgeList(G, 6);
+  checkGraphConArcList(G, 12);
+
+}
+
+template <class Graph>
+void checkGraphNodeSplit() {
+  TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+    Graph G;
+    Node n0 = G.addNode(),
+         n1 = G.addNode(),
+         n2 = G.addNode(),
+         n3 = G.addNode(),
+         n4 = G.addNode();
+    Edge a0 = G.addEdge(n4,n3,INVALID,INVALID),
+         a1 = G.addEdge(n4,n0,a0,INVALID),
+         a2 = G.addEdge(n1,n4,INVALID,a1),
+         a3 = G.addEdge(n4,n2,a2,INVALID),
+         a4 = G.addEdge(n2,n3,a3,a0);
+
+    checkGraphNodeList(G, 5);
+    checkGraphEdgeList(G, 5);
+    checkGraphArcList(G, 10);
+    checkGraphFaceList(G, 2);
+
+    checkGraphIncEdgeArcLists(G, n0, 1);
+    checkGraphIncEdgeArcLists(G, n1, 1);
+    checkGraphIncEdgeArcLists(G, n2, 2);
+    checkGraphIncEdgeArcLists(G, n3, 2);
+    checkGraphIncEdgeArcLists(G, n4, 4);
+
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a0,false)), 3);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a0,true)), 7);
+
+    Node n5 = G.split(n4,a3,a1,true);
+
+    checkGraphNodeList(G, 6);
+    checkGraphEdgeList(G, 6);
+    checkGraphArcList(G, 12);
+    checkGraphFaceList(G, 2);
+
+    checkGraphIncEdgeArcLists(G, n0, 1);
+    checkGraphIncEdgeArcLists(G, n1, 1);
+    checkGraphIncEdgeArcLists(G, n2, 2);
+    checkGraphIncEdgeArcLists(G, n3, 2);
+    checkGraphIncEdgeArcLists(G, n4, 3);
+    checkGraphIncEdgeArcLists(G, n5, 3);
+
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a0,false)), 3);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a0,true)), 9);
+}
+
+template <class Graph>
+void checkGraphContract() {
+  TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+    Graph G;
+    Node n0 = G.addNode(),
+         n1 = G.addNode(),
+         n2 = G.addNode(),
+         n3 = G.addNode();
+    Edge a0 = G.addEdge(n0,n1,INVALID,INVALID),
+         a1 = G.addEdge(n0,n1,a0,a0),
+         a2 = G.addEdge(n0,n1,a1,a0),
+         a3 = G.addEdge(n0,n2,a1,INVALID),
+         a4 = G.addEdge(n0,n3,a3,INVALID),
+         a5 = G.addEdge(n2,n1,a3,a2),
+         a6 = G.addEdge(n3,n1,a4,a2);
+
+    checkGraphNodeList(G, 4);
+    checkGraphEdgeList(G, 7);
+    checkGraphArcList(G, 14);
+    checkGraphFaceList(G, 5);
+
+    checkGraphIncEdgeArcLists(G, n0, 5);
+    checkGraphIncEdgeArcLists(G, n1, 5);
+    checkGraphIncEdgeArcLists(G, n2, 2);
+    checkGraphIncEdgeArcLists(G, n3, 2);
+
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a0,false)), 2);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a0,true)), 2);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a1,true)), 3);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a2,false)), 3);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a3,true)), 4);
+
+    G.contract(a0);
+
+    checkGraphNodeList(G, 3);
+    checkGraphEdgeList(G, 6);
+    checkGraphArcList(G, 12);
+    checkGraphFaceList(G, 5);
+
+    checkGraphIncEdgeArcLists(G, n0, 8);
+    checkGraphIncEdgeArcLists(G, n2, 2);
+    checkGraphIncEdgeArcLists(G, n3, 2);
+
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a1,false)), 1);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a2,true)), 1);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a1,true)), 3);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a2,false)), 3);
+    checkGraphBoundaryArcList(G, G.leftFace(G.direct(a3,true)), 4);
+
+    checkGraphCommonFaceList(G, n0, n2, 2);
+    checkGraphCommonFaceList(G, n0, n3, 2);
+    checkGraphCommonFaceList(G, n2, n3, 1);
+
+    checkGraphCommonNodeList(G, G.leftFace(G.direct(a3,true)), G.leftFace(G.
+            direct(a3,false)), 2);
+    checkGraphCommonNodeList(G, G.leftFace(G.direct(a1,false)), G.leftFace(G.
+            direct(a2,true)), 1);
+
+}
+
+template <class Graph>
+void checkGraphErase() {
+  TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+  Graph G;
+  Node n1 = G.addNode(),
+       n2 = G.addNode(),
+       n3 = G.addNode(),
+       n4 = G.addNode();
+  Edge e1 = G.addEdge(n1, n2, INVALID, INVALID),
+       e2 = G.addEdge(n2, n1, e1, e1),
+       e3 = G.addEdge(n2, n3, e1, INVALID),
+       e4 = G.addEdge(n1, n4, e2, INVALID),
+       e5 = G.addEdge(n4, n3, e4, e3);
+
+  // Check edge deletion
+  G.erase(e2);
+
+  checkGraphNodeList(G, 4);
+  checkGraphEdgeList(G, 4);
+  checkGraphArcList(G, 8);
+  checkGraphFaceList(G, 2);
+
+  checkGraphIncEdgeArcLists(G, n1, 2);
+  checkGraphIncEdgeArcLists(G, n2, 2);
+  checkGraphIncEdgeArcLists(G, n3, 2);
+  checkGraphIncEdgeArcLists(G, n4, 2);
+
+  checkGraphBoundaryArcList(G, G.leftFace(G.direct(e1,true)), 4);
+  checkGraphBoundaryArcList(G, G.leftFace(G.direct(e1,false)), 4);
+
+  checkGraphConEdgeList(G, 4);
+  checkGraphConArcList(G, 8);
+
+  // Check node deletion
+  G.erase(n3);
+
+  checkGraphNodeList(G, 3);
+  checkGraphEdgeList(G, 2);
+  checkGraphArcList(G, 4);
+  checkGraphFaceList(G, 1);
+
+  checkGraphIncEdgeArcLists(G, n1, 2);
+  checkGraphIncEdgeArcLists(G, n2, 1);
+  checkGraphIncEdgeArcLists(G, n4, 1);
+
+  checkGraphBoundaryArcList(G, G.leftFace(G.direct(e1,true)), 4);
+
+  checkGraphConEdgeList(G, 2);
+  checkGraphConArcList(G, 4);
+}
+
+
+template <class Graph>
+void checkGraphSnapshot() {
+  TEMPLATE_GRAPH_TYPEDEFS(Graph);
+
+  Graph G;
+  Node n1 = G.addNode(),
+       n2 = G.addNode(),
+       n3 = G.addNode();
+  Edge e1 = G.addEdge(n1, n2, INVALID, INVALID),
+       e2 = G.addEdge(n2, n1, e1, e1),
+       e3 = G.addEdge(n2, n3, e1, INVALID);
+
+  checkGraphNodeList(G, 3);
+  checkGraphEdgeList(G, 3);
+  checkGraphArcList(G, 6);
+  checkGraphFaceList(G, 2);
+
+  typename Graph::Snapshot snapshot(G);
+
+  Node n = G.addNode();
+  Edge ea4 = G.addEdge(n3, n, e3, INVALID),
+       ea5 = G.addEdge(n, n3, ea4, ea4),
+       ea6 = G.addEdge(n3, n2, ea5, e3);
+
+  checkGraphNodeList(G, 4);
+  checkGraphEdgeList(G, 6);
+  checkGraphArcList(G, 12);
+  checkGraphFaceList(G, 4);
+
+  snapshot.restore();
+
+  checkGraphNodeList(G, 3);
+  checkGraphEdgeList(G, 3);
+  checkGraphArcList(G, 6);
+  checkGraphFaceList(G, 2);
+
+  checkGraphIncEdgeArcLists(G, n1, 2);
+  checkGraphIncEdgeArcLists(G, n2, 3);
+  checkGraphIncEdgeArcLists(G, n3, 1);
+  checkGraphBoundaryArcList(G, G.leftFace(G.direct(e1,true)), 2);
+  checkGraphBoundaryArcList(G, G.leftFace(G.direct(e1,false)), 4);
+
+  checkGraphConEdgeList(G, 3);
+  checkGraphConArcList(G, 6);
+
+  checkNodeIds(G);
+  checkEdgeIds(G);
+  checkArcIds(G);
+  checkFaceIds(G);
+  checkGraphNodeMap(G);
+  checkGraphEdgeMap(G);
+  checkGraphArcMap(G);
+  checkGraphFaceMap(G);
+
+  G.addNode();
+  snapshot.save(G);
+
+  G.addEdge(G.addNode(), G.addNode(), INVALID, INVALID);
+
+  snapshot.restore();
+  snapshot.save(G);
+
+  checkGraphNodeList(G, 4);
+  checkGraphEdgeList(G, 3);
+  checkGraphArcList(G, 6);
+  checkGraphFaceList(G, 3);
+
+  G.addEdge(G.addNode(), G.addNode(), INVALID, INVALID);
+
+  snapshot.restore();
+
+  checkGraphNodeList(G, 4);
+  checkGraphEdgeList(G, 3);
+  checkGraphArcList(G, 6);
+  checkGraphFaceList(G, 3);
+}
+
+template <typename Graph>
+void checkGraphValidity() {
+  TEMPLATE_GRAPH_TYPEDEFS(Graph);
+  Graph g;
+
+  Node
+    n1 = g.addNode(),
+    n2 = g.addNode(),
+    n3 = g.addNode();
+
+  Edge
+    e1 = g.addEdge(n1, n2, INVALID, INVALID),
+    e2 = g.addEdge(n2, n3, e1, INVALID);
+
+  check(g.valid(n1), "Wrong validity check");
+  check(g.valid(e1), "Wrong validity check");
+  check(g.valid(g.direct(e1, true)), "Wrong validity check");
+
+  check(!g.valid(g.nodeFromId(-1)), "Wrong validity check");
+  check(!g.valid(g.edgeFromId(-1)), "Wrong validity check");
+  check(!g.valid(g.arcFromId(-1)), "Wrong validity check");
+  check(!g.valid(g.faceFromId(-1)), "Wrong validity check");
+}
+
+template <typename Graph>
+void checkGraphValidityErase() {
+  TEMPLATE_GRAPH_TYPEDEFS(Graph);
+  Graph g;
+
+  Node
+    n1 = g.addNode(),
+    n2 = g.addNode(),
+    n3 = g.addNode();
+
+  Edge
+    e1 = g.addEdge(n1, n2, INVALID, INVALID),
+    e2 = g.addEdge(n2, n3, e1, INVALID);
+
+  check(g.valid(n1), "Wrong validity check");
+  check(g.valid(e1), "Wrong validity check");
+  check(g.valid(g.direct(e1, true)), "Wrong validity check");
+
+  g.erase(n1);
+
+  check(!g.valid(n1), "Wrong validity check");
+  check(g.valid(n2), "Wrong validity check");
+  check(g.valid(n3), "Wrong validity check");
+  check(!g.valid(e1), "Wrong validity check");
+  check(g.valid(e2), "Wrong validity check");
+
+  check(!g.valid(g.nodeFromId(-1)), "Wrong validity check");
+  check(!g.valid(g.edgeFromId(-1)), "Wrong validity check");
+  check(!g.valid(g.arcFromId(-1)), "Wrong validity check");
+  check(!g.valid(g.faceFromId(-1)), "Wrong validity check");
+}
+
+void checkGraphs() {
+  { // Checking PlanarGraph
+    checkGraphBuild<PlanarGraph>();
+    checkPlaneGraphBuild<PlaneGraph<> >();
+    checkGraphArcSplit<PlanarGraph>();
+    checkGraphNodeSplit<PlanarGraph>();
+    checkGraphContract<PlanarGraph>();
+    checkGraphErase<PlanarGraph>();
+    checkGraphSnapshot<PlanarGraph>();
+    checkGraphValidityErase<PlanarGraph>();
+  }
+}
+
+int main() {
+  checkGraphs();
+  return 0;
+}

new file test/planar_graph_test.h

@@ +1 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2003-2009
+ * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
+ * (Egervary Research Group on Combinatorial Optimization, EGRES).
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_TEST_GRAPH_TEST_H
+#define LEMON_TEST_GRAPH_TEST_H
+
+#include <set>
+
+#include <lemon/core.h>
+#include <lemon/maps.h>
+
+#include "test_tools.h"
+
+namespace lemon {
+
+  template<class Graph>
+  void checkGraphNodeList(const Graph &G, int cnt)
+  {
+    typename Graph::NodeIt n(G);
+    for(int i=0;i<cnt;i++) {
+      check(n!=INVALID,"Wrong Node list linking.");
+      ++n;
+    }
+    check(n==INVALID,"Wrong Node list linking.");
+    check(countNodes(G)==cnt,"Wrong Node number.");
+  }
+
+  template<class Graph>
+  void checkGraphFaceList(const Graph &G, int cnt)
+  {
+    typename Graph::FaceIt n(G);
+    for(int i=0;i<cnt;i++) {
+      check(n!=INVALID,"Wrong Face list linking.");
+      ++n;
+    }
+    check(n==INVALID,"Wrong Face list linking.");
+    check(countFaces(G)==cnt,"Wrong Face number.");
+  }
+
+  template<class Graph>
+  void checkGraphArcList(const Graph &G, int cnt)
+  {
+    typename Graph::ArcIt e(G);
+    for(int i=0;i<cnt;i++) {
+      check(e!=INVALID,"Wrong Arc list linking.");
+      check(G.oppositeNode(G.source(e), e) == G.target(e),
+            "Wrong opposite node");
+      check(G.oppositeNode(G.target(e), e) == G.source(e),
+            "Wrong opposite node");
+      ++e;
+    }
+    check(e==INVALID,"Wrong Arc list linking.");
+    check(countArcs(G)==cnt,"Wrong Arc number.");
+  }
+
+  template<class Graph>
+  void checkGraphOutArcList(const Graph &G, typename Graph::Node n, int cnt)
+  {
+    typename Graph::OutArcIt e(G,n);
+    for(int i=0;i<cnt;i++) {
+      check(e!=INVALID,"Wrong OutArc list linking.");
+      check(n==G.source(e),"Wrong OutArc list linking.");
+      check(n==G.baseNode(e),"Wrong OutArc list linking.");
+      check(G.target(e)==G.runningNode(e),"Wrong OutArc list linking.");
+      ++e;
+    }
+    check(e==INVALID,"Wrong OutArc list linking.");
+    check(countOutArcs(G,n)==cnt,"Wrong OutArc number.");
+  }
+
+  template<class Graph>
+  void checkGraphInArcList(const Graph &G, typename Graph::Node n, int cnt)
+  {
+    typename Graph::InArcIt e(G,n);
+    for(int i=0;i<cnt;i++) {
+      check(e!=INVALID,"Wrong InArc list linking.");
+      check(n==G.target(e),"Wrong InArc list linking.");
+      check(n==G.baseNode(e),"Wrong OutArc list linking.");
+      check(G.source(e)==G.runningNode(e),"Wrong OutArc list linking.");
+      ++e;
+    }
+    check(e==INVALID,"Wrong InArc list linking.");
+    check(countInArcs(G,n)==cnt,"Wrong InArc number.");
+  }
+
+  template<class Graph>
+  void checkGraphEdgeList(const Graph &G, int cnt)
+  {
+    typename Graph::EdgeIt e(G);
+    for(int i=0;i<cnt;i++) {
+      check(e!=INVALID,"Wrong Edge list linking.");
+      check(G.oppositeNode(G.u(e), e) == G.v(e), "Wrong opposite node");
+      check(G.oppositeNode(G.v(e), e) == G.u(e), "Wrong opposite node");
+      ++e;
+    }
+    check(e==INVALID,"Wrong Edge list linking.");
+    check(countEdges(G)==cnt,"Wrong Edge number.");
+  }
+
+  template<class Graph>
+  void checkGraphIncEdgeList(const Graph &G, typename Graph::Node n, int cnt)
+  {
+    typename Graph::IncEdgeIt e(G,n);
+    for(int i=0;i<cnt;i++) {
+      check(e!=INVALID,"Wrong IncEdge list linking.");
+      check(n==G.u(e) || n==G.v(e),"Wrong IncEdge list linking.");
+      check(n==G.baseNode(e),"Wrong OutArc list linking.");
+      check(G.u(e)==G.runningNode(e) || G.v(e)==G.runningNode(e),
+            "Wrong OutArc list linking.");
+      ++e;
+    }
+    check(e==INVALID,"Wrong IncEdge list linking.");
+    check(countIncEdges(G,n)==cnt,"Wrong IncEdge number.");
+  }
+
+  template <class Graph>
+  void checkGraphIncEdgeArcLists(const Graph &G, typename Graph::Node n,
+                                 int cnt)
+  {
+    checkGraphIncEdgeList(G, n, cnt);
+    checkGraphOutArcList(G, n, cnt);
+    checkGraphInArcList(G, n, cnt);
+  }
+
+  template<class Graph>
+  void checkGraphBoundaryArcList(const Graph &G, typename Graph::Face n, int
+    cnt)
+  {
+    typename Graph::CwBoundaryArcIt e(G,n);
+    for(int i=0;i<cnt;i++) {
+      check(e!=INVALID,"Wrong CwBoundaryArc list linking.");
+      check(n==G.w1(e) || n==G.w2(e),"Wrong CwBoundaryArc list linking.");
+      ++e;
+    }
+    check(e==INVALID,"Wrong BoundaryArc list linking.");
+    check(countBoundaryArcs(G,n)==cnt,"Wrong IncEdge number.");
+  }
+
+  template <class Graph>
+  void checkGraphConArcList(const Graph &G, int cnt) {
+    int i = 0;
+    for (typename Graph::NodeIt u(G); u != INVALID; ++u) {
+      for (typename Graph::NodeIt v(G); v != INVALID; ++v) {
+        for (ConArcIt<Graph> a(G, u, v); a != INVALID; ++a) {
+          check(G.source(a) == u, "Wrong iterator.");
+          check(G.target(a) == v, "Wrong iterator.");
+          ++i;
+        }
+      }
+    }
+    check(cnt == i, "Wrong iterator.");
+  }
+
+  template <class Graph>
+  void checkGraphConEdgeList(const Graph &G, int cnt) {
+    int i = 0;
+    for (typename Graph::NodeIt u(G); u != INVALID; ++u) {
+      for (typename Graph::NodeIt v(G); v != INVALID; ++v) {
+        for (ConEdgeIt<Graph> e(G, u, v); e != INVALID; ++e) {
+          check((G.u(e) == u && G.v(e) == v) ||
+                (G.u(e) == v && G.v(e) == u), "Wrong iterator.");
+          i += u == v ? 2 : 1;
+        }
+      }
+    }
+    check(2 * cnt == i, "Wrong iterator.");
+  }
+
+  template<class Graph>
+  void checkGraphCommonNodeList(const Graph &G, typename Graph::Face f1,
+      typename Graph::Face f2, int cnt)
+  {
+    typename Graph::CommonNodesIt e(G,f1,f2);
+    for(int i=0;i<cnt;i++) {
+      check(e!=INVALID,"Wrong CommonNodeIt.");
+      ++e;
+    }
+    check(e==INVALID,"Wrong CommonNodeIt.");
+  }
+
+  template<class Graph>
+  void checkGraphCommonFaceList(const Graph &G, typename Graph::Node n1,
+      typename Graph::Node n2, int cnt)
+  {
+    typename Graph::CommonFacesIt e(G,n1,n2);
+    for(int i=0;i<cnt;i++) {
+      check(e!=INVALID,"Wrong CommonNodeIt.");
+      ++e;
+    }
+    check(e==INVALID,"Wrong CommonNodeIt.");
+  }
+
+  template <typename Graph>
+  void checkArcDirections(const Graph& G) {
+    for (typename Graph::ArcIt a(G); a != INVALID; ++a) {
+      check(G.source(a) == G.target(G.oppositeArc(a)), "Wrong direction");
+      check(G.target(a) == G.source(G.oppositeArc(a)), "Wrong direction");
+      check(G.direct(a, G.direction(a)) == a, "Wrong direction");
+    }
+  }
+
+  template <typename Graph>
+  void checkNodeIds(const Graph& G) {
+    std::set<int> values;
+    for (typename Graph::NodeIt n(G); n != INVALID; ++n) {
+      check(G.nodeFromId(G.id(n)) == n, "Wrong id");
+      check(values.find(G.id(n)) == values.end(), "Wrong id");
+      check(G.id(n) <= G.maxNodeId(), "Wrong maximum id");
+      values.insert(G.id(n));
+    }
+  }
+
+  template <typename Graph>
+  void checkArcIds(const Graph& G) {
+    std::set<int> values;
+    for (typename Graph::ArcIt a(G); a != INVALID; ++a) {
+      check(G.arcFromId(G.id(a)) == a, "Wrong id");
+      check(values.find(G.id(a)) == values.end(), "Wrong id");
+      check(G.id(a) <= G.maxArcId(), "Wrong maximum id");
+      values.insert(G.id(a));
+    }
+  }
+
+  template <typename Graph>
+  void checkEdgeIds(const Graph& G) {
+    std::set<int> values;
+    for (typename Graph::EdgeIt e(G); e != INVALID; ++e) {
+      check(G.edgeFromId(G.id(e)) == e, "Wrong id");
+      check(values.find(G.id(e)) == values.end(), "Wrong id");
+      check(G.id(e) <= G.maxEdgeId(), "Wrong maximum id");
+      values.insert(G.id(e));
+    }
+  }
+
+  template <typename Graph>
+  void checkFaceIds(const Graph& G) {
+    std::set<int> values;
+    for (typename Graph::FaceIt n(G); n != INVALID; ++n) {
+      check(G.faceFromId(G.id(n)) == n, "Wrong id");
+      check(values.find(G.id(n)) == values.end(), "Wrong id");
+      check(G.id(n) <= G.maxFaceId(), "Wrong maximum id");
+      values.insert(G.id(n));
+    }
+  }
+
+  template <typename Graph>
+  void checkGraphNodeMap(const Graph& G) {
+    typedef typename Graph::Node Node;
+    typedef typename Graph::NodeIt NodeIt;
+
+    typedef typename Graph::template NodeMap<int> IntNodeMap;
+    IntNodeMap map(G, 42);
+    for (NodeIt it(G); it != INVALID; ++it) {
+      check(map[it] == 42, "Wrong map constructor.");
+    }
+    int s = 0;
+    for (NodeIt it(G); it != INVALID; ++it) {
+      map[it] = 0;
+      check(map[it] == 0, "Wrong operator[].");
+      map.set(it, s);
+      check(map[it] == s, "Wrong set.");
+      ++s;
+    }
+    s = s * (s - 1) / 2;
+    for (NodeIt it(G); it != INVALID; ++it) {
+      s -= map[it];
+    }
+    check(s == 0, "Wrong sum.");
+
+    // map = constMap<Node>(12);
+    // for (NodeIt it(G); it != INVALID; ++it) {
+    //   check(map[it] == 12, "Wrong operator[].");
+    // }
+  }
+
+  template <typename Graph>
+  void checkGraphArcMap(const Graph& G) {
+    typedef typename Graph::Arc Arc;
+    typedef typename Graph::ArcIt ArcIt;
+
+    typedef typename Graph::template ArcMap<int> IntArcMap;
+    IntArcMap map(G, 42);
+    for (ArcIt it(G); it != INVALID; ++it) {
+      check(map[it] == 42, "Wrong map constructor.");
+    }
+    int s = 0;
+    for (ArcIt it(G); it != INVALID; ++it) {
+      map[it] = 0;
+      check(map[it] == 0, "Wrong operator[].");
+      map.set(it, s);
+      check(map[it] == s, "Wrong set.");
+      ++s;
+    }
+    s = s * (s - 1) / 2;
+    for (ArcIt it(G); it != INVALID; ++it) {
+      s -= map[it];
+    }
+    check(s == 0, "Wrong sum.");
+
+    // map = constMap<Arc>(12);
+    // for (ArcIt it(G); it != INVALID; ++it) {
+    //   check(map[it] == 12, "Wrong operator[].");
+    // }
+  }
+
+  template <typename Graph>
+  void checkGraphEdgeMap(const Graph& G) {
+    typedef typename Graph::Edge Edge;
+    typedef typename Graph::EdgeIt EdgeIt;
+
+    typedef typename Graph::template EdgeMap<int> IntEdgeMap;
+    IntEdgeMap map(G, 42);
+    for (EdgeIt it(G); it != INVALID; ++it) {
+      check(map[it] == 42, "Wrong map constructor.");
+    }
+    int s = 0;
+    for (EdgeIt it(G); it != INVALID; ++it) {
+      map[it] = 0;
+      check(map[it] == 0, "Wrong operator[].");
+      map.set(it, s);
+      check(map[it] == s, "Wrong set.");
+      ++s;
+    }
+    s = s * (s - 1) / 2;
+    for (EdgeIt it(G); it != INVALID; ++it) {
+      s -= map[it];
+    }
+    check(s == 0, "Wrong sum.");
+
+    // map = constMap<Edge>(12);
+    // for (EdgeIt it(G); it != INVALID; ++it) {
+    //   check(map[it] == 12, "Wrong operator[].");
+    // }
+  }
+
+
+  template <typename Graph>
+  void checkGraphFaceMap(const Graph& G) {
+    typedef typename Graph::Face Face;
+    typedef typename Graph::FaceIt FaceIt;
+
+    typedef typename Graph::template FaceMap<int> IntFaceMap;
+    IntFaceMap map(G, 42);
+    for (FaceIt it(G); it != INVALID; ++it) {
+      check(map[it] == 42, "Wrong map constructor.");
+    }
+    int s = 0;
+    for (FaceIt it(G); it != INVALID; ++it) {
+      map[it] = 0;
+      check(map[it] == 0, "Wrong operator[].");
+      map.set(it, s);
+      check(map[it] == s, "Wrong set.");
+      ++s;
+    }
+    s = s * (s - 1) / 2;
+    for (FaceIt it(G); it != INVALID; ++it) {
+      s -= map[it];
+    }
+    check(s == 0, "Wrong sum.");
+
+    // map = constMap<Node>(12);
+    // for (NodeIt it(G); it != INVALID; ++it) {
+    //   check(map[it] == 12, "Wrong operator[].");
+    // }
+  }
+
+} //namespace lemon
+
+#endif