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