tesseract 4.1.1
Loading...
Searching...
No Matches
quadlsq.cpp
Go to the documentation of this file.
1/**********************************************************************
2 * File: quadlsq.cpp (Formerly qlsq.c)
3 * Description: Code for least squares approximation of quadratics.
4 * Author: Ray Smith
5 *
6 * (C) Copyright 1993, Hewlett-Packard Ltd.
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 <cstdio>
20#include <cmath>
21#include "quadlsq.h"
22#include "tprintf.h"
23
24// Minimum variance in least squares before backing off to a lower degree.
25const long double kMinVariance = 1.0L / 1024;
26
27/**********************************************************************
28 * QLSQ::clear
29 *
30 * Function to initialize a QLSQ.
31 **********************************************************************/
32
33void QLSQ::clear() { // initialize
34 a = 0.0;
35 b = 0.0;
36 c = 0.0;
37 n = 0; // No elements.
38 sigx = 0.0; // Zero accumulators.
39 sigy = 0.0;
40 sigxx = 0.0;
41 sigxy = 0.0;
42 sigyy = 0.0;
43 sigxxx = 0.0;
44 sigxxy = 0.0;
45 sigxxxx = 0.0;
46}
47
48
49/**********************************************************************
50 * QLSQ::add
51 *
52 * Add an element to the accumulator.
53 **********************************************************************/
54
55void QLSQ::add(double x, double y) {
56 n++; // Count elements.
57 sigx += x; // Update accumulators.
58 sigy += y;
59 sigxx += x * x;
60 sigxy += x * y;
61 sigyy += y * y;
62 sigxxx += static_cast<long double>(x) * x * x;
63 sigxxy += static_cast<long double>(x) * x * y;
64 sigxxxx += static_cast<long double>(x) * x * x * x;
65}
66
67
68/**********************************************************************
69 * QLSQ::remove
70 *
71 * Delete an element from the accumulator.
72 **********************************************************************/
73
74void QLSQ::remove(double x, double y) {
75 if (n <= 0) {
76 tprintf("Can't remove an element from an empty QLSQ accumulator!\n");
77 return;
78 }
79 n--; // Count elements.
80 sigx -= x; // Update accumulators.
81 sigy -= y;
82 sigxx -= x * x;
83 sigxy -= x * y;
84 sigyy -= y * y;
85 sigxxx -= static_cast<long double>(x) * x * x;
86 sigxxy -= static_cast<long double>(x) * x * y;
87 sigxxxx -= static_cast<long double>(x) * x * x * x;
88}
89
90
91/**********************************************************************
92 * QLSQ::fit
93 *
94 * Fit the given degree of polynomial and store the result.
95 * This creates a quadratic of the form axx + bx + c, but limited to
96 * the given degree.
97 **********************************************************************/
98
99void QLSQ::fit(int degree) {
100 long double x_variance = static_cast<long double>(sigxx) * n -
101 static_cast<long double>(sigx) * sigx;
102
103 // Note: for computational efficiency, we do not normalize the variance,
104 // covariance and cube variance here as they are in the same order in both
105 // nominators and denominators. However, we need be careful in value range
106 // check.
107 if (x_variance < kMinVariance * n * n || degree < 1 || n < 2) {
108 // We cannot calculate b reliably so forget a and b, and just work on c.
109 a = b = 0.0;
110 if (n >= 1 && degree >= 0) {
111 c = sigy / n;
112 } else {
113 c = 0.0;
114 }
115 return;
116 }
117 long double top96 = 0.0; // Accurate top.
118 long double bottom96 = 0.0; // Accurate bottom.
119 long double cubevar = sigxxx * n - static_cast<long double>(sigxx) * sigx;
120 long double covariance = static_cast<long double>(sigxy) * n -
121 static_cast<long double>(sigx) * sigy;
122
123 if (n >= 4 && degree >= 2) {
124 top96 = cubevar * covariance;
125 top96 += x_variance * (static_cast<long double>(sigxx) * sigy - sigxxy * n);
126
127 bottom96 = cubevar * cubevar;
128 bottom96 -= x_variance *
129 (sigxxxx * n - static_cast<long double>(sigxx) * sigxx);
130 }
131 if (bottom96 >= kMinVariance * n * n * n * n) {
132 // Denominators looking good
133 a = top96 / bottom96;
134 top96 = covariance - cubevar * a;
135 b = top96 / x_variance;
136 } else {
137 // Forget a, and concentrate on b.
138 a = 0.0;
139 b = covariance / x_variance;
140 }
141 c = (sigy - a * sigxx - b * sigx) / n;
142}
const long double kMinVariance
Definition: quadlsq.cpp:25
DLLSYM void tprintf(const char *format,...)
Definition: tprintf.cpp:35
void remove(double x, double y)
Definition: quadlsq.cpp:74
void fit(int degree)
Definition: quadlsq.cpp:99
void clear()
Definition: quadlsq.cpp:33
void add(double x, double y)
Definition: quadlsq.cpp:55