Engauge Digitizer 2
Loading...
Searching...
No Matches
TestSpline.cpp
Go to the documentation of this file.
1#include <iostream>
2#include "Logger.h"
3#include "MainWindow.h"
4#include <map>
5#include <qmath.h>
6#include <QtTest/QtTest>
7#include "Spline.h"
8#include "SplinePair.h"
9#include <sstream>
10#include "Test/TestSpline.h"
11
12QTEST_MAIN (TestSpline)
13
14using namespace std;
15
16const QString WEBPAGE ("https://tools.timodenk.com/cubic-spline-interpolation");
17
18// Flags to enable extra information for investigating splines
19//#define SHOWCOEFFICIENTS 1
20//#define GNUPLOT 1
21
22const int NUM_ITERATIONS = 24;
23
24TestSpline::TestSpline(QObject *parent) :
25 QObject(parent)
26{
27}
28
29void TestSpline::cleanupTestCase ()
30{
31
32}
33
34bool TestSpline::coefCheckX (const vector<double> &t,
35#ifdef SHOWCOEFFICIENTS
36 const vector<SplinePair> &xy,
37#else
38 const vector<SplinePair> & /* xy */,
39#endif
40 const Spline &s) const
41{
42 unsigned int i;
43 double aUntranslated = 0, bUntranslated = 0, cUntranslated = 0, dUntranslated = 0;
44
45 // X coordinate fit
46#ifdef SHOWCOEFFICIENTS
47 cout << endl
48 << "(t,x) inputs to be copied to " << WEBPAGE.toLatin1().data()
49 << endl;
50 for (i = 0; i < t.size(); i++) {
51 cout << t[i] << " " << xy[i].x() << endl;
52 }
53 cout << endl
54 << "x=d*(t-ti)^3+c*(t-ti)^2+b*(t-ti)+a natural cubic spline results to be used in this code"
55 << endl;
56 for (i = 0; i < t.size() - 1; i++) {
57 coefShow ("x =",
58 "(t-ti)",
59 t[i],
60 t[i+1],
61 s.m_elements[i].a().x(),
62 s.m_elements[i].b().x(),
63 s.m_elements[i].c().x(),
64 s.m_elements[i].d().x());
65 }
66 cout << endl
67 << "x=d*t^3+c*t^2+b*t+a outputs to be compared to results from " << WEBPAGE.toLatin1().data()
68 << endl;
69#endif
70 for (i = 0; i < t.size() - 1; i++) {
71 s.computeUntranslatedCoefficients (s.m_elements[i].a().x(),
72 s.m_elements[i].b().x(),
73 s.m_elements[i].c().x(),
74 s.m_elements[i].d().x(),
75 t[i],
76 aUntranslated,
77 bUntranslated,
78 cUntranslated,
79 dUntranslated);
80#ifdef SHOWCOEFFICIENTS
81 coefShow ("x =",
82 "t",
83 t[i],
84 t[i+1],
85 aUntranslated,
86 bUntranslated,
87 cUntranslated,
88 dUntranslated);
89#endif
90 }
91
92 // Spot check the (arbitrarily chosen) final row with results from the WEBPAGE
93 bool success = true;
94 success &= (qAbs (aUntranslated - -8.3608) < 8.3608 / 10000.0);
95 success &= (qAbs (bUntranslated - 4.2505) < 4.2505 / 10000.0);
96 success &= (qAbs (cUntranslated - -0.63092) < 0.63092 / 10000.0);
97 success &= (qAbs (dUntranslated - 0.035051) < 0.035051 / 10000.0);
98
99 return success;
100}
101
102bool TestSpline::coefCheckY (const vector<double> &t,
103#ifdef SHOWCOEFFICIENTS
104 const vector<SplinePair> &xy,
105#else
106 const vector<SplinePair> & /* xy */,
107#endif
108 const Spline &s) const
109{
110 unsigned int i;
111 double aUntranslated = 0, bUntranslated = 0, cUntranslated = 0, dUntranslated = 0;
112
113 // Y coordinate fit
114#ifdef SHOWCOEFFICIENTS
115 cout << endl
116 << "(t,y) inputs to be copied to " << WEBPAGE.toLatin1().data()
117 << endl;
118 for (i = 0; i < xy.size(); i++) {
119 cout << t[i] << " " << xy[i].y() << endl;
120 }
121 cout << endl
122 << "y=d*(t-ti)^3+c*(t-ti)^2+b*(t-ti)+a natural cubic spline results to be used in this code"
123 << endl;
124 for (i = 0; i < xy.size() - 1; i++) {
125 coefShow ("y =",
126 "(t-ti)",
127 t[i],
128 t[i+1],
129 s.m_elements[i].a().y(),
130 s.m_elements[i].b().y(),
131 s.m_elements[i].c().y(),
132 s.m_elements[i].d().y());
133 }
134 cout << endl
135 << "y=d*t^3+c*t^2+b*t+a outputs to be compared to results from " << WEBPAGE.toLatin1().data()
136 << endl;
137#endif
138 for (i = 0; i < t.size() - 1; i++) {
139 s.computeUntranslatedCoefficients (s.m_elements[i].a().y(),
140 s.m_elements[i].b().y(),
141 s.m_elements[i].c().y(),
142 s.m_elements[i].d().y(),
143 t[i],
144 aUntranslated,
145 bUntranslated,
146 cUntranslated,
147 dUntranslated);
148#ifdef SHOWCOEFFICIENTS
149 coefShow ("y =",
150 "t",
151 t[i],
152 t[i+1],
153 aUntranslated,
154 bUntranslated,
155 cUntranslated,
156 dUntranslated);
157#endif
158 }
159
160 // Spot check the (arbitrarily chosen) final row with results from the WEBPAGE
161 bool success = true;
162 success &= (qAbs (aUntranslated - -7.0) < 7.0 / 10000.0);
163 success &= (qAbs (bUntranslated - 3.5667) < 3.5667 / 10000.0);
164 success &= (qAbs (cUntranslated - -0.6) < 0.6 / 10000.0);
165 success &= (qAbs (dUntranslated - 0.033333) < 0.033333 / 10000.0);
166
167 return success;
168}
169
170void TestSpline::coefShow (const QString &leftHandSide,
171 const QString &independentVariable,
172 double tLow,
173 double tHigh,
174 double a,
175 double b,
176 double c,
177 double d) const
178{
179 cout << leftHandSide.toLatin1().data() << scientific
180 << d << "*" << independentVariable.toLatin1().data() << "^3 + "
181 << c << "*" << independentVariable.toLatin1().data() << "^2 + "
182 << b << "*" << independentVariable.toLatin1().data() << " + "
183 << a << " (" << tLow << "<t<" << tHigh << ")" << endl;
184
185}
186
187void TestSpline::initTestCase ()
188{
189 const bool NO_DROP_REGRESSION = false;
190 const QString NO_ERROR_REPORT_LOG_FILE;
191 const QString NO_REGRESSION_OPEN_FILE;
192 const bool NO_GNUPLOT_LOG_FILES = false;
193 const bool NO_REGRESSION_IMPORT = false;
194 const bool NO_RESET = false;
195 const bool NO_EXPORT_ONLY = false;
196 const bool NO_EXPORT_IMAGE_ONLY = false;
197 const QString NO_EXPORT_IMAGE_EXTENSION;
198 const bool DEBUG_FLAG = false;
199 const QStringList NO_LOAD_STARTUP_FILES;
200 const QStringList NO_COMMAND_LINE;
201
202 initializeLogging ("engauge_test",
203 "engauge_test.log",
204 DEBUG_FLAG);
205
206 MainWindow w (NO_ERROR_REPORT_LOG_FILE,
211 NO_RESET,
213 NO_EXPORT_IMAGE_ONLY,
214 NO_EXPORT_IMAGE_EXTENSION,
217 w.show ();
218}
219
220void TestSpline::testCoefficientsFromOrdinals ()
221{
222 bool success = true;
223 vector<double> t;
224 vector<SplinePair> xy;
225
226 // Input data for testSharpTransition
227 xy.push_back (SplinePair (0.1, 1.0));
228 xy.push_back (SplinePair (0.5, 1.0));
229 xy.push_back (SplinePair (0.8, 1.0));
230 xy.push_back (SplinePair (1.0, 0.5));
231 xy.push_back (SplinePair (1.01, 0));
232 xy.push_back (SplinePair (1.5, 0.0));
233 xy.push_back (SplinePair (2.0, 0.0));
234
235 int counter = 0;
236 vector<SplinePair>::const_iterator itr;
237 for (itr = xy.begin(); itr != xy.end(); itr++) {
238 t.push_back (counter++);
239 }
240
241 // Generate the spline
242 Spline s (t, xy);
243
244 success &= coefCheckX (t,
245 xy,
246 s);
247 success &= coefCheckY (t,
248 xy,
249 s);
250
251 QVERIFY(success);
252}
253
254void TestSpline::testSharpTransition ()
255{
256 const int NUM_T = 60;
257 const double SPLINE_EPSILON = 0.01;
258
259 map<double, bool> xMerged; // Preserve order
260
261 bool success = true;
262
263 vector<double> t;
264 vector<SplinePair> xyBefore;
265
266 // Values with a sharp transition that can trigger overlap (=not a function)
267 xyBefore.push_back (SplinePair (0.1, 1.0));
268 xyBefore.push_back (SplinePair (0.5, 1.0));
269 xyBefore.push_back (SplinePair (0.8, 1.0));
270 xyBefore.push_back (SplinePair (1.0, 0.5));
271 xyBefore.push_back (SplinePair (1.01, 0));
272 xyBefore.push_back (SplinePair (1.5, 0.0));
273 xyBefore.push_back (SplinePair (2.0, 0.0));
274
275 // Copy x values into t vector and initial version of xMerged vector
276 vector<SplinePair>::const_iterator itrB;
277 for (itrB = xyBefore.begin(); itrB != xyBefore.end(); itrB++) {
278 const SplinePair &pair = *itrB;
279 t.push_back (pair.x());
280 xMerged [pair.x ()] = true;
281 }
282
283 // Merge many more plotting points into xMerged
284 double tStart = t[0];
285 double tStop = t[t.size() - 1];
286 for (int i = 0; i <= NUM_T; i++) {
287 double t = tStart + (double) i * (tStop - tStart) / (double) NUM_T;
288
289 if (xMerged.find (t) == xMerged.end ()) {
290 xMerged [t] = true;
291 }
292 }
293
294 // Generate the spline
295 Spline s (t, xyBefore, SPLINE_DISABLE_T_CHECK);
296
297 // Plot the points after generating the spline
298 vector<SplinePair> xyAfter;
299 map<double, bool>::const_iterator itrX;
300 for (itrX = xMerged.begin(); itrX != xMerged.end(); itrX++) {
301 double x = itrX->first;
302 SplinePair pair = s.interpolateCoeff (x);
303 xyAfter.push_back (pair);
304 }
305
306#ifdef GNUPLOT
307 cout << "set datafile missing \"?\"" << endl;
308 cout << "plot \"gnuplot.in\" using 1:2 with linespoints, \"gnuplot.in\" using 1:3 with lines" << endl;
309#endif
310
311 // Output the merged before/after curves
312 map<double, bool>::const_iterator itrM;
313 for (itrM = xMerged.begin(); itrM != xMerged.end(); itrM++) {
314 double x = itrM->first;
315
316 string yB = "?", yA = "?";
317
318 vector<SplinePair>::iterator itrB;
319 for (itrB = xyBefore.begin(); itrB != xyBefore.end(); itrB++) {
320 if (itrB->x() == x) {
321 ostringstream osB;
322 osB << itrB->y();
323 yB = osB.str();
324 break;
325 }
326 }
327
328 vector<SplinePair>::iterator itrA;
329 for (itrA = xyAfter.begin(); itrA != xyAfter.end(); itrA++) {
330 if (itrA->x() == x) {
331 ostringstream osA;
332 osA << itrA->y();
333 yA = osA.str();
334 break;
335 }
336 }
337
338 if (itrB != xyBefore.end() &&
339 itrA != xyAfter.end()) {
340
341 // This x value has y values from both before and after that we can compare for success
342 double yBefore = itrB->y();
343 double yAfter = itrA->y();
344 if (qAbs (yBefore - yAfter) > SPLINE_EPSILON) {
345 success = false;
346 }
347 }
348
349#ifdef GNUPLOT
350 cout << x << ", " << yB << ", " << yA << endl;
351#endif
352 }
353
354 QVERIFY (success);
355}
356
357void TestSpline::testSplinesAsControlPoints ()
358{
359 const int T_START = 1, T_STOP = 7;
360 const double SPLINE_EPSILON = 0.01;
361 const int NUM_T = 60;
362
363 bool success = true;
364
365 vector<double> t;
366 vector<SplinePair> xy;
367
368 // Independent variable are evenly spaced
369 t.push_back (T_START);
370 t.push_back (2);
371 t.push_back (3);
372 t.push_back (4);
373 t.push_back (5);
374 t.push_back (6);
375 t.push_back (T_STOP);
376
377 // Simple curve, with x values tweaked slightly (from even spacing) to make the test data more stressing
378 xy.push_back (SplinePair (1, 0.22));
379 xy.push_back (SplinePair (1.8, 0.04));
380 xy.push_back (SplinePair (3.2, -0.13));
381 xy.push_back (SplinePair (4.3, -0.17));
382 xy.push_back (SplinePair (5, -0.04));
383 xy.push_back (SplinePair (5.8, 0.09));
384 xy.push_back (SplinePair (7, 0.11));
385
386 Spline s (t, xy);
387
388 for (int i = 0; i <= NUM_T; i++) {
389 double t = T_START + (double) i * (T_STOP - T_START) / (double) NUM_T;
390 SplinePair spCoeff = s.interpolateCoeff (t);
391 SplinePair spBezier = s.interpolateControlPoints (t);
392
393 double xCoeff = spCoeff.x();
394 double yCoeff = spCoeff.y();
395 double xControl = spBezier.x();
396 double yControl = spBezier.y();
397
398 if (qAbs (xCoeff - xControl) > SPLINE_EPSILON) {
399 success = false;
400 }
401
402 if (qAbs (yCoeff - yControl) > SPLINE_EPSILON) {
403 success = false;
404 }
405 }
406
407 QVERIFY (success);
408}
void initializeLogging(const QString &name, const QString &filename, bool isDebug)
Definition Logger.cpp:21
@ SPLINE_DISABLE_T_CHECK
Definition Spline.h:16
const bool NO_EXPORT_ONLY
const QStringList NO_COMMAND_LINE
const QString NO_ERROR_REPORT_LOG_FILE
const bool NO_GNUPLOT_LOG_FILES
const QString NO_REGRESSION_OPEN_FILE
const QStringList NO_LOAD_STARTUP_FILES
const bool NO_REGRESSION_IMPORT
const bool NO_DROP_REGRESSION
const bool DEBUG_FLAG
const int NUM_ITERATIONS
const QString WEBPAGE("https://tools.timodenk.com/cubic-spline-interpolation")
double y() const
Get method for y.
double x() const
Get method for x.
Cubic interpolation given independent and dependent value vectors.
Definition Spline.h:30
Unit test of spline library.
Definition TestSpline.h:13
TestSpline(QObject *parent=0)
Single constructor.