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