libsidplayfp  1.0.3
spline.h
00001 //  ---------------------------------------------------------------------------
00002 //  This file is part of reSID, a MOS6581 SID emulator engine.
00003 //  Copyright (C) 2010  Dag Lem <resid@nimrod.no>
00004 //
00005 //  This program is free software; you can redistribute it and/or modify
00006 //  it under the terms of the GNU General Public License as published by
00007 //  the Free Software Foundation; either version 2 of the License, or
00008 //  (at your option) any later version.
00009 //
00010 //  This program is distributed in the hope that it will be useful,
00011 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013 //  GNU General Public License for more details.
00014 //
00015 //  You should have received a copy of the GNU General Public License
00016 //  along with this program; if not, write to the Free Software
00017 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00018 //  ---------------------------------------------------------------------------
00019 
00020 #ifndef RESID_SPLINE_H
00021 #define RESID_SPLINE_H
00022 
00023 namespace reSID
00024 {
00025 
00026 // Our objective is to construct a smooth interpolating single-valued function
00027 // y = f(x).
00028 //
00029 // Catmull-Rom splines are widely used for interpolation, however these are
00030 // parametric curves [x(t) y(t) ...] and can not be used to directly calculate
00031 // y = f(x).
00032 // For a discussion of Catmull-Rom splines see Catmull, E., and R. Rom,
00033 // "A Class of Local Interpolating Splines", Computer Aided Geometric Design.
00034 //
00035 // Natural cubic splines are single-valued functions, and have been used in
00036 // several applications e.g. to specify gamma curves for image display.
00037 // These splines do not afford local control, and a set of linear equations
00038 // including all interpolation points must be solved before any point on the
00039 // curve can be calculated. The lack of local control makes the splines
00040 // more difficult to handle than e.g. Catmull-Rom splines, and real-time
00041 // interpolation of a stream of data points is not possible.
00042 // For a discussion of natural cubic splines, see e.g. Kreyszig, E., "Advanced
00043 // Engineering Mathematics".
00044 //
00045 // Our approach is to approximate the properties of Catmull-Rom splines for
00046 // piecewice cubic polynomials f(x) = ax^3 + bx^2 + cx + d as follows:
00047 // Each curve segment is specified by four interpolation points,
00048 // p0, p1, p2, p3.
00049 // The curve between p1 and p2 must interpolate both p1 and p2, and in addition
00050 //   f'(p1.x) = k1 = (p2.y - p0.y)/(p2.x - p0.x) and
00051 //   f'(p2.x) = k2 = (p3.y - p1.y)/(p3.x - p1.x).
00052 //
00053 // The constraints are expressed by the following system of linear equations
00054 //
00055 //   [ 1  xi    xi^2    xi^3 ]   [ d ]    [ yi ]
00056 //   [     1  2*xi    3*xi^2 ] * [ c ] =  [ ki ]
00057 //   [ 1  xj    xj^2    xj^3 ]   [ b ]    [ yj ]
00058 //   [     1  2*xj    3*xj^2 ]   [ a ]    [ kj ]
00059 //
00060 // Solving using Gaussian elimination and back substitution, setting
00061 // dy = yj - yi, dx = xj - xi, we get
00062 // 
00063 //   a = ((ki + kj) - 2*dy/dx)/(dx*dx);
00064 //   b = ((kj - ki)/dx - 3*(xi + xj)*a)/2;
00065 //   c = ki - (3*xi*a + 2*b)*xi;
00066 //   d = yi - ((xi*a + b)*xi + c)*xi;
00067 //
00068 // Having calculated the coefficients of the cubic polynomial we have the
00069 // choice of evaluation by brute force
00070 //
00071 //   for (x = x1; x <= x2; x += res) {
00072 //     y = ((a*x + b)*x + c)*x + d;
00073 //     plot(x, y);
00074 //   }
00075 //
00076 // or by forward differencing
00077 //
00078 //   y = ((a*x1 + b)*x1 + c)*x1 + d;
00079 //   dy = (3*a*(x1 + res) + 2*b)*x1*res + ((a*res + b)*res + c)*res;
00080 //   d2y = (6*a*(x1 + res) + 2*b)*res*res;
00081 //   d3y = 6*a*res*res*res;
00082 //     
00083 //   for (x = x1; x <= x2; x += res) {
00084 //     plot(x, y);
00085 //     y += dy; dy += d2y; d2y += d3y;
00086 //   }
00087 //
00088 // See Foley, Van Dam, Feiner, Hughes, "Computer Graphics, Principles and
00089 // Practice" for a discussion of forward differencing.
00090 //
00091 // If we have a set of interpolation points p0, ..., pn, we may specify
00092 // curve segments between p0 and p1, and between pn-1 and pn by using the
00093 // following constraints:
00094 //   f''(p0.x) = 0 and
00095 //   f''(pn.x) = 0.
00096 //
00097 // Substituting the results for a and b in
00098 //
00099 //   2*b + 6*a*xi = 0
00100 //
00101 // we get
00102 //
00103 //   ki = (3*dy/dx - kj)/2;
00104 //
00105 // or by substituting the results for a and b in
00106 //
00107 //   2*b + 6*a*xj = 0
00108 //
00109 // we get
00110 //
00111 //   kj = (3*dy/dx - ki)/2;
00112 //
00113 // Finally, if we have only two interpolation points, the cubic polynomial
00114 // will degenerate to a straight line if we set
00115 //
00116 //   ki = kj = dy/dx;
00117 //
00118 
00119 
00120 #if SPLINE_BRUTE_FORCE
00121 #define interpolate_segment interpolate_brute_force
00122 #else
00123 #define interpolate_segment interpolate_forward_difference
00124 #endif
00125 
00126 
00127 // ----------------------------------------------------------------------------
00128 // Calculation of coefficients.
00129 // ----------------------------------------------------------------------------
00130 inline
00131 void cubic_coefficients(double x1, double y1, double x2, double y2,
00132                         double k1, double k2,
00133                         double& a, double& b, double& c, double& d)
00134 {
00135   double dx = x2 - x1, dy = y2 - y1;
00136 
00137   a = ((k1 + k2) - 2*dy/dx)/(dx*dx);
00138   b = ((k2 - k1)/dx - 3*(x1 + x2)*a)/2;
00139   c = k1 - (3*x1*a + 2*b)*x1;
00140   d = y1 - ((x1*a + b)*x1 + c)*x1;
00141 }
00142 
00143 // ----------------------------------------------------------------------------
00144 // Evaluation of cubic polynomial by brute force.
00145 // ----------------------------------------------------------------------------
00146 template<class PointPlotter>
00147 inline
00148 void interpolate_brute_force(double x1, double y1, double x2, double y2,
00149                              double k1, double k2,
00150                              PointPlotter plot, double res)
00151 {
00152   double a, b, c, d;
00153   cubic_coefficients(x1, y1, x2, y2, k1, k2, a, b, c, d);
00154   
00155   // Calculate each point.
00156   for (double x = x1; x <= x2; x += res) {
00157     double y = ((a*x + b)*x + c)*x + d;
00158     plot(x, y);
00159   }
00160 }
00161 
00162 // ----------------------------------------------------------------------------
00163 // Evaluation of cubic polynomial by forward differencing.
00164 // ----------------------------------------------------------------------------
00165 template<class PointPlotter>
00166 inline
00167 void interpolate_forward_difference(double x1, double y1, double x2, double y2,
00168                                     double k1, double k2,
00169                                     PointPlotter plot, double res)
00170 {
00171   double a, b, c, d;
00172   cubic_coefficients(x1, y1, x2, y2, k1, k2, a, b, c, d);
00173   
00174   double y = ((a*x1 + b)*x1 + c)*x1 + d;
00175   double dy = (3*a*(x1 + res) + 2*b)*x1*res + ((a*res + b)*res + c)*res;
00176   double d2y = (6*a*(x1 + res) + 2*b)*res*res;
00177   double d3y = 6*a*res*res*res;
00178     
00179   // Calculate each point.
00180   for (double x = x1; x <= x2; x += res) {
00181     plot(x, y);
00182     y += dy; dy += d2y; d2y += d3y;
00183   }
00184 }
00185 
00186 template<class PointIter>
00187 inline
00188 double x(PointIter p)
00189 {
00190   return (*p)[0];
00191 }
00192 
00193 template<class PointIter>
00194 inline
00195 double y(PointIter p)
00196 {
00197   return (*p)[1];
00198 }
00199 
00200 // ----------------------------------------------------------------------------
00201 // Evaluation of complete interpolating function.
00202 // Note that since each curve segment is controlled by four points, the
00203 // end points will not be interpolated. If extra control points are not
00204 // desirable, the end points can simply be repeated to ensure interpolation.
00205 // Note also that points of non-differentiability and discontinuity can be
00206 // introduced by repeating points.
00207 // ----------------------------------------------------------------------------
00208 template<class PointIter, class PointPlotter>
00209 inline
00210 void interpolate(PointIter p0, PointIter pn, PointPlotter plot, double res)
00211 {
00212   double k1, k2;
00213 
00214   // Set up points for first curve segment.
00215   PointIter p1 = p0; ++p1;
00216   PointIter p2 = p1; ++p2;
00217   PointIter p3 = p2; ++p3;
00218 
00219   // Draw each curve segment.
00220   for (; p2 != pn; ++p0, ++p1, ++p2, ++p3) {
00221     // p1 and p2 equal; single point.
00222     if (x(p1) == x(p2)) {
00223       continue;
00224     }
00225     // Both end points repeated; straight line.
00226     if (x(p0) == x(p1) && x(p2) == x(p3)) {
00227       k1 = k2 = (y(p2) - y(p1))/(x(p2) - x(p1));
00228     }
00229     // p0 and p1 equal; use f''(x1) = 0.
00230     else if (x(p0) == x(p1)) {
00231       k2 = (y(p3) - y(p1))/(x(p3) - x(p1));
00232       k1 = (3*(y(p2) - y(p1))/(x(p2) - x(p1)) - k2)/2;
00233     }
00234     // p2 and p3 equal; use f''(x2) = 0.
00235     else if (x(p2) == x(p3)) {
00236       k1 = (y(p2) - y(p0))/(x(p2) - x(p0));
00237       k2 = (3*(y(p2) - y(p1))/(x(p2) - x(p1)) - k1)/2;
00238     }
00239     // Normal curve.
00240     else {
00241       k1 = (y(p2) - y(p0))/(x(p2) - x(p0));
00242       k2 = (y(p3) - y(p1))/(x(p3) - x(p1));
00243     }
00244 
00245     interpolate_segment(x(p1), y(p1), x(p2), y(p2), k1, k2, plot, res);
00246   }
00247 }
00248 
00249 // ----------------------------------------------------------------------------
00250 // Class for plotting integers into an array.
00251 // ----------------------------------------------------------------------------
00252 template<class F>
00253 class PointPlotter
00254 {
00255  protected:
00256   F* f;
00257 
00258  public:
00259   PointPlotter(F* arr) : f(arr)
00260   {
00261   }
00262 
00263   void operator ()(double x, double y)
00264   {
00265     // Clamp negative values to zero.
00266     if (y < 0) {
00267       y = 0;
00268     }
00269 
00270     f[int(x)] = F(y + 0.5);
00271   }
00272 };
00273 
00274 } // namespace reSID
00275 
00276 #endif // not RESID_SPLINE_H