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
clear()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
add(double x,double y)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
remove(double x,double y)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
fit(int degree)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 }
148