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