Home | History | Annotate | Download | only in test 
      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