include bin version of point-in-polygon

R/RcppExports.R

@@ -85,6 +85,10 @@ sample_points_grid <- function(polymatrix, n_sam_pol) {
     .Call(`_sits_sample_points_grid`, polymatrix, n_sam_pol)
 }
 
+sample_points_bin <- function(polymatrix, n_sam_pol) {
+    .Call(`_sits_sample_points_bin`, polymatrix, n_sam_pol)
+}
+
 C_max_sampling <- function(x, nrows, ncols, window_size) {
     .Call(`_sits_C_max_sampling`, x, nrows, ncols, window_size)
 }

R/api_segments.R

@@ -521,7 +521,7 @@
                 purrr::map_dfr(function(i){
                     polymatrix <- sf::st_coordinates(sf_pols[i,])
                     polymatrix <- polymatrix[, 1:2]
-                    points_mx <- sample_points_inclusion(polymatrix,
+                    points_mx <- sample_points_bin(polymatrix,
                                                 sf_pols[i,]$n_sam_pol)
                     colnames(points_mx) <- c("longitude", "latitude")
                     points_row <- tibble::as_tibble(points_mx)
@@ -532,3 +532,5 @@
         })
     return(samples)
 }
+
+# sf::st_sample(sf_object[i, ], size = n_sam_pol)

src/RcppExports.cpp

@@ -285,6 +285,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// sample_points_bin
+NumericMatrix sample_points_bin(const NumericMatrix& polymatrix, const int n_sam_pol);
+RcppExport SEXP _sits_sample_points_bin(SEXP polymatrixSEXP, SEXP n_sam_polSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const NumericMatrix& >::type polymatrix(polymatrixSEXP);
+    Rcpp::traits::input_parameter< const int >::type n_sam_pol(n_sam_polSEXP);
+    rcpp_result_gen = Rcpp::wrap(sample_points_bin(polymatrix, n_sam_pol));
+    return rcpp_result_gen;
+END_RCPP
+}
 // C_max_sampling
 DataFrame C_max_sampling(const NumericVector& x, int nrows, int ncols, int window_size);
 RcppExport SEXP _sits_C_max_sampling(SEXP xSEXP, SEXP nrowsSEXP, SEXP ncolsSEXP, SEXP window_sizeSEXP) {
@@ -456,6 +468,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_sits_sample_points_inclusion", (DL_FUNC) &_sits_sample_points_inclusion, 2},
     {"_sits_sample_points_crossings", (DL_FUNC) &_sits_sample_points_crossings, 2},
     {"_sits_sample_points_grid", (DL_FUNC) &_sits_sample_points_grid, 2},
+    {"_sits_sample_points_bin", (DL_FUNC) &_sits_sample_points_bin, 2},
     {"_sits_C_max_sampling", (DL_FUNC) &_sits_C_max_sampling, 4},
     {"_sits_bayes_smoother", (DL_FUNC) &_sits_bayes_smoother, 6},
     {"_sits_bayes_var", (DL_FUNC) &_sits_bayes_var, 5},

src/sample_points.cpp

@@ -1,6 +1,7 @@
 #include <Rcpp.h>
 #include <iostream>
 #include <vector>
+#include <list>
 #include <cstdlib>
 #include <ctime>
 using namespace Rcpp;
@@ -119,6 +120,83 @@ std::vector<Point> polygon_grid_points (int n_intervals,
     }
     return points;
 }
+// ======= Trapezoid (bin) algorithm =======================================
+
+// Split polygons along set of y bins and sorts the edge fragments.  Testing
+// is done against these fragments.
+
+struct Edge {
+    Point p1, p2;
+};
+
+class BinnedPolygon {
+public:
+    BinnedPolygon(const std::vector<Point>& vertices, int binSize)
+        : binSize(binSize) {
+        // Find the bounding box of the polygon
+        for (const auto& v : vertices) {
+            minX = std::min(minX, v.x);
+            minY = std::min(minY, v.y);
+            maxX = std::max(maxX, v.x);
+            maxY = std::max(maxY, v.y);
+        }
+
+        // Create the bins
+        int binsX = (maxX - minX) / binSize + 1;
+        int binsY = (maxY - minY) / binSize + 1;
+        bins.resize(binsX, std::vector<std::list<Edge>>(binsY));
+
+        // Insert the edges into the bins
+        for (size_t i = 0; i < vertices.size(); ++i) {
+            size_t j = (i + 1) % vertices.size();
+            insertEdge({vertices[i], vertices[j]});
+        }
+    }
+
+    bool isInside(const Point& point) const {
+        int binX = (point.x - minX) / binSize;
+        int binY = (point.y - minY) / binSize;
+
+        if (binX < 0 || binX >= bins.size() || binY < 0 || binY >= bins[0].size())
+            return false;
+
+        // Ray casting algorithm
+        int intersections = 0;
+        for (const auto& edge : bins[binX][binY]) {
+            if ((edge.p1.y > point.y) != (edge.p2.y > point.y) &&
+                point.x < (edge.p2.x - edge.p1.x) * (point.y - edge.p1.y) /
+                    (edge.p2.y - edge.p1.y) + edge.p1.x) {
+                ++intersections;
+            }
+        }
+        return (intersections % 2) == 1;
+    }
+
+private:
+    void insertEdge(const Edge& edge) {
+        // Find the bins that this edge intersects with
+        int minXBin = std::min(edge.p1.x, edge.p2.x);
+        int maxXBin = std::max(edge.p1.x, edge.p2.x);
+        int minYBin = std::min(edge.p1.y, edge.p2.y);
+        int maxYBin = std::max(edge.p1.y, edge.p2.y);
+
+        for (int x = minXBin; x <= maxXBin; x += binSize) {
+            for (int y = minYBin; y <= maxYBin; y += binSize) {
+                int binX = (x - minX) / binSize;
+                int binY = (y - minY) / binSize;
+                bins[binX][binY].push_back(edge);
+            }
+        }
+    }
+
+    double minX = std::numeric_limits<double>::infinity();
+    double minY = std::numeric_limits<double>::infinity();
+    double maxX = -std::numeric_limits<double>::infinity();
+    double maxY = -std::numeric_limits<double>::infinity();
+    int binSize;
+    std::vector<std::vector<std::list<Edge>>> bins;
+};
+
 // [[Rcpp::export]]
 NumericMatrix sample_points_inclusion(const NumericMatrix& polymatrix,
                             const int n_sam_pol){
@@ -221,3 +299,43 @@ NumericMatrix sample_points_grid(const NumericMatrix& polymatrix,
     }
     return points_mx;
 }
+// [[Rcpp::export]]
+NumericMatrix sample_points_bin(const NumericMatrix& polymatrix,
+                                 const int n_sam_pol){
+
+    int n_bins = 20;
+    NumericMatrix points(n_sam_pol, 2); // output
+    // internal data structures
+    int n_vertices = polymatrix.nrow();
+    std::vector<Point> polygon(n_vertices);
+    // build polygon from NumericMatrix
+    for (size_t i = 0; i < n_vertices; ++i){
+        polygon[i].x = polymatrix(i, 0);
+        polygon[i].y = polymatrix(i, 1);
+    }
+    BinnedPolygon binnedPolygon(polygon, n_bins);
+
+    // find max and min
+    double minX = polygon[0].x, maxX = polygon[0].x;
+    double minY = polygon[0].y, maxY = polygon[0].y;
+    for (const auto& vertex : polygon) {
+        minX = std::min(minX, vertex.x);
+        maxX = std::max(maxX, vertex.x);
+        minY = std::min(minY, vertex.y);
+        maxY = std::max(maxY, vertex.y);
+    }
+
+    // Generate points inside the polygon
+    int i = 0;
+    while (i < n_sam_pol) {
+        Point p;
+        p.x = minX + (maxX - minX) * (rand() / static_cast<double>(RAND_MAX));
+        p.y = minY + (maxY - minY) * (rand() / static_cast<double>(RAND_MAX));
+        if (binnedPolygon.isInside(p)) {
+            points(i, 0) = p.x;
+            points(i, 1) = p.y;
+            i++;
+        }
+    }
+    return points;
+}