Sample Code for Euclidean MST using Boost and CGAL

From Wikipedia, the free encyclopedia

The Boost library and CGAL library together provide the building blocks needed to compute Euclidean MST using Delauney Triangulation and Prim MST

The code is based on the examples in the docs for CGAL and Boost

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_euclidean_traits_xy_3.h>
#include <CGAL/Delaunay_triangulation_2.h>
 
#include <stdio.h>
#include <fstream>
#include <cstdlib>
 
#include <boost/config.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/prim_minimum_spanning_tree.hpp>

struct K : CGAL::Exact_predicates_inexact_constructions_kernel {};

typedef CGAL::Triangulation_euclidean_traits_xy_3<K>  Gt;
typedef CGAL::Delaunay_triangulation_2<Gt> Delaunay;

typedef K::Point_3   Point;

typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS,
    boost::property<boost::vertex_distance_t, int>, boost::property < boost::edge_weight_t, double > > Graph;
typedef std::pair < int, int > Edge;


double euclid_dist(Point pt1,
                   Point pt2)
{
  double x = pt1.x()-pt2.x();
  double y = pt1.y()-pt2.y();
  return sqrt(x*x+y*y);
}

void
fast_mst(std::vector<Point> &nodes)
{
  std::cout << "fast mst pts\n";

  std::map<Point, int> pt_vert_map;
  
  Delaunay dt;
  
  int nodeid = 0;
  std::vector<Point>::iterator ptiter;
  for(ptiter=nodes.begin(); ptiter!=nodes.end(); ++ptiter) {
    Point pt = *ptiter;
    pt_vert_map[pt] = nodeid;
    nodeid++;
    dt.insert(pt);
  } 
  
  std::vector<Edge> edges;
  std::vector<double> weights;
  Delaunay::Finite_edges_iterator iter;
  for(iter = dt.finite_edges_begin(); iter!=dt.finite_edges_end(); iter++) {
    Delaunay::Edge pp = *iter;
    Delaunay::Face_handle f = pp.first;
    int fedge = pp.second;

    Point pt1 = f->vertex((fedge+1)%3)->point();
    Point pt2 = f->vertex((fedge+2)%3)->point();

    int id1 = pt_vert_map[pt1];
    int id2 = pt_vert_map[pt2];

    edges.push_back(Edge(id1,id2));
    weights.push_back(euclid_dist(pt1, pt2));
  }

  Graph g(edges.begin(), edges.end(), weights.begin(), nodes.size());
  boost::property_map<Graph, boost::edge_weight_t>::type weightmap = get(boost::edge_weight, g);

  std::vector < boost::graph_traits < Graph >::vertex_descriptor >
    p(num_vertices(g));

  prim_minimum_spanning_tree(g, &p[0]);

  for (std::size_t i = 0; i != p.size(); ++i) {
    if (p[i] != i) {
      std::cout << "edge between point: " << nodes[i] << "(" << i << " and point: ";
      std::cout << nodes[p[i]] << "(" << p[i] << ")" << std::endl;
    } 
  }

}

int main()
{
  std::vector<Point> nodes;
  for(int i=0; i<100; i++){
    Point pt(double(rand()%1000), double(rand()%1000), 0.0);
    std::cout << "pt: " << pt << std::endl;
    nodes.push_back(pt);
  } /* end of for(int i=0; i<100; i++) */
  
  fast_mst(nodes);
 }