1
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
2
<html xmlns="http://www.w3.org/1999/xhtml">
4
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
5
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
6
<title>tesseract: /usr/local/google/home/jbreiden/tesseract-ocr-read-only/ccstruct/linlsq.cpp Source File</title>
8
<link href="tabs.css" rel="stylesheet" type="text/css"/>
9
<link href="doxygen.css" rel="stylesheet" type="text/css" />
10
<link href="navtree.css" rel="stylesheet" type="text/css"/>
11
<script type="text/javascript" src="jquery.js"></script>
12
<script type="text/javascript" src="resize.js"></script>
13
<script type="text/javascript" src="navtree.js"></script>
14
<script type="text/javascript">
15
$(document).ready(initResizable);
17
<link href="search/search.css" rel="stylesheet" type="text/css"/>
18
<script type="text/javascript" src="search/search.js"></script>
19
<script type="text/javascript">
20
$(document).ready(function() { searchBox.OnSelectItem(0); });
25
<div id="top"><!-- do not remove this div! -->
29
<table cellspacing="0" cellpadding="0">
31
<tr style="height: 56px;">
34
<td style="padding-left: 0.5em;">
35
<div id="projectname">tesseract
36
 <span id="projectnumber">3.03</span>
48
<!-- Generated by Doxygen 1.7.6.1 -->
49
<script type="text/javascript">
50
var searchBox = new SearchBox("searchBox", "search",false,'Search');
52
<div id="navrow1" class="tabs">
54
<li><a href="index.html"><span>Main Page</span></a></li>
55
<li><a href="pages.html"><span>Related Pages</span></a></li>
56
<li><a href="modules.html"><span>Modules</span></a></li>
57
<li><a href="namespaces.html"><span>Namespaces</span></a></li>
58
<li><a href="annotated.html"><span>Classes</span></a></li>
59
<li class="current"><a href="files.html"><span>Files</span></a></li>
61
<div id="MSearchBox" class="MSearchBoxInactive">
63
<img id="MSearchSelect" src="search/mag_sel.png"
64
onmouseover="return searchBox.OnSearchSelectShow()"
65
onmouseout="return searchBox.OnSearchSelectHide()"
67
<input type="text" id="MSearchField" value="Search" accesskey="S"
68
onfocus="searchBox.OnSearchFieldFocus(true)"
69
onblur="searchBox.OnSearchFieldFocus(false)"
70
onkeyup="searchBox.OnSearchFieldChange(event)"/>
71
</span><span class="right">
72
<a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
78
<div id="navrow2" class="tabs2">
80
<li><a href="files.html"><span>File List</span></a></li>
81
<li><a href="globals.html"><span>File Members</span></a></li>
85
<div id="side-nav" class="ui-resizable side-nav-resizable">
87
<div id="nav-tree-contents">
90
<div id="splitbar" style="-moz-user-select:none;"
91
class="ui-resizable-handle">
94
<script type="text/javascript">
95
initNavTree('a00752.html','');
97
<div id="doc-content">
99
<div class="headertitle">
100
<div class="title">/usr/local/google/home/jbreiden/tesseract-ocr-read-only/ccstruct/linlsq.cpp</div> </div>
102
<div class="contents">
103
<a href="a00752.html">Go to the documentation of this file.</a><div class="fragment"><pre class="fragment"><a name="l00001"></a>00001 <span class="comment">/**********************************************************************</span>
104
<a name="l00002"></a>00002 <span class="comment"> * File: linlsq.cpp (Formerly llsq.c)</span>
105
<a name="l00003"></a>00003 <span class="comment"> * Description: Linear Least squares fitting code.</span>
106
<a name="l00004"></a>00004 <span class="comment"> * Author: Ray Smith</span>
107
<a name="l00005"></a>00005 <span class="comment"> * Created: Thu Sep 12 08:44:51 BST 1991</span>
108
<a name="l00006"></a>00006 <span class="comment"> *</span>
109
<a name="l00007"></a>00007 <span class="comment"> * (C) Copyright 1991, Hewlett-Packard Ltd.</span>
110
<a name="l00008"></a>00008 <span class="comment"> ** Licensed under the Apache License, Version 2.0 (the "License");</span>
111
<a name="l00009"></a>00009 <span class="comment"> ** you may not use this file except in compliance with the License.</span>
112
<a name="l00010"></a>00010 <span class="comment"> ** You may obtain a copy of the License at</span>
113
<a name="l00011"></a>00011 <span class="comment"> ** http://www.apache.org/licenses/LICENSE-2.0</span>
114
<a name="l00012"></a>00012 <span class="comment"> ** Unless required by applicable law or agreed to in writing, software</span>
115
<a name="l00013"></a>00013 <span class="comment"> ** distributed under the License is distributed on an "AS IS" BASIS,</span>
116
<a name="l00014"></a>00014 <span class="comment"> ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.</span>
117
<a name="l00015"></a>00015 <span class="comment"> ** See the License for the specific language governing permissions and</span>
118
<a name="l00016"></a>00016 <span class="comment"> ** limitations under the License.</span>
119
<a name="l00017"></a>00017 <span class="comment"> *</span>
120
<a name="l00018"></a>00018 <span class="comment"> **********************************************************************/</span>
121
<a name="l00019"></a>00019
122
<a name="l00020"></a>00020 <span class="preprocessor">#include <stdio.h></span>
123
<a name="l00021"></a>00021 <span class="preprocessor">#include <math.h></span>
124
<a name="l00022"></a>00022 <span class="preprocessor">#include "<a class="code" href="a00823.html">errcode.h</a>"</span>
125
<a name="l00023"></a>00023 <span class="preprocessor">#include "<a class="code" href="a00753.html">linlsq.h</a>"</span>
126
<a name="l00024"></a>00024
127
<a name="l00025"></a><a class="code" href="a00752.html#a39d4e21e57d29d2cbd0031886ba5822a">00025</a> <span class="keyword">const</span> <a class="code" href="a00372.html">ERRCODE</a> <a class="code" href="a00752.html#a39d4e21e57d29d2cbd0031886ba5822a">EMPTY_LLSQ</a> = <span class="stringliteral">"Can't delete from an empty LLSQ"</span>;
128
<a name="l00026"></a>00026
129
<a name="l00027"></a>00027 <span class="comment">/**********************************************************************</span>
130
<a name="l00028"></a>00028 <span class="comment"> * LLSQ::clear</span>
131
<a name="l00029"></a>00029 <span class="comment"> *</span>
132
<a name="l00030"></a>00030 <span class="comment"> * Function to initialize a LLSQ.</span>
133
<a name="l00031"></a>00031 <span class="comment"> **********************************************************************/</span>
134
<a name="l00032"></a>00032
135
<a name="l00033"></a><a class="code" href="a00454.html#a9f7b92005193802bc03539a6ae48a27a">00033</a> <span class="keywordtype">void</span> <a class="code" href="a00454.html#a9f7b92005193802bc03539a6ae48a27a">LLSQ::clear</a>() { <span class="comment">// initialize</span>
136
<a name="l00034"></a>00034 total_weight = 0.0; <span class="comment">// no elements</span>
137
<a name="l00035"></a>00035 sigx = 0.0; <span class="comment">// update accumulators</span>
138
<a name="l00036"></a>00036 sigy = 0.0;
139
<a name="l00037"></a>00037 sigxx = 0.0;
140
<a name="l00038"></a>00038 sigxy = 0.0;
141
<a name="l00039"></a>00039 sigyy = 0.0;
142
<a name="l00040"></a>00040 }
143
<a name="l00041"></a>00041
144
<a name="l00042"></a>00042
145
<a name="l00043"></a>00043 <span class="comment">/**********************************************************************</span>
146
<a name="l00044"></a>00044 <span class="comment"> * LLSQ::add</span>
147
<a name="l00045"></a>00045 <span class="comment"> *</span>
148
<a name="l00046"></a>00046 <span class="comment"> * Add an element to the accumulator.</span>
149
<a name="l00047"></a>00047 <span class="comment"> **********************************************************************/</span>
150
<a name="l00048"></a>00048
151
<a name="l00049"></a><a class="code" href="a00454.html#aa439294be0539828cf042f8941125b3c">00049</a> <span class="keywordtype">void</span> <a class="code" href="a00454.html#aa439294be0539828cf042f8941125b3c">LLSQ::add</a>(<span class="keywordtype">double</span> x, <span class="keywordtype">double</span> y) { <span class="comment">// add an element</span>
152
<a name="l00050"></a>00050 total_weight++; <span class="comment">// count elements</span>
153
<a name="l00051"></a>00051 sigx += x; <span class="comment">// update accumulators</span>
154
<a name="l00052"></a>00052 sigy += y;
155
<a name="l00053"></a>00053 sigxx += x * x;
156
<a name="l00054"></a>00054 sigxy += x * y;
157
<a name="l00055"></a>00055 sigyy += y * y;
158
<a name="l00056"></a>00056 }
159
<a name="l00057"></a>00057 <span class="comment">// Adds an element with a specified weight.</span>
160
<a name="l00058"></a><a class="code" href="a00454.html#a9a52d110b019972ea2d6009352c4c592">00058</a> <span class="keywordtype">void</span> <a class="code" href="a00454.html#aa439294be0539828cf042f8941125b3c">LLSQ::add</a>(<span class="keywordtype">double</span> x, <span class="keywordtype">double</span> y, <span class="keywordtype">double</span> weight) {
161
<a name="l00059"></a>00059 total_weight += weight;
162
<a name="l00060"></a>00060 sigx += x * weight; <span class="comment">// update accumulators</span>
163
<a name="l00061"></a>00061 sigy += y * weight;
164
<a name="l00062"></a>00062 sigxx += x * x * weight;
165
<a name="l00063"></a>00063 sigxy += x * y * weight;
166
<a name="l00064"></a>00064 sigyy += y * y * weight;
167
<a name="l00065"></a>00065 }
168
<a name="l00066"></a>00066 <span class="comment">// Adds a whole LLSQ.</span>
169
<a name="l00067"></a><a class="code" href="a00454.html#ae06248347a9cea5045ff8d5b92d76962">00067</a> <span class="keywordtype">void</span> <a class="code" href="a00454.html#aa439294be0539828cf042f8941125b3c">LLSQ::add</a>(<span class="keyword">const</span> <a class="code" href="a00454.html">LLSQ</a>& other) {
170
<a name="l00068"></a>00068 total_weight += other.total_weight;
171
<a name="l00069"></a>00069 sigx += other.sigx; <span class="comment">// update accumulators</span>
172
<a name="l00070"></a>00070 sigy += other.sigy;
173
<a name="l00071"></a>00071 sigxx += other.sigxx;
174
<a name="l00072"></a>00072 sigxy += other.sigxy;
175
<a name="l00073"></a>00073 sigyy += other.sigyy;
176
<a name="l00074"></a>00074 }
177
<a name="l00075"></a>00075
178
<a name="l00076"></a>00076
179
<a name="l00077"></a>00077 <span class="comment">/**********************************************************************</span>
180
<a name="l00078"></a>00078 <span class="comment"> * LLSQ::remove</span>
181
<a name="l00079"></a>00079 <span class="comment"> *</span>
182
<a name="l00080"></a>00080 <span class="comment"> * Delete an element from the acculuator.</span>
183
<a name="l00081"></a>00081 <span class="comment"> **********************************************************************/</span>
184
<a name="l00082"></a>00082
185
<a name="l00083"></a><a class="code" href="a00454.html#a2fc29873c4ec3ae881a26c54813d36ff">00083</a> <span class="keywordtype">void</span> <a class="code" href="a00454.html#a2fc29873c4ec3ae881a26c54813d36ff">LLSQ::remove</a>(<span class="keywordtype">double</span> x, <span class="keywordtype">double</span> y) { <span class="comment">// delete an element</span>
186
<a name="l00084"></a>00084 <span class="keywordflow">if</span> (total_weight <= 0.0) <span class="comment">// illegal</span>
187
<a name="l00085"></a>00085 EMPTY_LLSQ.<a class="code" href="a00372.html#abc4cc64aad0ddb74c02523b69c95a1c1">error</a>(<span class="stringliteral">"LLSQ::remove"</span>, <a class="code" href="a00823.html#a8a3a2e51383909b357fcdf8eda803637a781ad2788df9e25c59a70894c7832096">ABORT</a>, NULL);
188
<a name="l00086"></a>00086 total_weight--; <span class="comment">// count elements</span>
189
<a name="l00087"></a>00087 sigx -= x; <span class="comment">// update accumulators</span>
190
<a name="l00088"></a>00088 sigy -= y;
191
<a name="l00089"></a>00089 sigxx -= x * x;
192
<a name="l00090"></a>00090 sigxy -= x * y;
193
<a name="l00091"></a>00091 sigyy -= y * y;
194
<a name="l00092"></a>00092 }
195
<a name="l00093"></a>00093
196
<a name="l00094"></a>00094
197
<a name="l00095"></a>00095 <span class="comment">/**********************************************************************</span>
198
<a name="l00096"></a>00096 <span class="comment"> * LLSQ::m</span>
199
<a name="l00097"></a>00097 <span class="comment"> *</span>
200
<a name="l00098"></a>00098 <span class="comment"> * Return the gradient of the line fit.</span>
201
<a name="l00099"></a>00099 <span class="comment"> **********************************************************************/</span>
202
<a name="l00100"></a>00100
203
<a name="l00101"></a><a class="code" href="a00454.html#ad3aed3a5a46cf180dd80893c1034aa9e">00101</a> <span class="keywordtype">double</span> <a class="code" href="a00454.html#ad3aed3a5a46cf180dd80893c1034aa9e">LLSQ::m</a>()<span class="keyword"> const </span>{ <span class="comment">// get gradient</span>
204
<a name="l00102"></a>00102 <span class="keywordtype">double</span> covar = <a class="code" href="a00454.html#a36d6bf1e469425fc4f9b97521aee481e">covariance</a>();
205
<a name="l00103"></a>00103 <span class="keywordtype">double</span> x_var = <a class="code" href="a00454.html#a27994b02b1ea16d8e3a6c951c90205dd">x_variance</a>();
206
<a name="l00104"></a>00104 <span class="keywordflow">if</span> (x_var != 0.0)
207
<a name="l00105"></a>00105 <span class="keywordflow">return</span> covar / x_var;
208
<a name="l00106"></a>00106 <span class="keywordflow">else</span>
209
<a name="l00107"></a>00107 <span class="keywordflow">return</span> 0.0; <span class="comment">// too little</span>
210
<a name="l00108"></a>00108 }
211
<a name="l00109"></a>00109
212
<a name="l00110"></a>00110
213
<a name="l00111"></a>00111 <span class="comment">/**********************************************************************</span>
214
<a name="l00112"></a>00112 <span class="comment"> * LLSQ::c</span>
215
<a name="l00113"></a>00113 <span class="comment"> *</span>
216
<a name="l00114"></a>00114 <span class="comment"> * Return the constant of the line fit.</span>
217
<a name="l00115"></a>00115 <span class="comment"> **********************************************************************/</span>
218
<a name="l00116"></a>00116
219
<a name="l00117"></a><a class="code" href="a00454.html#af66af81c0276216acfa659178189019f">00117</a> <span class="keywordtype">double</span> <a class="code" href="a00454.html#af66af81c0276216acfa659178189019f">LLSQ::c</a>(<span class="keywordtype">double</span> m)<span class="keyword"> const </span>{ <span class="comment">// get constant</span>
220
<a name="l00118"></a>00118 <span class="keywordflow">if</span> (total_weight > 0.0)
221
<a name="l00119"></a>00119 <span class="keywordflow">return</span> (sigy - m * sigx) / total_weight;
222
<a name="l00120"></a>00120 <span class="keywordflow">else</span>
223
<a name="l00121"></a>00121 <span class="keywordflow">return</span> 0; <span class="comment">// too little</span>
224
<a name="l00122"></a>00122 }
225
<a name="l00123"></a>00123
226
<a name="l00124"></a>00124
227
<a name="l00125"></a>00125 <span class="comment">/**********************************************************************</span>
228
<a name="l00126"></a>00126 <span class="comment"> * LLSQ::rms</span>
229
<a name="l00127"></a>00127 <span class="comment"> *</span>
230
<a name="l00128"></a>00128 <span class="comment"> * Return the rms error of the fit.</span>
231
<a name="l00129"></a>00129 <span class="comment"> **********************************************************************/</span>
232
<a name="l00130"></a>00130
233
<a name="l00131"></a><a class="code" href="a00454.html#adb169c39429f6f206130e463c557af67">00131</a> <span class="keywordtype">double</span> <a class="code" href="a00454.html#adb169c39429f6f206130e463c557af67">LLSQ::rms</a>(<span class="keywordtype">double</span> m, <span class="keywordtype">double</span> c)<span class="keyword"> const </span>{ <span class="comment">// get error</span>
234
<a name="l00132"></a>00132 <span class="keywordtype">double</span> error; <span class="comment">// total error</span>
235
<a name="l00133"></a>00133
236
<a name="l00134"></a>00134 <span class="keywordflow">if</span> (total_weight > 0) {
237
<a name="l00135"></a>00135 error = sigyy + m * (m * sigxx + 2 * (c * sigx - sigxy)) + c *
238
<a name="l00136"></a>00136 (total_weight * c - 2 * sigy);
239
<a name="l00137"></a>00137 <span class="keywordflow">if</span> (error >= 0)
240
<a name="l00138"></a>00138 error = sqrt(error / total_weight); <span class="comment">// sqrt of mean</span>
241
<a name="l00139"></a>00139 <span class="keywordflow">else</span>
242
<a name="l00140"></a>00140 error = 0;
243
<a name="l00141"></a>00141 } <span class="keywordflow">else</span> {
244
<a name="l00142"></a>00142 error = 0; <span class="comment">// too little</span>
245
<a name="l00143"></a>00143 }
246
<a name="l00144"></a>00144 <span class="keywordflow">return</span> error;
247
<a name="l00145"></a>00145 }
248
<a name="l00146"></a>00146
249
<a name="l00147"></a>00147
250
<a name="l00148"></a>00148 <span class="comment">/**********************************************************************</span>
251
<a name="l00149"></a>00149 <span class="comment"> * LLSQ::pearson</span>
252
<a name="l00150"></a>00150 <span class="comment"> *</span>
253
<a name="l00151"></a>00151 <span class="comment"> * Return the pearson product moment correlation coefficient.</span>
254
<a name="l00152"></a>00152 <span class="comment"> **********************************************************************/</span>
255
<a name="l00153"></a>00153
256
<a name="l00154"></a><a class="code" href="a00454.html#a8537176180a1cd3ff9c3b1f97d1ffb78">00154</a> <span class="keywordtype">double</span> <a class="code" href="a00454.html#a8537176180a1cd3ff9c3b1f97d1ffb78">LLSQ::pearson</a>()<span class="keyword"> const </span>{ <span class="comment">// get correlation</span>
257
<a name="l00155"></a>00155 <span class="keywordtype">double</span> r = 0.0; <span class="comment">// Correlation is 0 if insufficent data.</span>
258
<a name="l00156"></a>00156
259
<a name="l00157"></a>00157 <span class="keywordtype">double</span> covar = <a class="code" href="a00454.html#a36d6bf1e469425fc4f9b97521aee481e">covariance</a>();
260
<a name="l00158"></a>00158 <span class="keywordflow">if</span> (covar != 0.0) {
261
<a name="l00159"></a>00159 <span class="keywordtype">double</span> var_product = <a class="code" href="a00454.html#a27994b02b1ea16d8e3a6c951c90205dd">x_variance</a>() * <a class="code" href="a00454.html#a3e386f3db4343cca32ac4a4f028c6ea4">y_variance</a>();
262
<a name="l00160"></a>00160 <span class="keywordflow">if</span> (var_product > 0.0)
263
<a name="l00161"></a>00161 r = covar / sqrt(var_product);
264
<a name="l00162"></a>00162 }
265
<a name="l00163"></a>00163 <span class="keywordflow">return</span> r;
266
<a name="l00164"></a>00164 }
267
<a name="l00165"></a>00165
268
<a name="l00166"></a>00166 <span class="comment">// Returns the x,y means as an FCOORD.</span>
269
<a name="l00167"></a><a class="code" href="a00454.html#ac786685f62b271ca20aa719a8a03aa3d">00167</a> <a class="code" href="a00375.html">FCOORD</a> <a class="code" href="a00454.html#ac786685f62b271ca20aa719a8a03aa3d">LLSQ::mean_point</a>()<span class="keyword"> const </span>{
270
<a name="l00168"></a>00168 <span class="keywordflow">if</span> (total_weight > 0.0) {
271
<a name="l00169"></a>00169 <span class="keywordflow">return</span> <a class="code" href="a00375.html">FCOORD</a>(sigx / total_weight, sigy / total_weight);
272
<a name="l00170"></a>00170 } <span class="keywordflow">else</span> {
273
<a name="l00171"></a>00171 <span class="keywordflow">return</span> <a class="code" href="a00375.html">FCOORD</a>(0.0f, 0.0f);
274
<a name="l00172"></a>00172 }
275
<a name="l00173"></a>00173 }
276
<a name="l00174"></a>00174
277
<a name="l00175"></a>00175 <span class="comment">// Returns the sqrt of the mean squared error measured perpendicular from the</span>
278
<a name="l00176"></a>00176 <span class="comment">// line through mean_point() in the direction dir.</span>
279
<a name="l00177"></a>00177 <span class="comment">//</span>
280
<a name="l00178"></a>00178 <span class="comment">// Derivation:</span>
281
<a name="l00179"></a>00179 <span class="comment">// Lemma: Let v and x_i (i=1..N) be a k-dimensional vectors (1xk matrices).</span>
282
<a name="l00180"></a>00180 <span class="comment">// Let % be dot product and ' be transpose. Note that:</span>
283
<a name="l00181"></a>00181 <span class="comment">// Sum[i=1..N] (v % x_i)^2</span>
284
<a name="l00182"></a>00182 <span class="comment">// = v * [x_1' x_2' ... x_N'] * [x_1' x_2' .. x_N']' * v'</span>
285
<a name="l00183"></a>00183 <span class="comment">// If x_i have average 0 we have:</span>
286
<a name="l00184"></a>00184 <span class="comment">// = v * (N * COVARIANCE_MATRIX(X)) * v'</span>
287
<a name="l00185"></a>00185 <span class="comment">// Expanded for the case that k = 2, where we treat the dimensions</span>
288
<a name="l00186"></a>00186 <span class="comment">// as x_i and y_i, this is:</span>
289
<a name="l00187"></a>00187 <span class="comment">// = v * (N * [VAR(X), COV(X,Y); COV(X,Y) VAR(Y)]) * v'</span>
290
<a name="l00188"></a>00188 <span class="comment">// Now, we are trying to calculate the mean squared error, where v is</span>
291
<a name="l00189"></a>00189 <span class="comment">// perpendicular to our line of interest:</span>
292
<a name="l00190"></a>00190 <span class="comment">// Mean squared error</span>
293
<a name="l00191"></a>00191 <span class="comment">// = E [ (v % (x_i - x_avg))) ^2 ]</span>
294
<a name="l00192"></a>00192 <span class="comment">// = Sum (v % (x_i - x_avg))^2 / N</span>
295
<a name="l00193"></a>00193 <span class="comment">// = v * N * [VAR(X) COV(X,Y); COV(X,Y) VAR(Y)] / N * v'</span>
296
<a name="l00194"></a>00194 <span class="comment">// = v * [VAR(X) COV(X,Y); COV(X,Y) VAR(Y)] * v'</span>
297
<a name="l00195"></a>00195 <span class="comment">// = code below</span>
298
<a name="l00196"></a><a class="code" href="a00454.html#af047401ce847d873866f6018d15d4a95">00196</a> <span class="keywordtype">double</span> <a class="code" href="a00454.html#af047401ce847d873866f6018d15d4a95">LLSQ::rms_orth</a>(<span class="keyword">const</span> <a class="code" href="a00375.html">FCOORD</a> &dir)<span class="keyword"> const </span>{
299
<a name="l00197"></a>00197 <a class="code" href="a00375.html">FCOORD</a> v = !dir;
300
<a name="l00198"></a>00198 v.<a class="code" href="a00375.html#ae038bea84e90d493366d15c9578e4be9" title="Convert to unit vec.">normalise</a>();
301
<a name="l00199"></a>00199 <span class="keywordflow">return</span> sqrt(v.<a class="code" href="a00375.html#a9edadbfc357c730f62fa3c92677b0d58">x</a>() * v.<a class="code" href="a00375.html#a9edadbfc357c730f62fa3c92677b0d58">x</a>() * <a class="code" href="a00454.html#a27994b02b1ea16d8e3a6c951c90205dd">x_variance</a>() +
302
<a name="l00200"></a>00200 2 * v.<a class="code" href="a00375.html#a9edadbfc357c730f62fa3c92677b0d58">x</a>() * v.<a class="code" href="a00375.html#a397f7ce997b246df18e7bd5a20cc422e">y</a>() * <a class="code" href="a00454.html#a36d6bf1e469425fc4f9b97521aee481e">covariance</a>() +
303
<a name="l00201"></a>00201 v.<a class="code" href="a00375.html#a397f7ce997b246df18e7bd5a20cc422e">y</a>() * v.<a class="code" href="a00375.html#a397f7ce997b246df18e7bd5a20cc422e">y</a>() * <a class="code" href="a00454.html#a3e386f3db4343cca32ac4a4f028c6ea4">y_variance</a>());
304
<a name="l00202"></a>00202 }
305
<a name="l00203"></a>00203
306
<a name="l00204"></a>00204 <span class="comment">// Returns the direction of the fitted line as a unit vector, using the</span>
307
<a name="l00205"></a>00205 <span class="comment">// least mean squared perpendicular distance. The line runs through the</span>
308
<a name="l00206"></a>00206 <span class="comment">// mean_point, i.e. a point p on the line is given by:</span>
309
<a name="l00207"></a>00207 <span class="comment">// p = mean_point() + lambda * vector_fit() for some real number lambda.</span>
310
<a name="l00208"></a>00208 <span class="comment">// Note that the result (0<=x<=1, -1<=y<=-1) is directionally ambiguous</span>
311
<a name="l00209"></a>00209 <span class="comment">// and may be negated without changing its meaning.</span>
312
<a name="l00210"></a>00210 <span class="comment">// Fitting a line m + 𝜆v to a set of N points Pi = (xi, yi), where</span>
313
<a name="l00211"></a>00211 <span class="comment">// m is the mean point (𝝁, 𝝂) and</span>
314
<a name="l00212"></a>00212 <span class="comment">// v is the direction vector (cos𝜃, sin𝜃)</span>
315
<a name="l00213"></a>00213 <span class="comment">// The perpendicular distance of each Pi from the line is:</span>
316
<a name="l00214"></a>00214 <span class="comment">// (Pi - m) x v, where x is the scalar cross product.</span>
317
<a name="l00215"></a>00215 <span class="comment">// Total squared error is thus:</span>
318
<a name="l00216"></a>00216 <span class="comment">// E = ∑((xi - 𝝁)sin𝜃 - (yi - 𝝂)cos𝜃)²</span>
319
<a name="l00217"></a>00217 <span class="comment">// = ∑(xi - 𝝁)²sin²𝜃 - 2∑(xi - 𝝁)(yi - 𝝂)sin𝜃 cos𝜃 + ∑(yi - 𝝂)²cos²𝜃</span>
320
<a name="l00218"></a>00218 <span class="comment">// = NVar(xi)sin²𝜃 - 2NCovar(xi, yi)sin𝜃 cos𝜃 + NVar(yi)cos²𝜃 (Eq 1)</span>
321
<a name="l00219"></a>00219 <span class="comment">// where Var(xi) is the variance of xi,</span>
322
<a name="l00220"></a>00220 <span class="comment">// and Covar(xi, yi) is the covariance of xi, yi.</span>
323
<a name="l00221"></a>00221 <span class="comment">// Taking the derivative wrt 𝜃 and setting to 0 to obtain the min/max:</span>
324
<a name="l00222"></a>00222 <span class="comment">// 0 = 2NVar(xi)sin𝜃 cos𝜃 -2NCovar(xi, yi)(cos²𝜃 - sin²𝜃) -2NVar(yi)sin𝜃 cos𝜃</span>
325
<a name="l00223"></a>00223 <span class="comment">// => Covar(xi, yi)(cos²𝜃 - sin²𝜃) = (Var(xi) - Var(yi))sin𝜃 cos𝜃</span>
326
<a name="l00224"></a>00224 <span class="comment">// Using double angles:</span>
327
<a name="l00225"></a>00225 <span class="comment">// 2Covar(xi, yi)cos2𝜃 = (Var(xi) - Var(yi))sin2𝜃 (Eq 2)</span>
328
<a name="l00226"></a>00226 <span class="comment">// So 𝜃 = 0.5 atan2(2Covar(xi, yi), Var(xi) - Var(yi)) (Eq 3)</span>
329
<a name="l00227"></a>00227
330
<a name="l00228"></a>00228 <span class="comment">// Because it involves 2𝜃 , Eq 2 has 2 solutions 90 degrees apart, but which</span>
331
<a name="l00229"></a>00229 <span class="comment">// is the min and which is the max? From Eq1:</span>
332
<a name="l00230"></a>00230 <span class="comment">// E/N = Var(xi)sin²𝜃 - 2Covar(xi, yi)sin𝜃 cos𝜃 + Var(yi)cos²𝜃</span>
333
<a name="l00231"></a>00231 <span class="comment">// and 90 degrees away, using sin/cos equivalences:</span>
334
<a name="l00232"></a>00232 <span class="comment">// E'/N = Var(xi)cos²𝜃 + 2Covar(xi, yi)sin𝜃 cos𝜃 + Var(yi)sin²𝜃</span>
335
<a name="l00233"></a>00233 <span class="comment">// The second error is smaller (making it the minimum) iff</span>
336
<a name="l00234"></a>00234 <span class="comment">// E'/N < E/N ie:</span>
337
<a name="l00235"></a>00235 <span class="comment">// (Var(xi) - Var(yi))(cos²𝜃 - sin²𝜃) < -4Covar(xi, yi)sin𝜃 cos𝜃</span>
338
<a name="l00236"></a>00236 <span class="comment">// Using double angles:</span>
339
<a name="l00237"></a>00237 <span class="comment">// (Var(xi) - Var(yi))cos2𝜃 < -2Covar(xi, yi)sin2𝜃 (InEq 1)</span>
340
<a name="l00238"></a>00238 <span class="comment">// But atan2(2Covar(xi, yi), Var(xi) - Var(yi)) picks 2𝜃 such that:</span>
341
<a name="l00239"></a>00239 <span class="comment">// sgn(cos2𝜃) = sgn(Var(xi) - Var(yi)) and sgn(sin2𝜃) = sgn(Covar(xi, yi))</span>
342
<a name="l00240"></a>00240 <span class="comment">// so InEq1 can *never* be true, making the atan2 result *always* the min!</span>
343
<a name="l00241"></a>00241 <span class="comment">// In the degenerate case, where Covar(xi, yi) = 0 AND Var(xi) = Var(yi),</span>
344
<a name="l00242"></a>00242 <span class="comment">// the 2 solutions have equal error and the inequality is still false.</span>
345
<a name="l00243"></a>00243 <span class="comment">// Therefore the solution really is as trivial as Eq 3.</span>
346
<a name="l00244"></a>00244
347
<a name="l00245"></a>00245 <span class="comment">// This is equivalent to returning the Principal Component in PCA, or the</span>
348
<a name="l00246"></a>00246 <span class="comment">// eigenvector corresponding to the largest eigenvalue in the covariance</span>
349
<a name="l00247"></a>00247 <span class="comment">// matrix. However, atan2 is much simpler! The one reference I found that</span>
350
<a name="l00248"></a>00248 <span class="comment">// uses this formula is http://web.mit.edu/18.06/www/Essays/tlsfit.pdf but</span>
351
<a name="l00249"></a>00249 <span class="comment">// that is still a much more complex derivation. It seems Pearson had already</span>
352
<a name="l00250"></a>00250 <span class="comment">// found this simple solution in 1901.</span>
353
<a name="l00251"></a>00251 <span class="comment">// http://books.google.com/books?id=WXwvAQAAIAAJ&pg=PA559</span>
354
<a name="l00252"></a><a class="code" href="a00454.html#a1e15590443e86960c79ae5785f850d6d">00252</a> <a class="code" href="a00375.html">FCOORD</a> <a class="code" href="a00454.html#a1e15590443e86960c79ae5785f850d6d">LLSQ::vector_fit</a>()<span class="keyword"> const </span>{
355
<a name="l00253"></a>00253 <span class="keywordtype">double</span> x_var = <a class="code" href="a00454.html#a27994b02b1ea16d8e3a6c951c90205dd">x_variance</a>();
356
<a name="l00254"></a>00254 <span class="keywordtype">double</span> y_var = <a class="code" href="a00454.html#a3e386f3db4343cca32ac4a4f028c6ea4">y_variance</a>();
357
<a name="l00255"></a>00255 <span class="keywordtype">double</span> covar = <a class="code" href="a00454.html#a36d6bf1e469425fc4f9b97521aee481e">covariance</a>();
358
<a name="l00256"></a>00256 <span class="keywordtype">double</span> theta = 0.5 * atan2(2.0 * covar, x_var - y_var);
359
<a name="l00257"></a>00257 <a class="code" href="a00375.html">FCOORD</a> result(cos(theta), sin(theta));
360
<a name="l00258"></a>00258 <span class="keywordflow">return</span> result;
361
<a name="l00259"></a>00259 }
362
</pre></div></div><!-- contents -->
364
<!-- window showing the filter options -->
365
<div id="MSearchSelectWindow"
366
onmouseover="return searchBox.OnSearchSelectShow()"
367
onmouseout="return searchBox.OnSearchSelectHide()"
368
onkeydown="return searchBox.OnSearchSelectKey(event)">
369
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark"> </span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark"> </span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark"> </span>Namespaces</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark"> </span>Files</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark"> </span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(5)"><span class="SelectionMark"> </span>Variables</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(6)"><span class="SelectionMark"> </span>Typedefs</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(7)"><span class="SelectionMark"> </span>Enumerations</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(8)"><span class="SelectionMark"> </span>Enumerator</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(9)"><span class="SelectionMark"> </span>Friends</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(10)"><span class="SelectionMark"> </span>Defines</a></div>
371
<!-- iframe showing the search results (closed by default) -->
372
<div id="MSearchResultsWindow">
373
<iframe src="javascript:void(0)" frameborder="0"
374
name="MSearchResults" id="MSearchResults">
378
<div id="nav-path" class="navpath">
380
<li class="navelem"><a class="el" href="a00752.html">linlsq.cpp</a> </li>
382
<li class="footer">Generated on Mon Feb 3 2014 10:59:07 for tesseract by
383
<a href="http://www.doxygen.org/index.html">
384
<img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.7.6.1 </li>
added
removed
doc/html/a00752_source.html
