Distances on Directed Graphs in R
mpadge.tngl.sh
1#include "match-points.h"
2
3//' Determine which side of intersecting line a point lies on.
4//'
5//' @param (ax, ay) Coordinates of one end of line
6//' @param (bx, by) Coordinates of other end of line
7//' @param (x, y) Coordinates of point
8//' @return 0 if point on line, -1 if to the left; +1 if to the right.
9//'
10//' @noRd
11int which_side_of_line (const double ax, const double ay,
12 const double bx, const double by, const double x, const double y)
13{
14 double val = (bx - ax) * (y - ay) - (by - ay) * (x - ax);
15 return (val > 0) ? 1 : ((val < 0) ? -1 : 0);
16}
17
18//' Simple match of points to nearest vertices
19//' @noRd
20struct OnePointIndex : public RcppParallel::Worker
21{
22 const RcppParallel::RVector <double> xy_x, xy_y, pt_x, pt_y;
23 const size_t nxy;
24 RcppParallel::RVector <int> index;
25
26 // constructor
27 OnePointIndex (
28 const RcppParallel::RVector <double> xy_x_in,
29 const RcppParallel::RVector <double> xy_y_in,
30 const RcppParallel::RVector <double> pt_x_in,
31 const RcppParallel::RVector <double> pt_y_in,
32 const size_t nxy_in,
33 Rcpp::IntegerVector index_in) :
34 xy_x (xy_x_in), xy_y (xy_y_in), pt_x (pt_x_in), pt_y (pt_y_in),
35 nxy (nxy_in), index (index_in)
36 {
37 }
38
39 // Parallel function operator
40 void operator() (std::size_t begin, std::size_t end)
41 {
42 for (std::size_t i = begin; i < end; i++)
43 {
44 double dmin = INFINITE_DOUBLE;
45 long int jmin = INFINITE_INT;
46 for (size_t j = 0; j < nxy; j++)
47 {
48 double dij = (xy_x [j] - pt_x [i]) * (xy_x [j] - pt_x [i]) +
49 (xy_y [j] - pt_y [i]) * (xy_y [j] - pt_y [i]);
50 if (dij < dmin)
51 {
52 dmin = dij;
53 jmin = static_cast <long int> (j);
54 }
55 }
56 index [i] = static_cast <int> (jmin);
57 }
58 }
59
60};
61
62//' Match points to nearest edge of graph at which perpendicular from point
63//' bisects edges. Uses psuedo-code from
64//' https://stackoverflow.com/a/6853926
65//' @noRd
66struct OneEdgeIndex : public RcppParallel::Worker
67{
68 const RcppParallel::RVector <double> pt_x, pt_y,
69 xfr, yfr, xto, yto;
70 const size_t nxy, npts;
71 RcppParallel::RVector <double> index;
72
73 // constructor
74 OneEdgeIndex (
75 const RcppParallel::RVector <double> pt_x_in,
76 const RcppParallel::RVector <double> pt_y_in,
77 const RcppParallel::RVector <double> xfr_in,
78 const RcppParallel::RVector <double> yfr_in,
79 const RcppParallel::RVector <double> xto_in,
80 const RcppParallel::RVector <double> yto_in,
81 const size_t nxy_in,
82 const size_t npts_in,
83 Rcpp::NumericVector index_in) :
84 pt_x (pt_x_in), pt_y (pt_y_in),
85 xfr (xfr_in), yfr (yfr_in), xto (xto_in), yto (yto_in),
86 nxy (nxy_in), npts (npts_in), index (index_in)
87 {
88 }
89
90 // Parallel function operator
91 void operator() (std::size_t begin, std::size_t end)
92 {
93 for (std::size_t i = begin; i < end; i++)
94 {
95 double dmin = INFINITE_DOUBLE;
96 double x_intersect = INFINITE_DOUBLE;
97 double y_intersect = INFINITE_DOUBLE;
98 long int jmin = INFINITE_INT;
99 int which_side = INFINITE_INT;
100
101 for (size_t j = 0; j < nxy; j++)
102 {
103 const double x1 = xfr [j], y1 = yfr [j];
104 const double x2 = xto [j], y2 = yto [j];
105
106 const double px = x2 - x1;
107 const double py = y2 - y1;
108
109 const double norm = px * px + py * py;
110
111 double u = ((pt_x [i] - x1) * px + (pt_y [i] - y1) * py) / norm;
112 u = (u > 1) ? 1 : ((u < 0) ? 0 : u);
113
114 const double xx = x1 + u * px;
115 const double yy = y1 + u * py;
116
117 const double dx = xx - pt_x [i];
118 const double dy = yy - pt_y [i];
119
120 const double dij = sqrt (dx * dx + dy * dy);
121
122 if (dij < dmin)
123 {
124 dmin = dij;
125 jmin = static_cast <long int> (j);
126 which_side = which_side_of_line (xfr [j], yfr [j],
127 xto [j], yto [j], pt_x [i], pt_y [i]);
128 x_intersect = xx;
129 y_intersect = yy;
130 }
131 }
132 index [i] = static_cast <double> (jmin);
133 index [i + npts] = x_intersect;
134 index [i + 2L * npts] = y_intersect;
135 }
136 }
137
138};
139
140//' rcpp_points_index_par
141//'
142//' Get index of nearest vertices to list of points
143//'
144//' @param xy Rcpp::DataFrame containing the vertex coordinates of the graph
145//' @param pts Rcpp::DataFrame containing the points to be matched
146//'
147//' @return 0-indexed Rcpp::NumericVector index into graph of nearest points
148//'
149//' @noRd
150// [[Rcpp::export]]
151Rcpp::IntegerVector rcpp_points_index_par (const Rcpp::DataFrame &xy,
152 Rcpp::DataFrame &pts)
153{
154 Rcpp::NumericVector ptx = pts ["x"];
155 Rcpp::NumericVector pty = pts ["y"];
156
157 Rcpp::NumericVector vtx = xy ["x"];
158 Rcpp::NumericVector vty = xy ["y"];
159
160 size_t npts = static_cast <size_t> (pts.nrow ()),
161 nxy = static_cast <size_t> (xy.nrow ());
162
163 //Rcpp::IntegerVector index (n, Rcpp::IntegerVector::get_na ());
164 Rcpp::IntegerVector index (npts);
165 // Create parallel worker
166 OnePointIndex one_pt_indx (RcppParallel::RVector <double> (vtx),
167 RcppParallel::RVector <double> (vty),
168 RcppParallel::RVector <double> (ptx),
169 RcppParallel::RVector <double> (pty), nxy, index);
170
171 RcppParallel::parallelFor (0, npts, one_pt_indx);
172
173 return index;
174}
175
176//' rcpp_points_to_edges_par
177//'
178//' Get index of nearest edges to list of points
179//'
180//' @param graph Rcpp::DataFrame containing the full edge-based graph
181//' @param pts Rcpp::DataFrame containing the points to be matched
182//'
183//' @return 0-indexed Rcpp::NumericVector index into graph of nearest points
184//'
185//' @noRd
186// [[Rcpp::export]]
187Rcpp::NumericVector rcpp_points_to_edges_par (const Rcpp::DataFrame &graph,
188 Rcpp::DataFrame &pts)
189{
190 Rcpp::NumericVector ptx = pts ["x"];
191 Rcpp::NumericVector pty = pts ["y"];
192
193 Rcpp::NumericVector xfr = graph ["xfr"];
194 Rcpp::NumericVector yfr = graph ["yfr"];
195 Rcpp::NumericVector xto = graph ["xto"];
196 Rcpp::NumericVector yto = graph ["yto"];
197
198 const size_t nxy = static_cast <size_t> (graph.nrow ()),
199 npts = static_cast <size_t> (pts.nrow ());
200
201 // index holds three vectors:
202 // 1. the integer index
203 // 2. the x coordinates of the intersection points
204 // 3. the y coordinates of the intersection points
205 Rcpp::NumericVector index (npts * 3L);
206
207 // Create parallel worker
208 OneEdgeIndex one_edge_indx (RcppParallel::RVector <double> (ptx),
209 RcppParallel::RVector <double> (pty),
210 RcppParallel::RVector <double> (xfr),
211 RcppParallel::RVector <double> (yfr),
212 RcppParallel::RVector <double> (xto),
213 RcppParallel::RVector <double> (yto),
214 nxy, npts, index);
215
216 RcppParallel::parallelFor (0, npts, one_edge_indx);
217
218 return index;
219}