lemon/planarity.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 19 Feb 2010 14:08:32 +0100
changeset 844 a6eb9698c321
parent 798 58c330ad0b5c
child 877 141f9c0db4a3
permissions -rw-r--r--
Support tolerance technique for BellmanFord (#51)

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