home *** CD-ROM | disk | FTP | other *** search
/ Large Pack of OldSkool DOS MOD Trackers / goattracker_2.73.zip / src / resid / spline.h < prev    next >
Encoding:
C/C++ Source or Header  |  2014-07-23  |  8.8 KB  |  273 lines
  1. //  ---------------------------------------------------------------------------
  2. //  This file is part of reSID, a MOS6581 SID emulator engine.
  3. //  Copyright (C) 2004  Dag Lem <resid@nimrod.no>
  4. //
  5. //  This program is free software; you can redistribute it and/or modify
  6. //  it under the terms of the GNU General Public License as published by
  7. //  the Free Software Foundation; either version 2 of the License, or
  8. //  (at your option) any later version.
  9. //
  10. //  This program is distributed in the hope that it will be useful,
  11. //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  12. //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13. //  GNU General Public License for more details.
  14. //
  15. //  You should have received a copy of the GNU General Public License
  16. //  along with this program; if not, write to the Free Software
  17. //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  18. //  ---------------------------------------------------------------------------
  19.  
  20. #ifndef __SPLINE_H__
  21. #define __SPLINE_H__
  22.  
  23. // Our objective is to construct a smooth interpolating single-valued function
  24. // y = f(x).
  25. //
  26. // Catmull-Rom splines are widely used for interpolation, however these are
  27. // parametric curves [x(t) y(t) ...] and can not be used to directly calculate
  28. // y = f(x).
  29. // For a discussion of Catmull-Rom splines see Catmull, E., and R. Rom,
  30. // "A Class of Local Interpolating Splines", Computer Aided Geometric Design.
  31. //
  32. // Natural cubic splines are single-valued functions, and have been used in
  33. // several applications e.g. to specify gamma curves for image display.
  34. // These splines do not afford local control, and a set of linear equations
  35. // including all interpolation points must be solved before any point on the
  36. // curve can be calculated. The lack of local control makes the splines
  37. // more difficult to handle than e.g. Catmull-Rom splines, and real-time
  38. // interpolation of a stream of data points is not possible.
  39. // For a discussion of natural cubic splines, see e.g. Kreyszig, E., "Advanced
  40. // Engineering Mathematics".
  41. //
  42. // Our approach is to approximate the properties of Catmull-Rom splines for
  43. // piecewice cubic polynomials f(x) = ax^3 + bx^2 + cx + d as follows:
  44. // Each curve segment is specified by four interpolation points,
  45. // p0, p1, p2, p3.
  46. // The curve between p1 and p2 must interpolate both p1 and p2, and in addition
  47. //   f'(p1.x) = k1 = (p2.y - p0.y)/(p2.x - p0.x) and
  48. //   f'(p2.x) = k2 = (p3.y - p1.y)/(p3.x - p1.x).
  49. //
  50. // The constraints are expressed by the following system of linear equations
  51. //
  52. //   [ 1  xi    xi^2    xi^3 ]   [ d ]    [ yi ]
  53. //   [     1  2*xi    3*xi^2 ] * [ c ] =  [ ki ]
  54. //   [ 1  xj    xj^2    xj^3 ]   [ b ]    [ yj ]
  55. //   [     1  2*xj    3*xj^2 ]   [ a ]    [ kj ]
  56. //
  57. // Solving using Gaussian elimination and back substitution, setting
  58. // dy = yj - yi, dx = xj - xi, we get
  59. // 
  60. //   a = ((ki + kj) - 2*dy/dx)/(dx*dx);
  61. //   b = ((kj - ki)/dx - 3*(xi + xj)*a)/2;
  62. //   c = ki - (3*xi*a + 2*b)*xi;
  63. //   d = yi - ((xi*a + b)*xi + c)*xi;
  64. //
  65. // Having calculated the coefficients of the cubic polynomial we have the
  66. // choice of evaluation by brute force
  67. //
  68. //   for (x = x1; x <= x2; x += res) {
  69. //     y = ((a*x + b)*x + c)*x + d;
  70. //     plot(x, y);
  71. //   }
  72. //
  73. // or by forward differencing
  74. //
  75. //   y = ((a*x1 + b)*x1 + c)*x1 + d;
  76. //   dy = (3*a*(x1 + res) + 2*b)*x1*res + ((a*res + b)*res + c)*res;
  77. //   d2y = (6*a*(x1 + res) + 2*b)*res*res;
  78. //   d3y = 6*a*res*res*res;
  79. //     
  80. //   for (x = x1; x <= x2; x += res) {
  81. //     plot(x, y);
  82. //     y += dy; dy += d2y; d2y += d3y;
  83. //   }
  84. //
  85. // See Foley, Van Dam, Feiner, Hughes, "Computer Graphics, Principles and
  86. // Practice" for a discussion of forward differencing.
  87. //
  88. // If we have a set of interpolation points p0, ..., pn, we may specify
  89. // curve segments between p0 and p1, and between pn-1 and pn by using the
  90. // following constraints:
  91. //   f''(p0.x) = 0 and
  92. //   f''(pn.x) = 0.
  93. //
  94. // Substituting the results for a and b in
  95. //
  96. //   2*b + 6*a*xi = 0
  97. //
  98. // we get
  99. //
  100. //   ki = (3*dy/dx - kj)/2;
  101. //
  102. // or by substituting the results for a and b in
  103. //
  104. //   2*b + 6*a*xj = 0
  105. //
  106. // we get
  107. //
  108. //   kj = (3*dy/dx - ki)/2;
  109. //
  110. // Finally, if we have only two interpolation points, the cubic polynomial
  111. // will degenerate to a straight line if we set
  112. //
  113. //   ki = kj = dy/dx;
  114. //
  115.  
  116.  
  117. #if SPLINE_BRUTE_FORCE
  118. #define interpolate_segment interpolate_brute_force
  119. #else
  120. #define interpolate_segment interpolate_forward_difference
  121. #endif
  122.  
  123.  
  124. // ----------------------------------------------------------------------------
  125. // Calculation of coefficients.
  126. // ----------------------------------------------------------------------------
  127. inline
  128. void cubic_coefficients(double x1, double y1, double x2, double y2,
  129.       double k1, double k2,
  130.       double& a, double& b, double& c, double& d)
  131. {
  132.   double dx = x2 - x1, dy = y2 - y1;
  133.  
  134.   a = ((k1 + k2) - 2*dy/dx)/(dx*dx);
  135.   b = ((k2 - k1)/dx - 3*(x1 + x2)*a)/2;
  136.   c = k1 - (3*x1*a + 2*b)*x1;
  137.   d = y1 - ((x1*a + b)*x1 + c)*x1;
  138. }
  139.  
  140. // ----------------------------------------------------------------------------
  141. // Evaluation of cubic polynomial by brute force.
  142. // ----------------------------------------------------------------------------
  143. template<class PointPlotter>
  144. inline
  145. void interpolate_brute_force(double x1, double y1, double x2, double y2,
  146.            double k1, double k2,
  147.            PointPlotter plot, double res)
  148. {
  149.   double a, b, c, d;
  150.   cubic_coefficients(x1, y1, x2, y2, k1, k2, a, b, c, d);
  151.   
  152.   // Calculate each point.
  153.   for (double x = x1; x <= x2; x += res) {
  154.     double y = ((a*x + b)*x + c)*x + d;
  155.     plot(x, y);
  156.   }
  157. }
  158.  
  159. // ----------------------------------------------------------------------------
  160. // Evaluation of cubic polynomial by forward differencing.
  161. // ----------------------------------------------------------------------------
  162. template<class PointPlotter>
  163. inline
  164. void interpolate_forward_difference(double x1, double y1, double x2, double y2,
  165.             double k1, double k2,
  166.             PointPlotter plot, double res)
  167. {
  168.   double a, b, c, d;
  169.   cubic_coefficients(x1, y1, x2, y2, k1, k2, a, b, c, d);
  170.   
  171.   double y = ((a*x1 + b)*x1 + c)*x1 + d;
  172.   double dy = (3*a*(x1 + res) + 2*b)*x1*res + ((a*res + b)*res + c)*res;
  173.   double d2y = (6*a*(x1 + res) + 2*b)*res*res;
  174.   double d3y = 6*a*res*res*res;
  175.     
  176.   // Calculate each point.
  177.   for (double x = x1; x <= x2; x += res) {
  178.     plot(x, y);
  179.     y += dy; dy += d2y; d2y += d3y;
  180.   }
  181. }
  182.  
  183. template<class PointIter>
  184. inline
  185. double x(PointIter p)
  186. {
  187.   return (*p)[0];
  188. }
  189.  
  190. template<class PointIter>
  191. inline
  192. double y(PointIter p)
  193. {
  194.   return (*p)[1];
  195. }
  196.  
  197. // ----------------------------------------------------------------------------
  198. // Evaluation of complete interpolating function.
  199. // Note that since each curve segment is controlled by four points, the
  200. // end points will not be interpolated. If extra control points are not
  201. // desirable, the end points can simply be repeated to ensure interpolation.
  202. // Note also that points of non-differentiability and discontinuity can be
  203. // introduced by repeating points.
  204. // ----------------------------------------------------------------------------
  205. template<class PointIter, class PointPlotter>
  206. inline
  207. void interpolate(PointIter p0, PointIter pn, PointPlotter plot, double res)
  208. {
  209.   double k1, k2;
  210.  
  211.   // Set up points for first curve segment.
  212.   PointIter p1 = p0; ++p1;
  213.   PointIter p2 = p1; ++p2;
  214.   PointIter p3 = p2; ++p3;
  215.  
  216.   // Draw each curve segment.
  217.   for (; p2 != pn; ++p0, ++p1, ++p2, ++p3) {
  218.     // p1 and p2 equal; single point.
  219.     if (x(p1) == x(p2)) {
  220.       continue;
  221.     }
  222.     // Both end points repeated; straight line.
  223.     if (x(p0) == x(p1) && x(p2) == x(p3)) {
  224.       k1 = k2 = (y(p2) - y(p1))/(x(p2) - x(p1));
  225.     }
  226.     // p0 and p1 equal; use f''(x1) = 0.
  227.     else if (x(p0) == x(p1)) {
  228.       k2 = (y(p3) - y(p1))/(x(p3) - x(p1));
  229.       k1 = (3*(y(p2) - y(p1))/(x(p2) - x(p1)) - k2)/2;
  230.     }
  231.     // p2 and p3 equal; use f''(x2) = 0.
  232.     else if (x(p2) == x(p3)) {
  233.       k1 = (y(p2) - y(p0))/(x(p2) - x(p0));
  234.       k2 = (3*(y(p2) - y(p1))/(x(p2) - x(p1)) - k1)/2;
  235.     }
  236.     // Normal curve.
  237.     else {
  238.       k1 = (y(p2) - y(p0))/(x(p2) - x(p0));
  239.       k2 = (y(p3) - y(p1))/(x(p3) - x(p1));
  240.     }
  241.  
  242.     interpolate_segment(x(p1), y(p1), x(p2), y(p2), k1, k2, plot, res);
  243.   }
  244. }
  245.  
  246. // ----------------------------------------------------------------------------
  247. // Class for plotting integers into an array.
  248. // ----------------------------------------------------------------------------
  249. template<class F>
  250. class PointPlotter
  251. {
  252.  protected:
  253.   F* f;
  254.  
  255.  public:
  256.   PointPlotter(F* arr) : f(arr)
  257.   {
  258.   }
  259.  
  260.   void operator ()(double x, double y)
  261.   {
  262.     // Clamp negative values to zero.
  263.     if (y < 0) {
  264.       y = 0;
  265.     }
  266.  
  267.     f[F(x)] = F(y);
  268.   }
  269. };
  270.  
  271.  
  272. #endif // not __SPLINE_H__
  273.