libs/math/test/catmull_rom_test.cpp

   1 /*
   2  * Copyright Nick Thompson, 2017
   3  * Use, modification and distribution are subject to the
   4  * Boost Software License, Version 1.0. (See accompanying file
   5  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6  */
   7 #define BOOST_TEST_MODULE catmull_rom_test
   8
   9 #include <array>
  10 #include <random>
  11 #include <boost/cstdfloat.hpp>
  12 #include <boost/type_index.hpp>
  13 #include <boost/test/included/unit_test.hpp>
  14 #include <boost/test/tools/floating_point_comparison.hpp>
  15 #include <boost/math/constants/constants.hpp>
  16 #include <boost/math/interpolators/catmull_rom.hpp>
  17 #include <boost/multiprecision/cpp_bin_float.hpp>
  18 #include <boost/multiprecision/cpp_dec_float.hpp>
  19 #include <boost/numeric/ublas/vector.hpp>
  20
  21 using std::abs;
  22 using boost::multiprecision::cpp_bin_float_50;
  23 using boost::math::catmull_rom;
  24
  25 template<class Real>
  26 void test_alpha_distance()
  27 {
  28     Real tol = std::numeric_limits<Real>::epsilon();
  29     std::array<Real, 3> v1 = {0,0,0};
  30     std::array<Real, 3> v2 = {1,0,0};
  31     Real alpha = 0.5;
  32     Real d = boost::math::detail::alpha_distance<std::array<Real, 3>>(v1, v2, alpha);
  33     BOOST_CHECK_CLOSE_FRACTION(d, 1, tol);
  34
  35     d = boost::math::detail::alpha_distance<std::array<Real, 3>>(v1, v2, 0.0);
  36     BOOST_CHECK_CLOSE_FRACTION(d, 1, tol);
  37
  38     d = boost::math::detail::alpha_distance<std::array<Real, 3>>(v1, v2, 1.0);
  39     BOOST_CHECK_CLOSE_FRACTION(d, 1, tol);
  40
  41     v2[0] = 2;
  42     d = boost::math::detail::alpha_distance<std::array<Real, 3>>(v1, v2, alpha);
  43     BOOST_CHECK_CLOSE_FRACTION(d, pow(2, (Real)1/ (Real) 2), tol);
  44
  45     d = boost::math::detail::alpha_distance<std::array<Real, 3>>(v1, v2, 0.0);
  46     BOOST_CHECK_CLOSE_FRACTION(d, 1, tol);
  47
  48     d = boost::math::detail::alpha_distance<std::array<Real, 3>>(v1, v2, 1.0);
  49     BOOST_CHECK_CLOSE_FRACTION(d, 2, tol);
  50 }
  51
  52
  53 template<class Real>
  54 void test_linear()
  55 {
  56     std::cout << "Testing that the Catmull-Rom spline interpolates linear functions correctly on type "
  57               << boost::typeindex::type_id<Real>().pretty_name() << "\n";
  58
  59     Real tol = 10*std::numeric_limits<Real>::epsilon();
  60     std::vector<std::array<Real, 3>> v(4);
  61     v[0] = {0,0,0};
  62     v[1] = {1,0,0};
  63     v[2] = {2,0,0};
  64     v[3] = {3,0,0};
  65     catmull_rom<std::array<Real, 3>> cr(std::move(v));
  66
  67     // Test that the interpolation condition is obeyed:
  68     BOOST_CHECK_CLOSE_FRACTION(cr.max_parameter(), 3, tol);
  69     auto p0 = cr(0.0);
  70     BOOST_CHECK_SMALL(p0[0], tol);
  71     BOOST_CHECK_SMALL(p0[1], tol);
  72     BOOST_CHECK_SMALL(p0[2], tol);
  73     auto p1 = cr(1.0);
  74     BOOST_CHECK_CLOSE_FRACTION(p1[0], 1, tol);
  75     BOOST_CHECK_SMALL(p1[1], tol);
  76     BOOST_CHECK_SMALL(p1[2], tol);
  77
  78     auto p2 = cr(2.0);
  79     BOOST_CHECK_CLOSE_FRACTION(p2[0], 2, tol);
  80     BOOST_CHECK_SMALL(p2[1], tol);
  81     BOOST_CHECK_SMALL(p2[2], tol);
  82
  83
  84     auto p3 = cr(3.0);
  85     BOOST_CHECK_CLOSE_FRACTION(p3[0], 3, tol);
  86     BOOST_CHECK_SMALL(p3[1], tol);
  87     BOOST_CHECK_SMALL(p3[2], tol);
  88
  89     Real s = cr.parameter_at_point(0);
  90     BOOST_CHECK_SMALL(s, tol);
  91
  92     s = cr.parameter_at_point(1);
  93     BOOST_CHECK_CLOSE_FRACTION(s, 1, tol);
  94
  95     s = cr.parameter_at_point(2);
  96     BOOST_CHECK_CLOSE_FRACTION(s, 2, tol);
  97
  98     s = cr.parameter_at_point(3);
  99     BOOST_CHECK_CLOSE_FRACTION(s, 3, tol);
 100
 101     // Test that the function is linear on the interval [1,2]:
 102     for (double s = 1; s < 2; s += 0.01)
 103     {
 104         auto p = cr(s);
 105         BOOST_CHECK_CLOSE_FRACTION(p[0], s, tol);
 106         BOOST_CHECK_SMALL(p[1], tol);
 107         BOOST_CHECK_SMALL(p[2], tol);
 108
 109         auto tangent = cr.prime(s);
 110         BOOST_CHECK_CLOSE_FRACTION(tangent[0], 1, tol);
 111         BOOST_CHECK_SMALL(tangent[1], tol);
 112         BOOST_CHECK_SMALL(tangent[2], tol);
 113     }
 114
 115 }
 116
 117 template<class Real>
 118 void test_circle()
 119 {
 120     using boost::math::constants::pi;
 121     using std::cos;
 122     using std::sin;
 123
 124     std::cout << "Testing that the Catmull-Rom spline interpolates circles correctly on type "
 125               << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 126
 127     Real tol = 10*std::numeric_limits<Real>::epsilon();
 128     std::vector<std::array<Real, 2>> v(20*sizeof(Real));
 129     std::vector<std::array<Real, 2>> u(20*sizeof(Real));
 130     for (size_t i = 0; i < v.size(); ++i)
 131     {
 132         Real theta = ((Real) i/ (Real) v.size())*2*pi<Real>();
 133         v[i] = {cos(theta), sin(theta)};
 134         u[i] = v[i];
 135     }
 136     catmull_rom<std::array<Real, 2>> circle(std::move(v), true);
 137
 138     // Interpolation condition:
 139     for (size_t i = 0; i < v.size(); ++i)
 140     {
 141         Real s = circle.parameter_at_point(i);
 142         auto p = circle(s);
 143         Real x = p[0];
 144         Real y = p[1];
 145         if (abs(x) < std::numeric_limits<Real>::epsilon())
 146         {
 147             BOOST_CHECK_SMALL(u[i][0], tol);
 148         }
 149         if (abs(y) < std::numeric_limits<Real>::epsilon())
 150         {
 151             BOOST_CHECK_SMALL(u[i][1], tol);
 152         }
 153         else
 154         {
 155             BOOST_CHECK_CLOSE_FRACTION(x, u[i][0], tol);
 156             BOOST_CHECK_CLOSE_FRACTION(y, u[i][1], tol);
 157         }
 158     }
 159
 160     Real max_s = circle.max_parameter();
 161     for(Real s = 0; s < max_s; s += 0.01)
 162     {
 163         auto p = circle(s);
 164         Real x = p[0];
 165         Real y = p[1];
 166         BOOST_CHECK_CLOSE_FRACTION(x*x+y*y, 1, 0.001);
 167     }
 168 }
 169
 170
 171 template<class Real, size_t dimension>
 172 void test_affine_invariance()
 173 {
 174     std::cout << "Testing that the Catmull-Rom spline is affine invariant in dimension "
 175               << dimension << " on type "
 176               << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 177
 178     Real tol = 1000*std::numeric_limits<Real>::epsilon();
 179     std::vector<std::array<Real, dimension>> v(100);
 180     std::vector<std::array<Real, dimension>> u(100);
 181     std::mt19937_64 gen(438232);
 182     Real inv_denom = (Real) 100/( (Real) (gen.max)() + (Real) 2);
 183     for(size_t j = 0; j < dimension; ++j)
 184     {
 185         v[0][j] = gen()*inv_denom;
 186         u[0][j] = v[0][j];
 187     }
 188
 189     for (size_t i = 1; i < v.size(); ++i)
 190     {
 191         for(size_t j = 0; j < dimension; ++j)
 192         {
 193             v[i][j] = v[i-1][j] + gen()*inv_denom;
 194             u[i][j] = v[i][j];
 195         }
 196     }
 197     std::array<Real, dimension> affine_shift;
 198     for (size_t j = 0; j < dimension; ++j)
 199     {
 200         affine_shift[j] = gen()*inv_denom;
 201     }
 202
 203     catmull_rom<std::array<Real, dimension>> cr1(std::move(v));
 204
 205     for(size_t i = 0; i< u.size(); ++i)
 206     {
 207         for(size_t j = 0; j < dimension; ++j)
 208         {
 209             u[i][j] += affine_shift[j];
 210         }
 211     }
 212
 213     catmull_rom<std::array<Real, dimension>> cr2(std::move(u));
 214
 215     BOOST_CHECK_CLOSE_FRACTION(cr1.max_parameter(), cr2.max_parameter(), tol);
 216
 217     Real ds = cr1.max_parameter()/1024;
 218     for (Real s = 0; s < cr1.max_parameter(); s += ds)
 219     {
 220         auto p0 = cr1(s);
 221         auto p1 = cr2(s);
 222         auto tangent0 = cr1.prime(s);
 223         auto tangent1 = cr2.prime(s);
 224         for (size_t j = 0; j < dimension; ++j)
 225         {
 226             BOOST_CHECK_CLOSE_FRACTION(p0[j] + affine_shift[j], p1[j], tol);
 227             if (abs(tangent0[j]) > 5000*tol)
 228             {
 229                 BOOST_CHECK_CLOSE_FRACTION(tangent0[j], tangent1[j], 5000*tol);
 230             }
 231         }
 232     }
 233 }
 234
 235 template<class Real>
 236 void test_helix()
 237 {
 238     using boost::math::constants::pi;
 239     std::cout << "Testing that the Catmull-Rom spline interpolates helices correctly on type "
 240               << boost::typeindex::type_id<Real>().pretty_name() << "\n";
 241
 242     Real tol = 0.001;
 243     std::vector<std::array<Real, 3>> v(400*sizeof(Real));
 244     for (size_t i = 0; i < v.size(); ++i)
 245     {
 246         Real theta = ((Real) i/ (Real) v.size())*2*pi<Real>();
 247         v[i] = {cos(theta), sin(theta), theta};
 248     }
 249     catmull_rom<std::array<Real, 3>> helix(std::move(v));
 250
 251     // Interpolation condition:
 252     for (size_t i = 0; i < v.size(); ++i)
 253     {
 254         Real s = helix.parameter_at_point(i);
 255         auto p = helix(s);
 256         Real t = p[2];
 257
 258         Real x = p[0];
 259         Real y = p[1];
 260         if (abs(x) < tol)
 261         {
 262             BOOST_CHECK_SMALL(cos(t), tol);
 263         }
 264         if (abs(y) < tol)
 265         {
 266             BOOST_CHECK_SMALL(sin(t), tol);
 267         }
 268         else
 269         {
 270             BOOST_CHECK_CLOSE_FRACTION(x, cos(t), tol);
 271             BOOST_CHECK_CLOSE_FRACTION(y, sin(t), tol);
 272         }
 273     }
 274
 275     Real max_s = helix.max_parameter();
 276     for(Real s = helix.parameter_at_point(1); s < max_s; s += 0.01)
 277     {
 278         auto p = helix(s);
 279         Real x = p[0];
 280         Real y = p[1];
 281         Real t = p[2];
 282         BOOST_CHECK_CLOSE_FRACTION(x*x+y*y, (Real) 1, (Real) 0.01);
 283         if (abs(x) < 0.01)
 284         {
 285             BOOST_CHECK_SMALL(cos(t),  (Real) 0.05);
 286         }
 287         if (abs(y) < 0.01)
 288         {
 289             BOOST_CHECK_SMALL(sin(t), (Real) 0.05);
 290         }
 291         else
 292         {
 293             BOOST_CHECK_CLOSE_FRACTION(x, cos(t), (Real) 0.05);
 294             BOOST_CHECK_CLOSE_FRACTION(y, sin(t), (Real) 0.05);
 295         }
 296     }
 297 }
 298
 299
 300 template<class Real>
 301 class mypoint3d
 302 {
 303 public:
 304     // Must define a value_type:
 305     typedef Real value_type;
 306
 307     // Regular constructor:
 308     mypoint3d(Real x, Real y, Real z)
 309     {
 310         m_vec[0] = x;
 311         m_vec[1] = y;
 312         m_vec[2] = z;
 313     }
 314
 315     // Must define a default constructor:
 316     mypoint3d() {}
 317
 318     // Must define array access:
 319     Real operator[](size_t i) const
 320     {
 321         return m_vec[i];
 322     }
 323
 324     // Array element assignment:
 325     Real& operator[](size_t i)
 326     {
 327         return m_vec[i];
 328     }
 329
 330
 331 private:
 332     std::array<Real, 3>  m_vec;
 333 };
 334
 335
 336 // Must define the free function "size()":
 337 template<class Real>
 338 BOOST_CONSTEXPR std::size_t size(const mypoint3d<Real>& c)
 339 {
 340     return 3;
 341 }
 342
 343 template<class Real>
 344 void test_data_representations()
 345 {
 346     std::cout << "Testing that the Catmull-Rom spline works with multiple data representations.\n";
 347     mypoint3d<Real> p0(0.1, 0.2, 0.3);
 348     mypoint3d<Real> p1(0.2, 0.3, 0.4);
 349     mypoint3d<Real> p2(0.3, 0.4, 0.5);
 350     mypoint3d<Real> p3(0.4, 0.5, 0.6);
 351     mypoint3d<Real> p4(0.5, 0.6, 0.7);
 352     mypoint3d<Real> p5(0.6, 0.7, 0.8);
 353
 354
 355     // Tests initializer_list:
 356     catmull_rom<mypoint3d<Real>> cat({p0, p1, p2, p3, p4, p5});
 357
 358     Real tol = 0.001;
 359     auto p = cat(cat.parameter_at_point(0));
 360     BOOST_CHECK_CLOSE_FRACTION(p[0], p0[0], tol);
 361     BOOST_CHECK_CLOSE_FRACTION(p[1], p0[1], tol);
 362     BOOST_CHECK_CLOSE_FRACTION(p[2], p0[2], tol);
 363     p = cat(cat.parameter_at_point(1));
 364     BOOST_CHECK_CLOSE_FRACTION(p[0], p1[0], tol);
 365     BOOST_CHECK_CLOSE_FRACTION(p[1], p1[1], tol);
 366     BOOST_CHECK_CLOSE_FRACTION(p[2], p1[2], tol);
 367 }
 368
 369 template<class Real>
 370 void test_random_access_container()
 371 {
 372     std::cout << "Testing that the Catmull-Rom spline works with multiple data representations.\n";
 373     mypoint3d<Real> p0(0.1, 0.2, 0.3);
 374     mypoint3d<Real> p1(0.2, 0.3, 0.4);
 375     mypoint3d<Real> p2(0.3, 0.4, 0.5);
 376     mypoint3d<Real> p3(0.4, 0.5, 0.6);
 377     mypoint3d<Real> p4(0.5, 0.6, 0.7);
 378     mypoint3d<Real> p5(0.6, 0.7, 0.8);
 379
 380     boost::numeric::ublas::vector<mypoint3d<Real>> u(6);
 381     u[0] = p0;
 382     u[1] = p1;
 383     u[2] = p2;
 384     u[3] = p3;
 385     u[4] = p4;
 386     u[5] = p5;
 387
 388     // Tests initializer_list:
 389     catmull_rom<mypoint3d<Real>, decltype(u)> cat(std::move(u));
 390
 391     Real tol = 0.001;
 392     auto p = cat(cat.parameter_at_point(0));
 393     BOOST_CHECK_CLOSE_FRACTION(p[0], p0[0], tol);
 394     BOOST_CHECK_CLOSE_FRACTION(p[1], p0[1], tol);
 395     BOOST_CHECK_CLOSE_FRACTION(p[2], p0[2], tol);
 396     p = cat(cat.parameter_at_point(1));
 397     BOOST_CHECK_CLOSE_FRACTION(p[0], p1[0], tol);
 398     BOOST_CHECK_CLOSE_FRACTION(p[1], p1[1], tol);
 399     BOOST_CHECK_CLOSE_FRACTION(p[2], p1[2], tol);
 400 }
 401
 402 BOOST_AUTO_TEST_CASE(catmull_rom_test)
 403 {
 404 #if !defined(TEST) || (TEST == 1)
 405     test_data_representations<float>();
 406     test_alpha_distance<double>();
 407
 408     test_linear<double>();
 409     test_linear<long double>();
 410
 411     test_circle<float>();
 412     test_circle<double>();
 413 #endif
 414 #if !defined(TEST) || (TEST == 2)
 415     test_helix<double>();
 416
 417     test_affine_invariance<double, 1>();
 418     test_affine_invariance<double, 2>();
 419     test_affine_invariance<double, 3>();
 420     test_affine_invariance<double, 4>();
 421
 422     test_random_access_container<double>();
 423 #endif
 424 #if !defined(TEST) || (TEST == 3)
 425     test_affine_invariance<cpp_bin_float_50, 4>();
 426 #endif
 427 }