13
13
14
14
using namespace cpp11 ::literals;
15
15
16
- class Point {
16
+ class Point2 {
17
17
public:
18
18
double x;
19
19
double y;
20
20
21
- Point () : x(0.0 ), y(0.0 ) {};
22
- Point (double _x, double _y) : x(_x), y(_y) {};
23
- Point (const Point & p) : x(p.x), y(p.y) {};
21
+ Point2 () : x(0.0 ), y(0.0 ) {};
22
+ Point2 (double _x, double _y) : x(_x), y(_y) {};
23
+ Point2 (const Point2 & p) : x(p.x), y(p.y) {};
24
24
25
- Point operator +(const Point & p) const {
25
+ Point2 operator +(const Point2 & p) const {
26
26
return {x + p.x , y + p.y };
27
27
};
28
28
29
- Point & operator +=(const Point & p) {
29
+ Point2 & operator +=(const Point2 & p) {
30
30
x += p.x ;
31
31
y += p.y ;
32
32
return *this ;
33
33
};
34
34
35
- Point operator -(const Point & p) const {
35
+ Point2 operator -(const Point2 & p) const {
36
36
return {x - p.x , y - p.y };
37
37
};
38
38
39
- Point & operator -=(const Point & p) {
39
+ Point2 & operator -=(const Point2 & p) {
40
40
x -= p.x ;
41
41
y -= p.y ;
42
42
return *this ;
43
43
};
44
44
45
- Point operator *(double fac) const {
45
+ Point2 operator *(double fac) const {
46
46
return {x * fac, y * fac};
47
47
};
48
48
49
- Point & operator *=(double fac) {
49
+ Point2 & operator *=(double fac) {
50
50
x *= fac;
51
51
y *= fac;
52
52
return *this ;
53
53
};
54
54
55
- double operator *(const Point & p) const {
55
+ double operator *(const Point2 & p) const {
56
56
return dot (p);
57
57
};
58
58
59
- Point operator /(double fac) const {
59
+ Point2 operator /(double fac) const {
60
60
return {x / fac, y / fac};
61
61
};
62
62
63
- Point & operator /=(double fac) {
63
+ Point2 & operator /=(double fac) {
64
64
x /= fac;
65
65
y /= fac;
66
66
return *this ;
67
67
};
68
68
69
- bool operator ==(const Point & p) const {
69
+ bool operator ==(const Point2 & p) const {
70
70
return x == p.x && y == p.y ;
71
71
};
72
72
73
- bool operator !=(const Point & p) const {
73
+ bool operator !=(const Point2 & p) const {
74
74
return x != p.x || y != p.y ;
75
75
};
76
76
77
- double distance_to (const Point & p, double eps = 0 ) const {
77
+ double distance_to (const Point2 & p, double eps = 0 ) const {
78
78
double dx = x - p.x ;
79
79
double dy = y - p.y ;
80
80
if (std::abs (dx) < eps && std::abs (dy) < eps) return eps;
81
81
82
82
return std::sqrt (dx * dx + dy * dy);
83
83
}
84
84
85
- double dot (const Point & p) const {
85
+ double dot (const Point2 & p) const {
86
86
return x * p.x + y * p.y ;
87
87
}
88
88
};
89
89
class Segment {
90
90
public:
91
- Point a;
92
- Point b;
91
+ Point2 a;
92
+ Point2 b;
93
93
94
94
Segment () : a(), b() {};
95
- Segment (Point _a, Point _b) : a(_a), b(_b) {};
95
+ Segment (Point2 _a, Point2 _b) : a(_a), b(_b) {};
96
96
Segment (double x0, double y0, double x1, double y1) : a(x0, y0), b(x1, y1) {};
97
97
Segment (const Segment& s) : a(s.a), b(s.b) {};
98
98
99
99
double length (double eps = 0.0 ) const {
100
100
return a.distance_to (b, eps);
101
101
}
102
102
103
- Point as_vec () const {
103
+ Point2 as_vec () const {
104
104
return b - a;
105
105
}
106
106
107
- Point project (const Point & p, double eps = 0.0 ) const {
108
- Point vec = as_vec ();
107
+ Point2 project (const Point2 & p, double eps = 0.0 ) const {
108
+ Point2 vec = as_vec ();
109
109
110
110
double l2 = vec * vec;
111
111
double u = ((p.x - a.x ) * vec.x + (p.y - a.y ) * vec.y ) / l2;
@@ -117,7 +117,7 @@ class Segment {
117
117
return {project (s.a , eps), project (s.b , eps)};
118
118
}
119
119
120
- Point midpoint () const {
120
+ Point2 midPoint2 () const {
121
121
return (a + b) / 2.0 ;
122
122
}
123
123
@@ -128,7 +128,7 @@ class Segment {
128
128
private:
129
129
double visibility (const Segment& s, double eps = 0.0 ) const {
130
130
Segment proj_s = project (s, eps);
131
- return std::max (1.0 - 2.0 * midpoint ().distance_to (proj_s.midpoint (), eps) / proj_s.length (eps), 0.0 );
131
+ return std::max (1.0 - 2.0 * midPoint2 ().distance_to (proj_s.midPoint2 (), eps) / proj_s.length (eps), 0.0 );
132
132
}
133
133
double angle_comp (const Segment& Q, double eps = 0.0 ) const {
134
134
double dot_PQ = as_vec () * Q.as_vec ();
@@ -148,7 +148,7 @@ class Segment {
148
148
149
149
double lavg = (euc_P + euc_Q) / 2.0 ;
150
150
151
- double euc_mid = midpoint ().distance_to (Q.midpoint (), eps);
151
+ double euc_mid = midPoint2 ().distance_to (Q.midPoint2 (), eps);
152
152
return lavg / (lavg + euc_mid);
153
153
}
154
154
double visibility_comp (const Segment& Q, double eps = 0.0 ) const {
@@ -159,7 +159,7 @@ class Segment {
159
159
}
160
160
};
161
161
162
- typedef std::vector<Point > Path;
162
+ typedef std::vector<Point2 > Path;
163
163
164
164
165
165
double compute_divided_edge_length (Path& edge) {
@@ -174,7 +174,7 @@ double compute_divided_edge_length(Path& edge) {
174
174
void update_edge_divisions (std::vector<Path>& elist, int P) {
175
175
for (size_t e_idx = 0 ; e_idx < elist.size (); ++e_idx) {
176
176
if (P == 1 ) {
177
- elist[e_idx].insert (elist[e_idx].begin () + 1 , Segment (elist[e_idx][0 ], elist[e_idx][1 ]).midpoint ());
177
+ elist[e_idx].insert (elist[e_idx].begin () + 1 , Segment (elist[e_idx][0 ], elist[e_idx][1 ]).midPoint2 ());
178
178
} else {
179
179
Path edge = elist[e_idx];
180
180
double current_edge_length = compute_divided_edge_length (edge);
@@ -184,19 +184,19 @@ void update_edge_divisions(std::vector<Path>& elist, int P) {
184
184
elist[e_idx].reserve (P + 2 );
185
185
elist[e_idx].push_back (edge[0 ]);
186
186
for (size_t i = 1 ; i < edge.size (); ++i) {
187
- Point& start_point = edge[i - 1 ];
188
- Point vec = edge[i] - start_point ;
189
- double segment_length = start_point .distance_to (edge[i]);
187
+ Point2& start_Point2 = edge[i - 1 ];
188
+ Point2 vec = edge[i] - start_Point2 ;
189
+ double segment_length = start_Point2 .distance_to (edge[i]);
190
190
while (segment_length > current_segment_length) {
191
191
double percent_position = current_segment_length / segment_length;
192
192
193
- elist[e_idx].push_back (start_point + vec * percent_position);
193
+ elist[e_idx].push_back (start_Point2 + vec * percent_position);
194
194
195
195
current_segment_length += new_segment_length;
196
196
}
197
197
current_segment_length -= segment_length;
198
198
}
199
- while (elist[e_idx].size () > P + 1 ) {
199
+ while (elist[e_idx].size () > size_t ( P + 1 ) ) {
200
200
elist[e_idx].pop_back ();
201
201
}
202
202
elist[e_idx].push_back (edge.back ());
@@ -220,16 +220,16 @@ std::vector< std::vector<int> > compute_compatibility_lists(std::vector<Path>& e
220
220
return comp;
221
221
}
222
222
223
- Point apply_spring_force (const std::vector<Path>& edges, int e_idx, int i, double kP ) {
223
+ Point2 apply_spring_force (const std::vector<Path>& edges, int e_idx, int i, double kP ) {
224
224
const Path& edge = edges[e_idx];
225
225
226
226
return (edge[i - 1 ] + edge[i + 1 ] - edge[i] * 2 ) * kP ;
227
227
}
228
228
229
- Point apply_electrostatic_force (const std::vector<Path>& edges,
229
+ Point2 apply_electrostatic_force (const std::vector<Path>& edges,
230
230
const std::vector< std::vector<int > >& comp,
231
231
int e_idx, int i, double eps) {
232
- Point sum_of_forces;
232
+ Point2 sum_of_forces;
233
233
if (comp[e_idx].empty ()) {
234
234
return sum_of_forces;
235
235
}
@@ -238,7 +238,7 @@ Point apply_electrostatic_force(const std::vector<Path>& edges,
238
238
239
239
for (size_t oe = 0 ; oe < compatible_edges.size (); ++oe) {
240
240
const Path& edge2 = edges[compatible_edges[oe]];
241
- Point force = edge2[i] - edge[i];
241
+ Point2 force = edge2[i] - edge[i];
242
242
if ((std::abs (force.x ) > eps) || (std::abs (force.y ) > eps)) {
243
243
double euc = edge2[i].distance_to (edge[i]);
244
244
double diff = std::pow (euc, -1.0 );
@@ -249,7 +249,7 @@ Point apply_electrostatic_force(const std::vector<Path>& edges,
249
249
return sum_of_forces;
250
250
}
251
251
252
- std::vector<Point> apply_resulting_forces_on_subdivision_points (const std::vector<Path>& edges,
252
+ std::vector<Point2> apply_resulting_forces_on_subdivision_Point2s (const std::vector<Path>& edges,
253
253
const std::vector< std::vector<int > >& comp,
254
254
int e_idx, int P,
255
255
double S, double K,
@@ -258,16 +258,16 @@ std::vector<Point> apply_resulting_forces_on_subdivision_points(const std::vecto
258
258
259
259
double kP = K / (edge[0 ].distance_to (edge[P + 1 ], eps) * (P + 1 ));
260
260
261
- std::vector<Point> resulting_forces_for_subdivision_points (P + 2 );
261
+ std::vector<Point2> resulting_forces_for_subdivision_Point2s (P + 2 );
262
262
for (int i = 1 ; i < (P + 1 ); ++i) {
263
- Point spring_force = apply_spring_force (edges, e_idx, i, kP );
263
+ Point2 spring_force = apply_spring_force (edges, e_idx, i, kP );
264
264
265
- Point electrostatic_force = apply_electrostatic_force (edges, comp, e_idx, i, eps);
265
+ Point2 electrostatic_force = apply_electrostatic_force (edges, comp, e_idx, i, eps);
266
266
267
- resulting_forces_for_subdivision_points [i] = (spring_force + electrostatic_force) * S;
267
+ resulting_forces_for_subdivision_Point2s [i] = (spring_force + electrostatic_force) * S;
268
268
}
269
269
270
- return resulting_forces_for_subdivision_points ;
270
+ return resulting_forces_for_subdivision_Point2s ;
271
271
}
272
272
273
273
@@ -279,8 +279,8 @@ cpp11::writable::data_frame force_bundle_iter(cpp11::doubles_matrix<> edges_xy,
279
279
// first division
280
280
std::vector<Path> edges (m);
281
281
for (R_xlen_t i = 0 ; i < m; ++i) {
282
- edges[i].push_back (Point (edges_xy (i, 0 ), edges_xy (i, 1 )));
283
- edges[i].push_back (Point (edges_xy (i, 2 ), edges_xy (i, 3 )));
282
+ edges[i].push_back (Point2 (edges_xy (i, 0 ), edges_xy (i, 1 )));
283
+ edges[i].push_back (Point2 (edges_xy (i, 2 ), edges_xy (i, 3 )));
284
284
}
285
285
286
286
// compute compatibility list
@@ -291,9 +291,9 @@ cpp11::writable::data_frame force_bundle_iter(cpp11::doubles_matrix<> edges_xy,
291
291
// main loop
292
292
for (int cycle = 0 ; cycle < C; ++cycle) {
293
293
for (int iteration = 0 ; iteration < I; ++iteration) {
294
- std::vector< std::vector<Point > > forces (m);
294
+ std::vector< std::vector<Point2 > > forces (m);
295
295
for (R_xlen_t e = 0 ; e < m; ++e) {
296
- forces[e] = apply_resulting_forces_on_subdivision_points (edges, comp, e, P, S, K, eps);
296
+ forces[e] = apply_resulting_forces_on_subdivision_Point2s (edges, comp, e, P, S, K, eps);
297
297
}
298
298
for (int e = 0 ; e < m; ++e) {
299
299
for (int i = 0 ; i < (P + 1 ); ++i) {
0 commit comments