[WIP] setting up Inst/Include header

[dcooley]

I'm making a work-in-progress PR so you can see what I'm doing as and when I commit changes.

I'll let you know when I think it's ready

---

TODO

• [x] move header to /inst/include and compile
• [x] test LinkingTo from Rcpp code / packages
• [x] SEXP API
• [ ] expose earcut( Rcpp::List polygon )
• [x] test exposed function (for out-of-range errors, unsupported types, the usual.)

[mdsumner]

Oh sweet thanks so much. Hoping to level up and contribute more soonish!

[dcooley]

I've set up the header, and made Rcpp wrappers for the Point and Polygon objects.

I've written a proof-of-concept earcht_sfc() which (when it works) will return each coordinate pair for each triangle, in a list[[ matrix ]] structure. At the moment it doesn't work when you include holes.

x <- c(0, 0, 0.75, 1, 0.5, 0.8, 0.69)
y <- c(0, 1, 1, 0.8, 0.7, 0.6, 0)

(ind <- earcut(cbind(x, y)))
# [1] 2 1 7 7 6 5 5 4 3 2 7 5 5 3 2
plot_ears(cbind(x, y), ind)

mat <- cbind(x,y)
sfg <- sfheaders::sfg_polygon( mat, close = F )
(ind <- decido:::earcut_sfc( sfg )) + 1
# [1] 2 1 7 7 6 5 5 4 3 2 7 5 5 3 2
plot_ears(cbind(x, y), ind + 1)

[dcooley]

Here's a proof-of-concept on how it could be called by sfheaders to return sf objects

library(Rcpp)

cppFunction(
  depends = c("decido","sfheaders","geometries") ## geometries because I'm using dev version of sfheaders
  , includes = c(
    '#include "decido/decido.hpp"'
    , '#include "sfheaders/df/sfg.hpp"'
    , '#include "sfheaders/sfheaders.hpp"'
    )
  , code = '
    SEXP triangles( SEXP sfg ) {
      Polygons polyrings = Rcpp::as< Polygons >( sfg );
      std::vector< uint32_t > indices = mapbox::earcut< uint32_t >( polyrings );
      Rcpp::IntegerVector idx = Rcpp::wrap( indices );

      Rcpp::DataFrame df = sfheaders::df::sfg_to_df( sfg );
      Rcpp::NumericVector x = df["x"];
      Rcpp::NumericVector y = df["y"];
      Rcpp::NumericVector xx = x[ idx ];
      Rcpp::NumericVector yy = y[ idx ];
      int n_triangles = idx.length() / 3;
      Rcpp::IntegerVector triangle_idx = Rcpp::seq(1, n_triangles );
      Rcpp::IntegerVector triangle_ids = Rcpp::rep_each( triangle_idx, 3 );

      Rcpp::DataFrame df_tri = Rcpp::DataFrame::create(
        Rcpp::_["triangle_id"] = triangle_ids,
        Rcpp::_["x"] = xx,
        Rcpp::_["y"] = yy
      );

      Rcpp::StringVector geometry_cols {"x","y"};
      Rcpp::String polygon_id = "triangle_id";
      std::string xyzm = "XY";
      bool close = false;
      Rcpp::List sfg_tri = sfheaders::sfc::sfc_polygon(
        df_tri, geometry_cols, polygon_id, polygon_id, xyzm, close
      );
      return sfg_tri;
    }
  '
)

library(sf)

library(mapdeck)
set_token( secret::get_secret("mapbox"))

nc <- sf::st_read(system.file("./shape/nc.shp", package = "sf"))

sf_poly <- sfheaders::sf_cast( nc, "POLYGON" )
sfg <- sf_poly[1, ]$geometry[[1]]
sf_tri <- triangles( sfg )

sf <- sf::st_as_sf( sf_tri ) 
sf$id <- 1:nrow( sf )

mapdeck() %>%
  add_polygon(
    data = sf
    , fill_colour = "id"
    , stroke_width = 20
    , tooltip = "id"
    , legend = T
  )

[dcooley]

Note to me for future reference.

I tried making a Decido class inherting Earcut, which had Rcpp-specific functions, but it was slower

sf <- geojsonsf::geojson_sf("https://symbolixau.github.io/data/geojson/SA2_2016_VIC.json")
sf_poly <- sfheaders::sf_cast( sf, "POLYGON" )
sfc <- sf_poly$geometry

lst <- decido:::earcut_sfc( sfc )
lst2 <- decido:::earcut_sfc2( sfc )

all.equal(lst, lst2)
# [1] TRUE

microbenchmark::microbenchmark(
  wrap = { lst <- decido:::earcut_sfc( sfc ) },
  rcpp = { lst2 <- decido:::earcut_sfc2( sfc ) }
)

# Unit: milliseconds
# expr      min       lq     mean   median       uq      max neval
# wrap 203.1208 210.5301 225.9748 213.9351 218.7850 311.4982   100
# rcpp 285.7782 305.3584 337.2406 311.4472 385.1729 473.1197   100

(on branch inherited)

[dcooley]

@mdsumner I think with this last commit I've finished what I needed to do so is ready for your review.

Things to note

• I've made a decido::api:: in /inst/include/decido/decido.hpp with two overloaded earcut() functions
• I've moved the internals of your earcut_cpp() code into this API, and added a second earcut() function which works on a List[[ Matrix ]] structure
• I've not exposed this second function through R (yet). But we can if you think it's worth having it?
• I have not changed any of the R code.

Examples

sf <- geojsonsf::geojson_sf("https://symbolixau.github.io/data/geojson/SA2_2016_VIC.json")
sf_poly <- sfheaders::sf_cast( sf, "POLYGON" )

library(Rcpp)

## Single Polygon 
cppFunction(
  depends = "decido"
  , includes = '#include "decido/decido.hpp"'
  , code = '
    Rcpp::IntegerVector earcut( Rcpp::List polygon ) {
      return decido::api::earcut( polygon );
    }
  '
)

## Single Polygon
which( sapply( sf_poly$geometry, length ) > 1 )
r <- 17
idx <- earcut( sf_poly[r, ]$geometry[[1]] )
df <- sfheaders::sf_to_df( sf[r,  ])

ids <- rep(1:(length(idx)/3), each = 3)

df <- df[ idx, ]
df$id <- ids

sf_triangles <- sfheaders::sf_polygon(
  obj = df
  , polygon_id = "id"
  , x = "x"
  , y = "y"
  )

mapdeck() %>%
  add_polygon(
    data = sf_triangles
    , fill_colour = "id"
    , fill_opacity = 0.6
    , tooltip = "id"
  )

### Multiple Polygons
cppFunction(
  depends = "decido"
  , includes = '#include "decido/decido.hpp"'
  , code = '
    Rcpp::List earcuts( Rcpp::List polygons ) {
      R_xlen_t i;
      R_xlen_t n = polygons.length();
      Rcpp::List res( n );
      for( i = 0; i < n; ++i ) {
        Rcpp::List polygon = polygons[ i ];
        res[i] = decido::api::earcut( polygon );
      }
      return res;
    }
  '
)

lst <- earcuts( sf_poly$geometry )

str( lst )
# List of 616
# $ : int [1:270] 92 91 90 90 89 88 88 87 86 86 ...
# $ : int [1:669] 225 224 223 223 222 221 220 219 218 218 ...
# $ : int [1:999] 333 332 331 331 330 329 329 328 327 327 ...
# $ : int [1:858] 1 288 287 287 286 285 285 284 283 281 ...
# $ : int [1:1626] 544 543 542 542 541 540 539 538 537 537 ...
# $ : int [1:258] 1 88 87 87 86 85 85 84 83 83 ...

[mdsumner]

Well, great day's work (!) - I can't see any problems. Thanks! I definitely think it's worth exposing that new function. Do you want to suggest a NEWS item? I did have one drafted earlier today but you might as well write it. It's Saturday night and I'm unlikely to delve much but really excited with what you've done here. I'll merge now if you want, it's passing check and hasn't changed anything that I was doing before afaics.

[dcooley]

It's Saturday night and I'm unlikely to delve much

yeah me either now I'm two glasses of wine down. But I'll write the NEWs entry tomorrow.

[mdsumner]

yeah let's step away haha

[dcooley]

I've updated the NEWS and added some tests.

I think the task of "exposing to R" is better suited to the discussion in #9 , as I think it needs a philosophical discussion on what the inputs and outputs are going to be, and which other (if any) packages are required as dependencies.

Therefore, I'm inclined to say this PR is complete.

[mdsumner]

Brilliant, thanks so much! No excuses for me now I have to get into this headers stuff :)