COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/planarity.h @ 2553:bfced05fa852

Last change on this file since 2553:bfced05fa852 was 2553:bfced05fa852, checked in by Alpar Juttner, 12 years ago

Happy New Year to LEMON (+ better update-copyright-header script)

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