COIN-OR::LEMON - Graph Library

source: lemon-main/lemon/planarity.h

Last change on this file was 1182:6b79d93e812f, checked in by Balazs Dezso <deba@…>, 13 months ago

Planar drawing algorithm now works for less than 3 nodes (#611)

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