COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/planarity.h @ 2521:05c0ba99cc27

Last change on this file since 2521:05c0ba99cc27 was 2508:c86db0f7f917, checked in by Balazs Dezso, 16 years ago

Planar graph coloring

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