]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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> | |
1e59de90 | 338 | constexpr std::size_t size(const mypoint3d<Real>& c) |
92f5a8d4 TL |
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>(); | |
1e59de90 | 409 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS |
92f5a8d4 | 410 | test_linear<long double>(); |
1e59de90 | 411 | #endif |
92f5a8d4 TL |
412 | |
413 | test_circle<float>(); | |
414 | test_circle<double>(); | |
415 | #endif | |
416 | #if !defined(TEST) || (TEST == 2) | |
417 | test_helix<double>(); | |
418 | ||
419 | test_affine_invariance<double, 1>(); | |
420 | test_affine_invariance<double, 2>(); | |
421 | test_affine_invariance<double, 3>(); | |
422 | test_affine_invariance<double, 4>(); | |
423 | ||
424 | test_random_access_container<double>(); | |
425 | #endif | |
426 | #if !defined(TEST) || (TEST == 3) | |
427 | test_affine_invariance<cpp_bin_float_50, 4>(); | |
428 | #endif | |
429 | } |