Convex Hull of Polygon using Boost.Geometry

February 23, 2014
By

(This article was first published on Rcpp Gallery, and kindly contributed to R-bloggers)

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

To leave a comment for the author, please follow the link and comment on his blog: Rcpp Gallery.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.