lemon-main: lemon/planarity.h@ddd3c0d3d9bf

     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-2010

     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  */

    19 #ifndef LEMON_PLANARITY_H

    20 #define LEMON_PLANARITY_H

    22 /// \ingroup planar

    23 /// \file

    24 /// \brief Planarity checking, embedding, drawing and coloring

    26 #include <vector>

    27 #include <list>

    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>

    40 namespace lemon {

    42   namespace _planarity_bits {

    44     template <typename Graph>

    45     struct PlanarityVisitor : DfsVisitor<Graph> {

    47       TEMPLATE_GRAPH_TYPEDEFS(Graph);

    49       typedef typename Graph::template NodeMap<Arc> PredMap;

    51       typedef typename Graph::template EdgeMap<bool> TreeMap;

    53       typedef typename Graph::template NodeMap<int> OrderMap;

    54       typedef std::vector<Node> OrderList;

    56       typedef typename Graph::template NodeMap<int> LowMap;

    57       typedef typename Graph::template NodeMap<int> AncestorMap;

    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) {}

    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       }

    74       void discover(const Arc& arc) {

    75         Node source = _graph.source(arc);

    76         Node target = _graph.target(arc);

    78         _tree_map[arc] = true;

    79         _pred_map[target] = arc;

    80       }

    82       void examine(const Arc& arc) {

    83         Node source = _graph.source(arc);

    84         Node target = _graph.target(arc);

    86         if (_order_map[target] < _order_map[source] && !_tree_map[arc]) {

    87           if (_low_map[source] > _order_map[target]) {

    88             _low_map[source] = _order_map[target];

    89           }

    90           if (_ancestor_map[source] > _order_map[target]) {

    91             _ancestor_map[source] = _order_map[target];

    92           }

    93         }

    94       }

    96       void backtrack(const Arc& arc) {

    97         Node source = _graph.source(arc);

    98         Node target = _graph.target(arc);

   100         if (_low_map[source] > _low_map[target]) {

   101           _low_map[source] = _low_map[target];

   102         }

   103       }

   105       const Graph& _graph;

   106       PredMap& _pred_map;

   107       TreeMap& _tree_map;

   108       OrderMap& _order_map;

   109       OrderList& _order_list;

   110       AncestorMap& _ancestor_map;

   111       LowMap& _low_map;

   112     };

   114     template <typename Graph, bool embedding = true>

   115     struct NodeDataNode {

   116       int prev, next;

   117       int visited;

   118       typename Graph::Arc first;

   119       bool inverted;

   120     };

   122     template <typename Graph>

   123     struct NodeDataNode<Graph, false> {

   124       int prev, next;

   125       int visited;

   126     };

   128     template <typename Graph>

   129     struct ChildListNode {

   130       typedef typename Graph::Node Node;

   131       Node first;

   132       Node prev, next;

   133     };

   135     template <typename Graph>

   136     struct ArcListNode {

   137       typename Graph::Arc prev, next;

   138     };

   140     template <typename Graph>

   141     class PlanarityChecking {

   142     private:

   144       TEMPLATE_GRAPH_TYPEDEFS(Graph);

   146       const Graph& _graph;

   148     private:

   150       typedef typename Graph::template NodeMap<Arc> PredMap;

   152       typedef typename Graph::template EdgeMap<bool> TreeMap;

   154       typedef typename Graph::template NodeMap<int> OrderMap;

   155       typedef std::vector<Node> OrderList;

   157       typedef typename Graph::template NodeMap<int> LowMap;

   158       typedef typename Graph::template NodeMap<int> AncestorMap;

   160       typedef _planarity_bits::NodeDataNode<Graph> NodeDataNode;

   161       typedef std::vector<NodeDataNode> NodeData;

   163       typedef _planarity_bits::ChildListNode<Graph> ChildListNode;

   164       typedef typename Graph::template NodeMap<ChildListNode> ChildLists;

   166       typedef typename Graph::template NodeMap<std::list<int> > MergeRoots;

   168       typedef typename Graph::template NodeMap<bool> EmbedArc;

   170     public:

   172       PlanarityChecking(const Graph& graph) : _graph(graph) {}

   174       bool run() {

   175         typedef _planarity_bits::PlanarityVisitor<Graph> Visitor;

   177         PredMap pred_map(_graph, INVALID);

   178         TreeMap tree_map(_graph, false);

   180         OrderMap order_map(_graph, -1);

   181         OrderList order_list;

   183         AncestorMap ancestor_map(_graph, -1);

   184         LowMap low_map(_graph, -1);

   186         Visitor visitor(_graph, pred_map, tree_map,

   187                         order_map, order_list, ancestor_map, low_map);

   188         DfsVisit<Graph, Visitor> visit(_graph, visitor);

   189         visit.run();

   191         ChildLists child_lists(_graph);

   192         createChildLists(tree_map, order_map, low_map, child_lists);

   194         NodeData node_data(2 * order_list.size());

   196         EmbedArc embed_arc(_graph, false);

   198         MergeRoots merge_roots(_graph);

   200         for (int i = order_list.size() - 1; i >= 0; --i) {

   202           Node node = order_list[i];

   204           Node source = node;

   205           for (OutArcIt e(_graph, node); e != INVALID; ++e) {

   206             Node target = _graph.target(e);

   208             if (order_map[source] < order_map[target] && tree_map[e]) {

   209               initFace(target, node_data, order_map, order_list);

   210             }

   211           }

   213           for (OutArcIt e(_graph, node); e != INVALID; ++e) {

   214             Node target = _graph.target(e);

   216             if (order_map[source] < order_map[target] && !tree_map[e]) {

   217               embed_arc[target] = true;

   218               walkUp(target, source, i, pred_map, low_map,

   219                      order_map, order_list, node_data, merge_roots);

   220             }

   221           }

   223           for (typename MergeRoots::Value::iterator it =

   224                  merge_roots[node].begin();

   225                it != merge_roots[node].end(); ++it) {

   226             int rn = *it;

   227             walkDown(rn, i, node_data, order_list, child_lists,

   228                      ancestor_map, low_map, embed_arc, merge_roots);

   229           }

   230           merge_roots[node].clear();

   232           for (OutArcIt e(_graph, node); e != INVALID; ++e) {

   233             Node target = _graph.target(e);

   235             if (order_map[source] < order_map[target] && !tree_map[e]) {

   236               if (embed_arc[target]) {

   237                 return false;

   238               }

   239             }

   240           }

   241         }

   243         return true;

   244       }

   246     private:

   248       void createChildLists(const TreeMap& tree_map, const OrderMap& order_map,

   249                             const LowMap& low_map, ChildLists& child_lists) {

   251         for (NodeIt n(_graph); n != INVALID; ++n) {

   252           Node source = n;

   254           std::vector<Node> targets;

   255           for (OutArcIt e(_graph, n); e != INVALID; ++e) {

   256             Node target = _graph.target(e);

   258             if (order_map[source] < order_map[target] && tree_map[e]) {

   259               targets.push_back(target);

   260             }

   261           }

   263           if (targets.size() == 0) {

   264             child_lists[source].first = INVALID;

   265           } else if (targets.size() == 1) {

   266             child_lists[source].first = targets[0];

   267             child_lists[targets[0]].prev = INVALID;

   268             child_lists[targets[0]].next = INVALID;

   269           } else {

   270             radixSort(targets.begin(), targets.end(), mapToFunctor(low_map));

   271             for (int i = 1; i < int(targets.size()); ++i) {

   272               child_lists[targets[i]].prev = targets[i - 1];

   273               child_lists[targets[i - 1]].next = targets[i];

   274             }

   275             child_lists[targets.back()].next = INVALID;

   276             child_lists[targets.front()].prev = INVALID;

   277             child_lists[source].first = targets.front();

   278           }

   279         }

   280       }

   282       void walkUp(const Node& node, Node root, int rorder,

   283                   const PredMap& pred_map, const LowMap& low_map,

   284                   const OrderMap& order_map, const OrderList& order_list,

   285                   NodeData& node_data, MergeRoots& merge_roots) {

   287         int na, nb;

   288         bool da, db;

   290         na = nb = order_map[node];

   291         da = true; db = false;

   293         while (true) {

   295           if (node_data[na].visited == rorder) break;

   296           if (node_data[nb].visited == rorder) break;

   298           node_data[na].visited = rorder;

   299           node_data[nb].visited = rorder;

   301           int rn = -1;

   303           if (na >= int(order_list.size())) {

   304             rn = na;

   305           } else if (nb >= int(order_list.size())) {

   306             rn = nb;

   307           }

   309           if (rn == -1) {

   310             int nn;

   312             nn = da ? node_data[na].prev : node_data[na].next;

   313             da = node_data[nn].prev != na;

   314             na = nn;

   316             nn = db ? node_data[nb].prev : node_data[nb].next;

   317             db = node_data[nn].prev != nb;

   318             nb = nn;

   320           } else {

   322             Node rep = order_list[rn - order_list.size()];

   323             Node parent = _graph.source(pred_map[rep]);

   325             if (low_map[rep] < rorder) {

   326               merge_roots[parent].push_back(rn);

   327             } else {

   328               merge_roots[parent].push_front(rn);

   329             }

   331             if (parent != root) {

   332               na = nb = order_map[parent];

   333               da = true; db = false;

   334             } else {

   335               break;

   336             }

   337           }

   338         }

   339       }

   341       void walkDown(int rn, int rorder, NodeData& node_data,

   342                     OrderList& order_list, ChildLists& child_lists,

   343                     AncestorMap& ancestor_map, LowMap& low_map,

   344                     EmbedArc& embed_arc, MergeRoots& merge_roots) {

   346         std::vector<std::pair<int, bool> > merge_stack;

   348         for (int di = 0; di < 2; ++di) {

   349           bool rd = di == 0;

   350           int pn = rn;

   351           int n = rd ? node_data[rn].next : node_data[rn].prev;

   353           while (n != rn) {

   355             Node node = order_list[n];

   357             if (embed_arc[node]) {

   359               // Merging components on the critical path

   360               while (!merge_stack.empty()) {

   362                 // Component root

   363                 int cn = merge_stack.back().first;

   364                 bool cd = merge_stack.back().second;

   365                 merge_stack.pop_back();

   367                 // Parent of component

   368                 int dn = merge_stack.back().first;

   369                 bool dd = merge_stack.back().second;

   370                 merge_stack.pop_back();

   372                 Node parent = order_list[dn];

   374                 // Erasing from merge_roots

   375                 merge_roots[parent].pop_front();

   377                 Node child = order_list[cn - order_list.size()];

   379                 // Erasing from child_lists

   380                 if (child_lists[child].prev != INVALID) {

   381                   child_lists[child_lists[child].prev].next =

   382                     child_lists[child].next;

   383                 } else {

   384                   child_lists[parent].first = child_lists[child].next;

   385                 }

   387                 if (child_lists[child].next != INVALID) {

   388                   child_lists[child_lists[child].next].prev =

   389                     child_lists[child].prev;

   390                 }

   392                 // Merging external faces

   393                 {

   394                   int en = cn;

   395                   cn = cd ? node_data[cn].prev : node_data[cn].next;

   396                   cd = node_data[cn].next == en;

   398                 }

   400                 if (cd) node_data[cn].next = dn; else node_data[cn].prev = dn;

   401                 if (dd) node_data[dn].prev = cn; else node_data[dn].next = cn;

   403               }

   405               bool d = pn == node_data[n].prev;

   407               if (node_data[n].prev == node_data[n].next &&

   408                   node_data[n].inverted) {

   409                 d = !d;

   410               }

   412               // Embedding arc into external face

   413               if (rd) node_data[rn].next = n; else node_data[rn].prev = n;

   414               if (d) node_data[n].prev = rn; else node_data[n].next = rn;

   415               pn = rn;

   417               embed_arc[order_list[n]] = false;

   418             }

   420             if (!merge_roots[node].empty()) {

   422               bool d = pn == node_data[n].prev;

   424               merge_stack.push_back(std::make_pair(n, d));

   426               int rn = merge_roots[node].front();

   428               int xn = node_data[rn].next;

   429               Node xnode = order_list[xn];

   431               int yn = node_data[rn].prev;

   432               Node ynode = order_list[yn];

   434               bool rd;

   435               if (!external(xnode, rorder, child_lists,

   436                             ancestor_map, low_map)) {

   437                 rd = true;

   438               } else if (!external(ynode, rorder, child_lists,

   439                                    ancestor_map, low_map)) {

   440                 rd = false;

   441               } else if (pertinent(xnode, embed_arc, merge_roots)) {

   442                 rd = true;

   443               } else {

   444                 rd = false;

   445               }

   447               merge_stack.push_back(std::make_pair(rn, rd));

   449               pn = rn;

   450               n = rd ? xn : yn;

   452             } else if (!external(node, rorder, child_lists,

   453                                  ancestor_map, low_map)) {

   454               int nn = (node_data[n].next != pn ?

   455                         node_data[n].next : node_data[n].prev);

   457               bool nd = n == node_data[nn].prev;

   459               if (nd) node_data[nn].prev = pn;

   460               else node_data[nn].next = pn;

   462               if (n == node_data[pn].prev) node_data[pn].prev = nn;

   463               else node_data[pn].next = nn;

   465               node_data[nn].inverted =

   466                 (node_data[nn].prev == node_data[nn].next && nd != rd);

   468               n = nn;

   469             }

   470             else break;

   472           }

   474           if (!merge_stack.empty() || n == rn) {

   475             break;

   476           }

   477         }

   478       }

   480       void initFace(const Node& node, NodeData& node_data,

   481                     const OrderMap& order_map, const OrderList& order_list) {

   482         int n = order_map[node];

   483         int rn = n + order_list.size();

   485         node_data[n].next = node_data[n].prev = rn;

   486         node_data[rn].next = node_data[rn].prev = n;

   488         node_data[n].visited = order_list.size();

   489         node_data[rn].visited = order_list.size();

   491       }

   493       bool external(const Node& node, int rorder,

   494                     ChildLists& child_lists, AncestorMap& ancestor_map,

   495                     LowMap& low_map) {

   496         Node child = child_lists[node].first;

   498         if (child != INVALID) {

   499           if (low_map[child] < rorder) return true;

   500         }

   502         if (ancestor_map[node] < rorder) return true;

   504         return false;

   505       }

   507       bool pertinent(const Node& node, const EmbedArc& embed_arc,

   508                      const MergeRoots& merge_roots) {

   509         return !merge_roots[node].empty() || embed_arc[node];

   510       }

   512     };

   514   }

   516   /// \ingroup planar

   517   ///

   518   /// \brief Planarity checking of an undirected simple graph

   519   ///

   520   /// This function implements the Boyer-Myrvold algorithm for

   521   /// planarity checking of an undirected simple graph. It is a simplified

   522   /// version of the PlanarEmbedding algorithm class because neither

   523   /// the embedding nor the Kuratowski subdivisons are computed.

   524   template <typename GR>

   525   bool checkPlanarity(const GR& graph) {

   526     _planarity_bits::PlanarityChecking<GR> pc(graph);

   527     return pc.run();

   528   }

   530   /// \ingroup planar

   531   ///

   532   /// \brief Planar embedding of an undirected simple graph

   533   ///

   534   /// This class implements the Boyer-Myrvold algorithm for planar

   535   /// embedding of an undirected simple graph. The planar embedding is an

   536   /// ordering of the outgoing edges of the nodes, which is a possible

   537   /// configuration to draw the graph in the plane. If there is not

   538   /// such ordering then the graph contains a K<sub>5</sub> (full graph

   539   /// with 5 nodes) or a K<sub>3,3</sub> (complete bipartite graph on

   540   /// 3 Red and 3 Blue nodes) subdivision.

   541   ///

   542   /// The current implementation calculates either an embedding or a

   543   /// Kuratowski subdivision. The running time of the algorithm is O(n).

   544   ///

   545   /// \see PlanarDrawing, checkPlanarity()

   546   template <typename Graph>

   547   class PlanarEmbedding {

   548   private:

   550     TEMPLATE_GRAPH_TYPEDEFS(Graph);

   552     const Graph& _graph;

   553     typename Graph::template ArcMap<Arc> _embedding;

   555     typename Graph::template EdgeMap<bool> _kuratowski;

   557   private:

   559     typedef typename Graph::template NodeMap<Arc> PredMap;

   561     typedef typename Graph::template EdgeMap<bool> TreeMap;

   563     typedef typename Graph::template NodeMap<int> OrderMap;

   564     typedef std::vector<Node> OrderList;

   566     typedef typename Graph::template NodeMap<int> LowMap;

   567     typedef typename Graph::template NodeMap<int> AncestorMap;

   569     typedef _planarity_bits::NodeDataNode<Graph> NodeDataNode;

   570     typedef std::vector<NodeDataNode> NodeData;

   572     typedef _planarity_bits::ChildListNode<Graph> ChildListNode;

   573     typedef typename Graph::template NodeMap<ChildListNode> ChildLists;

   575     typedef typename Graph::template NodeMap<std::list<int> > MergeRoots;

   577     typedef typename Graph::template NodeMap<Arc> EmbedArc;

   579     typedef _planarity_bits::ArcListNode<Graph> ArcListNode;

   580     typedef typename Graph::template ArcMap<ArcListNode> ArcLists;

   582     typedef typename Graph::template NodeMap<bool> FlipMap;

   584     typedef typename Graph::template NodeMap<int> TypeMap;

   586     enum IsolatorNodeType {

   587       HIGHX = 6, LOWX = 7,

   588       HIGHY = 8, LOWY = 9,

   589       ROOT = 10, PERTINENT = 11,

   590       INTERNAL = 12

   591     };

   593   public:

   595     /// \brief The map type for storing the embedding

   596     ///

   597     /// The map type for storing the embedding.

   598     /// \see embeddingMap()

   599     typedef typename Graph::template ArcMap<Arc> EmbeddingMap;

   601     /// \brief Constructor

   602     ///

   603     /// Constructor.

   604     /// \pre The graph must be simple, i.e. it should not

   605     /// contain parallel or loop arcs.

   606     PlanarEmbedding(const Graph& graph)

   607       : _graph(graph), _embedding(_graph), _kuratowski(graph, false) {}

   609     /// \brief Run the algorithm.

   610     ///

   611     /// This function runs the algorithm.

   612     /// \param kuratowski If this parameter is set to \c false, then the

   613     /// algorithm does not compute a Kuratowski subdivision.

   614     /// \return \c true if the graph is planar.

   615     bool run(bool kuratowski = true) {

   616       typedef _planarity_bits::PlanarityVisitor<Graph> Visitor;

   618       PredMap pred_map(_graph, INVALID);

   619       TreeMap tree_map(_graph, false);

   621       OrderMap order_map(_graph, -1);

   622       OrderList order_list;

   624       AncestorMap ancestor_map(_graph, -1);

   625       LowMap low_map(_graph, -1);

   627       Visitor visitor(_graph, pred_map, tree_map,

   628                       order_map, order_list, ancestor_map, low_map);

   629       DfsVisit<Graph, Visitor> visit(_graph, visitor);

   630       visit.run();

   632       ChildLists child_lists(_graph);

   633       createChildLists(tree_map, order_map, low_map, child_lists);

   635       NodeData node_data(2 * order_list.size());

   637       EmbedArc embed_arc(_graph, INVALID);

   639       MergeRoots merge_roots(_graph);

   641       ArcLists arc_lists(_graph);

   643       FlipMap flip_map(_graph, false);

   645       for (int i = order_list.size() - 1; i >= 0; --i) {

   647         Node node = order_list[i];

   649         node_data[i].first = INVALID;

   651         Node source = node;

   652         for (OutArcIt e(_graph, node); e != INVALID; ++e) {

   653           Node target = _graph.target(e);

   655           if (order_map[source] < order_map[target] && tree_map[e]) {

   656             initFace(target, arc_lists, node_data,

   657                      pred_map, order_map, order_list);

   658           }

   659         }

   661         for (OutArcIt e(_graph, node); e != INVALID; ++e) {

   662           Node target = _graph.target(e);

   664           if (order_map[source] < order_map[target] && !tree_map[e]) {

   665             embed_arc[target] = e;

   666             walkUp(target, source, i, pred_map, low_map,

   667                    order_map, order_list, node_data, merge_roots);

   668           }

   669         }

   671         for (typename MergeRoots::Value::iterator it =

   672                merge_roots[node].begin(); it != merge_roots[node].end(); ++it) {

   673           int rn = *it;

   674           walkDown(rn, i, node_data, arc_lists, flip_map, order_list,

   675                    child_lists, ancestor_map, low_map, embed_arc, merge_roots);

   676         }

   677         merge_roots[node].clear();

   679         for (OutArcIt e(_graph, node); e != INVALID; ++e) {

   680           Node target = _graph.target(e);

   682           if (order_map[source] < order_map[target] && !tree_map[e]) {

   683             if (embed_arc[target] != INVALID) {

   684               if (kuratowski) {

   685                 isolateKuratowski(e, node_data, arc_lists, flip_map,

   686                                   order_map, order_list, pred_map, child_lists,

   687                                   ancestor_map, low_map,

   688                                   embed_arc, merge_roots);

   689               }

   690               return false;

   691             }

   692           }

   693         }

   694       }

   696       for (int i = 0; i < int(order_list.size()); ++i) {

   698         mergeRemainingFaces(order_list[i], node_data, order_list, order_map,

   699                             child_lists, arc_lists);

   700         storeEmbedding(order_list[i], node_data, order_map, pred_map,

   701                        arc_lists, flip_map);

   702       }

   704       return true;

   705     }

   707     /// \brief Give back the successor of an arc

   708     ///

   709     /// This function gives back the successor of an arc. It makes

   710     /// possible to query the cyclic order of the outgoing arcs from

   711     /// a node.

   712     Arc next(const Arc& arc) const {

   713       return _embedding[arc];

   714     }

   716     /// \brief Give back the calculated embedding map

   717     ///

   718     /// This function gives back the calculated embedding map, which

   719     /// contains the successor of each arc in the cyclic order of the

   720     /// outgoing arcs of its source node.

   721     const EmbeddingMap& embeddingMap() const {

   722       return _embedding;

   723     }

   725     /// \brief Give back \c true if the given edge is in the Kuratowski

   726     /// subdivision

   727     ///

   728     /// This function gives back \c true if the given edge is in the found

   729     /// Kuratowski subdivision.

   730     /// \pre The \c run() function must be called with \c true parameter

   731     /// before using this function.

   732     bool kuratowski(const Edge& edge) const {

   733       return _kuratowski[edge];

   734     }

   736   private:

   738     void createChildLists(const TreeMap& tree_map, const OrderMap& order_map,

   739                           const LowMap& low_map, ChildLists& child_lists) {

   741       for (NodeIt n(_graph); n != INVALID; ++n) {

   742         Node source = n;

   744         std::vector<Node> targets;

   745         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

   746           Node target = _graph.target(e);

   748           if (order_map[source] < order_map[target] && tree_map[e]) {

   749             targets.push_back(target);

   750           }

   751         }

   753         if (targets.size() == 0) {

   754           child_lists[source].first = INVALID;

   755         } else if (targets.size() == 1) {

   756           child_lists[source].first = targets[0];

   757           child_lists[targets[0]].prev = INVALID;

   758           child_lists[targets[0]].next = INVALID;

   759         } else {

   760           radixSort(targets.begin(), targets.end(), mapToFunctor(low_map));

   761           for (int i = 1; i < int(targets.size()); ++i) {

   762             child_lists[targets[i]].prev = targets[i - 1];

   763             child_lists[targets[i - 1]].next = targets[i];

   764           }

   765           child_lists[targets.back()].next = INVALID;

   766           child_lists[targets.front()].prev = INVALID;

   767           child_lists[source].first = targets.front();

   768         }

   769       }

   770     }

   772     void walkUp(const Node& node, Node root, int rorder,

   773                 const PredMap& pred_map, const LowMap& low_map,

   774                 const OrderMap& order_map, const OrderList& order_list,

   775                 NodeData& node_data, MergeRoots& merge_roots) {

   777       int na, nb;

   778       bool da, db;

   780       na = nb = order_map[node];

   781       da = true; db = false;

   783       while (true) {

   785         if (node_data[na].visited == rorder) break;

   786         if (node_data[nb].visited == rorder) break;

   788         node_data[na].visited = rorder;

   789         node_data[nb].visited = rorder;

   791         int rn = -1;

   793         if (na >= int(order_list.size())) {

   794           rn = na;

   795         } else if (nb >= int(order_list.size())) {

   796           rn = nb;

   797         }

   799         if (rn == -1) {

   800           int nn;

   802           nn = da ? node_data[na].prev : node_data[na].next;

   803           da = node_data[nn].prev != na;

   804           na = nn;

   806           nn = db ? node_data[nb].prev : node_data[nb].next;

   807           db = node_data[nn].prev != nb;

   808           nb = nn;

   810         } else {

   812           Node rep = order_list[rn - order_list.size()];

   813           Node parent = _graph.source(pred_map[rep]);

   815           if (low_map[rep] < rorder) {

   816             merge_roots[parent].push_back(rn);

   817           } else {

   818             merge_roots[parent].push_front(rn);

   819           }

   821           if (parent != root) {

   822             na = nb = order_map[parent];

   823             da = true; db = false;

   824           } else {

   825             break;

   826           }

   827         }

   828       }

   829     }

   831     void walkDown(int rn, int rorder, NodeData& node_data,

   832                   ArcLists& arc_lists, FlipMap& flip_map,

   833                   OrderList& order_list, ChildLists& child_lists,

   834                   AncestorMap& ancestor_map, LowMap& low_map,

   835                   EmbedArc& embed_arc, MergeRoots& merge_roots) {

   837       std::vector<std::pair<int, bool> > merge_stack;

   839       for (int di = 0; di < 2; ++di) {

   840         bool rd = di == 0;

   841         int pn = rn;

   842         int n = rd ? node_data[rn].next : node_data[rn].prev;

   844         while (n != rn) {

   846           Node node = order_list[n];

   848           if (embed_arc[node] != INVALID) {

   850             // Merging components on the critical path

   851             while (!merge_stack.empty()) {

   853               // Component root

   854               int cn = merge_stack.back().first;

   855               bool cd = merge_stack.back().second;

   856               merge_stack.pop_back();

   858               // Parent of component

   859               int dn = merge_stack.back().first;

   860               bool dd = merge_stack.back().second;

   861               merge_stack.pop_back();

   863               Node parent = order_list[dn];

   865               // Erasing from merge_roots

   866               merge_roots[parent].pop_front();

   868               Node child = order_list[cn - order_list.size()];

   870               // Erasing from child_lists

   871               if (child_lists[child].prev != INVALID) {

   872                 child_lists[child_lists[child].prev].next =

   873                   child_lists[child].next;

   874               } else {

   875                 child_lists[parent].first = child_lists[child].next;

   876               }

   878               if (child_lists[child].next != INVALID) {

   879                 child_lists[child_lists[child].next].prev =

   880                   child_lists[child].prev;

   881               }

   883               // Merging arcs + flipping

   884               Arc de = node_data[dn].first;

   885               Arc ce = node_data[cn].first;

   887               flip_map[order_list[cn - order_list.size()]] = cd != dd;

   888               if (cd != dd) {

   889                 std::swap(arc_lists[ce].prev, arc_lists[ce].next);

   890                 ce = arc_lists[ce].prev;

   891                 std::swap(arc_lists[ce].prev, arc_lists[ce].next);

   892               }

   894               {

   895                 Arc dne = arc_lists[de].next;

   896                 Arc cne = arc_lists[ce].next;

   898                 arc_lists[de].next = cne;

   899                 arc_lists[ce].next = dne;

   901                 arc_lists[dne].prev = ce;

   902                 arc_lists[cne].prev = de;

   903               }

   905               if (dd) {

   906                 node_data[dn].first = ce;

   907               }

   909               // Merging external faces

   910               {

   911                 int en = cn;

   912                 cn = cd ? node_data[cn].prev : node_data[cn].next;

   913                 cd = node_data[cn].next == en;

   915                  if (node_data[cn].prev == node_data[cn].next &&

   916                     node_data[cn].inverted) {

   917                    cd = !cd;

   918                  }

   919               }

   921               if (cd) node_data[cn].next = dn; else node_data[cn].prev = dn;

   922               if (dd) node_data[dn].prev = cn; else node_data[dn].next = cn;

   924             }

   926             bool d = pn == node_data[n].prev;

   928             if (node_data[n].prev == node_data[n].next &&

   929                 node_data[n].inverted) {

   930               d = !d;

   931             }

   933             // Add new arc

   934             {

   935               Arc arc = embed_arc[node];

   936               Arc re = node_data[rn].first;

   938               arc_lists[arc_lists[re].next].prev = arc;

   939               arc_lists[arc].next = arc_lists[re].next;

   940               arc_lists[arc].prev = re;

   941               arc_lists[re].next = arc;

   943               if (!rd) {

   944                 node_data[rn].first = arc;

   945               }

   947               Arc rev = _graph.oppositeArc(arc);

   948               Arc e = node_data[n].first;

   950               arc_lists[arc_lists[e].next].prev = rev;

   951               arc_lists[rev].next = arc_lists[e].next;

   952               arc_lists[rev].prev = e;

   953               arc_lists[e].next = rev;

   955               if (d) {

   956                 node_data[n].first = rev;

   957               }

   959             }

   961             // Embedding arc into external face

   962             if (rd) node_data[rn].next = n; else node_data[rn].prev = n;

   963             if (d) node_data[n].prev = rn; else node_data[n].next = rn;

   964             pn = rn;

   966             embed_arc[order_list[n]] = INVALID;

   967           }

   969           if (!merge_roots[node].empty()) {

   971             bool d = pn == node_data[n].prev;

   972             if (node_data[n].prev == node_data[n].next &&

   973                 node_data[n].inverted) {

   974               d = !d;

   975             }

   977             merge_stack.push_back(std::make_pair(n, d));

   979             int rn = merge_roots[node].front();

   981             int xn = node_data[rn].next;

   982             Node xnode = order_list[xn];

   984             int yn = node_data[rn].prev;

   985             Node ynode = order_list[yn];

   987             bool rd;

   988             if (!external(xnode, rorder, child_lists, ancestor_map, low_map)) {

   989               rd = true;

   990             } else if (!external(ynode, rorder, child_lists,

   991                                  ancestor_map, low_map)) {

   992               rd = false;

   993             } else if (pertinent(xnode, embed_arc, merge_roots)) {

   994               rd = true;

   995             } else {

   996               rd = false;

   997             }

   999             merge_stack.push_back(std::make_pair(rn, rd));

  1001             pn = rn;

  1002             n = rd ? xn : yn;

  1004           } else if (!external(node, rorder, child_lists,

  1005                                ancestor_map, low_map)) {

  1006             int nn = (node_data[n].next != pn ?

  1007                       node_data[n].next : node_data[n].prev);

  1009             bool nd = n == node_data[nn].prev;

  1011             if (nd) node_data[nn].prev = pn;

  1012             else node_data[nn].next = pn;

  1014             if (n == node_data[pn].prev) node_data[pn].prev = nn;

  1015             else node_data[pn].next = nn;

  1017             node_data[nn].inverted =

  1018               (node_data[nn].prev == node_data[nn].next && nd != rd);

  1020             n = nn;

  1021           }

  1022           else break;

  1024         }

  1026         if (!merge_stack.empty() || n == rn) {

  1027           break;

  1028         }

  1029       }

  1030     }

  1032     void initFace(const Node& node, ArcLists& arc_lists,

  1033                   NodeData& node_data, const PredMap& pred_map,

  1034                   const OrderMap& order_map, const OrderList& order_list) {

  1035       int n = order_map[node];

  1036       int rn = n + order_list.size();

  1038       node_data[n].next = node_data[n].prev = rn;

  1039       node_data[rn].next = node_data[rn].prev = n;

  1041       node_data[n].visited = order_list.size();

  1042       node_data[rn].visited = order_list.size();

  1044       node_data[n].inverted = false;

  1045       node_data[rn].inverted = false;

  1047       Arc arc = pred_map[node];

  1048       Arc rev = _graph.oppositeArc(arc);

  1050       node_data[rn].first = arc;

  1051       node_data[n].first = rev;

  1053       arc_lists[arc].prev = arc;

  1054       arc_lists[arc].next = arc;

  1056       arc_lists[rev].prev = rev;

  1057       arc_lists[rev].next = rev;

  1059     }

  1061     void mergeRemainingFaces(const Node& node, NodeData& node_data,

  1062                              OrderList& order_list, OrderMap& order_map,

  1063                              ChildLists& child_lists, ArcLists& arc_lists) {

  1064       while (child_lists[node].first != INVALID) {

  1065         int dd = order_map[node];

  1066         Node child = child_lists[node].first;

  1067         int cd = order_map[child] + order_list.size();

  1068         child_lists[node].first = child_lists[child].next;

  1070         Arc de = node_data[dd].first;

  1071         Arc ce = node_data[cd].first;

  1073         if (de != INVALID) {

  1074           Arc dne = arc_lists[de].next;

  1075           Arc cne = arc_lists[ce].next;

  1077           arc_lists[de].next = cne;

  1078           arc_lists[ce].next = dne;

  1080           arc_lists[dne].prev = ce;

  1081           arc_lists[cne].prev = de;

  1082         }

  1084         node_data[dd].first = ce;

  1086       }

  1087     }

  1089     void storeEmbedding(const Node& node, NodeData& node_data,

  1090                         OrderMap& order_map, PredMap& pred_map,

  1091                         ArcLists& arc_lists, FlipMap& flip_map) {

  1093       if (node_data[order_map[node]].first == INVALID) return;

  1095       if (pred_map[node] != INVALID) {

  1096         Node source = _graph.source(pred_map[node]);

  1097         flip_map[node] = flip_map[node] != flip_map[source];

  1098       }

  1100       Arc first = node_data[order_map[node]].first;

  1101       Arc prev = first;

  1103       Arc arc = flip_map[node] ?

  1104         arc_lists[prev].prev : arc_lists[prev].next;

  1106       _embedding[prev] = arc;

  1108       while (arc != first) {

  1109         Arc next = arc_lists[arc].prev == prev ?

  1110           arc_lists[arc].next : arc_lists[arc].prev;

  1111         prev = arc; arc = next;

  1112         _embedding[prev] = arc;

  1113       }

  1114     }

  1117     bool external(const Node& node, int rorder,

  1118                   ChildLists& child_lists, AncestorMap& ancestor_map,

  1119                   LowMap& low_map) {

  1120       Node child = child_lists[node].first;

  1122       if (child != INVALID) {

  1123         if (low_map[child] < rorder) return true;

  1124       }

  1126       if (ancestor_map[node] < rorder) return true;

  1128       return false;

  1129     }

  1131     bool pertinent(const Node& node, const EmbedArc& embed_arc,

  1132                    const MergeRoots& merge_roots) {

  1133       return !merge_roots[node].empty() || embed_arc[node] != INVALID;

  1134     }

  1136     int lowPoint(const Node& node, OrderMap& order_map, ChildLists& child_lists,

  1137                  AncestorMap& ancestor_map, LowMap& low_map) {

  1138       int low_point;

  1140       Node child = child_lists[node].first;

  1142       if (child != INVALID) {

  1143         low_point = low_map[child];

  1144       } else {

  1145         low_point = order_map[node];

  1146       }

  1148       if (low_point > ancestor_map[node]) {

  1149         low_point = ancestor_map[node];

  1150       }

  1152       return low_point;

  1153     }

  1155     int findComponentRoot(Node root, Node node, ChildLists& child_lists,

  1156                           OrderMap& order_map, OrderList& order_list) {

  1158       int order = order_map[root];

  1159       int norder = order_map[node];

  1161       Node child = child_lists[root].first;

  1162       while (child != INVALID) {

  1163         int corder = order_map[child];

  1164         if (corder > order && corder < norder) {

  1165           order = corder;

  1166         }

  1167         child = child_lists[child].next;

  1168       }

  1169       return order + order_list.size();

  1170     }

  1172     Node findPertinent(Node node, OrderMap& order_map, NodeData& node_data,

  1173                        EmbedArc& embed_arc, MergeRoots& merge_roots) {

  1174       Node wnode =_graph.target(node_data[order_map[node]].first);

  1175       while (!pertinent(wnode, embed_arc, merge_roots)) {

  1176         wnode = _graph.target(node_data[order_map[wnode]].first);

  1177       }

  1178       return wnode;

  1179     }

  1182     Node findExternal(Node node, int rorder, OrderMap& order_map,

  1183                       ChildLists& child_lists, AncestorMap& ancestor_map,

  1184                       LowMap& low_map, NodeData& node_data) {

  1185       Node wnode =_graph.target(node_data[order_map[node]].first);

  1186       while (!external(wnode, rorder, child_lists, ancestor_map, low_map)) {

  1187         wnode = _graph.target(node_data[order_map[wnode]].first);

  1188       }

  1189       return wnode;

  1190     }

  1192     void markCommonPath(Node node, int rorder, Node& wnode, Node& znode,

  1193                         OrderList& order_list, OrderMap& order_map,

  1194                         NodeData& node_data, ArcLists& arc_lists,

  1195                         EmbedArc& embed_arc, MergeRoots& merge_roots,

  1196                         ChildLists& child_lists, AncestorMap& ancestor_map,

  1197                         LowMap& low_map) {

  1199       Node cnode = node;

  1200       Node pred = INVALID;

  1202       while (true) {

  1204         bool pert = pertinent(cnode, embed_arc, merge_roots);

  1205         bool ext = external(cnode, rorder, child_lists, ancestor_map, low_map);

  1207         if (pert && ext) {

  1208           if (!merge_roots[cnode].empty()) {

  1209             int cn = merge_roots[cnode].back();

  1211             if (low_map[order_list[cn - order_list.size()]] < rorder) {

  1212               Arc arc = node_data[cn].first;

  1213               _kuratowski.set(arc, true);

  1215               pred = cnode;

  1216               cnode = _graph.target(arc);

  1218               continue;

  1219             }

  1220           }

  1221           wnode = znode = cnode;

  1222           return;

  1224         } else if (pert) {

  1225           wnode = cnode;

  1227           while (!external(cnode, rorder, child_lists, ancestor_map, low_map)) {

  1228             Arc arc = node_data[order_map[cnode]].first;

  1230             if (_graph.target(arc) == pred) {

  1231               arc = arc_lists[arc].next;

  1232             }

  1233             _kuratowski.set(arc, true);

  1235             Node next = _graph.target(arc);

  1236             pred = cnode; cnode = next;

  1237           }

  1239           znode = cnode;

  1240           return;

  1242         } else if (ext) {

  1243           znode = cnode;

  1245           while (!pertinent(cnode, embed_arc, merge_roots)) {

  1246             Arc arc = node_data[order_map[cnode]].first;

  1248             if (_graph.target(arc) == pred) {

  1249               arc = arc_lists[arc].next;

  1250             }

  1251             _kuratowski.set(arc, true);

  1253             Node next = _graph.target(arc);

  1254             pred = cnode; cnode = next;

  1255           }

  1257           wnode = cnode;

  1258           return;

  1260         } else {

  1261           Arc arc = node_data[order_map[cnode]].first;

  1263           if (_graph.target(arc) == pred) {

  1264             arc = arc_lists[arc].next;

  1265           }

  1266           _kuratowski.set(arc, true);

  1268           Node next = _graph.target(arc);

  1269           pred = cnode; cnode = next;

  1270         }

  1272       }

  1274     }

  1276     void orientComponent(Node root, int rn, OrderMap& order_map,

  1277                          PredMap& pred_map, NodeData& node_data,

  1278                          ArcLists& arc_lists, FlipMap& flip_map,

  1279                          TypeMap& type_map) {

  1280       node_data[order_map[root]].first = node_data[rn].first;

  1281       type_map[root] = 1;

  1283       std::vector<Node> st, qu;

  1285       st.push_back(root);

  1286       while (!st.empty()) {

  1287         Node node = st.back();

  1288         st.pop_back();

  1289         qu.push_back(node);

  1291         Arc arc = node_data[order_map[node]].first;

  1293         if (type_map[_graph.target(arc)] == 0) {

  1294           st.push_back(_graph.target(arc));

  1295           type_map[_graph.target(arc)] = 1;

  1296         }

  1298         Arc last = arc, pred = arc;

  1299         arc = arc_lists[arc].next;

  1300         while (arc != last) {

  1302           if (type_map[_graph.target(arc)] == 0) {

  1303             st.push_back(_graph.target(arc));

  1304             type_map[_graph.target(arc)] = 1;

  1305           }

  1307           Arc next = arc_lists[arc].next != pred ?

  1308             arc_lists[arc].next : arc_lists[arc].prev;

  1309           pred = arc; arc = next;

  1310         }

  1312       }

  1314       type_map[root] = 2;

  1315       flip_map[root] = false;

  1317       for (int i = 1; i < int(qu.size()); ++i) {

  1319         Node node = qu[i];

  1321         while (type_map[node] != 2) {

  1322           st.push_back(node);

  1323           type_map[node] = 2;

  1324           node = _graph.source(pred_map[node]);

  1325         }

  1327         bool flip = flip_map[node];

  1329         while (!st.empty()) {

  1330           node = st.back();

  1331           st.pop_back();

  1333           flip_map[node] = flip != flip_map[node];

  1334           flip = flip_map[node];

  1336           if (flip) {

  1337             Arc arc = node_data[order_map[node]].first;

  1338             std::swap(arc_lists[arc].prev, arc_lists[arc].next);

  1339             arc = arc_lists[arc].prev;

  1340             std::swap(arc_lists[arc].prev, arc_lists[arc].next);

  1341             node_data[order_map[node]].first = arc;

  1342           }

  1343         }

  1344       }

  1346       for (int i = 0; i < int(qu.size()); ++i) {

  1348         Arc arc = node_data[order_map[qu[i]]].first;

  1349         Arc last = arc, pred = arc;

  1351         arc = arc_lists[arc].next;

  1352         while (arc != last) {

  1354           if (arc_lists[arc].next == pred) {

  1355             std::swap(arc_lists[arc].next, arc_lists[arc].prev);

  1356           }

  1357           pred = arc; arc = arc_lists[arc].next;

  1358         }

  1360       }

  1361     }

  1363     void setFaceFlags(Node root, Node wnode, Node ynode, Node xnode,

  1364                       OrderMap& order_map, NodeData& node_data,

  1365                       TypeMap& type_map) {

  1366       Node node = _graph.target(node_data[order_map[root]].first);

  1368       while (node != ynode) {

  1369         type_map[node] = HIGHY;

  1370         node = _graph.target(node_data[order_map[node]].first);

  1371       }

  1373       while (node != wnode) {

  1374         type_map[node] = LOWY;

  1375         node = _graph.target(node_data[order_map[node]].first);

  1376       }

  1378       node = _graph.target(node_data[order_map[wnode]].first);

  1380       while (node != xnode) {

  1381         type_map[node] = LOWX;

  1382         node = _graph.target(node_data[order_map[node]].first);

  1383       }

  1384       type_map[node] = LOWX;

  1386       node = _graph.target(node_data[order_map[xnode]].first);

  1387       while (node != root) {

  1388         type_map[node] = HIGHX;

  1389         node = _graph.target(node_data[order_map[node]].first);

  1390       }

  1392       type_map[wnode] = PERTINENT;

  1393       type_map[root] = ROOT;

  1394     }

  1396     void findInternalPath(std::vector<Arc>& ipath,

  1397                           Node wnode, Node root, TypeMap& type_map,

  1398                           OrderMap& order_map, NodeData& node_data,

  1399                           ArcLists& arc_lists) {

  1400       std::vector<Arc> st;

  1402       Node node = wnode;

  1404       while (node != root) {

  1405         Arc arc = arc_lists[node_data[order_map[node]].first].next;

  1406         st.push_back(arc);

  1407         node = _graph.target(arc);

  1408       }

  1410       while (true) {

  1411         Arc arc = st.back();

  1412         if (type_map[_graph.target(arc)] == LOWX ||

  1413             type_map[_graph.target(arc)] == HIGHX) {

  1414           break;

  1415         }

  1416         if (type_map[_graph.target(arc)] == 2) {

  1417           type_map[_graph.target(arc)] = 3;

  1419           arc = arc_lists[_graph.oppositeArc(arc)].next;

  1420           st.push_back(arc);

  1421         } else {

  1422           st.pop_back();

  1423           arc = arc_lists[arc].next;

  1425           while (_graph.oppositeArc(arc) == st.back()) {

  1426             arc = st.back();

  1427             st.pop_back();

  1428             arc = arc_lists[arc].next;

  1429           }

  1430           st.push_back(arc);

  1431         }

  1432       }

  1434       for (int i = 0; i < int(st.size()); ++i) {

  1435         if (type_map[_graph.target(st[i])] != LOWY &&

  1436             type_map[_graph.target(st[i])] != HIGHY) {

  1437           for (; i < int(st.size()); ++i) {

  1438             ipath.push_back(st[i]);

  1439           }

  1440         }

  1441       }

  1442     }

  1444     void setInternalFlags(std::vector<Arc>& ipath, TypeMap& type_map) {

  1445       for (int i = 1; i < int(ipath.size()); ++i) {

  1446         type_map[_graph.source(ipath[i])] = INTERNAL;

  1447       }

  1448     }

  1450     void findPilePath(std::vector<Arc>& ppath,

  1451                       Node root, TypeMap& type_map, OrderMap& order_map,

  1452                       NodeData& node_data, ArcLists& arc_lists) {

  1453       std::vector<Arc> st;

  1455       st.push_back(_graph.oppositeArc(node_data[order_map[root]].first));

  1456       st.push_back(node_data[order_map[root]].first);

  1458       while (st.size() > 1) {

  1459         Arc arc = st.back();

  1460         if (type_map[_graph.target(arc)] == INTERNAL) {

  1461           break;

  1462         }

  1463         if (type_map[_graph.target(arc)] == 3) {

  1464           type_map[_graph.target(arc)] = 4;

  1466           arc = arc_lists[_graph.oppositeArc(arc)].next;

  1467           st.push_back(arc);

  1468         } else {

  1469           st.pop_back();

  1470           arc = arc_lists[arc].next;

  1472           while (!st.empty() && _graph.oppositeArc(arc) == st.back()) {

  1473             arc = st.back();

  1474             st.pop_back();

  1475             arc = arc_lists[arc].next;

  1476           }

  1477           st.push_back(arc);

  1478         }

  1479       }

  1481       for (int i = 1; i < int(st.size()); ++i) {

  1482         ppath.push_back(st[i]);

  1483       }

  1484     }

  1487     int markExternalPath(Node node, OrderMap& order_map,

  1488                          ChildLists& child_lists, PredMap& pred_map,

  1489                          AncestorMap& ancestor_map, LowMap& low_map) {

  1490       int lp = lowPoint(node, order_map, child_lists,

  1491                         ancestor_map, low_map);

  1493       if (ancestor_map[node] != lp) {

  1494         node = child_lists[node].first;

  1495         _kuratowski[pred_map[node]] = true;

  1497         while (ancestor_map[node] != lp) {

  1498           for (OutArcIt e(_graph, node); e != INVALID; ++e) {

  1499             Node tnode = _graph.target(e);

  1500             if (order_map[tnode] > order_map[node] && low_map[tnode] == lp) {

  1501               node = tnode;

  1502               _kuratowski[e] = true;

  1503               break;

  1504             }

  1505           }

  1506         }

  1507       }

  1509       for (OutArcIt e(_graph, node); e != INVALID; ++e) {

  1510         if (order_map[_graph.target(e)] == lp) {

  1511           _kuratowski[e] = true;

  1512           break;

  1513         }

  1514       }

  1516       return lp;

  1517     }

  1519     void markPertinentPath(Node node, OrderMap& order_map,

  1520                            NodeData& node_data, ArcLists& arc_lists,

  1521                            EmbedArc& embed_arc, MergeRoots& merge_roots) {

  1522       while (embed_arc[node] == INVALID) {

  1523         int n = merge_roots[node].front();

  1524         Arc arc = node_data[n].first;

  1526         _kuratowski.set(arc, true);

  1528         Node pred = node;

  1529         node = _graph.target(arc);

  1530         while (!pertinent(node, embed_arc, merge_roots)) {

  1531           arc = node_data[order_map[node]].first;

  1532           if (_graph.target(arc) == pred) {

  1533             arc = arc_lists[arc].next;

  1534           }

  1535           _kuratowski.set(arc, true);

  1536           pred = node;

  1537           node = _graph.target(arc);

  1538         }

  1539       }

  1540       _kuratowski.set(embed_arc[node], true);

  1541     }

  1543     void markPredPath(Node node, Node snode, PredMap& pred_map) {

  1544       while (node != snode) {

  1545         _kuratowski.set(pred_map[node], true);

  1546         node = _graph.source(pred_map[node]);

  1547       }

  1548     }

  1550     void markFacePath(Node ynode, Node xnode,

  1551                       OrderMap& order_map, NodeData& node_data) {

  1552       Arc arc = node_data[order_map[ynode]].first;

  1553       Node node = _graph.target(arc);

  1554       _kuratowski.set(arc, true);

  1556       while (node != xnode) {

  1557         arc = node_data[order_map[node]].first;

  1558         _kuratowski.set(arc, true);

  1559         node = _graph.target(arc);

  1560       }

  1561     }

  1563     void markInternalPath(std::vector<Arc>& path) {

  1564       for (int i = 0; i < int(path.size()); ++i) {

  1565         _kuratowski.set(path[i], true);

  1566       }

  1567     }

  1569     void markPilePath(std::vector<Arc>& path) {

  1570       for (int i = 0; i < int(path.size()); ++i) {

  1571         _kuratowski.set(path[i], true);

  1572       }

  1573     }

  1575     void isolateKuratowski(Arc arc, NodeData& node_data,

  1576                            ArcLists& arc_lists, FlipMap& flip_map,

  1577                            OrderMap& order_map, OrderList& order_list,

  1578                            PredMap& pred_map, ChildLists& child_lists,

  1579                            AncestorMap& ancestor_map, LowMap& low_map,

  1580                            EmbedArc& embed_arc, MergeRoots& merge_roots) {

  1582       Node root = _graph.source(arc);

  1583       Node enode = _graph.target(arc);

  1585       int rorder = order_map[root];

  1587       TypeMap type_map(_graph, 0);

  1589       int rn = findComponentRoot(root, enode, child_lists,

  1590                                  order_map, order_list);

  1592       Node xnode = order_list[node_data[rn].next];

  1593       Node ynode = order_list[node_data[rn].prev];

  1595       // Minor-A

  1596       {

  1597         while (!merge_roots[xnode].empty() || !merge_roots[ynode].empty()) {

  1599           if (!merge_roots[xnode].empty()) {

  1600             root = xnode;

  1601             rn = merge_roots[xnode].front();

  1602           } else {

  1603             root = ynode;

  1604             rn = merge_roots[ynode].front();

  1605           }

  1607           xnode = order_list[node_data[rn].next];

  1608           ynode = order_list[node_data[rn].prev];

  1609         }

  1611         if (root != _graph.source(arc)) {

  1612           orientComponent(root, rn, order_map, pred_map,

  1613                           node_data, arc_lists, flip_map, type_map);

  1614           markFacePath(root, root, order_map, node_data);

  1615           int xlp = markExternalPath(xnode, order_map, child_lists,

  1616                                      pred_map, ancestor_map, low_map);

  1617           int ylp = markExternalPath(ynode, order_map, child_lists,

  1618                                      pred_map, ancestor_map, low_map);

  1619           markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

  1620           Node lwnode = findPertinent(ynode, order_map, node_data,

  1621                                       embed_arc, merge_roots);

  1623           markPertinentPath(lwnode, order_map, node_data, arc_lists,

  1624                             embed_arc, merge_roots);

  1626           return;

  1627         }

  1628       }

  1630       orientComponent(root, rn, order_map, pred_map,

  1631                       node_data, arc_lists, flip_map, type_map);

  1633       Node wnode = findPertinent(ynode, order_map, node_data,

  1634                                  embed_arc, merge_roots);

  1635       setFaceFlags(root, wnode, ynode, xnode, order_map, node_data, type_map);

  1638       //Minor-B

  1639       if (!merge_roots[wnode].empty()) {

  1640         int cn = merge_roots[wnode].back();

  1641         Node rep = order_list[cn - order_list.size()];

  1642         if (low_map[rep] < rorder) {

  1643           markFacePath(root, root, order_map, node_data);

  1644           int xlp = markExternalPath(xnode, order_map, child_lists,

  1645                                      pred_map, ancestor_map, low_map);

  1646           int ylp = markExternalPath(ynode, order_map, child_lists,

  1647                                      pred_map, ancestor_map, low_map);

  1649           Node lwnode, lznode;

  1650           markCommonPath(wnode, rorder, lwnode, lznode, order_list,

  1651                          order_map, node_data, arc_lists, embed_arc,

  1652                          merge_roots, child_lists, ancestor_map, low_map);

  1654           markPertinentPath(lwnode, order_map, node_data, arc_lists,

  1655                             embed_arc, merge_roots);

  1656           int zlp = markExternalPath(lznode, order_map, child_lists,

  1657                                      pred_map, ancestor_map, low_map);

  1659           int minlp = xlp < ylp ? xlp : ylp;

  1660           if (zlp < minlp) minlp = zlp;

  1662           int maxlp = xlp > ylp ? xlp : ylp;

  1663           if (zlp > maxlp) maxlp = zlp;

  1665           markPredPath(order_list[maxlp], order_list[minlp], pred_map);

  1667           return;

  1668         }

  1669       }

  1671       Node pxnode, pynode;

  1672       std::vector<Arc> ipath;

  1673       findInternalPath(ipath, wnode, root, type_map, order_map,

  1674                        node_data, arc_lists);

  1675       setInternalFlags(ipath, type_map);

  1676       pynode = _graph.source(ipath.front());

  1677       pxnode = _graph.target(ipath.back());

  1679       wnode = findPertinent(pynode, order_map, node_data,

  1680                             embed_arc, merge_roots);

  1682       // Minor-C

  1683       {

  1684         if (type_map[_graph.source(ipath.front())] == HIGHY) {

  1685           if (type_map[_graph.target(ipath.back())] == HIGHX) {

  1686             markFacePath(xnode, pxnode, order_map, node_data);

  1687           }

  1688           markFacePath(root, xnode, order_map, node_data);

  1689           markPertinentPath(wnode, order_map, node_data, arc_lists,

  1690                             embed_arc, merge_roots);

  1691           markInternalPath(ipath);

  1692           int xlp = markExternalPath(xnode, order_map, child_lists,

  1693                                      pred_map, ancestor_map, low_map);

  1694           int ylp = markExternalPath(ynode, order_map, child_lists,

  1695                                      pred_map, ancestor_map, low_map);

  1696           markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

  1697           return;

  1698         }

  1700         if (type_map[_graph.target(ipath.back())] == HIGHX) {

  1701           markFacePath(ynode, root, order_map, node_data);

  1702           markPertinentPath(wnode, order_map, node_data, arc_lists,

  1703                             embed_arc, merge_roots);

  1704           markInternalPath(ipath);

  1705           int xlp = markExternalPath(xnode, order_map, child_lists,

  1706                                      pred_map, ancestor_map, low_map);

  1707           int ylp = markExternalPath(ynode, order_map, child_lists,

  1708                                      pred_map, ancestor_map, low_map);

  1709           markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

  1710           return;

  1711         }

  1712       }

  1714       std::vector<Arc> ppath;

  1715       findPilePath(ppath, root, type_map, order_map, node_data, arc_lists);

  1717       // Minor-D

  1718       if (!ppath.empty()) {

  1719         markFacePath(ynode, xnode, order_map, node_data);

  1720         markPertinentPath(wnode, order_map, node_data, arc_lists,

  1721                           embed_arc, merge_roots);

  1722         markPilePath(ppath);

  1723         markInternalPath(ipath);

  1724         int xlp = markExternalPath(xnode, order_map, child_lists,

  1725                                    pred_map, ancestor_map, low_map);

  1726         int ylp = markExternalPath(ynode, order_map, child_lists,

  1727                                    pred_map, ancestor_map, low_map);

  1728         markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

  1729         return;

  1730       }

  1732       // Minor-E*

  1733       {

  1735         if (!external(wnode, rorder, child_lists, ancestor_map, low_map)) {

  1736           Node znode = findExternal(pynode, rorder, order_map,

  1737                                     child_lists, ancestor_map,

  1738                                     low_map, node_data);

  1740           if (type_map[znode] == LOWY) {

  1741             markFacePath(root, xnode, order_map, node_data);

  1742             markPertinentPath(wnode, order_map, node_data, arc_lists,

  1743                               embed_arc, merge_roots);

  1744             markInternalPath(ipath);

  1745             int xlp = markExternalPath(xnode, order_map, child_lists,

  1746                                        pred_map, ancestor_map, low_map);

  1747             int zlp = markExternalPath(znode, order_map, child_lists,

  1748                                        pred_map, ancestor_map, low_map);

  1749             markPredPath(root, order_list[xlp < zlp ? xlp : zlp], pred_map);

  1750           } else {

  1751             markFacePath(ynode, root, order_map, node_data);

  1752             markPertinentPath(wnode, order_map, node_data, arc_lists,

  1753                               embed_arc, merge_roots);

  1754             markInternalPath(ipath);

  1755             int ylp = markExternalPath(ynode, order_map, child_lists,

  1756                                        pred_map, ancestor_map, low_map);

  1757             int zlp = markExternalPath(znode, order_map, child_lists,

  1758                                        pred_map, ancestor_map, low_map);

  1759             markPredPath(root, order_list[ylp < zlp ? ylp : zlp], pred_map);

  1760           }

  1761           return;

  1762         }

  1764         int xlp = markExternalPath(xnode, order_map, child_lists,

  1765                                    pred_map, ancestor_map, low_map);

  1766         int ylp = markExternalPath(ynode, order_map, child_lists,

  1767                                    pred_map, ancestor_map, low_map);

  1768         int wlp = markExternalPath(wnode, order_map, child_lists,

  1769                                    pred_map, ancestor_map, low_map);

  1771         if (wlp > xlp && wlp > ylp) {

  1772           markFacePath(root, root, order_map, node_data);

  1773           markPredPath(root, order_list[xlp < ylp ? xlp : ylp], pred_map);

  1774           return;

  1775         }

  1777         markInternalPath(ipath);

  1778         markPertinentPath(wnode, order_map, node_data, arc_lists,

  1779                           embed_arc, merge_roots);

  1781         if (xlp > ylp && xlp > wlp) {

  1782           markFacePath(root, pynode, order_map, node_data);

  1783           markFacePath(wnode, xnode, order_map, node_data);

  1784           markPredPath(root, order_list[ylp < wlp ? ylp : wlp], pred_map);

  1785           return;

  1786         }

  1788         if (ylp > xlp && ylp > wlp) {

  1789           markFacePath(pxnode, root, order_map, node_data);

  1790           markFacePath(ynode, wnode, order_map, node_data);

  1791           markPredPath(root, order_list[xlp < wlp ? xlp : wlp], pred_map);

  1792           return;

  1793         }

  1795         if (pynode != ynode) {

  1796           markFacePath(pxnode, wnode, order_map, node_data);

  1798           int minlp = xlp < ylp ? xlp : ylp;

  1799           if (wlp < minlp) minlp = wlp;

  1801           int maxlp = xlp > ylp ? xlp : ylp;

  1802           if (wlp > maxlp) maxlp = wlp;

  1804           markPredPath(order_list[maxlp], order_list[minlp], pred_map);

  1805           return;

  1806         }

  1808         if (pxnode != xnode) {

  1809           markFacePath(wnode, pynode, order_map, node_data);

  1811           int minlp = xlp < ylp ? xlp : ylp;

  1812           if (wlp < minlp) minlp = wlp;

  1814           int maxlp = xlp > ylp ? xlp : ylp;

  1815           if (wlp > maxlp) maxlp = wlp;

  1817           markPredPath(order_list[maxlp], order_list[minlp], pred_map);

  1818           return;

  1819         }

  1821         markFacePath(root, root, order_map, node_data);

  1822         int minlp = xlp < ylp ? xlp : ylp;

  1823         if (wlp < minlp) minlp = wlp;

  1824         markPredPath(root, order_list[minlp], pred_map);

  1825         return;

  1826       }

  1828     }

  1830   };

  1832   namespace _planarity_bits {

  1834     template <typename Graph, typename EmbeddingMap>

  1835     void makeConnected(Graph& graph, EmbeddingMap& embedding) {

  1836       DfsVisitor<Graph> null_visitor;

  1837       DfsVisit<Graph, DfsVisitor<Graph> > dfs(graph, null_visitor);

  1838       dfs.init();

  1840       typename Graph::Node u = INVALID;

  1841       for (typename Graph::NodeIt n(graph); n != INVALID; ++n) {

  1842         if (!dfs.reached(n)) {

  1843           dfs.addSource(n);

  1844           dfs.start();

  1845           if (u == INVALID) {

  1846             u = n;

  1847           } else {

  1848             typename Graph::Node v = n;

  1850             typename Graph::Arc ue = typename Graph::OutArcIt(graph, u);

  1851             typename Graph::Arc ve = typename Graph::OutArcIt(graph, v);

  1853             typename Graph::Arc e = graph.direct(graph.addEdge(u, v), true);

  1855             if (ue != INVALID) {

  1856               embedding[e] = embedding[ue];

  1857               embedding[ue] = e;

  1858             } else {

  1859               embedding[e] = e;

  1860             }

  1862             if (ve != INVALID) {

  1863               embedding[graph.oppositeArc(e)] = embedding[ve];

  1864               embedding[ve] = graph.oppositeArc(e);

  1865             } else {

  1866               embedding[graph.oppositeArc(e)] = graph.oppositeArc(e);

  1867             }

  1868           }

  1869         }

  1870       }

  1871     }

  1873     template <typename Graph, typename EmbeddingMap>

  1874     void makeBiNodeConnected(Graph& graph, EmbeddingMap& embedding) {

  1875       typename Graph::template ArcMap<bool> processed(graph);

  1877       std::vector<typename Graph::Arc> arcs;

  1878       for (typename Graph::ArcIt e(graph); e != INVALID; ++e) {

  1879         arcs.push_back(e);

  1880       }

  1882       IterableBoolMap<Graph, typename Graph::Node> visited(graph, false);

  1884       for (int i = 0; i < int(arcs.size()); ++i) {

  1885         typename Graph::Arc pp = arcs[i];

  1886         if (processed[pp]) continue;

  1888         typename Graph::Arc e = embedding[graph.oppositeArc(pp)];

  1889         processed[e] = true;

  1890         visited.set(graph.source(e), true);

  1892         typename Graph::Arc p = e, l = e;

  1893         e = embedding[graph.oppositeArc(e)];

  1895         while (e != l) {

  1896           processed[e] = true;

  1898           if (visited[graph.source(e)]) {

  1900             typename Graph::Arc n =

  1901               graph.direct(graph.addEdge(graph.source(p),

  1902                                            graph.target(e)), true);

  1903             embedding[n] = p;

  1904             embedding[graph.oppositeArc(pp)] = n;

  1906             embedding[graph.oppositeArc(n)] =

  1907               embedding[graph.oppositeArc(e)];

  1908             embedding[graph.oppositeArc(e)] =

  1909               graph.oppositeArc(n);

  1911             p = n;

  1912             e = embedding[graph.oppositeArc(n)];

  1913           } else {

  1914             visited.set(graph.source(e), true);

  1915             pp = p;

  1916             p = e;

  1917             e = embedding[graph.oppositeArc(e)];

  1918           }

  1919         }

  1920         visited.setAll(false);

  1921       }

  1922     }

  1925     template <typename Graph, typename EmbeddingMap>

  1926     void makeMaxPlanar(Graph& graph, EmbeddingMap& embedding) {

  1928       typename Graph::template NodeMap<int> degree(graph);

  1930       for (typename Graph::NodeIt n(graph); n != INVALID; ++n) {

  1931         degree[n] = countIncEdges(graph, n);

  1932       }

  1934       typename Graph::template ArcMap<bool> processed(graph);

  1935       IterableBoolMap<Graph, typename Graph::Node> visited(graph, false);

  1937       std::vector<typename Graph::Arc> arcs;

  1938       for (typename Graph::ArcIt e(graph); e != INVALID; ++e) {

  1939         arcs.push_back(e);

  1940       }

  1942       for (int i = 0; i < int(arcs.size()); ++i) {

  1943         typename Graph::Arc e = arcs[i];

  1945         if (processed[e]) continue;

  1946         processed[e] = true;

  1948         typename Graph::Arc mine = e;

  1949         int mind = degree[graph.source(e)];

  1951         int face_size = 1;

  1953         typename Graph::Arc l = e;

  1954         e = embedding[graph.oppositeArc(e)];

  1955         while (l != e) {

  1956           processed[e] = true;

  1958           ++face_size;

  1960           if (degree[graph.source(e)] < mind) {

  1961             mine = e;

  1962             mind = degree[graph.source(e)];

  1963           }

  1965           e = embedding[graph.oppositeArc(e)];

  1966         }

  1968         if (face_size < 4) {

  1969           continue;

  1970         }

  1972         typename Graph::Node s = graph.source(mine);

  1973         for (typename Graph::OutArcIt e(graph, s); e != INVALID; ++e) {

  1974           visited.set(graph.target(e), true);

  1975         }

  1977         typename Graph::Arc oppe = INVALID;

  1979         e = embedding[graph.oppositeArc(mine)];

  1980         e = embedding[graph.oppositeArc(e)];

  1981         while (graph.target(e) != s) {

  1982           if (visited[graph.source(e)]) {

  1983             oppe = e;

  1984             break;

  1985           }

  1986           e = embedding[graph.oppositeArc(e)];

  1987         }

  1988         visited.setAll(false);

  1990         if (oppe == INVALID) {

  1992           e = embedding[graph.oppositeArc(mine)];

  1993           typename Graph::Arc pn = mine, p = e;

  1995           e = embedding[graph.oppositeArc(e)];

  1996           while (graph.target(e) != s) {

  1997             typename Graph::Arc n =

  1998               graph.direct(graph.addEdge(s, graph.source(e)), true);

  2000             embedding[n] = pn;

  2001             embedding[graph.oppositeArc(n)] = e;

  2002             embedding[graph.oppositeArc(p)] = graph.oppositeArc(n);

  2004             pn = n;

  2006             p = e;

  2007             e = embedding[graph.oppositeArc(e)];

  2008           }

  2010           embedding[graph.oppositeArc(e)] = pn;

  2012         } else {

  2014           mine = embedding[graph.oppositeArc(mine)];

  2015           s = graph.source(mine);

  2016           oppe = embedding[graph.oppositeArc(oppe)];

  2017           typename Graph::Node t = graph.source(oppe);

  2019           typename Graph::Arc ce = graph.direct(graph.addEdge(s, t), true);

  2020           embedding[ce] = mine;

  2021           embedding[graph.oppositeArc(ce)] = oppe;

  2023           typename Graph::Arc pn = ce, p = oppe;

  2024           e = embedding[graph.oppositeArc(oppe)];

  2025           while (graph.target(e) != s) {

  2026             typename Graph::Arc n =

  2027               graph.direct(graph.addEdge(s, graph.source(e)), true);

  2029             embedding[n] = pn;

  2030             embedding[graph.oppositeArc(n)] = e;

  2031             embedding[graph.oppositeArc(p)] = graph.oppositeArc(n);

  2033             pn = n;

  2035             p = e;

  2036             e = embedding[graph.oppositeArc(e)];

  2038           }

  2039           embedding[graph.oppositeArc(e)] = pn;

  2041           pn = graph.oppositeArc(ce), p = mine;

  2042           e = embedding[graph.oppositeArc(mine)];

  2043           while (graph.target(e) != t) {

  2044             typename Graph::Arc n =

  2045               graph.direct(graph.addEdge(t, graph.source(e)), true);

  2047             embedding[n] = pn;

  2048             embedding[graph.oppositeArc(n)] = e;

  2049             embedding[graph.oppositeArc(p)] = graph.oppositeArc(n);

  2051             pn = n;

  2053             p = e;

  2054             e = embedding[graph.oppositeArc(e)];

  2056           }

  2057           embedding[graph.oppositeArc(e)] = pn;

  2058         }

  2059       }

  2060     }

  2062   }

  2064   /// \ingroup planar

  2065   ///

  2066   /// \brief Schnyder's planar drawing algorithm

  2067   ///

  2068   /// The planar drawing algorithm calculates positions for the nodes

  2069   /// in the plane. These coordinates satisfy that if the edges are

  2070   /// represented with straight lines, then they will not intersect

  2071   /// each other.

  2072   ///

  2073   /// Scnyder's algorithm embeds the graph on an \c (n-2)x(n-2) size grid,

  2074   /// i.e. each node will be located in the \c [0..n-2]x[0..n-2] square.

  2075   /// The time complexity of the algorithm is O(n).

  2076   ///

  2077   /// \see PlanarEmbedding

  2078   template <typename Graph>

  2079   class PlanarDrawing {

  2080   public:

  2082     TEMPLATE_GRAPH_TYPEDEFS(Graph);

  2084     /// \brief The point type for storing coordinates

  2085     typedef dim2::Point<int> Point;

  2086     /// \brief The map type for storing the coordinates of the nodes

  2087     typedef typename Graph::template NodeMap<Point> PointMap;

  2090     /// \brief Constructor

  2091     ///

  2092     /// Constructor

  2093     /// \pre The graph must be simple, i.e. it should not

  2094     /// contain parallel or loop arcs.

  2095     PlanarDrawing(const Graph& graph)

  2096       : _graph(graph), _point_map(graph) {}

  2098   private:

  2100     template <typename AuxGraph, typename AuxEmbeddingMap>

  2101     void drawing(const AuxGraph& graph,

  2102                  const AuxEmbeddingMap& next,

  2103                  PointMap& point_map) {

  2104       TEMPLATE_GRAPH_TYPEDEFS(AuxGraph);

  2106       typename AuxGraph::template ArcMap<Arc> prev(graph);

  2108       for (NodeIt n(graph); n != INVALID; ++n) {

  2109         Arc e = OutArcIt(graph, n);

  2111         Arc p = e, l = e;

  2113         e = next[e];

  2114         while (e != l) {

  2115           prev[e] = p;

  2116           p = e;

  2117           e = next[e];

  2118         }

  2119         prev[e] = p;

  2120       }

  2122       Node anode, bnode, cnode;

  2124       {

  2125         Arc e = ArcIt(graph);

  2126         anode = graph.source(e);

  2127         bnode = graph.target(e);

  2128         cnode = graph.target(next[graph.oppositeArc(e)]);

  2129       }

  2131       IterableBoolMap<AuxGraph, Node> proper(graph, false);

  2132       typename AuxGraph::template NodeMap<int> conn(graph, -1);

  2134       conn[anode] = conn[bnode] = -2;

  2135       {

  2136         for (OutArcIt e(graph, anode); e != INVALID; ++e) {

  2137           Node m = graph.target(e);

  2138           if (conn[m] == -1) {

  2139             conn[m] = 1;

  2140           }

  2141         }

  2142         conn[cnode] = 2;

  2144         for (OutArcIt e(graph, bnode); e != INVALID; ++e) {

  2145           Node m = graph.target(e);

  2146           if (conn[m] == -1) {

  2147             conn[m] = 1;

  2148           } else if (conn[m] != -2) {

  2149             conn[m] += 1;

  2150             Arc pe = graph.oppositeArc(e);

  2151             if (conn[graph.target(next[pe])] == -2) {

  2152               conn[m] -= 1;

  2153             }

  2154             if (conn[graph.target(prev[pe])] == -2) {

  2155               conn[m] -= 1;

  2156             }

  2158             proper.set(m, conn[m] == 1);

  2159           }

  2160         }

  2161       }

  2164       typename AuxGraph::template ArcMap<int> angle(graph, -1);

  2166       while (proper.trueNum() != 0) {

  2167         Node n = typename IterableBoolMap<AuxGraph, Node>::TrueIt(proper);

  2168         proper.set(n, false);

  2169         conn[n] = -2;

  2171         for (OutArcIt e(graph, n); e != INVALID; ++e) {

  2172           Node m = graph.target(e);

  2173           if (conn[m] == -1) {

  2174             conn[m] = 1;

  2175           } else if (conn[m] != -2) {

  2176             conn[m] += 1;

  2177             Arc pe = graph.oppositeArc(e);

  2178             if (conn[graph.target(next[pe])] == -2) {

  2179               conn[m] -= 1;

  2180             }

  2181             if (conn[graph.target(prev[pe])] == -2) {

  2182               conn[m] -= 1;

  2183             }

  2185             proper.set(m, conn[m] == 1);

  2186           }

  2187         }

  2189         {

  2190           Arc e = OutArcIt(graph, n);

  2191           Arc p = e, l = e;

  2193           e = next[e];

  2194           while (e != l) {

  2196             if (conn[graph.target(e)] == -2 && conn[graph.target(p)] == -2) {

  2197               Arc f = e;

  2198               angle[f] = 0;

  2199               f = next[graph.oppositeArc(f)];

  2200               angle[f] = 1;

  2201               f = next[graph.oppositeArc(f)];

  2202               angle[f] = 2;

  2203             }

  2205             p = e;

  2206             e = next[e];

  2207           }

  2209           if (conn[graph.target(e)] == -2 && conn[graph.target(p)] == -2) {

  2210             Arc f = e;

  2211             angle[f] = 0;

  2212             f = next[graph.oppositeArc(f)];

  2213             angle[f] = 1;

  2214             f = next[graph.oppositeArc(f)];

  2215             angle[f] = 2;

  2216           }

  2217         }

  2218       }

  2220       typename AuxGraph::template NodeMap<Node> apred(graph, INVALID);

  2221       typename AuxGraph::template NodeMap<Node> bpred(graph, INVALID);

  2222       typename AuxGraph::template NodeMap<Node> cpred(graph, INVALID);

  2224       typename AuxGraph::template NodeMap<int> apredid(graph, -1);

  2225       typename AuxGraph::template NodeMap<int> bpredid(graph, -1);

  2226       typename AuxGraph::template NodeMap<int> cpredid(graph, -1);

  2228       for (ArcIt e(graph); e != INVALID; ++e) {

  2229         if (angle[e] == angle[next[e]]) {

  2230           switch (angle[e]) {

  2231           case 2:

  2232             apred[graph.target(e)] = graph.source(e);

  2233             apredid[graph.target(e)] = graph.id(graph.source(e));

  2234             break;

  2235           case 1:

  2236             bpred[graph.target(e)] = graph.source(e);

  2237             bpredid[graph.target(e)] = graph.id(graph.source(e));

  2238             break;

  2239           case 0:

  2240             cpred[graph.target(e)] = graph.source(e);

  2241             cpredid[graph.target(e)] = graph.id(graph.source(e));

  2242             break;

  2243           }

  2244         }

  2245       }

  2247       cpred[anode] = INVALID;

  2248       cpred[bnode] = INVALID;

  2250       std::vector<Node> aorder, border, corder;

  2252       {

  2253         typename AuxGraph::template NodeMap<bool> processed(graph, false);

  2254         std::vector<Node> st;

  2255         for (NodeIt n(graph); n != INVALID; ++n) {

  2256           if (!processed[n] && n != bnode && n != cnode) {

  2257             st.push_back(n);

  2258             processed[n] = true;

  2259             Node m = apred[n];

  2260             while (m != INVALID && !processed[m]) {

  2261               st.push_back(m);

  2262               processed[m] = true;

  2263               m = apred[m];

  2264             }

  2265             while (!st.empty()) {

  2266               aorder.push_back(st.back());

  2267               st.pop_back();

  2268             }

  2269           }

  2270         }

  2271       }

  2273       {

  2274         typename AuxGraph::template NodeMap<bool> processed(graph, false);

  2275         std::vector<Node> st;

  2276         for (NodeIt n(graph); n != INVALID; ++n) {

  2277           if (!processed[n] && n != cnode && n != anode) {

  2278             st.push_back(n);

  2279             processed[n] = true;

  2280             Node m = bpred[n];

  2281             while (m != INVALID && !processed[m]) {

  2282               st.push_back(m);

  2283               processed[m] = true;

  2284               m = bpred[m];

  2285             }

  2286             while (!st.empty()) {

  2287               border.push_back(st.back());

  2288               st.pop_back();

  2289             }

  2290           }

  2291         }

  2292       }

  2294       {

  2295         typename AuxGraph::template NodeMap<bool> processed(graph, false);

  2296         std::vector<Node> st;

  2297         for (NodeIt n(graph); n != INVALID; ++n) {

  2298           if (!processed[n] && n != anode && n != bnode) {

  2299             st.push_back(n);

  2300             processed[n] = true;

  2301             Node m = cpred[n];

  2302             while (m != INVALID && !processed[m]) {

  2303               st.push_back(m);

  2304               processed[m] = true;

  2305               m = cpred[m];

  2306             }

  2307             while (!st.empty()) {

  2308               corder.push_back(st.back());

  2309               st.pop_back();

  2310             }

  2311           }

  2312         }

  2313       }

  2315       typename AuxGraph::template NodeMap<int> atree(graph, 0);

  2316       for (int i = aorder.size() - 1; i >= 0; --i) {

  2317         Node n = aorder[i];

  2318         atree[n] = 1;

  2319         for (OutArcIt e(graph, n); e != INVALID; ++e) {

  2320           if (apred[graph.target(e)] == n) {

  2321             atree[n] += atree[graph.target(e)];

  2322           }

  2323         }

  2324       }

  2326       typename AuxGraph::template NodeMap<int> btree(graph, 0);

  2327       for (int i = border.size() - 1; i >= 0; --i) {

  2328         Node n = border[i];

  2329         btree[n] = 1;

  2330         for (OutArcIt e(graph, n); e != INVALID; ++e) {

  2331           if (bpred[graph.target(e)] == n) {

  2332             btree[n] += btree[graph.target(e)];

  2333           }

  2334         }

  2335       }

  2337       typename AuxGraph::template NodeMap<int> apath(graph, 0);

  2338       apath[bnode] = apath[cnode] = 1;

  2339       typename AuxGraph::template NodeMap<int> apath_btree(graph, 0);

  2340       apath_btree[bnode] = btree[bnode];

  2341       for (int i = 1; i < int(aorder.size()); ++i) {

  2342         Node n = aorder[i];

  2343         apath[n] = apath[apred[n]] + 1;

  2344         apath_btree[n] = btree[n] + apath_btree[apred[n]];

  2345       }

  2347       typename AuxGraph::template NodeMap<int> bpath_atree(graph, 0);

  2348       bpath_atree[anode] = atree[anode];

  2349       for (int i = 1; i < int(border.size()); ++i) {

  2350         Node n = border[i];

  2351         bpath_atree[n] = atree[n] + bpath_atree[bpred[n]];

  2352       }

  2354       typename AuxGraph::template NodeMap<int> cpath(graph, 0);

  2355       cpath[anode] = cpath[bnode] = 1;

  2356       typename AuxGraph::template NodeMap<int> cpath_atree(graph, 0);

  2357       cpath_atree[anode] = atree[anode];

  2358       typename AuxGraph::template NodeMap<int> cpath_btree(graph, 0);

  2359       cpath_btree[bnode] = btree[bnode];

  2360       for (int i = 1; i < int(corder.size()); ++i) {

  2361         Node n = corder[i];

  2362         cpath[n] = cpath[cpred[n]] + 1;

  2363         cpath_atree[n] = atree[n] + cpath_atree[cpred[n]];

  2364         cpath_btree[n] = btree[n] + cpath_btree[cpred[n]];

  2365       }

  2367       typename AuxGraph::template NodeMap<int> third(graph);

  2368       for (NodeIt n(graph); n != INVALID; ++n) {

  2369         point_map[n].x =

  2370           bpath_atree[n] + cpath_atree[n] - atree[n] - cpath[n] + 1;

  2371         point_map[n].y =

  2372           cpath_btree[n] + apath_btree[n] - btree[n] - apath[n] + 1;

  2373       }

  2375     }

  2377   public:

  2379     /// \brief Calculate the node positions

  2380     ///

  2381     /// This function calculates the node positions on the plane.

  2382     /// \return \c true if the graph is planar.

  2383     bool run() {

  2384       PlanarEmbedding<Graph> pe(_graph);

  2385       if (!pe.run()) return false;

  2387       run(pe);

  2388       return true;

  2389     }

  2391     /// \brief Calculate the node positions according to a

  2392     /// combinatorical embedding

  2393     ///

  2394     /// This function calculates the node positions on the plane.

  2395     /// The given \c embedding map should contain a valid combinatorical

  2396     /// embedding, i.e. a valid cyclic order of the arcs.

  2397     /// It can be computed using PlanarEmbedding.

  2398     template <typename EmbeddingMap>

  2399     void run(const EmbeddingMap& embedding) {

  2400       typedef SmartEdgeSet<Graph> AuxGraph;

  2402       if (3 * countNodes(_graph) - 6 == countEdges(_graph)) {

  2403         drawing(_graph, embedding, _point_map);

  2404         return;

  2405       }

  2407       AuxGraph aux_graph(_graph);

  2408       typename AuxGraph::template ArcMap<typename AuxGraph::Arc>

  2409         aux_embedding(aux_graph);

  2411       {

  2413         typename Graph::template EdgeMap<typename AuxGraph::Edge>

  2414           ref(_graph);

  2416         for (EdgeIt e(_graph); e != INVALID; ++e) {

  2417           ref[e] = aux_graph.addEdge(_graph.u(e), _graph.v(e));

  2418         }

  2420         for (EdgeIt e(_graph); e != INVALID; ++e) {

  2421           Arc ee = embedding[_graph.direct(e, true)];

  2422           aux_embedding[aux_graph.direct(ref[e], true)] =

  2423             aux_graph.direct(ref[ee], _graph.direction(ee));

  2424           ee = embedding[_graph.direct(e, false)];

  2425           aux_embedding[aux_graph.direct(ref[e], false)] =

  2426             aux_graph.direct(ref[ee], _graph.direction(ee));

  2427         }

  2428       }

  2429       _planarity_bits::makeConnected(aux_graph, aux_embedding);

  2430       _planarity_bits::makeBiNodeConnected(aux_graph, aux_embedding);

  2431       _planarity_bits::makeMaxPlanar(aux_graph, aux_embedding);

  2432       drawing(aux_graph, aux_embedding, _point_map);

  2433     }

  2435     /// \brief The coordinate of the given node

  2436     ///

  2437     /// This function returns the coordinate of the given node.

  2438     Point operator[](const Node& node) const {

  2439       return _point_map[node];

  2440     }

  2442     /// \brief Return the grid embedding in a node map

  2443     ///

  2444     /// This function returns the grid embedding in a node map of

  2445     /// \c dim2::Point<int> coordinates.

  2446     const PointMap& coords() const {

  2447       return _point_map;

  2448     }

  2450   private:

  2452     const Graph& _graph;

  2453     PointMap _point_map;

  2455   };

  2457   namespace _planarity_bits {

  2459     template <typename ColorMap>

  2460     class KempeFilter {

  2461     public:

  2462       typedef typename ColorMap::Key Key;

  2463       typedef bool Value;

  2465       KempeFilter(const ColorMap& color_map,

  2466                   const typename ColorMap::Value& first,

  2467                   const typename ColorMap::Value& second)

  2468         : _color_map(color_map), _first(first), _second(second) {}

  2470       Value operator[](const Key& key) const {

  2471         return _color_map[key] == _first || _color_map[key] == _second;

  2472       }

  2474     private:

  2475       const ColorMap& _color_map;

  2476       typename ColorMap::Value _first, _second;

  2477     };

  2478   }

  2480   /// \ingroup planar

  2481   ///

  2482   /// \brief Coloring planar graphs

  2483   ///

  2484   /// The graph coloring problem is the coloring of the graph nodes

  2485   /// so that there are no adjacent nodes with the same color. The

  2486   /// planar graphs can always be colored with four colors, which is

  2487   /// proved by Appel and Haken. Their proofs provide a quadratic

  2488   /// time algorithm for four coloring, but it could not be used to

  2489   /// implement an efficient algorithm. The five and six coloring can be

  2490   /// made in linear time, but in this class, the five coloring has

  2491   /// quadratic worst case time complexity. The two coloring (if

  2492   /// possible) is solvable with a graph search algorithm and it is

  2493   /// implemented in \ref bipartitePartitions() function in LEMON. To

  2494   /// decide whether a planar graph is three colorable is NP-complete.

  2495   ///

  2496   /// This class contains member functions for calculate colorings

  2497   /// with five and six colors. The six coloring algorithm is a simple

  2498   /// greedy coloring on the backward minimum outgoing order of nodes.

  2499   /// This order can be computed by selecting the node with least

  2500   /// outgoing arcs to unprocessed nodes in each phase. This order

  2501   /// guarantees that when a node is chosen for coloring it has at

  2502   /// most five already colored adjacents. The five coloring algorithm

  2503   /// use the same method, but if the greedy approach fails to color

  2504   /// with five colors, i.e. the node has five already different

  2505   /// colored neighbours, it swaps the colors in one of the connected

  2506   /// two colored sets with the Kempe recoloring method.

  2507   template <typename Graph>

  2508   class PlanarColoring {

  2509   public:

  2511     TEMPLATE_GRAPH_TYPEDEFS(Graph);

  2513     /// \brief The map type for storing color indices

  2514     typedef typename Graph::template NodeMap<int> IndexMap;

  2515     /// \brief The map type for storing colors

  2516     ///

  2517     /// The map type for storing colors.

  2518     /// \see Palette, Color

  2519     typedef ComposeMap<Palette, IndexMap> ColorMap;

  2521     /// \brief Constructor

  2522     ///

  2523     /// Constructor.

  2524     /// \pre The graph must be simple, i.e. it should not

  2525     /// contain parallel or loop arcs.

  2526     PlanarColoring(const Graph& graph)

  2527       : _graph(graph), _color_map(graph), _palette(0) {

  2528       _palette.add(Color(1,0,0));

  2529       _palette.add(Color(0,1,0));

  2530       _palette.add(Color(0,0,1));

  2531       _palette.add(Color(1,1,0));

  2532       _palette.add(Color(1,0,1));

  2533       _palette.add(Color(0,1,1));

  2534     }

  2536     /// \brief Return the node map of color indices

  2537     ///

  2538     /// This function returns the node map of color indices. The values are

  2539     /// in the range \c [0..4] or \c [0..5] according to the coloring method.

  2540     IndexMap colorIndexMap() const {

  2541       return _color_map;

  2542     }

  2544     /// \brief Return the node map of colors

  2545     ///

  2546     /// This function returns the node map of colors. The values are among

  2547     /// five or six distinct \ref lemon::Color "colors".

  2548     ColorMap colorMap() const {

  2549       return composeMap(_palette, _color_map);

  2550     }

  2552     /// \brief Return the color index of the node

  2553     ///

  2554     /// This function returns the color index of the given node. The value is

  2555     /// in the range \c [0..4] or \c [0..5] according to the coloring method.

  2556     int colorIndex(const Node& node) const {

  2557       return _color_map[node];

  2558     }

  2560     /// \brief Return the color of the node

  2561     ///

  2562     /// This function returns the color of the given node. The value is among

  2563     /// five or six distinct \ref lemon::Color "colors".

  2564     Color color(const Node& node) const {

  2565       return _palette[_color_map[node]];

  2566     }

  2569     /// \brief Calculate a coloring with at most six colors

  2570     ///

  2571     /// This function calculates a coloring with at most six colors. The time

  2572     /// complexity of this variant is linear in the size of the graph.

  2573     /// \return \c true if the algorithm could color the graph with six colors.

  2574     /// If the algorithm fails, then the graph is not planar.

  2575     /// \note This function can return \c true if the graph is not

  2576     /// planar, but it can be colored with at most six colors.

  2577     bool runSixColoring() {

  2579       typename Graph::template NodeMap<int> heap_index(_graph, -1);

  2580       BucketHeap<typename Graph::template NodeMap<int> > heap(heap_index);

  2582       for (NodeIt n(_graph); n != INVALID; ++n) {

  2583         _color_map[n] = -2;

  2584         heap.push(n, countOutArcs(_graph, n));

  2585       }

  2587       std::vector<Node> order;

  2589       while (!heap.empty()) {

  2590         Node n = heap.top();

  2591         heap.pop();

  2592         _color_map[n] = -1;

  2593         order.push_back(n);

  2594         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

  2595           Node t = _graph.runningNode(e);

  2596           if (_color_map[t] == -2) {

  2597             heap.decrease(t, heap[t] - 1);

  2598           }

  2599         }

  2600       }

  2602       for (int i = order.size() - 1; i >= 0; --i) {

  2603         std::vector<bool> forbidden(6, false);

  2604         for (OutArcIt e(_graph, order[i]); e != INVALID; ++e) {

  2605           Node t = _graph.runningNode(e);

  2606           if (_color_map[t] != -1) {

  2607             forbidden[_color_map[t]] = true;

  2608           }

  2609         }

  2610                for (int k = 0; k < 6; ++k) {

  2611           if (!forbidden[k]) {

  2612             _color_map[order[i]] = k;

  2613             break;

  2614           }

  2615         }

  2616         if (_color_map[order[i]] == -1) {

  2617           return false;

  2618         }

  2619       }

  2620       return true;

  2621     }

  2623   private:

  2625     bool recolor(const Node& u, const Node& v) {

  2626       int ucolor = _color_map[u];

  2627       int vcolor = _color_map[v];

  2628       typedef _planarity_bits::KempeFilter<IndexMap> KempeFilter;

  2629       KempeFilter filter(_color_map, ucolor, vcolor);

  2631       typedef FilterNodes<const Graph, const KempeFilter> KempeGraph;

  2632       KempeGraph kempe_graph(_graph, filter);

  2634       std::vector<Node> comp;

  2635       Bfs<KempeGraph> bfs(kempe_graph);

  2636       bfs.init();

  2637       bfs.addSource(u);

  2638       while (!bfs.emptyQueue()) {

  2639         Node n = bfs.nextNode();

  2640         if (n == v) return false;

  2641         comp.push_back(n);

  2642         bfs.processNextNode();

  2643       }

  2645       int scolor = ucolor + vcolor;

  2646       for (int i = 0; i < static_cast<int>(comp.size()); ++i) {

  2647         _color_map[comp[i]] = scolor - _color_map[comp[i]];

  2648       }

  2650       return true;

  2651     }

  2653     template <typename EmbeddingMap>

  2654     void kempeRecoloring(const Node& node, const EmbeddingMap& embedding) {

  2655       std::vector<Node> nodes;

  2656       nodes.reserve(4);

  2658       for (Arc e = OutArcIt(_graph, node); e != INVALID; e = embedding[e]) {

  2659         Node t = _graph.target(e);

  2660         if (_color_map[t] != -1) {

  2661           nodes.push_back(t);

  2662           if (nodes.size() == 4) break;

  2663         }

  2664       }

  2666       int color = _color_map[nodes[0]];

  2667       if (recolor(nodes[0], nodes[2])) {

  2668         _color_map[node] = color;

  2669       } else {

  2670         color = _color_map[nodes[1]];

  2671         recolor(nodes[1], nodes[3]);

  2672         _color_map[node] = color;

  2673       }

  2674     }

  2676   public:

  2678     /// \brief Calculate a coloring with at most five colors

  2679     ///

  2680     /// This function calculates a coloring with at most five

  2681     /// colors. The worst case time complexity of this variant is

  2682     /// quadratic in the size of the graph.

  2683     /// \param embedding This map should contain a valid combinatorical

  2684     /// embedding, i.e. a valid cyclic order of the arcs.

  2685     /// It can be computed using PlanarEmbedding.

  2686     template <typename EmbeddingMap>

  2687     void runFiveColoring(const EmbeddingMap& embedding) {

  2689       typename Graph::template NodeMap<int> heap_index(_graph, -1);

  2690       BucketHeap<typename Graph::template NodeMap<int> > heap(heap_index);

  2692       for (NodeIt n(_graph); n != INVALID; ++n) {

  2693         _color_map[n] = -2;

  2694         heap.push(n, countOutArcs(_graph, n));

  2695       }

  2697       std::vector<Node> order;

  2699       while (!heap.empty()) {

  2700         Node n = heap.top();

  2701         heap.pop();

  2702         _color_map[n] = -1;

  2703         order.push_back(n);

  2704         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

  2705           Node t = _graph.runningNode(e);

  2706           if (_color_map[t] == -2) {

  2707             heap.decrease(t, heap[t] - 1);

  2708           }

  2709         }

  2710       }

  2712       for (int i = order.size() - 1; i >= 0; --i) {

  2713         std::vector<bool> forbidden(5, false);

  2714         for (OutArcIt e(_graph, order[i]); e != INVALID; ++e) {

  2715           Node t = _graph.runningNode(e);

  2716           if (_color_map[t] != -1) {

  2717             forbidden[_color_map[t]] = true;

  2718           }

  2719         }

  2720         for (int k = 0; k < 5; ++k) {

  2721           if (!forbidden[k]) {

  2722             _color_map[order[i]] = k;

  2723             break;

  2724           }

  2725         }

  2726         if (_color_map[order[i]] == -1) {

  2727           kempeRecoloring(order[i], embedding);

  2728         }

  2729       }

  2730     }

  2732     /// \brief Calculate a coloring with at most five colors

  2733     ///

  2734     /// This function calculates a coloring with at most five

  2735     /// colors. The worst case time complexity of this variant is

  2736     /// quadratic in the size of the graph.

  2737     /// \return \c true if the graph is planar.

  2738     bool runFiveColoring() {

  2739       PlanarEmbedding<Graph> pe(_graph);

  2740       if (!pe.run()) return false;

  2742       runFiveColoring(pe.embeddingMap());

  2743       return true;

  2744     }

  2746   private:

  2748     const Graph& _graph;

  2749     IndexMap _color_map;

  2750     Palette _palette;

  2751   };

  2753 }

  2755 #endif

author	Peter Kovacs <kpeter@inf.elte.hu>
	Tue, 15 Mar 2011 19:32:21 +0100
changeset 936	ddd3c0d3d9bf
parent 828	5fd7fafc4470
child 999	00f8d9f9920d
permissions	-rw-r--r--