add rcpp_points_to_edges_par fn for #103

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: dodgr
 Title: Distances on Directed Graphs
-Version: 0.2.14.081
+Version: 0.2.14.082
 Authors@R: c(
     person("Mark", "Padgham", , "[email protected]", role = c("aut", "cre")),
     person("Andreas", "Petutschnig", role = "aut"),

R/RcppExports.R

@@ -249,12 +249,22 @@ rcpp_unique_rownames <- function(xyfrom, xyto, precision = 10L) {
     .Call(`_dodgr_rcpp_unique_rownames`, xyfrom, xyto, precision)
 }
 
+#' Simple match of points to nearest vertices
+#' @noRd
+NULL
+
+#' Match points to nearest edge of graph at which perpendicular from point
+#' bisects edges. Uses psuedo-code from
+#' https://stackoverflow.com/a/6853926
+#' @noRd
+NULL
+
 #' rcpp_points_index_par
 #'
 #' Get index of nearest vertices to list of points
 #'
-#' @param graph Rcpp::DataFrame containing the graph
-#' @param pts Rcpp::DataFrame containing the routing points
+#' @param xy Rcpp::DataFrame containing the vertex coordinates of the graph
+#' @param pts Rcpp::DataFrame containing the points to be matched
 #'
 #' @return 0-indexed Rcpp::NumericVector index into graph of nearest points
 #'
@@ -263,6 +273,20 @@ rcpp_points_index_par <- function(xy, pts) {
     .Call(`_dodgr_rcpp_points_index_par`, xy, pts)
 }
 
+#' rcpp_points_to_edges_par
+#'
+#' Get index of nearest edges to list of points
+#'
+#' @param graph Rcpp::DataFrame containing the full edge-based graph
+#' @param pts Rcpp::DataFrame containing the points to be matched
+#'
+#' @return 0-indexed Rcpp::NumericVector index into graph of nearest points
+#'
+#' @noRd
+rcpp_points_to_edges_par <- function(graph, pts) {
+    .Call(`_dodgr_rcpp_points_to_edges_par`, graph, pts)
+}
+
 #' rcpp_get_sp_dists_par
 #'
 #' @noRd

codemeta.json

@@ -7,7 +7,7 @@
   "codeRepository": "https://github.com/ATFutures/dodgr",
   "issueTracker": "https://github.com/ATFutures/dodgr/issues",
   "license": "https://spdx.org/licenses/GPL-3.0",
-  "version": "0.2.14.081",
+  "version": "0.2.14.082",
   "programmingLanguage": {
     "@type": "ComputerLanguage",
     "name": "R",

src/RcppExports.cpp

@@ -178,6 +178,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// rcpp_points_to_edges_par
+Rcpp::IntegerVector rcpp_points_to_edges_par(const Rcpp::DataFrame& graph, Rcpp::DataFrame& pts);
+RcppExport SEXP _dodgr_rcpp_points_to_edges_par(SEXP graphSEXP, SEXP ptsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const Rcpp::DataFrame& >::type graph(graphSEXP);
+    Rcpp::traits::input_parameter< Rcpp::DataFrame& >::type pts(ptsSEXP);
+    rcpp_result_gen = Rcpp::wrap(rcpp_points_to_edges_par(graph, pts));
+    return rcpp_result_gen;
+END_RCPP
+}
 // rcpp_get_sp_dists_par
 Rcpp::NumericMatrix rcpp_get_sp_dists_par(const Rcpp::DataFrame graph, const Rcpp::DataFrame vert_map_in, Rcpp::IntegerVector fromi, Rcpp::IntegerVector toi_in, const std::string& heap_type, const bool is_spatial);
 RcppExport SEXP _dodgr_rcpp_get_sp_dists_par(SEXP graphSEXP, SEXP vert_map_inSEXP, SEXP fromiSEXP, SEXP toi_inSEXP, SEXP heap_typeSEXP, SEXP is_spatialSEXP) {
@@ -352,6 +364,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_dodgr_rcpp_get_component_vector", (DL_FUNC) &_dodgr_rcpp_get_component_vector, 1},
     {"_dodgr_rcpp_unique_rownames", (DL_FUNC) &_dodgr_rcpp_unique_rownames, 3},
     {"_dodgr_rcpp_points_index_par", (DL_FUNC) &_dodgr_rcpp_points_index_par, 2},
+    {"_dodgr_rcpp_points_to_edges_par", (DL_FUNC) &_dodgr_rcpp_points_to_edges_par, 2},
     {"_dodgr_rcpp_get_sp_dists_par", (DL_FUNC) &_dodgr_rcpp_get_sp_dists_par, 6},
     {"_dodgr_rcpp_get_sp_dists_paired_par", (DL_FUNC) &_dodgr_rcpp_get_sp_dists_paired_par, 6},
     {"_dodgr_rcpp_get_iso", (DL_FUNC) &_dodgr_rcpp_get_iso, 5},

src/match-points.cpp

@@ -1,5 +1,7 @@
 #include "match-points.h"
 
+//' Simple match of points to nearest vertices
+//' @noRd
 struct OnePointIndex : public RcppParallel::Worker
 {
     const RcppParallel::RVector <double> xy_x, xy_y, pt_x, pt_y;
@@ -42,6 +44,91 @@ struct OnePointIndex : public RcppParallel::Worker
 
 };
 
+//' Match points to nearest edge of graph at which perpendicular from point
+//' bisects edges. Uses psuedo-code from
+//' https://stackoverflow.com/a/6853926
+//' @noRd
+struct OneEdgeIndex : public RcppParallel::Worker
+{
+    const RcppParallel::RVector <double> pt_x, pt_y,
+          xfr, yfr, xto, yto;
+    const size_t nxy;
+    RcppParallel::RVector <int> index;
+
+    // constructor
+    OneEdgeIndex (
+            const RcppParallel::RVector <double> pt_x_in,
+            const RcppParallel::RVector <double> pt_y_in,
+            const RcppParallel::RVector <double> xfr_in,
+            const RcppParallel::RVector <double> yfr_in,
+            const RcppParallel::RVector <double> xto_in,
+            const RcppParallel::RVector <double> yto_in,
+            const size_t nxy_in,
+            Rcpp::IntegerVector index_in) :
+        pt_x (pt_x_in), pt_y (pt_y_in),
+        xfr (xfr_in), yfr (yfr_in), xto (xto_in), yto (yto_in),
+        nxy (nxy_in), index (index_in)
+    {
+    }
+
+    // Parallel function operator
+    void operator() (std::size_t begin, std::size_t end)
+    {
+        for (std::size_t i = begin; i < end; i++)
+        {
+            double dmin = INFINITE_DOUBLE;
+            long int jmin = INFINITE_INT;
+
+            for (size_t j = 0; j < nxy; j++)
+            {
+                const double x1 = xfr [j], y1 = yfr [j];
+                const double x2 = xto [j], y2 = yto [j];
+
+                const double A = pt_x [i] - x1;
+                const double B = pt_y [i] - y1;
+                const double C = x2 - x1;
+                const double D = y2 - y1;
+
+                const double dot = A * C + B * D;
+                const double len_sq = C * C + D * D;
+                double param = -1.0;
+                if (fabs (len_sq) < 1.0e-12)
+                {
+                    param = dot / len_sq;
+                }
+
+                double xx, yy;
+                if (param < 0.0)
+                {
+                    xx = x1;
+                    yy = y1;
+                } else if (param > 1.0)
+                {
+                    xx = x2;
+                    yy = y2;
+                } else
+                {
+                    xx = x1 + param * C;
+                    yy = y1 + param * D;
+                }
+
+                const double dx = pt_x [i] - xx;
+                const double dy = pt_y [i] - yy;
+
+                const double dij = sqrt (dx * dx + dy * dy);
+
+                if (dij < dmin)
+                {
+                    dmin = dij;
+                    jmin = static_cast <long int> (j);
+                }
+            }
+            index [i] = static_cast <int> (jmin);
+        }
+    }
+
+};
+
 //' rcpp_points_index_par
 //'
 //' Get index of nearest vertices to list of points
@@ -77,3 +164,44 @@ Rcpp::IntegerVector rcpp_points_index_par (const Rcpp::DataFrame &xy,
 
     return index;
 }
+
+//' rcpp_points_to_edges_par
+//'
+//' Get index of nearest edges to list of points
+//'
+//' @param graph Rcpp::DataFrame containing the full edge-based graph
+//' @param pts Rcpp::DataFrame containing the points to be matched
+//'
+//' @return 0-indexed Rcpp::NumericVector index into graph of nearest points
+//'
+//' @noRd
+// [[Rcpp::export]]
+Rcpp::IntegerVector rcpp_points_to_edges_par (const Rcpp::DataFrame &graph,
+        Rcpp::DataFrame &pts)
+{
+    Rcpp::NumericVector ptx = pts ["x"];
+    Rcpp::NumericVector pty = pts ["y"];
+
+    Rcpp::NumericVector xfr = graph ["xfr"];
+    Rcpp::NumericVector yfr = graph ["yfr"];
+    Rcpp::NumericVector xto = graph ["xto"];
+    Rcpp::NumericVector yto = graph ["yto"];
+
+    const size_t npts = static_cast <size_t> (pts.nrow ()),
+           nxy = static_cast <size_t> (graph.nrow ());
+
+    Rcpp::IntegerVector index (npts);
+
+    // Create parallel worker
+    OneEdgeIndex one_edge_indx (RcppParallel::RVector <double> (ptx),
+            RcppParallel::RVector <double> (pty),
+            RcppParallel::RVector <double> (xfr),
+            RcppParallel::RVector <double> (yfr),
+            RcppParallel::RVector <double> (xto),
+            RcppParallel::RVector <double> (yto),
+            nxy, index);
+
+    RcppParallel::parallelFor (0, npts, one_edge_indx);
+
+    return index;
+}

src/match-points.h

@@ -11,5 +11,5 @@ constexpr int INFINITE_INT =  std::numeric_limits<int>::max ();
 
 Rcpp::IntegerVector rcpp_points_index (const Rcpp::DataFrame &xy,
         Rcpp::DataFrame &pts);
-Rcpp::IntegerVector rcpp_points_index_par (const Rcpp::DataFrame &xy,
+Rcpp::IntegerVector rcpp_points_to_edges_par (const Rcpp::DataFrame &graph,
         Rcpp::DataFrame &pts);