1[/ 2 Copyright 2017 Nick Thompson 3 4 Distributed under the Boost Software License, Version 1.0. 5 (See accompanying file LICENSE_1_0.txt or copy at 6 http://www.boost.org/LICENSE_1_0.txt). 7] 8 9[section:catmull_rom Catmull-Rom Splines] 10 11[heading Synopsis] 12 13`` 14#include <boost/math/interpolators/catmull_rom.hpp> 15 16namespace boost{ namespace math{ 17 18 template<class Point, class RandomAccessContainer = std::vector<Point> > 19 class catmull_rom 20 { 21 public: 22 23 catmull_rom(RandomAccessContainer&& points, bool closed = false, Real alpha = (Real) 1/ (Real) 2) 24 25 catmull_rom(std::initializer_list<Point> l, bool closed = false, typename Point::value_type alpha = (typename Point::value_type) 1/ (typename Point::value_type) 2); 26 27 Real operator()(Real s) const; 28 29 Real max_parameter() const; 30 31 Real parameter_at_point(size_t i) const; 32 33 Point prime(Real s) const; 34 }; 35 36}} 37`` 38 39[heading Description] 40 41Catmull-Rom splines are a family of interpolating curves which are commonly used in computer graphics and animation. 42Catmull-Rom splines enjoy the following properties: 43 44* Affine invariance: The interpolant commutes with affine transformations. 45* Local support of the basis functions: This gives stability and fast evaluation. 46* /C/[super 2]-smoothness 47* Interpolation of control points-this means the curve passes through the control points. 48Many curves (such as B[eacute]zier) are /approximating/ - they do not pass through their control points. 49This makes them more difficult to use than interpolating splines. 50 51The `catmull_rom` class provided by Boost.Math creates a cubic Catmull-Rom spline from an array of points in any dimension. 52Since there are numerous ways to represent a point in /n/-dimensional space, 53the class attempts to be flexible by templating on the point type. 54The requirements on the point type are discussing in more detail below, but roughly, it must have a dereference operator defined (e.g., `p[0]` is not a syntax error), 55it must be able to be dereferenced up to `dimension -1`, and `p[i]` is of type `Real`, define `value_type`, and the free function `size()`. 56These requirements are met by `std::vector` and `std::array`. 57The basic usage is shown here: 58 59 std::vector<std::array<Real, 3>> points(4); 60 points[0] = {0,0,0}; 61 points[1] = {1,0,0}; 62 points[2] = {0,1,0}; 63 points[3] = {0,0,1}; 64 catmull_rom<std::array<Real, 3>> cr(std::move(points)); 65 // Interpolate at s = 0.1: 66 auto point = cr(0.1); 67 68The spline can be either open or /closed/, closed meaning that there is some /s > 0/ such that /P(s) = P(0)/. 69The default is open, but this can be easily changed: 70 71 // closed = true 72 catmull_rom<std::array<Real, 3>> cr(std::move(points), true); 73 74In either case, evaluating the interpolator at /s=0/ returns the first point in the list. 75 76If the curve is open, then the first and last segments may have strange behavior. 77The traditional solution is to prepend a carefully selected control point to the data so that the first data segment (second interpolator segment) has reasonable tangent vectors, and simply ignore the first interpolator segment. 78A control point is appended to the data using similar criteria. 79However, we recommend not going through this effort until it proves to be necessary: For most use-cases, the curve is good enough without prepending and appending control points, and responsible selection of non-data control points is difficult. 80 81Inside `catmull_rom`, the curve is represented as closed. 82This is because an open Catmull-Rom curve is /implicitly closed/, but the closing point is the zero vector. 83There is no reason to suppose that the zero vector is a better closing point than the endpoint (or any other point, for that matter), 84so traditionally Catmull-Rom splines leave the segment between the first and second point undefined, 85as well as the segment between the second-to-last and last point. 86We find this property of the traditional implementation of Catmull-Rom splines annoying and confusing to the user. 87Hence internally, we close the curve so that the first and last segments are defined. 88Of course, this causes the /tangent/ vectors to the first and last points to be bizarre. 89This is a "pick your poison" design decision-either the curve cannot interpolate in its first and last segments, 90or the tangents along the first and last segments are meaningless. 91In the vast majority of cases, this will be no problem to the user. 92However, if it becomes a problem, then the user should add one extra point in a position they believe is reasonable and close the curve. 93 94Since the routine internally represents the curve as closed, 95a question arises: Why does the user have to specify if the curve is open or closed? 96The answer is that the parameterization is chosen by the routine, 97so it is of interest to the user to understand the values where a meaningful result is returned. 98 99 Real max_s = cr.max_parameter(); 100 101If you attempt to interpolate for `s > max_s`, an exception is thrown. 102If the curve is closed, then `cr(max_s) = p0`, where `p0` is the first point on the curve. 103If the curve is open, then `cr(max_s) = pf`, where `pf` is the final point on the curve. 104 105 106The Catmull-Rom curve admits an infinite number of parameterizations. 107The default parameterization of the `catmull_rom` class is the so-called /centripedal/ parameterization. 108This parameterization has been shown to be the only parameterization that does not form cusps or self-intersections within segments. 109However, for advanced users, other parameterizations can be chosen using the /alpha/ parameter: 110 111 // alpha = 1 is the "chordal" parameterization. 112 catmull_rom<std::array<double, 3>> cr(std::move(points), false, 1.0); 113 114The alpha parameter must always be in the range `[0,1]`. 115 116Finally, the tangent vector to any point of the curve can be computed via 117 118 double s = 0.1; 119 Point tangent = cr.prime(s); 120 121Since the magnitude of the tangent vector is dependent on the parameterization, 122it is not meaningful (unless the user chooses the chordal parameterization /alpha = 1/ which parameterizes by Euclidean distance between points.) 123However, its direction is meaningful no matter the parameterization, so the user may wish to normalize this result. 124 125[heading Examples] 126 127[import ../../example/catmull_rom_example.cpp] 128 129[heading Performance] 130 131The following performance numbers were generated for a call to the Catmull-Rom interpolation method. 132The number that follows the slash is the number of points passed to the interpolant. 133We see that evaluation of the interpolant is [bigo](/log/(/N/)). 134 135 136 Run on 2700 MHz CPU 137 CPU Caches: 138 L1 Data 32K (x2) 139 L1 Instruction 32K (x2) 140 L2 Unified 262K (x2) 141 L3 Unified 3145K (x1) 142 --------------------------------------------------------- 143 Benchmark Time CPU 144 --------------------------------------------------------- 145 BM_CatmullRom<double>/4 20 ns 20 ns 146 BM_CatmullRom<double>/8 21 ns 21 ns 147 BM_CatmullRom<double>/16 23 ns 23 ns 148 BM_CatmullRom<double>/32 24 ns 24 ns 149 BM_CatmullRom<double>/64 27 ns 27 ns 150 BM_CatmullRom<double>/128 27 ns 27 ns 151 BM_CatmullRom<double>/256 30 ns 30 ns 152 BM_CatmullRom<double>/512 32 ns 31 ns 153 BM_CatmullRom<double>/1024 33 ns 33 ns 154 BM_CatmullRom<double>/2048 34 ns 34 ns 155 BM_CatmullRom<double>/4096 36 ns 36 ns 156 BM_CatmullRom<double>/8192 38 ns 38 ns 157 BM_CatmullRom<double>/16384 39 ns 39 ns 158 BM_CatmullRom<double>/32768 40 ns 40 ns 159 BM_CatmullRom<double>/65536 45 ns 44 ns 160 BM_CatmullRom<double>/131072 46 ns 46 ns 161 BM_CatmullRom<double>/262144 50 ns 50 ns 162 BM_CatmullRom<double>/524288 53 ns 52 ns 163 BM_CatmullRom<double>/1048576 58 ns 57 ns 164 BM_CatmullRom<double>_BigO 2.97 lgN 2.97 lgN 165 BM_CatmullRom<double>_RMS 19 % 19 % 166 167 168[heading Point types] 169 170We have already discussed that certain conditions on the `Point` type template argument must be obeyed. 171The following shows a custom point type in 3D which can be used as a template argument to Catmull-Rom: 172 173 template<class Real> 174 class mypoint3d 175 { 176 public: 177 // Must define a value_type: 178 typedef Real value_type; 179 180 // Regular constructor--need not be of this form. 181 mypoint3d(Real x, Real y, Real z) {m_vec[0] = x; m_vec[1] = y; m_vec[2] = z; } 182 183 // Must define a default constructor: 184 mypoint3d() {} 185 186 // Must define array access: 187 Real operator[](size_t i) const 188 { 189 return m_vec[i]; 190 } 191 192 // Must define array element assignment: 193 Real& operator[](size_t i) 194 { 195 return m_vec[i]; 196 } 197 198 private: 199 std::array<Real, 3> m_vec; 200 }; 201 202 203 // Must define the free function "size()": 204 template<class Real> 205 constexpr size_t size(const mypoint3d<Real>& c) 206 { 207 return 3; 208 } 209 210These conditions are satisfied by both `std::array` and `std::vector`, but it may nonetheless be useful to define your own point class so that (say) you can define geometric distance between them. 211 212 213[heading Caveats] 214 215The Catmull-Rom interpolator requires memory for three more points than is provided by the user. 216This causes the class to call a `resize()` on the input vector. 217If `v.capacity() >= v.size() + 3`, then no problems arise; there are no reallocs, and in practice this condition is almost always satisfied. 218However, if `v.capacity() < v.size() + 3`, the `realloc` causes a performance penalty of roughly 20%. 219 220[heading Generic Containers] 221 222The `Point` type may be stored in a different container than `std::vector`. 223For example, here is how to store the points in a Boost.uBLAS vector: 224 225 mypoint3d<Real> p0(0.1, 0.2, 0.3); 226 mypoint3d<Real> p1(0.2, 0.3, 0.4); 227 mypoint3d<Real> p2(0.3, 0.4, 0.5); 228 mypoint3d<Real> p3(0.4, 0.5, 0.6); 229 mypoint3d<Real> p4(0.5, 0.6, 0.7); 230 mypoint3d<Real> p5(0.6, 0.7, 0.8); 231 232 boost::numeric::ublas::vector<mypoint3d<Real>> u(6); 233 u[0] = p0; 234 u[1] = p1; 235 u[2] = p2; 236 u[3] = p3; 237 u[4] = p4; 238 u[5] = p5; 239 240 // Tests initializer_list: 241 catmull_rom<mypoint3d<Real>, decltype(u)> cat(std::move(u)); 242 243 244[heading References] 245 246* Cem Yuksel, Scott Schaefer, and John Keyser, ['Parameterization and applications of Catmull–Rom curves], Computer-Aided Design 43 (2011) 747–755. 247* Phillip J. Barry and Ronald N. Goldman, ['A Recursive Evaluation Algorithm for a Class of Catmull-Rom Splines], Computer Graphics, Volume 22, Number 4, August 1988 248 249[endsect] [/section:catmull_rom Catmull-Rom Splines] 250