Tesseract  3.02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
linlsq.cpp
Go to the documentation of this file.
1 /**********************************************************************
2  * File: linlsq.cpp (Formerly llsq.c)
3  * Description: Linear Least squares fitting code.
4  * Author: Ray Smith
5  * Created: Thu Sep 12 08:44:51 BST 1991
6  *
7  * (C) Copyright 1991, 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" // Must be first include for windows.
21 #include <stdio.h>
22 #include <math.h>
23 #include "errcode.h"
24 #include "linlsq.h"
25 
26 const ERRCODE EMPTY_LLSQ = "Can't delete from an empty LLSQ";
27 
28 /**********************************************************************
29  * LLSQ::clear
30  *
31  * Function to initialize a LLSQ.
32  **********************************************************************/
33 
34 void LLSQ::clear() { // initialize
35  total_weight = 0.0; // no elements
36  sigx = 0.0; // update accumulators
37  sigy = 0.0;
38  sigxx = 0.0;
39  sigxy = 0.0;
40  sigyy = 0.0;
41 }
42 
43 
44 /**********************************************************************
45  * LLSQ::add
46  *
47  * Add an element to the accumulator.
48  **********************************************************************/
49 
50 void LLSQ::add(double x, double y) { // add an element
51  total_weight++; // count elements
52  sigx += x; // update accumulators
53  sigy += y;
54  sigxx += x * x;
55  sigxy += x * y;
56  sigyy += y * y;
57 }
58 // Adds an element with a specified weight.
59 void LLSQ::add(double x, double y, double weight) {
60  total_weight += weight;
61  sigx += x * weight; // update accumulators
62  sigy += y * weight;
63  sigxx += x * x * weight;
64  sigxy += x * y * weight;
65  sigyy += y * y * weight;
66 }
67 // Adds a whole LLSQ.
68 void LLSQ::add(const LLSQ& other) {
69  total_weight += other.total_weight;
70  sigx += other.sigx; // update accumulators
71  sigy += other.sigy;
72  sigxx += other.sigxx;
73  sigxy += other.sigxy;
74  sigyy += other.sigyy;
75 }
76 
77 
78 /**********************************************************************
79  * LLSQ::remove
80  *
81  * Delete an element from the acculuator.
82  **********************************************************************/
83 
84 void LLSQ::remove(double x, double y) { // delete an element
85  if (total_weight <= 0.0) // illegal
86  EMPTY_LLSQ.error("LLSQ::remove", ABORT, NULL);
87  total_weight--; // count elements
88  sigx -= x; // update accumulators
89  sigy -= y;
90  sigxx -= x * x;
91  sigxy -= x * y;
92  sigyy -= y * y;
93 }
94 
95 
96 /**********************************************************************
97  * LLSQ::m
98  *
99  * Return the gradient of the line fit.
100  **********************************************************************/
101 
102 double LLSQ::m() const { // get gradient
103  double covar = covariance();
104  double x_var = x_variance();
105  if (x_var != 0.0)
106  return covar / x_var;
107  else
108  return 0.0; // too little
109 }
110 
111 
112 /**********************************************************************
113  * LLSQ::c
114  *
115  * Return the constant of the line fit.
116  **********************************************************************/
117 
118 double LLSQ::c(double m) const { // get constant
119  if (total_weight > 0.0)
120  return (sigy - m * sigx) / total_weight;
121  else
122  return 0; // too little
123 }
124 
125 
126 /**********************************************************************
127  * LLSQ::rms
128  *
129  * Return the rms error of the fit.
130  **********************************************************************/
131 
132 double LLSQ::rms(double m, double c) const { // get error
133  double error; // total error
134 
135  if (total_weight > 0) {
136  error = sigyy + m * (m * sigxx + 2 * (c * sigx - sigxy)) + c *
137  (total_weight * c - 2 * sigy);
138  if (error >= 0)
139  error = sqrt(error / total_weight); // sqrt of mean
140  else
141  error = 0;
142  } else {
143  error = 0; // too little
144  }
145  return error;
146 }
147 
148 
149 /**********************************************************************
150  * LLSQ::pearson
151  *
152  * Return the pearson product moment correlation coefficient.
153  **********************************************************************/
154 
155 double LLSQ::pearson() const { // get correlation
156  double r = 0.0; // Correlation is 0 if insufficent data.
157 
158  double covar = covariance();
159  if (covar != 0.0) {
160  double var_product = x_variance() * y_variance();
161  if (var_product > 0.0)
162  r = covar / sqrt(var_product);
163  }
164  return r;
165 }
166 
167 // Returns the x,y means as an FCOORD.
169  if (total_weight > 0.0) {
170  return FCOORD(sigx / total_weight, sigy / total_weight);
171  } else {
172  return FCOORD(0.0f, 0.0f);
173  }
174 }
175 
176 // Returns the direction of the fitted line as a unit vector, using the
177 // least mean squared perpendicular distance. The line runs through the
178 // mean_point, i.e. a point p on the line is given by:
179 // p = mean_point() + lambda * vector_fit() for some real number lambda.
180 // Note that the result (0<=x<=1, -1<=y<=-1) is directionally ambiguous
181 // and may be negated without changing its meaning.
183  double x_var = x_variance();
184  double y_var = y_variance();
185  double covar = covariance();
186  FCOORD result;
187  if (x_var >= y_var) {
188  if (x_var == 0.0)
189  return FCOORD(0.0f, 0.0f);
190  result.set_x(x_var / sqrt(x_var * x_var + covar * covar));
191  result.set_y(sqrt(1.0 - result.x() * result.x()));
192  } else {
193  result.set_y(y_var / sqrt(y_var * y_var + covar * covar));
194  result.set_x(sqrt(1.0 - result.y() * result.y()));
195  }
196  if (covar < 0.0)
197  result.set_y(-result.y());
198  return result;
199 }