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
14 namespace Eigen {
15
16 // lets do some explicit instantiations and thus
17 // force the compilation of all spline functions...
18 template class Spline<double, 2, Dynamic>;
19 template class Spline<double, 3, Dynamic>;
20
21 template class Spline<double, 2, 2>;
22 template class Spline<double, 2, 3>;
23 template class Spline<double, 2, 4>;
24 template class Spline<double, 2, 5>;
25
26 template class Spline<float, 2, Dynamic>;
27 template class Spline<float, 3, Dynamic>;
28
29 template class Spline<float, 3, 2>;
30 template class Spline<float, 3, 3>;
31 template class Spline<float, 3, 4>;
32 template class Spline<float, 3, 5>;
33
34 }
35
closed_spline2d()36 Spline<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 */
spline3d()67 Spline<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 */
eval_spline3d()97 void 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 */
eval_spline3d_onbrks()134 void 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
eval_closed_spline2d()161 void 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
check_global_interpolation2d()201 void 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
check_global_interpolation_with_derivatives2d()237 void check_global_interpolation_with_derivatives2d()
238 {
239 typedef Spline2d::PointType PointType;
240 typedef Spline2d::KnotVectorType KnotVectorType;
241
242 const Eigen::DenseIndex numPoints = 100;
243 const unsigned int dimension = 2;
244 const unsigned int degree = 3;
245
246 ArrayXXd points = ArrayXXd::Random(dimension, numPoints);
247
248 KnotVectorType knots;
249 Eigen::ChordLengths(points, knots);
250
251 ArrayXXd derivatives = ArrayXXd::Random(dimension, numPoints);
252 VectorXd derivativeIndices(numPoints);
253
254 for (Eigen::DenseIndex i = 0; i < numPoints; ++i)
255 derivativeIndices(i) = static_cast<double>(i);
256
257 const Spline2d spline = SplineFitting<Spline2d>::InterpolateWithDerivatives(
258 points, derivatives, derivativeIndices, degree);
259
260 for (Eigen::DenseIndex i = 0; i < points.cols(); ++i)
261 {
262 PointType point = spline(knots(i));
263 PointType referencePoint = points.col(i);
264 VERIFY_IS_APPROX(point, referencePoint);
265 PointType derivative = spline.derivatives(knots(i), 1).col(1);
266 PointType referenceDerivative = derivatives.col(i);
267 VERIFY_IS_APPROX(derivative, referenceDerivative);
268 }
269 }
270
test_splines()271 void test_splines()
272 {
273 for (int i = 0; i < g_repeat; ++i)
274 {
275 CALL_SUBTEST( eval_spline3d() );
276 CALL_SUBTEST( eval_spline3d_onbrks() );
277 CALL_SUBTEST( eval_closed_spline2d() );
278 CALL_SUBTEST( check_global_interpolation2d() );
279 CALL_SUBTEST( check_global_interpolation_with_derivatives2d() );
280 }
281 }
282