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)
7 #define BOOST_TEST_MODULE catmull_rom_test
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>
22 using boost::multiprecision::cpp_bin_float_50
;
23 using boost::math::catmull_rom
;
26 void test_alpha_distance()
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};
32 Real d
= boost::math::detail::alpha_distance
<std::array
<Real
, 3>>(v1
, v2
, alpha
);
33 BOOST_CHECK_CLOSE_FRACTION(d
, 1, tol
);
35 d
= boost::math::detail::alpha_distance
<std::array
<Real
, 3>>(v1
, v2
, 0.0);
36 BOOST_CHECK_CLOSE_FRACTION(d
, 1, tol
);
38 d
= boost::math::detail::alpha_distance
<std::array
<Real
, 3>>(v1
, v2
, 1.0);
39 BOOST_CHECK_CLOSE_FRACTION(d
, 1, tol
);
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
);
45 d
= boost::math::detail::alpha_distance
<std::array
<Real
, 3>>(v1
, v2
, 0.0);
46 BOOST_CHECK_CLOSE_FRACTION(d
, 1, tol
);
48 d
= boost::math::detail::alpha_distance
<std::array
<Real
, 3>>(v1
, v2
, 1.0);
49 BOOST_CHECK_CLOSE_FRACTION(d
, 2, tol
);
56 std::cout
<< "Testing that the Catmull-Rom spline interpolates linear functions correctly on type "
57 << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
59 Real tol
= 10*std::numeric_limits
<Real
>::epsilon();
60 std::vector
<std::array
<Real
, 3>> v(4);
65 catmull_rom
<std::array
<Real
, 3>> cr(std::move(v
));
67 // Test that the interpolation condition is obeyed:
68 BOOST_CHECK_CLOSE_FRACTION(cr
.max_parameter(), 3, tol
);
70 BOOST_CHECK_SMALL(p0
[0], tol
);
71 BOOST_CHECK_SMALL(p0
[1], tol
);
72 BOOST_CHECK_SMALL(p0
[2], tol
);
74 BOOST_CHECK_CLOSE_FRACTION(p1
[0], 1, tol
);
75 BOOST_CHECK_SMALL(p1
[1], tol
);
76 BOOST_CHECK_SMALL(p1
[2], tol
);
79 BOOST_CHECK_CLOSE_FRACTION(p2
[0], 2, tol
);
80 BOOST_CHECK_SMALL(p2
[1], tol
);
81 BOOST_CHECK_SMALL(p2
[2], tol
);
85 BOOST_CHECK_CLOSE_FRACTION(p3
[0], 3, tol
);
86 BOOST_CHECK_SMALL(p3
[1], tol
);
87 BOOST_CHECK_SMALL(p3
[2], tol
);
89 Real s
= cr
.parameter_at_point(0);
90 BOOST_CHECK_SMALL(s
, tol
);
92 s
= cr
.parameter_at_point(1);
93 BOOST_CHECK_CLOSE_FRACTION(s
, 1, tol
);
95 s
= cr
.parameter_at_point(2);
96 BOOST_CHECK_CLOSE_FRACTION(s
, 2, tol
);
98 s
= cr
.parameter_at_point(3);
99 BOOST_CHECK_CLOSE_FRACTION(s
, 3, tol
);
101 // Test that the function is linear on the interval [1,2]:
102 for (double s
= 1; s
< 2; s
+= 0.01)
105 BOOST_CHECK_CLOSE_FRACTION(p
[0], s
, tol
);
106 BOOST_CHECK_SMALL(p
[1], tol
);
107 BOOST_CHECK_SMALL(p
[2], tol
);
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
);
120 using boost::math::constants::pi
;
124 std::cout
<< "Testing that the Catmull-Rom spline interpolates circles correctly on type "
125 << boost::typeindex::type_id
<Real
>().pretty_name() << "\n";
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
)
132 Real theta
= ((Real
) i
/ (Real
) v
.size())*2*pi
<Real
>();
133 v
[i
] = {cos(theta
), sin(theta
)};
136 catmull_rom
<std::array
<Real
, 2>> circle(std::move(v
), true);
138 // Interpolation condition:
139 for (size_t i
= 0; i
< v
.size(); ++i
)
141 Real s
= circle
.parameter_at_point(i
);
145 if (abs(x
) < std::numeric_limits
<Real
>::epsilon())
147 BOOST_CHECK_SMALL(u
[i
][0], tol
);
149 if (abs(y
) < std::numeric_limits
<Real
>::epsilon())
151 BOOST_CHECK_SMALL(u
[i
][1], tol
);
155 BOOST_CHECK_CLOSE_FRACTION(x
, u
[i
][0], tol
);
156 BOOST_CHECK_CLOSE_FRACTION(y
, u
[i
][1], tol
);
160 Real max_s
= circle
.max_parameter();
161 for(Real s
= 0; s
< max_s
; s
+= 0.01)
166 BOOST_CHECK_CLOSE_FRACTION(x
*x
+y
*y
, 1, 0.001);
171 template<class Real
, size_t dimension
>
172 void test_affine_invariance()
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";
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
)
185 v
[0][j
] = gen()*inv_denom
;
189 for (size_t i
= 1; i
< v
.size(); ++i
)
191 for(size_t j
= 0; j
< dimension
; ++j
)
193 v
[i
][j
] = v
[i
-1][j
] + gen()*inv_denom
;
197 std::array
<Real
, dimension
> affine_shift
;
198 for (size_t j
= 0; j
< dimension
; ++j
)
200 affine_shift
[j
] = gen()*inv_denom
;
203 catmull_rom
<std::array
<Real
, dimension
>> cr1(std::move(v
));
205 for(size_t i
= 0; i
< u
.size(); ++i
)
207 for(size_t j
= 0; j
< dimension
; ++j
)
209 u
[i
][j
] += affine_shift
[j
];
213 catmull_rom
<std::array
<Real
, dimension
>> cr2(std::move(u
));
215 BOOST_CHECK_CLOSE_FRACTION(cr1
.max_parameter(), cr2
.max_parameter(), tol
);
217 Real ds
= cr1
.max_parameter()/1024;
218 for (Real s
= 0; s
< cr1
.max_parameter(); s
+= ds
)
222 auto tangent0
= cr1
.prime(s
);
223 auto tangent1
= cr2
.prime(s
);
224 for (size_t j
= 0; j
< dimension
; ++j
)
226 BOOST_CHECK_CLOSE_FRACTION(p0
[j
] + affine_shift
[j
], p1
[j
], tol
);
227 if (abs(tangent0
[j
]) > 5000*tol
)
229 BOOST_CHECK_CLOSE_FRACTION(tangent0
[j
], tangent1
[j
], 5000*tol
);
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";
243 std::vector
<std::array
<Real
, 3>> v(400*sizeof(Real
));
244 for (size_t i
= 0; i
< v
.size(); ++i
)
246 Real theta
= ((Real
) i
/ (Real
) v
.size())*2*pi
<Real
>();
247 v
[i
] = {cos(theta
), sin(theta
), theta
};
249 catmull_rom
<std::array
<Real
, 3>> helix(std::move(v
));
251 // Interpolation condition:
252 for (size_t i
= 0; i
< v
.size(); ++i
)
254 Real s
= helix
.parameter_at_point(i
);
262 BOOST_CHECK_SMALL(cos(t
), tol
);
266 BOOST_CHECK_SMALL(sin(t
), tol
);
270 BOOST_CHECK_CLOSE_FRACTION(x
, cos(t
), tol
);
271 BOOST_CHECK_CLOSE_FRACTION(y
, sin(t
), tol
);
275 Real max_s
= helix
.max_parameter();
276 for(Real s
= helix
.parameter_at_point(1); s
< max_s
; s
+= 0.01)
282 BOOST_CHECK_CLOSE_FRACTION(x
*x
+y
*y
, (Real
) 1, (Real
) 0.01);
285 BOOST_CHECK_SMALL(cos(t
), (Real
) 0.05);
289 BOOST_CHECK_SMALL(sin(t
), (Real
) 0.05);
293 BOOST_CHECK_CLOSE_FRACTION(x
, cos(t
), (Real
) 0.05);
294 BOOST_CHECK_CLOSE_FRACTION(y
, sin(t
), (Real
) 0.05);
304 // Must define a value_type:
305 typedef Real value_type
;
307 // Regular constructor:
308 mypoint3d(Real x
, Real y
, Real z
)
315 // Must define a default constructor:
318 // Must define array access:
319 Real
operator[](size_t i
) const
324 // Array element assignment:
325 Real
& operator[](size_t i
)
332 std::array
<Real
, 3> m_vec
;
336 // Must define the free function "size()":
338 constexpr std::size_t size(const mypoint3d
<Real
>& c
)
344 void test_data_representations()
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);
355 // Tests initializer_list:
356 catmull_rom
<mypoint3d
<Real
>> cat({p0
, p1
, p2
, p3
, p4
, p5
});
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
);
370 void test_random_access_container()
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);
380 boost::numeric::ublas::vector
<mypoint3d
<Real
>> u(6);
388 // Tests initializer_list:
389 catmull_rom
<mypoint3d
<Real
>, decltype(u
)> cat(std::move(u
));
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
);
402 BOOST_AUTO_TEST_CASE(catmull_rom_test
)
404 #if !defined(TEST) || (TEST == 1)
405 test_data_representations
<float>();
406 test_alpha_distance
<double>();
408 test_linear
<double>();
409 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
410 test_linear
<long double>();
413 test_circle
<float>();
414 test_circle
<double>();
416 #if !defined(TEST) || (TEST == 2)
417 test_helix
<double>();
419 test_affine_invariance
<double, 1>();
420 test_affine_invariance
<double, 2>();
421 test_affine_invariance
<double, 3>();
422 test_affine_invariance
<double, 4>();
424 test_random_access_container
<double>();
426 #if !defined(TEST) || (TEST == 3)
427 test_affine_invariance
<cpp_bin_float_50
, 4>();