1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
use crate::common::Result;
use crate::{Exceptions, Point};
use super::RegressionLineTrait;
#[derive(Clone)]
pub struct DMRegressionLine {
points: Vec<Point>,
direction_inward: Point,
pub(super) a: f32,
pub(super) b: f32,
pub(super) c: f32,
// std::vector<PointF> _points;
// PointF _directionInward;
// PointF::value_t a = NAN, b = NAN, c = NAN;
}
impl Default for DMRegressionLine {
fn default() -> Self {
Self {
points: Default::default(),
direction_inward: Default::default(),
a: f32::NAN,
b: f32::NAN,
c: f32::NAN,
}
}
}
impl RegressionLineTrait for DMRegressionLine {
fn points(&self) -> &[Point] {
&self.points
}
fn length(&self) -> u32 {
if self.points.len() >= 2 {
Point::distance(*self.points.first().unwrap(), *self.points.last().unwrap()) as u32
} else {
0
}
}
fn isValid(&self) -> bool {
!self.a.is_nan()
}
fn normal(&self) -> Point {
if self.isValid() {
Point {
x: self.a,
y: self.b,
}
} else {
self.direction_inward
}
}
fn signedDistance(&self, p: Point) -> f32 {
Point::dot(self.normal(), p) - self.c
}
fn distance_single(&self, p: Point) -> f32 {
(self.signedDistance(p)).abs()
}
fn reset(&mut self) {
self.points.clear();
self.direction_inward = Point { x: 0.0, y: 0.0 };
self.a = f32::NAN;
self.b = f32::NAN;
self.c = f32::NAN;
}
fn add(&mut self, p: Point) -> Result<()> {
if self.direction_inward == Point::default() {
return Err(Exceptions::ILLEGAL_STATE);
}
self.points.push(p);
if self.points.len() == 1 {
self.c = Point::dot(self.normal(), p);
}
Ok(())
}
fn pop_back(&mut self) {
self.points.pop();
}
fn setDirectionInward(&mut self, d: Point) {
self.direction_inward = Point::normalized(d);
}
fn evaluate_max_distance(
&mut self,
maxSignedDist: Option<f64>,
updatePoints: Option<bool>,
) -> bool {
let maxSignedDist = if let Some(m) = maxSignedDist { m } else { -1.0 };
let updatePoints = if let Some(u) = updatePoints { u } else { false };
let mut ret = self.evaluateSelf();
if maxSignedDist > 0.0 {
let mut points = self.points.clone();
loop {
let old_points_size = points.len();
// remove points that are further 'inside' than maxSignedDist or further 'outside' than 2 x maxSignedDist
// auto end = std::remove_if(points.begin(), points.end(), [this, maxSignedDist](auto p) {
// auto sd = this->signedDistance(p);
// return sd > maxSignedDist || sd < -2 * maxSignedDist;
// });
// points.erase(end, points.end());
points.retain(|&p| {
let sd = self.signedDistance(p) as f64;
!(sd > maxSignedDist || sd < -2.0 * maxSignedDist)
});
if old_points_size == points.len() {
break;
}
// #ifdef PRINT_DEBUG
// printf("removed %zu points\n", old_points_size - points.size());
// #endif
ret = self.evaluate(&points);
}
if updatePoints {
self.points = points;
}
}
ret
}
fn isHighRes(&self) -> bool {
let Some(mut min) = self.points.first().copied() else {
return false;
};
let Some(mut max) = self.points.first().copied() else {
return false;
};
for p in &self.points {
min.x = f32::min(min.x, p.x);
min.y = f32::min(min.y, p.y);
max.x = f32::max(max.x, p.x);
max.y = f32::max(max.y, p.y);
}
let diff = max - min;
let len = diff.maxAbsComponent();
let steps = f32::min(diff.x.abs(), diff.y.abs());
// due to aliasing we get bad extrapolations if the line is short and too close to vertical/horizontal
steps > 2.0 || len > 50.0
}
fn evaluate(&mut self, points: &[Point]) -> bool {
let mean = points.iter().sum::<Point>() / points.len() as f32;
let mut sumXX = 0.0;
let mut sumYY = 0.0;
let mut sumXY = 0.0;
for p in points {
// for (auto p = begin; p != end; ++p) {
let d = *p - mean;
sumXX += d.x * d.x;
sumYY += d.y * d.y;
sumXY += d.x * d.y;
}
if sumYY >= sumXX {
let l = (sumYY * sumYY + sumXY * sumXY).sqrt();
self.a = sumYY / l;
self.b = -sumXY / l;
} else {
let l = (sumXX * sumXX + sumXY * sumXY).sqrt();
self.a = sumXY / l;
self.b = -sumXX / l;
}
if Point::dot(self.direction_inward, self.normal()) < 0.0 {
// if (dot(_directionInward, normal()) < 0) {
self.a = -self.a;
self.b = -self.b;
}
self.c = Point::dot(self.normal(), mean); // (a*mean.x + b*mean.y);
Point::dot(self.direction_inward, self.normal()) > 0.5
// angle between original and new direction is at most 60 degree
}
fn evaluateSelf(&mut self) -> bool {
let mean = self.points.iter().sum::<Point>() / self.points.len() as f32;
let mut sumXX = 0.0;
let mut sumYY = 0.0;
let mut sumXY = 0.0;
for p in &self.points {
// for (auto p = begin; p != end; ++p) {
let d = *p - mean;
sumXX += d.x * d.x;
sumYY += d.y * d.y;
sumXY += d.x * d.y;
}
if sumYY >= sumXX {
let l = (sumYY * sumYY + sumXY * sumXY).sqrt();
self.a = sumYY / l;
self.b = -sumXY / l;
} else {
let l = (sumXX * sumXX + sumXY * sumXY).sqrt();
self.a = sumXY / l;
self.b = -sumXX / l;
}
if Point::dot(self.direction_inward, self.normal()) < 0.0 {
// if (dot(_directionInward, normal()) < 0) {
self.a = -self.a;
self.b = -self.b;
}
self.c = Point::dot(self.normal(), mean); // (a*mean.x + b*mean.y);
Point::dot(self.direction_inward, self.normal()) > 0.5
// angle between original and new direction is at most 60 degree
}
fn a(&self) -> f32 {
self.a
}
fn b(&self) -> f32 {
self.b
}
fn c(&self) -> f32 {
self.c
}
}
impl DMRegressionLine {
pub fn new(point_1: Point, point_2: Point) -> Self {
let mut new = Self::default();
RegressionLineTrait::evaluate(&mut new, &[point_1, point_2]);
new
}
// template <typename Container, typename Filter>
fn average<T>(c: &[f64], f: T) -> f64
where
T: Fn(f64) -> bool,
{
let mut sum: f64 = 0.0;
let mut num = 0;
for v in c {
// for (const auto& v : c)
if f(*v) {
sum += *v;
num += 1;
}
}
sum / num as f64
}
pub fn reverse(&mut self) {
self.points.reverse();
}
pub fn modules(&mut self, beg: Point, end: Point) -> Result<f64> {
if self.points.len() <= 3 {
return Err(Exceptions::ILLEGAL_STATE);
}
// re-evaluate and filter out all points too far away. required for the gapSizes calculation.
self.evaluate_max_distance(Some(1.0), Some(true));
// std::vector<double> gapSizes, modSizes;
let mut gapSizes: Vec<f64> = Vec::new();
let mut modSizes = Vec::new();
gapSizes.reserve(self.points.len());
// calculate the distance between the points projected onto the regression line
for i in 1..self.points.len() {
// for (size_t i = 1; i < _points.size(); ++i)
gapSizes.push(Point::distance(
self.project(self.points[i]),
self.project(self.points[i - 1]),
) as f64);
}
// calculate the (expected average) distance of two adjacent pixels
let unitPixelDist = Point::length(Point::bresenhamDirection(
self.points
.last()
.copied()
.ok_or(Exceptions::INDEX_OUT_OF_BOUNDS)?
- self
.points
.first()
.copied()
.ok_or(Exceptions::INDEX_OUT_OF_BOUNDS)?,
)) as f64;
// calculate the width of 2 modules (first black pixel to first black pixel)
let mut sumFront: f64 =
Point::distance(beg, self.project(self.points[0])) as f64 - unitPixelDist;
let mut sumBack: f64 = 0.0; // (last black pixel to last black pixel)
for dist in gapSizes {
// for (auto dist : gapSizes) {
if dist > 1.9 * unitPixelDist {
modSizes.push(std::mem::take(&mut sumBack));
}
sumFront += dist;
sumBack += dist;
if dist > 1.9 * unitPixelDist {
modSizes.push(std::mem::take(&mut sumFront));
}
}
modSizes.push(
sumFront
+ Point::distance(
end,
self.project(
self.points
.last()
.copied()
.ok_or(Exceptions::INDEX_OUT_OF_BOUNDS)?,
),
) as f64,
);
modSizes[0] = 0.0; // the first element is an invalid sumBack value, would be pop_front() if vector supported this
let lineLength = Point::distance(beg, end) as f64 - unitPixelDist;
let mut meanModSize = Self::average(&modSizes, |_: f64| true);
// let meanModSize = average(modSizes, [](double){ return true; });
// #ifdef PRINT_DEBUG
// printf("unit pixel dist: %.1f\n", unitPixelDist);
// printf("lineLength: %.1f, meanModSize: %.1f, gaps: %lu\n", lineLength, meanModSize, modSizes.size());
// #endif
for i in 0..2 {
// for (int i = 0; i < 2; ++i)
meanModSize = Self::average(&modSizes, |dist: f64| {
(dist - meanModSize).abs() < meanModSize / (2 + i) as f64
});
// meanModSize = average(modSizes, [=](double dist) { return std::abs(dist - meanModSize) < meanModSize / (2 + i); });
}
// #ifdef PRINT_DEBUG
// printf("post filter meanModSize: %.1f\n", meanModSize);
// #endif
Ok(lineLength / meanModSize)
}
}