libs/math/doc/interpolators/catmull_rom.qbk

   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
  16 namespace 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
  41 Catmull-Rom splines are a family of interpolating curves which are commonly used in computer graphics and animation.
  42 Catmull-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.
  48 Many curves (such as B[eacute]zier) are /approximating/ - they do not pass through their control points.
  49 This makes them more difficult to use than interpolating splines.
  50
  51 The `catmull_rom` class provided by Boost.Math creates a cubic Catmull-Rom spline from an array of points in any dimension.
  52 Since there are numerous ways to represent a point in /n/-dimensional space,
  53 the class attempts to be flexible by templating on the point type.
  54 The 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),
  55 it 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()`.
  56 These requirements are met by `std::vector` and `std::array`.
  57 The 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
  68 The spline can be either open or /closed/, closed meaning that there is some /s > 0/ such that /P(s) = P(0)/.
  69 The 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
  74 In either case, evaluating the interpolator at /s=0/ returns the first point in the list.
  75
  76 If the curve is open, then the first and last segments may have strange behavior.
  77 The 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.
  78 A control point is appended to the data using similar criteria.
  79 However, 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
  81 Inside `catmull_rom`, the curve is represented as closed.
  82 This is because an open Catmull-Rom curve is /implicitly closed/, but the closing point is the zero vector.
  83 There is no reason to suppose that the zero vector is a better closing point than the endpoint (or any other point, for that matter),
  84 so traditionally Catmull-Rom splines leave the segment between the first and second point undefined,
  85 as well as the segment between the second-to-last and last point.
  86 We find this property of the traditional implementation of Catmull-Rom splines annoying and confusing to the user.
  87 Hence internally, we close the curve so that the first and last segments are defined.
  88 Of course, this causes the /tangent/ vectors to the first and last points to be bizarre.
  89 This is a "pick your poison" design decision-either the curve cannot interpolate in its first and last segments,
  90 or the tangents along the first and last segments are meaningless.
  91 In the vast majority of cases, this will be no problem to the user.
  92 However, 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
  94 Since the routine internally represents the curve as closed,
  95 a question arises: Why does the user have to specify if the curve is open or closed?
  96 The answer is that the parameterization is chosen by the routine,
  97 so 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
 101 If you attempt to interpolate for `s > max_s`, an exception is thrown.
 102 If the curve is closed, then `cr(max_s) = p0`, where `p0` is the first point on the curve.
 103 If the curve is open, then `cr(max_s) = pf`, where `pf` is the final point on the curve.
 104
 105
 106 The Catmull-Rom curve admits an infinite number of parameterizations.
 107 The default parameterization of the `catmull_rom` class is the so-called /centripedal/ parameterization.
 108 This parameterization has been shown to be the only parameterization that does not form cusps or self-intersections within segments.
 109 However, 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
 114 The alpha parameter must always be in the range `[0,1]`.
 115
 116 Finally, 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
 121 Since the magnitude of the tangent vector is dependent on the parameterization,
 122 it is not meaningful (unless the user chooses the chordal parameterization /alpha = 1/ which parameterizes by Euclidean distance between points.)
 123 However, 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
 131 The following performance numbers were generated for a call to the Catmull-Rom interpolation method.
 132 The number that follows the slash is the number of points passed to the interpolant.
 133 We 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
 170 We have already discussed that certain conditions on the `Point` type template argument must be obeyed.
 171 The 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
 210 These 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
 215 The Catmull-Rom interpolator requires memory for three more points than is provided by the user.
 216 This causes the class to call a `resize()` on the input vector.
 217 If `v.capacity() >= v.size() + 3`, then no problems arise; there are no reallocs, and in practice this condition is almost always satisfied.
 218 However, if `v.capacity() < v.size() + 3`, the `realloc` causes a performance penalty of roughly 20%.
 219
 220 [heading Generic Containers]
 221
 222 The `Point` type may be stored in a different container than `std::vector`.
 223 For 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]