/Src/Dependencies/Boost/libs/geometry/example/07_a_graph_route_example.cpp

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