Search code examples
c++boostconvex-hullboost-geometry

Boost convex hull with longitude and latitude


I am trying to use boost's convex_hull algorithm with longitude/latitude coordinates.

From here: https://archive.fosdem.org/2017/schedule/event/geo_boost_geography/attachments/slides/1748/export/events/attachments/geo_boost_geography/slides/1748/FOSDEM17_vissarion.pdf

I can see that we can calculate the distance between two points and even find the area using longitude/latitude coordinates (see page 19 and 22 of the PDF document).

Combining that with https://www.boost.org/doc/libs/1_75_0/libs/geometry/doc/html/geometry/reference/algorithms/convex_hull.html

I came up with this: https://wandbox.org/permlink/2AGPUtHPWrlGFMTf, but it does not compile, code here for convenience:

#include <iostream>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/adapted/boost_tuple.hpp>

namespace bg = boost::geometry;

int main()
{
    typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree>> point;
    typedef boost::geometry::model::polygon<point> polygon;
    polygon poly;
    bg::read_wkt(" POLYGON ((4.346693 50.858306, 4.367945 50.852455, 4.366227 50.840809, 4.344961 50.833264, 4.338074 50.848677,4.346693 50.858306))",
    poly );
    
    polygon hull;
    boost::geometry::convex_hull(poly, hull);

    using boost::geometry::dsv;
    std::cout << "polygon: " << dsv(poly) << std::endl << "hull: " << dsv(hull) << std::endl;
}

Any help is much appreciated.


Solution

  • Yeah, though you may be right that a strategy /can/ be made, that is not implemented.

    A little side-by-side tester clearly shows that the strategy is not implemented for the geographical coordinate system:

    template <typename cs> void test() {
        using point   = bg::model::point<double, 2, cs>;
        using polygon = bg::model::polygon<point>;
        polygon poly;
        bg::read_wkt("POLYGON((4.346693 50.858306, 4.367945 50.852455, 4.366227 "
                     "50.840809, 4.344961 50.833264, 4.338074 50.848677,4.346693 "
                     "50.858306))",
                     poly);
    
        std::cout << std::fixed;
        std::cout << "Polygon:   " << bg::dsv(poly)       << std::endl;
        std::cout << "Perimeter: " << bg::perimeter(poly) << std::endl;
        std::cout << "Area:      " << bg::area(poly)      << std::endl;
    
        using Strategy = typename bg::strategy_convex_hull<polygon, point>::type;
        std::cout << "Strategy " << boost::core::demangle(typeid(Strategy).name()) << "\n";
    
        if constexpr (not std::is_same_v<Strategy, bg::strategy::not_implemented>) {
            polygon hull;
            bg::convex_hull(poly, hull);
            std::cout << "Hull: " << bg::dsv(hull) << std::endl;
        }
    }
    

    It also demonstrates that some other strategies are implemented (e.g. distance).

    See it Live On Coliru

    int main() {
        std::cout << "Cartesian:\n";
        std::cout << "----------\n";
        test<bg::cs::cartesian>();
    
        std::cout << "\nGeographic:\n";
        std::cout <<   "-----------\n";
        test<bg::cs::geographic<bg::degree>>();
    }
    

    Simplifying the the typenames in the output:

        Cartesian:
        ----------
        Polygon:   (((4.346693, 50.858306), (4.367945, 50.852455), (4.366227, 50.840809), (4.344961, 50.833264), (4.338074, 50.848677), (4.346693, 50.858306)))
        Perimeter: 0.086184
        Area:      0.000488
        Strategy bg::strategy::convex_hull::graham_andrew<polygon, point>
        Hull: (((4.338074, 50.848677), (4.346693, 50.858306), (4.367945, 50.852455), (4.366227, 50.840809), (4.344961, 50.833264), (4.338074, 50.848677)))
    
        Geographic:
        -----------
        Polygon:   (((4.346693, 50.858306), (4.367945, 50.852455), (4.366227, 50.840809), (4.344961, 50.833264), (4.338074, 50.848677), (4.346693, 50.858306)))
        Perimeter: 7663.398262
        Area:      3848183.734567
        Strategy bg::strategy::not_implemented
    

    A look at the documented strategies suggests that graham_andrew is in fact the only one available.

    You should probably find out what tweaks are required to get things to work. it is technically possible to force convex_hull to use the Graham/Andrew strategy, but that seems ill-advised as the trait implies that the strategy is specifically deselected based on the coordinate system:

    /*!
        \brief Traits class binding a convex hull calculation strategy to a coordinate system
        \ingroup convex_hull
        \tparam Tag tag of coordinate system
        \tparam Geometry the geometry type (hull operates internally per hull over geometry)
        \tparam Point point-type of output points
    */
    template
    <
        typename Geometry1,
        typename Point,
        typename CsTag = typename cs_tag<Point>::type
    >
    struct strategy_convex_hull
    {
        typedef strategy::not_implemented type;
    };
    

    Digging into the implementation of the strategy here's a hopeful hint:

        // TODO: User-defiend CS-specific side strategy
        typename strategy::side::services::default_strategy<cs_tag>::type side;
    

    Perhaps we could be "done" with "just" specializing the Side Strategy for your coordinate systems? And more interestingly: a strategy::side::geographic exists. I'm out of my depth understanding the parameters (e.g. what the geodetic solution policy means?), but maybe yourself can take it from there?

    I'm convinced that if you know what needs to be done, the helpful devs over at the mailing list will be very willing to guide the technical questions on how to fit it into the library.