Remove dependency to boost

NEWS.md

@@ -7,6 +7,7 @@ If you are viewing this file on CRAN, please check [the latest news on GitHub](h
 - Fix: #757
 - Fix: algorithm `random_per_voxel` that was not working
 - Fix: `readLAScatalog` if `sp` is missing
+- Internal: remove dependency to `boost` for `point_in_polygon` #763
 
 ## lidR v4.1.1 (Release date: 2024-02-03)
 

R/RcppExports.R

@@ -49,8 +49,8 @@ C_lowest <- function(las, layout) {
     .Call(`_lidR_C_lowest`, las, layout)
 }
 
-C_in_polygon <- function(las, wkts, by_poly) {
-    .Call(`_lidR_C_in_polygon`, las, wkts, by_poly)
+C_in_polygon <- function(las, polygons, by_poly) {
+    .Call(`_lidR_C_in_polygon`, las, polygons, by_poly)
 }
 
 C_lasdetectshape <- function(las, method, th, k, filter, ncpu) {

R/utils_geometry.R

@@ -8,7 +8,29 @@ point_in_polygons <- function(las, sfc, by_poly = FALSE)
   id  <- as.integer(id)
   sfc2 <- sfc[id]
 
-  wkt <- sf::st_as_text(sfc2, digit = 10)
-  res <- C_in_polygon(las, wkt, by_poly)
+  rings = extract_polygons(sfc2)
+
+  res <- C_in_polygon(las, rings, by_poly)
   return(res)
 }
+
+extract_polygons = function(sfc)
+{
+  # Fix #763. Complex code to avoid using boost
+  coordinates = lapply(sfc, sf::st_coordinates)
+  rings = lapply(coordinates, extract_rings)
+}
+
+extract_rings = function(coordinates)
+{
+  if (ncol(coordinates) == 5L)
+    features = lapply( split( coordinates[,1:5], coordinates[,5] ), matrix, ncol = 5)
+  else
+    features = list(coordinates)
+
+  polygons = lapply(features,  function (x) lapply( split( x[,1:3], x[,4] ), matrix, ncol = 3))
+  rings = lapply(polygons,  function (x) lapply(x, function(x) lapply( split( x[,1:3], x[,3] ), matrix, ncol = 3)))
+  rings = lapply(rings, function(x) Reduce(c, x))
+
+  return(rings[[1]])
+}

src/LAS.cpp

@@ -3,10 +3,24 @@
 #include "myomp.h"
 #include "SpatialIndex.h"
 #include <limits>
-#include <boost/geometry.hpp>
 
 using namespace lidR;
 
+bool pnpoly(NumericMatrix polygon, double x, double y)
+{
+  int nvert = polygon.nrow();
+  bool c = false;
+  int i, j;
+  for (i = 0, j = nvert - 1; i < nvert; j = i++)
+  {
+    if (((polygon(i, 1) > y) != (polygon(j, 1) > y)) &&
+        (x < (polygon(j, 0) - polygon(i, 0)) * (y - polygon(i, 1)) / (polygon(j, 1) - polygon(i, 1)) + polygon(i, 0))) {
+      c = !c;
+    }
+  }
+  return c;
+}
+
 LAS::LAS(S4 las, int ncpu)
 {
   Rcpp::List index = las.slot("index");
@@ -586,18 +600,13 @@ void LAS::filter_with_grid(List layout, bool max)
   return;
 }
 
-SEXP LAS::find_polygon_ids(CharacterVector wkts, bool by_poly)
+SEXP LAS::find_polygon_ids(Rcpp::List polygons, bool by_poly)
 {
-  typedef boost::geometry::model::point<double, 2, boost::geometry::cs::cartesian> Point;
-  typedef boost::geometry::model::polygon<Point> Polygon;
-  typedef boost::geometry::model::multi_polygon<Polygon> MultiPolygon;
-  typedef boost::geometry::model::box<Point> Bbox;
-
   std::vector<std::vector<int>> res;
   std::vector<int> poly_id;
   if (by_poly)
   {
-    res.resize(wkts.size());
+    res.resize(polygons.size());
   }
   else
   {
@@ -607,76 +616,77 @@ SEXP LAS::find_polygon_ids(CharacterVector wkts, bool by_poly)
 
   SpatialIndex tree(las);
 
-  for (unsigned int j = 0 ; j < wkts.size() ; j++)
+  for (unsigned int i = 0 ; i < polygons.size() ; i++)
   {
     bool is_in_polygon = false;
-    std::string wkt = as<std::string>(wkts[j]);
-
-    if (wkt.find("MULTIPOLYGON") != std::string::npos)
-    {
-      Point p;
-      Bbox bbox;
-      MultiPolygon polygons;
 
-      boost::geometry::read_wkt(wkt, polygons);
-      boost::geometry::envelope(polygons, bbox);
+    // This list can be made of several rings (MULTIPOLYGON) and interior rings
+    Rcpp::List rings = polygons[i];
 
-      double min_x = bbox.min_corner().get<0>();
-      double min_y = bbox.min_corner().get<1>();
-      double max_x = bbox.max_corner().get<0>();
-      double max_y = bbox.max_corner().get<1>();
+    // Find the bbox of the polygons
+    double min_x = INFINITY;
+    double min_y = INFINITY;
+    double max_x = -INFINITY;
+    double max_y = -INFINITY;
+    for (int j = 0 ; j < rings.size() ; j++)
+    {
+      Rcpp::NumericMatrix ring = rings[j];
+      Rcpp::NumericVector x = ring(_, 0);
+      Rcpp::NumericVector y = ring(_, 1);
+
+      double maxx = max(x);
+      double maxy = max(y);
+      double minx = min(x);
+      double miny = min(y);
+
+      if (min_x > minx) min_x = minx;
+      if (min_y > miny) min_y = miny;
+      if (max_x < maxx) max_x = maxx;
+      if (max_y < maxy) max_y = maxy;
+    }
 
-      lidR::Rectangle rect(min_x, max_x, min_y, max_y);
-      std::vector<PointXYZ> pts;
-      tree.lookup(rect, pts);
+    // Spatial query of the point that are in the bounding box of the polygon
+    lidR::Rectangle rect(min_x, max_x, min_y, max_y);
+    std::vector<PointXYZ> pts;
+    tree.lookup(rect, pts);
 
-      for(unsigned int i = 0 ; i < pts.size() ; i++)
-      {
-        p.set<0>(pts[i].x);
-        p.set<1>(pts[i].y);
-        if (boost::geometry::covered_by(p, polygons))
-        {
-          if (by_poly)
-            res[j].push_back(pts[i].id+1);
-          else
-            poly_id[pts[i].id] = j+1;
-        }
-      }
-    }
-    else if (wkt.find("POLYGON") != std::string::npos)
+    // For each point we check if it is in the potential multi part polygon
+    for (unsigned int k = 0 ; k < pts.size() ; k++)
     {
-      Point p;
-      Bbox bbox;
-      Polygon polygon;
+      bool inpoly = false;
 
-      boost::geometry::read_wkt(wkt, polygon);
-      boost::geometry::envelope(polygon, bbox);
+      // Loop through sub polygons (ring)
+      for (int j = 0 ; j < rings.size() ; j++)
+      {
+        Rcpp::NumericMatrix ring = rings[j];
 
-      double min_x = bbox.min_corner().get<0>();
-      double min_y = bbox.min_corner().get<1>();
-      double max_x = bbox.max_corner().get<0>();
-      double max_y = bbox.max_corner().get<1>();
+        // We need to know if the ring is an exterior/interior ring (hole)
+        bool exterior_ring = ring(0,2) == 1;
 
-      lidR::Rectangle rect(min_x, max_x, min_y, max_y);
-      std::vector<PointXYZ> pts;
-      tree.lookup(rect, pts);
+        bool b = pnpoly(ring, pts[k].x, pts[k].y);
 
-      for(unsigned int i = 0 ; i < pts.size() ; i++)
-      {
-        p.set<0>(pts[i].x);
-        p.set<1>(pts[i].y);
-        if (boost::geometry::covered_by(p, polygon))
+        if (b)
         {
-          if (by_poly)
-            res[j].push_back(pts[i].id+1);
+          if (exterior_ring)
+          {
+            inpoly = true;
+          }
           else
-            poly_id[pts[i].id] = j+1;
+          {
+            inpoly = false;
+            break;
+          }
         }
       }
-    }
-    else
-      Rcpp::stop("Unexpected error in point_in_polygon: WKT is not a POLYGON or MULTIPOLYGON"); // # nocov
 
+      if (inpoly)
+      {
+        if (by_poly)
+          res[i].push_back(pts[k].id+1);
+        else
+          poly_id[pts[k].id] = i+1;
+      }
+    }
   }
 
   if (by_poly)

src/LAS.h

@@ -41,7 +41,7 @@ class LAS
 
     NumericVector rasterize(List layout, double subcircle, int method);
     NumericVector compute_range(DataFrame flightlines);
-    SEXP find_polygon_ids(CharacterVector wkts, bool by_poly);
+    SEXP find_polygon_ids(Rcpp::List polygons, bool by_poly);
     IntegerVector segment_snags(NumericVector neigh_radii, double low_int_thrsh, double uppr_int_thrsh, int pt_den_req, NumericMatrix BBPRthrsh_mat);
     IntegerVector segment_trees(double dt1, double dt2, double Zu, double R, double th_tree, double radius);
     List point_metrics(unsigned int k, double r, DataFrame data, int nalloc, SEXP call, SEXP env);

src/RcppExports.cpp

@@ -172,14 +172,14 @@ BEGIN_RCPP
 END_RCPP
 }
 // C_in_polygon
-SEXP C_in_polygon(S4 las, CharacterVector wkts, bool by_poly);
-RcppExport SEXP _lidR_C_in_polygon(SEXP lasSEXP, SEXP wktsSEXP, SEXP by_polySEXP) {
+SEXP C_in_polygon(S4 las, Rcpp::List polygons, bool by_poly);
+RcppExport SEXP _lidR_C_in_polygon(SEXP lasSEXP, SEXP polygonsSEXP, SEXP by_polySEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::traits::input_parameter< S4 >::type las(lasSEXP);
-    Rcpp::traits::input_parameter< CharacterVector >::type wkts(wktsSEXP);
+    Rcpp::traits::input_parameter< Rcpp::List >::type polygons(polygonsSEXP);
     Rcpp::traits::input_parameter< bool >::type by_poly(by_polySEXP);
-    rcpp_result_gen = Rcpp::wrap(C_in_polygon(las, wkts, by_poly));
+    rcpp_result_gen = Rcpp::wrap(C_in_polygon(las, polygons, by_poly));
     return rcpp_result_gen;
 END_RCPP
 }

src/RcppFunction.cpp

@@ -66,10 +66,10 @@ LogicalVector C_lowest(S4 las, List layout)
 }
 
 // [[Rcpp::export(rng = false)]]
-SEXP C_in_polygon(S4 las, CharacterVector wkts, bool by_poly)
+SEXP C_in_polygon(S4 las, Rcpp::List polygons, bool by_poly)
 {
   LAS pt(las);
-  return pt.find_polygon_ids(wkts, by_poly);
+  return pt.find_polygon_ids(polygons, by_poly);
 }
 
 // [[Rcpp::export(rng = false)]]