libs/geometry/example/07_b_graph_route_example.cpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
   4 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
   5 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
   6
   7 // Use, modification and distribution is subject to the Boost Software License,
   8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
   9 // http://www.boost.org/LICENSE_1_0.txt)
  10 //
  11 // Example showing Boost.Geometry combined with Boost.Graph, calculating shortest routes
  12 // input: two WKT's, provided in subfolder data
  13 // output: text, + an SVG, displayable in e.g. Firefox)
  14
  15
  16 // Second variant, with bundled properties
  17
  18
  19 #include <iostream>
  20 #include <fstream>
  21 #include <iomanip>
  22 #include <limits>
  23
  24 #include <boost/tuple/tuple.hpp>
  25 #include <boost/foreach.hpp>
  26
  27 #include <boost/graph/adjacency_list.hpp>
  28 #include <boost/graph/dijkstra_shortest_paths.hpp>
  29
  30 #include <boost/geometry/geometry.hpp>
  31 #include <boost/geometry/geometries/linestring.hpp>
  32 #include <boost/geometry/io/wkt/read.hpp>
  33
  34 // For output:
  35 #include <boost/geometry/io/svg/svg_mapper.hpp>
  36
  37 // For distance-calculations over the Earth:
  38 //#include <boost/geometry/extensions/gis/geographic/strategies/andoyer.hpp>
  39
  40
  41
  42 // Read an ASCII file containing WKT's, fill a vector of tuples
  43 // The tuples consist of at least <0> a geometry and <1> an identifying string
  44 template <typename Geometry, typename Tuple, typename Box>
  45 void read_wkt(std::string const& filename, std::vector<Tuple>& tuples, Box& box)
  46 {
  47     std::ifstream cpp_file(filename.c_str());
  48     if (cpp_file.is_open())
  49     {
  50         while (! cpp_file.eof() )
  51         {
  52             std::string line;
  53             std::getline(cpp_file, line);
  54             Geometry geometry;
  55             boost::trim(line);
  56             if (! line.empty() && ! boost::starts_with(line, "#"))
  57             {
  58                 std::string name;
  59
  60                 // Split at ';', if any
  61                 std::string::size_type pos = line.find(";");
  62                 if (pos != std::string::npos)
  63                 {
  64                     name = line.substr(pos + 1);
  65                     line.erase(pos);
  66
  67                     boost::trim(line);
  68                     boost::trim(name);
  69                 }
  70
  71                 Geometry geometry;
  72                 boost::geometry::read_wkt(line, geometry);
  73
  74                 Tuple tuple(geometry, name);
  75
  76                 tuples.push_back(tuple);
  77                 boost::geometry::expand(box, boost::geometry::return_envelope<Box>(geometry));
  78             }
  79         }
  80     }
  81 }
  82
  83 // To calculate distance, declare and construct a strategy with average earth radius
  84 boost::geometry::strategy::distance::haversine<double> const haversine(6372795.0);
  85
  86 // Define properties for vertex
  87 template <typename Point>
  88 struct bg_vertex_property
  89 {
  90     bg_vertex_property()
  91     {
  92         boost::geometry::assign_zero(location);
  93     }
  94     bg_vertex_property(Point const& loc)
  95         : location(loc)
  96     {
  97     }
  98
  99     Point location;
 100 };
 101
 102 // Define properties for edge
 103 template <typename Linestring>
 104 struct bg_edge_property
 105 {
 106     bg_edge_property(Linestring const& line)
 107         : length(boost::geometry::length(line, haversine))
 108         , m_line(line)
 109     {
 110     }
 111
 112     inline Linestring const& line() const
 113     {
 114         return m_line;
 115     }
 116
 117     double length;
 118 private :
 119     Linestring m_line;
 120 };
 121
 122 // Utility function to add a vertex to a graph. It might exist already. Then do not insert,
 123 // but return vertex descriptor back. It might not exist. Then add it (and return).
 124 // To efficiently handle this, a std::map is used.
 125 template <typename M, typename K, typename G>
 126 inline typename boost::graph_traits<G>::vertex_descriptor find_or_insert(M& map, K const& key, G& graph)
 127 {
 128     typename M::const_iterator it = map.find(key);
 129     if (it == map.end())
 130     {
 131         // Add a vertex to the graph
 132         typename boost::graph_traits<G>::vertex_descriptor new_vertex
 133             = boost::add_vertex(bg_vertex_property<typename M::key_type>(key), graph);
 134
 135         // Add to the map, using POINT as key
 136         map[key] = new_vertex;
 137         return new_vertex;
 138     }
 139     return it->second;
 140 }
 141
 142 template
 143 <
 144     typename Graph,
 145     typename RoadTupleVector,
 146     typename CityTupleVector
 147 >
 148 void add_roads_and_connect_cities(Graph& graph,
 149             RoadTupleVector const& roads,
 150             CityTupleVector& cities)
 151 {
 152     typedef typename boost::range_value<RoadTupleVector>::type road_type;
 153     typedef typename boost::tuples::element<0, road_type>::type line_type;
 154     typedef typename boost::geometry::point_type<line_type>::type point_type;
 155
 156     typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_type;
 157
 158     // Define a map to be used during graph filling
 159     // Maps from point to vertex-id's
 160     typedef std::map<point_type, vertex_type, boost::geometry::less<point_type> > map_type;
 161     map_type map;
 162
 163
 164     // Fill the graph
 165     BOOST_FOREACH(road_type const& road, roads)
 166     {
 167         line_type const& line = road.template get<0>();
 168         // Find or add begin/end point of these line
 169         vertex_type from = find_or_insert(map, line.front(), graph);
 170         vertex_type to = find_or_insert(map, line.back(), graph);
 171         boost::add_edge(from, to, bg_edge_property<line_type>(line), graph);
 172     }
 173
 174     // Find nearest graph vertex for each city, using the map
 175     typedef typename boost::range_value<CityTupleVector>::type city_type;
 176     BOOST_FOREACH(city_type& city, cities)
 177     {
 178         double min_distance = 1e300;
 179         for(typename map_type::const_iterator it = map.begin(); it != map.end(); ++it)
 180         {
 181             double dist = boost::geometry::distance(it->first, city.template get<0>());
 182             if (dist < min_distance)
 183             {
 184                 min_distance = dist;
 185                 // Set the vertex
 186                 city.template get<2>() = it->second;
 187             }
 188         }
 189     }
 190 }
 191
 192 template <typename Graph, typename Route>
 193 inline void add_edge_to_route(Graph const& graph,
 194             typename boost::graph_traits<Graph>::vertex_descriptor vertex1,
 195             typename boost::graph_traits<Graph>::vertex_descriptor vertex2,
 196             Route& route)
 197 {
 198     std::pair
 199         <
 200             typename boost::graph_traits<Graph>::edge_descriptor,
 201             bool
 202         > opt_edge = boost::edge(vertex1, vertex2, graph);
 203     if (opt_edge.second)
 204     {
 205         // Get properties of edge and of vertex
 206         bg_edge_property<Route> const& edge_prop = graph[opt_edge.first];
 207         bg_vertex_property
 208             <
 209                 typename boost::geometry::point_type<Route>::type
 210             > const& vertex_prop = graph[vertex2];
 211
 212         // Depending on how edge connects to vertex, copy it forward or backward
 213         if (boost::geometry::equals(edge_prop.line().front(), vertex_prop.location))
 214         {
 215             std::copy(edge_prop.line().begin(), edge_prop.line().end(),
 216                 std::back_inserter(route));
 217         }
 218         else
 219         {
 220             std::reverse_copy(edge_prop.line().begin(), edge_prop.line().end(),
 221                 std::back_inserter(route));
 222         }
 223     }
 224 }
 225
 226
 227 template <typename Graph, typename Route>
 228 inline void build_route(Graph const& graph,
 229             std::vector<typename boost::graph_traits<Graph>::vertex_descriptor> const& predecessors,
 230             typename boost::graph_traits<Graph>::vertex_descriptor vertex1,
 231             typename boost::graph_traits<Graph>::vertex_descriptor vertex2,
 232             Route& route)
 233 {
 234     typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex_type;
 235     vertex_type pred = predecessors[vertex2];
 236
 237     add_edge_to_route(graph, vertex2, pred, route);
 238     while (pred != vertex1)
 239     {
 240         add_edge_to_route(graph, predecessors[pred], pred, route);
 241         pred = predecessors[pred];
 242     }
 243 }
 244
 245
 246 int main()
 247 {
 248     // Define a point in the Geographic coordinate system (currently spherical-equatorial)
 249     // (geographic calculations are in an extension; for sample it makes no difference)
 250     typedef boost::geometry::model::point
 251         <
 252             double, 2, boost::geometry::cs::spherical_equatorial<boost::geometry::degree>
 253         > point_type;
 254
 255     typedef boost::geometry::model::linestring<point_type> line_type;
 256
 257     // Define the graph, lateron containing the road network
 258     // With this, specify bundled properties for vertex and edge,
 259     // as in http://www.boost.org/doc/libs/1_43_0/libs/graph/doc/bundles.html
 260     typedef boost::adjacency_list
 261         <
 262             boost::vecS, boost::vecS, boost::undirectedS
 263             , bg_vertex_property<point_type> // bundled
 264             , bg_edge_property<line_type>
 265         > graph_type;
 266
 267     typedef boost::graph_traits<graph_type>::vertex_descriptor vertex_type;
 268
 269
 270     // Init a bounding box, lateron used to define SVG map
 271     boost::geometry::model::box<point_type> box;
 272     boost::geometry::assign_inverse(box);
 273
 274     graph_type graph;
 275
 276     // Read the cities
 277     typedef boost::tuple<point_type, std::string, vertex_type> city_type;
 278     std::vector<city_type> cities;
 279     read_wkt<point_type>("data/cities.wkt", cities, box);
 280
 281     // Read the road network
 282     typedef boost::tuple<line_type, std::string> road_type;
 283     std::vector<road_type> roads;
 284     read_wkt<line_type>("data/roads.wkt", roads, box);
 285
 286
 287     // Add roads and connect cities
 288     add_roads_and_connect_cities(graph, roads, cities);
 289
 290     double const km = 1000.0;
 291     std::cout << "distances, all in KM" << std::endl
 292         << std::fixed << std::setprecision(0);
 293
 294     // Main functionality: calculate shortest routes from/to all cities
 295
 296     // For the first one, the complete route is stored as a linestring
 297     bool first = true;
 298     line_type route;
 299
 300     int const n = boost::num_vertices(graph);
 301     BOOST_FOREACH(city_type const& city1, cities)
 302     {
 303         std::vector<vertex_type> predecessors(n);
 304         std::vector<double> costs(n);
 305
 306         // Call Dijkstra (without named-parameter to be compatible with all VC)
 307         boost::dijkstra_shortest_paths(graph, city1.get<2>(),
 308                 &predecessors[0], &costs[0],
 309                 boost::get(&bg_edge_property<line_type>::length, graph),
 310                 boost::get(boost::vertex_index, graph),
 311                 std::less<double>(), std::plus<double>(),
 312                 (std::numeric_limits<double>::max)(), double(),
 313                 boost::dijkstra_visitor<boost::null_visitor>());
 314
 315         BOOST_FOREACH(city_type const& city2, cities)
 316         {
 317             if (! boost::equals(city1.get<1>(), city2.get<1>()))
 318             {
 319                 double distance = costs[city2.get<2>()] / km;
 320                 double acof = boost::geometry::distance(city1.get<0>(), city2.get<0>(), haversine) / km;
 321
 322                 std::cout
 323                     << std::setiosflags (std::ios_base::left) << std::setw(15)
 324                         << city1.get<1>() << " - "
 325                     << std::setiosflags (std::ios_base::left) << std::setw(15)
 326                         << city2.get<1>()
 327                     << " -> through the air: " << std::setw(4) << acof
 328                     << " , over the road: " << std::setw(4) << distance
 329                     << std::endl;
 330
 331                 if (first)
 332                 {
 333                     build_route(graph, predecessors,
 334                             city1.get<2>(), city2.get<2>(),
 335                             route);
 336                     first = false;
 337                 }
 338             }
 339         }
 340     }
 341
 342 #if defined(HAVE_SVG)
 343     // Create the SVG
 344     std::ofstream stream("routes.svg");
 345     boost::geometry::svg_mapper<point_type> mapper(stream, 600, 600);
 346
 347     // Map roads
 348     BOOST_FOREACH(road_type const& road, roads)
 349     {
 350         mapper.add(road.get<0>());
 351     }
 352
 353     BOOST_FOREACH(road_type const& road, roads)
 354     {
 355         mapper.map(road.get<0>(),
 356                 "stroke:rgb(128,128,128);stroke-width:1");
 357     }
 358
 359     mapper.map(route,
 360             "stroke:rgb(0, 255, 0);stroke-width:6;opacity:0.5");
 361
 362     // Map cities
 363     BOOST_FOREACH(city_type const& city, cities)
 364     {
 365         mapper.map(city.get<0>(),
 366                 "fill:rgb(255,255,0);stroke:rgb(0,0,0);stroke-width:1");
 367         mapper.text(city.get<0>(), city.get<1>(),
 368                 "fill:rgb(0,0,0);font-family:Arial;font-size:10px", 5, 5);
 369     }
 370 #endif
 371
 372     return 0;
 373 }