Ticket #363: lemon_rev964.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 \

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-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_PLANAR_GRAPH_H
+#define LEMON_PLANAR_GRAPH_H
+
+///\ingroup graphs
+///\file
+///\brief ListDigraph and ListGraph classes.
+
+#include <lemon/core.h>
+#include <lemon/error.h>
+#include <lemon/bits/graph_extender.h>
+
+#include <vector>
+#include <list>
+
+#ifdef REMOVE_BEFORE_RELEASE
+#include <iostream>
+#endif
+
+namespace lemon {
+
+    class PlanarDigraph;
+
+    class PlanarDigraphBase {
+
+    protected:
+        struct NodeT {
+            int first_in, first_out;
+            int first_ccw, last_ccw;
+            int prev, next;
+            int component;
+        };
+
+        struct ArcT {
+            int target, source;
+            int prev_in, prev_out;
+            int next_in, next_out;
+
+            int prev_s, prev_t;
+            int next_s, next_t;
+            int left_face, right_face;
+
+            int target_at(bool d) const {
+                return d?target:source;
+            }
+            int &target_at(bool d) {
+                return d?target:source;
+            }
+            int source_at(bool d) const {
+                return d?source:target;
+            }
+            int &source_at(bool d) {
+                return d?source:target;
+            }
+            int prev_at(bool d) const {
+                return d?prev_s:prev_t;
+            }
+            int &prev_at(bool d) {
+                return d?prev_s:prev_t;
+            }
+            int next_at(bool d) const {
+                return d?next_s:next_t;
+            }
+            int &next_at(bool d) {
+                return d?next_s:next_t;
+            }
+            int left_face_at(bool d) const {
+                return d?left_face:right_face;
+            }
+            int &left_face_at(bool d) {
+                return d?left_face:right_face;
+            }
+            int right_face_at(bool d) const {
+                return d?right_face:left_face;
+            }
+            int &right_face_at(bool d) {
+                return d?right_face:left_face;
+            }
+        };
+
+        struct FaceT {
+            int prev, next;
+            int first_arc;
+            bool first_arc_dir;
+        };
+
+        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;
+
+        int component_id;
+
+    public:
+
+        typedef PlanarDigraphBase Digraph;
+
+        class Node {
+            friend class PlanarDigraphBase;
+            friend class PlanarDigraph;
+        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 Arc {
+            friend class PlanarDigraphBase;
+            friend class PlanarDigraph;
+        protected:
+
+            int id;
+            bool dir;   //travel direction during iteration,
+                        //true means we are travelling from source to target
+            explicit Arc(int pid) { id = pid; dir = true; }
+
+        public:
+            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 PlanarDigraphBase;
+            friend class PlanarDigraph;
+        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;}
+        };
+
+
+
+        PlanarDigraphBase()
+            : nodes(), first_node(-1), first_free_node(-1),
+              arcs(), first_free_arc(-1),
+              faces(), first_face(-1), first_free_face(-1),
+              component_id(0) {
+        }
+
+
+        int maxNodeId() const { return nodes.size()-1; }
+        int maxArcId() const { return arcs.size()-1; }
+        int maxFaceId() const { return faces.size()-1; }
+
+        Node source(Arc e) const { return Node(arcs[e.id].source); }
+        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].right_face); }
+
+
+        void first(Node& node) const {
+            node.id = first_node;
+        }
+
+        void next(Node& node) const {
+            node.id = nodes[node.id].next;
+        }
+
+
+        void first(Arc& arc) const {
+            int n;
+            for(n = first_node;
+                n != -1 && nodes[n].first_out == -1;
+                n = nodes[n].next) {}
+            arc.id = (n == -1) ? -1 : nodes[n].first_out;
+            arc.dir = true;
+        }
+
+        void next(Arc& arc) const {
+            if (arcs[arc.id].next_out != -1) {
+                arc.id = arcs[arc.id].next_out;
+            } else {
+                int n;
+                for(n = nodes[arcs[arc.id].source].next;
+                    n != -1 && nodes[n].first_out == -1;
+                    n = nodes[n].next) {}
+                arc.id = (n == -1) ? -1 : nodes[n].first_out;
+            }
+        }
+
+        void firstOut(Arc &e, const Node& v) const {
+            e.id = nodes[v.id].first_out;
+            e.dir = true;
+        }
+        void nextOut(Arc &e) const {
+            e.id=arcs[e.id].next_out;
+        }
+
+        void firstIn(Arc &e, const Node& v) const {
+            e.id = nodes[v.id].first_in;
+            e.dir = false;
+        }
+        void nextIn(Arc &e) const {
+            e.id=arcs[e.id].next_in;
+        }
+
+
+        void firstCcw(Arc &e, const Node &v) const {
+            e.id = nodes[v.id].first_ccw;
+            if (e.id > -1) e.dir = v.id == arcs[e.id].source;
+        }
+        void nextCcw(Arc &e) const {
+            int n = arcs[e.id].source_at(e.dir);
+            e.id = arcs[e.id].next_at(e.dir);
+            if (e.id > -1) e.dir = arcs[e.id].source == n;
+        }
+
+        void firstCw(Arc &e, const Node &v) const {
+            e.id = nodes[v.id].last_ccw;
+            if (e.id > -1) e.dir = v.id == arcs[e.id].source;
+        }
+
+        void nextCw(Arc &e) const {
+            int n = arcs[e.id].source_at(e.dir);
+            e.id = arcs[e.id].prev_at(e.dir);
+            if (e.id > -1) e.dir = arcs[e.id].source == n;
+        }
+
+        void turnLeft(Arc &e) const {
+            int n = arcs[e.id].source_at(e.dir);
+            e.id = arcs[e.id].next_at(e.dir);
+            if (e.id == -1) e.id = nodes[n].first_ccw;
+            if (e.id > -1) e.dir = arcs[e.id].source == n;
+        }
+        void turnRight(Arc &e) const {
+            int n = arcs[e.id].source_at(e.dir);
+            e.id = arcs[e.id].prev_at(e.dir);
+            if (e.id == -1) e.id = nodes[n].last_ccw;
+            if (e.id > -1) e.dir = arcs[e.id].source == n;
+        }
+
+        void first(Face& face) const {
+            face.id = first_face;
+        }
+
+        void next(Face& face) const {
+            face.id = faces[face.id].next;
+        }
+
+        void firstCcwF(Arc &e, const Face &f) const {
+            e.id = faces[f.id].first_arc;
+            e.dir = faces[f.id].first_arc_dir;
+        }
+        void nextCcwF(Arc &e) const {
+            int f = arcs[e.id].left_face_at(e.dir);
+            turnLeftF(e);
+            if (e.id == faces[f].first_arc && e.dir == faces[f].first_arc_dir) e.id = -1;
+        }
+
+        void firstCwF(Arc &e, const Face &f) const {
+            e.id = faces[f.id].first_arc;
+            e.dir = faces[f.id].first_arc_dir;
+        }
+        void nextCwF(Arc &e) const {
+            int f = arcs[e.id].left_face_at(e.dir);
+            turnRightF(e);
+            if (e.id == faces[f].first_arc && e.dir == faces[f].first_arc_dir) e = INVALID;
+        }
+
+        void turnLeftF(Arc &e) const {
+            e.dir = !e.dir;
+            turnRight(e);
+        }
+
+        void turnRightF(Arc &e) const {
+            turnLeft(e);
+            e.dir = !e.dir;
+        }
+
+        void firstInF(Arc &e, const Face &f) const {
+            e.id = faces[f.id].first_arc;
+            e.dir = faces[f.id].first_arc_dir;
+            while (e.id > -1 && arcs[e.id].right_face != f.id) {
+                turnLeft(e);
+            }
+        }
+        void nextInF(Arc &e) const {
+            int f = arcs[e.id].right_face_at(e.dir);
+            do {
+                turnLeft(e);
+                if (e.id == faces[f].first_arc && e.dir == faces[f].first_arc_dir) e.id = -1;
+            } while (e.id > -1 && arcs[e.id].right_face != f);
+        }
+
+        void firstOutF(Arc &e, const Face &f) const {
+            e.id = faces[f.id].first_arc;
+            e.dir = faces[f.id].first_arc_dir;
+            while (e.id > -1 && arcs[e.id].left_face != f.id) {
+                turnRight(e);
+            }
+        }
+        void nextOutF(Arc &e) const {
+            int f = arcs[e.id].left_face_at(e.dir);
+            do {
+                turnLeft(e);
+                if (e.id == faces[f].first_arc && e.dir == faces[f].first_arc_dir) e.id = -1;
+            } while (e.id > -1 && arcs[e.id].left_face != f);
+        }
+
+
+        static int id(Node v) { return v.id; }
+        static int id(Arc 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 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_in != -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;
+
+            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;
+            nodes[n].component = component_id++;
+            if(first_node != -1) nodes[first_node].prev = n;
+            first_node = n;
+            nodes[n].prev = -1;
+
+            nodes[n].first_in = nodes[n].first_out = -1;
+            nodes[n].first_ccw = -1;
+            nodes[n].last_ccw = -1;
+
+            return Node(n);
+        }
+
+        Arc addArc(Node u, Node v, Arc p_u, Arc p_v) {
+            bool ud = p_u.id > -1 && u.id == arcs[p_u.id].source;
+            bool vd = p_v.id > -1 && v.id == arcs[p_v.id].source;
+
+            if (p_u.id > -1 && p_v.id > -1 && arcs[p_u.id].left_face_at(ud) != arcs[p_v.id].left_face_at(vd)
+                && nodes[u.id].component == nodes[v.id].component) return INVALID;
+            int n = addBlankArc();
+
+            arcs[n].source = u.id;
+            arcs[n].target = v.id;
+
+            arcs[n].next_out = nodes[u.id].first_out;
+            if(nodes[u.id].first_out != -1) {
+                arcs[nodes[u.id].first_out].prev_out = n;
+            }
+
+            arcs[n].next_in = nodes[v.id].first_in;
+            if(nodes[v.id].first_in != -1) {
+                arcs[nodes[v.id].first_in].prev_in = n;
+            }
+
+            arcs[n].prev_in = arcs[n].prev_out = -1;
+
+            nodes[u.id].first_out = nodes[v.id].first_in = n;
+
+            arcs[n].prev_s = p_u.id;
+            if (p_u.id == -1) {
+                arcs[n].next_s = nodes[u.id].first_ccw;
+                nodes[u.id].first_ccw = n;
+            } else {
+                arcs[n].next_s = arcs[p_u.id].next_at(ud);
+                arcs[p_u.id].next_at(ud) = n;
+            }
+            if (arcs[n].next_s != -1) {
+                bool pvd = u.id == arcs[arcs[n].next_s].source;
+                arcs[n].prev_s = arcs[arcs[n].next_s].prev_at(pvd);
+                arcs[arcs[n].next_s].prev_at(pvd) = n;
+            } else {
+                nodes[u.id].last_ccw = n;
+            }
+
+            arcs[n].prev_t = p_v.id;
+            if (p_v.id == -1) {
+                arcs[n].next_t = nodes[v.id].first_ccw;
+                nodes[v.id].first_ccw = n;
+            } else {
+                arcs[n].next_t = arcs[p_v.id].next_at(vd);
+                arcs[p_v.id].next_at(vd) = n;
+            }
+            if (arcs[n].next_t != -1) {
+                bool nvd = v.id == arcs[arcs[n].next_t].source;
+                arcs[n].prev_t = arcs[arcs[n].next_t].prev_at(nvd);
+                arcs[arcs[n].next_t].prev_at(nvd) = n;
+            } else {
+                nodes[v.id].last_ccw = n;
+            }
+
+            //Add the extra face, if needed
+            if (p_u.id > -1 && p_v.id > -1) {
+                int oldf = arcs[p_u.id].left_face_at(ud);
+                int oldfb = arcs[p_v.id].left_face_at(vd);
+                arcs[n].left_face = arcs[n].right_face = oldf;
+                Face f = addFace();
+                faces[f.id].first_arc = n;
+                faces[f.id].first_arc_dir = true;
+                faces[oldf].first_arc = n;
+                faces[oldf].first_arc_dir = false;
+                Arc arc(n);
+                wall_paint(arc,f.id,arc);
+                if (nodes[v.id].component != nodes[u.id].component) {
+                    erase(Face(oldf));
+                    erase(Face(oldfb));
+                    int ca = nodes[u.id].component;
+                    int cb = nodes[v.id].component;
+                    int k = first_node;
+                    while (k != -1) {
+                        if (nodes[k].component == cb)
+                            nodes[k].component = ca;
+                        k = nodes[k].next;
+                    }
+                }
+            } else if (p_u.id > -1) {
+                arcs[n].left_face = arcs[n].right_face = arcs[p_u.id].left_face_at(ud);
+                faces[arcs[n].left_face].first_arc = n;
+                faces[arcs[n].left_face].first_arc_dir = true;
+                nodes[v.id].component = nodes[u.id].component;
+            } else if (p_v.id > -1) {
+                arcs[n].left_face = arcs[n].right_face = arcs[p_v.id].left_face_at(vd);
+                faces[arcs[n].left_face].first_arc = n;
+                faces[arcs[n].left_face].first_arc_dir = true;
+                nodes[u.id].component = nodes[v.id].component;
+            } else {    //both prevs are INVALID
+                Face f = addFace();
+                arcs[n].left_face = arcs[n].right_face = f.id;
+                faces[f.id].first_arc = n;
+                faces[f.id].first_arc_dir = true;
+                nodes[v.id].component = nodes[u.id].component;
+            }
+
+            return Arc(n);
+        }
+
+    public:
+
+        void erase(const Node& node) {  //breaks graph unless node is isolated
+            int n = node.id;
+
+            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;
+
+        }
+
+        void erase(const Arc& arc) {
+            int n = arc.id;
+
+            //"retreat" the incident faces' first arcs
+            int fl = arcs[n].left_face;
+            if (faces[fl].first_arc == n) {
+                Arc e(faces[fl].first_arc);
+                e.dir = faces[fl].first_arc;
+                turnLeftF(e);
+                if (e.id == n) turnLeftF(e);
+                faces[fl].first_arc = e.id;
+                faces[fl].first_arc_dir = e.dir;
+            }
+
+            int fr = arcs[n].right_face;
+
+            bool comp_split = false;
+            if (fr != fl) {
+                Arc arc(faces[fr].first_arc);
+                arc.dir = faces[fr].first_arc_dir;
+                wall_paint(arc,fl,arc);
+                erase(Face(fr));
+            } else if ((arcs[n].next_s > -1 || arcs[n].prev_s > -1) && (arcs[n].next_t > -1 || arcs[n].prev_t > -1)) {
+                comp_split = true;
+                Arc arc(n);
+                Arc ed = arc;
+                ed.dir = false;
+                turnRightF(arc);
+                Face f = addFace();
+                wall_paint(arc,f.id,ed);
+                faces[f.id].first_arc = arc.id;
+                faces[f.id].first_arc_dir = arc.dir;
+            }
+
+            unlink_arc_inout(n);
+
+            if (arcs[n].next_s != -1) {
+                arcs[arcs[n].next_s].prev_at(arcs[n].source == arcs[arcs[n].next_s].source) = arcs[n].prev_s;
+            } else {
+                nodes[arcs[n].source].last_ccw = arcs[n].prev_s;
+            }
+            if (arcs[n].prev_s != -1) {
+                arcs[arcs[n].prev_s].next_at(arcs[n].source == arcs[arcs[n].prev_s].source) = arcs[n].next_s;
+            } else {
+                nodes[arcs[n].source].first_ccw = arcs[n].next_s;
+            }
+
+            if (arcs[n].next_t != -1) {
+                arcs[arcs[n].next_t].prev_at(arcs[n].target == arcs[arcs[n].next_t].source) = arcs[n].prev_t;
+            } else {
+                nodes[arcs[n].target].last_ccw = arcs[n].prev_t;
+            }
+            if (arcs[n].prev_t != -1) {
+                arcs[arcs[n].prev_t].next_at(arcs[n].target == arcs[arcs[n].prev_t].source) = arcs[n].next_t;
+            } else {
+                nodes[arcs[n].target].first_ccw = arcs[n].next_t;
+            }
+
+            unlink_arc_id(n);
+
+            if (comp_split) component_relabel(Node(arcs[n].target), component_id++);
+        }
+
+
+        void clear() {
+            arcs.clear();
+            nodes.clear();
+            first_node = first_free_node = first_free_arc = first_face = first_free_face = -1;
+            faces.clear();
+        }
+
+        void reverseArc(Arc a) {
+            int n = a.id;
+            int tmp;
+            if (faces[arcs[n].left_face].first_arc == n)
+                faces[arcs[n].left_face].first_arc_dir = !faces[arcs[n].left_face].first_arc_dir;
+            if (faces[arcs[n].right_face].first_arc == n)
+                faces[arcs[n].right_face].first_arc_dir = !faces[arcs[n].right_face].first_arc_dir;
+
+            if (arcs[n].prev_out > -1)
+                arcs[arcs[n].prev_out].next_out = arcs[n].next_out;
+            else
+                nodes[arcs[n].source].first_out = arcs[n].next_out;
+            if (arcs[n].next_out > -1)
+                arcs[arcs[n].next_out].prev_out = arcs[n].prev_out;
+
+            if (arcs[n].prev_in > -1)
+                arcs[arcs[n].prev_in].next_in = arcs[n].next_in;
+            else
+                nodes[arcs[n].source].first_in = arcs[n].next_in;
+            if (arcs[n].next_in > -1)
+                arcs[arcs[n].next_in].prev_in = arcs[n].prev_in;
+
+            tmp = arcs[n].left_face; arcs[n].left_face = arcs[n].right_face; arcs[n].right_face = tmp;
+            tmp = arcs[n].prev_s; arcs[n].prev_s = arcs[n].prev_t; arcs[n].prev_t = tmp;
+            tmp = arcs[n].next_s; arcs[n].next_s = arcs[n].next_t; arcs[n].next_t = tmp;
+            tmp = arcs[n].source; arcs[n].source = arcs[n].target; arcs[n].target = tmp;
+
+            arcs[n].prev_out = -1;
+            arcs[n].next_out = nodes[arcs[n].source].first_out;
+            nodes[arcs[n].source].first_out = n;
+
+            arcs[n].prev_in = -1;
+            arcs[n].next_in = nodes[arcs[n].target].first_in;
+            nodes[arcs[n].target].first_in = n;
+        }
+
+        Node split(Arc a) {
+            Node v = addNode();
+            nodes[v.id].component = nodes[arcs[a.id].target].component;
+            int b = addBlankArc();
+            arcs[b] = arcs[a.id];
+            arcs[b].source = v.id;
+            arcs[b].next_out = -1;
+            arcs[b].prev_out = -1;
+            arcs[b].next_s = -1;
+            arcs[b].prev_s = a.id;
+
+            if (arcs[b].next_in != -1) {
+                arcs[arcs[b].next_in].prev_in = b;
+            }
+            if (arcs[b].prev_in != -1) {
+                arcs[arcs[b].prev_in].next_in = b;
+            } else {
+                nodes[arcs[b].target].first_in = b;
+            }
+
+            if (arcs[b].next_t != -1) {
+                arcs[arcs[b].next_t].prev_at(arcs[b].target) = b;
+            } else {
+                nodes[arcs[b].target].last_ccw = b;
+            }
+            if (arcs[b].prev_t != -1) {
+                arcs[arcs[b].prev_t].next_at(arcs[b].target) = b;
+            } else {
+                nodes[arcs[b].target].first_ccw = b;
+            }
+
+            arcs[a.id].target = v.id;
+            arcs[a.id].next_in = -1;
+            arcs[a.id].prev_in = -1;
+            arcs[a.id].next_t = b;
+            arcs[a.id].prev_t = -1;
+
+            nodes[v.id].first_in = a.id;
+            nodes[v.id].first_out = b;
+            nodes[v.id].first_ccw = a.id;
+            nodes[v.id].last_ccw = b;
+
+            return v;
+        }
+
+        bool contract(Node n1, Node n2) {
+            Arc fs;
+            firstCcw(fs,n1);
+            while (fs != INVALID && arcs[fs.id].target_at(fs.dir) != n2.id)
+                nextCcw(fs);
+            if (fs == INVALID) return false;
+            Arc bs = fs;
+            turnRight(bs);
+            while (bs != fs && arcs[bs.id].target_at(bs.dir) == n2.id)
+                turnRight(bs);
+            if (bs == fs) {
+                while (nodes[n1.id].first_ccw != -1)
+                    erase(Arc(nodes[n1.id].first_ccw));
+                erase(n1);
+                return true;
+            }
+            turnLeft(bs);
+            while (bs != fs) {
+                Arc b2s = bs;
+                turnLeft(b2s);
+                erase(bs);
+                bs = b2s;
+            }
+
+            int fl = arcs[fs.id].left_face;
+            if (faces[fl].first_arc == fs.id) {
+                Arc e(faces[fl].first_arc);
+                e.dir = faces[fl].first_arc;
+                turnLeftF(e);
+                if (e.id == fs.id) turnLeftF(e);
+                faces[fl].first_arc = e.id;
+                faces[fl].first_arc_dir = e.dir;
+            }
+            int fr = arcs[fs.id].right_face;
+            if (faces[fr].first_arc == fs.id) {
+                Arc e(faces[fr].first_arc);
+                e.dir = faces[fr].first_arc;
+                turnLeftF(e);
+                if (e.id == fs.id) turnLeftF(e);
+                faces[fr].first_arc = e.id;
+                faces[fr].first_arc_dir = e.dir;
+            }
+
+            Arc a;
+            firstIn(a,n2);
+            while (arcs[a.id].next_in != -1) {
+                if (a.id != fs.id) arcs[a.id].target = n1.id;
+                nextIn(a);
+            }
+            if (a.id != fs.id) arcs[a.id].target = n1.id;
+            arcs[a.id].next_in = nodes[n1.id].first_in;
+            nodes[n1.id].first_in = nodes[n2.id].first_in;
+
+            firstOut(a,n2);
+            while (arcs[a.id].next_out != -1) {
+                if (a.id != fs.id) arcs[a.id].source = n1.id;
+                nextOut(a);
+            }
+            if (a.id != fs.id) arcs[a.id].source = n1.id;
+            arcs[a.id].next_out = nodes[n1.id].first_out;
+            nodes[n1.id].first_out = nodes[n2.id].first_out;
+            unlink_arc_inout(fs.id);
+
+            int la = nodes[n2.id].last_ccw;
+            arcs[la].next_at(arcs[la].source == n2.id) = nodes[n2.id].first_ccw;
+            la = nodes[n2.id].first_ccw;
+            arcs[la].prev_at(arcs[la].source == n2.id) = nodes[n2.id].last_ccw;
+            nodes[n2.id].first_ccw = -1;
+            nodes[n2.id].last_ccw = -1;
+
+            Arc ta = fs;
+            turnRight(ta);
+            if (nodes[n1.id].last_ccw == fs.id) nodes[n1.id].last_ccw = ta.id;
+            arcs[ta.id].next_at(ta.dir) = arcs[fs.id].next_at(!fs.dir);
+            ta = fs;
+            turnLeft(ta);
+            if (nodes[n1.id].first_ccw == fs.id) nodes[n1.id].first_ccw = ta.id;
+            arcs[ta.id].prev_at(ta.dir) = arcs[fs.id].prev_at(!fs.dir);
+            ta = fs;
+            ta.dir = !ta.dir;
+            turnLeft(ta);
+            arcs[ta.id].prev_at(ta.dir) = arcs[fs.id].next_at(fs.dir);
+            ta = fs;
+            ta.dir = !ta.dir;
+            turnRight(ta);
+            arcs[ta.id].next_at(ta.dir) = arcs[fs.id].prev_at(fs.dir);
+
+            unlink_arc_id(fs.id);
+            erase(n2);
+            return true;
+        }
+
+
+    protected:
+
+        //traverses a face clockwise, relabelling the arcs with the face id
+        void wall_paint(Arc arc, int f_id, Arc ed) {
+            do {
+                arcs[arc.id].left_face_at(arc.dir) = f_id;
+                turnRightF(arc);
+            } while (arc.id != ed.id || arc.dir != ed.dir);
+        }
+
+        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;
+            while (!q.empty()) {
+                int n = q.front();
+                ns[n] = 2;
+                nodes[n].component = comp_id;
+                q.pop_front();
+                Arc arc;
+                firstCcw(arc,Node(n));
+                while (arc.id > -1) {
+                    int m = arcs[arc.id].target_at(arcs[arc.id].source == n);
+                    if (ns[m] == 0) {
+                        ns[m] = 1;
+                        q.push_back(m);
+                    }
+                    nextCcw(arc);
+                }
+            }
+        }
+
+        void unlink_arc_inout(int n) {
+            if(arcs[n].next_in!=-1) {
+                arcs[arcs[n].next_in].prev_in = arcs[n].prev_in;
+            }
+
+            if(arcs[n].prev_in!=-1) {
+                arcs[arcs[n].prev_in].next_in = arcs[n].next_in;
+            } else {
+                nodes[arcs[n].target].first_in = arcs[n].next_in;
+            }
+
+
+            if(arcs[n].next_out!=-1) {
+                arcs[arcs[n].next_out].prev_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].source].first_out = arcs[n].next_out;
+            }
+
+        }
+
+        void unlink_arc_id(int n) {
+            arcs[n].next_in = first_free_arc;
+            first_free_arc = n;
+            arcs[n].prev_in = -2;
+        }
+
+        int addBlankArc() {
+            int n;
+            if (first_free_arc == -1) {
+                n = arcs.size();
+                arcs.push_back(ArcT());
+            } else {
+                n = first_free_arc;
+                first_free_arc = arcs[n].next_in;
+            }
+            return 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;
+
+        }
+
+
+        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;
+
+            return Face(n);
+        }
+
+#ifdef REMOVE_BEFORE_RELEASE
+    public:
+        void print() {
+            std::cout << "Nodes: " << std::endl;
+            for (int i=0; i<nodes.size(); ++i)
+                std::cout << i << ": fi=" << nodes[i].first_in
+                          << " fo=" << nodes[i].first_out
+                          << " fc=" << nodes[i].first_ccw
+                          << " lc=" << nodes[i].last_ccw
+                          << " pr=" << nodes[i].prev
+                          << " nx=" << nodes[i].next
+                          << " co=" << nodes[i].component
+                          <<std::endl;
+            std::cout << "Arcs: " << std::endl;
+            for (int i=0; i<arcs.size(); ++i) {
+                if (arcs[i].prev_in > -2) {
+                    std::cout << i << ": sc=" << arcs[i].source
+                          << " tg=" << arcs[i].target
+                          << " pi=" << arcs[i].prev_in
+                          << " ni=" << arcs[i].next_in
+                          << " po=" << arcs[i].prev_out
+                          << " no=" << arcs[i].next_out
+                          << " ps=" << arcs[i].prev_s
+                          << " ns=" << arcs[i].next_s
+                          << " pt=" << arcs[i].prev_t
+                          << " nt=" << arcs[i].next_t
+                          << " lf=" << arcs[i].left_face
+                          << " rf=" << arcs[i].right_face
+                          <<std::endl;
+                } else std::cout << i << ": (deleted)" << std::endl;
+            }
+            std::cout << "Faces: " << std::endl;
+            for (int i=0; i<faces.size(); ++i)
+                std::cout << i
+                          << " pr=" << faces[i].prev
+                          << " nx=" << faces[i].next
+                          << " fa=" << faces[i].first_arc
+                          << " fd=" << faces[i].first_arc_dir
+                          << (faces[i].prev > -2? "" : " (deleted)")
+                          <<std::endl;
+        }
+#endif
+
+    };
+
+    template<typename Base>
+    class PlanarDigraphExtender : public Base{
+
+        typedef Base Parent;
+
+            public:
+        typedef PlanarDigraphExtender Digraph;
+
+        PlanarDigraphExtender() {}
+
+        typedef typename Parent::Node Node;
+        typedef typename Parent::Arc Arc;
+        typedef typename Parent::Face Face;
+
+        class FaceIt : public Face {
+            const Digraph* _digraph;
+        public:
+
+            FaceIt() {}
+
+            FaceIt(Invalid i) : Face(i) { }
+
+            explicit FaceIt(const Digraph& digraph) : _digraph(&digraph) {
+                _digraph->first(static_cast<Face&>(*this));
+            }
+
+            FaceIt(const Digraph& digraph, const Face& face)
+                : Face(face), _digraph(&digraph) {}
+
+            FaceIt& operator++() {
+                _digraph->next(*this);
+                return *this;
+            }
+
+        };
+
+
+        class CwPerimeterArcIt : public Arc {
+          const Digraph* _digraph;
+          Face _face;
+        public:
+
+          CwPerimeterArcIt() { }
+
+          CwPerimeterArcIt(Invalid i) : Arc(i) { }
+
+          CwPerimeterArcIt(const Digraph& digraph, const Face& face)
+            : _digraph(&digraph), _face(face) {
+            _digraph->firstCwF(*this, _face);
+          }
+
+          CwPerimeterArcIt(const Digraph& digraph, const Arc& arc)
+            : Arc(arc), _digraph(&digraph) {}
+
+          CwPerimeterArcIt& operator++() {
+            _digraph->nextCwF(*this);
+            return *this;
+          }
+
+        };
+
+        class CcwArcIt : public Arc {
+          const Digraph* _digraph;
+          const Node _node;
+        public:
+
+          CcwArcIt() { }
+
+          CcwArcIt(Invalid i) : Arc(i) { }
+
+          CcwArcIt(const Digraph& digraph, const Node& node)
+            : _digraph(&digraph), _node(node) {
+            _digraph->firstCcw(*this, node);
+          }
+
+          CcwArcIt(const Digraph& digraph, const Arc& arc)
+            : Arc(arc), _digraph(&digraph) {}
+
+          CcwArcIt& operator++() {
+            _digraph->nextCcw(*this);
+            return *this;
+          }
+
+        };
+
+    };
+
+    typedef PlanarDigraphExtender<DigraphExtender<PlanarDigraphBase> > ExtendedPlanarDigraphBase;
+
+    /// \addtogroup graphs
+    /// @{
+
+    ///A general directed graph structure.
+
+    ///\ref ListDigraph is a versatile and fast directed graph
+    ///implementation based on linked lists that are stored in
+    ///\c std::vector structures.
+    ///
+    ///This type fully conforms to the \ref concepts::Digraph "Digraph concept"
+    ///and it also provides several useful additional functionalities.
+    ///Most of its member functions and nested classes are documented
+    ///only in the concept class.
+    ///
+    ///This class provides only linear time counting for nodes and arcs.
+    ///
+    ///\sa concepts::Digraph
+    ///\sa ListGraph
+    class PlanarDigraph : public ExtendedPlanarDigraphBase {
+        typedef ExtendedPlanarDigraphBase Parent;
+
+    private:
+        /// Digraphs are \e not copy constructible. Use DigraphCopy instead.
+        PlanarDigraph(const PlanarDigraph &) :ExtendedPlanarDigraphBase() {};
+        /// \brief Assignment of a digraph to another one is \e not allowed.
+        /// Use DigraphCopy instead.
+        void operator=(const PlanarDigraph &) {}
+    public:
+
+        /// Constructor
+
+        /// Constructor.
+        ///
+        PlanarDigraph() {}
+
+        ///Add a new node to the digraph.
+
+        ///This function adds a new node to the digraph.
+        ///\return The new node.
+        Node addNode() { return Parent::addNode(); }
+
+        ///Add a new arc to the digraph.
+
+        ///This function adds a new arc to the digraph with source node \c s
+        ///and target node \c t.
+        ///\return The new arc.
+        Arc addArc(Node s, Node t, Arc p_u, Arc p_v) {
+            Arc arc = PlanarDigraphBase::addArc(s, t, p_u, p_v);
+            if (arc != INVALID)
+                notifier(Arc()).add(arc);
+            return arc;
+        }
+
+        ///\brief Erase a node from the digraph.
+        ///
+        ///This function erases the given node along with its outgoing and
+        ///incoming arcs from the digraph.
+        ///
+        ///\note All iterators referencing the removed node or the connected
+        ///arcs are invalidated, of course.
+        void erase(Node n) { Parent::erase(n); }
+
+        ///\brief Erase an arc from the digraph.
+        ///
+        ///This function erases the given arc from the digraph.
+        ///
+        ///\note All iterators referencing the removed arc are invalidated,
+        ///of course.
+        void erase(Arc a) { Parent::erase(a); }
+
+        /// Node validity check
+
+        /// This function gives back \c true if the given node is valid,
+        /// i.e. it is a real node of the digraph.
+        ///
+        /// \warning A removed node could become valid again if new nodes are
+        /// added to the digraph.
+        bool valid(Node n) const { return Parent::valid(n); }
+
+        /// Arc validity check
+
+        /// This function gives back \c true if the given arc is valid,
+        /// i.e. it is a real arc of the digraph.
+        ///
+        /// \warning A removed arc could become valid again if new arcs are
+        /// added to the digraph.
+        bool valid(Arc a) const { return Parent::valid(a); }
+
+        /// Change the target node of an arc
+
+        /// Planar graphs don't support changing endpoints of arcs.
+        void changeTarget(Arc a, Node n) = delete;
+        /// Change the source node of an arc
+
+        ///Planar graphs don't support changing endpoints of arcs.
+        void changeSource(Arc a, Node n) = delete;
+
+        /// Reverse the direction of an arc.
+
+        /// This function reverses the direction of the given arc.
+        ///\note \c ArcIt, \c OutArcIt and \c InArcIt iterators referencing
+        ///the changed arc are invalidated.
+        ///
+        ///\warning This functionality cannot be used together with the Snapshot
+        ///feature.
+        void reverseArc(Arc a) {
+            Parent::reverseArc(a);
+        }
+
+        ///Contract two nodes.
+
+        ///This function contracts the given two nodes.
+        ///Node \c v is removed, but instead of deleting its
+        ///incident arcs, they are joined to node \c u.
+        ///If the last parameter \c r is \c true (this is the default value),
+        ///then the newly created loops are removed.
+        ///
+        ///\note The moved arcs are joined to node \c u using changeSource()
+        ///or changeTarget(), thus \c ArcIt and \c OutArcIt iterators are
+        ///invalidated for the outgoing arcs of node \c v and \c InArcIt
+        ///iterators are invalidated for the incomming arcs of \c v.
+        ///Moreover all iterators referencing node \c v or the removed
+        ///loops are also invalidated. Other iterators remain valid.
+        ///
+        ///\warning This functionality cannot be used together with the Snapshot
+        ///feature.
+        bool contract(Node u, Node v)
+        {
+            return Parent::contract(u,v);
+        }
+
+        ///Split a node.
+
+        ///Planar digraphs don't support splitting of nodes.
+        Node split(Node n, bool connect = true) = delete;
+/*        Node split(Node n, bool connect = true) {
+            Node b = addNode();
+            nodes[b.id].first_out=nodes[n.id].first_out;
+            nodes[n.id].first_out=-1;
+            for(int i=nodes[b.id].first_out; i!=-1; i=arcs[i].next_out) {
+                arcs[i].source=b.id;
+            }
+            if (connect) addArc(n,b);
+            return b;
+        }*/
+
+        ///Split an arc.
+
+        ///This function splits the given arc. First, a new node \c v is
+        ///added to the digraph, then the target node of the original arc
+        ///is set to \c v. Finally, an arc from \c v to the original target
+        ///is added.
+        ///\return The newly created node.
+        ///
+        ///\note \c InArcIt iterators referencing the original arc are
+        ///invalidated. Other iterators remain valid.
+        ///
+        ///\warning This functionality cannot be used together with the
+        ///Snapshot feature.
+        Node split(Arc a) {
+            return Parent::split(a);
+        }
+
+        ///Clear the digraph.
+
+        ///This function erases all nodes and arcs from the digraph.
+        ///
+        ///\note All iterators of the digraph 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 digraph you want to build will
+        /// be large (e.g. it will contain millions of nodes and/or arcs),
+        /// then it is worth reserving space for this amount before starting
+        /// to build the digraph.
+        /// \sa reserveArc()
+        void reserveNode(int n) { nodes.reserve(n); };
+
+        /// Reserve memory for arcs.
+
+        /// Using this function, it is possible to avoid superfluous memory
+        /// allocation: if you know that the digraph you want to build will
+        /// be large (e.g. it will contain millions of nodes and/or arcs),
+        /// then it is worth reserving space for this amount before starting
+        /// to build the digraph.
+        /// \sa reserveNode()
+        void reserveArc(int m) { arcs.reserve(m); };
+
+        class DualBase {
+            const Digraph *_digraph;
+        protected:
+            void initialize(const Digraph *digraph) { _digraph = digraph; }
+        public:
+
+            typedef PlanarDigraph::Face Node;
+            typedef PlanarDigraph::Arc Arc;
+            typedef PlanarDigraph::Node Face;
+
+            int maxNodeId() const { return _digraph->maxFaceId(); }
+            int maxArcId() const { return _digraph->maxArcId(); }
+            int maxFaceId() const { return _digraph->maxNodeId(); }
+
+            Node source(Arc e) const { return _digraph->leftFace(e); }
+            Node target(Arc e) const { return _digraph->rightFace(e); }
+            Face leftFace(Arc e) const { return _digraph->target(e); }
+            Face rightFace(Arc e) const { return _digraph->source(e); }
+
+            void first(Node &i) const { _digraph->first(i); }
+            void next(Node &i) const { _digraph->next(i); }
+            void first(Arc &i) const { _digraph->first(i); }
+            void next(Arc &i) const { _digraph->next(i); }
+            void firstIn(Arc& i, const Node& n) const { _digraph->firstInF(i, n); }
+            void nextIn(Arc& i) const { _digraph->nextInF(i); }
+            void firstOut(Arc& i, const Node& n ) const { _digraph->firstOutF(i, n); }
+            void nextOut(Arc& i) const { _digraph->nextOutF(i); }
+            void firstCcw(Arc &e, const Node &v) const { _digraph->firstCcwF(e, v); }
+            void nextCcw(Arc &e) const { _digraph->nextCcwF(e); }
+            void turnLeft(Arc &e, const Node &v) const { _digraph->turnLeftF(e); }
+            void turnRight(Arc &e, const Node &v) const { _digraph->turnRightF(e); }
+            void first(Face &i) const { _digraph->first(i); }
+            void next(Face &i) const { _digraph->next(i); }
+            void firstCwF(Arc &arc, const Face& face) const { _digraph->firstCw(arc, face); }
+            void nextCwF(Arc &arc) const { _digraph->nextCw(arc); }
+
+            static int id(Node v) { return v.id; }
+            static int id(Arc 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 Face faceFromId(int id) { return Face(id);}
+
+            bool valid(Node n) const { return _digraph->valid(n); }
+            bool valid(Arc n) const { return _digraph->valid(n); }
+            bool valid(Face n) const { return _digraph->valid(n); }
+
+        };
+
+        typedef PlanarDigraphExtender<DigraphExtender<DualBase> > ExtendedDualBase;
+        /// \brief Adaptor class for the dual of a planar graph.
+        ///
+        /// Adaptor class for the dual of a planar graph.
+        class Dual : public ExtendedDualBase {
+                public:
+            Dual(const PlanarDigraph &digraph) { initialize(&digraph); }
+
+        };
+
+        /// \brief Class to make a snapshot of the digraph and restore
+        /// it later.
+        ///
+        /// Class to make a snapshot of the digraph and restore it later.
+        ///
+        /// The newly added nodes and arcs 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 arcs again using
+        /// another Snapshot instance.
+        ///
+        /// \warning Node and arc deletions and other modifications (e.g.
+        /// reversing, contracting, splitting arcs or nodes) cannot be
+        /// restored. These events invalidate the snapshot.
+        /// However, the arcs and nodes that were added to the digraph 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 ArcObserverProxy : public ArcNotifier::ObserverBase {
+            public:
+
+                ArcObserverProxy(Snapshot& _snapshot)
+                    : snapshot(_snapshot) {}
+
+                using ArcNotifier::ObserverBase::attach;
+                using ArcNotifier::ObserverBase::detach;
+                using ArcNotifier::ObserverBase::attached;
+
+            protected:
+
+                virtual void add(const Arc& arc) {
+                    snapshot.addArc(arc);
+                }
+                virtual void add(const std::vector<Arc>& arcs) {
+                    for (int i = arcs.size() - 1; i >= 0; ++i) {
+                        snapshot.addArc(arcs[i]);
+                    }
+                }
+                virtual void erase(const Arc& arc) {
+                    snapshot.eraseArc(arc);
+                }
+                virtual void erase(const std::vector<Arc>& arcs) {
+                    for (int i = 0; i < int(arcs.size()); ++i) {
+                        snapshot.eraseArc(arcs[i]);
+                    }
+                }
+                virtual void build() {
+                    Arc arc;
+                    std::vector<Arc> arcs;
+                    for (notifier()->first(arc); arc != INVALID;
+                    notifier()->next(arc)) {
+                        arcs.push_back(arc);
+                    }
+                    for (int i = arcs.size() - 1; i >= 0; --i) {
+                        snapshot.addArc(arcs[i]);
+                    }
+                }
+                virtual void clear() {
+                    Arc arc;
+                    for (notifier()->first(arc); arc != INVALID;
+                    notifier()->next(arc)) {
+                        snapshot.eraseArc(arc);
+                    }
+                }
+
+                Snapshot& snapshot;
+            };
+
+            PlanarDigraph *digraph;
+
+            NodeObserverProxy node_observer_proxy;
+            ArcObserverProxy arc_observer_proxy;
+
+            std::list<Node> added_nodes;
+            std::list<Arc> added_arcs;
+
+
+            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();
+                    arc_observer_proxy.detach();
+                    throw NodeNotifier::ImmediateDetach();
+                } else {
+                    added_nodes.erase(it);
+                }
+            }
+
+            void addArc(const Arc& arc) {
+                added_arcs.push_front(arc);
+            }
+            void eraseArc(const Arc& arc) {
+                std::list<Arc>::iterator it =
+                        std::find(added_arcs.begin(), added_arcs.end(), arc);
+                if (it == added_arcs.end()) {
+                    clear();
+                    node_observer_proxy.detach();
+                    throw ArcNotifier::ImmediateDetach();
+                } else {
+                    added_arcs.erase(it);
+                }
+            }
+
+            void attach(PlanarDigraph &_digraph) {
+                digraph = &_digraph;
+                node_observer_proxy.attach(digraph->notifier(Node()));
+                arc_observer_proxy.attach(digraph->notifier(Arc()));
+            }
+
+            void detach() {
+                node_observer_proxy.detach();
+                arc_observer_proxy.detach();
+            }
+
+            bool attached() const {
+                return node_observer_proxy.attached();
+            }
+
+            void clear() {
+                added_nodes.clear();
+                added_arcs.clear();
+            }
+
+        public:
+
+            /// \brief Default constructor.
+            ///
+            /// Default constructor.
+            /// You have to call save() to actually make a snapshot.
+            Snapshot()
+                : digraph(0), node_observer_proxy(*this),
+                arc_observer_proxy(*this) {}
+
+            /// \brief Constructor that immediately makes a snapshot.
+            ///
+            /// This constructor immediately makes a snapshot of the given digraph.
+            Snapshot(PlanarDigraph &gr)
+                : node_observer_proxy(*this),
+                arc_observer_proxy(*this) {
+                attach(gr);
+            }
+
+            /// \brief Make a snapshot.
+            ///
+            /// This function makes a snapshot of the given digraph.
+            /// It can be called more than once. In case of a repeated
+            /// call, the previous snapshot gets lost.
+            void save(PlanarDigraph &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(ListDigraph&).
+            ///
+            /// \warning This method invalidates the snapshot, i.e. repeated
+            /// restoring is not supported unless you call save() again.
+            void restore() {
+                detach();
+                for(std::list<Arc>::iterator it = added_arcs.begin();
+                it != added_arcs.end(); ++it) {
+                    digraph->erase(*it);
+                }
+                for(std::list<Node>::iterator it = added_nodes.begin();
+                it != added_nodes.end(); ++it) {
+                    digraph->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();
+            }
+        };
+
+    };
+
+    ///@}
+
+    class PlanarGraphBase {
+
+    protected:
+
+        struct NodeT {
+            int first_out, last_out;
+            int prev, next;
+            int component;
+        };
+
+        struct ArcT {
+            int target;
+            int left_face;
+            int prev_out, next_out;
+        };
+
+        struct FaceT {
+            int first_arc;
+            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;
+
+        int component_id;
+
+    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),
+              component_id(0) {
+        }
+
+
+        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; }
+
+        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); }
+
+        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 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 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;
+
+            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;
+            nodes[n].component = component_id++;
+            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 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 = addBlankEdge();
+
+            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;
+            }
+
+            //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));
+                    int ca = nodes[u.id].component;
+                    int cb = nodes[v.id].component;
+                    int k = first_node;
+                    while (k != -1) {
+                        if (nodes[k].component == cb)
+                            nodes[k].component = ca;
+                        k = nodes[k].next;
+                    }
+                }
+            } 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;
+            } 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;
+                nodes[u.id].component = nodes[v.id].component;
+            } else {    //both prevs are INVALID
+                Face f = addFace();
+                arcs[n].left_face = arcs[n | 1].left_face = f.id;
+                faces[f.id].first_arc = n | 1;
+                nodes[v.id].component = nodes[u.id].component;
+            }
+
+            return Edge(n / 2);
+        }
+
+        void erase(const Node& node) {
+            int n = node.id;
+
+            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;
+        }
+
+        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);
+                erase(Face(fr));
+            } else if ((arcs[n].next_out > -1 || arcs[n].prev_out > -1) &&
+                (arcs[n | 1].next_out > -1 || arcs[n | 1].prev_out > -1)) {
+                comp_split = true;
+                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;
+            }
+
+            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;
+            }
+
+            arcs[n].next_out = first_free_arc;
+            first_free_arc = n;
+            arcs[n].prev_out = -2;
+            arcs[n | 1].prev_out = -2;
+
+            if (comp_split) component_relabel(Node(arcs[n | 1].target), component_id++);
+
+        }
+
+        void clear() {
+            arcs.clear();
+            nodes.clear();
+            faces.clear();
+            first_node = first_free_node = first_free_arc = first_face = first_free_face = -1;
+        }
+
+        Node split(Edge e) {
+            Node v = addNode();
+            Arc a(e.id*2);
+            int b = addBlankEdge();
+
+            nodes[v.id].component = nodes[arcs[a.id].target].component;
+            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;
+        }
+
+    protected:
+
+        void wall_paint(Arc arc, int f_id, Arc ed) {
+            do {
+                arcs[arc.id].left_face = f_id;
+                turnRightF(arc);
+            } while (arc != ed);
+        }
+
+        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;
+            while (!q.empty()) {
+                int n = q.front();
+                ns[n] = 2;
+                nodes[n].component = comp_id;
+                q.pop_front();
+                Arc arc;
+                firstOut(arc,Node(n));
+                while (arc.id > -1) {
+                    int m = arcs[arc.id].target;
+                    if (ns[m] == 0) {
+                        ns[m] = 1;
+                        q.push_back(m);
+                    }
+                    nextOut(arc);
+                }
+            }
+        }
+
+        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 addBlankEdge() {
+            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;
+        }
+
+#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
+                          << " pr=" << nodes[i].prev
+                          << " nx=" << nodes[i].next
+                          << " co=" << nodes[i].component
+                          <<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)
+                std::cout << i
+                          << " pr=" << faces[i].prev
+                          << " nx=" << faces[i].next
+                          << " fa=" << faces[i].first_arc
+                          <<std::endl;
+        }
+#endif
+
+    };
+
+    template<typename Base>
+    class PlanarGraphExtender : public Base{
+
+        typedef Base Parent;
+
+            public:
+        typedef PlanarGraphExtender Graph;
+
+        PlanarGraphExtender() {}
+
+        typedef typename Parent::Node Node;
+        typedef typename Parent::Arc Arc;
+        typedef typename Parent::Edge Edge;
+        typedef typename Parent::Face Face;
+
+        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;
+            }
+
+        };
+
+
+        class CwPerimeterArcIt : public Arc {
+          const Graph* _graph;
+          Face _face;
+          Arc f_arc;
+        public:
+
+          CwPerimeterArcIt() { }
+
+          CwPerimeterArcIt(Invalid i) : Arc(i) { }
+
+          CwPerimeterArcIt(const Graph& graph, const Face& face)
+            : _graph(&graph), _face(face) {
+            _graph->firstCwF(static_cast<Arc&>(*this),face);
+            f_arc = *this;
+          }
+
+          CwPerimeterArcIt(const Graph& graph, const Arc& arc)
+            : Arc(arc), _graph(&graph) {}
+
+          CwPerimeterArcIt& operator++() {
+            _graph->nextCwF(*this);
+            return *this;
+          }
+
+        };
+
+        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;
+          }
+
+        };
+
+    };
+
+    typedef PlanarGraphExtender<GraphExtender<PlanarGraphBase> > ExtendedPlanarGraphBase;
+
+
+    /// \addtogroup graphs
+    /// @{
+
+    ///A general undirected graph structure.
+
+    ///\ref ListGraph is a versatile and fast undirected graph
+    ///implementation based on linked lists that are stored in
+    ///\c std::vector structures.
+    ///
+    ///This type fully conforms to the \ref concepts::Graph "Graph concept"
+    ///and it also provides several useful additional functionalities.
+    ///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 and arcs.
+    ///
+    ///\sa concepts::Graph
+    ///\sa ListDigraph
+    class PlanarGraph : public ExtendedPlanarGraphBase {
+        typedef ExtendedPlanarGraphBase Parent;
+
+    private:
+        /// Graphs are \e not copy constructible. Use GraphCopy instead.
+        PlanarGraph(const PlanarGraph &) :ExtendedPlanarGraphBase()  {};
+        /// \brief Assignment of a graph to another one is \e not allowed.
+        /// Use GraphCopy instead.
+        void operator=(const PlanarGraph &) {}
+    public:
+        /// Constructor
+
+        /// Constructor.
+        ///
+        PlanarGraph() {}
+
+        typedef Parent::OutArcIt IncEdgeIt;
+
+        /// \brief Add a new node to the graph.
+        ///
+        /// This function adds a new node to the graph.
+        /// \return The new node.
+        Node addNode() { return Parent::addNode(); }
+
+        /// \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.
+        /// \return The new edge.
+        Edge addEdge(Node u, Node v, Edge p_u, Edge p_v) {
+            return PlanarGraphBase::addEdge(u, v, p_u, p_v);
+        }
+
+        ///\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) { Parent::erase(n); }
+
+        ///\brief Erase an edge from the graph.
+        ///
+        /// This function erases the given edge from the graph.
+        ///
+        /// \note All iterators referencing the removed edge are invalidated,
+        /// of course.
+        void erase(Edge e) { Parent::erase(e); }
+        /// 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); }
+
+        /// \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 given two nodes.
+        /// Node \c b is removed, but instead of deleting
+        /// its incident edges, they are joined to node \c a.
+        /// If the last parameter \c r is \c true (this is the default value),
+        /// then the newly created loops are removed.
+        ///
+        /// \note The moved edges are joined to node \c a using changeU()
+        /// or changeV(), thus all edge and arc iterators whose base node is
+        /// \c b are invalidated.
+        /// Moreover all iterators referencing node \c b or the removed
+        /// loops are also invalidated. Other iterators remain valid.
+        ///
+        ///\warning This functionality cannot be used together with the
+        ///Snapshot feature.
+/*TODO: rewrite this function
+        void contract(Node a, Node b, bool r = true) {
+            for(IncEdgeIt e(*this, b); e!=INVALID;) {
+                IncEdgeIt f = e; ++f;
+                if (r && runningNode(e) == a) {
+                    erase(e);
+                } else if (u(e) == b) {
+                    changeU(e, a);
+                } else {
+                    changeV(e, a);
+                }
+                e = f;
+            }
+            erase(b);
+        }
+*/
+
+        ///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); };
+
+        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;
+        /// \brief Adaptor class for the dual of a planar graph.
+        ///
+        /// 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(ListGraph&).
+            ///
+            /// \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-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_PLANE_GRAPH_H
+#define LEMON_PLANE_GRAPH_H
+
+///\ingroup graphs
+///\file
+///\brief PlaneDigraph and PlaneGraph classes.
+
+#include <lemon/planar_graph.h>
+#include <lemon/dim2.h>
+#include <cmath>
+#include <lemon/math.h>
+
+namespace lemon {
+
+    class PlaneDigraph : public PlanarDigraph {
+
+    private:
+        NodeMap<dim2::Point<double> > nodepos;
+
+        static double angle(double x, double y) {
+            return atan2(y,x);
+        }
+    public:
+
+        PlaneDigraph() : nodepos(*this) {}
+
+        Node addNode() = delete;
+        Arc addArc(Node, Node, Arc, Arc) = delete;
+
+        Node addNode(const dim2::Point<double> &point) {
+            Node n = PlanarDigraph::addNode();
+            nodepos[n] = point;
+            return n;
+        }
+
+        Arc addArc(Node n1, Node n2) {
+            Arc p_u, p_v;
+            firstCcw(p_u,n1);
+            if (p_u != INVALID) {
+                double ad = angle(nodepos[n2].x-nodepos[n1].x,nodepos[n2].y-nodepos[n1].y);
+                double pd, nd;
+                pd = angle(nodepos[oppositeNode(n1,p_u)].x-nodepos[n1].x,nodepos[oppositeNode(n1,p_u)].y-nodepos[n1].y);
+                nd = pd;
+                if (ad < pd) ad += 2*PI;
+                Arc n_u = p_u;
+                while (n_u != INVALID && ad > nd) {
+                    p_u = n_u;
+                    pd = nd;
+                    nextCcw(n_u,n1);
+                    if (n_u != INVALID)
+                        nd = angle(nodepos[oppositeNode(n1,n_u)].x-nodepos[n1].x,nodepos[oppositeNode(n1,n_u)].y-nodepos[n1].y);
+                }
+            }
+
+            firstCcw(p_v,n2);
+            if (p_v != INVALID) {
+                double ad = angle(nodepos[n1].x-nodepos[n2].x,nodepos[n1].y-nodepos[n2].y);
+                double pd, nd;
+                pd = angle(nodepos[oppositeNode(n2,p_v)].x-nodepos[n2].x,nodepos[oppositeNode(n2,p_v)].y-nodepos[n2].y);
+                nd = pd;
+                if (ad < pd) ad += 2*PI;
+                Arc n_v = p_v;
+                while (n_v != INVALID && ad > nd) {
+                    p_v = n_v;
+                    pd = nd;
+                    nextCcw(n_v,n2);
+                    if (n_v != INVALID)
+                        nd = angle(nodepos[oppositeNode(n2,n_v)].x-nodepos[n2].x,nodepos[oppositeNode(n2,n_v)].y-nodepos[n2].y);
+                }
+            }
+
+#ifdef REMOVE_BEFORE_RELEASE
+//            std::cout << "AddArc: " << id(n1) << "->" << id(n2) << ", p_u: " << id(p_u) << ", p_v: " << id(p_v) << std::endl;
+#endif
+            return PlanarDigraph::addArc(n1,n2,p_u,p_v);
+
+        }
+
+    };
+
+    class PlaneGraph : public PlanarGraph {
+
+    private:
+        NodeMap<dim2::Point<double> > nodepos;
+
+        static double angle(double x, double y) {
+            return atan2(y,x);
+        }
+    public:
+
+        PlaneGraph() : nodepos(*this) {}
+
+        Node addNode() = delete;
+        Arc addArc(Node, Node, Arc, Arc) = delete;
+
+        Node addNode(const dim2::Point<double> &point) {
+            Node n = PlanarGraph::addNode();
+            nodepos[n] = point;
+            return n;
+        }
+
+        Edge addEdge(Node n1, Node n2) {
+            Arc p_u, p_v;
+            firstOut(p_u,n1);
+            if (p_u != INVALID) {
+                double ad = angle(nodepos[n2].x-nodepos[n1].x,nodepos[n2].y-nodepos[n1].y);
+                double pd, nd;
+                pd = angle(nodepos[target(p_u)].x-nodepos[n1].x,nodepos[target(p_u)].y-nodepos[n1].y);
+                nd = pd;
+                if (ad < pd) ad += 2*PI;
+                Arc n_u = p_u;
+                while (n_u != INVALID && ad > nd) {
+                    p_u = n_u;
+                    pd = nd;
+                    nextOut(n_u);
+                    if (n_u != INVALID)
+                        nd = angle(nodepos[target(n_u)].x-nodepos[n1].x,nodepos[target(n_u)].y-nodepos[n1].y);
+                }
+            }
+
+            firstOut(p_v,n2);
+            if (p_v != INVALID) {
+                double ad = angle(nodepos[n1].x-nodepos[n2].x,nodepos[n1].y-nodepos[n2].y);
+                double pd, nd;
+                pd = angle(nodepos[target(p_v)].x-nodepos[n2].x,nodepos[target(p_v)].y-nodepos[n2].y);
+                nd = pd;
+                if (ad < pd) ad += 2*PI;
+                Arc n_v = p_v;
+                while (n_v != INVALID && ad > nd) {
+                    p_v = n_v;
+                    pd = nd;
+                    nextOut(n_v);
+                    if (n_v != INVALID)
+                        nd = angle(nodepos[target(n_v)].x-nodepos[n2].x,nodepos[target(n_v)].y-nodepos[n2].y);
+                }
+            }
+
+#ifdef REMOVE_BEFORE_RELEASE
+//            std::cout << "AddArc: " << id(n1) << "->" << id(n2) << ", p_u: " << id(p_u) << ", p_v: " << id(p_v) << std::endl;
+#endif
+            return PlanarGraph::addEdge(n1,n2,p_u,p_v);
+
+        }
+
+    };
+
+} //namespace lemon
+
+#endif