lemon-main: lemon/planarity.h@2e959a5a0c2d

     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 target = _graph.target(arc);

    77         _tree_map[arc] = true;

    78         _pred_map[target] = arc;

    79       }

    81       void examine(const Arc& arc) {

    82         Node source = _graph.source(arc);

    83         Node target = _graph.target(arc);

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

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

    87             _low_map[source] = _order_map[target];

    88           }

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

    90             _ancestor_map[source] = _order_map[target];

    91           }

    92         }

    93       }

    95       void backtrack(const Arc& arc) {

    96         Node source = _graph.source(arc);

    97         Node target = _graph.target(arc);

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

   100           _low_map[source] = _low_map[target];

   101         }

   102       }

   104       const Graph& _graph;

   105       PredMap& _pred_map;

   106       TreeMap& _tree_map;

   107       OrderMap& _order_map;

   108       OrderList& _order_list;

   109       AncestorMap& _ancestor_map;

   110       LowMap& _low_map;

   111     };

   113     template <typename Graph, bool embedding = true>

   114     struct NodeDataNode {

   115       int prev, next;

   116       int visited;

   117       typename Graph::Arc first;

   118       bool inverted;

   119     };

   121     template <typename Graph>

   122     struct NodeDataNode<Graph, false> {

   123       int prev, next;

   124       int visited;

   125     };

   127     template <typename Graph>

   128     struct ChildListNode {

   129       typedef typename Graph::Node Node;

   130       Node first;

   131       Node prev, next;

   132     };

   134     template <typename Graph>

   135     struct ArcListNode {

   136       typename Graph::Arc prev, next;

   137     };

   139     template <typename Graph>

   140     class PlanarityChecking {

   141     private:

   143       TEMPLATE_GRAPH_TYPEDEFS(Graph);

   145       const Graph& _graph;

   147     private:

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

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

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

   154       typedef std::vector<Node> OrderList;

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

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

   159       typedef _planarity_bits::NodeDataNode<Graph> NodeDataNode;

   160       typedef std::vector<NodeDataNode> NodeData;

   162       typedef _planarity_bits::ChildListNode<Graph> ChildListNode;

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

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

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

   169     public:

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

   173       bool run() {

   174         typedef _planarity_bits::PlanarityVisitor<Graph> Visitor;

   176         PredMap pred_map(_graph, INVALID);

   177         TreeMap tree_map(_graph, false);

   179         OrderMap order_map(_graph, -1);

   180         OrderList order_list;

   182         AncestorMap ancestor_map(_graph, -1);

   183         LowMap low_map(_graph, -1);

   185         Visitor visitor(_graph, pred_map, tree_map,

   186                         order_map, order_list, ancestor_map, low_map);

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

   188         visit.run();

   190         ChildLists child_lists(_graph);

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

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

   195         EmbedArc embed_arc(_graph, false);

   197         MergeRoots merge_roots(_graph);

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

   201           Node node = order_list[i];

   203           Node source = node;

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

   205             Node target = _graph.target(e);

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

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

   209             }

   210           }

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

   213             Node target = _graph.target(e);

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

   216               embed_arc[target] = true;

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

   218                      order_map, order_list, node_data, merge_roots);

   219             }

   220           }

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

   223                  merge_roots[node].begin();

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

   225             int rn = *it;

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

   227                      ancestor_map, low_map, embed_arc, merge_roots);

   228           }

   229           merge_roots[node].clear();

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

   232             Node target = _graph.target(e);

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

   235               if (embed_arc[target]) {

   236                 return false;

   237               }

   238             }

   239           }

   240         }

   242         return true;

   243       }

   245     private:

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

   248                             const LowMap& low_map, ChildLists& child_lists) {

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

   251           Node source = n;

   253           std::vector<Node> targets;

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

   255             Node target = _graph.target(e);

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

   258               targets.push_back(target);

   259             }

   260           }

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

   263             child_lists[source].first = INVALID;

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

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

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

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

   268           } else {

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

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

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

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

   273             }

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

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

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

   277           }

   278         }

   279       }

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

   282                   const PredMap& pred_map, const LowMap& low_map,

   283                   const OrderMap& order_map, const OrderList& order_list,

   284                   NodeData& node_data, MergeRoots& merge_roots) {

   286         int na, nb;

   287         bool da, db;

   289         na = nb = order_map[node];

   290         da = true; db = false;

   292         while (true) {

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

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

   297           node_data[na].visited = rorder;

   298           node_data[nb].visited = rorder;

   300           int rn = -1;

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

   303             rn = na;

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

   305             rn = nb;

   306           }

   308           if (rn == -1) {

   309             int nn;

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

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

   313             na = nn;

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

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

   317             nb = nn;

   319           } else {

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

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

   324             if (low_map[rep] < rorder) {

   325               merge_roots[parent].push_back(rn);

   326             } else {

   327               merge_roots[parent].push_front(rn);

   328             }

   330             if (parent != root) {

   331               na = nb = order_map[parent];

   332               da = true; db = false;

   333             } else {

   334               break;

   335             }

   336           }

   337         }

   338       }

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

   341                     OrderList& order_list, ChildLists& child_lists,

   342                     AncestorMap& ancestor_map, LowMap& low_map,

   343                     EmbedArc& embed_arc, MergeRoots& merge_roots) {

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

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

   348           bool rd = di == 0;

   349           int pn = rn;

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

   352           while (n != rn) {

   354             Node node = order_list[n];

   356             if (embed_arc[node]) {

   358               // Merging components on the critical path

   359               while (!merge_stack.empty()) {

   361                 // Component root

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

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

   364                 merge_stack.pop_back();

   366                 // Parent of component

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

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

   369                 merge_stack.pop_back();

   371                 Node parent = order_list[dn];

   373                 // Erasing from merge_roots

   374                 merge_roots[parent].pop_front();

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

   378                 // Erasing from child_lists

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

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

   381                     child_lists[child].next;

   382                 } else {

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

   384                 }

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

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

   388                     child_lists[child].prev;

   389                 }

   391                 // Merging external faces

   392                 {

   393                   int en = cn;

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

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

   397                 }

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

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

   402               }

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

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

   407                   node_data[n].inverted) {

   408                 d = !d;

   409               }

   411               // Embedding arc into external face

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

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

   414               pn = rn;

   416               embed_arc[order_list[n]] = false;

   417             }

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

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

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

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

   427               int xn = node_data[rn].next;

   428               Node xnode = order_list[xn];

   430               int yn = node_data[rn].prev;

   431               Node ynode = order_list[yn];

   433               bool rd;

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

   435                             ancestor_map, low_map)) {

   436                 rd = true;

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

   438                                    ancestor_map, low_map)) {

   439                 rd = false;

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

   441                 rd = true;

   442               } else {

   443                 rd = false;

   444               }

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

   448               pn = rn;

   449               n = rd ? xn : yn;

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

   452                                  ancestor_map, low_map)) {

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

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

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

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

   459               else node_data[nn].next = pn;

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

   462               else node_data[pn].next = nn;

   464               node_data[nn].inverted =

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

   467               n = nn;

   468             }

   469             else break;

   471           }

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

   474             break;

   475           }

   476         }

   477       }

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

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

   481         int n = order_map[node];

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

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

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

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

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

   490       }

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

   493                     ChildLists& child_lists, AncestorMap& ancestor_map,

   494                     LowMap& low_map) {

   495         Node child = child_lists[node].first;

   497         if (child != INVALID) {

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

   499         }

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

   503         return false;

   504       }

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

   507                      const MergeRoots& merge_roots) {

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

   509       }

   511     };

   513   }

   515   /// \ingroup planar

   516   ///

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

   518   ///

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

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

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

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

   523   template <typename GR>

   524   bool checkPlanarity(const GR& graph) {

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

   526     return pc.run();

   527   }

   529   /// \ingroup planar

   530   ///

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

   532   ///

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

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

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

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

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

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

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

   540   ///

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

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

   543   ///

   544   /// \see PlanarDrawing, checkPlanarity()

   545   template <typename Graph>

   546   class PlanarEmbedding {

   547   private:

   549     TEMPLATE_GRAPH_TYPEDEFS(Graph);

   551     const Graph& _graph;

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

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

   556   private:

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

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

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

   563     typedef std::vector<Node> OrderList;

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

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

   568     typedef _planarity_bits::NodeDataNode<Graph> NodeDataNode;

   569     typedef std::vector<NodeDataNode> NodeData;

   571     typedef _planarity_bits::ChildListNode<Graph> ChildListNode;

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

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

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

   578     typedef _planarity_bits::ArcListNode<Graph> ArcListNode;

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

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

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

   585     enum IsolatorNodeType {

   586       HIGHX = 6, LOWX = 7,

   587       HIGHY = 8, LOWY = 9,

   588       ROOT = 10, PERTINENT = 11,

   589       INTERNAL = 12

   590     };

   592   public:

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

   595     ///

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

   597     /// \see embeddingMap()

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

   600     /// \brief Constructor

   601     ///

   602     /// Constructor.

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

   604     /// contain parallel or loop arcs.

   605     PlanarEmbedding(const Graph& graph)

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

   608     /// \brief Run the algorithm.

   609     ///

   610     /// This function runs the algorithm.

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

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

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

   614     bool run(bool kuratowski = true) {

   615       typedef _planarity_bits::PlanarityVisitor<Graph> Visitor;

   617       PredMap pred_map(_graph, INVALID);

   618       TreeMap tree_map(_graph, false);

   620       OrderMap order_map(_graph, -1);

   621       OrderList order_list;

   623       AncestorMap ancestor_map(_graph, -1);

   624       LowMap low_map(_graph, -1);

   626       Visitor visitor(_graph, pred_map, tree_map,

   627                       order_map, order_list, ancestor_map, low_map);

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

   629       visit.run();

   631       ChildLists child_lists(_graph);

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

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

   636       EmbedArc embed_arc(_graph, INVALID);

   638       MergeRoots merge_roots(_graph);

   640       ArcLists arc_lists(_graph);

   642       FlipMap flip_map(_graph, false);

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

   646         Node node = order_list[i];

   648         node_data[i].first = INVALID;

   650         Node source = node;

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

   652           Node target = _graph.target(e);

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

   655             initFace(target, arc_lists, node_data,

   656                      pred_map, order_map, order_list);

   657           }

   658         }

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

   661           Node target = _graph.target(e);

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

   664             embed_arc[target] = e;

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

   666                    order_map, order_list, node_data, merge_roots);

   667           }

   668         }

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

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

   672           int rn = *it;

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

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

   675         }

   676         merge_roots[node].clear();

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

   679           Node target = _graph.target(e);

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

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

   683               if (kuratowski) {

   684                 isolateKuratowski(e, node_data, arc_lists, flip_map,

   685                                   order_map, order_list, pred_map, child_lists,

   686                                   ancestor_map, low_map,

   687                                   embed_arc, merge_roots);

   688               }

   689               return false;

   690             }

   691           }

   692         }

   693       }

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

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

   698                             child_lists, arc_lists);

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

   700                        arc_lists, flip_map);

   701       }

   703       return true;

   704     }

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

   707     ///

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

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

   710     /// a node.

   711     Arc next(const Arc& arc) const {

   712       return _embedding[arc];

   713     }

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

   716     ///

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

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

   719     /// outgoing arcs of its source node.

   720     const EmbeddingMap& embeddingMap() const {

   721       return _embedding;

   722     }

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

   725     /// subdivision

   726     ///

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

   728     /// Kuratowski subdivision.

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

   730     /// before using this function.

   731     bool kuratowski(const Edge& edge) const {

   732       return _kuratowski[edge];

   733     }

   735   private:

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

   738                           const LowMap& low_map, ChildLists& child_lists) {

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

   741         Node source = n;

   743         std::vector<Node> targets;

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

   745           Node target = _graph.target(e);

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

   748             targets.push_back(target);

   749           }

   750         }

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

   753           child_lists[source].first = INVALID;

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

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

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

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

   758         } else {

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

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

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

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

   763           }

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

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

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

   767         }

   768       }

   769     }

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

   772                 const PredMap& pred_map, const LowMap& low_map,

   773                 const OrderMap& order_map, const OrderList& order_list,

   774                 NodeData& node_data, MergeRoots& merge_roots) {

   776       int na, nb;

   777       bool da, db;

   779       na = nb = order_map[node];

   780       da = true; db = false;

   782       while (true) {

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

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

   787         node_data[na].visited = rorder;

   788         node_data[nb].visited = rorder;

   790         int rn = -1;

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

   793           rn = na;

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

   795           rn = nb;

   796         }

   798         if (rn == -1) {

   799           int nn;

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

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

   803           na = nn;

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

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

   807           nb = nn;

   809         } else {

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

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

   814           if (low_map[rep] < rorder) {

   815             merge_roots[parent].push_back(rn);

   816           } else {

   817             merge_roots[parent].push_front(rn);

   818           }

   820           if (parent != root) {

   821             na = nb = order_map[parent];

   822             da = true; db = false;

   823           } else {

   824             break;

   825           }

   826         }

   827       }

   828     }

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

   831                   ArcLists& arc_lists, FlipMap& flip_map,

   832                   OrderList& order_list, ChildLists& child_lists,

   833                   AncestorMap& ancestor_map, LowMap& low_map,

   834                   EmbedArc& embed_arc, MergeRoots& merge_roots) {

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

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

   839         bool rd = di == 0;

   840         int pn = rn;

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

   843         while (n != rn) {

   845           Node node = order_list[n];

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

   849             // Merging components on the critical path

   850             while (!merge_stack.empty()) {

   852               // Component root

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

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

   855               merge_stack.pop_back();

   857               // Parent of component

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

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

   860               merge_stack.pop_back();

   862               Node parent = order_list[dn];

   864               // Erasing from merge_roots

   865               merge_roots[parent].pop_front();

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

   869               // Erasing from child_lists

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

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

   872                   child_lists[child].next;

   873               } else {

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

   875               }

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

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

   879                   child_lists[child].prev;

   880               }

   882               // Merging arcs + flipping

   883               Arc de = node_data[dn].first;

   884               Arc ce = node_data[cn].first;

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

   887               if (cd != dd) {

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

   889                 ce = arc_lists[ce].prev;

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

   891               }

   893               {

   894                 Arc dne = arc_lists[de].next;

   895                 Arc cne = arc_lists[ce].next;

   897                 arc_lists[de].next = cne;

   898                 arc_lists[ce].next = dne;

   900                 arc_lists[dne].prev = ce;

   901                 arc_lists[cne].prev = de;

   902               }

   904               if (dd) {

   905                 node_data[dn].first = ce;

   906               }

   908               // Merging external faces

   909               {

   910                 int en = cn;

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

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

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

   915                     node_data[cn].inverted) {

   916                    cd = !cd;

   917                  }

   918               }

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

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

   923             }

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

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

   928                 node_data[n].inverted) {

   929               d = !d;

   930             }

   932             // Add new arc

   933             {

   934               Arc arc = embed_arc[node];

   935               Arc re = node_data[rn].first;

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

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

   939               arc_lists[arc].prev = re;

   940               arc_lists[re].next = arc;

   942               if (!rd) {

   943                 node_data[rn].first = arc;

   944               }

   946               Arc rev = _graph.oppositeArc(arc);

   947               Arc e = node_data[n].first;

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

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

   951               arc_lists[rev].prev = e;

   952               arc_lists[e].next = rev;

   954               if (d) {

   955                 node_data[n].first = rev;

   956               }

   958             }

   960             // Embedding arc into external face

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

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

   963             pn = rn;

   965             embed_arc[order_list[n]] = INVALID;

   966           }

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

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

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

   972                 node_data[n].inverted) {

   973               d = !d;

   974             }

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

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

   980             int xn = node_data[rn].next;

   981             Node xnode = order_list[xn];

   983             int yn = node_data[rn].prev;

   984             Node ynode = order_list[yn];

   986             bool rd;

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

   988               rd = true;

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

   990                                  ancestor_map, low_map)) {

   991               rd = false;

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

   993               rd = true;

   994             } else {

   995               rd = false;

   996             }

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

  1000             pn = rn;

  1001             n = rd ? xn : yn;

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

  1004                                ancestor_map, low_map)) {

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

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

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

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

  1011             else node_data[nn].next = pn;

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

  1014             else node_data[pn].next = nn;

  1016             node_data[nn].inverted =

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

  1019             n = nn;

  1020           }

  1021           else break;

  1023         }

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

  1026           break;

  1027         }

  1028       }

  1029     }

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

  1032                   NodeData& node_data, const PredMap& pred_map,

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

  1034       int n = order_map[node];

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

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

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

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

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

  1043       node_data[n].inverted = false;

  1044       node_data[rn].inverted = false;

  1046       Arc arc = pred_map[node];

  1047       Arc rev = _graph.oppositeArc(arc);

  1049       node_data[rn].first = arc;

  1050       node_data[n].first = rev;

  1052       arc_lists[arc].prev = arc;

  1053       arc_lists[arc].next = arc;

  1055       arc_lists[rev].prev = rev;

  1056       arc_lists[rev].next = rev;

  1058     }

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

  1061                              OrderList& order_list, OrderMap& order_map,

  1062                              ChildLists& child_lists, ArcLists& arc_lists) {

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

  1064         int dd = order_map[node];

  1065         Node child = child_lists[node].first;

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

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

  1069         Arc de = node_data[dd].first;

  1070         Arc ce = node_data[cd].first;

  1072         if (de != INVALID) {

  1073           Arc dne = arc_lists[de].next;

  1074           Arc cne = arc_lists[ce].next;

  1076           arc_lists[de].next = cne;

  1077           arc_lists[ce].next = dne;

  1079           arc_lists[dne].prev = ce;

  1080           arc_lists[cne].prev = de;

  1081         }

  1083         node_data[dd].first = ce;

  1085       }

  1086     }

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

  1089                         OrderMap& order_map, PredMap& pred_map,

  1090                         ArcLists& arc_lists, FlipMap& flip_map) {

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

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

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

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

  1097       }

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

  1100       Arc prev = first;

  1102       Arc arc = flip_map[node] ?

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

  1105       _embedding[prev] = arc;

  1107       while (arc != first) {

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

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

  1110         prev = arc; arc = next;

  1111         _embedding[prev] = arc;

  1112       }

  1113     }

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

  1117                   ChildLists& child_lists, AncestorMap& ancestor_map,

  1118                   LowMap& low_map) {

  1119       Node child = child_lists[node].first;

  1121       if (child != INVALID) {

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

  1123       }

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

  1127       return false;

  1128     }

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

  1131                    const MergeRoots& merge_roots) {

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

  1133     }

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

  1136                  AncestorMap& ancestor_map, LowMap& low_map) {

  1137       int low_point;

  1139       Node child = child_lists[node].first;

  1141       if (child != INVALID) {

  1142         low_point = low_map[child];

  1143       } else {

  1144         low_point = order_map[node];

  1145       }

  1147       if (low_point > ancestor_map[node]) {

  1148         low_point = ancestor_map[node];

  1149       }

  1151       return low_point;

  1152     }

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

  1155                           OrderMap& order_map, OrderList& order_list) {

  1157       int order = order_map[root];

  1158       int norder = order_map[node];

  1160       Node child = child_lists[root].first;

  1161       while (child != INVALID) {

  1162         int corder = order_map[child];

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

  1164           order = corder;

  1165         }

  1166         child = child_lists[child].next;

  1167       }

  1168       return order + order_list.size();

  1169     }

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

  1172                        EmbedArc& embed_arc, MergeRoots& merge_roots) {

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

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

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

  1176       }

  1177       return wnode;

  1178     }

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

  1182                       ChildLists& child_lists, AncestorMap& ancestor_map,

  1183                       LowMap& low_map, NodeData& node_data) {

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

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

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

  1187       }

  1188       return wnode;

  1189     }

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

  1192                         OrderList& order_list, OrderMap& order_map,

  1193                         NodeData& node_data, ArcLists& arc_lists,

  1194                         EmbedArc& embed_arc, MergeRoots& merge_roots,

  1195                         ChildLists& child_lists, AncestorMap& ancestor_map,

  1196                         LowMap& low_map) {

  1198       Node cnode = node;

  1199       Node pred = INVALID;

  1201       while (true) {

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

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

  1206         if (pert && ext) {

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

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

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

  1211               Arc arc = node_data[cn].first;

  1212               _kuratowski.set(arc, true);

  1214               pred = cnode;

  1215               cnode = _graph.target(arc);

  1217               continue;

  1218             }

  1219           }

  1220           wnode = znode = cnode;

  1221           return;

  1223         } else if (pert) {

  1224           wnode = cnode;

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

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

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

  1230               arc = arc_lists[arc].next;

  1231             }

  1232             _kuratowski.set(arc, true);

  1234             Node next = _graph.target(arc);

  1235             pred = cnode; cnode = next;

  1236           }

  1238           znode = cnode;

  1239           return;

  1241         } else if (ext) {

  1242           znode = cnode;

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

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

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

  1248               arc = arc_lists[arc].next;

  1249             }

  1250             _kuratowski.set(arc, true);

  1252             Node next = _graph.target(arc);

  1253             pred = cnode; cnode = next;

  1254           }

  1256           wnode = cnode;

  1257           return;

  1259         } else {

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

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

  1263             arc = arc_lists[arc].next;

  1264           }

  1265           _kuratowski.set(arc, true);

  1267           Node next = _graph.target(arc);

  1268           pred = cnode; cnode = next;

  1269         }

  1271       }

  1273     }

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

  1276                          PredMap& pred_map, NodeData& node_data,

  1277                          ArcLists& arc_lists, FlipMap& flip_map,

  1278                          TypeMap& type_map) {

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

  1280       type_map[root] = 1;

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

  1284       st.push_back(root);

  1285       while (!st.empty()) {

  1286         Node node = st.back();

  1287         st.pop_back();

  1288         qu.push_back(node);

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

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

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

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

  1295         }

  1297         Arc last = arc, pred = arc;

  1298         arc = arc_lists[arc].next;

  1299         while (arc != last) {

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

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

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

  1304           }

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

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

  1308           pred = arc; arc = next;

  1309         }

  1311       }

  1313       type_map[root] = 2;

  1314       flip_map[root] = false;

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

  1318         Node node = qu[i];

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

  1321           st.push_back(node);

  1322           type_map[node] = 2;

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

  1324         }

  1326         bool flip = flip_map[node];

  1328         while (!st.empty()) {

  1329           node = st.back();

  1330           st.pop_back();

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

  1333           flip = flip_map[node];

  1335           if (flip) {

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

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

  1338             arc = arc_lists[arc].prev;

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

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

  1341           }

  1342         }

  1343       }

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

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

  1348         Arc last = arc, pred = arc;

  1350         arc = arc_lists[arc].next;

  1351         while (arc != last) {

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

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

  1355           }

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

  1357         }

  1359       }

  1360     }

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

  1363                       OrderMap& order_map, NodeData& node_data,

  1364                       TypeMap& type_map) {

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

  1367       while (node != ynode) {

  1368         type_map[node] = HIGHY;

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

  1370       }

  1372       while (node != wnode) {

  1373         type_map[node] = LOWY;

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

  1375       }

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

  1379       while (node != xnode) {

  1380         type_map[node] = LOWX;

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

  1382       }

  1383       type_map[node] = LOWX;

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

  1386       while (node != root) {

  1387         type_map[node] = HIGHX;

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

  1389       }

  1391       type_map[wnode] = PERTINENT;

  1392       type_map[root] = ROOT;

  1393     }

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

  1396                           Node wnode, Node root, TypeMap& type_map,

  1397                           OrderMap& order_map, NodeData& node_data,

  1398                           ArcLists& arc_lists) {

  1399       std::vector<Arc> st;

  1401       Node node = wnode;

  1403       while (node != root) {

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

  1405         st.push_back(arc);

  1406         node = _graph.target(arc);

  1407       }

  1409       while (true) {

  1410         Arc arc = st.back();

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

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

  1413           break;

  1414         }

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

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

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

  1419           st.push_back(arc);

  1420         } else {

  1421           st.pop_back();

  1422           arc = arc_lists[arc].next;

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

  1425             arc = st.back();

  1426             st.pop_back();

  1427             arc = arc_lists[arc].next;

  1428           }

  1429           st.push_back(arc);

  1430         }

  1431       }

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

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

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

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

  1437             ipath.push_back(st[i]);

  1438           }

  1439         }

  1440       }

  1441     }

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

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

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

  1446       }

  1447     }

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

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

  1451                       NodeData& node_data, ArcLists& arc_lists) {

  1452       std::vector<Arc> st;

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

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

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

  1458         Arc arc = st.back();

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

  1460           break;

  1461         }

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

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

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

  1466           st.push_back(arc);

  1467         } else {

  1468           st.pop_back();

  1469           arc = arc_lists[arc].next;

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

  1472             arc = st.back();

  1473             st.pop_back();

  1474             arc = arc_lists[arc].next;

  1475           }

  1476           st.push_back(arc);

  1477         }

  1478       }

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

  1481         ppath.push_back(st[i]);

  1482       }

  1483     }

  1486     int markExternalPath(Node node, OrderMap& order_map,

  1487                          ChildLists& child_lists, PredMap& pred_map,

  1488                          AncestorMap& ancestor_map, LowMap& low_map) {

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

  1490                         ancestor_map, low_map);

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

  1493         node = child_lists[node].first;

  1494         _kuratowski[pred_map[node]] = true;

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

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

  1498             Node tnode = _graph.target(e);

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

  1500               node = tnode;

  1501               _kuratowski[e] = true;

  1502               break;

  1503             }

  1504           }

  1505         }

  1506       }

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

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

  1510           _kuratowski[e] = true;

  1511           break;

  1512         }

  1513       }

  1515       return lp;

  1516     }

  1518     void markPertinentPath(Node node, OrderMap& order_map,

  1519                            NodeData& node_data, ArcLists& arc_lists,

  1520                            EmbedArc& embed_arc, MergeRoots& merge_roots) {

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

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

  1523         Arc arc = node_data[n].first;

  1525         _kuratowski.set(arc, true);

  1527         Node pred = node;

  1528         node = _graph.target(arc);

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

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

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

  1532             arc = arc_lists[arc].next;

  1533           }

  1534           _kuratowski.set(arc, true);

  1535           pred = node;

  1536           node = _graph.target(arc);

  1537         }

  1538       }

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

  1540     }

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

  1543       while (node != snode) {

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

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

  1546       }

  1547     }

  1549     void markFacePath(Node ynode, Node xnode,

  1550                       OrderMap& order_map, NodeData& node_data) {

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

  1552       Node node = _graph.target(arc);

  1553       _kuratowski.set(arc, true);

  1555       while (node != xnode) {

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

  1557         _kuratowski.set(arc, true);

  1558         node = _graph.target(arc);

  1559       }

  1560     }

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

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

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

  1565       }

  1566     }

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

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

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

  1571       }

  1572     }

  1574     void isolateKuratowski(Arc arc, NodeData& node_data,

  1575                            ArcLists& arc_lists, FlipMap& flip_map,

  1576                            OrderMap& order_map, OrderList& order_list,

  1577                            PredMap& pred_map, ChildLists& child_lists,

  1578                            AncestorMap& ancestor_map, LowMap& low_map,

  1579                            EmbedArc& embed_arc, MergeRoots& merge_roots) {

  1581       Node root = _graph.source(arc);

  1582       Node enode = _graph.target(arc);

  1584       int rorder = order_map[root];

  1586       TypeMap type_map(_graph, 0);

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

  1589                                  order_map, order_list);

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

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

  1594       // Minor-A

  1595       {

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

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

  1599             root = xnode;

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

  1601           } else {

  1602             root = ynode;

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

  1604           }

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

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

  1608         }

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

  1611           orientComponent(root, rn, order_map, pred_map,

  1612                           node_data, arc_lists, flip_map, type_map);

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

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

  1615                                      pred_map, ancestor_map, low_map);

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

  1617                                      pred_map, ancestor_map, low_map);

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

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

  1620                                       embed_arc, merge_roots);

  1622           markPertinentPath(lwnode, order_map, node_data, arc_lists,

  1623                             embed_arc, merge_roots);

  1625           return;

  1626         }

  1627       }

  1629       orientComponent(root, rn, order_map, pred_map,

  1630                       node_data, arc_lists, flip_map, type_map);

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

  1633                                  embed_arc, merge_roots);

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

  1637       //Minor-B

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

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

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

  1641         if (low_map[rep] < rorder) {

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

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

  1644                                      pred_map, ancestor_map, low_map);

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

  1646                                      pred_map, ancestor_map, low_map);

  1648           Node lwnode, lznode;

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

  1650                          order_map, node_data, arc_lists, embed_arc,

  1651                          merge_roots, child_lists, ancestor_map, low_map);

  1653           markPertinentPath(lwnode, order_map, node_data, arc_lists,

  1654                             embed_arc, merge_roots);

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

  1656                                      pred_map, ancestor_map, low_map);

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

  1659           if (zlp < minlp) minlp = zlp;

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

  1662           if (zlp > maxlp) maxlp = zlp;

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

  1666           return;

  1667         }

  1668       }

  1670       Node pxnode, pynode;

  1671       std::vector<Arc> ipath;

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

  1673                        node_data, arc_lists);

  1674       setInternalFlags(ipath, type_map);

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

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

  1678       wnode = findPertinent(pynode, order_map, node_data,

  1679                             embed_arc, merge_roots);

  1681       // Minor-C

  1682       {

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

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

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

  1686           }

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

  1688           markPertinentPath(wnode, order_map, node_data, arc_lists,

  1689                             embed_arc, merge_roots);

  1690           markInternalPath(ipath);

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

  1692                                      pred_map, ancestor_map, low_map);

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

  1694                                      pred_map, ancestor_map, low_map);

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

  1696           return;

  1697         }

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

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

  1701           markPertinentPath(wnode, order_map, node_data, arc_lists,

  1702                             embed_arc, merge_roots);

  1703           markInternalPath(ipath);

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

  1705                                      pred_map, ancestor_map, low_map);

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

  1707                                      pred_map, ancestor_map, low_map);

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

  1709           return;

  1710         }

  1711       }

  1713       std::vector<Arc> ppath;

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

  1716       // Minor-D

  1717       if (!ppath.empty()) {

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

  1719         markPertinentPath(wnode, order_map, node_data, arc_lists,

  1720                           embed_arc, merge_roots);

  1721         markPilePath(ppath);

  1722         markInternalPath(ipath);

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

  1724                                    pred_map, ancestor_map, low_map);

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

  1726                                    pred_map, ancestor_map, low_map);

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

  1728         return;

  1729       }

  1731       // Minor-E*

  1732       {

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

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

  1736                                     child_lists, ancestor_map,

  1737                                     low_map, node_data);

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

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

  1741             markPertinentPath(wnode, order_map, node_data, arc_lists,

  1742                               embed_arc, merge_roots);

  1743             markInternalPath(ipath);

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

  1745                                        pred_map, ancestor_map, low_map);

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

  1747                                        pred_map, ancestor_map, low_map);

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

  1749           } else {

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

  1751             markPertinentPath(wnode, order_map, node_data, arc_lists,

  1752                               embed_arc, merge_roots);

  1753             markInternalPath(ipath);

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

  1755                                        pred_map, ancestor_map, low_map);

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

  1757                                        pred_map, ancestor_map, low_map);

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

  1759           }

  1760           return;

  1761         }

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

  1764                                    pred_map, ancestor_map, low_map);

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

  1766                                    pred_map, ancestor_map, low_map);

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

  1768                                    pred_map, ancestor_map, low_map);

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

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

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

  1773           return;

  1774         }

  1776         markInternalPath(ipath);

  1777         markPertinentPath(wnode, order_map, node_data, arc_lists,

  1778                           embed_arc, merge_roots);

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

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

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

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

  1784           return;

  1785         }

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

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

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

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

  1791           return;

  1792         }

  1794         if (pynode != ynode) {

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

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

  1798           if (wlp < minlp) minlp = wlp;

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

  1801           if (wlp > maxlp) maxlp = wlp;

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

  1804           return;

  1805         }

  1807         if (pxnode != xnode) {

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

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

  1811           if (wlp < minlp) minlp = wlp;

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

  1814           if (wlp > maxlp) maxlp = wlp;

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

  1817           return;

  1818         }

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

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

  1822         if (wlp < minlp) minlp = wlp;

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

  1824         return;

  1825       }

  1827     }

  1829   };

  1831   namespace _planarity_bits {

  1833     template <typename Graph, typename EmbeddingMap>

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

  1835       DfsVisitor<Graph> null_visitor;

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

  1837       dfs.init();

  1839       typename Graph::Node u = INVALID;

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

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

  1842           dfs.addSource(n);

  1843           dfs.start();

  1844           if (u == INVALID) {

  1845             u = n;

  1846           } else {

  1847             typename Graph::Node v = n;

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

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

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

  1854             if (ue != INVALID) {

  1855               embedding[e] = embedding[ue];

  1856               embedding[ue] = e;

  1857             } else {

  1858               embedding[e] = e;

  1859             }

  1861             if (ve != INVALID) {

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

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

  1864             } else {

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

  1866             }

  1867           }

  1868         }

  1869       }

  1870     }

  1872     template <typename Graph, typename EmbeddingMap>

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

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

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

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

  1878         arcs.push_back(e);

  1879       }

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

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

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

  1885         if (processed[pp]) continue;

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

  1888         processed[e] = true;

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

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

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

  1894         while (e != l) {

  1895           processed[e] = true;

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

  1899             typename Graph::Arc n =

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

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

  1902             embedding[n] = p;

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

  1905             embedding[graph.oppositeArc(n)] =

  1906               embedding[graph.oppositeArc(e)];

  1907             embedding[graph.oppositeArc(e)] =

  1908               graph.oppositeArc(n);

  1910             p = n;

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

  1912           } else {

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

  1914             pp = p;

  1915             p = e;

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

  1917           }

  1918         }

  1919         visited.setAll(false);

  1920       }

  1921     }

  1924     template <typename Graph, typename EmbeddingMap>

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

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

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

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

  1931       }

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

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

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

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

  1938         arcs.push_back(e);

  1939       }

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

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

  1944         if (processed[e]) continue;

  1945         processed[e] = true;

  1947         typename Graph::Arc mine = e;

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

  1950         int face_size = 1;

  1952         typename Graph::Arc l = e;

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

  1954         while (l != e) {

  1955           processed[e] = true;

  1957           ++face_size;

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

  1960             mine = e;

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

  1962           }

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

  1965         }

  1967         if (face_size < 4) {

  1968           continue;

  1969         }

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

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

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

  1974         }

  1976         typename Graph::Arc oppe = INVALID;

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

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

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

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

  1982             oppe = e;

  1983             break;

  1984           }

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

  1986         }

  1987         visited.setAll(false);

  1989         if (oppe == INVALID) {

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

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

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

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

  1996             typename Graph::Arc n =

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

  1999             embedding[n] = pn;

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

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

  2003             pn = n;

  2005             p = e;

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

  2007           }

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

  2011         } else {

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

  2014           s = graph.source(mine);

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

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

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

  2019           embedding[ce] = mine;

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

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

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

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

  2025             typename Graph::Arc n =

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

  2028             embedding[n] = pn;

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

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

  2032             pn = n;

  2034             p = e;

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

  2037           }

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

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

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

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

  2043             typename Graph::Arc n =

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

  2046             embedding[n] = pn;

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

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

  2050             pn = n;

  2052             p = e;

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

  2055           }

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

  2057         }

  2058       }

  2059     }

  2061   }

  2063   /// \ingroup planar

  2064   ///

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

  2066   ///

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

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

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

  2070   /// each other.

  2071   ///

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

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

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

  2075   ///

  2076   /// \see PlanarEmbedding

  2077   template <typename Graph>

  2078   class PlanarDrawing {

  2079   public:

  2081     TEMPLATE_GRAPH_TYPEDEFS(Graph);

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

  2084     typedef dim2::Point<int> Point;

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

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

  2089     /// \brief Constructor

  2090     ///

  2091     /// Constructor

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

  2093     /// contain parallel or loop arcs.

  2094     PlanarDrawing(const Graph& graph)

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

  2097   private:

  2099     template <typename AuxGraph, typename AuxEmbeddingMap>

  2100     void drawing(const AuxGraph& graph,

  2101                  const AuxEmbeddingMap& next,

  2102                  PointMap& point_map) {

  2103       TEMPLATE_GRAPH_TYPEDEFS(AuxGraph);

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

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

  2108         Arc e = OutArcIt(graph, n);

  2110         Arc p = e, l = e;

  2112         e = next[e];

  2113         while (e != l) {

  2114           prev[e] = p;

  2115           p = e;

  2116           e = next[e];

  2117         }

  2118         prev[e] = p;

  2119       }

  2121       Node anode, bnode, cnode;

  2123       {

  2124         Arc e = ArcIt(graph);

  2125         anode = graph.source(e);

  2126         bnode = graph.target(e);

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

  2128       }

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

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

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

  2134       {

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

  2136           Node m = graph.target(e);

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

  2138             conn[m] = 1;

  2139           }

  2140         }

  2141         conn[cnode] = 2;

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

  2144           Node m = graph.target(e);

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

  2146             conn[m] = 1;

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

  2148             conn[m] += 1;

  2149             Arc pe = graph.oppositeArc(e);

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

  2151               conn[m] -= 1;

  2152             }

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

  2154               conn[m] -= 1;

  2155             }

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

  2158           }

  2159         }

  2160       }

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

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

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

  2167         proper.set(n, false);

  2168         conn[n] = -2;

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

  2171           Node m = graph.target(e);

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

  2173             conn[m] = 1;

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

  2175             conn[m] += 1;

  2176             Arc pe = graph.oppositeArc(e);

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

  2178               conn[m] -= 1;

  2179             }

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

  2181               conn[m] -= 1;

  2182             }

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

  2185           }

  2186         }

  2188         {

  2189           Arc e = OutArcIt(graph, n);

  2190           Arc p = e, l = e;

  2192           e = next[e];

  2193           while (e != l) {

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

  2196               Arc f = e;

  2197               angle[f] = 0;

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

  2199               angle[f] = 1;

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

  2201               angle[f] = 2;

  2202             }

  2204             p = e;

  2205             e = next[e];

  2206           }

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

  2209             Arc f = e;

  2210             angle[f] = 0;

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

  2212             angle[f] = 1;

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

  2214             angle[f] = 2;

  2215           }

  2216         }

  2217       }

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

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

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

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

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

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

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

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

  2229           switch (angle[e]) {

  2230           case 2:

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

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

  2233             break;

  2234           case 1:

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

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

  2237             break;

  2238           case 0:

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

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

  2241             break;

  2242           }

  2243         }

  2244       }

  2246       cpred[anode] = INVALID;

  2247       cpred[bnode] = INVALID;

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

  2251       {

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

  2253         std::vector<Node> st;

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

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

  2256             st.push_back(n);

  2257             processed[n] = true;

  2258             Node m = apred[n];

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

  2260               st.push_back(m);

  2261               processed[m] = true;

  2262               m = apred[m];

  2263             }

  2264             while (!st.empty()) {

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

  2266               st.pop_back();

  2267             }

  2268           }

  2269         }

  2270       }

  2272       {

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

  2274         std::vector<Node> st;

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

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

  2277             st.push_back(n);

  2278             processed[n] = true;

  2279             Node m = bpred[n];

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

  2281               st.push_back(m);

  2282               processed[m] = true;

  2283               m = bpred[m];

  2284             }

  2285             while (!st.empty()) {

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

  2287               st.pop_back();

  2288             }

  2289           }

  2290         }

  2291       }

  2293       {

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

  2295         std::vector<Node> st;

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

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

  2298             st.push_back(n);

  2299             processed[n] = true;

  2300             Node m = cpred[n];

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

  2302               st.push_back(m);

  2303               processed[m] = true;

  2304               m = cpred[m];

  2305             }

  2306             while (!st.empty()) {

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

  2308               st.pop_back();

  2309             }

  2310           }

  2311         }

  2312       }

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

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

  2316         Node n = aorder[i];

  2317         atree[n] = 1;

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

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

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

  2321           }

  2322         }

  2323       }

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

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

  2327         Node n = border[i];

  2328         btree[n] = 1;

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

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

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

  2332           }

  2333         }

  2334       }

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

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

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

  2339       apath_btree[bnode] = btree[bnode];

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

  2341         Node n = aorder[i];

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

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

  2344       }

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

  2347       bpath_atree[anode] = atree[anode];

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

  2349         Node n = border[i];

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

  2351       }

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

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

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

  2356       cpath_atree[anode] = atree[anode];

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

  2358       cpath_btree[bnode] = btree[bnode];

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

  2360         Node n = corder[i];

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

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

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

  2364       }

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

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

  2368         point_map[n].x =

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

  2370         point_map[n].y =

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

  2372       }

  2374     }

  2376   public:

  2378     /// \brief Calculate the node positions

  2379     ///

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

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

  2382     bool run() {

  2383       PlanarEmbedding<Graph> pe(_graph);

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

  2386       run(pe);

  2387       return true;

  2388     }

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

  2391     /// combinatorical embedding

  2392     ///

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

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

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

  2396     /// It can be computed using PlanarEmbedding.

  2397     template <typename EmbeddingMap>

  2398     void run(const EmbeddingMap& embedding) {

  2399       typedef SmartEdgeSet<Graph> AuxGraph;

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

  2402         drawing(_graph, embedding, _point_map);

  2403         return;

  2404       }

  2406       AuxGraph aux_graph(_graph);

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

  2408         aux_embedding(aux_graph);

  2410       {

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

  2413           ref(_graph);

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

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

  2417         }

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

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

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

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

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

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

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

  2426         }

  2427       }

  2428       _planarity_bits::makeConnected(aux_graph, aux_embedding);

  2429       _planarity_bits::makeBiNodeConnected(aux_graph, aux_embedding);

  2430       _planarity_bits::makeMaxPlanar(aux_graph, aux_embedding);

  2431       drawing(aux_graph, aux_embedding, _point_map);

  2432     }

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

  2435     ///

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

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

  2438       return _point_map[node];

  2439     }

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

  2442     ///

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

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

  2445     const PointMap& coords() const {

  2446       return _point_map;

  2447     }

  2449   private:

  2451     const Graph& _graph;

  2452     PointMap _point_map;

  2454   };

  2456   namespace _planarity_bits {

  2458     template <typename ColorMap>

  2459     class KempeFilter {

  2460     public:

  2461       typedef typename ColorMap::Key Key;

  2462       typedef bool Value;

  2464       KempeFilter(const ColorMap& color_map,

  2465                   const typename ColorMap::Value& first,

  2466                   const typename ColorMap::Value& second)

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

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

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

  2471       }

  2473     private:

  2474       const ColorMap& _color_map;

  2475       typename ColorMap::Value _first, _second;

  2476     };

  2477   }

  2479   /// \ingroup planar

  2480   ///

  2481   /// \brief Coloring planar graphs

  2482   ///

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

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

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

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

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

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

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

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

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

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

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

  2494   ///

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

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

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

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

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

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

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

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

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

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

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

  2506   template <typename Graph>

  2507   class PlanarColoring {

  2508   public:

  2510     TEMPLATE_GRAPH_TYPEDEFS(Graph);

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

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

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

  2515     ///

  2516     /// The map type for storing colors.

  2517     /// \see Palette, Color

  2518     typedef ComposeMap<Palette, IndexMap> ColorMap;

  2520     /// \brief Constructor

  2521     ///

  2522     /// Constructor.

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

  2524     /// contain parallel or loop arcs.

  2525     PlanarColoring(const Graph& graph)

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

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

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

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

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

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

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

  2533     }

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

  2536     ///

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

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

  2539     IndexMap colorIndexMap() const {

  2540       return _color_map;

  2541     }

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

  2544     ///

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

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

  2547     ColorMap colorMap() const {

  2548       return composeMap(_palette, _color_map);

  2549     }

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

  2552     ///

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

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

  2555     int colorIndex(const Node& node) const {

  2556       return _color_map[node];

  2557     }

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

  2560     ///

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

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

  2563     Color color(const Node& node) const {

  2564       return _palette[_color_map[node]];

  2565     }

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

  2569     ///

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

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

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

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

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

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

  2576     bool runSixColoring() {

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

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

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

  2582         _color_map[n] = -2;

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

  2584       }

  2586       std::vector<Node> order;

  2588       while (!heap.empty()) {

  2589         Node n = heap.top();

  2590         heap.pop();

  2591         _color_map[n] = -1;

  2592         order.push_back(n);

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

  2594           Node t = _graph.runningNode(e);

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

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

  2597           }

  2598         }

  2599       }

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

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

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

  2604           Node t = _graph.runningNode(e);

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

  2606             forbidden[_color_map[t]] = true;

  2607           }

  2608         }

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

  2610           if (!forbidden[k]) {

  2611             _color_map[order[i]] = k;

  2612             break;

  2613           }

  2614         }

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

  2616           return false;

  2617         }

  2618       }

  2619       return true;

  2620     }

  2622   private:

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

  2625       int ucolor = _color_map[u];

  2626       int vcolor = _color_map[v];

  2627       typedef _planarity_bits::KempeFilter<IndexMap> KempeFilter;

  2628       KempeFilter filter(_color_map, ucolor, vcolor);

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

  2631       KempeGraph kempe_graph(_graph, filter);

  2633       std::vector<Node> comp;

  2634       Bfs<KempeGraph> bfs(kempe_graph);

  2635       bfs.init();

  2636       bfs.addSource(u);

  2637       while (!bfs.emptyQueue()) {

  2638         Node n = bfs.nextNode();

  2639         if (n == v) return false;

  2640         comp.push_back(n);

  2641         bfs.processNextNode();

  2642       }

  2644       int scolor = ucolor + vcolor;

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

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

  2647       }

  2649       return true;

  2650     }

  2652     template <typename EmbeddingMap>

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

  2654       std::vector<Node> nodes;

  2655       nodes.reserve(4);

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

  2658         Node t = _graph.target(e);

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

  2660           nodes.push_back(t);

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

  2662         }

  2663       }

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

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

  2667         _color_map[node] = color;

  2668       } else {

  2669         color = _color_map[nodes[1]];

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

  2671         _color_map[node] = color;

  2672       }

  2673     }

  2675   public:

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

  2678     ///

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

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

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

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

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

  2684     /// It can be computed using PlanarEmbedding.

  2685     template <typename EmbeddingMap>

  2686     void runFiveColoring(const EmbeddingMap& embedding) {

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

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

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

  2692         _color_map[n] = -2;

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

  2694       }

  2696       std::vector<Node> order;

  2698       while (!heap.empty()) {

  2699         Node n = heap.top();

  2700         heap.pop();

  2701         _color_map[n] = -1;

  2702         order.push_back(n);

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

  2704           Node t = _graph.runningNode(e);

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

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

  2707           }

  2708         }

  2709       }

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

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

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

  2714           Node t = _graph.runningNode(e);

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

  2716             forbidden[_color_map[t]] = true;

  2717           }

  2718         }

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

  2720           if (!forbidden[k]) {

  2721             _color_map[order[i]] = k;

  2722             break;

  2723           }

  2724         }

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

  2726           kempeRecoloring(order[i], embedding);

  2727         }

  2728       }

  2729     }

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

  2732     ///

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

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

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

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

  2737     bool runFiveColoring() {

  2738       PlanarEmbedding<Graph> pe(_graph);

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

  2741       runFiveColoring(pe.embeddingMap());

  2742       return true;

  2743     }

  2745   private:

  2747     const Graph& _graph;

  2748     IndexMap _color_map;

  2749     Palette _palette;

  2750   };

  2752 }

  2754 #endif

author	Balazs Dezso <deba@inf.elte.hu>
	Sun, 14 Nov 2010 16:35:31 +0100
changeset 1018	2e959a5a0c2d
parent 877	141f9c0db4a3
child 1092	dceba191c00d
permissions	-rw-r--r--