COIN-OR::LEMON - Graph Library

source: lemon/lemon/planarity.h @ 862:58c330ad0b5c

Last change on this file since 862:58c330ad0b5c was 862:58c330ad0b5c, checked in by Balazs Dezso <deba@…>, 10 years ago

Planarity checking function instead of class (#62)

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