COIN-OR::LEMON - Graph Library

Changeset 2447:260ce674cc65 in lemon-0.x for tools


Ignore:
Timestamp:
06/05/07 13:49:19 (17 years ago)
Author:
Balazs Dezso
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3284
Message:

Delaunay triangulation
Faster geometric minimum spanning tree

File:
1 edited

Legend:

Unmodified
Added
Removed
  • tools/lgf-gen.cc

    r2446 r2447  
    3030#include <cmath>
    3131#include <algorithm>
    32 #include <lemon/unionfind.h>
     32#include <lemon/kruskal.h>
    3333#include <lemon/time_measure.h>
    3434
     
    147147std::vector<UEdge> edges;
    148148
    149 void triangle()
    150 {
     149namespace _delaunay_bits {
     150
     151  struct Part {
     152    int prev, curr, next;
     153
     154    Part(int p, int c, int n) : prev(p), curr(c), next(n) {}
     155  };
     156
     157  inline std::ostream& operator<<(std::ostream& os, const Part& part) {
     158    os << '(' << part.prev << ',' << part.curr << ',' << part.next << ')';
     159    return os;
     160  }
     161
     162  inline double circle_point(const Point& p, const Point& q, const Point& r) {
     163    double a = p.x * (q.y - r.y) + q.x * (r.y - p.y) + r.x * (p.y - q.y);
     164    if (a == 0) return std::numeric_limits<double>::quiet_NaN();
     165
     166    double d = (p.x * p.x + p.y * p.y) * (q.y - r.y) +
     167      (q.x * q.x + q.y * q.y) * (r.y - p.y) +
     168      (r.x * r.x + r.y * r.y) * (p.y - q.y);
     169
     170    double e = (p.x * p.x + p.y * p.y) * (q.x - r.x) +
     171      (q.x * q.x + q.y * q.y) * (r.x - p.x) +
     172      (r.x * r.x + r.y * r.y) * (p.x - q.x);
     173
     174    double f = (p.x * p.x + p.y * p.y) * (q.x * r.y - r.x * q.y) +
     175      (q.x * q.x + q.y * q.y) * (r.x * p.y - p.x * r.y) +
     176      (r.x * r.x + r.y * r.y) * (p.x * q.y - q.x * p.y);
     177
     178    return d / (2 * a) + sqrt((d * d + e * e) / (4 * a * a) + f / a);
     179  }
     180
     181  inline bool circle_form(const Point& p, const Point& q, const Point& r) {
     182    return rot90(q - p) * (r - q) < 0.0;
     183  }
     184
     185  inline double intersection(const Point& p, const Point& q, double sx) {
     186    const double epsilon = 1e-8;
     187
     188    if (p.x == q.x) return (p.y + q.y) / 2.0;
     189
     190    if (sx < p.x + epsilon) return p.y;
     191    if (sx < q.x + epsilon) return q.y;
     192   
     193    double a = q.x - p.x;
     194    double b = (q.x - sx) * p.y - (p.x - sx) * q.y;   
     195    double d = (q.x - sx) * (p.x - sx) * (p - q).normSquare();
     196    return (b - sqrt(d)) / a;
     197  }
     198
     199  struct YLess {
     200
     201
     202    YLess(const std::vector<Point>& points, double& sweep)
     203      : _points(points), _sweep(sweep) {}
     204
     205    bool operator()(const Part& l, const Part& r) const {
     206      const double epsilon = 1e-8;
     207
     208      //      std::cerr << l << " vs " << r << std::endl;
     209      double lbx = l.prev != -1 ?
     210        intersection(_points[l.prev], _points[l.curr], _sweep) :
     211        - std::numeric_limits<double>::infinity();
     212      double rbx = r.prev != -1 ?
     213        intersection(_points[r.prev], _points[r.curr], _sweep) :
     214        - std::numeric_limits<double>::infinity();
     215      double lex = l.next != -1 ?
     216        intersection(_points[l.curr], _points[l.next], _sweep) :
     217        std::numeric_limits<double>::infinity();
     218      double rex = r.next != -1 ?
     219        intersection(_points[r.curr], _points[r.next], _sweep) :
     220        std::numeric_limits<double>::infinity();
     221
     222      if (lbx > lex) std::swap(lbx, lex);
     223      if (rbx > rex) std::swap(rbx, rex);
     224
     225      if (lex < epsilon + rex && lbx + epsilon < rex) return true;
     226      if (rex < epsilon + lex && rbx + epsilon < lex) return false;
     227      return lex < rex;
     228    }
     229   
     230    const std::vector<Point>& _points;
     231    double& _sweep;
     232  };
     233 
     234  struct BeachIt;
     235 
     236  typedef std::multimap<double, BeachIt> SpikeHeap;
     237
     238  typedef std::multimap<Part, SpikeHeap::iterator, YLess> Beach;
     239
     240  struct BeachIt {
     241    Beach::iterator it;
     242
     243    BeachIt(Beach::iterator iter) : it(iter) {}
     244  };
     245
     246}
     247
     248inline void delaunay() {
    151249  Counter cnt("Number of edges added: ");
    152   std::vector<Pedge> pedges;
    153   for(NodeIt n(g);n!=INVALID;++n)
    154     for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
     250 
     251  using namespace _delaunay_bits;
     252
     253  typedef _delaunay_bits::Part Part;
     254  typedef std::vector<std::pair<double, int> > SiteHeap;
     255
     256
     257  std::vector<Point> points;
     258  std::vector<Node> nodes;
     259
     260  for (NodeIt it(g); it != INVALID; ++it) {
     261    nodes.push_back(it);
     262    points.push_back(coords[it]);
     263  }
     264
     265  SiteHeap siteheap(points.size());
     266
     267  double sweep;
     268
     269
     270  for (int i = 0; i < int(siteheap.size()); ++i) {
     271    siteheap[i] = std::make_pair(points[i].x, i);
     272  }
     273 
     274  std::sort(siteheap.begin(), siteheap.end());
     275  sweep = siteheap.front().first;
     276 
     277  YLess yless(points, sweep);
     278  Beach beach(yless);
     279
     280  SpikeHeap spikeheap;
     281
     282  std::set<std::pair<int, int> > edges;
     283
     284  beach.insert(std::make_pair(Part(-1, siteheap[0].second, -1),
     285                              spikeheap.end()));
     286  int siteindex = 1;
     287
     288  while (siteindex < int(points.size()) || !spikeheap.empty()) {
     289
     290    SpikeHeap::iterator spit = spikeheap.begin();
     291
     292    if (siteindex < int(points.size()) &&
     293        (spit == spikeheap.end() || siteheap[siteindex].first < spit->first)) {
     294      int site = siteheap[siteindex].second;
     295      sweep = siteheap[siteindex].first;
     296         
     297      Beach::iterator bit = beach.upper_bound(Part(site, site, site));
     298     
     299      if (bit->second != spikeheap.end()) {
     300        spikeheap.erase(bit->second);   
     301      }
     302
     303      int prev = bit->first.prev;
     304      int curr = bit->first.curr;
     305      int next = bit->first.next;
     306
     307      beach.erase(bit);
     308     
     309      SpikeHeap::iterator pit = spikeheap.end();
     310      if (prev != -1 &&
     311          circle_form(points[prev], points[curr], points[site])) {
     312        double x = circle_point(points[prev], points[curr], points[site]);
     313        pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
     314        pit->second.it =
     315          beach.insert(std::make_pair(Part(prev, curr, site), pit));
     316      } else {
     317        beach.insert(std::make_pair(Part(prev, curr, site), pit));
     318      }
     319
     320      beach.insert(std::make_pair(Part(curr, site, curr), spikeheap.end()));
     321     
     322      SpikeHeap::iterator nit = spikeheap.end();
     323      if (next != -1 &&
     324          circle_form(points[site], points[curr],points[next])) {
     325        double x = circle_point(points[site], points[curr], points[next]);
     326        nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
     327        nit->second.it =
     328          beach.insert(std::make_pair(Part(site, curr, next), nit));
     329      } else {
     330        beach.insert(std::make_pair(Part(site, curr, next), nit));
     331      }
     332     
     333      ++siteindex;
     334    } else {
     335      sweep = spit->first;     
     336
     337      Beach::iterator bit = spit->second.it;
     338
     339      int prev = bit->first.prev;
     340      int curr = bit->first.curr;
     341      int next = bit->first.next;
     342
    155343      {
    156         Pedge p;
    157         p.a=n;
    158         p.b=m;
    159         p.len=(coords[m]-coords[n]).normSquare();
    160         pedges.push_back(p);
    161       }
    162   std::sort(pedges.begin(),pedges.end(),pedgeLess);
    163   for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
    164     {
    165       Line li(pi->a,pi->b);
    166       UEdgeIt e(g);
    167       for(;e!=INVALID && !cross(e,li);++e) ;
    168       UEdge ne;
    169       if(e==INVALID) {
    170         ne=g.addEdge(pi->a,pi->b);
    171         edges.push_back(ne);
    172         cnt++;
    173       }
    174     }
     344        std::pair<int, int> edge;
     345
     346        edge = prev < curr ?
     347          std::make_pair(prev, curr) : std::make_pair(curr, prev);
     348       
     349        if (edges.find(edge) == edges.end()) {
     350          edges.insert(edge);
     351          g.addEdge(nodes[prev], nodes[curr]);
     352          ++cnt;
     353        }
     354
     355        edge = curr < next ?
     356          std::make_pair(curr, next) : std::make_pair(next, curr);
     357       
     358        if (edges.find(edge) == edges.end()) {
     359          edges.insert(edge);
     360          g.addEdge(nodes[curr], nodes[next]);
     361          ++cnt;
     362        }
     363      }
     364     
     365      Beach::iterator pbit = bit; --pbit;
     366      int ppv = pbit->first.prev;
     367      Beach::iterator nbit = bit; ++nbit;
     368      int nnt = nbit->first.next;
     369
     370      if (bit->second != spikeheap.end()) spikeheap.erase(bit->second);
     371      if (pbit->second != spikeheap.end()) spikeheap.erase(pbit->second);
     372      if (nbit->second != spikeheap.end()) spikeheap.erase(nbit->second);
     373
     374      beach.erase(nbit);
     375      beach.erase(bit);
     376      beach.erase(pbit);
     377
     378      SpikeHeap::iterator pit = spikeheap.end();
     379      if (ppv != -1 && ppv != next &&
     380          circle_form(points[ppv], points[prev], points[next])) {
     381        double x = circle_point(points[ppv], points[prev], points[next]);
     382        if (x < sweep) x = sweep;
     383        pit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
     384        pit->second.it =
     385          beach.insert(std::make_pair(Part(ppv, prev, next), pit));
     386      } else {
     387        beach.insert(std::make_pair(Part(ppv, prev, next), pit));
     388      }
     389
     390      SpikeHeap::iterator nit = spikeheap.end();
     391      if (nnt != -1 && prev != nnt &&
     392          circle_form(points[prev], points[next], points[nnt])) {
     393        double x = circle_point(points[prev], points[next], points[nnt]);
     394        if (x < sweep) x = sweep;
     395        nit = spikeheap.insert(std::make_pair(x, BeachIt(beach.end())));
     396        nit->second.it =
     397          beach.insert(std::make_pair(Part(prev, next, nnt), nit));
     398      } else {
     399        beach.insert(std::make_pair(Part(prev, next, nnt), nit));
     400      }
     401     
     402    }
     403  }
     404
     405  for (Beach::iterator it = beach.begin(); it != beach.end(); ++it) {
     406    int curr = it->first.curr;
     407    int next = it->first.next;
     408
     409    if (next == -1) continue;
     410
     411    std::pair<int, int> edge;
     412
     413    edge = curr < next ?
     414      std::make_pair(curr, next) : std::make_pair(next, curr);
     415   
     416    if (edges.find(edge) == edges.end()) {
     417      edges.insert(edge);
     418      g.addEdge(nodes[curr], nodes[next]);
     419      ++cnt;
     420    }
     421  }
    175422}
    176423
     
    246493}
    247494
     495template <typename UGraph, typename CoordMap>
     496class LengthSquareMap {
     497public:
     498  typedef typename UGraph::UEdge Key;
     499  typedef typename CoordMap::Value::Value Value;
     500
     501  LengthSquareMap(const UGraph& ugraph, const CoordMap& coords)
     502    : _ugraph(ugraph), _coords(coords) {}
     503
     504  Value operator[](const Key& key) const {
     505    return (_coords[_ugraph.target(key)] -
     506            _coords[_ugraph.source(key)]).normSquare();
     507  }
     508
     509private:
     510
     511  const UGraph& _ugraph;
     512  const CoordMap& _coords;
     513};
     514
    248515void minTree() {
    249   int en=0;
    250   int pr=0;
    251516  std::vector<Pedge> pedges;
    252517  Timer T;
    253   std::cout << T.realTime() << "s: Setting up the edges...\n";
    254   for(NodeIt n(g);n!=INVALID;++n)
    255     {
    256       for(NodeIt m=++(NodeIt(n));m!=INVALID;++m)
    257         {
    258           Pedge p;
    259           p.a=n;
    260           p.b=m;
    261           p.len=(coords[m]-coords[n]).normSquare();
    262           pedges.push_back(p);
    263         }
    264       if(progress && en>=pr*double(N)/100)
    265         {
    266           std::cout << pr << "%  \r" << std::flush;
    267           pr++;
    268         }
    269     }
    270   std::cout << T.realTime() << "s: Sorting the edges...\n";
    271   std::sort(pedges.begin(),pedges.end(),pedgeLess);
    272   std::cout << T.realTime() << "s: Creating the tree...\n";
    273   ListUGraph::NodeMap<int> comp(g);
    274   UnionFind<ListUGraph::NodeMap<int> > uf(comp);
    275   for (NodeIt it(g); it != INVALID; ++it)
    276     uf.insert(it); 
    277   for(std::vector<Pedge>::iterator pi=pedges.begin();pi!=pedges.end();++pi)
    278     {
    279       if ( uf.join(pi->a,pi->b) ) {
    280         g.addEdge(pi->a,pi->b);
    281         en++;
    282         if(en>=N-1)
    283           break;
    284       }
    285     }
     518  std::cout << T.realTime() << "s: Creating delaunay triangulation...\n";
     519  delaunay();
     520  std::cout << T.realTime() << "s: Calculating spanning tree...\n";
     521  LengthSquareMap<ListUGraph, ListUGraph::NodeMap<Point> > ls(g, coords);
     522  ListUGraph::UEdgeMap<bool> tree(g);
     523  kruskal(g, ls, tree);
     524  std::cout << T.realTime() << "s: Removing non tree edges...\n";
     525  std::vector<UEdge> remove;
     526  for (UEdgeIt e(g); e != INVALID; ++e) {
     527    if (!tree[e]) remove.push_back(e);
     528  }
     529  for(int i = 0; i < int(remove.size()); ++i) {
     530    g.erase(remove[i]);
     531  }
    286532  std::cout << T.realTime() << "s: Done\n";
    287533}
     
    409655    .boolOption("tsp2", "Create a TSP tour (tree based)")
    410656    .optionGroup("alg","tsp2")
     657    .boolOption("dela", "Delaunay triangulation graph")
     658    .optionGroup("alg","dela")
    411659    .onlyOneGroup("alg")
     660    .boolOption("rand", "Use time seed for random number generator")
     661    .optionGroup("rand", "rand")
     662    .intOption("seed", "Random seed", -1)
     663    .optionGroup("rand", "seed")
     664    .onlyOneGroup("rand")
    412665    .other("[prefix]","Prefix of the output files. Default is 'lgf-gen-out'")
    413666    .run();
     667
     668  if (ap["rand"]) {
     669    int seed = time(0);
     670    std::cout << "Random number seed: " << seed << std::endl;
     671    rnd = Random(seed);
     672  }
     673  if (ap.given("seed")) {
     674    int seed = ap["seed"];
     675    std::cout << "Random number seed: " << seed << std::endl;
     676    rnd = Random(seed);
     677  }
    414678 
    415679  std::string prefix;
     
    464728        }
    465729    }
     730
     731//   for (ListUGraph::NodeIt n(g); n != INVALID; ++n) {
     732//     std::cerr << coords[n] << std::endl;
     733//   }
    466734 
    467735  if(ap["tsp"]) {
     
    482750  else if(ap["tree"]) {
    483751    minTree();
     752  }
     753  else if(ap["dela"]) {
     754    delaunay();
    484755  }
    485756 
Note: See TracChangeset for help on using the changeset viewer.