From c90491d4fd69daed9499a00d8eb0c8fec877edb8 Mon Sep 17 00:00:00 2001 From: mpadge Date: Wed, 10 Aug 2022 16:05:06 +0200 Subject: [PATCH] add rcpp_points_to_edges_par fn for #103 --- DESCRIPTION | 2 +- R/RcppExports.R | 28 +++++++++- codemeta.json | 2 +- src/RcppExports.cpp | 13 +++++ src/match-points.cpp | 128 +++++++++++++++++++++++++++++++++++++++++++ src/match-points.h | 2 +- 6 files changed, 170 insertions(+), 5 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5f7090900..42d89d21c 100644 --- a/DESCRIPTION +++ b/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", , "mark.padgham@email.com", role = c("aut", "cre")), person("Andreas", "Petutschnig", role = "aut"), diff --git a/R/RcppExports.R b/R/RcppExports.R index 83690af21..24c5cdc1f 100644 --- a/R/RcppExports.R +++ b/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 diff --git a/codemeta.json b/codemeta.json index 5198b94ea..617df0586 100644 --- a/codemeta.json +++ b/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", diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 7946ddd5a..e8041d3d2 100644 --- a/src/RcppExports.cpp +++ b/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}, diff --git a/src/match-points.cpp b/src/match-points.cpp index 0dd7f6e03..ad32d2537 100644 --- a/src/match-points.cpp +++ b/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 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 pt_x, pt_y, + xfr, yfr, xto, yto; + const size_t nxy; + RcppParallel::RVector index; + + // constructor + OneEdgeIndex ( + const RcppParallel::RVector pt_x_in, + const RcppParallel::RVector pt_y_in, + const RcppParallel::RVector xfr_in, + const RcppParallel::RVector yfr_in, + const RcppParallel::RVector xto_in, + const RcppParallel::RVector 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 (j); + } + } + index [i] = static_cast (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 (pts.nrow ()), + nxy = static_cast (graph.nrow ()); + + Rcpp::IntegerVector index (npts); + + // Create parallel worker + OneEdgeIndex one_edge_indx (RcppParallel::RVector (ptx), + RcppParallel::RVector (pty), + RcppParallel::RVector (xfr), + RcppParallel::RVector (yfr), + RcppParallel::RVector (xto), + RcppParallel::RVector (yto), + nxy, index); + + RcppParallel::parallelFor (0, npts, one_edge_indx); + + return index; +} diff --git a/src/match-points.h b/src/match-points.h index 5bc76794c..95c66206e 100644 --- a/src/match-points.h +++ b/src/match-points.h @@ -11,5 +11,5 @@ constexpr int INFINITE_INT = std::numeric_limits::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);