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