Engauge Digitizer  2
 All Classes Functions Variables Typedefs Enumerations Friends Pages
TestSpline.cpp
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 
12 QTEST_MAIN (TestSpline)
13 
14 using namespace std;
15 
16 const 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 
22 const int NUM_ITERATIONS = 24;
23 
24 TestSpline::TestSpline(QObject *parent) :
25  QObject(parent)
26 {
27 }
28 
29 void TestSpline::cleanupTestCase ()
30 {
31 
32 }
33 
34 bool 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 
102 bool 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 
170 void 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 
187 void TestSpline::initTestCase ()
188 {
189  const QString NO_ERROR_REPORT_LOG_FILE;
190  const QString NO_REGRESSION_OPEN_FILE;
191  const bool NO_GNUPLOT_LOG_FILES = false;
192  const bool NO_REGRESSION_IMPORT = false;
193  const bool NO_RESET = false;
194  const bool NO_EXPORT_ONLY = false;
195  const bool NO_EXPORT_IMAGE_ONLY = false;
196  const QString NO_EXPORT_IMAGE_EXTENSION;
197  const bool DEBUG_FLAG = false;
198  const QStringList NO_LOAD_STARTUP_FILES;
199  const QStringList NO_COMMAND_LINE;
200 
201  initializeLogging ("engauge_test",
202  "engauge_test.log",
203  DEBUG_FLAG);
204 
205  MainWindow w (NO_ERROR_REPORT_LOG_FILE,
206  NO_REGRESSION_OPEN_FILE,
207  NO_REGRESSION_IMPORT,
208  NO_GNUPLOT_LOG_FILES,
209  NO_RESET,
210  NO_EXPORT_ONLY,
211  NO_EXPORT_IMAGE_ONLY,
212  NO_EXPORT_IMAGE_EXTENSION,
213  NO_LOAD_STARTUP_FILES,
214  NO_COMMAND_LINE);
215  w.show ();
216 }
217 
218 void TestSpline::testCoefficientsFromOrdinals ()
219 {
220  bool success = true;
221  vector<double> t;
222  vector<SplinePair> xy;
223 
224  // Input data for testSharpTransition
225  xy.push_back (SplinePair (0.1, 1.0));
226  xy.push_back (SplinePair (0.5, 1.0));
227  xy.push_back (SplinePair (0.8, 1.0));
228  xy.push_back (SplinePair (1.0, 0.5));
229  xy.push_back (SplinePair (1.01, 0));
230  xy.push_back (SplinePair (1.5, 0.0));
231  xy.push_back (SplinePair (2.0, 0.0));
232 
233  int counter = 0;
234  vector<SplinePair>::const_iterator itr;
235  for (itr = xy.begin(); itr != xy.end(); itr++) {
236  t.push_back (counter++);
237  }
238 
239  // Generate the spline
240  Spline s (t, xy);
241 
242  success &= coefCheckX (t,
243  xy,
244  s);
245  success &= coefCheckY (t,
246  xy,
247  s);
248 
249  QVERIFY(success);
250 }
251 
252 void TestSpline::testSharpTransition ()
253 {
254  const int NUM_T = 60;
255  const double SPLINE_EPSILON = 0.01;
256 
257  map<double, bool> xMerged; // Preserve order
258 
259  bool success = true;
260 
261  vector<double> t;
262  vector<SplinePair> xyBefore;
263 
264  // Values with a sharp transition that can trigger overlap (=not a function)
265  xyBefore.push_back (SplinePair (0.1, 1.0));
266  xyBefore.push_back (SplinePair (0.5, 1.0));
267  xyBefore.push_back (SplinePair (0.8, 1.0));
268  xyBefore.push_back (SplinePair (1.0, 0.5));
269  xyBefore.push_back (SplinePair (1.01, 0));
270  xyBefore.push_back (SplinePair (1.5, 0.0));
271  xyBefore.push_back (SplinePair (2.0, 0.0));
272 
273  // Copy x values into t vector and initial version of xMerged vector
274  vector<SplinePair>::const_iterator itrB;
275  for (itrB = xyBefore.begin(); itrB != xyBefore.end(); itrB++) {
276  const SplinePair &pair = *itrB;
277  t.push_back (pair.x());
278  xMerged [pair.x ()] = true;
279  }
280 
281  // Merge many more plotting points into xMerged
282  double tStart = t[0];
283  double tStop = t[t.size() - 1];
284  for (int i = 0; i <= NUM_T; i++) {
285  double t = tStart + (double) i * (tStop - tStart) / (double) NUM_T;
286 
287  if (xMerged.find (t) == xMerged.end ()) {
288  xMerged [t] = true;
289  }
290  }
291 
292  // Generate the spline
293  Spline s (t, xyBefore, SPLINE_DISABLE_T_CHECK);
294 
295  // Plot the points after generating the spline
296  vector<SplinePair> xyAfter;
297  map<double, bool>::const_iterator itrX;
298  for (itrX = xMerged.begin(); itrX != xMerged.end(); itrX++) {
299  double x = itrX->first;
300  SplinePair pair = s.interpolateCoeff (x);
301  xyAfter.push_back (pair);
302  }
303 
304 #ifdef GNUPLOT
305  cout << "set datafile missing \"?\"" << endl;
306  cout << "plot \"gnuplot.in\" using 1:2 with linespoints, \"gnuplot.in\" using 1:3 with lines" << endl;
307 #endif
308 
309  // Output the merged before/after curves
310  map<double, bool>::const_iterator itrM;
311  for (itrM = xMerged.begin(); itrM != xMerged.end(); itrM++) {
312  double x = itrM->first;
313 
314  string yB = "?", yA = "?";
315 
316  vector<SplinePair>::iterator itrB;
317  for (itrB = xyBefore.begin(); itrB != xyBefore.end(); itrB++) {
318  if (itrB->x() == x) {
319  ostringstream osB;
320  osB << itrB->y();
321  yB = osB.str();
322  break;
323  }
324  }
325 
326  vector<SplinePair>::iterator itrA;
327  for (itrA = xyAfter.begin(); itrA != xyAfter.end(); itrA++) {
328  if (itrA->x() == x) {
329  ostringstream osA;
330  osA << itrA->y();
331  yA = osA.str();
332  break;
333  }
334  }
335 
336  if (itrB != xyBefore.end() &&
337  itrA != xyAfter.end()) {
338 
339  // This x value has y values from both before and after that we can compare for success
340  double yBefore = itrB->y();
341  double yAfter = itrA->y();
342  if (qAbs (yBefore - yAfter) > SPLINE_EPSILON) {
343  success = false;
344  }
345  }
346 
347 #ifdef GNUPLOT
348  cout << x << ", " << yB << ", " << yA << endl;
349 #endif
350  }
351 
352  QVERIFY (success);
353 }
354 
355 void TestSpline::testSplinesAsControlPoints ()
356 {
357  const int T_START = 1, T_STOP = 7;
358  const double SPLINE_EPSILON = 0.01;
359  const int NUM_T = 60;
360 
361  bool success = true;
362 
363  vector<double> t;
364  vector<SplinePair> xy;
365 
366  // Independent variable are evenly spaced
367  t.push_back (T_START);
368  t.push_back (2);
369  t.push_back (3);
370  t.push_back (4);
371  t.push_back (5);
372  t.push_back (6);
373  t.push_back (T_STOP);
374 
375  // Simple curve, with x values tweaked slightly (from even spacing) to make the test data more stressing
376  xy.push_back (SplinePair (1, 0.22));
377  xy.push_back (SplinePair (1.8, 0.04));
378  xy.push_back (SplinePair (3.2, -0.13));
379  xy.push_back (SplinePair (4.3, -0.17));
380  xy.push_back (SplinePair (5, -0.04));
381  xy.push_back (SplinePair (5.8, 0.09));
382  xy.push_back (SplinePair (7, 0.11));
383 
384  Spline s (t, xy);
385 
386  for (int i = 0; i <= NUM_T; i++) {
387  double t = T_START + (double) i * (T_STOP - T_START) / (double) NUM_T;
388  SplinePair spCoeff = s.interpolateCoeff (t);
389  SplinePair spBezier = s.interpolateControlPoints (t);
390 
391  double xCoeff = spCoeff.x();
392  double yCoeff = spCoeff.y();
393  double xControl = spBezier.x();
394  double yControl = spBezier.y();
395 
396  if (qAbs (xCoeff - xControl) > SPLINE_EPSILON) {
397  success = false;
398  }
399 
400  if (qAbs (yCoeff - yControl) > SPLINE_EPSILON) {
401  success = false;
402  }
403  }
404 
405  QVERIFY (success);
406 }
Cubic interpolation given independent and dependent value vectors.
Definition: Spline.h:29
double y() const
Get method for y.
Definition: SplinePair.cpp:88
double x() const
Get method for x.
Definition: SplinePair.cpp:83
Unit test of spline library.
Definition: TestSpline.h:12
Main window consisting of menu, graphics scene, status bar and optional toolbars as a Single Document...
Definition: MainWindow.h:91
Single X/Y pair for cubic spline interpolation initialization and calculations.
Definition: SplinePair.h:13