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 
    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