From 182db1b901ce1d081dff3046c00ce99231991de7 Mon Sep 17 00:00:00 2001 From: mpadge Date: Fri, 2 Sep 2022 11:19:27 +0200 Subject: [PATCH] calc geodesic dists to geometric edge intersections for #103 --- DESCRIPTION | 2 +- R/match-points.R | 47 ++++++++++++++++++++++++++++++++++++++--- codemeta.json | 2 +- src/match-points.cpp | 50 +++++++++++++++++++------------------------- 4 files changed, 67 insertions(+), 34 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 971427016..a19bd56a3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: dodgr Title: Distances on Directed Graphs -Version: 0.2.15.009 +Version: 0.2.15.010 Authors@R: c( person("Mark", "Padgham", , "mark.padgham@email.com", role = c("aut", "cre")), person("Andreas", "Petutschnig", role = "aut"), diff --git a/R/match-points.R b/R/match-points.R index a6fccce1f..b7568124c 100644 --- a/R/match-points.R +++ b/R/match-points.R @@ -155,10 +155,14 @@ match_pts_to_graph <- function (graph, xy, connected = FALSE) { # rcpp_points_index is 0-indexed, so ... graph_index <- as.integer (res [index]) + 1L - # edge_dist not yet used here; see #103 - # edge_dist <- res [index + nrow (xy)] - return (graph_index) + # ret <- data.frame ( + # index = graph_index, + # d_signed = signed_intersection_dists (graph, xy, res) + # ) + ret <- graph_index + + return (ret) } #' match_points_to_graph @@ -172,6 +176,43 @@ match_points_to_graph <- function (graph, xy, connected = FALSE) { match_pts_to_graph (graph, xy, connected = connected) } +#' Get geodesic distances to intersection points for match_pts_to_graph. +#' +#' @param res Output of 'rcpp_points_index' function +#' @noRd +signed_intersection_dists <- function (graph, xy, res) { + + n <- nrow (xy) + index <- seq (n) + + # rcpp_points_index is 0-indexed, so ... + graph_index <- as.integer (res [index]) + 1L + + # coordinates not yet used here; see #103 + xy_intersect <- data.frame ( + x = res [index + nrow (xy)], + y = res [index + 2L * nrow (xy)] + ) + + d <- geodist::geodist ( + xy, + xy_intersect, + paired = TRUE, + measure = "geodesic" + ) + + # Then coordinates of graph edges for sign calculation + xy_cols <- c ("xfr", "yfr", "xto", "yto") + gxy <- graph [graph_index, xy_cols] + + which_side <- (gxy$xto - gxy$xfr) * (xy_intersect$y - gxy$yfr) - + (gxy$yto - gxy$yfr) * (xy_intersect$x - gxy$xfr) + which_side [which_side < 0.0] <- -1 + which_side [which_side > 0.0] <- 1 + + return (d * which_side) +} + #' Insert new nodes into a graph, breaking edges at point of nearest #' intersection. #' diff --git a/codemeta.json b/codemeta.json index b9552dc1c..a599fba2b 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.15.009", + "version": "0.2.15.010", "programmingLanguage": { "@type": "ComputerLanguage", "name": "R", diff --git a/src/match-points.cpp b/src/match-points.cpp index 6bf63b94c..72609b04c 100644 --- a/src/match-points.cpp +++ b/src/match-points.cpp @@ -93,6 +93,8 @@ struct OneEdgeIndex : public RcppParallel::Worker for (std::size_t i = begin; i < end; i++) { double dmin = INFINITE_DOUBLE; + double x_intersect = INFINITE_DOUBLE; + double y_intersect = INFINITE_DOUBLE; long int jmin = INFINITE_INT; int which_side = INFINITE_INT; @@ -101,36 +103,19 @@ struct OneEdgeIndex : public RcppParallel::Worker 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 px = x2 - x1; + const double py = 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; - } + const double norm = px * px + py * py; - 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; - } + double u = ((pt_x [i] - x1) * px + (pt_y [i] - y1) * py) / norm; + u = (u > 1) ? 1 : ((u < 0) ? 0 : u); + + const double xx = x1 + u * px; + const double yy = y1 + u * py; - const double dx = pt_x [i] - xx; - const double dy = pt_y [i] - yy; + const double dx = xx - pt_x [i]; + const double dy = yy - pt_y [i]; const double dij = sqrt (dx * dx + dy * dy); @@ -140,10 +125,13 @@ struct OneEdgeIndex : public RcppParallel::Worker jmin = static_cast (j); which_side = which_side_of_line (xfr [j], yfr [j], xto [j], yto [j], pt_x [i], pt_y [i]); + x_intersect = xx; + y_intersect = yy; } } index [i] = static_cast (jmin); - index [i + npts] = which_side * dmin; + index [i + npts] = x_intersect; + index [i + 2L * npts] = y_intersect; } } @@ -210,7 +198,11 @@ Rcpp::NumericVector rcpp_points_to_edges_par (const Rcpp::DataFrame &graph, const size_t nxy = static_cast (graph.nrow ()), npts = static_cast (pts.nrow ()); - Rcpp::NumericVector index (npts * 2L); + // index holds three vectors: + // 1. the integer index + // 2. the x coordinates of the intersection points + // 3. the y coordinates of the intersection points + Rcpp::NumericVector index (npts * 3L); // Create parallel worker OneEdgeIndex one_edge_indx (RcppParallel::RVector (ptx),