calc geodesic dists to geometric edge intersections for #103

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", , "[email protected]", role = c("aut", "cre")),
     person("Andreas", "Petutschnig", role = "aut"),

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.
 #'

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",

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 <long int> (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 <double> (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 <size_t> (graph.nrow ()),
         npts = static_cast <size_t> (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 <double> (ptx),