Index: lemon/Makefile.am
===================================================================
--- lemon/Makefile.am	(revision 2462)
+++ lemon/Makefile.am	(revision 2480)
@@ -100,4 +100,5 @@
 	lemon/path.h \
 	lemon/path_utils.h \
+	lemon/planarity.h \
 	lemon/polynomial.h \
 	lemon/preflow.h \
Index: lemon/planarity.h
===================================================================
--- lemon/planarity.h	(revision 2480)
+++ lemon/planarity.h	(revision 2480)
@@ -0,0 +1,1815 @@
+/* -*- C++ -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library
+ *
+ * Copyright (C) 2003-2007
+ * 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_PLANARITY_H
+#define LEMON_PLANARITY_H
+
+/// \ingroup graph_prop
+/// \file
+/// \brief Planarity checking, embedding
+
+#include <vector>
+#include <list>
+
+#include <lemon/dfs.h>
+#include <lemon/radix_sort.h>
+#include <lemon/maps.h>
+#include <lemon/path.h>
+
+
+namespace lemon {
+
+  namespace _planarity_bits {
+
+    template <typename UGraph>
+    struct PlanarityVisitor : DfsVisitor<UGraph> {
+
+      typedef typename UGraph::Node Node;
+      typedef typename UGraph::Edge Edge;
+
+      typedef typename UGraph::template NodeMap<Edge> PredMap;
+
+      typedef typename UGraph::template UEdgeMap<bool> TreeMap;
+
+      typedef typename UGraph::template NodeMap<int> OrderMap;
+      typedef std::vector<Node> OrderList;
+
+      typedef typename UGraph::template NodeMap<int> LowMap;
+      typedef typename UGraph::template NodeMap<int> AncestorMap;
+
+      PlanarityVisitor(const UGraph& ugraph,
+		       PredMap& pred_map, TreeMap& tree_map,
+		       OrderMap& order_map, OrderList& order_list,
+		       AncestorMap& ancestor_map, LowMap& low_map)
+	: _ugraph(ugraph), _pred_map(pred_map), _tree_map(tree_map),
+	  _order_map(order_map), _order_list(order_list),
+	  _ancestor_map(ancestor_map), _low_map(low_map) {}
+      
+      void reach(const Node& node) {
+	_order_map[node] = _order_list.size();
+	_low_map[node] = _order_list.size();
+	_ancestor_map[node] = _order_list.size();
+	_order_list.push_back(node);
+      }
+
+      void discover(const Edge& edge) {
+	Node source = _ugraph.source(edge);
+	Node target = _ugraph.target(edge);
+
+	_tree_map[edge] = true;
+	_pred_map[target] = edge;
+      }
+
+      void examine(const Edge& edge) {
+	Node source = _ugraph.source(edge);
+	Node target = _ugraph.target(edge);
+	
+	if (_order_map[target] < _order_map[source] && !_tree_map[edge]) {
+	  if (_low_map[source] > _order_map[target]) {
+	    _low_map[source] = _order_map[target];
+	  }
+	  if (_ancestor_map[source] > _order_map[target]) {
+	    _ancestor_map[source] = _order_map[target];
+	  }
+	}
+      }
+
+      void backtrack(const Edge& edge) {
+	Node source = _ugraph.source(edge);
+	Node target = _ugraph.target(edge);
+
+	if (_low_map[source] > _low_map[target]) {
+	  _low_map[source] = _low_map[target];
+	}
+      }
+
+      const UGraph& _ugraph;
+      PredMap& _pred_map;
+      TreeMap& _tree_map;
+      OrderMap& _order_map;
+      OrderList& _order_list;
+      AncestorMap& _ancestor_map;
+      LowMap& _low_map;
+    };
+
+    template <typename UGraph, bool embedding = true>
+    struct NodeDataNode {
+      int prev, next;
+      int visited;
+      typename UGraph::Edge first;
+      bool inverted;
+    };
+
+    template <typename UGraph>
+    struct NodeDataNode<UGraph, false> {
+      int prev, next;
+      int visited;
+    };
+
+    template <typename UGraph>
+    struct ChildListNode {
+      typedef typename UGraph::Node Node;
+      Node first;
+      Node prev, next;
+    };
+
+    template <typename UGraph>
+    struct EdgeListNode {
+      typename UGraph::Edge prev, next;
+    };
+
+  }
+
+  /// \ingroup  graph_prop
+  ///
+  /// \brief Planarity checking of an undirected simple graph
+  ///
+  /// This class implements the Boyer-Myrvold algorithm for planar
+  /// checking of an undirected graph. This class is a simplified
+  /// version of the PlanarEmbedding algorithm class, and it does
+  /// provide neither embedding nor kuratowski subdivisons.
+  template <typename UGraph>
+  class PlanarityChecking {
+  private:
+    
+    UGRAPH_TYPEDEFS(typename UGraph)
+      
+    const UGraph& _ugraph;
+
+  private:
+
+    typedef typename UGraph::template NodeMap<Edge> PredMap;
+    
+    typedef typename UGraph::template UEdgeMap<bool> TreeMap;
+
+    typedef typename UGraph::template NodeMap<int> OrderMap;
+    typedef std::vector<Node> OrderList;
+
+    typedef typename UGraph::template NodeMap<int> LowMap;
+    typedef typename UGraph::template NodeMap<int> AncestorMap;
+
+    typedef _planarity_bits::NodeDataNode<UGraph> NodeDataNode;
+    typedef std::vector<NodeDataNode> NodeData;
+    
+    typedef _planarity_bits::ChildListNode<UGraph> ChildListNode;
+    typedef typename UGraph::template NodeMap<ChildListNode> ChildLists;
+
+    typedef typename UGraph::template NodeMap<std::list<int> > MergeRoots;
+ 
+    typedef typename UGraph::template NodeMap<bool> EmbedEdge;
+
+  public:
+
+    /// \brief Constructor
+    ///
+    /// \warining The graph should be simple, i.e. parallel and loop edge
+    /// free.
+    PlanarityChecking(const UGraph& ugraph) : _ugraph(ugraph) {}
+
+    /// \brief Runs the algorithm.
+    ///
+    /// Runs the algorithm.  
+    /// \param kuratowski If the parameter is false, then the
+    /// algorithm does not calculate the isolate Kuratowski
+    /// subdivisions.
+    /// \return %True when the graph is planar.
+    bool run(bool kuratowski = true) {
+      typedef _planarity_bits::PlanarityVisitor<UGraph> Visitor;
+
+      PredMap pred_map(_ugraph, INVALID);
+      TreeMap tree_map(_ugraph, false);
+
+      OrderMap order_map(_ugraph, -1);
+      OrderList order_list;
+
+      AncestorMap ancestor_map(_ugraph, -1);
+      LowMap low_map(_ugraph, -1);
+
+      Visitor visitor(_ugraph, pred_map, tree_map,
+		      order_map, order_list, ancestor_map, low_map);
+      DfsVisit<UGraph, Visitor> visit(_ugraph, visitor);
+      visit.run();
+
+      ChildLists child_lists(_ugraph);
+      createChildLists(tree_map, order_map, low_map, child_lists);
+
+      NodeData node_data(2 * order_list.size());
+      
+      EmbedEdge embed_edge(_ugraph, false);
+
+      MergeRoots merge_roots(_ugraph);
+      
+      for (int i = order_list.size() - 1; i >= 0; --i) {
+
+	Node node = order_list[i];
+
+	Node source = node;
+	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
+	  Node target = _ugraph.target(e);
+	  
+	  if (order_map[source] < order_map[target] && tree_map[e]) {
+	    initFace(target, node_data, pred_map, order_map, order_list);
+	  }
+	}
+	
+	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
+	  Node target = _ugraph.target(e);
+	  
+	  if (order_map[source] < order_map[target] && !tree_map[e]) {
+	    embed_edge[target] = true;
+	    walkUp(target, source, i, pred_map, low_map,
+		   order_map, order_list, node_data, merge_roots);
+	  }
+	}
+
+	for (typename MergeRoots::Value::iterator it = 
+	       merge_roots[node].begin(); it != merge_roots[node].end(); ++it) {
+	  int rn = *it;
+	  walkDown(rn, i, node_data, order_list, child_lists, 
+		   ancestor_map, low_map, embed_edge, merge_roots);
+	}
+	merge_roots[node].clear();
+	
+	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
+	  Node target = _ugraph.target(e);
+	  
+	  if (order_map[source] < order_map[target] && !tree_map[e]) {
+	    if (embed_edge[target]) {
+	      return false;
+	    }
+	  }
+	}
+      }
+
+      return true;
+    }
+    
+  private:
+
+    void createChildLists(const TreeMap& tree_map, const OrderMap& order_map,
+			  const LowMap& low_map, ChildLists& child_lists) {
+
+      for (NodeIt n(_ugraph); n != INVALID; ++n) {
+	Node source = n;
+	
+	std::vector<Node> targets;  
+	for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node target = _ugraph.target(e);
+
+	  if (order_map[source] < order_map[target] && tree_map[e]) {
+	    targets.push_back(target);
+	  }
+	}	
+
+	if (targets.size() == 0) {
+	  child_lists[source].first = INVALID;
+	} else if (targets.size() == 1) {
+	  child_lists[source].first = targets[0];
+	  child_lists[targets[0]].prev = INVALID;
+	  child_lists[targets[0]].next = INVALID;
+	} else {
+	  radixSort(targets.begin(), targets.end(), mapFunctor(low_map));
+	  for (int i = 1; i < int(targets.size()); ++i) {
+	    child_lists[targets[i]].prev = targets[i - 1];
+	    child_lists[targets[i - 1]].next = targets[i];
+	  }
+	  child_lists[targets.back()].next = INVALID; 
+	  child_lists[targets.front()].prev = INVALID;
+	  child_lists[source].first = targets.front();
+	}
+      }
+    }
+
+    void walkUp(const Node& node, Node root, int rorder,
+		const PredMap& pred_map, const LowMap& low_map,
+		const OrderMap& order_map, const OrderList& order_list,
+		NodeData& node_data, MergeRoots& merge_roots) {
+
+      int na, nb;
+      bool da, db;
+      
+      na = nb = order_map[node];
+      da = true; db = false;
+      
+      while (true) {
+	
+	if (node_data[na].visited == rorder) break;
+	if (node_data[nb].visited == rorder) break;
+
+	node_data[na].visited = rorder;
+	node_data[nb].visited = rorder;
+
+	int rn = -1;
+
+	if (na >= int(order_list.size())) {
+	  rn = na;
+	} else if (nb >= int(order_list.size())) {
+	  rn = nb;
+	}
+
+	if (rn == -1) {
+	  int nn;
+	  
+	  nn = da ? node_data[na].prev : node_data[na].next;
+	  da = node_data[nn].prev != na;
+	  na = nn;
+	  
+	  nn = db ? node_data[nb].prev : node_data[nb].next;
+	  db = node_data[nn].prev != nb;
+	  nb = nn;
+
+	} else {
+
+	  Node rep = order_list[rn - order_list.size()];
+	  Node parent = _ugraph.source(pred_map[rep]);
+
+	  if (low_map[rep] < rorder) {
+	    merge_roots[parent].push_back(rn);
+	  } else {
+	    merge_roots[parent].push_front(rn);
+	  }
+	  
+	  if (parent != root) {  
+	    na = nb = order_map[parent];
+	    da = true; db = false;
+	  } else {
+	    break;
+	  }
+	}	
+      }
+    }
+
+    void walkDown(int rn, int rorder, NodeData& node_data,
+		  OrderList& order_list, ChildLists& child_lists,
+		  AncestorMap& ancestor_map, LowMap& low_map,
+		  EmbedEdge& embed_edge, MergeRoots& merge_roots) {
+
+      std::vector<std::pair<int, bool> > merge_stack;
+
+      for (int di = 0; di < 2; ++di) {
+	bool rd = di == 0;
+	int pn = rn;
+	int n = rd ? node_data[rn].next : node_data[rn].prev;
+	
+	while (n != rn) {
+	  
+	  Node node = order_list[n];
+	  
+	  if (embed_edge[node]) {
+
+	    // Merging components on the critical path
+	    while (!merge_stack.empty()) {
+
+	      // Component root
+	      int cn = merge_stack.back().first;
+	      bool cd = merge_stack.back().second;
+	      merge_stack.pop_back();
+
+	      // Parent of component
+	      int dn = merge_stack.back().first;
+	      bool dd = merge_stack.back().second;
+	      merge_stack.pop_back();
+
+	      Node parent = order_list[dn];
+
+	      // Erasing from merge_roots
+	      merge_roots[parent].pop_front();
+	      
+	      Node child = order_list[cn - order_list.size()];
+
+	      // Erasing from child_lists
+	      if (child_lists[child].prev != INVALID) {
+		child_lists[child_lists[child].prev].next =
+		  child_lists[child].next;
+	      } else {
+		child_lists[parent].first = child_lists[child].next;
+	      }
+	      
+	      if (child_lists[child].next != INVALID) {
+		child_lists[child_lists[child].next].prev =
+		  child_lists[child].prev;
+	      }
+	      
+	      // Merging external faces
+	      {
+		int en = cn;
+		cn = cd ? node_data[cn].prev : node_data[cn].next;
+		cd = node_data[cn].next == en;
+
+	      }
+
+	      if (cd) node_data[cn].next = dn; else node_data[cn].prev = dn;
+	      if (dd) node_data[dn].prev = cn; else node_data[dn].next = cn;
+
+	    }
+
+	    bool d = pn == node_data[n].prev;
+
+	    if (node_data[n].prev == node_data[n].next && 
+		node_data[n].inverted) {
+	      d = !d;
+	    }
+
+	    // Embedding edge into external face
+	    if (rd) node_data[rn].next = n; else node_data[rn].prev = n;
+	    if (d) node_data[n].prev = rn; else node_data[n].next = rn;
+	    pn = rn;
+
+	    embed_edge[order_list[n]] = false;
+	  }
+
+	  if (!merge_roots[node].empty()) {
+
+	    bool d = pn == node_data[n].prev;
+
+	    merge_stack.push_back(std::make_pair(n, d));
+
+	    int rn = merge_roots[node].front();
+	    
+	    int xn = node_data[rn].next;
+	    Node xnode = order_list[xn];
+	    
+	    int yn = node_data[rn].prev;
+	    Node ynode = order_list[yn];
+	    
+	    bool rd;
+	    if (!external(xnode, rorder, child_lists, ancestor_map, low_map)) {
+	      rd = true;
+	    } else if (!external(ynode, rorder, child_lists, 
+				 ancestor_map, low_map)) {
+	      rd = false;
+	    } else if (pertinent(xnode, embed_edge, merge_roots)) {
+	      rd = true;
+	    } else {
+	      rd = false;
+	    }
+	    
+	    merge_stack.push_back(std::make_pair(rn, rd));
+	    
+	    pn = rn;
+	    n = rd ? xn : yn;	      
+	    	    
+	  } else if (!external(node, rorder, child_lists,
+			       ancestor_map, low_map)) {
+	    int nn = (node_data[n].next != pn ? 
+		      node_data[n].next : node_data[n].prev);
+
+	    bool nd = n == node_data[nn].prev;
+
+	    if (nd) node_data[nn].prev = pn;
+	    else node_data[nn].next = pn; 
+
+	    if (n == node_data[pn].prev) node_data[pn].prev = nn;
+	    else node_data[pn].next = nn;
+
+	    node_data[nn].inverted = 
+	      (node_data[nn].prev == node_data[nn].next && nd != rd);
+
+	    n = nn;
+	  }
+	  else break;
+	  
+	}
+
+	if (!merge_stack.empty() || n == rn) {
+	  break;
+	}
+      }
+    }
+
+    void initFace(const Node& node, NodeData& node_data, 
+		  const PredMap& pred_map, const OrderMap& order_map, 
+		  const OrderList& order_list) {
+      int n = order_map[node];
+      int rn = n + order_list.size();
+      
+      node_data[n].next = node_data[n].prev = rn;
+      node_data[rn].next = node_data[rn].prev = n;
+      
+      node_data[n].visited = order_list.size();
+      node_data[rn].visited = order_list.size();
+      
+    }
+
+    bool external(const Node& node, int rorder,
+		  ChildLists& child_lists, AncestorMap& ancestor_map, 
+		  LowMap& low_map) {
+      Node child = child_lists[node].first;
+
+      if (child != INVALID) {
+	if (low_map[child] < rorder) return true;
+      }
+
+      if (ancestor_map[node] < rorder) return true;
+
+      return false;
+    }
+
+    bool pertinent(const Node& node, const EmbedEdge& embed_edge,
+		   const MergeRoots& merge_roots) {
+      return !merge_roots[node].empty() || embed_edge[node];
+    }
+
+  };
+
+  /// \ingroup graph_prop
+  ///
+  /// \brief Planar embedding of an undirected simple graph
+  ///
+  /// This class implements the Boyer-Myrvold algorithm for planar
+  /// embedding of an undirected graph. The planar embeding is an
+  /// ordering of the outgoing edges in each node, which is a possible
+  /// configuration to draw the graph in the plane. If there is not
+  /// such ordering then the graph contains a \f$ K_5 \f$ (full graph
+  /// with 5 nodes) or an \f$ K_{3,3} \f$ (complete bipartite graph on
+  /// 3 ANode and 3 BNode) subdivision.
+  ///
+  /// The current implementation calculates an embedding or an
+  /// Kuratowski subdivision if the graph is not planar. The running
+  /// time of the algorithm is \f$ O(n) \f$.
+  template <typename UGraph>
+  class PlanarEmbedding {
+  private:
+    
+    UGRAPH_TYPEDEFS(typename UGraph)
+      
+    const UGraph& _ugraph;
+    typename UGraph::template EdgeMap<Edge> _embedding;
+
+    typename UGraph::template UEdgeMap<bool> _kuratowski;
+
+  private:
+
+    typedef typename UGraph::template NodeMap<Edge> PredMap;
+    
+    typedef typename UGraph::template UEdgeMap<bool> TreeMap;
+
+    typedef typename UGraph::template NodeMap<int> OrderMap;
+    typedef std::vector<Node> OrderList;
+
+    typedef typename UGraph::template NodeMap<int> LowMap;
+    typedef typename UGraph::template NodeMap<int> AncestorMap;
+
+    typedef _planarity_bits::NodeDataNode<UGraph> NodeDataNode;
+    typedef std::vector<NodeDataNode> NodeData;
+    
+    typedef _planarity_bits::ChildListNode<UGraph> ChildListNode;
+    typedef typename UGraph::template NodeMap<ChildListNode> ChildLists;
+
+    typedef typename UGraph::template NodeMap<std::list<int> > MergeRoots;
+ 
+    typedef typename UGraph::template NodeMap<Edge> EmbedEdge;
+
+    typedef _planarity_bits::EdgeListNode<UGraph> EdgeListNode;
+    typedef typename UGraph::template EdgeMap<EdgeListNode> EdgeLists;
+
+    typedef typename UGraph::template NodeMap<bool> FlipMap;
+
+    typedef typename UGraph::template NodeMap<int> TypeMap;
+
+    enum IsolatorNodeType {
+      HIGHX = 6, LOWX = 7,
+      HIGHY = 8, LOWY = 9,
+      ROOT = 10, PERTINENT = 11,
+      INTERNAL = 12
+    };
+
+  public:
+
+    /// \brief Constructor
+    ///
+    /// \warining The graph should be simple, i.e. parallel and loop edge
+    /// free.
+    PlanarEmbedding(const UGraph& ugraph)
+      : _ugraph(ugraph), _embedding(_ugraph), _kuratowski(ugraph, false) {}
+
+    /// \brief Runs the algorithm.
+    ///
+    /// Runs the algorithm.  
+    /// \param kuratowski If the parameter is false, then the
+    /// algorithm does not calculate the isolate Kuratowski
+    /// subdivisions.
+    ///\return %True when the graph is planar.
+    bool run(bool kuratowski = true) {
+      typedef _planarity_bits::PlanarityVisitor<UGraph> Visitor;
+
+      PredMap pred_map(_ugraph, INVALID);
+      TreeMap tree_map(_ugraph, false);
+
+      OrderMap order_map(_ugraph, -1);
+      OrderList order_list;
+
+      AncestorMap ancestor_map(_ugraph, -1);
+      LowMap low_map(_ugraph, -1);
+
+      Visitor visitor(_ugraph, pred_map, tree_map,
+		      order_map, order_list, ancestor_map, low_map);
+      DfsVisit<UGraph, Visitor> visit(_ugraph, visitor);
+      visit.run();
+
+      ChildLists child_lists(_ugraph);
+      createChildLists(tree_map, order_map, low_map, child_lists);
+
+      NodeData node_data(2 * order_list.size());
+      
+      EmbedEdge embed_edge(_ugraph, INVALID);
+
+      MergeRoots merge_roots(_ugraph);
+
+      EdgeLists edge_lists(_ugraph);
+
+      FlipMap flip_map(_ugraph, false);
+
+      for (int i = order_list.size() - 1; i >= 0; --i) {
+
+	Node node = order_list[i];
+
+	node_data[i].first = INVALID;
+	
+	Node source = node;
+	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
+	  Node target = _ugraph.target(e);
+	  
+	  if (order_map[source] < order_map[target] && tree_map[e]) {
+	    initFace(target, edge_lists, node_data,
+		      pred_map, order_map, order_list);
+	  }
+	}
+	
+	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
+	  Node target = _ugraph.target(e);
+	  
+	  if (order_map[source] < order_map[target] && !tree_map[e]) {
+	    embed_edge[target] = e;
+	    walkUp(target, source, i, pred_map, low_map,
+		   order_map, order_list, node_data, merge_roots);
+	  }
+	}
+
+	for (typename MergeRoots::Value::iterator it = 
+	       merge_roots[node].begin(); it != merge_roots[node].end(); ++it) {
+	  int rn = *it;
+	  walkDown(rn, i, node_data, edge_lists, flip_map, order_list, 
+		   child_lists, ancestor_map, low_map, embed_edge, merge_roots);
+	}
+	merge_roots[node].clear();
+	
+	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
+	  Node target = _ugraph.target(e);
+	  
+	  if (order_map[source] < order_map[target] && !tree_map[e]) {
+	    if (embed_edge[target] != INVALID) {
+	      if (kuratowski) {
+		isolateKuratowski(e, node_data, edge_lists, flip_map,
+				  order_map, order_list, pred_map, child_lists,
+				  ancestor_map, low_map, 
+				  embed_edge, merge_roots);
+	      }
+	      return false;
+	    }
+	  }
+	}
+      }
+
+      for (int i = 0; i < int(order_list.size()); ++i) {
+
+	mergeRemainingFaces(order_list[i], node_data, order_list, order_map,
+			    child_lists, edge_lists);
+	storeEmbedding(order_list[i], node_data, order_map, pred_map,
+		       edge_lists, flip_map);
+      }
+
+      return true;
+    }
+
+    /// \brief Gives back the successor of an edge
+    ///
+    /// Gives back the successor of an edge. This function makes
+    /// possible to query the cyclic order of the outgoing edges from
+    /// a node.
+    Edge next(const Edge& edge) const {
+      return _embedding[edge];
+    }
+
+    /// \brief Gives back true when the undirected edge is in the
+    /// kuratowski subdivision
+    ///
+    /// Gives back true when the undirected edge is in the kuratowski
+    /// subdivision
+    bool kuratowski(const UEdge& uedge) {
+      return _kuratowski[uedge];
+    }
+
+  private:
+
+    void createChildLists(const TreeMap& tree_map, const OrderMap& order_map,
+			  const LowMap& low_map, ChildLists& child_lists) {
+
+      for (NodeIt n(_ugraph); n != INVALID; ++n) {
+	Node source = n;
+	
+	std::vector<Node> targets;  
+	for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
+	  Node target = _ugraph.target(e);
+
+	  if (order_map[source] < order_map[target] && tree_map[e]) {
+	    targets.push_back(target);
+	  }
+	}	
+
+	if (targets.size() == 0) {
+	  child_lists[source].first = INVALID;
+	} else if (targets.size() == 1) {
+	  child_lists[source].first = targets[0];
+	  child_lists[targets[0]].prev = INVALID;
+	  child_lists[targets[0]].next = INVALID;
+	} else {
+	  radixSort(targets.begin(), targets.end(), mapFunctor(low_map));
+	  for (int i = 1; i < int(targets.size()); ++i) {
+	    child_lists[targets[i]].prev = targets[i - 1];
+	    child_lists[targets[i - 1]].next = targets[i];
+	  }
+	  child_lists[targets.back()].next = INVALID; 
+	  child_lists[targets.front()].prev = INVALID;
+	  child_lists[source].first = targets.front();
+	}
+      }
+    }
+
+    void walkUp(const Node& node, Node root, int rorder,
+		const PredMap& pred_map, const LowMap& low_map,
+		const OrderMap& order_map, const OrderList& order_list,
+		NodeData& node_data, MergeRoots& merge_roots) {
+
+      int na, nb;
+      bool da, db;
+      
+      na = nb = order_map[node];
+      da = true; db = false;
+      
+      while (true) {
+	
+	if (node_data[na].visited == rorder) break;
+	if (node_data[nb].visited == rorder) break;
+
+	node_data[na].visited = rorder;
+	node_data[nb].visited = rorder;
+
+	int rn = -1;
+
+	if (na >= int(order_list.size())) {
+	  rn = na;
+	} else if (nb >= int(order_list.size())) {
+	  rn = nb;
+	}
+
+	if (rn == -1) {
+	  int nn;
+	  
+	  nn = da ? node_data[na].prev : node_data[na].next;
+	  da = node_data[nn].prev != na;
+	  na = nn;
+	  
+	  nn = db ? node_data[nb].prev : node_data[nb].next;
+	  db = node_data[nn].prev != nb;
+	  nb = nn;
+
+	} else {
+
+	  Node rep = order_list[rn - order_list.size()];
+	  Node parent = _ugraph.source(pred_map[rep]);
+
+	  if (low_map[rep] < rorder) {
+	    merge_roots[parent].push_back(rn);
+	  } else {
+	    merge_roots[parent].push_front(rn);
+	  }
+	  
+	  if (parent != root) {  
+	    na = nb = order_map[parent];
+	    da = true; db = false;
+	  } else {
+	    break;
+	  }
+	}	
+      }
+    }
+
+    void walkDown(int rn, int rorder, NodeData& node_data,
+		  EdgeLists& edge_lists, FlipMap& flip_map, 
+		  OrderList& order_list, ChildLists& child_lists,
+		  AncestorMap& ancestor_map, LowMap& low_map,
+		  EmbedEdge& embed_edge, MergeRoots& merge_roots) {
+
+      std::vector<std::pair<int, bool> > merge_stack;
+
+      for (int di = 0; di < 2; ++di) {
+	bool rd = di == 0;
+	int pn = rn;
+	int n = rd ? node_data[rn].next : node_data[rn].prev;
+	
+	while (n != rn) {
+	  
+	  Node node = order_list[n];
+	  
+	  if (embed_edge[node] != INVALID) {
+
+	    // Merging components on the critical path
+	    while (!merge_stack.empty()) {
+
+	      // Component root
+	      int cn = merge_stack.back().first;
+	      bool cd = merge_stack.back().second;
+	      merge_stack.pop_back();
+
+	      // Parent of component
+	      int dn = merge_stack.back().first;
+	      bool dd = merge_stack.back().second;
+	      merge_stack.pop_back();
+
+	      Node parent = order_list[dn];
+
+	      // Erasing from merge_roots
+	      merge_roots[parent].pop_front();
+	      
+	      Node child = order_list[cn - order_list.size()];
+
+	      // Erasing from child_lists
+	      if (child_lists[child].prev != INVALID) {
+		child_lists[child_lists[child].prev].next =
+		  child_lists[child].next;
+	      } else {
+		child_lists[parent].first = child_lists[child].next;
+	      }
+	      
+	      if (child_lists[child].next != INVALID) {
+		child_lists[child_lists[child].next].prev =
+		  child_lists[child].prev;
+	      }
+
+	      // Merging edges + flipping
+	      Edge de = node_data[dn].first;
+	      Edge ce = node_data[cn].first;
+
+	      flip_map[order_list[cn - order_list.size()]] = cd != dd;
+	      if (cd != dd) {
+		std::swap(edge_lists[ce].prev, edge_lists[ce].next);
+		ce = edge_lists[ce].prev;
+		std::swap(edge_lists[ce].prev, edge_lists[ce].next);
+	      }
+
+	      {
+		Edge dne = edge_lists[de].next; 
+		Edge cne = edge_lists[ce].next; 
+
+		edge_lists[de].next = cne;
+		edge_lists[ce].next = dne;
+	      
+		edge_lists[dne].prev = ce;
+		edge_lists[cne].prev = de;
+	      }
+	      	      
+	      if (dd) {
+		node_data[dn].first = ce;
+	      }
+	      
+	      // Merging external faces
+	      {
+		int en = cn;
+		cn = cd ? node_data[cn].prev : node_data[cn].next;
+		cd = node_data[cn].next == en;
+
+ 		if (node_data[cn].prev == node_data[cn].next && 
+		    node_data[cn].inverted) {
+ 		  cd = !cd;
+ 		}
+	      }
+
+	      if (cd) node_data[cn].next = dn; else node_data[cn].prev = dn;
+	      if (dd) node_data[dn].prev = cn; else node_data[dn].next = cn;
+
+	    }
+
+	    bool d = pn == node_data[n].prev;
+
+	    if (node_data[n].prev == node_data[n].next && 
+		node_data[n].inverted) {
+	      d = !d;
+	    }
+
+	    // Add new edge
+	    {
+	      Edge edge = embed_edge[node];
+	      Edge re = node_data[rn].first;
+
+	      edge_lists[edge_lists[re].next].prev = edge;
+	      edge_lists[edge].next = edge_lists[re].next;
+	      edge_lists[edge].prev = re;
+	      edge_lists[re].next = edge;
+
+	      if (!rd) {
+		node_data[rn].first = edge;
+	      }
+
+	      Edge rev = _ugraph.oppositeEdge(edge);
+	      Edge e = node_data[n].first;
+
+	      edge_lists[edge_lists[e].next].prev = rev;
+	      edge_lists[rev].next = edge_lists[e].next;
+	      edge_lists[rev].prev = e;
+	      edge_lists[e].next = rev;
+
+	      if (d) {
+		node_data[n].first = rev;
+	      }
+	      
+	    }
+
+	    // Embedding edge into external face
+	    if (rd) node_data[rn].next = n; else node_data[rn].prev = n;
+	    if (d) node_data[n].prev = rn; else node_data[n].next = rn;
+	    pn = rn;
+
+	    embed_edge[order_list[n]] = INVALID;
+	  }
+
+	  if (!merge_roots[node].empty()) {
+
+	    bool d = pn == node_data[n].prev;
+	    if (node_data[n].prev == node_data[n].next && 
+		node_data[n].inverted) {
+	      d = !d;
+	    }
+
+	    merge_stack.push_back(std::make_pair(n, d));
+
+	    int rn = merge_roots[node].front();
+	    
+	    int xn = node_data[rn].next;
+	    Node xnode = order_list[xn];
+	    
+	    int yn = node_data[rn].prev;
+	    Node ynode = order_list[yn];
+	    
+	    bool rd;
+	    if (!external(xnode, rorder, child_lists, ancestor_map, low_map)) {
+	      rd = true;
+	    } else if (!external(ynode, rorder, child_lists, 
+				 ancestor_map, low_map)) {
+	      rd = false;
+	    } else if (pertinent(xnode, embed_edge, merge_roots)) {
+	      rd = true;
+	    } else {
+	      rd = false;
+	    }
+	    
+	    merge_stack.push_back(std::make_pair(rn, rd));
+	    
+	    pn = rn;
+	    n = rd ? xn : yn;	      
+	    	    
+	  } else if (!external(node, rorder, child_lists,
+			       ancestor_map, low_map)) {
+	    int nn = (node_data[n].next != pn ? 
+		      node_data[n].next : node_data[n].prev);
+
+	    bool nd = n == node_data[nn].prev;
+
+	    if (nd) node_data[nn].prev = pn;
+	    else node_data[nn].next = pn; 
+
+	    if (n == node_data[pn].prev) node_data[pn].prev = nn;
+	    else node_data[pn].next = nn;
+
+	    node_data[nn].inverted = 
+	      (node_data[nn].prev == node_data[nn].next && nd != rd);
+
+	    n = nn;
+	  }
+	  else break;
+	  
+	}
+
+	if (!merge_stack.empty() || n == rn) {
+	  break;
+	}
+      }
+    }
+
+    void initFace(const Node& node, EdgeLists& edge_lists,
+		   NodeData& node_data, const PredMap& pred_map,
+		   const OrderMap& order_map, const OrderList& order_list) {
+      int n = order_map[node];
+      int rn = n + order_list.size();
+      
+      node_data[n].next = node_data[n].prev = rn;
+      node_data[rn].next = node_data[rn].prev = n;
+      
+      node_data[n].visited = order_list.size();
+      node_data[rn].visited = order_list.size();
+
+      node_data[n].inverted = false;
+      node_data[rn].inverted = false;
+
+      Edge edge = pred_map[node];
+      Edge rev = _ugraph.oppositeEdge(edge);
+
+      node_data[rn].first = edge;
+      node_data[n].first = rev;
+
+      edge_lists[edge].prev = edge;
+      edge_lists[edge].next = edge;
+
+      edge_lists[rev].prev = rev;
+      edge_lists[rev].next = rev;
+
+    }
+
+    void mergeRemainingFaces(const Node& node, NodeData& node_data,
+			     OrderList& order_list, OrderMap& order_map,
+			     ChildLists& child_lists, EdgeLists& edge_lists) {
+      while (child_lists[node].first != INVALID) {
+	int dd = order_map[node];
+	Node child = child_lists[node].first; 
+	int cd = order_map[child] + order_list.size();
+	child_lists[node].first = child_lists[child].next;
+
+	Edge de = node_data[dd].first;
+	Edge ce = node_data[cd].first;
+
+	if (de != INVALID) {
+	  Edge dne = edge_lists[de].next; 
+	  Edge cne = edge_lists[ce].next; 
+	  
+	  edge_lists[de].next = cne;
+	  edge_lists[ce].next = dne;
+	  
+	  edge_lists[dne].prev = ce;
+	  edge_lists[cne].prev = de;
+	}
+	
+	node_data[dd].first = ce;
+
+      }
+    }
+
+    void storeEmbedding(const Node& node, NodeData& node_data,
+			OrderMap& order_map, PredMap& pred_map,
+			EdgeLists& edge_lists, FlipMap& flip_map) {
+
+      if (node_data[order_map[node]].first == INVALID) return;
+
+      if (pred_map[node] != INVALID) {
+	Node source = _ugraph.source(pred_map[node]);
+	flip_map[node] = flip_map[node] != flip_map[source];
+      }
+      
+      Edge first = node_data[order_map[node]].first;
+      Edge prev = first;
+
+      Edge edge = flip_map[node] ?
+	edge_lists[prev].prev : edge_lists[prev].next;
+
+      _embedding[prev] = edge;
+      
+      while (edge != first) {
+	Edge next = edge_lists[edge].prev == prev ?
+	  edge_lists[edge].next : edge_lists[edge].prev;
+	prev = edge; edge = next;
+	_embedding[prev] = edge;
+      }
+    }
+
+
+    bool external(const Node& node, int rorder,
+		  ChildLists& child_lists, AncestorMap& ancestor_map, 
+		  LowMap& low_map) {
+      Node child = child_lists[node].first;
+
+      if (child != INVALID) {
+	if (low_map[child] < rorder) return true;
+      }
+
+      if (ancestor_map[node] < rorder) return true;
+
+      return false;
+    }
+
+    bool pertinent(const Node& node, const EmbedEdge& embed_edge,
+		   const MergeRoots& merge_roots) {
+      return !merge_roots[node].empty() || embed_edge[node] != INVALID;
+    }
+
+    int lowPoint(const Node& node, OrderMap& order_map, ChildLists& child_lists,
+		 AncestorMap& ancestor_map, LowMap& low_map) {
+      int low_point;
+
+      Node child = child_lists[node].first;
+
+      if (child != INVALID) {
+	low_point = low_map[child];
+      } else {
+	low_point = order_map[node];
+      }
+
+      if (low_point > ancestor_map[node]) {
+	low_point = ancestor_map[node];
+      }
+
+      return low_point;
+    }
+
+    int findComponentRoot(Node root, Node node, ChildLists& child_lists, 
+			  OrderMap& order_map, OrderList& order_list) {
+
+      int order = order_map[root];
+      int norder = order_map[node];
+
+      Node child = child_lists[root].first;
+      while (child != INVALID) {
+	int corder = order_map[child];
+	if (corder > order && corder < norder) {
+	  order = corder;
+	}
+	child = child_lists[child].next;
+      }
+      return order + order_list.size();
+    }
+
+    Node findPertinent(Node node, OrderMap& order_map, NodeData& node_data,
+		       EmbedEdge& embed_edge, MergeRoots& merge_roots) {
+      Node wnode =_ugraph.target(node_data[order_map[node]].first);
+      while (!pertinent(wnode, embed_edge, merge_roots)) {
+	wnode = _ugraph.target(node_data[order_map[wnode]].first);
+      }
+      return wnode;
+    }
+
+
+    Node findExternal(Node node, int rorder, OrderMap& order_map, 
+		      ChildLists& child_lists, AncestorMap& ancestor_map,
+		      LowMap& low_map, NodeData& node_data) {
+      Node wnode =_ugraph.target(node_data[order_map[node]].first);
+      while (!external(wnode, rorder, child_lists, ancestor_map, low_map)) {
+	wnode = _ugraph.target(node_data[order_map[wnode]].first);
+      }
+      return wnode;
+    }
+
+    void markCommonPath(Node node, int rorder, Node& wnode, Node& znode, 
+			OrderList& order_list, OrderMap& order_map, 
+			NodeData& node_data, EdgeLists& edge_lists, 
+			EmbedEdge& embed_edge, MergeRoots& merge_roots, 
+			ChildLists& child_lists, AncestorMap& ancestor_map, 
+			LowMap& low_map) {
+      
+      Node cnode = node;
+      Node pred = INVALID;
+      
+      while (true) {
+
+	bool pert = pertinent(cnode, embed_edge, merge_roots);
+	bool ext = external(cnode, rorder, child_lists, ancestor_map, low_map);
+
+	if (pert && ext) {
+	  if (!merge_roots[cnode].empty()) {
+	    int cn = merge_roots[cnode].back();
+	    
+	    if (low_map[order_list[cn - order_list.size()]] < rorder) {
+	      Edge edge = node_data[cn].first;
+	      _kuratowski.set(edge, true);
+	      
+	      pred = cnode;
+	      cnode = _ugraph.target(edge);
+	    
+	      continue;
+	    }
+	  }
+	  wnode = znode = cnode;
+	  return;
+
+	} else if (pert) {
+	  wnode = cnode;
+	  
+	  while (!external(cnode, rorder, child_lists, ancestor_map, low_map)) {
+	    Edge edge = node_data[order_map[cnode]].first;
+	  
+	    if (_ugraph.target(edge) == pred) {
+	      edge = edge_lists[edge].next;
+	    }
+	    _kuratowski.set(edge, true);
+	    
+	    Node next = _ugraph.target(edge);
+	    pred = cnode; cnode = next;
+	  }
+	  
+	  znode = cnode;
+	  return;
+
+	} else if (ext) {
+	  znode = cnode;
+	  
+	  while (!pertinent(cnode, embed_edge, merge_roots)) {
+	    Edge edge = node_data[order_map[cnode]].first;
+	  
+	    if (_ugraph.target(edge) == pred) {
+	      edge = edge_lists[edge].next;
+	    }
+	    _kuratowski.set(edge, true);
+	    
+	    Node next = _ugraph.target(edge);
+	    pred = cnode; cnode = next;
+	  }
+	  
+	  wnode = cnode;
+	  return;
+	  
+	} else {
+	  Edge edge = node_data[order_map[cnode]].first;
+	  
+	  if (_ugraph.target(edge) == pred) {
+	    edge = edge_lists[edge].next;
+	  }
+	  _kuratowski.set(edge, true);
+
+	  Node next = _ugraph.target(edge);
+	  pred = cnode; cnode = next;
+	}
+	
+      }
+
+    }
+
+    void orientComponent(Node root, int rn, OrderMap& order_map,
+			 PredMap& pred_map, NodeData& node_data, 
+			 EdgeLists& edge_lists, FlipMap& flip_map, 
+			 TypeMap& type_map) {
+      node_data[order_map[root]].first = node_data[rn].first;
+      type_map[root] = 1;
+
+      std::vector<Node> st, qu;
+
+      st.push_back(root);
+      while (!st.empty()) {
+	Node node = st.back();
+	st.pop_back();
+	qu.push_back(node);
+	
+	Edge edge = node_data[order_map[node]].first;
+	
+	if (type_map[_ugraph.target(edge)] == 0) {
+	  st.push_back(_ugraph.target(edge));
+	  type_map[_ugraph.target(edge)] = 1;
+	} 
+	
+	Edge last = edge, pred = edge;
+	edge = edge_lists[edge].next;
+	while (edge != last) {
+
+	  if (type_map[_ugraph.target(edge)] == 0) {
+	    st.push_back(_ugraph.target(edge));
+	    type_map[_ugraph.target(edge)] = 1;
+	  } 
+	  
+	  Edge next = edge_lists[edge].next != pred ? 
+	    edge_lists[edge].next : edge_lists[edge].prev;
+	  pred = edge; edge = next;
+	}
+
+      }
+
+      type_map[root] = 2;
+      flip_map[root] = false;
+
+      for (int i = 1; i < int(qu.size()); ++i) {
+
+	Node node = qu[i];
+
+	while (type_map[node] != 2) {
+	  st.push_back(node);
+	  type_map[node] = 2;
+	  node = _ugraph.source(pred_map[node]);
+	}
+
+	bool flip = flip_map[node];
+
+	while (!st.empty()) {
+	  node = st.back();
+	  st.pop_back();
+	  
+	  flip_map[node] = flip != flip_map[node];
+	  flip = flip_map[node];
+
+	  if (flip) {
+	    Edge edge = node_data[order_map[node]].first;
+	    std::swap(edge_lists[edge].prev, edge_lists[edge].next);
+	    edge = edge_lists[edge].prev;
+	    std::swap(edge_lists[edge].prev, edge_lists[edge].next);
+	    node_data[order_map[node]].first = edge;
+	  }
+	}
+      }
+
+      for (int i = 0; i < int(qu.size()); ++i) {
+
+	Edge edge = node_data[order_map[qu[i]]].first;
+	Edge last = edge, pred = edge;
+
+	edge = edge_lists[edge].next;
+	while (edge != last) {
+
+	  if (edge_lists[edge].next == pred) {
+	    std::swap(edge_lists[edge].next, edge_lists[edge].prev);
+	  } 
+	  pred = edge; edge = edge_lists[edge].next;
+	}
+	
+      }
+    }
+
+    void setFaceFlags(Node root, Node wnode, Node ynode, Node xnode,
+		      OrderMap& order_map, NodeData& node_data, 
+		      TypeMap& type_map) {
+      Node node = _ugraph.target(node_data[order_map[root]].first);
+
+      while (node != ynode) {
+	type_map[node] = HIGHY;
+	node = _ugraph.target(node_data[order_map[node]].first);
+      }
+
+      while (node != wnode) {
+	type_map[node] = LOWY;
+	node = _ugraph.target(node_data[order_map[node]].first);
+      }
+      
+      node = _ugraph.target(node_data[order_map[wnode]].first);
+      
+      while (node != xnode) {
+	type_map[node] = LOWX;
+	node = _ugraph.target(node_data[order_map[node]].first);
+      }
+      type_map[node] = LOWX;
+
+      node = _ugraph.target(node_data[order_map[xnode]].first);
+      while (node != root) {
+	type_map[node] = HIGHX;
+	node = _ugraph.target(node_data[order_map[node]].first);
+      }
+
+      type_map[wnode] = PERTINENT;
+      type_map[root] = ROOT;
+    }
+
+    void findInternalPath(std::vector<Edge>& ipath,
+			  Node wnode, Node root, TypeMap& type_map, 
+			  OrderMap& order_map, NodeData& node_data, 
+			  EdgeLists& edge_lists) {
+      std::vector<Edge> st;
+
+      Node node = wnode;
+      
+      while (node != root) {
+	Edge edge = edge_lists[node_data[order_map[node]].first].next;
+	st.push_back(edge);
+	node = _ugraph.target(edge);
+      }
+      
+      while (true) {
+	Edge edge = st.back();
+	if (type_map[_ugraph.target(edge)] == LOWX ||
+	    type_map[_ugraph.target(edge)] == HIGHX) {
+	  break;
+	}
+	if (type_map[_ugraph.target(edge)] == 2) {
+	  type_map[_ugraph.target(edge)] = 3;
+	  
+	  edge = edge_lists[_ugraph.oppositeEdge(edge)].next;
+	  st.push_back(edge);
+	} else {
+	  st.pop_back();
+	  edge = edge_lists[edge].next;
+	  
+	  while (_ugraph.oppositeEdge(edge) == st.back()) {
+	    edge = st.back();
+	    st.pop_back();
+	    edge = edge_lists[edge].next;
+	  }
+	  st.push_back(edge);
+	}
+      }
+      
+      for (int i = 0; i < int(st.size()); ++i) {
+	if (type_map[_ugraph.target(st[i])] != LOWY &&
+	    type_map[_ugraph.target(st[i])] != HIGHY) {
+	  for (; i < int(st.size()); ++i) {
+	    ipath.push_back(st[i]);
+	  }
+	}
+      }
+    }
+
+    void setInternalFlags(std::vector<Edge>& ipath, TypeMap& type_map) {
+      for (int i = 1; i < int(ipath.size()); ++i) {
+	type_map[_ugraph.source(ipath[i])] = INTERNAL;
+      }
+    }
+
+    void findPilePath(std::vector<Edge>& ppath,
+		      Node root, TypeMap& type_map, OrderMap& order_map, 
+		      NodeData& node_data, EdgeLists& edge_lists) {
+      std::vector<Edge> st;
+
+      st.push_back(_ugraph.oppositeEdge(node_data[order_map[root]].first));
+      st.push_back(node_data[order_map[root]].first);
+      
+      while (st.size() > 1) {
+	Edge edge = st.back();
+	if (type_map[_ugraph.target(edge)] == INTERNAL) {
+	  break;
+	}
+	if (type_map[_ugraph.target(edge)] == 3) {
+	  type_map[_ugraph.target(edge)] = 4;
+	  
+	  edge = edge_lists[_ugraph.oppositeEdge(edge)].next;
+	  st.push_back(edge);
+	} else {
+	  st.pop_back();
+	  edge = edge_lists[edge].next;
+	  
+	  while (!st.empty() && _ugraph.oppositeEdge(edge) == st.back()) {
+	    edge = st.back();
+	    st.pop_back();
+	    edge = edge_lists[edge].next;
+	  }
+	  st.push_back(edge);
+	}
+      }
+      
+      for (int i = 1; i < int(st.size()); ++i) {
+	ppath.push_back(st[i]);
+      }
+    }
+
+
+    int markExternalPath(Node node, OrderMap& order_map,
+			 ChildLists& child_lists, PredMap& pred_map,
+			 AncestorMap& ancestor_map, LowMap& low_map) {
+      int lp = lowPoint(node, order_map, child_lists,
+			ancestor_map, low_map);
+      
+      if (ancestor_map[node] != lp) {
+	node = child_lists[node].first;
+	_kuratowski[pred_map[node]] = true;
+
+	while (ancestor_map[node] != lp) {
+	  for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
+	    Node tnode = _ugraph.target(e); 
+	    if (order_map[tnode] > order_map[node] && low_map[tnode] == lp) {
+	      node = tnode;
+	      _kuratowski[e] = true;
+	      break;
+	    }
+	  }
+	}
+      }
+
+      for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
+	if (order_map[_ugraph.target(e)] == lp) {
+	  _kuratowski[e] = true;
+	  break;
+	}
+      }
+      
+      return lp;
+    }
+
+    void markPertinentPath(Node node, OrderMap& order_map, 
+			   NodeData& node_data, EdgeLists& edge_lists,
+			   EmbedEdge& embed_edge, MergeRoots& merge_roots) {
+      while (embed_edge[node] == INVALID) {
+	int n = merge_roots[node].front();
+	Edge edge = node_data[n].first;
+
+	_kuratowski.set(edge, true);
+	
+	Node pred = node;
+	node = _ugraph.target(edge);
+	while (!pertinent(node, embed_edge, merge_roots)) {
+	  edge = node_data[order_map[node]].first;
+	  if (_ugraph.target(edge) == pred) {
+	    edge = edge_lists[edge].next;
+	  }
+	  _kuratowski.set(edge, true);
+	  pred = node;
+	  node = _ugraph.target(edge);
+	}
+      }
+      _kuratowski.set(embed_edge[node], true);
+    } 
+
+    void markPredPath(Node node, Node snode, PredMap& pred_map) {
+      while (node != snode) {
+	_kuratowski.set(pred_map[node], true);
+	node = _ugraph.source(pred_map[node]);
+      }
+    }
+
+    void markFacePath(Node ynode, Node xnode, 
+		      OrderMap& order_map, NodeData& node_data) {
+      Edge edge = node_data[order_map[ynode]].first;
+      Node node = _ugraph.target(edge);
+      _kuratowski.set(edge, true);
+	
+      while (node != xnode) {
+	edge = node_data[order_map[node]].first;
+	_kuratowski.set(edge, true);
+	node = _ugraph.target(edge);
+      }
+    }
+
+    void markInternalPath(std::vector<Edge>& path) {
+      for (int i = 0; i < int(path.size()); ++i) {
+	_kuratowski.set(path[i], true);
+      }
+    }
+
+    void markPilePath(std::vector<Edge>& path) {
+      for (int i = 0; i < int(path.size()); ++i) {
+	_kuratowski.set(path[i], true);
+      }
+    }
+
+    void isolateKuratowski(Edge edge, NodeData& node_data, 
+			   EdgeLists& edge_lists, FlipMap& flip_map,
+			   OrderMap& order_map, OrderList& order_list, 
+			   PredMap& pred_map, ChildLists& child_lists,
+			   AncestorMap& ancestor_map, LowMap& low_map, 
+			   EmbedEdge& embed_edge, MergeRoots& merge_roots) {
+
+      Node root = _ugraph.source(edge);
+      Node enode = _ugraph.target(edge);
+
+      int rorder = order_map[root];
+
+      TypeMap type_map(_ugraph, 0);
+
+      int rn = findComponentRoot(root, enode, child_lists, 
+				 order_map, order_list);
+
+      Node xnode = order_list[node_data[rn].next];
+      Node ynode = order_list[node_data[rn].prev];
+
+      // Minor-A
+      {
+	while (!merge_roots[xnode].empty() || !merge_roots[ynode].empty()) {
+	  
+	  if (!merge_roots[xnode].empty()) {
+	    root = xnode;
+	    rn = merge_roots[xnode].front();
+	  } else {
+	    root = ynode;
+	    rn = merge_roots[ynode].front();
+	  }
+	  
+	  xnode = order_list[node_data[rn].next];
+	  ynode = order_list[node_data[rn].prev];
+	}
+	
+	if (root != _ugraph.source(edge)) {
+	  orientComponent(root, rn, order_map, pred_map, 
+			  node_data, edge_lists, flip_map, type_map);
+	  markFacePath(root, root, order_map, node_data);
+	  int xlp = markExternalPath(xnode, order_map, child_lists, 
+				     pred_map, ancestor_map, low_map);
+	  int ylp = markExternalPath(ynode, order_map, child_lists, 
+				     pred_map, ancestor_map, low_map);
+	  markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
+	  Node lwnode = findPertinent(ynode, order_map, node_data,
+				     embed_edge, merge_roots);
+	  
+	  markPertinentPath(lwnode, order_map, node_data, edge_lists,
+			    embed_edge, merge_roots);
+	  
+	  return;
+	}
+      }
+      
+      orientComponent(root, rn, order_map, pred_map, 
+		      node_data, edge_lists, flip_map, type_map);
+
+      Node wnode = findPertinent(ynode, order_map, node_data,
+				 embed_edge, merge_roots);
+      setFaceFlags(root, wnode, ynode, xnode, order_map, node_data, type_map);
+
+      
+      //Minor-B
+      if (!merge_roots[wnode].empty()) {
+	int cn = merge_roots[wnode].back();
+	Node rep = order_list[cn - order_list.size()];
+	if (low_map[rep] < rorder) {
+	  markFacePath(root, root, order_map, node_data);
+	  int xlp = markExternalPath(xnode, order_map, child_lists, 
+				     pred_map, ancestor_map, low_map);
+	  int ylp = markExternalPath(ynode, order_map, child_lists, 
+				     pred_map, ancestor_map, low_map);
+
+	  Node lwnode, lznode;
+	  markCommonPath(wnode, rorder, lwnode, lznode, order_list, 
+			 order_map, node_data, edge_lists, embed_edge, 
+			 merge_roots, child_lists, ancestor_map, low_map);
+	  	  
+	  markPertinentPath(lwnode, order_map, node_data, edge_lists,
+			    embed_edge, merge_roots);
+	  int zlp = markExternalPath(lznode, order_map, child_lists, 
+				     pred_map, ancestor_map, low_map);
+
+	  int minlp = xlp < ylp ? xlp : ylp;
+	  if (zlp < minlp) minlp = zlp;
+
+	  int maxlp = xlp > ylp ? xlp : ylp;
+	  if (zlp > maxlp) maxlp = zlp;
+	  
+	  markPredPath(order_list[maxlp], order_list[minlp], pred_map);
+	  
+	  return;
+	}
+      }
+
+      Node pxnode, pynode;
+      std::vector<Edge> ipath;
+      findInternalPath(ipath, wnode, root, type_map, order_map,
+		       node_data, edge_lists);
+      setInternalFlags(ipath, type_map);
+      pynode = _ugraph.source(ipath.front());
+      pxnode = _ugraph.target(ipath.back());
+
+      wnode = findPertinent(pynode, order_map, node_data,
+			    embed_edge, merge_roots);
+      
+      // Minor-C
+      {
+	if (type_map[_ugraph.source(ipath.front())] == HIGHY) {
+	  if (type_map[_ugraph.target(ipath.back())] == HIGHX) {
+	    markFacePath(xnode, pxnode, order_map, node_data);
+	  }
+	  markFacePath(root, xnode, order_map, node_data);
+	  markPertinentPath(wnode, order_map, node_data, edge_lists,
+			    embed_edge, merge_roots);
+	  markInternalPath(ipath);
+	  int xlp = markExternalPath(xnode, order_map, child_lists, 
+				     pred_map, ancestor_map, low_map);
+	  int ylp = markExternalPath(ynode, order_map, child_lists, 
+				     pred_map, ancestor_map, low_map);
+	  markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
+	  return;
+	}
+
+	if (type_map[_ugraph.target(ipath.back())] == HIGHX) {
+	  markFacePath(ynode, root, order_map, node_data);
+	  markPertinentPath(wnode, order_map, node_data, edge_lists,
+			    embed_edge, merge_roots);
+	  markInternalPath(ipath);
+	  int xlp = markExternalPath(xnode, order_map, child_lists, 
+				     pred_map, ancestor_map, low_map);
+	  int ylp = markExternalPath(ynode, order_map, child_lists, 
+				     pred_map, ancestor_map, low_map);
+	  markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
+	  return;
+	}
+      }
+
+      std::vector<Edge> ppath;
+      findPilePath(ppath, root, type_map, order_map, node_data, edge_lists);
+      
+      // Minor-D
+      if (!ppath.empty()) {
+	markFacePath(ynode, xnode, order_map, node_data);
+	markPertinentPath(wnode, order_map, node_data, edge_lists,
+			  embed_edge, merge_roots);
+	markPilePath(ppath);
+	markInternalPath(ipath);
+	int xlp = markExternalPath(xnode, order_map, child_lists, 
+				   pred_map, ancestor_map, low_map);
+	int ylp = markExternalPath(ynode, order_map, child_lists, 
+				   pred_map, ancestor_map, low_map);
+	markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
+	return;
+      }
+
+      // Minor-E*
+      {
+
+	if (!external(wnode, rorder, child_lists, ancestor_map, low_map)) {
+	  Node znode = findExternal(pynode, rorder, order_map, 
+				    child_lists, ancestor_map,
+				    low_map, node_data);
+	  
+	  if (type_map[znode] == LOWY) {
+	    markFacePath(root, xnode, order_map, node_data);
+	    markPertinentPath(wnode, order_map, node_data, edge_lists,
+			      embed_edge, merge_roots);
+	    markInternalPath(ipath);
+	    int xlp = markExternalPath(xnode, order_map, child_lists, 
+				       pred_map, ancestor_map, low_map);
+	    int zlp = markExternalPath(znode, order_map, child_lists, 
+				       pred_map, ancestor_map, low_map);
+	    markPredPath(root, order_list[xlp < zlp ? xlp : zlp], pred_map);
+	  } else {
+	    markFacePath(ynode, root, order_map, node_data);
+	    markPertinentPath(wnode, order_map, node_data, edge_lists,
+			      embed_edge, merge_roots);
+	    markInternalPath(ipath);
+	    int ylp = markExternalPath(ynode, order_map, child_lists, 
+				       pred_map, ancestor_map, low_map);
+	    int zlp = markExternalPath(znode, order_map, child_lists, 
+				       pred_map, ancestor_map, low_map);
+	    markPredPath(root, order_list[ylp < zlp ? ylp : zlp], pred_map);
+	  }
+	  return;
+	}
+
+	int xlp = markExternalPath(xnode, order_map, child_lists, 
+				   pred_map, ancestor_map, low_map);
+	int ylp = markExternalPath(ynode, order_map, child_lists, 
+				   pred_map, ancestor_map, low_map);
+	int wlp = markExternalPath(wnode, order_map, child_lists, 
+				   pred_map, ancestor_map, low_map);
+
+	if (wlp > xlp && wlp > ylp) {
+	  markFacePath(root, root, order_map, node_data);
+	  markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
+	  return;
+	}
+
+	markInternalPath(ipath);
+	markPertinentPath(wnode, order_map, node_data, edge_lists,
+			  embed_edge, merge_roots);
+
+	if (xlp > ylp && xlp > wlp) {
+	  markFacePath(root, pynode, order_map, node_data);
+	  markFacePath(wnode, xnode, order_map, node_data);
+	  markPredPath(root, order_list[ylp < wlp ? ylp : wlp], pred_map);
+	  return;
+	}
+
+	if (ylp > xlp && ylp > wlp) {
+	  markFacePath(pxnode, root, order_map, node_data);
+	  markFacePath(ynode, wnode, order_map, node_data);
+	  markPredPath(root, order_list[xlp < wlp ? xlp : wlp], pred_map);
+	  return;
+	}
+
+	if (pynode != ynode) {
+	  markFacePath(pxnode, wnode, order_map, node_data);
+
+	  int minlp = xlp < ylp ? xlp : ylp;
+	  if (wlp < minlp) minlp = wlp;
+
+	  int maxlp = xlp > ylp ? xlp : ylp;
+	  if (wlp > maxlp) maxlp = wlp;
+	  
+	  markPredPath(order_list[maxlp], order_list[minlp], pred_map);
+	  return;
+	}
+
+	if (pxnode != xnode) {
+	  markFacePath(wnode, pynode, order_map, node_data);
+
+	  int minlp = xlp < ylp ? xlp : ylp;
+	  if (wlp < minlp) minlp = wlp;
+
+	  int maxlp = xlp > ylp ? xlp : ylp;
+	  if (wlp > maxlp) maxlp = wlp;
+	  
+	  markPredPath(order_list[maxlp], order_list[minlp], pred_map);
+	  return;
+	}
+
+	markFacePath(root, root, order_map, node_data);
+	int minlp = xlp < ylp ? xlp : ylp;
+	if (wlp < minlp) minlp = wlp;
+	markPredPath(root, order_list[minlp], pred_map);
+	return;
+      }
+      
+    }
+
+  };
+
+}
+
+#endif
