COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/planarity.h @ 875:07ec2b52e53d

Last change on this file since 875:07ec2b52e53d was 828:5fd7fafc4470, checked in by Peter Kovacs <kpeter@…>, 14 years ago

Doc improvements for planarity related tools (#62)

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