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 
237 
test_splines()238 void 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