1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2010-2011 Hauke Heibel <heibel@gmail.com>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#include "main.h"
11
12#include <unsupported/Eigen/Splines>
13
14namespace Eigen {
15
16// lets do some explicit instantiations and thus
17// force the compilation of all spline functions...
18template class Spline<double, 2, Dynamic>;
19template class Spline<double, 3, Dynamic>;
20
21template class Spline<double, 2, 2>;
22template class Spline<double, 2, 3>;
23template class Spline<double, 2, 4>;
24template class Spline<double, 2, 5>;
25
26template class Spline<float, 2, Dynamic>;
27template class Spline<float, 3, Dynamic>;
28
29template class Spline<float, 3, 2>;
30template class Spline<float, 3, 3>;
31template class Spline<float, 3, 4>;
32template class Spline<float, 3, 5>;
33
34}
35
36Spline<double, 2, Dynamic> closed_spline2d()
37{
38  RowVectorXd knots(12);
39  knots << 0,
40    0,
41    0,
42    0,
43    0.867193179093898,
44    1.660330955342408,
45    2.605084834823134,
46    3.484154586374428,
47    4.252699478956276,
48    4.252699478956276,
49    4.252699478956276,
50    4.252699478956276;
51
52  MatrixXd ctrls(8,2);
53  ctrls << -0.370967741935484,   0.236842105263158,
54    -0.231401860693277,   0.442245185027632,
55    0.344361228532831,   0.773369994120753,
56    0.828990216203802,   0.106550882647595,
57    0.407270163678382,  -1.043452922172848,
58    -0.488467813584053,  -0.390098582530090,
59    -0.494657189446427,   0.054804824897884,
60    -0.370967741935484,   0.236842105263158;
61  ctrls.transposeInPlace();
62
63  return Spline<double, 2, Dynamic>(knots, ctrls);
64}
65
66/* create a reference spline */
67Spline<double, 3, Dynamic> spline3d()
68{
69  RowVectorXd knots(11);
70  knots << 0,
71    0,
72    0,
73    0.118997681558377,
74    0.162611735194631,
75    0.498364051982143,
76    0.655098003973841,
77    0.679702676853675,
78    1.000000000000000,
79    1.000000000000000,
80    1.000000000000000;
81
82  MatrixXd ctrls(8,3);
83  ctrls <<    0.959743958516081,   0.340385726666133,   0.585267750979777,
84    0.223811939491137,   0.751267059305653,   0.255095115459269,
85    0.505957051665142,   0.699076722656686,   0.890903252535799,
86    0.959291425205444,   0.547215529963803,   0.138624442828679,
87    0.149294005559057,   0.257508254123736,   0.840717255983663,
88    0.254282178971531,   0.814284826068816,   0.243524968724989,
89    0.929263623187228,   0.349983765984809,   0.196595250431208,
90    0.251083857976031,   0.616044676146639,   0.473288848902729;
91  ctrls.transposeInPlace();
92
93  return Spline<double, 3, Dynamic>(knots, ctrls);
94}
95
96/* compares evaluations against known results */
97void eval_spline3d()
98{
99  Spline3d spline = spline3d();
100
101  RowVectorXd u(10);
102  u << 0.351659507062997,
103    0.830828627896291,
104    0.585264091152724,
105    0.549723608291140,
106    0.917193663829810,
107    0.285839018820374,
108    0.757200229110721,
109    0.753729094278495,
110    0.380445846975357,
111    0.567821640725221;
112
113  MatrixXd pts(10,3);
114  pts << 0.707620811535916,   0.510258911240815,   0.417485437023409,
115    0.603422256426978,   0.529498282727551,   0.270351549348981,
116    0.228364197569334,   0.423745615677815,   0.637687289287490,
117    0.275556796335168,   0.350856706427970,   0.684295784598905,
118    0.514519311047655,   0.525077224890754,   0.351628308305896,
119    0.724152914315666,   0.574461155457304,   0.469860285484058,
120    0.529365063753288,   0.613328702656816,   0.237837040141739,
121    0.522469395136878,   0.619099658652895,   0.237139665242069,
122    0.677357023849552,   0.480655768435853,   0.422227610314397,
123    0.247046593173758,   0.380604672404750,   0.670065791405019;
124  pts.transposeInPlace();
125
126  for (int i=0; i<u.size(); ++i)
127  {
128    Vector3d pt = spline(u(i));
129    VERIFY( (pt - pts.col(i)).norm() < 1e-14 );
130  }
131}
132
133/* compares evaluations on corner cases */
134void eval_spline3d_onbrks()
135{
136  Spline3d spline = spline3d();
137
138  RowVectorXd u = spline.knots();
139
140  MatrixXd pts(11,3);
141  pts <<    0.959743958516081,   0.340385726666133,   0.585267750979777,
142    0.959743958516081,   0.340385726666133,   0.585267750979777,
143    0.959743958516081,   0.340385726666133,   0.585267750979777,
144    0.430282980289940,   0.713074680056118,   0.720373307943349,
145    0.558074875553060,   0.681617921034459,   0.804417124839942,
146    0.407076008291750,   0.349707710518163,   0.617275937419545,
147    0.240037008286602,   0.738739390398014,   0.324554153129411,
148    0.302434111480572,   0.781162443963899,   0.240177089094644,
149    0.251083857976031,   0.616044676146639,   0.473288848902729,
150    0.251083857976031,   0.616044676146639,   0.473288848902729,
151    0.251083857976031,   0.616044676146639,   0.473288848902729;
152  pts.transposeInPlace();
153
154  for (int i=0; i<u.size(); ++i)
155  {
156    Vector3d pt = spline(u(i));
157    VERIFY( (pt - pts.col(i)).norm() < 1e-14 );
158  }
159}
160
161void eval_closed_spline2d()
162{
163  Spline2d spline = closed_spline2d();
164
165  RowVectorXd u(12);
166  u << 0,
167    0.332457030395796,
168    0.356467130532952,
169    0.453562180176215,
170    0.648017921874804,
171    0.973770235555003,
172    1.882577647219307,
173    2.289408593930498,
174    3.511951429883045,
175    3.884149321369450,
176    4.236261590369414,
177    4.252699478956276;
178
179  MatrixXd pts(12,2);
180  pts << -0.370967741935484,   0.236842105263158,
181    -0.152576775123250,   0.448975001279334,
182    -0.133417538277668,   0.461615613865667,
183    -0.053199060826740,   0.507630360006299,
184    0.114249591147281,   0.570414135097409,
185    0.377810316891987,   0.560497102875315,
186    0.665052120135908,  -0.157557441109611,
187    0.516006487053228,  -0.559763292174825,
188    -0.379486035348887,  -0.331959640488223,
189    -0.462034726249078,  -0.039105670080824,
190    -0.378730600917982,   0.225127015099919,
191    -0.370967741935484,   0.236842105263158;
192  pts.transposeInPlace();
193
194  for (int i=0; i<u.size(); ++i)
195  {
196    Vector2d pt = spline(u(i));
197    VERIFY( (pt - pts.col(i)).norm() < 1e-14 );
198  }
199}
200
201void check_global_interpolation2d()
202{
203  typedef Spline2d::PointType PointType;
204  typedef Spline2d::KnotVectorType KnotVectorType;
205  typedef Spline2d::ControlPointVectorType ControlPointVectorType;
206
207  ControlPointVectorType points = ControlPointVectorType::Random(2,100);
208
209  KnotVectorType chord_lengths; // knot parameters
210  Eigen::ChordLengths(points, chord_lengths);
211
212  // interpolation without knot parameters
213  {
214    const Spline2d spline = SplineFitting<Spline2d>::Interpolate(points,3);
215
216    for (Eigen::DenseIndex i=0; i<points.cols(); ++i)
217    {
218      PointType pt = spline( chord_lengths(i) );
219      PointType ref = points.col(i);
220      VERIFY( (pt - ref).matrix().norm() < 1e-14 );
221    }
222  }
223
224  // interpolation with given knot parameters
225  {
226    const Spline2d spline = SplineFitting<Spline2d>::Interpolate(points,3,chord_lengths);
227
228    for (Eigen::DenseIndex i=0; i<points.cols(); ++i)
229    {
230      PointType pt = spline( chord_lengths(i) );
231      PointType ref = points.col(i);
232      VERIFY( (pt - ref).matrix().norm() < 1e-14 );
233    }
234  }
235}
236
237
238void test_splines()
239{
240  CALL_SUBTEST( eval_spline3d() );
241  CALL_SUBTEST( eval_spline3d_onbrks() );
242  CALL_SUBTEST( eval_closed_spline2d() );
243  CALL_SUBTEST( check_global_interpolation2d() );
244}
245