Convex Hull of Polygon using Boost.Geometry

Rcpp can be used to convert basic R data types to and from Boost.Geometry models.

In this example we take a matrix of 2d-points and convert it into a Boost.Geometry polygon. We then compute the convex hull of this polygon using a Boost.Geometry function boost::geometry::convex_hull. The convex hull is then converted back to an R matrix.

The conversions to and from R and Boost.Geometry types are are done using two custom as() and wrap() convenience converter functions.

#include <Rcpp.h>
using namespace Rcpp;

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

BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian)

typedef boost::tuple<double, double> point;
typedef boost::geometry::model::polygon<point, true, true> polygon; 

namespace Rcpp {

    // `as` converter from R to Boost.Geometry's polygon type
    template <> polygon as(SEXP pointsMatrixSEXP) {
        // the coordinates are the rows of the (n x 2) matrix
        NumericMatrix pointsMatrix(pointsMatrixSEXP);
        polygon poly;
        for (int i = 0; i < pointsMatrix.nrow(); ++i) {
            double x = pointsMatrix(i,0);
            double y = pointsMatrix(i,1);
            point p(x,y);
            poly.outer().push_back(p); 
        }
        return (poly);
    } 
    
    // `wrap` converter from Boost.Geometry's polygon to an R(cpp) matrix
    // The Rcpp NumericMatrix can be converted to/from a SEXP
    template <> SEXP wrap(const polygon& poly) {
        const std::vector<point>& points = poly.outer();
        NumericMatrix rmat(points.size(), 2);
        for(int i = 0; i < points.size(); ++i) {
            const point& p = points[i];
            rmat(i,0) = p.get<0>();
            rmat(i,1) = p.get<1>();
        }
        return Rcpp::wrap(rmat);
    }
}


// [[Rcpp::export]]
NumericMatrix convexHullRcpp(SEXP pointsMatrixSEXP){

    // Conversion of pointsMatrix here to boost::geometry polygon
    polygon poly = as<polygon>(pointsMatrixSEXP);

    polygon hull;
    // Compute the convex hull
    boost::geometry::convex_hull(poly, hull);

    // Convert hull into a NumericMatrixsomething that Rcpp can hand back to R
    return wrap(hull);
}

Now we use R to replicate the convex_hull example in the Boost.Geometry documentation.

points <- c(2.0, 1.3, 2.4, 1.7, 2.8, 1.8, 3.4, 1.2, 3.7, 1.6, 3.4, 2.0, 4.1, 3.0, 5.3, 2.6, 5.4, 1.2, 4.9, 0.8, 2.9, 0.7, 2.0, 1.3)
points.matrix <- matrix(points, ncol=2, byrow=TRUE)
points.matrix


      [,1] [,2]
 [1,]  2.0  1.3
 [2,]  2.4  1.7
 [3,]  2.8  1.8
 [4,]  3.4  1.2
 [5,]  3.7  1.6
 [6,]  3.4  2.0
 [7,]  4.1  3.0
 [8,]  5.3  2.6
 [9,]  5.4  1.2
[10,]  4.9  0.8
[11,]  2.9  0.7
[12,]  2.0  1.3

convexHullRcpp(points.matrix)


     [,1] [,2]
[1,]  2.0  1.3
[2,]  2.4  1.7
[3,]  4.1  3.0
[4,]  5.3  2.6
[5,]  5.4  1.2
[6,]  4.9  0.8
[7,]  2.9  0.7
[8,]  2.0  1.3