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 (at) 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 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 */ 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 */ 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 */ 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 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 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 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 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