Tesseract
3.02
|
00001 /********************************************************************** 00002 * File: linlsq.cpp (Formerly llsq.c) 00003 * Description: Linear Least squares fitting code. 00004 * Author: Ray Smith 00005 * Created: Thu Sep 12 08:44:51 BST 1991 00006 * 00007 * (C) Copyright 1991, Hewlett-Packard Ltd. 00008 ** Licensed under the Apache License, Version 2.0 (the "License"); 00009 ** you may not use this file except in compliance with the License. 00010 ** You may obtain a copy of the License at 00011 ** http://www.apache.org/licenses/LICENSE-2.0 00012 ** Unless required by applicable law or agreed to in writing, software 00013 ** distributed under the License is distributed on an "AS IS" BASIS, 00014 ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 00015 ** See the License for the specific language governing permissions and 00016 ** limitations under the License. 00017 * 00018 **********************************************************************/ 00019 00020 #include "mfcpch.h" // Must be first include for windows. 00021 #include <stdio.h> 00022 #include <math.h> 00023 #include "errcode.h" 00024 #include "linlsq.h" 00025 00026 const ERRCODE EMPTY_LLSQ = "Can't delete from an empty LLSQ"; 00027 00028 /********************************************************************** 00029 * LLSQ::clear 00030 * 00031 * Function to initialize a LLSQ. 00032 **********************************************************************/ 00033 00034 void LLSQ::clear() { // initialize 00035 total_weight = 0.0; // no elements 00036 sigx = 0.0; // update accumulators 00037 sigy = 0.0; 00038 sigxx = 0.0; 00039 sigxy = 0.0; 00040 sigyy = 0.0; 00041 } 00042 00043 00044 /********************************************************************** 00045 * LLSQ::add 00046 * 00047 * Add an element to the accumulator. 00048 **********************************************************************/ 00049 00050 void LLSQ::add(double x, double y) { // add an element 00051 total_weight++; // count elements 00052 sigx += x; // update accumulators 00053 sigy += y; 00054 sigxx += x * x; 00055 sigxy += x * y; 00056 sigyy += y * y; 00057 } 00058 // Adds an element with a specified weight. 00059 void LLSQ::add(double x, double y, double weight) { 00060 total_weight += weight; 00061 sigx += x * weight; // update accumulators 00062 sigy += y * weight; 00063 sigxx += x * x * weight; 00064 sigxy += x * y * weight; 00065 sigyy += y * y * weight; 00066 } 00067 // Adds a whole LLSQ. 00068 void LLSQ::add(const LLSQ& other) { 00069 total_weight += other.total_weight; 00070 sigx += other.sigx; // update accumulators 00071 sigy += other.sigy; 00072 sigxx += other.sigxx; 00073 sigxy += other.sigxy; 00074 sigyy += other.sigyy; 00075 } 00076 00077 00078 /********************************************************************** 00079 * LLSQ::remove 00080 * 00081 * Delete an element from the acculuator. 00082 **********************************************************************/ 00083 00084 void LLSQ::remove(double x, double y) { // delete an element 00085 if (total_weight <= 0.0) // illegal 00086 EMPTY_LLSQ.error("LLSQ::remove", ABORT, NULL); 00087 total_weight--; // count elements 00088 sigx -= x; // update accumulators 00089 sigy -= y; 00090 sigxx -= x * x; 00091 sigxy -= x * y; 00092 sigyy -= y * y; 00093 } 00094 00095 00096 /********************************************************************** 00097 * LLSQ::m 00098 * 00099 * Return the gradient of the line fit. 00100 **********************************************************************/ 00101 00102 double LLSQ::m() const { // get gradient 00103 double covar = covariance(); 00104 double x_var = x_variance(); 00105 if (x_var != 0.0) 00106 return covar / x_var; 00107 else 00108 return 0.0; // too little 00109 } 00110 00111 00112 /********************************************************************** 00113 * LLSQ::c 00114 * 00115 * Return the constant of the line fit. 00116 **********************************************************************/ 00117 00118 double LLSQ::c(double m) const { // get constant 00119 if (total_weight > 0.0) 00120 return (sigy - m * sigx) / total_weight; 00121 else 00122 return 0; // too little 00123 } 00124 00125 00126 /********************************************************************** 00127 * LLSQ::rms 00128 * 00129 * Return the rms error of the fit. 00130 **********************************************************************/ 00131 00132 double LLSQ::rms(double m, double c) const { // get error 00133 double error; // total error 00134 00135 if (total_weight > 0) { 00136 error = sigyy + m * (m * sigxx + 2 * (c * sigx - sigxy)) + c * 00137 (total_weight * c - 2 * sigy); 00138 if (error >= 0) 00139 error = sqrt(error / total_weight); // sqrt of mean 00140 else 00141 error = 0; 00142 } else { 00143 error = 0; // too little 00144 } 00145 return error; 00146 } 00147 00148 00149 /********************************************************************** 00150 * LLSQ::pearson 00151 * 00152 * Return the pearson product moment correlation coefficient. 00153 **********************************************************************/ 00154 00155 double LLSQ::pearson() const { // get correlation 00156 double r = 0.0; // Correlation is 0 if insufficent data. 00157 00158 double covar = covariance(); 00159 if (covar != 0.0) { 00160 double var_product = x_variance() * y_variance(); 00161 if (var_product > 0.0) 00162 r = covar / sqrt(var_product); 00163 } 00164 return r; 00165 } 00166 00167 // Returns the x,y means as an FCOORD. 00168 FCOORD LLSQ::mean_point() const { 00169 if (total_weight > 0.0) { 00170 return FCOORD(sigx / total_weight, sigy / total_weight); 00171 } else { 00172 return FCOORD(0.0f, 0.0f); 00173 } 00174 } 00175 00176 // Returns the direction of the fitted line as a unit vector, using the 00177 // least mean squared perpendicular distance. The line runs through the 00178 // mean_point, i.e. a point p on the line is given by: 00179 // p = mean_point() + lambda * vector_fit() for some real number lambda. 00180 // Note that the result (0<=x<=1, -1<=y<=-1) is directionally ambiguous 00181 // and may be negated without changing its meaning. 00182 FCOORD LLSQ::vector_fit() const { 00183 double x_var = x_variance(); 00184 double y_var = y_variance(); 00185 double covar = covariance(); 00186 FCOORD result; 00187 if (x_var >= y_var) { 00188 if (x_var == 0.0) 00189 return FCOORD(0.0f, 0.0f); 00190 result.set_x(x_var / sqrt(x_var * x_var + covar * covar)); 00191 result.set_y(sqrt(1.0 - result.x() * result.x())); 00192 } else { 00193 result.set_y(y_var / sqrt(y_var * y_var + covar * covar)); 00194 result.set_x(sqrt(1.0 - result.y() * result.y())); 00195 } 00196 if (covar < 0.0) 00197 result.set_y(-result.y()); 00198 return result; 00199 }