CODEWITHSUNDEEP
c++rboostrcppboost-geometry
politicalEconomist
politicalEconomist

Reputation: 1041

Converting Rcpp NumericMatrix for use with Boost Geometry

I am finding that I am lost without the nice <as> and <wrap> commands that Rcpp and their related packages provide for conversion between different object types.

I have a matrix of points for which the rows represent points in a two dimensional cartesian space:

 pointsMatrix <- matrix(runif(100,-1,1),50,50)

I then want to use the convex_hull algorithm from boost geometry to find the convex hull of the points.

However, I am not sure how to convert the NumericMatrix into one of the data types that convex_hull understands. Further I am not sure how to convert the output from Boost Geometry back into something Rcpp can hand back to R.

 #include <Rcpp.h>
 #include <boost/geometry.hpp>
 #include <boost/geometry/geometries/polygon.hpp>
 using namespace Rcpp;

 BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian)

 // [[Rcpp::export]]
 NumericMatrix convexHullRcpp(NumericMatrix pointsMatrix){

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

 // Some sort of conversion of pointsMatrix here to pointsMatrixBG//

     polygon hull;
     boost::geometry::convex_hull(pointsMatrixBG, hull);

     //Now to convert hull into something that Rcpp can hand back to R.//

  return hullToR;
 }

It looks like boost.tuple may be the best choice

Upvotes: 1

Views: 822

Answers (2)

Sameer
Sameer

Reputation: 1807

Here is a little test.cpp file

#include <Rcpp.h>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/point_xy.hpp>

using namespace Rcpp;

typedef boost::geometry::model::d2::point_xy<double, boost::geometry::cs::cartesian> Point;
typedef boost::geometry::model::polygon<Point, true, true> Polygon; 

namespace Rcpp {
  template <> Polygon as(SEXP pointsMatrixSEXP) {
    NumericMatrix pointsMatrix(pointsMatrixSEXP);
    Polygon pointsMatrixBG;
    for (int i = 0; i < pointsMatrix.nrow(); ++i) {
      double x = pointsMatrix(i,0);
      double y = pointsMatrix(i,1);
      Point p(x,y);
      pointsMatrixBG.outer().push_back(p); 
    }
    return (pointsMatrixBG);
  } 

  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.x();
      rmat(i,1) = p.y();
    }
    return Rcpp::wrap(rmat);
  }
}

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

  // Conversion of pointsMatrix here to pointsMatrixBG
  Polygon pointsMatrixBG = as<Polygon>(pointsMatrixSEXP);

  Polygon hull;
  boost::geometry::convex_hull(pointsMatrixBG, hull);

  //Now to convert hull into something that Rcpp can hand back to R.//
  return wrap(hull);
}

Then you can test it out in R.

library(Rcpp)
sourceCpp("test.cpp")
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)
> 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

Upvotes: 5

Dirk is no longer here
Dirk is no longer here

Reputation: 368301

In general, once you go to types not known by R, you need to build your own custom converter functions as<>() and wrap().

As noted in the comment, there is an entire vignette devoted to it. There are also examples in packages as well as eg this article in the Rcpp Gallery about custom as<>() and wrap(). I would start from the Boost Geometry examples in C++ to understand how their datastructures are filled, and then fill them from my C++. There is no shortcut. The No Free Lunch theorem still holds.

Upvotes: 0

Related Questions