Convex Hull of Polygon using Boost.Geometry

Sameer D'Costa — written Feb 23, 2014 — source

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 the 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 which are implemented below, followed by a function using them in order to take data from R, call a Boost.Geometry function, and return the result to R.

#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 (unsigned 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 can 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

In fact, because of the interplay because providing as<>() and wrap() for the compiler, and how compileAttributes() works, we can write the function more cleanly with the desired new types in the interface.

// [[Rcpp::export]]
polygon convexHullRcpp2(polygon pointsMatrixBG){
    polygon hull;
    boost::geometry::convex_hull(pointsMatrixBG, hull);
    return hull;
}

We can run this second function as well:

convexHullRcpp2(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

tags: boost

Extending R with C++ and Fortran — Dirk Eddelbuettel and JBrandon Duck-Mayr
Benchmarking Rcpp code with RcppClock — Zach DeBruine
Simulation Smoother using RcppArmadillo — Tomasz Woźniak

Convex Hull of Polygon using Boost.Geometry

Related Articles