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