1 #pragma once
2 /*
3 * Copyright 2020 Axel Waggershauser
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 #include "Point.h"
19
20 #include <algorithm>
21 #include <cmath>
22 #include <numeric>
23 #include <vector>
24
25 namespace ZXing {
26
27 class RegressionLine
28 {
29 protected:
30 std::vector<PointF> _points;
31 PointF _directionInward;
32 PointF::value_t a = NAN, b = NAN, c = NAN;
33
34 friend PointF intersect(const RegressionLine& l1, const RegressionLine& l2);
35
evaluate(const std::vector<PointF> & ps)36 bool evaluate(const std::vector<PointF>& ps)
37 {
38 auto mean = std::accumulate(ps.begin(), ps.end(), PointF()) / ps.size();
39 PointF::value_t sumXX = 0, sumYY = 0, sumXY = 0;
40 for (auto& p : ps) {
41 auto d = p - mean;
42 sumXX += d.x * d.x;
43 sumYY += d.y * d.y;
44 sumXY += d.x * d.y;
45 }
46 if (sumYY >= sumXX) {
47 auto l = std::sqrt(sumYY * sumYY + sumXY * sumXY);
48 a = +sumYY / l;
49 b = -sumXY / l;
50 } else {
51 auto l = std::sqrt(sumXX * sumXX + sumXY * sumXY);
52 a = +sumXY / l;
53 b = -sumXX / l;
54 }
55 if (dot(_directionInward, normal()) < 0) {
56 a = -a;
57 b = -b;
58 }
59 c = dot(normal(), mean); // (a*mean.x + b*mean.y);
60 return dot(_directionInward, normal()) > 0.5f; // angle between original and new direction is at most 60 degree
61 }
62
63 public:
RegressionLine()64 RegressionLine() { _points.reserve(16); } // arbitrary but plausible start size (tiny performance improvement)
65
points()66 const auto& points() const { return _points; }
length()67 int length() const { return _points.size() >= 2 ? int(distance(_points.front(), _points.back())) : 0; }
isValid()68 bool isValid() const { return !std::isnan(a); }
normal()69 PointF normal() const { return isValid() ? PointF(a, b) : _directionInward; }
signedDistance(PointF p)70 auto signedDistance(PointF p) const { return dot(normal(), p) - c; }
project(PointF p)71 PointF project(PointF p) const { return p - signedDistance(p) * normal(); }
72
reset()73 void reset()
74 {
75 _points.clear();
76 _directionInward = {};
77 a = b = c = NAN;
78 }
79
add(PointF p)80 void add(PointF p) {
81 assert(_directionInward != PointF());
82 _points.push_back(p);
83 if (_points.size() == 1)
84 c = dot(normal(), p);
85 }
86
pop_back()87 void pop_back() { _points.pop_back(); }
88
setDirectionInward(PointF d)89 void setDirectionInward(PointF d) { _directionInward = normalized(d); }
90
91 bool evaluate(double maxSignedDist = -1, bool updatePoints = false)
92 {
93 bool ret = evaluate(_points);
94 if (maxSignedDist > 0) {
95 auto points = _points;
96 while (true) {
97 auto old_points_size = points.size();
98 points.erase(
99 std::remove_if(points.begin(), points.end(),
100 [this, maxSignedDist](auto p) { return this->signedDistance(p) > maxSignedDist; }),
101 points.end());
102 if (old_points_size == points.size())
103 break;
104 #ifdef PRINT_DEBUG
105 printf("removed %zu points\n", old_points_size - points.size());
106 #endif
107 ret = evaluate(points);
108 }
109
110 if (updatePoints)
111 _points = std::move(points);
112 }
113 return ret;
114 }
115
isHighRes()116 bool isHighRes() const
117 {
118 PointF min = _points.front(), max = _points.front();
119 for (auto p : _points) {
120 min.x = std::min(min.x, p.x);
121 min.y = std::min(min.y, p.y);
122 max.x = std::max(max.x, p.x);
123 max.y = std::max(max.y, p.y);
124 }
125 auto diff = max - min;
126 auto len = maxAbsComponent(diff);
127 auto steps = std::min(std::abs(diff.x), std::abs(diff.y));
128 // due to aliasing we get bad extrapolations if the line is short and too close to vertical/horizontal
129 return steps > 2 || len > 50;
130 }
131 };
132
intersect(const RegressionLine & l1,const RegressionLine & l2)133 inline PointF intersect(const RegressionLine& l1, const RegressionLine& l2)
134 {
135 assert(l1.isValid() && l2.isValid());
136 auto d = l1.a * l2.b - l1.b * l2.a;
137 auto x = (l1.c * l2.b - l1.b * l2.c) / d;
138 auto y = (l1.a * l2.c - l1.c * l2.a) / d;
139 return {x, y};
140 }
141
142 } // ZXing
143
144