/src/tesseract/src/ccstruct/detlinefit.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | /////////////////////////////////////////////////////////////////////// |
2 | | // File: detlinefit.cpp |
3 | | // Description: Deterministic least median squares line fitting. |
4 | | // Author: Ray Smith |
5 | | // |
6 | | // (C) Copyright 2008, Google Inc. |
7 | | // Licensed under the Apache License, Version 2.0 (the "License"); |
8 | | // you may not use this file except in compliance with the License. |
9 | | // You may obtain a copy of the License at |
10 | | // http://www.apache.org/licenses/LICENSE-2.0 |
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 | | |
19 | | #include "detlinefit.h" |
20 | | #include "helpers.h" // for IntCastRounded |
21 | | #include "statistc.h" |
22 | | #include "tesserrstream.h" // for tesserr |
23 | | |
24 | | #include <algorithm> |
25 | | #include <cfloat> // for FLT_MAX |
26 | | |
27 | | namespace tesseract { |
28 | | |
29 | | // The number of points to consider at each end. |
30 | | const int kNumEndPoints = 3; |
31 | | // The minimum number of points at which to switch to number of points |
32 | | // for badly fitted lines. |
33 | | // To ensure a sensible error metric, kMinPointsForErrorCount should be at |
34 | | // least kMaxRealDistance / (1 - %ile) where %ile is the fractile used in |
35 | | // ComputeUpperQuartileError. |
36 | | const int kMinPointsForErrorCount = 16; |
37 | | // The maximum real distance to use before switching to number of |
38 | | // mis-fitted points, which will get square-rooted for true distance. |
39 | | const int kMaxRealDistance = 2.0; |
40 | | |
41 | 228k | DetLineFit::DetLineFit() : square_length_(0.0) {} |
42 | | |
43 | | // Delete all Added points. |
44 | 56.3k | void DetLineFit::Clear() { |
45 | 56.3k | pts_.clear(); |
46 | 56.3k | distances_.clear(); |
47 | 56.3k | } |
48 | | |
49 | | // Add a new point. Takes a copy - the pt doesn't need to stay in scope. |
50 | 644k | void DetLineFit::Add(const ICOORD &pt) { |
51 | 644k | pts_.emplace_back(pt, 0); |
52 | 644k | } |
53 | | // Associates a half-width with the given point if a point overlaps the |
54 | | // previous point by more than half the width, and its distance is further |
55 | | // than the previous point, then the more distant point is ignored in the |
56 | | // distance calculation. Useful for ignoring i dots and other diacritics. |
57 | 522k | void DetLineFit::Add(const ICOORD &pt, int halfwidth) { |
58 | 522k | pts_.emplace_back(pt, halfwidth); |
59 | 522k | } |
60 | | |
61 | | // Fits a line to the points, ignoring the skip_first initial points and the |
62 | | // skip_last final points, returning the fitted line as a pair of points, |
63 | | // and the upper quartile error. |
64 | 134k | double DetLineFit::Fit(int skip_first, int skip_last, ICOORD *pt1, ICOORD *pt2) { |
65 | | // Do something sensible with no points. |
66 | 134k | if (pts_.empty()) { |
67 | 0 | pt1->set_x(0); |
68 | 0 | pt1->set_y(0); |
69 | 0 | *pt2 = *pt1; |
70 | 0 | return 0.0; |
71 | 0 | } |
72 | | // Count the points and find the first and last kNumEndPoints. |
73 | 134k | int pt_count = pts_.size(); |
74 | 134k | ICOORD *starts[kNumEndPoints]; |
75 | 134k | if (skip_first >= pt_count) { |
76 | 0 | skip_first = pt_count - 1; |
77 | 0 | } |
78 | 134k | int start_count = 0; |
79 | 134k | int end_i = std::min(skip_first + kNumEndPoints, pt_count); |
80 | 428k | for (int i = skip_first; i < end_i; ++i) { |
81 | 294k | starts[start_count++] = &pts_[i].pt; |
82 | 294k | } |
83 | 134k | ICOORD *ends[kNumEndPoints]; |
84 | 134k | if (skip_last >= pt_count) { |
85 | 0 | skip_last = pt_count - 1; |
86 | 0 | } |
87 | 134k | int end_count = 0; |
88 | 134k | end_i = std::max(0, pt_count - kNumEndPoints - skip_last); |
89 | 428k | for (int i = pt_count - 1 - skip_last; i >= end_i; --i) { |
90 | 294k | ends[end_count++] = &pts_[i].pt; |
91 | 294k | } |
92 | | // 1 or 2 points need special treatment. |
93 | 134k | if (pt_count <= 2) { |
94 | 63.7k | *pt1 = *starts[0]; |
95 | 63.7k | if (pt_count > 1) { |
96 | 19.7k | *pt2 = *ends[0]; |
97 | 44.0k | } else { |
98 | 44.0k | *pt2 = *pt1; |
99 | 44.0k | } |
100 | 63.7k | return 0.0; |
101 | 63.7k | } |
102 | | // Although with between 2 and 2*kNumEndPoints-1 points, there will be |
103 | | // overlap in the starts, ends sets, this is OK and taken care of by the |
104 | | // if (*start != *end) test below, which also tests for equal input points. |
105 | 70.3k | double best_uq = -1.0; |
106 | | // Iterate each pair of points and find the best fitting line. |
107 | 281k | for (int i = 0; i < start_count; ++i) { |
108 | 211k | ICOORD *start = starts[i]; |
109 | 844k | for (int j = 0; j < end_count; ++j) { |
110 | 633k | ICOORD *end = ends[j]; |
111 | 633k | if (*start != *end) { |
112 | 572k | ComputeDistances(*start, *end); |
113 | | // Compute the upper quartile error from the line. |
114 | 572k | double dist = EvaluateLineFit(); |
115 | 572k | if (dist < best_uq || best_uq < 0.0) { |
116 | 144k | best_uq = dist; |
117 | 144k | *pt1 = *start; |
118 | 144k | *pt2 = *end; |
119 | 144k | } |
120 | 572k | } |
121 | 633k | } |
122 | 211k | } |
123 | | // Finally compute the square root to return the true distance. |
124 | 70.3k | return best_uq > 0.0 ? sqrt(best_uq) : best_uq; |
125 | 134k | } |
126 | | |
127 | | // Constrained fit with a supplied direction vector. Finds the best line_pt, |
128 | | // that is one of the supplied points having the median cross product with |
129 | | // direction, ignoring points that have a cross product outside of the range |
130 | | // [min_dist, max_dist]. Returns the resulting error metric using the same |
131 | | // reduced set of points. |
132 | | // *Makes use of floating point arithmetic* |
133 | | double DetLineFit::ConstrainedFit(const FCOORD &direction, double min_dist, double max_dist, |
134 | 218k | bool debug, ICOORD *line_pt) { |
135 | 218k | ComputeConstrainedDistances(direction, min_dist, max_dist); |
136 | | // Do something sensible with no points or computed distances. |
137 | 218k | if (pts_.empty() || distances_.empty()) { |
138 | 0 | line_pt->set_x(0); |
139 | 0 | line_pt->set_y(0); |
140 | 0 | return 0.0; |
141 | 0 | } |
142 | 218k | auto median_index = distances_.size() / 2; |
143 | 218k | std::nth_element(distances_.begin(), distances_.begin() + median_index, distances_.end()); |
144 | 218k | *line_pt = distances_[median_index].data(); |
145 | 218k | if (debug) { |
146 | 0 | tesserr << "Constrained fit to dir " << direction.x() << ", " |
147 | 0 | << direction.y() << " = " |
148 | 0 | << line_pt->x() << ", " << line_pt->y() |
149 | 0 | << " :" << distances_.size() << " distances:\n"; |
150 | 0 | for (unsigned i = 0; i < distances_.size(); ++i) { |
151 | 0 | tesserr << i << ": " |
152 | 0 | << distances_[i].data().x() << ", " |
153 | 0 | << distances_[i].data().y() << " -> " |
154 | 0 | << distances_[i].key() << '\n'; |
155 | 0 | } |
156 | 0 | tesserr << "Result = " << median_index << '\n'; |
157 | 0 | } |
158 | | // Center distances on the fitted point. |
159 | 218k | double dist_origin = direction * *line_pt; |
160 | 910k | for (auto &distance : distances_) { |
161 | 910k | distance.key() -= dist_origin; |
162 | 910k | } |
163 | 218k | return sqrt(EvaluateLineFit()); |
164 | 218k | } |
165 | | |
166 | | // Returns true if there were enough points at the last call to Fit or |
167 | | // ConstrainedFit for the fitted points to be used on a badly fitted line. |
168 | 116k | bool DetLineFit::SufficientPointsForIndependentFit() const { |
169 | 116k | return distances_.size() >= kMinPointsForErrorCount; |
170 | 116k | } |
171 | | |
172 | | // Backwards compatible fit returning a gradient and constant. |
173 | | // Deprecated. Prefer Fit(ICOORD*, ICOORD*) where possible, but use this |
174 | | // function in preference to the LMS class. |
175 | 69.9k | double DetLineFit::Fit(float *m, float *c) { |
176 | 69.9k | ICOORD start, end; |
177 | 69.9k | double error = Fit(&start, &end); |
178 | 69.9k | if (end.x() != start.x()) { |
179 | 36.8k | *m = static_cast<float>(end.y() - start.y()) / (end.x() - start.x()); |
180 | 36.8k | *c = start.y() - *m * start.x(); |
181 | 36.8k | } else { |
182 | 33.0k | *m = 0.0f; |
183 | 33.0k | *c = 0.0f; |
184 | 33.0k | } |
185 | 69.9k | return error; |
186 | 69.9k | } |
187 | | |
188 | | // Backwards compatible constrained fit with a supplied gradient. |
189 | | // Deprecated. Use ConstrainedFit(const FCOORD& direction) where possible |
190 | | // to avoid potential difficulties with infinite gradients. |
191 | 102k | double DetLineFit::ConstrainedFit(double m, float *c) { |
192 | | // Do something sensible with no points. |
193 | 102k | if (pts_.empty()) { |
194 | 0 | *c = 0.0f; |
195 | 0 | return 0.0; |
196 | 0 | } |
197 | 102k | double cos = 1.0 / sqrt(1.0 + m * m); |
198 | 102k | FCOORD direction(cos, m * cos); |
199 | 102k | ICOORD line_pt; |
200 | 102k | double error = ConstrainedFit(direction, -FLT_MAX, FLT_MAX, false, &line_pt); |
201 | 102k | *c = line_pt.y() - line_pt.x() * m; |
202 | 102k | return error; |
203 | 102k | } |
204 | | |
205 | | // Computes and returns the squared evaluation metric for a line fit. |
206 | 791k | double DetLineFit::EvaluateLineFit() { |
207 | | // Compute the upper quartile error from the line. |
208 | 791k | double dist = ComputeUpperQuartileError(); |
209 | 791k | if (distances_.size() >= kMinPointsForErrorCount && dist > kMaxRealDistance * kMaxRealDistance) { |
210 | | // Use the number of mis-fitted points as the error metric, as this |
211 | | // gives a better measure of fit for badly fitted lines where more |
212 | | // than a quarter are badly fitted. |
213 | 138k | double threshold = kMaxRealDistance * sqrt(square_length_); |
214 | 138k | dist = NumberOfMisfittedPoints(threshold); |
215 | 138k | } |
216 | 791k | return dist; |
217 | 791k | } |
218 | | |
219 | | // Computes the absolute error distances of the points from the line, |
220 | | // and returns the squared upper-quartile error distance. |
221 | 791k | double DetLineFit::ComputeUpperQuartileError() { |
222 | 791k | int num_errors = distances_.size(); |
223 | 791k | if (num_errors == 0) { |
224 | 0 | return 0.0; |
225 | 0 | } |
226 | | // Get the absolute values of the errors. |
227 | 9.72M | for (int i = 0; i < num_errors; ++i) { |
228 | 8.92M | if (distances_[i].key() < 0) { |
229 | 3.26M | distances_[i].key() = -distances_[i].key(); |
230 | 3.26M | } |
231 | 8.92M | } |
232 | | // Now get the upper quartile distance. |
233 | 791k | auto index = 3 * num_errors / 4; |
234 | 791k | std::nth_element(distances_.begin(), distances_.begin() + index, distances_.end()); |
235 | 791k | double dist = distances_[index].key(); |
236 | | // The true distance is the square root of the dist squared / square_length. |
237 | | // Don't bother with the square root. Just return the square distance. |
238 | 791k | return square_length_ > 0.0 ? dist * dist / square_length_ : 0.0; |
239 | 791k | } |
240 | | |
241 | | // Returns the number of sample points that have an error more than threshold. |
242 | 138k | int DetLineFit::NumberOfMisfittedPoints(double threshold) const { |
243 | 138k | int num_misfits = 0; |
244 | 138k | int num_dists = distances_.size(); |
245 | | // Get the absolute values of the errors. |
246 | 4.87M | for (int i = 0; i < num_dists; ++i) { |
247 | 4.73M | if (distances_[i].key() > threshold) { |
248 | 3.00M | ++num_misfits; |
249 | 3.00M | } |
250 | 4.73M | } |
251 | 138k | return num_misfits; |
252 | 138k | } |
253 | | |
254 | | // Computes all the cross product distances of the points from the line, |
255 | | // storing the actual (signed) cross products in distances. |
256 | | // Ignores distances of points that are further away than the previous point, |
257 | | // and overlaps the previous point by at least half. |
258 | 572k | void DetLineFit::ComputeDistances(const ICOORD &start, const ICOORD &end) { |
259 | 572k | distances_.clear(); |
260 | 572k | ICOORD line_vector = end; |
261 | 572k | line_vector -= start; |
262 | 572k | square_length_ = line_vector.sqlength(); |
263 | 572k | int line_length = IntCastRounded(sqrt(square_length_)); |
264 | | // Compute the distance of each point from the line. |
265 | 572k | int prev_abs_dist = 0; |
266 | 572k | int prev_dot = 0; |
267 | 8.96M | for (unsigned i = 0; i < pts_.size(); ++i) { |
268 | 8.39M | ICOORD pt_vector = pts_[i].pt; |
269 | 8.39M | pt_vector -= start; |
270 | 8.39M | int dot = line_vector % pt_vector; |
271 | | // Compute |line_vector||pt_vector|sin(angle between) |
272 | 8.39M | int dist = line_vector * pt_vector; |
273 | 8.39M | int abs_dist = dist < 0 ? -dist : dist; |
274 | 8.39M | if (abs_dist > prev_abs_dist && i > 0) { |
275 | | // Ignore this point if it overlaps the previous one. |
276 | 3.63M | int separation = abs(dot - prev_dot); |
277 | 3.63M | if (separation < line_length * pts_[i].halfwidth || |
278 | 3.63M | separation < line_length * pts_[i - 1].halfwidth) { |
279 | 371k | continue; |
280 | 371k | } |
281 | 3.63M | } |
282 | 8.01M | distances_.emplace_back(dist, pts_[i].pt); |
283 | 8.01M | prev_abs_dist = abs_dist; |
284 | 8.01M | prev_dot = dot; |
285 | 8.01M | } |
286 | 572k | } |
287 | | |
288 | | // Computes all the cross product distances of the points perpendicular to |
289 | | // the given direction, ignoring distances outside of the give distance range, |
290 | | // storing the actual (signed) cross products in distances_. |
291 | | void DetLineFit::ComputeConstrainedDistances(const FCOORD &direction, double min_dist, |
292 | 218k | double max_dist) { |
293 | 218k | distances_.clear(); |
294 | 218k | square_length_ = direction.sqlength(); |
295 | | // Compute the distance of each point from the line. |
296 | 1.49M | for (auto &pt : pts_) { |
297 | 1.49M | FCOORD pt_vector = pt.pt; |
298 | | // Compute |line_vector||pt_vector|sin(angle between) |
299 | 1.49M | double dist = direction * pt_vector; |
300 | 1.49M | if (min_dist <= dist && dist <= max_dist) { |
301 | 910k | distances_.emplace_back(dist, pt.pt); |
302 | 910k | } |
303 | 1.49M | } |
304 | 218k | } |
305 | | |
306 | | } // namespace tesseract. |