Tesseract  3.02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  * Created: Wed Oct 6 15:14:23 BST 1993
6  *
7  * (C) Copyright 1993, Hewlett-Packard Ltd.
8  ** Licensed under the Apache License, Version 2.0 (the "License");
9  ** you may not use this file except in compliance with the License.
10  ** You may obtain a copy of the License at
11  ** http://www.apache.org/licenses/LICENSE-2.0
12  ** Unless required by applicable law or agreed to in writing, software
13  ** distributed under the License is distributed on an "AS IS" BASIS,
14  ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15  ** See the License for the specific language governing permissions and
16  ** limitations under the License.
17  *
18  **********************************************************************/
19 
20 #include "mfcpch.h"
21 #include <stdio.h>
22 #include <math.h>
23 #include "errcode.h"
24 #include "quadlsq.h"
25 
26 const ERRCODE EMPTY_QLSQ = "Can't delete from an empty QLSQ";
27 
28 #define EXTERN
29 
30 /**********************************************************************
31  * QLSQ::clear
32  *
33  * Function to initialize a QLSQ.
34  **********************************************************************/
35 
36 void QLSQ::clear() { //initialize
37  a = 0;
38  b = 0;
39  c = 0;
40  n = 0; //no elements
41  sigx = 0; //update accumulators
42  sigy = 0;
43  sigxx = 0;
44  sigxy = 0;
45  sigyy = 0;
46  sigxxx = 0;
47  sigxxy = 0;
48  sigxxxx = 0;
49 }
50 
51 
52 /**********************************************************************
53  * QLSQ::add
54  *
55  * Add an element to the accumulator.
56  **********************************************************************/
57 
58 void QLSQ::add( //add an element
59  double x, //xcoord
60  double y //ycoord
61  ) {
62  n++; //count elements
63  sigx += x; //update accumulators
64  sigy += y;
65  sigxx += x * x;
66  sigxy += x * y;
67  sigyy += y * y;
68  sigxxx += (long double) x *x * x;
69  sigxxy += (long double) x *x * y;
70  sigxxxx += (long double) x *x * x * x;
71 }
72 
73 
74 /**********************************************************************
75  * QLSQ::remove
76  *
77  * Delete an element from the acculuator.
78  **********************************************************************/
79 
80 void QLSQ::remove( //delete an element
81  double x, //xcoord
82  double y //ycoord
83  ) {
84  if (n <= 0)
85  //illegal
86  EMPTY_QLSQ.error ("QLSQ::remove", ABORT, NULL);
87  n--; //count elements
88  sigx -= x; //update accumulators
89  sigy -= y;
90  sigxx -= x * x;
91  sigxy -= x * y;
92  sigyy -= y * y;
93  sigxxx -= (long double) x *x * x;
94  sigxxy -= (long double) x *x * y;
95  sigxxxx -= (long double) x *x * x * x;
96 }
97 
98 
99 /**********************************************************************
100  * QLSQ::fit
101  *
102  * Fit the given degree of polynomial and store the result.
103  **********************************************************************/
104 
105 void QLSQ::fit( //fit polynomial
106  int degree //degree to fit
107  ) {
108  long double cubetemp; //intermediates
109  long double squaretemp;
110  long double top96, bottom96; /*accurate top & bottom */
111 
112  if (n >= 4 && degree >= 2) {
113  cubetemp = sigxxx * n - (long double) sigxx *sigx;
114 
115  top96 =
116  cubetemp * ((long double) sigxy * n - (long double) sigx * sigy);
117 
118  squaretemp = (long double) sigxx *n - (long double) sigx *sigx;
119 
120  top96 += squaretemp * ((long double) sigxx * sigy - sigxxy * n);
121 
122  bottom96 = cubetemp * cubetemp;
123 
124  bottom96 -= squaretemp * (sigxxxx * n - (long double) sigxx * sigxx);
125 
126  a = top96 / bottom96;
127 
128  top96 = ((long double) sigxx * sigx - sigxxx * n) * a
129  + (long double) sigxy *n - (long double) sigx *sigy;
130  bottom96 = (long double) sigxx *n - (long double) sigx *sigx;
131  b = top96 / bottom96;
132 
133  c = (sigy - a * sigxx - b * sigx) / n;
134  }
135  else if (n == 0 || degree < 0) {
136  a = b = c = 0;
137  }
138  else {
139  a = 0;
140  if (n > 1 && degree > 0) {
141  b = (sigxy * n - sigx * sigy) / (sigxx * n - sigx * sigx);
142  }
143  else
144  b = 0;
145  c = (sigy - b * sigx) / n;
146  }
147 }