lemon/planarity.h
author deba
Tue, 30 Oct 2007 10:51:07 +0000
changeset 2504 46a82ce84cc6
parent 2499 c97596611d59
child 2508 c86db0f7f917
permissions -rw-r--r--
Bug fix
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2007
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 #ifndef LEMON_PLANARITY_H
    19 #define LEMON_PLANARITY_H
    20 
    21 /// \ingroup planar
    22 /// \file
    23 /// \brief Planarity checking, embedding
    24 
    25 #include <vector>
    26 #include <list>
    27 
    28 #include <lemon/dfs.h>
    29 #include <lemon/radix_sort.h>
    30 #include <lemon/maps.h>
    31 #include <lemon/path.h>
    32 #include <lemon/iterable_maps.h>
    33 #include <lemon/edge_set.h>
    34 
    35 
    36 namespace lemon {
    37 
    38   namespace _planarity_bits {
    39 
    40     template <typename UGraph>
    41     struct PlanarityVisitor : DfsVisitor<UGraph> {
    42 
    43       typedef typename UGraph::Node Node;
    44       typedef typename UGraph::Edge Edge;
    45 
    46       typedef typename UGraph::template NodeMap<Edge> PredMap;
    47 
    48       typedef typename UGraph::template UEdgeMap<bool> TreeMap;
    49 
    50       typedef typename UGraph::template NodeMap<int> OrderMap;
    51       typedef std::vector<Node> OrderList;
    52 
    53       typedef typename UGraph::template NodeMap<int> LowMap;
    54       typedef typename UGraph::template NodeMap<int> AncestorMap;
    55 
    56       PlanarityVisitor(const UGraph& ugraph,
    57 		       PredMap& pred_map, TreeMap& tree_map,
    58 		       OrderMap& order_map, OrderList& order_list,
    59 		       AncestorMap& ancestor_map, LowMap& low_map)
    60 	: _ugraph(ugraph), _pred_map(pred_map), _tree_map(tree_map),
    61 	  _order_map(order_map), _order_list(order_list),
    62 	  _ancestor_map(ancestor_map), _low_map(low_map) {}
    63       
    64       void reach(const Node& node) {
    65 	_order_map[node] = _order_list.size();
    66 	_low_map[node] = _order_list.size();
    67 	_ancestor_map[node] = _order_list.size();
    68 	_order_list.push_back(node);
    69       }
    70 
    71       void discover(const Edge& edge) {
    72 	Node source = _ugraph.source(edge);
    73 	Node target = _ugraph.target(edge);
    74 
    75 	_tree_map[edge] = true;
    76 	_pred_map[target] = edge;
    77       }
    78 
    79       void examine(const Edge& edge) {
    80 	Node source = _ugraph.source(edge);
    81 	Node target = _ugraph.target(edge);
    82 	
    83 	if (_order_map[target] < _order_map[source] && !_tree_map[edge]) {
    84 	  if (_low_map[source] > _order_map[target]) {
    85 	    _low_map[source] = _order_map[target];
    86 	  }
    87 	  if (_ancestor_map[source] > _order_map[target]) {
    88 	    _ancestor_map[source] = _order_map[target];
    89 	  }
    90 	}
    91       }
    92 
    93       void backtrack(const Edge& edge) {
    94 	Node source = _ugraph.source(edge);
    95 	Node target = _ugraph.target(edge);
    96 
    97 	if (_low_map[source] > _low_map[target]) {
    98 	  _low_map[source] = _low_map[target];
    99 	}
   100       }
   101 
   102       const UGraph& _ugraph;
   103       PredMap& _pred_map;
   104       TreeMap& _tree_map;
   105       OrderMap& _order_map;
   106       OrderList& _order_list;
   107       AncestorMap& _ancestor_map;
   108       LowMap& _low_map;
   109     };
   110 
   111     template <typename UGraph, bool embedding = true>
   112     struct NodeDataNode {
   113       int prev, next;
   114       int visited;
   115       typename UGraph::Edge first;
   116       bool inverted;
   117     };
   118 
   119     template <typename UGraph>
   120     struct NodeDataNode<UGraph, false> {
   121       int prev, next;
   122       int visited;
   123     };
   124 
   125     template <typename UGraph>
   126     struct ChildListNode {
   127       typedef typename UGraph::Node Node;
   128       Node first;
   129       Node prev, next;
   130     };
   131 
   132     template <typename UGraph>
   133     struct EdgeListNode {
   134       typename UGraph::Edge prev, next;
   135     };
   136 
   137   }
   138 
   139   /// \ingroup planar
   140   ///
   141   /// \brief Planarity checking of an undirected simple graph
   142   ///
   143   /// This class implements the Boyer-Myrvold algorithm for planarity
   144   /// checking of an undirected graph. This class is a simplified
   145   /// version of the PlanarEmbedding algorithm class, and it does
   146   /// provide neither embedding nor kuratowski subdivisons.
   147   template <typename UGraph>
   148   class PlanarityChecking {
   149   private:
   150     
   151     UGRAPH_TYPEDEFS(typename UGraph);
   152       
   153     const UGraph& _ugraph;
   154 
   155   private:
   156 
   157     typedef typename UGraph::template NodeMap<Edge> PredMap;
   158     
   159     typedef typename UGraph::template UEdgeMap<bool> TreeMap;
   160 
   161     typedef typename UGraph::template NodeMap<int> OrderMap;
   162     typedef std::vector<Node> OrderList;
   163 
   164     typedef typename UGraph::template NodeMap<int> LowMap;
   165     typedef typename UGraph::template NodeMap<int> AncestorMap;
   166 
   167     typedef _planarity_bits::NodeDataNode<UGraph> NodeDataNode;
   168     typedef std::vector<NodeDataNode> NodeData;
   169     
   170     typedef _planarity_bits::ChildListNode<UGraph> ChildListNode;
   171     typedef typename UGraph::template NodeMap<ChildListNode> ChildLists;
   172 
   173     typedef typename UGraph::template NodeMap<std::list<int> > MergeRoots;
   174  
   175     typedef typename UGraph::template NodeMap<bool> EmbedEdge;
   176 
   177   public:
   178 
   179     /// \brief Constructor
   180     ///
   181     /// \warining The graph should be simple, i.e. parallel and loop edge
   182     /// free.
   183     PlanarityChecking(const UGraph& ugraph) : _ugraph(ugraph) {}
   184 
   185     /// \brief Runs the algorithm.
   186     ///
   187     /// Runs the algorithm.  
   188     /// \return %True when the graph is planar.
   189     bool run() {
   190       typedef _planarity_bits::PlanarityVisitor<UGraph> Visitor;
   191 
   192       PredMap pred_map(_ugraph, INVALID);
   193       TreeMap tree_map(_ugraph, false);
   194 
   195       OrderMap order_map(_ugraph, -1);
   196       OrderList order_list;
   197 
   198       AncestorMap ancestor_map(_ugraph, -1);
   199       LowMap low_map(_ugraph, -1);
   200 
   201       Visitor visitor(_ugraph, pred_map, tree_map,
   202 		      order_map, order_list, ancestor_map, low_map);
   203       DfsVisit<UGraph, Visitor> visit(_ugraph, visitor);
   204       visit.run();
   205 
   206       ChildLists child_lists(_ugraph);
   207       createChildLists(tree_map, order_map, low_map, child_lists);
   208 
   209       NodeData node_data(2 * order_list.size());
   210       
   211       EmbedEdge embed_edge(_ugraph, false);
   212 
   213       MergeRoots merge_roots(_ugraph);
   214       
   215       for (int i = order_list.size() - 1; i >= 0; --i) {
   216 
   217 	Node node = order_list[i];
   218 
   219 	Node source = node;
   220 	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
   221 	  Node target = _ugraph.target(e);
   222 	  
   223 	  if (order_map[source] < order_map[target] && tree_map[e]) {
   224 	    initFace(target, node_data, order_map, order_list);
   225 	  }
   226 	}
   227 	
   228 	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
   229 	  Node target = _ugraph.target(e);
   230 	  
   231 	  if (order_map[source] < order_map[target] && !tree_map[e]) {
   232 	    embed_edge[target] = true;
   233 	    walkUp(target, source, i, pred_map, low_map,
   234 		   order_map, order_list, node_data, merge_roots);
   235 	  }
   236 	}
   237 
   238 	for (typename MergeRoots::Value::iterator it = 
   239 	       merge_roots[node].begin(); it != merge_roots[node].end(); ++it) {
   240 	  int rn = *it;
   241 	  walkDown(rn, i, node_data, order_list, child_lists, 
   242 		   ancestor_map, low_map, embed_edge, merge_roots);
   243 	}
   244 	merge_roots[node].clear();
   245 	
   246 	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
   247 	  Node target = _ugraph.target(e);
   248 	  
   249 	  if (order_map[source] < order_map[target] && !tree_map[e]) {
   250 	    if (embed_edge[target]) {
   251 	      return false;
   252 	    }
   253 	  }
   254 	}
   255       }
   256 
   257       return true;
   258     }
   259     
   260   private:
   261 
   262     void createChildLists(const TreeMap& tree_map, const OrderMap& order_map,
   263 			  const LowMap& low_map, ChildLists& child_lists) {
   264 
   265       for (NodeIt n(_ugraph); n != INVALID; ++n) {
   266 	Node source = n;
   267 	
   268 	std::vector<Node> targets;  
   269 	for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
   270 	  Node target = _ugraph.target(e);
   271 
   272 	  if (order_map[source] < order_map[target] && tree_map[e]) {
   273 	    targets.push_back(target);
   274 	  }
   275 	}	
   276 
   277 	if (targets.size() == 0) {
   278 	  child_lists[source].first = INVALID;
   279 	} else if (targets.size() == 1) {
   280 	  child_lists[source].first = targets[0];
   281 	  child_lists[targets[0]].prev = INVALID;
   282 	  child_lists[targets[0]].next = INVALID;
   283 	} else {
   284 	  radixSort(targets.begin(), targets.end(), mapFunctor(low_map));
   285 	  for (int i = 1; i < int(targets.size()); ++i) {
   286 	    child_lists[targets[i]].prev = targets[i - 1];
   287 	    child_lists[targets[i - 1]].next = targets[i];
   288 	  }
   289 	  child_lists[targets.back()].next = INVALID; 
   290 	  child_lists[targets.front()].prev = INVALID;
   291 	  child_lists[source].first = targets.front();
   292 	}
   293       }
   294     }
   295 
   296     void walkUp(const Node& node, Node root, int rorder,
   297 		const PredMap& pred_map, const LowMap& low_map,
   298 		const OrderMap& order_map, const OrderList& order_list,
   299 		NodeData& node_data, MergeRoots& merge_roots) {
   300 
   301       int na, nb;
   302       bool da, db;
   303       
   304       na = nb = order_map[node];
   305       da = true; db = false;
   306       
   307       while (true) {
   308 	
   309 	if (node_data[na].visited == rorder) break;
   310 	if (node_data[nb].visited == rorder) break;
   311 
   312 	node_data[na].visited = rorder;
   313 	node_data[nb].visited = rorder;
   314 
   315 	int rn = -1;
   316 
   317 	if (na >= int(order_list.size())) {
   318 	  rn = na;
   319 	} else if (nb >= int(order_list.size())) {
   320 	  rn = nb;
   321 	}
   322 
   323 	if (rn == -1) {
   324 	  int nn;
   325 	  
   326 	  nn = da ? node_data[na].prev : node_data[na].next;
   327 	  da = node_data[nn].prev != na;
   328 	  na = nn;
   329 	  
   330 	  nn = db ? node_data[nb].prev : node_data[nb].next;
   331 	  db = node_data[nn].prev != nb;
   332 	  nb = nn;
   333 
   334 	} else {
   335 
   336 	  Node rep = order_list[rn - order_list.size()];
   337 	  Node parent = _ugraph.source(pred_map[rep]);
   338 
   339 	  if (low_map[rep] < rorder) {
   340 	    merge_roots[parent].push_back(rn);
   341 	  } else {
   342 	    merge_roots[parent].push_front(rn);
   343 	  }
   344 	  
   345 	  if (parent != root) {  
   346 	    na = nb = order_map[parent];
   347 	    da = true; db = false;
   348 	  } else {
   349 	    break;
   350 	  }
   351 	}	
   352       }
   353     }
   354 
   355     void walkDown(int rn, int rorder, NodeData& node_data,
   356 		  OrderList& order_list, ChildLists& child_lists,
   357 		  AncestorMap& ancestor_map, LowMap& low_map,
   358 		  EmbedEdge& embed_edge, MergeRoots& merge_roots) {
   359 
   360       std::vector<std::pair<int, bool> > merge_stack;
   361 
   362       for (int di = 0; di < 2; ++di) {
   363 	bool rd = di == 0;
   364 	int pn = rn;
   365 	int n = rd ? node_data[rn].next : node_data[rn].prev;
   366 	
   367 	while (n != rn) {
   368 	  
   369 	  Node node = order_list[n];
   370 	  
   371 	  if (embed_edge[node]) {
   372 
   373 	    // Merging components on the critical path
   374 	    while (!merge_stack.empty()) {
   375 
   376 	      // Component root
   377 	      int cn = merge_stack.back().first;
   378 	      bool cd = merge_stack.back().second;
   379 	      merge_stack.pop_back();
   380 
   381 	      // Parent of component
   382 	      int dn = merge_stack.back().first;
   383 	      bool dd = merge_stack.back().second;
   384 	      merge_stack.pop_back();
   385 
   386 	      Node parent = order_list[dn];
   387 
   388 	      // Erasing from merge_roots
   389 	      merge_roots[parent].pop_front();
   390 	      
   391 	      Node child = order_list[cn - order_list.size()];
   392 
   393 	      // Erasing from child_lists
   394 	      if (child_lists[child].prev != INVALID) {
   395 		child_lists[child_lists[child].prev].next =
   396 		  child_lists[child].next;
   397 	      } else {
   398 		child_lists[parent].first = child_lists[child].next;
   399 	      }
   400 	      
   401 	      if (child_lists[child].next != INVALID) {
   402 		child_lists[child_lists[child].next].prev =
   403 		  child_lists[child].prev;
   404 	      }
   405 	      
   406 	      // Merging external faces
   407 	      {
   408 		int en = cn;
   409 		cn = cd ? node_data[cn].prev : node_data[cn].next;
   410 		cd = node_data[cn].next == en;
   411 
   412 	      }
   413 
   414 	      if (cd) node_data[cn].next = dn; else node_data[cn].prev = dn;
   415 	      if (dd) node_data[dn].prev = cn; else node_data[dn].next = cn;
   416 
   417 	    }
   418 
   419 	    bool d = pn == node_data[n].prev;
   420 
   421 	    if (node_data[n].prev == node_data[n].next && 
   422 		node_data[n].inverted) {
   423 	      d = !d;
   424 	    }
   425 
   426 	    // Embedding edge into external face
   427 	    if (rd) node_data[rn].next = n; else node_data[rn].prev = n;
   428 	    if (d) node_data[n].prev = rn; else node_data[n].next = rn;
   429 	    pn = rn;
   430 
   431 	    embed_edge[order_list[n]] = false;
   432 	  }
   433 
   434 	  if (!merge_roots[node].empty()) {
   435 
   436 	    bool d = pn == node_data[n].prev;
   437 
   438 	    merge_stack.push_back(std::make_pair(n, d));
   439 
   440 	    int rn = merge_roots[node].front();
   441 	    
   442 	    int xn = node_data[rn].next;
   443 	    Node xnode = order_list[xn];
   444 	    
   445 	    int yn = node_data[rn].prev;
   446 	    Node ynode = order_list[yn];
   447 	    
   448 	    bool rd;
   449 	    if (!external(xnode, rorder, child_lists, ancestor_map, low_map)) {
   450 	      rd = true;
   451 	    } else if (!external(ynode, rorder, child_lists, 
   452 				 ancestor_map, low_map)) {
   453 	      rd = false;
   454 	    } else if (pertinent(xnode, embed_edge, merge_roots)) {
   455 	      rd = true;
   456 	    } else {
   457 	      rd = false;
   458 	    }
   459 	    
   460 	    merge_stack.push_back(std::make_pair(rn, rd));
   461 	    
   462 	    pn = rn;
   463 	    n = rd ? xn : yn;	      
   464 	    	    
   465 	  } else if (!external(node, rorder, child_lists,
   466 			       ancestor_map, low_map)) {
   467 	    int nn = (node_data[n].next != pn ? 
   468 		      node_data[n].next : node_data[n].prev);
   469 
   470 	    bool nd = n == node_data[nn].prev;
   471 
   472 	    if (nd) node_data[nn].prev = pn;
   473 	    else node_data[nn].next = pn; 
   474 
   475 	    if (n == node_data[pn].prev) node_data[pn].prev = nn;
   476 	    else node_data[pn].next = nn;
   477 
   478 	    node_data[nn].inverted = 
   479 	      (node_data[nn].prev == node_data[nn].next && nd != rd);
   480 
   481 	    n = nn;
   482 	  }
   483 	  else break;
   484 	  
   485 	}
   486 
   487 	if (!merge_stack.empty() || n == rn) {
   488 	  break;
   489 	}
   490       }
   491     }
   492 
   493     void initFace(const Node& node, NodeData& node_data, 
   494 		  const OrderMap& order_map, const OrderList& order_list) {
   495       int n = order_map[node];
   496       int rn = n + order_list.size();
   497       
   498       node_data[n].next = node_data[n].prev = rn;
   499       node_data[rn].next = node_data[rn].prev = n;
   500       
   501       node_data[n].visited = order_list.size();
   502       node_data[rn].visited = order_list.size();
   503       
   504     }
   505 
   506     bool external(const Node& node, int rorder,
   507 		  ChildLists& child_lists, AncestorMap& ancestor_map, 
   508 		  LowMap& low_map) {
   509       Node child = child_lists[node].first;
   510 
   511       if (child != INVALID) {
   512 	if (low_map[child] < rorder) return true;
   513       }
   514 
   515       if (ancestor_map[node] < rorder) return true;
   516 
   517       return false;
   518     }
   519 
   520     bool pertinent(const Node& node, const EmbedEdge& embed_edge,
   521 		   const MergeRoots& merge_roots) {
   522       return !merge_roots[node].empty() || embed_edge[node];
   523     }
   524 
   525   };
   526 
   527   /// \ingroup planar
   528   ///
   529   /// \brief Planar embedding of an undirected simple graph
   530   ///
   531   /// This class implements the Boyer-Myrvold algorithm for planar
   532   /// embedding of an undirected graph. The planar embeding is an
   533   /// ordering of the outgoing edges in each node, which is a possible
   534   /// configuration to draw the graph in the plane. If there is not
   535   /// such ordering then the graph contains a \f$ K_5 \f$ (full graph
   536   /// with 5 nodes) or an \f$ K_{3,3} \f$ (complete bipartite graph on
   537   /// 3 ANode and 3 BNode) subdivision.
   538   ///
   539   /// The current implementation calculates an embedding or an
   540   /// Kuratowski subdivision if the graph is not planar. The running
   541   /// time of the algorithm is \f$ O(n) \f$.
   542   template <typename UGraph>
   543   class PlanarEmbedding {
   544   private:
   545     
   546     UGRAPH_TYPEDEFS(typename UGraph);
   547       
   548     const UGraph& _ugraph;
   549     typename UGraph::template EdgeMap<Edge> _embedding;
   550 
   551     typename UGraph::template UEdgeMap<bool> _kuratowski;
   552 
   553   private:
   554 
   555     typedef typename UGraph::template NodeMap<Edge> PredMap;
   556     
   557     typedef typename UGraph::template UEdgeMap<bool> TreeMap;
   558 
   559     typedef typename UGraph::template NodeMap<int> OrderMap;
   560     typedef std::vector<Node> OrderList;
   561 
   562     typedef typename UGraph::template NodeMap<int> LowMap;
   563     typedef typename UGraph::template NodeMap<int> AncestorMap;
   564 
   565     typedef _planarity_bits::NodeDataNode<UGraph> NodeDataNode;
   566     typedef std::vector<NodeDataNode> NodeData;
   567     
   568     typedef _planarity_bits::ChildListNode<UGraph> ChildListNode;
   569     typedef typename UGraph::template NodeMap<ChildListNode> ChildLists;
   570 
   571     typedef typename UGraph::template NodeMap<std::list<int> > MergeRoots;
   572  
   573     typedef typename UGraph::template NodeMap<Edge> EmbedEdge;
   574 
   575     typedef _planarity_bits::EdgeListNode<UGraph> EdgeListNode;
   576     typedef typename UGraph::template EdgeMap<EdgeListNode> EdgeLists;
   577 
   578     typedef typename UGraph::template NodeMap<bool> FlipMap;
   579 
   580     typedef typename UGraph::template NodeMap<int> TypeMap;
   581 
   582     enum IsolatorNodeType {
   583       HIGHX = 6, LOWX = 7,
   584       HIGHY = 8, LOWY = 9,
   585       ROOT = 10, PERTINENT = 11,
   586       INTERNAL = 12
   587     };
   588 
   589   public:
   590 
   591     /// \brief The map for store of embedding
   592     typedef typename UGraph::template EdgeMap<Edge> EmbeddingMap;
   593 
   594     /// \brief Constructor
   595     ///
   596     /// \warining The graph should be simple, i.e. parallel and loop edge
   597     /// free.
   598     PlanarEmbedding(const UGraph& ugraph)
   599       : _ugraph(ugraph), _embedding(_ugraph), _kuratowski(ugraph, false) {}
   600 
   601     /// \brief Runs the algorithm.
   602     ///
   603     /// Runs the algorithm.  
   604     /// \param kuratowski If the parameter is false, then the
   605     /// algorithm does not calculate the isolate Kuratowski
   606     /// subdivisions.
   607     ///\return %True when the graph is planar.
   608     bool run(bool kuratowski = true) {
   609       typedef _planarity_bits::PlanarityVisitor<UGraph> Visitor;
   610 
   611       PredMap pred_map(_ugraph, INVALID);
   612       TreeMap tree_map(_ugraph, false);
   613 
   614       OrderMap order_map(_ugraph, -1);
   615       OrderList order_list;
   616 
   617       AncestorMap ancestor_map(_ugraph, -1);
   618       LowMap low_map(_ugraph, -1);
   619 
   620       Visitor visitor(_ugraph, pred_map, tree_map,
   621 		      order_map, order_list, ancestor_map, low_map);
   622       DfsVisit<UGraph, Visitor> visit(_ugraph, visitor);
   623       visit.run();
   624 
   625       ChildLists child_lists(_ugraph);
   626       createChildLists(tree_map, order_map, low_map, child_lists);
   627 
   628       NodeData node_data(2 * order_list.size());
   629       
   630       EmbedEdge embed_edge(_ugraph, INVALID);
   631 
   632       MergeRoots merge_roots(_ugraph);
   633 
   634       EdgeLists edge_lists(_ugraph);
   635 
   636       FlipMap flip_map(_ugraph, false);
   637 
   638       for (int i = order_list.size() - 1; i >= 0; --i) {
   639 
   640 	Node node = order_list[i];
   641 
   642 	node_data[i].first = INVALID;
   643 	
   644 	Node source = node;
   645 	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
   646 	  Node target = _ugraph.target(e);
   647 	  
   648 	  if (order_map[source] < order_map[target] && tree_map[e]) {
   649 	    initFace(target, edge_lists, node_data,
   650 		      pred_map, order_map, order_list);
   651 	  }
   652 	}
   653 	
   654 	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
   655 	  Node target = _ugraph.target(e);
   656 	  
   657 	  if (order_map[source] < order_map[target] && !tree_map[e]) {
   658 	    embed_edge[target] = e;
   659 	    walkUp(target, source, i, pred_map, low_map,
   660 		   order_map, order_list, node_data, merge_roots);
   661 	  }
   662 	}
   663 
   664 	for (typename MergeRoots::Value::iterator it = 
   665 	       merge_roots[node].begin(); it != merge_roots[node].end(); ++it) {
   666 	  int rn = *it;
   667 	  walkDown(rn, i, node_data, edge_lists, flip_map, order_list, 
   668 		   child_lists, ancestor_map, low_map, embed_edge, merge_roots);
   669 	}
   670 	merge_roots[node].clear();
   671 	
   672 	for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
   673 	  Node target = _ugraph.target(e);
   674 	  
   675 	  if (order_map[source] < order_map[target] && !tree_map[e]) {
   676 	    if (embed_edge[target] != INVALID) {
   677 	      if (kuratowski) {
   678 		isolateKuratowski(e, node_data, edge_lists, flip_map,
   679 				  order_map, order_list, pred_map, child_lists,
   680 				  ancestor_map, low_map, 
   681 				  embed_edge, merge_roots);
   682 	      }
   683 	      return false;
   684 	    }
   685 	  }
   686 	}
   687       }
   688 
   689       for (int i = 0; i < int(order_list.size()); ++i) {
   690 
   691 	mergeRemainingFaces(order_list[i], node_data, order_list, order_map,
   692 			    child_lists, edge_lists);
   693 	storeEmbedding(order_list[i], node_data, order_map, pred_map,
   694 		       edge_lists, flip_map);
   695       }
   696 
   697       return true;
   698     }
   699 
   700     /// \brief Gives back the successor of an edge
   701     ///
   702     /// Gives back the successor of an edge. This function makes
   703     /// possible to query the cyclic order of the outgoing edges from
   704     /// a node.
   705     Edge next(const Edge& edge) const {
   706       return _embedding[edge];
   707     }
   708 
   709     /// \brief Gives back the calculated embedding map
   710     ///
   711     /// The returned map contains the successor of each edge in the
   712     /// graph.
   713     const EmbeddingMap& embedding() const {
   714       return _embedding;
   715     }
   716 
   717     /// \brief Gives back true when the undirected edge is in the
   718     /// kuratowski subdivision
   719     ///
   720     /// Gives back true when the undirected edge is in the kuratowski
   721     /// subdivision
   722     bool kuratowski(const UEdge& uedge) {
   723       return _kuratowski[uedge];
   724     }
   725 
   726   private:
   727 
   728     void createChildLists(const TreeMap& tree_map, const OrderMap& order_map,
   729 			  const LowMap& low_map, ChildLists& child_lists) {
   730 
   731       for (NodeIt n(_ugraph); n != INVALID; ++n) {
   732 	Node source = n;
   733 	
   734 	std::vector<Node> targets;  
   735 	for (OutEdgeIt e(_ugraph, n); e != INVALID; ++e) {
   736 	  Node target = _ugraph.target(e);
   737 
   738 	  if (order_map[source] < order_map[target] && tree_map[e]) {
   739 	    targets.push_back(target);
   740 	  }
   741 	}	
   742 
   743 	if (targets.size() == 0) {
   744 	  child_lists[source].first = INVALID;
   745 	} else if (targets.size() == 1) {
   746 	  child_lists[source].first = targets[0];
   747 	  child_lists[targets[0]].prev = INVALID;
   748 	  child_lists[targets[0]].next = INVALID;
   749 	} else {
   750 	  radixSort(targets.begin(), targets.end(), mapFunctor(low_map));
   751 	  for (int i = 1; i < int(targets.size()); ++i) {
   752 	    child_lists[targets[i]].prev = targets[i - 1];
   753 	    child_lists[targets[i - 1]].next = targets[i];
   754 	  }
   755 	  child_lists[targets.back()].next = INVALID; 
   756 	  child_lists[targets.front()].prev = INVALID;
   757 	  child_lists[source].first = targets.front();
   758 	}
   759       }
   760     }
   761 
   762     void walkUp(const Node& node, Node root, int rorder,
   763 		const PredMap& pred_map, const LowMap& low_map,
   764 		const OrderMap& order_map, const OrderList& order_list,
   765 		NodeData& node_data, MergeRoots& merge_roots) {
   766 
   767       int na, nb;
   768       bool da, db;
   769       
   770       na = nb = order_map[node];
   771       da = true; db = false;
   772       
   773       while (true) {
   774 	
   775 	if (node_data[na].visited == rorder) break;
   776 	if (node_data[nb].visited == rorder) break;
   777 
   778 	node_data[na].visited = rorder;
   779 	node_data[nb].visited = rorder;
   780 
   781 	int rn = -1;
   782 
   783 	if (na >= int(order_list.size())) {
   784 	  rn = na;
   785 	} else if (nb >= int(order_list.size())) {
   786 	  rn = nb;
   787 	}
   788 
   789 	if (rn == -1) {
   790 	  int nn;
   791 	  
   792 	  nn = da ? node_data[na].prev : node_data[na].next;
   793 	  da = node_data[nn].prev != na;
   794 	  na = nn;
   795 	  
   796 	  nn = db ? node_data[nb].prev : node_data[nb].next;
   797 	  db = node_data[nn].prev != nb;
   798 	  nb = nn;
   799 
   800 	} else {
   801 
   802 	  Node rep = order_list[rn - order_list.size()];
   803 	  Node parent = _ugraph.source(pred_map[rep]);
   804 
   805 	  if (low_map[rep] < rorder) {
   806 	    merge_roots[parent].push_back(rn);
   807 	  } else {
   808 	    merge_roots[parent].push_front(rn);
   809 	  }
   810 	  
   811 	  if (parent != root) {  
   812 	    na = nb = order_map[parent];
   813 	    da = true; db = false;
   814 	  } else {
   815 	    break;
   816 	  }
   817 	}	
   818       }
   819     }
   820 
   821     void walkDown(int rn, int rorder, NodeData& node_data,
   822 		  EdgeLists& edge_lists, FlipMap& flip_map, 
   823 		  OrderList& order_list, ChildLists& child_lists,
   824 		  AncestorMap& ancestor_map, LowMap& low_map,
   825 		  EmbedEdge& embed_edge, MergeRoots& merge_roots) {
   826 
   827       std::vector<std::pair<int, bool> > merge_stack;
   828 
   829       for (int di = 0; di < 2; ++di) {
   830 	bool rd = di == 0;
   831 	int pn = rn;
   832 	int n = rd ? node_data[rn].next : node_data[rn].prev;
   833 	
   834 	while (n != rn) {
   835 	  
   836 	  Node node = order_list[n];
   837 	  
   838 	  if (embed_edge[node] != INVALID) {
   839 
   840 	    // Merging components on the critical path
   841 	    while (!merge_stack.empty()) {
   842 
   843 	      // Component root
   844 	      int cn = merge_stack.back().first;
   845 	      bool cd = merge_stack.back().second;
   846 	      merge_stack.pop_back();
   847 
   848 	      // Parent of component
   849 	      int dn = merge_stack.back().first;
   850 	      bool dd = merge_stack.back().second;
   851 	      merge_stack.pop_back();
   852 
   853 	      Node parent = order_list[dn];
   854 
   855 	      // Erasing from merge_roots
   856 	      merge_roots[parent].pop_front();
   857 	      
   858 	      Node child = order_list[cn - order_list.size()];
   859 
   860 	      // Erasing from child_lists
   861 	      if (child_lists[child].prev != INVALID) {
   862 		child_lists[child_lists[child].prev].next =
   863 		  child_lists[child].next;
   864 	      } else {
   865 		child_lists[parent].first = child_lists[child].next;
   866 	      }
   867 	      
   868 	      if (child_lists[child].next != INVALID) {
   869 		child_lists[child_lists[child].next].prev =
   870 		  child_lists[child].prev;
   871 	      }
   872 
   873 	      // Merging edges + flipping
   874 	      Edge de = node_data[dn].first;
   875 	      Edge ce = node_data[cn].first;
   876 
   877 	      flip_map[order_list[cn - order_list.size()]] = cd != dd;
   878 	      if (cd != dd) {
   879 		std::swap(edge_lists[ce].prev, edge_lists[ce].next);
   880 		ce = edge_lists[ce].prev;
   881 		std::swap(edge_lists[ce].prev, edge_lists[ce].next);
   882 	      }
   883 
   884 	      {
   885 		Edge dne = edge_lists[de].next; 
   886 		Edge cne = edge_lists[ce].next; 
   887 
   888 		edge_lists[de].next = cne;
   889 		edge_lists[ce].next = dne;
   890 	      
   891 		edge_lists[dne].prev = ce;
   892 		edge_lists[cne].prev = de;
   893 	      }
   894 	      	      
   895 	      if (dd) {
   896 		node_data[dn].first = ce;
   897 	      }
   898 	      
   899 	      // Merging external faces
   900 	      {
   901 		int en = cn;
   902 		cn = cd ? node_data[cn].prev : node_data[cn].next;
   903 		cd = node_data[cn].next == en;
   904 
   905  		if (node_data[cn].prev == node_data[cn].next && 
   906 		    node_data[cn].inverted) {
   907  		  cd = !cd;
   908  		}
   909 	      }
   910 
   911 	      if (cd) node_data[cn].next = dn; else node_data[cn].prev = dn;
   912 	      if (dd) node_data[dn].prev = cn; else node_data[dn].next = cn;
   913 
   914 	    }
   915 
   916 	    bool d = pn == node_data[n].prev;
   917 
   918 	    if (node_data[n].prev == node_data[n].next && 
   919 		node_data[n].inverted) {
   920 	      d = !d;
   921 	    }
   922 
   923 	    // Add new edge
   924 	    {
   925 	      Edge edge = embed_edge[node];
   926 	      Edge re = node_data[rn].first;
   927 
   928 	      edge_lists[edge_lists[re].next].prev = edge;
   929 	      edge_lists[edge].next = edge_lists[re].next;
   930 	      edge_lists[edge].prev = re;
   931 	      edge_lists[re].next = edge;
   932 
   933 	      if (!rd) {
   934 		node_data[rn].first = edge;
   935 	      }
   936 
   937 	      Edge rev = _ugraph.oppositeEdge(edge);
   938 	      Edge e = node_data[n].first;
   939 
   940 	      edge_lists[edge_lists[e].next].prev = rev;
   941 	      edge_lists[rev].next = edge_lists[e].next;
   942 	      edge_lists[rev].prev = e;
   943 	      edge_lists[e].next = rev;
   944 
   945 	      if (d) {
   946 		node_data[n].first = rev;
   947 	      }
   948 	      
   949 	    }
   950 
   951 	    // Embedding edge into external face
   952 	    if (rd) node_data[rn].next = n; else node_data[rn].prev = n;
   953 	    if (d) node_data[n].prev = rn; else node_data[n].next = rn;
   954 	    pn = rn;
   955 
   956 	    embed_edge[order_list[n]] = INVALID;
   957 	  }
   958 
   959 	  if (!merge_roots[node].empty()) {
   960 
   961 	    bool d = pn == node_data[n].prev;
   962 	    if (node_data[n].prev == node_data[n].next && 
   963 		node_data[n].inverted) {
   964 	      d = !d;
   965 	    }
   966 
   967 	    merge_stack.push_back(std::make_pair(n, d));
   968 
   969 	    int rn = merge_roots[node].front();
   970 	    
   971 	    int xn = node_data[rn].next;
   972 	    Node xnode = order_list[xn];
   973 	    
   974 	    int yn = node_data[rn].prev;
   975 	    Node ynode = order_list[yn];
   976 	    
   977 	    bool rd;
   978 	    if (!external(xnode, rorder, child_lists, ancestor_map, low_map)) {
   979 	      rd = true;
   980 	    } else if (!external(ynode, rorder, child_lists, 
   981 				 ancestor_map, low_map)) {
   982 	      rd = false;
   983 	    } else if (pertinent(xnode, embed_edge, merge_roots)) {
   984 	      rd = true;
   985 	    } else {
   986 	      rd = false;
   987 	    }
   988 	    
   989 	    merge_stack.push_back(std::make_pair(rn, rd));
   990 	    
   991 	    pn = rn;
   992 	    n = rd ? xn : yn;	      
   993 	    	    
   994 	  } else if (!external(node, rorder, child_lists,
   995 			       ancestor_map, low_map)) {
   996 	    int nn = (node_data[n].next != pn ? 
   997 		      node_data[n].next : node_data[n].prev);
   998 
   999 	    bool nd = n == node_data[nn].prev;
  1000 
  1001 	    if (nd) node_data[nn].prev = pn;
  1002 	    else node_data[nn].next = pn; 
  1003 
  1004 	    if (n == node_data[pn].prev) node_data[pn].prev = nn;
  1005 	    else node_data[pn].next = nn;
  1006 
  1007 	    node_data[nn].inverted = 
  1008 	      (node_data[nn].prev == node_data[nn].next && nd != rd);
  1009 
  1010 	    n = nn;
  1011 	  }
  1012 	  else break;
  1013 	  
  1014 	}
  1015 
  1016 	if (!merge_stack.empty() || n == rn) {
  1017 	  break;
  1018 	}
  1019       }
  1020     }
  1021 
  1022     void initFace(const Node& node, EdgeLists& edge_lists,
  1023 		   NodeData& node_data, const PredMap& pred_map,
  1024 		   const OrderMap& order_map, const OrderList& order_list) {
  1025       int n = order_map[node];
  1026       int rn = n + order_list.size();
  1027       
  1028       node_data[n].next = node_data[n].prev = rn;
  1029       node_data[rn].next = node_data[rn].prev = n;
  1030       
  1031       node_data[n].visited = order_list.size();
  1032       node_data[rn].visited = order_list.size();
  1033 
  1034       node_data[n].inverted = false;
  1035       node_data[rn].inverted = false;
  1036 
  1037       Edge edge = pred_map[node];
  1038       Edge rev = _ugraph.oppositeEdge(edge);
  1039 
  1040       node_data[rn].first = edge;
  1041       node_data[n].first = rev;
  1042 
  1043       edge_lists[edge].prev = edge;
  1044       edge_lists[edge].next = edge;
  1045 
  1046       edge_lists[rev].prev = rev;
  1047       edge_lists[rev].next = rev;
  1048 
  1049     }
  1050 
  1051     void mergeRemainingFaces(const Node& node, NodeData& node_data,
  1052 			     OrderList& order_list, OrderMap& order_map,
  1053 			     ChildLists& child_lists, EdgeLists& edge_lists) {
  1054       while (child_lists[node].first != INVALID) {
  1055 	int dd = order_map[node];
  1056 	Node child = child_lists[node].first; 
  1057 	int cd = order_map[child] + order_list.size();
  1058 	child_lists[node].first = child_lists[child].next;
  1059 
  1060 	Edge de = node_data[dd].first;
  1061 	Edge ce = node_data[cd].first;
  1062 
  1063 	if (de != INVALID) {
  1064 	  Edge dne = edge_lists[de].next; 
  1065 	  Edge cne = edge_lists[ce].next; 
  1066 	  
  1067 	  edge_lists[de].next = cne;
  1068 	  edge_lists[ce].next = dne;
  1069 	  
  1070 	  edge_lists[dne].prev = ce;
  1071 	  edge_lists[cne].prev = de;
  1072 	}
  1073 	
  1074 	node_data[dd].first = ce;
  1075 
  1076       }
  1077     }
  1078 
  1079     void storeEmbedding(const Node& node, NodeData& node_data,
  1080 			OrderMap& order_map, PredMap& pred_map,
  1081 			EdgeLists& edge_lists, FlipMap& flip_map) {
  1082 
  1083       if (node_data[order_map[node]].first == INVALID) return;
  1084 
  1085       if (pred_map[node] != INVALID) {
  1086 	Node source = _ugraph.source(pred_map[node]);
  1087 	flip_map[node] = flip_map[node] != flip_map[source];
  1088       }
  1089       
  1090       Edge first = node_data[order_map[node]].first;
  1091       Edge prev = first;
  1092 
  1093       Edge edge = flip_map[node] ?
  1094 	edge_lists[prev].prev : edge_lists[prev].next;
  1095 
  1096       _embedding[prev] = edge;
  1097       
  1098       while (edge != first) {
  1099 	Edge next = edge_lists[edge].prev == prev ?
  1100 	  edge_lists[edge].next : edge_lists[edge].prev;
  1101 	prev = edge; edge = next;
  1102 	_embedding[prev] = edge;
  1103       }
  1104     }
  1105 
  1106 
  1107     bool external(const Node& node, int rorder,
  1108 		  ChildLists& child_lists, AncestorMap& ancestor_map, 
  1109 		  LowMap& low_map) {
  1110       Node child = child_lists[node].first;
  1111 
  1112       if (child != INVALID) {
  1113 	if (low_map[child] < rorder) return true;
  1114       }
  1115 
  1116       if (ancestor_map[node] < rorder) return true;
  1117 
  1118       return false;
  1119     }
  1120 
  1121     bool pertinent(const Node& node, const EmbedEdge& embed_edge,
  1122 		   const MergeRoots& merge_roots) {
  1123       return !merge_roots[node].empty() || embed_edge[node] != INVALID;
  1124     }
  1125 
  1126     int lowPoint(const Node& node, OrderMap& order_map, ChildLists& child_lists,
  1127 		 AncestorMap& ancestor_map, LowMap& low_map) {
  1128       int low_point;
  1129 
  1130       Node child = child_lists[node].first;
  1131 
  1132       if (child != INVALID) {
  1133 	low_point = low_map[child];
  1134       } else {
  1135 	low_point = order_map[node];
  1136       }
  1137 
  1138       if (low_point > ancestor_map[node]) {
  1139 	low_point = ancestor_map[node];
  1140       }
  1141 
  1142       return low_point;
  1143     }
  1144 
  1145     int findComponentRoot(Node root, Node node, ChildLists& child_lists, 
  1146 			  OrderMap& order_map, OrderList& order_list) {
  1147 
  1148       int order = order_map[root];
  1149       int norder = order_map[node];
  1150 
  1151       Node child = child_lists[root].first;
  1152       while (child != INVALID) {
  1153 	int corder = order_map[child];
  1154 	if (corder > order && corder < norder) {
  1155 	  order = corder;
  1156 	}
  1157 	child = child_lists[child].next;
  1158       }
  1159       return order + order_list.size();
  1160     }
  1161 
  1162     Node findPertinent(Node node, OrderMap& order_map, NodeData& node_data,
  1163 		       EmbedEdge& embed_edge, MergeRoots& merge_roots) {
  1164       Node wnode =_ugraph.target(node_data[order_map[node]].first);
  1165       while (!pertinent(wnode, embed_edge, merge_roots)) {
  1166 	wnode = _ugraph.target(node_data[order_map[wnode]].first);
  1167       }
  1168       return wnode;
  1169     }
  1170 
  1171 
  1172     Node findExternal(Node node, int rorder, OrderMap& order_map, 
  1173 		      ChildLists& child_lists, AncestorMap& ancestor_map,
  1174 		      LowMap& low_map, NodeData& node_data) {
  1175       Node wnode =_ugraph.target(node_data[order_map[node]].first);
  1176       while (!external(wnode, rorder, child_lists, ancestor_map, low_map)) {
  1177 	wnode = _ugraph.target(node_data[order_map[wnode]].first);
  1178       }
  1179       return wnode;
  1180     }
  1181 
  1182     void markCommonPath(Node node, int rorder, Node& wnode, Node& znode, 
  1183 			OrderList& order_list, OrderMap& order_map, 
  1184 			NodeData& node_data, EdgeLists& edge_lists, 
  1185 			EmbedEdge& embed_edge, MergeRoots& merge_roots, 
  1186 			ChildLists& child_lists, AncestorMap& ancestor_map, 
  1187 			LowMap& low_map) {
  1188       
  1189       Node cnode = node;
  1190       Node pred = INVALID;
  1191       
  1192       while (true) {
  1193 
  1194 	bool pert = pertinent(cnode, embed_edge, merge_roots);
  1195 	bool ext = external(cnode, rorder, child_lists, ancestor_map, low_map);
  1196 
  1197 	if (pert && ext) {
  1198 	  if (!merge_roots[cnode].empty()) {
  1199 	    int cn = merge_roots[cnode].back();
  1200 	    
  1201 	    if (low_map[order_list[cn - order_list.size()]] < rorder) {
  1202 	      Edge edge = node_data[cn].first;
  1203 	      _kuratowski.set(edge, true);
  1204 	      
  1205 	      pred = cnode;
  1206 	      cnode = _ugraph.target(edge);
  1207 	    
  1208 	      continue;
  1209 	    }
  1210 	  }
  1211 	  wnode = znode = cnode;
  1212 	  return;
  1213 
  1214 	} else if (pert) {
  1215 	  wnode = cnode;
  1216 	  
  1217 	  while (!external(cnode, rorder, child_lists, ancestor_map, low_map)) {
  1218 	    Edge edge = node_data[order_map[cnode]].first;
  1219 	  
  1220 	    if (_ugraph.target(edge) == pred) {
  1221 	      edge = edge_lists[edge].next;
  1222 	    }
  1223 	    _kuratowski.set(edge, true);
  1224 	    
  1225 	    Node next = _ugraph.target(edge);
  1226 	    pred = cnode; cnode = next;
  1227 	  }
  1228 	  
  1229 	  znode = cnode;
  1230 	  return;
  1231 
  1232 	} else if (ext) {
  1233 	  znode = cnode;
  1234 	  
  1235 	  while (!pertinent(cnode, embed_edge, merge_roots)) {
  1236 	    Edge edge = node_data[order_map[cnode]].first;
  1237 	  
  1238 	    if (_ugraph.target(edge) == pred) {
  1239 	      edge = edge_lists[edge].next;
  1240 	    }
  1241 	    _kuratowski.set(edge, true);
  1242 	    
  1243 	    Node next = _ugraph.target(edge);
  1244 	    pred = cnode; cnode = next;
  1245 	  }
  1246 	  
  1247 	  wnode = cnode;
  1248 	  return;
  1249 	  
  1250 	} else {
  1251 	  Edge edge = node_data[order_map[cnode]].first;
  1252 	  
  1253 	  if (_ugraph.target(edge) == pred) {
  1254 	    edge = edge_lists[edge].next;
  1255 	  }
  1256 	  _kuratowski.set(edge, true);
  1257 
  1258 	  Node next = _ugraph.target(edge);
  1259 	  pred = cnode; cnode = next;
  1260 	}
  1261 	
  1262       }
  1263 
  1264     }
  1265 
  1266     void orientComponent(Node root, int rn, OrderMap& order_map,
  1267 			 PredMap& pred_map, NodeData& node_data, 
  1268 			 EdgeLists& edge_lists, FlipMap& flip_map, 
  1269 			 TypeMap& type_map) {
  1270       node_data[order_map[root]].first = node_data[rn].first;
  1271       type_map[root] = 1;
  1272 
  1273       std::vector<Node> st, qu;
  1274 
  1275       st.push_back(root);
  1276       while (!st.empty()) {
  1277 	Node node = st.back();
  1278 	st.pop_back();
  1279 	qu.push_back(node);
  1280 	
  1281 	Edge edge = node_data[order_map[node]].first;
  1282 	
  1283 	if (type_map[_ugraph.target(edge)] == 0) {
  1284 	  st.push_back(_ugraph.target(edge));
  1285 	  type_map[_ugraph.target(edge)] = 1;
  1286 	} 
  1287 	
  1288 	Edge last = edge, pred = edge;
  1289 	edge = edge_lists[edge].next;
  1290 	while (edge != last) {
  1291 
  1292 	  if (type_map[_ugraph.target(edge)] == 0) {
  1293 	    st.push_back(_ugraph.target(edge));
  1294 	    type_map[_ugraph.target(edge)] = 1;
  1295 	  } 
  1296 	  
  1297 	  Edge next = edge_lists[edge].next != pred ? 
  1298 	    edge_lists[edge].next : edge_lists[edge].prev;
  1299 	  pred = edge; edge = next;
  1300 	}
  1301 
  1302       }
  1303 
  1304       type_map[root] = 2;
  1305       flip_map[root] = false;
  1306 
  1307       for (int i = 1; i < int(qu.size()); ++i) {
  1308 
  1309 	Node node = qu[i];
  1310 
  1311 	while (type_map[node] != 2) {
  1312 	  st.push_back(node);
  1313 	  type_map[node] = 2;
  1314 	  node = _ugraph.source(pred_map[node]);
  1315 	}
  1316 
  1317 	bool flip = flip_map[node];
  1318 
  1319 	while (!st.empty()) {
  1320 	  node = st.back();
  1321 	  st.pop_back();
  1322 	  
  1323 	  flip_map[node] = flip != flip_map[node];
  1324 	  flip = flip_map[node];
  1325 
  1326 	  if (flip) {
  1327 	    Edge edge = node_data[order_map[node]].first;
  1328 	    std::swap(edge_lists[edge].prev, edge_lists[edge].next);
  1329 	    edge = edge_lists[edge].prev;
  1330 	    std::swap(edge_lists[edge].prev, edge_lists[edge].next);
  1331 	    node_data[order_map[node]].first = edge;
  1332 	  }
  1333 	}
  1334       }
  1335 
  1336       for (int i = 0; i < int(qu.size()); ++i) {
  1337 
  1338 	Edge edge = node_data[order_map[qu[i]]].first;
  1339 	Edge last = edge, pred = edge;
  1340 
  1341 	edge = edge_lists[edge].next;
  1342 	while (edge != last) {
  1343 
  1344 	  if (edge_lists[edge].next == pred) {
  1345 	    std::swap(edge_lists[edge].next, edge_lists[edge].prev);
  1346 	  } 
  1347 	  pred = edge; edge = edge_lists[edge].next;
  1348 	}
  1349 	
  1350       }
  1351     }
  1352 
  1353     void setFaceFlags(Node root, Node wnode, Node ynode, Node xnode,
  1354 		      OrderMap& order_map, NodeData& node_data, 
  1355 		      TypeMap& type_map) {
  1356       Node node = _ugraph.target(node_data[order_map[root]].first);
  1357 
  1358       while (node != ynode) {
  1359 	type_map[node] = HIGHY;
  1360 	node = _ugraph.target(node_data[order_map[node]].first);
  1361       }
  1362 
  1363       while (node != wnode) {
  1364 	type_map[node] = LOWY;
  1365 	node = _ugraph.target(node_data[order_map[node]].first);
  1366       }
  1367       
  1368       node = _ugraph.target(node_data[order_map[wnode]].first);
  1369       
  1370       while (node != xnode) {
  1371 	type_map[node] = LOWX;
  1372 	node = _ugraph.target(node_data[order_map[node]].first);
  1373       }
  1374       type_map[node] = LOWX;
  1375 
  1376       node = _ugraph.target(node_data[order_map[xnode]].first);
  1377       while (node != root) {
  1378 	type_map[node] = HIGHX;
  1379 	node = _ugraph.target(node_data[order_map[node]].first);
  1380       }
  1381 
  1382       type_map[wnode] = PERTINENT;
  1383       type_map[root] = ROOT;
  1384     }
  1385 
  1386     void findInternalPath(std::vector<Edge>& ipath,
  1387 			  Node wnode, Node root, TypeMap& type_map, 
  1388 			  OrderMap& order_map, NodeData& node_data, 
  1389 			  EdgeLists& edge_lists) {
  1390       std::vector<Edge> st;
  1391 
  1392       Node node = wnode;
  1393       
  1394       while (node != root) {
  1395 	Edge edge = edge_lists[node_data[order_map[node]].first].next;
  1396 	st.push_back(edge);
  1397 	node = _ugraph.target(edge);
  1398       }
  1399       
  1400       while (true) {
  1401 	Edge edge = st.back();
  1402 	if (type_map[_ugraph.target(edge)] == LOWX ||
  1403 	    type_map[_ugraph.target(edge)] == HIGHX) {
  1404 	  break;
  1405 	}
  1406 	if (type_map[_ugraph.target(edge)] == 2) {
  1407 	  type_map[_ugraph.target(edge)] = 3;
  1408 	  
  1409 	  edge = edge_lists[_ugraph.oppositeEdge(edge)].next;
  1410 	  st.push_back(edge);
  1411 	} else {
  1412 	  st.pop_back();
  1413 	  edge = edge_lists[edge].next;
  1414 	  
  1415 	  while (_ugraph.oppositeEdge(edge) == st.back()) {
  1416 	    edge = st.back();
  1417 	    st.pop_back();
  1418 	    edge = edge_lists[edge].next;
  1419 	  }
  1420 	  st.push_back(edge);
  1421 	}
  1422       }
  1423       
  1424       for (int i = 0; i < int(st.size()); ++i) {
  1425 	if (type_map[_ugraph.target(st[i])] != LOWY &&
  1426 	    type_map[_ugraph.target(st[i])] != HIGHY) {
  1427 	  for (; i < int(st.size()); ++i) {
  1428 	    ipath.push_back(st[i]);
  1429 	  }
  1430 	}
  1431       }
  1432     }
  1433 
  1434     void setInternalFlags(std::vector<Edge>& ipath, TypeMap& type_map) {
  1435       for (int i = 1; i < int(ipath.size()); ++i) {
  1436 	type_map[_ugraph.source(ipath[i])] = INTERNAL;
  1437       }
  1438     }
  1439 
  1440     void findPilePath(std::vector<Edge>& ppath,
  1441 		      Node root, TypeMap& type_map, OrderMap& order_map, 
  1442 		      NodeData& node_data, EdgeLists& edge_lists) {
  1443       std::vector<Edge> st;
  1444 
  1445       st.push_back(_ugraph.oppositeEdge(node_data[order_map[root]].first));
  1446       st.push_back(node_data[order_map[root]].first);
  1447       
  1448       while (st.size() > 1) {
  1449 	Edge edge = st.back();
  1450 	if (type_map[_ugraph.target(edge)] == INTERNAL) {
  1451 	  break;
  1452 	}
  1453 	if (type_map[_ugraph.target(edge)] == 3) {
  1454 	  type_map[_ugraph.target(edge)] = 4;
  1455 	  
  1456 	  edge = edge_lists[_ugraph.oppositeEdge(edge)].next;
  1457 	  st.push_back(edge);
  1458 	} else {
  1459 	  st.pop_back();
  1460 	  edge = edge_lists[edge].next;
  1461 	  
  1462 	  while (!st.empty() && _ugraph.oppositeEdge(edge) == st.back()) {
  1463 	    edge = st.back();
  1464 	    st.pop_back();
  1465 	    edge = edge_lists[edge].next;
  1466 	  }
  1467 	  st.push_back(edge);
  1468 	}
  1469       }
  1470       
  1471       for (int i = 1; i < int(st.size()); ++i) {
  1472 	ppath.push_back(st[i]);
  1473       }
  1474     }
  1475 
  1476 
  1477     int markExternalPath(Node node, OrderMap& order_map,
  1478 			 ChildLists& child_lists, PredMap& pred_map,
  1479 			 AncestorMap& ancestor_map, LowMap& low_map) {
  1480       int lp = lowPoint(node, order_map, child_lists,
  1481 			ancestor_map, low_map);
  1482       
  1483       if (ancestor_map[node] != lp) {
  1484 	node = child_lists[node].first;
  1485 	_kuratowski[pred_map[node]] = true;
  1486 
  1487 	while (ancestor_map[node] != lp) {
  1488 	  for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
  1489 	    Node tnode = _ugraph.target(e); 
  1490 	    if (order_map[tnode] > order_map[node] && low_map[tnode] == lp) {
  1491 	      node = tnode;
  1492 	      _kuratowski[e] = true;
  1493 	      break;
  1494 	    }
  1495 	  }
  1496 	}
  1497       }
  1498 
  1499       for (OutEdgeIt e(_ugraph, node); e != INVALID; ++e) {
  1500 	if (order_map[_ugraph.target(e)] == lp) {
  1501 	  _kuratowski[e] = true;
  1502 	  break;
  1503 	}
  1504       }
  1505       
  1506       return lp;
  1507     }
  1508 
  1509     void markPertinentPath(Node node, OrderMap& order_map, 
  1510 			   NodeData& node_data, EdgeLists& edge_lists,
  1511 			   EmbedEdge& embed_edge, MergeRoots& merge_roots) {
  1512       while (embed_edge[node] == INVALID) {
  1513 	int n = merge_roots[node].front();
  1514 	Edge edge = node_data[n].first;
  1515 
  1516 	_kuratowski.set(edge, true);
  1517 	
  1518 	Node pred = node;
  1519 	node = _ugraph.target(edge);
  1520 	while (!pertinent(node, embed_edge, merge_roots)) {
  1521 	  edge = node_data[order_map[node]].first;
  1522 	  if (_ugraph.target(edge) == pred) {
  1523 	    edge = edge_lists[edge].next;
  1524 	  }
  1525 	  _kuratowski.set(edge, true);
  1526 	  pred = node;
  1527 	  node = _ugraph.target(edge);
  1528 	}
  1529       }
  1530       _kuratowski.set(embed_edge[node], true);
  1531     } 
  1532 
  1533     void markPredPath(Node node, Node snode, PredMap& pred_map) {
  1534       while (node != snode) {
  1535 	_kuratowski.set(pred_map[node], true);
  1536 	node = _ugraph.source(pred_map[node]);
  1537       }
  1538     }
  1539 
  1540     void markFacePath(Node ynode, Node xnode, 
  1541 		      OrderMap& order_map, NodeData& node_data) {
  1542       Edge edge = node_data[order_map[ynode]].first;
  1543       Node node = _ugraph.target(edge);
  1544       _kuratowski.set(edge, true);
  1545 	
  1546       while (node != xnode) {
  1547 	edge = node_data[order_map[node]].first;
  1548 	_kuratowski.set(edge, true);
  1549 	node = _ugraph.target(edge);
  1550       }
  1551     }
  1552 
  1553     void markInternalPath(std::vector<Edge>& path) {
  1554       for (int i = 0; i < int(path.size()); ++i) {
  1555 	_kuratowski.set(path[i], true);
  1556       }
  1557     }
  1558 
  1559     void markPilePath(std::vector<Edge>& path) {
  1560       for (int i = 0; i < int(path.size()); ++i) {
  1561 	_kuratowski.set(path[i], true);
  1562       }
  1563     }
  1564 
  1565     void isolateKuratowski(Edge edge, NodeData& node_data, 
  1566 			   EdgeLists& edge_lists, FlipMap& flip_map,
  1567 			   OrderMap& order_map, OrderList& order_list, 
  1568 			   PredMap& pred_map, ChildLists& child_lists,
  1569 			   AncestorMap& ancestor_map, LowMap& low_map, 
  1570 			   EmbedEdge& embed_edge, MergeRoots& merge_roots) {
  1571 
  1572       Node root = _ugraph.source(edge);
  1573       Node enode = _ugraph.target(edge);
  1574 
  1575       int rorder = order_map[root];
  1576 
  1577       TypeMap type_map(_ugraph, 0);
  1578 
  1579       int rn = findComponentRoot(root, enode, child_lists, 
  1580 				 order_map, order_list);
  1581 
  1582       Node xnode = order_list[node_data[rn].next];
  1583       Node ynode = order_list[node_data[rn].prev];
  1584 
  1585       // Minor-A
  1586       {
  1587 	while (!merge_roots[xnode].empty() || !merge_roots[ynode].empty()) {
  1588 	  
  1589 	  if (!merge_roots[xnode].empty()) {
  1590 	    root = xnode;
  1591 	    rn = merge_roots[xnode].front();
  1592 	  } else {
  1593 	    root = ynode;
  1594 	    rn = merge_roots[ynode].front();
  1595 	  }
  1596 	  
  1597 	  xnode = order_list[node_data[rn].next];
  1598 	  ynode = order_list[node_data[rn].prev];
  1599 	}
  1600 	
  1601 	if (root != _ugraph.source(edge)) {
  1602 	  orientComponent(root, rn, order_map, pred_map, 
  1603 			  node_data, edge_lists, flip_map, type_map);
  1604 	  markFacePath(root, root, order_map, node_data);
  1605 	  int xlp = markExternalPath(xnode, order_map, child_lists, 
  1606 				     pred_map, ancestor_map, low_map);
  1607 	  int ylp = markExternalPath(ynode, order_map, child_lists, 
  1608 				     pred_map, ancestor_map, low_map);
  1609 	  markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
  1610 	  Node lwnode = findPertinent(ynode, order_map, node_data,
  1611 				     embed_edge, merge_roots);
  1612 	  
  1613 	  markPertinentPath(lwnode, order_map, node_data, edge_lists,
  1614 			    embed_edge, merge_roots);
  1615 	  
  1616 	  return;
  1617 	}
  1618       }
  1619       
  1620       orientComponent(root, rn, order_map, pred_map, 
  1621 		      node_data, edge_lists, flip_map, type_map);
  1622 
  1623       Node wnode = findPertinent(ynode, order_map, node_data,
  1624 				 embed_edge, merge_roots);
  1625       setFaceFlags(root, wnode, ynode, xnode, order_map, node_data, type_map);
  1626 
  1627       
  1628       //Minor-B
  1629       if (!merge_roots[wnode].empty()) {
  1630 	int cn = merge_roots[wnode].back();
  1631 	Node rep = order_list[cn - order_list.size()];
  1632 	if (low_map[rep] < rorder) {
  1633 	  markFacePath(root, root, order_map, node_data);
  1634 	  int xlp = markExternalPath(xnode, order_map, child_lists, 
  1635 				     pred_map, ancestor_map, low_map);
  1636 	  int ylp = markExternalPath(ynode, order_map, child_lists, 
  1637 				     pred_map, ancestor_map, low_map);
  1638 
  1639 	  Node lwnode, lznode;
  1640 	  markCommonPath(wnode, rorder, lwnode, lznode, order_list, 
  1641 			 order_map, node_data, edge_lists, embed_edge, 
  1642 			 merge_roots, child_lists, ancestor_map, low_map);
  1643 	  	  
  1644 	  markPertinentPath(lwnode, order_map, node_data, edge_lists,
  1645 			    embed_edge, merge_roots);
  1646 	  int zlp = markExternalPath(lznode, order_map, child_lists, 
  1647 				     pred_map, ancestor_map, low_map);
  1648 
  1649 	  int minlp = xlp < ylp ? xlp : ylp;
  1650 	  if (zlp < minlp) minlp = zlp;
  1651 
  1652 	  int maxlp = xlp > ylp ? xlp : ylp;
  1653 	  if (zlp > maxlp) maxlp = zlp;
  1654 	  
  1655 	  markPredPath(order_list[maxlp], order_list[minlp], pred_map);
  1656 	  
  1657 	  return;
  1658 	}
  1659       }
  1660 
  1661       Node pxnode, pynode;
  1662       std::vector<Edge> ipath;
  1663       findInternalPath(ipath, wnode, root, type_map, order_map,
  1664 		       node_data, edge_lists);
  1665       setInternalFlags(ipath, type_map);
  1666       pynode = _ugraph.source(ipath.front());
  1667       pxnode = _ugraph.target(ipath.back());
  1668 
  1669       wnode = findPertinent(pynode, order_map, node_data,
  1670 			    embed_edge, merge_roots);
  1671       
  1672       // Minor-C
  1673       {
  1674 	if (type_map[_ugraph.source(ipath.front())] == HIGHY) {
  1675 	  if (type_map[_ugraph.target(ipath.back())] == HIGHX) {
  1676 	    markFacePath(xnode, pxnode, order_map, node_data);
  1677 	  }
  1678 	  markFacePath(root, xnode, order_map, node_data);
  1679 	  markPertinentPath(wnode, order_map, node_data, edge_lists,
  1680 			    embed_edge, merge_roots);
  1681 	  markInternalPath(ipath);
  1682 	  int xlp = markExternalPath(xnode, order_map, child_lists, 
  1683 				     pred_map, ancestor_map, low_map);
  1684 	  int ylp = markExternalPath(ynode, order_map, child_lists, 
  1685 				     pred_map, ancestor_map, low_map);
  1686 	  markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
  1687 	  return;
  1688 	}
  1689 
  1690 	if (type_map[_ugraph.target(ipath.back())] == HIGHX) {
  1691 	  markFacePath(ynode, root, order_map, node_data);
  1692 	  markPertinentPath(wnode, order_map, node_data, edge_lists,
  1693 			    embed_edge, merge_roots);
  1694 	  markInternalPath(ipath);
  1695 	  int xlp = markExternalPath(xnode, order_map, child_lists, 
  1696 				     pred_map, ancestor_map, low_map);
  1697 	  int ylp = markExternalPath(ynode, order_map, child_lists, 
  1698 				     pred_map, ancestor_map, low_map);
  1699 	  markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
  1700 	  return;
  1701 	}
  1702       }
  1703 
  1704       std::vector<Edge> ppath;
  1705       findPilePath(ppath, root, type_map, order_map, node_data, edge_lists);
  1706       
  1707       // Minor-D
  1708       if (!ppath.empty()) {
  1709 	markFacePath(ynode, xnode, order_map, node_data);
  1710 	markPertinentPath(wnode, order_map, node_data, edge_lists,
  1711 			  embed_edge, merge_roots);
  1712 	markPilePath(ppath);
  1713 	markInternalPath(ipath);
  1714 	int xlp = markExternalPath(xnode, order_map, child_lists, 
  1715 				   pred_map, ancestor_map, low_map);
  1716 	int ylp = markExternalPath(ynode, order_map, child_lists, 
  1717 				   pred_map, ancestor_map, low_map);
  1718 	markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
  1719 	return;
  1720       }
  1721 
  1722       // Minor-E*
  1723       {
  1724 
  1725 	if (!external(wnode, rorder, child_lists, ancestor_map, low_map)) {
  1726 	  Node znode = findExternal(pynode, rorder, order_map, 
  1727 				    child_lists, ancestor_map,
  1728 				    low_map, node_data);
  1729 	  
  1730 	  if (type_map[znode] == LOWY) {
  1731 	    markFacePath(root, xnode, order_map, node_data);
  1732 	    markPertinentPath(wnode, order_map, node_data, edge_lists,
  1733 			      embed_edge, merge_roots);
  1734 	    markInternalPath(ipath);
  1735 	    int xlp = markExternalPath(xnode, order_map, child_lists, 
  1736 				       pred_map, ancestor_map, low_map);
  1737 	    int zlp = markExternalPath(znode, order_map, child_lists, 
  1738 				       pred_map, ancestor_map, low_map);
  1739 	    markPredPath(root, order_list[xlp < zlp ? xlp : zlp], pred_map);
  1740 	  } else {
  1741 	    markFacePath(ynode, root, order_map, node_data);
  1742 	    markPertinentPath(wnode, order_map, node_data, edge_lists,
  1743 			      embed_edge, merge_roots);
  1744 	    markInternalPath(ipath);
  1745 	    int ylp = markExternalPath(ynode, order_map, child_lists, 
  1746 				       pred_map, ancestor_map, low_map);
  1747 	    int zlp = markExternalPath(znode, order_map, child_lists, 
  1748 				       pred_map, ancestor_map, low_map);
  1749 	    markPredPath(root, order_list[ylp < zlp ? ylp : zlp], pred_map);
  1750 	  }
  1751 	  return;
  1752 	}
  1753 
  1754 	int xlp = markExternalPath(xnode, order_map, child_lists, 
  1755 				   pred_map, ancestor_map, low_map);
  1756 	int ylp = markExternalPath(ynode, order_map, child_lists, 
  1757 				   pred_map, ancestor_map, low_map);
  1758 	int wlp = markExternalPath(wnode, order_map, child_lists, 
  1759 				   pred_map, ancestor_map, low_map);
  1760 
  1761 	if (wlp > xlp && wlp > ylp) {
  1762 	  markFacePath(root, root, order_map, node_data);
  1763 	  markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);
  1764 	  return;
  1765 	}
  1766 
  1767 	markInternalPath(ipath);
  1768 	markPertinentPath(wnode, order_map, node_data, edge_lists,
  1769 			  embed_edge, merge_roots);
  1770 
  1771 	if (xlp > ylp && xlp > wlp) {
  1772 	  markFacePath(root, pynode, order_map, node_data);
  1773 	  markFacePath(wnode, xnode, order_map, node_data);
  1774 	  markPredPath(root, order_list[ylp < wlp ? ylp : wlp], pred_map);
  1775 	  return;
  1776 	}
  1777 
  1778 	if (ylp > xlp && ylp > wlp) {
  1779 	  markFacePath(pxnode, root, order_map, node_data);
  1780 	  markFacePath(ynode, wnode, order_map, node_data);
  1781 	  markPredPath(root, order_list[xlp < wlp ? xlp : wlp], pred_map);
  1782 	  return;
  1783 	}
  1784 
  1785 	if (pynode != ynode) {
  1786 	  markFacePath(pxnode, wnode, order_map, node_data);
  1787 
  1788 	  int minlp = xlp < ylp ? xlp : ylp;
  1789 	  if (wlp < minlp) minlp = wlp;
  1790 
  1791 	  int maxlp = xlp > ylp ? xlp : ylp;
  1792 	  if (wlp > maxlp) maxlp = wlp;
  1793 	  
  1794 	  markPredPath(order_list[maxlp], order_list[minlp], pred_map);
  1795 	  return;
  1796 	}
  1797 
  1798 	if (pxnode != xnode) {
  1799 	  markFacePath(wnode, pynode, order_map, node_data);
  1800 
  1801 	  int minlp = xlp < ylp ? xlp : ylp;
  1802 	  if (wlp < minlp) minlp = wlp;
  1803 
  1804 	  int maxlp = xlp > ylp ? xlp : ylp;
  1805 	  if (wlp > maxlp) maxlp = wlp;
  1806 	  
  1807 	  markPredPath(order_list[maxlp], order_list[minlp], pred_map);
  1808 	  return;
  1809 	}
  1810 
  1811 	markFacePath(root, root, order_map, node_data);
  1812 	int minlp = xlp < ylp ? xlp : ylp;
  1813 	if (wlp < minlp) minlp = wlp;
  1814 	markPredPath(root, order_list[minlp], pred_map);
  1815 	return;
  1816       }
  1817       
  1818     }
  1819 
  1820   };
  1821 
  1822   namespace _planarity_bits {
  1823 
  1824     template <typename UGraph, typename EmbeddingMap>
  1825     void makeConnected(UGraph& ugraph, EmbeddingMap& embedding) {
  1826       DfsVisitor<UGraph> null_visitor;
  1827       DfsVisit<UGraph, DfsVisitor<UGraph> > dfs(ugraph, null_visitor);
  1828       dfs.init();
  1829 
  1830       typename UGraph::Node u = INVALID;
  1831       for (typename UGraph::NodeIt n(ugraph); n != INVALID; ++n) {
  1832 	if (!dfs.reached(n)) {
  1833 	  dfs.addSource(n);
  1834 	  dfs.start();
  1835 	  if (u == INVALID) {
  1836 	    u = n;
  1837 	  } else {
  1838 	    typename UGraph::Node v = n;
  1839 	    
  1840 	    typename UGraph::Edge ue = typename UGraph::OutEdgeIt(ugraph, u);
  1841 	    typename UGraph::Edge ve = typename UGraph::OutEdgeIt(ugraph, v);
  1842 
  1843 	    typename UGraph::Edge e = ugraph.direct(ugraph.addEdge(u, v), true);
  1844 	    
  1845 	    if (ue != INVALID) {
  1846 	      embedding[e] = embedding[ue];
  1847 	      embedding[ue] = e;
  1848 	    } else {
  1849 	      embedding[e] = e;
  1850 	    }
  1851 
  1852 	    if (ve != INVALID) {
  1853 	      embedding[ugraph.oppositeEdge(e)] = embedding[ve];
  1854 	      embedding[ve] = ugraph.oppositeEdge(e);
  1855 	    } else {
  1856 	      embedding[ugraph.oppositeEdge(e)] = ugraph.oppositeEdge(e);
  1857 	    }
  1858 	  }
  1859 	}
  1860       }
  1861     }
  1862 
  1863     template <typename UGraph, typename EmbeddingMap>
  1864     void makeBiNodeConnected(UGraph& ugraph, EmbeddingMap& embedding) {
  1865       typename UGraph::template EdgeMap<bool> processed(ugraph);
  1866 
  1867       std::vector<typename UGraph::Edge> edges;
  1868       for (typename UGraph::EdgeIt e(ugraph); e != INVALID; ++e) {
  1869 	edges.push_back(e);
  1870       }
  1871 
  1872       IterableBoolMap<UGraph, typename UGraph::Node> visited(ugraph, false);
  1873       
  1874       for (int i = 0; i < int(edges.size()); ++i) {
  1875 	typename UGraph::Edge pp = edges[i];
  1876 	if (processed[pp]) continue;
  1877 
  1878 	typename UGraph::Edge e = embedding[ugraph.oppositeEdge(pp)];
  1879 	processed[e] = true;
  1880 	visited.set(ugraph.source(e), true);
  1881 	
  1882 	typename UGraph::Edge p = e, l = e;
  1883 	e = embedding[ugraph.oppositeEdge(e)];
  1884 	
  1885 	while (e != l) {
  1886 	  processed[e] = true;
  1887 
  1888 	  if (visited[ugraph.source(e)]) {
  1889 	    
  1890 	    typename UGraph::Edge n = 
  1891 	      ugraph.direct(ugraph.addEdge(ugraph.source(p), 
  1892 					   ugraph.target(e)), true);
  1893 	    embedding[n] = p;
  1894 	    embedding[ugraph.oppositeEdge(pp)] = n;
  1895 
  1896 	    embedding[ugraph.oppositeEdge(n)] = 
  1897 	      embedding[ugraph.oppositeEdge(e)];
  1898 	    embedding[ugraph.oppositeEdge(e)] = 
  1899 	      ugraph.oppositeEdge(n);
  1900 	    
  1901 	    p = n;
  1902 	    e = embedding[ugraph.oppositeEdge(n)];
  1903 	  } else {
  1904 	    visited.set(ugraph.source(e), true);
  1905 	    pp = p;
  1906 	    p = e;
  1907 	    e = embedding[ugraph.oppositeEdge(e)];
  1908 	  }
  1909 	}
  1910 	visited.setAll(false);
  1911       }
  1912     }
  1913     
  1914     
  1915     template <typename UGraph, typename EmbeddingMap>
  1916     void makeMaxPlanar(UGraph& ugraph, EmbeddingMap& embedding) {
  1917       
  1918       typename UGraph::template NodeMap<int> degree(ugraph);
  1919 
  1920       for (typename UGraph::NodeIt n(ugraph); n != INVALID; ++n) {
  1921 	degree[n] = countIncEdges(ugraph, n);
  1922       }
  1923 
  1924       typename UGraph::template EdgeMap<bool> processed(ugraph);
  1925       IterableBoolMap<UGraph, typename UGraph::Node> visited(ugraph, false);
  1926 
  1927       std::vector<typename UGraph::Edge> edges;
  1928       for (typename UGraph::EdgeIt e(ugraph); e != INVALID; ++e) {
  1929 	edges.push_back(e);
  1930       }
  1931 
  1932       for (int i = 0; i < int(edges.size()); ++i) {
  1933 	typename UGraph::Edge e = edges[i];
  1934 
  1935 	if (processed[e]) continue;
  1936 	processed[e] = true;
  1937 
  1938 	typename UGraph::Edge mine = e;
  1939 	int mind = degree[ugraph.source(e)];
  1940 
  1941 	int face_size = 1;	
  1942 
  1943 	typename UGraph::Edge l = e;
  1944 	e = embedding[ugraph.oppositeEdge(e)];
  1945 	while (l != e) {
  1946 	  processed[e] = true;
  1947 
  1948 	  ++face_size;
  1949 
  1950 	  if (degree[ugraph.source(e)] < mind) {
  1951 	    mine = e;
  1952 	    mind = degree[ugraph.source(e)];
  1953 	  }
  1954 	  
  1955 	  e = embedding[ugraph.oppositeEdge(e)];	  
  1956 	}
  1957 	
  1958 	if (face_size < 4) {
  1959 	  continue;
  1960 	}
  1961 
  1962 	typename UGraph::Node s = ugraph.source(mine);
  1963 	for (typename UGraph::OutEdgeIt e(ugraph, s); e != INVALID; ++e) {
  1964 	  visited.set(ugraph.target(e), true); 
  1965 	}
  1966 
  1967 	typename UGraph::Edge oppe = INVALID;
  1968 
  1969 	e = embedding[ugraph.oppositeEdge(mine)];
  1970 	e = embedding[ugraph.oppositeEdge(e)];
  1971 	while (ugraph.target(e) != s) {
  1972 	  if (visited[ugraph.source(e)]) {
  1973 	    oppe = e;
  1974 	    break;
  1975 	  }
  1976 	  e = embedding[ugraph.oppositeEdge(e)];
  1977 	}
  1978 	visited.setAll(false);
  1979 	
  1980 	if (oppe == INVALID) {
  1981 
  1982 	  e = embedding[ugraph.oppositeEdge(mine)];
  1983 	  typename UGraph::Edge pn = mine, p = e;
  1984 
  1985 	  e = embedding[ugraph.oppositeEdge(e)];
  1986 	  while (ugraph.target(e) != s) {
  1987 	    typename UGraph::Edge n = 
  1988 	      ugraph.direct(ugraph.addEdge(s, ugraph.source(e)), true);
  1989 
  1990 	    embedding[n] = pn;
  1991 	    embedding[ugraph.oppositeEdge(n)] = e;
  1992 	    embedding[ugraph.oppositeEdge(p)] = ugraph.oppositeEdge(n);
  1993 
  1994 	    pn = n;
  1995 	    
  1996 	    p = e;
  1997 	    e = embedding[ugraph.oppositeEdge(e)];
  1998 	  }
  1999 
  2000 	  embedding[ugraph.oppositeEdge(e)] = pn;
  2001 
  2002 	} else {
  2003 
  2004 	  mine = embedding[ugraph.oppositeEdge(mine)];
  2005 	  s = ugraph.source(mine);
  2006 	  oppe = embedding[ugraph.oppositeEdge(oppe)];
  2007 	  typename UGraph::Node t = ugraph.source(oppe);
  2008 	  
  2009 	  typename UGraph::Edge ce = ugraph.direct(ugraph.addEdge(s, t), true);
  2010 	  embedding[ce] = mine;
  2011 	  embedding[ugraph.oppositeEdge(ce)] = oppe;
  2012 
  2013 	  typename UGraph::Edge pn = ce, p = oppe;	  
  2014 	  e = embedding[ugraph.oppositeEdge(oppe)];
  2015 	  while (ugraph.target(e) != s) {
  2016 	    typename UGraph::Edge n = 
  2017 	      ugraph.direct(ugraph.addEdge(s, ugraph.source(e)), true);
  2018 
  2019 	    embedding[n] = pn;
  2020 	    embedding[ugraph.oppositeEdge(n)] = e;
  2021 	    embedding[ugraph.oppositeEdge(p)] = ugraph.oppositeEdge(n);
  2022 
  2023 	    pn = n;
  2024 	    
  2025 	    p = e;
  2026 	    e = embedding[ugraph.oppositeEdge(e)];
  2027 	    
  2028 	  }
  2029 	  embedding[ugraph.oppositeEdge(e)] = pn;
  2030 
  2031 	  pn = ugraph.oppositeEdge(ce), p = mine;	  
  2032 	  e = embedding[ugraph.oppositeEdge(mine)];
  2033 	  while (ugraph.target(e) != t) {
  2034 	    typename UGraph::Edge n = 
  2035 	      ugraph.direct(ugraph.addEdge(t, ugraph.source(e)), true);
  2036 
  2037 	    embedding[n] = pn;
  2038 	    embedding[ugraph.oppositeEdge(n)] = e;
  2039 	    embedding[ugraph.oppositeEdge(p)] = ugraph.oppositeEdge(n);
  2040 
  2041 	    pn = n;
  2042 	    
  2043 	    p = e;
  2044 	    e = embedding[ugraph.oppositeEdge(e)];
  2045 	    
  2046 	  }
  2047 	  embedding[ugraph.oppositeEdge(e)] = pn;
  2048 	}
  2049       }
  2050     }
  2051 
  2052   }
  2053 
  2054   /// \ingroup planar
  2055   ///
  2056   /// \brief Schnyder's planar drawing algorithms
  2057   ///
  2058   /// The planar drawing algorithm calculates location for each node
  2059   /// in the plane, which coordinates satisfies that if each edge is
  2060   /// represented with a straight line then the edges will not
  2061   /// intersect each other.
  2062   ///
  2063   /// Scnyder's algorithm embeds the graph on \c (n-2,n-2) size grid,
  2064   /// ie. each node will be located in the \c [0,n-2]x[0,n-2] square.
  2065   /// The time complexity of the algorithm is O(n).
  2066   template <typename UGraph>
  2067   class PlanarDrawing {
  2068   public:
  2069 
  2070     UGRAPH_TYPEDEFS(typename UGraph);
  2071 
  2072     /// \brief The point type for store coordinates
  2073     typedef dim2::Point<int> Point;
  2074     /// \brief The map type for store coordinates
  2075     typedef typename UGraph::template NodeMap<Point> PointMap;
  2076 
  2077 
  2078     /// \brief Constructor
  2079     ///
  2080     /// Constructor
  2081     /// \pre The ugraph should be simple, ie. loop and parallel edge free. 
  2082     PlanarDrawing(const UGraph& ugraph)
  2083       : _ugraph(ugraph), _point_map(ugraph) {}
  2084 
  2085   private:
  2086 
  2087     template <typename AuxUGraph, typename AuxEmbeddingMap>
  2088     void drawing(const AuxUGraph& ugraph, 
  2089 		 const AuxEmbeddingMap& next, 
  2090 		 PointMap& point_map) {
  2091       UGRAPH_TYPEDEFS(typename AuxUGraph);
  2092 
  2093       typename AuxUGraph::template EdgeMap<Edge> prev(ugraph);
  2094 
  2095       for (NodeIt n(ugraph); n != INVALID; ++n) {
  2096 	Edge e = OutEdgeIt(ugraph, n);
  2097 	
  2098 	Edge p = e, l = e;
  2099 	
  2100 	e = next[e];
  2101 	while (e != l) {
  2102 	  prev[e] = p;
  2103 	  p = e;
  2104 	  e = next[e];
  2105 	}
  2106 	prev[e] = p;
  2107       }
  2108 
  2109       Node anode, bnode, cnode;
  2110 
  2111       {
  2112 	Edge e = EdgeIt(ugraph);
  2113 	anode = ugraph.source(e);
  2114 	bnode = ugraph.target(e);
  2115 	cnode = ugraph.target(next[ugraph.oppositeEdge(e)]);
  2116       }
  2117       
  2118       IterableBoolMap<AuxUGraph, Node> proper(ugraph, false);
  2119       typename AuxUGraph::template NodeMap<int> conn(ugraph, -1);
  2120 
  2121       conn[anode] = conn[bnode] = -2;      
  2122       {
  2123 	for (OutEdgeIt e(ugraph, anode); e != INVALID; ++e) {
  2124 	  Node m = ugraph.target(e);
  2125 	  if (conn[m] == -1) {
  2126 	    conn[m] = 1;
  2127 	  }
  2128 	}
  2129 	conn[cnode] = 2;
  2130 
  2131 	for (OutEdgeIt e(ugraph, bnode); e != INVALID; ++e) {
  2132 	  Node m = ugraph.target(e);
  2133 	  if (conn[m] == -1) {
  2134 	    conn[m] = 1;
  2135 	  } else if (conn[m] != -2) {
  2136 	    conn[m] += 1;	    
  2137 	    Edge pe = ugraph.oppositeEdge(e);
  2138 	    if (conn[ugraph.target(next[pe])] == -2) {
  2139 	      conn[m] -= 1;
  2140 	    }
  2141 	    if (conn[ugraph.target(prev[pe])] == -2) {
  2142 	      conn[m] -= 1;
  2143 	    }
  2144 
  2145 	    proper.set(m, conn[m] == 1);
  2146 	  }
  2147 	}
  2148       }
  2149       
  2150 
  2151       typename AuxUGraph::template EdgeMap<int> angle(ugraph, -1);
  2152 
  2153       while (proper.trueNum() != 0) {
  2154 	Node n = typename IterableBoolMap<AuxUGraph, Node>::TrueIt(proper);
  2155 	proper.set(n, false);
  2156 	conn[n] = -2;
  2157 
  2158 	for (OutEdgeIt e(ugraph, n); e != INVALID; ++e) {
  2159 	  Node m = ugraph.target(e);
  2160 	  if (conn[m] == -1) {
  2161 	    conn[m] = 1;
  2162 	  } else if (conn[m] != -2) {
  2163 	    conn[m] += 1;	    
  2164 	    Edge pe = ugraph.oppositeEdge(e);
  2165 	    if (conn[ugraph.target(next[pe])] == -2) {
  2166 	      conn[m] -= 1;
  2167 	    }
  2168 	    if (conn[ugraph.target(prev[pe])] == -2) {
  2169 	      conn[m] -= 1;
  2170 	    }
  2171 
  2172 	    proper.set(m, conn[m] == 1);
  2173 	  }
  2174 	}
  2175 
  2176 	{
  2177 	  Edge e = OutEdgeIt(ugraph, n);
  2178 	  Edge p = e, l = e;
  2179 	  
  2180 	  e = next[e];
  2181 	  while (e != l) {
  2182 	    
  2183 	    if (conn[ugraph.target(e)] == -2 && conn[ugraph.target(p)] == -2) {
  2184 	      Edge f = e;
  2185 	      angle[f] = 0;
  2186 	      f = next[ugraph.oppositeEdge(f)];
  2187 	      angle[f] = 1;
  2188 	      f = next[ugraph.oppositeEdge(f)];
  2189 	      angle[f] = 2;
  2190 	    }
  2191 	    
  2192 	    p = e;
  2193 	    e = next[e];
  2194 	  }
  2195 	  
  2196 	  if (conn[ugraph.target(e)] == -2 && conn[ugraph.target(p)] == -2) {
  2197 	    Edge f = e;
  2198 	    angle[f] = 0;
  2199 	    f = next[ugraph.oppositeEdge(f)];
  2200 	    angle[f] = 1;
  2201 	    f = next[ugraph.oppositeEdge(f)];
  2202 	    angle[f] = 2;
  2203 	  }
  2204 	}
  2205       }
  2206 
  2207       typename AuxUGraph::template NodeMap<Node> apred(ugraph, INVALID);
  2208       typename AuxUGraph::template NodeMap<Node> bpred(ugraph, INVALID);
  2209       typename AuxUGraph::template NodeMap<Node> cpred(ugraph, INVALID);
  2210 
  2211       typename AuxUGraph::template NodeMap<int> apredid(ugraph, -1);
  2212       typename AuxUGraph::template NodeMap<int> bpredid(ugraph, -1);
  2213       typename AuxUGraph::template NodeMap<int> cpredid(ugraph, -1);
  2214 
  2215       for (EdgeIt e(ugraph); e != INVALID; ++e) {
  2216 	if (angle[e] == angle[next[e]]) {
  2217 	  switch (angle[e]) {
  2218 	  case 2:
  2219 	    apred[ugraph.target(e)] = ugraph.source(e);
  2220 	    apredid[ugraph.target(e)] = ugraph.id(ugraph.source(e));
  2221 	    break;
  2222 	  case 1:
  2223 	    bpred[ugraph.target(e)] = ugraph.source(e);
  2224 	    bpredid[ugraph.target(e)] = ugraph.id(ugraph.source(e));
  2225 	    break;
  2226 	  case 0:
  2227 	    cpred[ugraph.target(e)] = ugraph.source(e);
  2228 	    cpredid[ugraph.target(e)] = ugraph.id(ugraph.source(e));
  2229 	    break;
  2230 	  }
  2231 	}
  2232       }
  2233 
  2234       cpred[anode] = INVALID;
  2235       cpred[bnode] = INVALID;
  2236 
  2237       std::vector<Node> aorder, border, corder; 
  2238 
  2239       {
  2240 	typename AuxUGraph::template NodeMap<bool> processed(ugraph, false);
  2241 	std::vector<Node> st;
  2242 	for (NodeIt n(ugraph); n != INVALID; ++n) {
  2243 	  if (!processed[n] && n != bnode && n != cnode) {
  2244 	    st.push_back(n);
  2245 	    processed[n] = true;
  2246 	    Node m = apred[n];
  2247 	    while (m != INVALID && !processed[m]) {
  2248 	      st.push_back(m);
  2249 	      processed[m] = true;
  2250 	      m = apred[m];
  2251 	    }
  2252 	    while (!st.empty()) {
  2253 	      aorder.push_back(st.back());
  2254 	      st.pop_back();
  2255 	    }
  2256 	  }
  2257 	}
  2258       }
  2259 
  2260       {
  2261 	typename AuxUGraph::template NodeMap<bool> processed(ugraph, false);
  2262 	std::vector<Node> st;
  2263 	for (NodeIt n(ugraph); n != INVALID; ++n) {
  2264 	  if (!processed[n] && n != cnode && n != anode) {
  2265 	    st.push_back(n);
  2266 	    processed[n] = true;
  2267 	    Node m = bpred[n];
  2268 	    while (m != INVALID && !processed[m]) {
  2269 	      st.push_back(m);
  2270 	      processed[m] = true;
  2271 	      m = bpred[m];
  2272 	    }
  2273 	    while (!st.empty()) {
  2274 	      border.push_back(st.back());
  2275 	      st.pop_back();
  2276 	    }
  2277 	  }
  2278 	}
  2279       }
  2280 
  2281       {
  2282 	typename AuxUGraph::template NodeMap<bool> processed(ugraph, false);
  2283 	std::vector<Node> st;
  2284 	for (NodeIt n(ugraph); n != INVALID; ++n) {
  2285 	  if (!processed[n] && n != anode && n != bnode) {
  2286 	    st.push_back(n);
  2287 	    processed[n] = true;
  2288 	    Node m = cpred[n];
  2289 	    while (m != INVALID && !processed[m]) {
  2290 	      st.push_back(m);
  2291 	      processed[m] = true;
  2292 	      m = cpred[m];
  2293 	    }
  2294 	    while (!st.empty()) {
  2295 	      corder.push_back(st.back());
  2296 	      st.pop_back();
  2297 	    }
  2298 	  }
  2299 	}
  2300       }
  2301 
  2302       typename AuxUGraph::template NodeMap<int> atree(ugraph, 0);
  2303       for (int i = aorder.size() - 1; i >= 0; --i) {
  2304 	Node n = aorder[i];
  2305 	atree[n] = 1;
  2306 	for (OutEdgeIt e(ugraph, n); e != INVALID; ++e) {
  2307 	  if (apred[ugraph.target(e)] == n) {
  2308 	    atree[n] += atree[ugraph.target(e)];
  2309 	  }
  2310 	} 
  2311       }
  2312 
  2313       typename AuxUGraph::template NodeMap<int> btree(ugraph, 0);
  2314       for (int i = border.size() - 1; i >= 0; --i) {
  2315 	Node n = border[i];
  2316 	btree[n] = 1;
  2317 	for (OutEdgeIt e(ugraph, n); e != INVALID; ++e) {
  2318 	  if (bpred[ugraph.target(e)] == n) {
  2319 	    btree[n] += btree[ugraph.target(e)];
  2320 	  }
  2321 	} 
  2322       }
  2323       
  2324       typename AuxUGraph::template NodeMap<int> apath(ugraph, 0);
  2325       apath[bnode] = apath[cnode] = 1;
  2326       typename AuxUGraph::template NodeMap<int> apath_btree(ugraph, 0);
  2327       apath_btree[bnode] = btree[bnode];
  2328       for (int i = 1; i < int(aorder.size()); ++i) {
  2329 	Node n = aorder[i];
  2330 	apath[n] = apath[apred[n]] + 1;
  2331 	apath_btree[n] = btree[n] + apath_btree[apred[n]];
  2332       }
  2333 
  2334       typename AuxUGraph::template NodeMap<int> bpath_atree(ugraph, 0);
  2335       bpath_atree[anode] = atree[anode];
  2336       for (int i = 1; i < int(border.size()); ++i) {
  2337 	Node n = border[i];
  2338 	bpath_atree[n] = atree[n] + bpath_atree[bpred[n]];
  2339       }
  2340 
  2341       typename AuxUGraph::template NodeMap<int> cpath(ugraph, 0);
  2342       cpath[anode] = cpath[bnode] = 1;
  2343       typename AuxUGraph::template NodeMap<int> cpath_atree(ugraph, 0);
  2344       cpath_atree[anode] = atree[anode];
  2345       typename AuxUGraph::template NodeMap<int> cpath_btree(ugraph, 0);
  2346       cpath_btree[bnode] = btree[bnode];
  2347       for (int i = 1; i < int(corder.size()); ++i) {
  2348 	Node n = corder[i];
  2349 	cpath[n] = cpath[cpred[n]] + 1;
  2350 	cpath_atree[n] = atree[n] + cpath_atree[cpred[n]];
  2351 	cpath_btree[n] = btree[n] + cpath_btree[cpred[n]];
  2352       }
  2353 
  2354       typename AuxUGraph::template NodeMap<int> third(ugraph);
  2355       for (NodeIt n(ugraph); n != INVALID; ++n) {
  2356 	point_map[n].x = 
  2357 	  bpath_atree[n] + cpath_atree[n] - atree[n] - cpath[n] + 1;
  2358 	point_map[n].y = 
  2359 	  cpath_btree[n] + apath_btree[n] - btree[n] - apath[n] + 1;
  2360       }
  2361       
  2362     }
  2363 
  2364   public:
  2365 
  2366     /// \brief Calculates the node locations
  2367     ///
  2368     /// This function calculates the node locations.
  2369     bool run() {
  2370       PlanarEmbedding<UGraph> pe(_ugraph);
  2371       if (!pe.run()) return false;
  2372       
  2373       run(pe);
  2374       return true;
  2375     }
  2376 
  2377     /// \brief Calculates the node locations according to a
  2378     /// combinatorical embedding
  2379     ///
  2380     /// This function calculates the node locations. The \c embedding
  2381     /// parameter should contain a valid combinatorical embedding, ie.
  2382     /// a valid cyclic order of the edges.
  2383     template <typename EmbeddingMap>
  2384     void run(const EmbeddingMap& embedding) {
  2385       typedef SmartUEdgeSet<UGraph> AuxUGraph; 
  2386 
  2387       if (3 * countNodes(_ugraph) - 6 == countUEdges(_ugraph)) {
  2388 	drawing(_ugraph, embedding, _point_map);
  2389 	return;
  2390       }
  2391 
  2392       AuxUGraph aux_ugraph(_ugraph);
  2393       typename AuxUGraph::template EdgeMap<typename AuxUGraph::Edge> 
  2394 	aux_embedding(aux_ugraph);
  2395 
  2396       {
  2397 
  2398 	typename UGraph::template UEdgeMap<typename AuxUGraph::UEdge> 
  2399 	  ref(_ugraph);
  2400 	
  2401 	for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
  2402 	  ref[e] = aux_ugraph.addEdge(_ugraph.source(e), _ugraph.target(e));
  2403 	}
  2404 
  2405 	for (UEdgeIt e(_ugraph); e != INVALID; ++e) {
  2406 	  Edge ee = embedding[_ugraph.direct(e, true)];
  2407 	  aux_embedding[aux_ugraph.direct(ref[e], true)] = 
  2408 	    aux_ugraph.direct(ref[ee], _ugraph.direction(ee));
  2409 	  ee = embedding[_ugraph.direct(e, false)];
  2410 	  aux_embedding[aux_ugraph.direct(ref[e], false)] = 
  2411 	    aux_ugraph.direct(ref[ee], _ugraph.direction(ee));
  2412 	}
  2413       }
  2414       _planarity_bits::makeConnected(aux_ugraph, aux_embedding);
  2415       _planarity_bits::makeBiNodeConnected(aux_ugraph, aux_embedding);
  2416       _planarity_bits::makeMaxPlanar(aux_ugraph, aux_embedding);
  2417       drawing(aux_ugraph, aux_embedding, _point_map);
  2418     }
  2419 
  2420     /// \brief The coordinate of the given node
  2421     ///
  2422     /// The coordinate of the given node.
  2423     Point operator[](const Node& node) {
  2424       return _point_map[node];
  2425     }
  2426 
  2427     /// \brief Returns the grid embedding in a \e NodeMap.
  2428     ///
  2429     /// Returns the grid embedding in a \e NodeMap of \c dim2::Point<int> .
  2430     const PointMap& coords() const {
  2431       return _point_map;
  2432     }
  2433 
  2434   private:
  2435     
  2436     const UGraph& _ugraph;
  2437     PointMap _point_map;
  2438     
  2439   };
  2440 
  2441 }
  2442 
  2443 #endif