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