Draw a smooth curve through a set of 2D points with Cubic Spline

I would like to provide you with the code to draw a smooth curve through a set of 2D points with cubic spline. If we have some tabulated function y_i=f(x_i) it's easy to get its cubic spline interpolant with some library code. For example, you could use the code from "Numerical Recipes in C, 2-nd Edition" book - proved source of a lot of math algorithms. Cubic spline gives an excellent interpolation in the most cases.

Cubic spline is comprised from a sequence of cubic polynomials, so to draw the curve we have to approximate each partial cubic polynomial with the polyline.

Let we have a cubic polynomial defined at [x1, x2] interval.

To approximate it with polyline we should do the following:

1. Get the deviation polynomial, i.e. the difference between the initial cubic polynomial and the straight line passing through its left and right bound points. This polynomial is either identically equal to zero or has one or two extremum(s) at [x1, x2].
2. Evaluate the values of deviation polynomial at extremum points. It its absolute values are lower than the tolerance then the initial cubic polynomial can be approximated with a straight line passing through points (x1, y1) and (x2, y2). Otherwise
3. Split the initial interval  [x1, x2] on two or three subintervals (depending on extremum count) and repeat the procedure recursively from (1) for each of subintervals.

    /// 
    /// Approximating Cubic Polynomial with PolyLine.
    /// 
    public static class CubicPolynomialPolylineApproximation
    {
        /// 
        /// Gets the approximation of the polynomial with polyline.
        /// 
        /// The polynomial.
        /// The abscissas start.
        /// The abscissas stop.
        /// The tolerance is the maximum distance from the cubic 
        /// polynomial to the approximating polyline.
        /// 
        public static Collection Approximate(Polynomial polynomial, double x1, double x2, double tolerance)
        {
            Debug.Assert(x1 <= x2, "x1 <= x2");
            Debug.Assert(polynomial.Order == 3, "polynomial.Order == 3");
 
            Collection points = new Collection();
 
            // Get difference between given polynomial and the straight line passing its node points.
            Polynomial deviation = DeviationPolynomial(polynomial, x1, x2);
            Debug.Assert(deviation.Order == 3, "diff.Order == 3");
            if (deviation[0] == 0 && deviation[1] == 0 && deviation[2] == 0 && deviation[3] == 0)
            {
                points.Add(new Point(x1, polynomial.GetValue(x1)));
                points.Add(new Point(x2, polynomial.GetValue(x2)));
                return points;
            }
 
            // Get previouse polynomial first derivative
            Polynomial firstDerivative = new Polynomial(new double[] { deviation[1], 2 * deviation[2], 3 * deviation[3] });
 
            // Difference polinomial extremums.
            // Fing first derivative roots.
            Complex[] complexRoots = firstDerivative.Solve();
            // Get real roots in [x1, x2].
            List roots = new List();
            foreach (Complex complexRoot in complexRoots)
            {
                if (complexRoot.Imaginary == 0)
                {
                    double r = complexRoot.Real;
                    if (r > x1 && r < x2)
                        roots.Add(r);
                }
            }
            Debug.Assert(roots.Count > 0, "roots.Count > 0");
            Debug.Assert(roots.Count <= 2, "roots.Count <= 2");
 
            // Check difference polynomial extremal values.
            bool approximates = true;
            foreach (double x in roots)
            {
                if (Math.Abs(deviation.GetValue(x)) > tolerance)
                {
                    approximates = false;
                    break;
                }
            }
            if (approximates)
            {// Approximation is good enough.
                points.Add(new Point(x1, polynomial.GetValue(x1)));
                points.Add(new Point(x2, polynomial.GetValue(x2)));
                return points;
            }
 
            if (roots.Count == 2)
            {
                if (roots[0] == roots[1])
                    roots.RemoveAt(1);
                else if (roots[0] > roots[1])
                {// Sort the roots
                    // Swap roots
                    double x = roots[0];
                    roots[0] = roots[1];
                    roots[1] = x;
                }
            }
            // Add the end abscissas.
            roots.Add(x2);
 
            // First subinterval.
            Collection pts = Approximate(polynomial, x1, roots[0], tolerance);
            // Copy all points.
            foreach (Point pt in pts)
            {
                points.Add(pt);
            }
            // The remnant of subintervals.
            for (int i = 0; i < roots.Count - 1; ++i)
            {
                pts = Approximate(polynomial, roots[i], roots[i + 1], tolerance);
                // Copy all points but the first one.
                for (int j = 1; j < pts.Count; ++j)
                {
                    points.Add(pts[j]);
                }
            }
            return points;
        }
 
        /// 
        /// Gets the difference between given polynomial and the straight line passing through its node points.
        /// 
        /// The polynomial.
        /// The abscissas start.
        /// The abscissas stop.
        /// 
        static Polynomial DeviationPolynomial(Polynomial polynomial, double x1, double x2)
        {
            double y1 = polynomial.GetValue(x1);
            double y2 = polynomial.GetValue(x2);
            double a = (y2 - y1) / (x2 - x1);
            double b = y1 - a * x1;
            if (a != 0)
                return polynomial.Subtract(new Polynomial(new double[] { b, a }));
            else if (b != 0)
                return polynomial.Subtract(new Polynomial(new double[] { b }));
            else
                return polynomial;
        }
    }

In the code above I'm using the helper class Polynomial encapsulating operations on polynomials including addition, subtraction, dividing, root finding, etc. It's ported from "Numerical Recipes in C, 2-nd Edition" book with some additions and bug fixes.

The sample supplied with this article is Visual Studio 2008 solution targeted to .NET 3.5. It contains WPF Windows Application project designed to demonstrate some curves drawn with cubic spline. You can select one of the curves from Combo Box at the top of the Window, experiment with point counts, tolerance and set appropriate XY Scales. You can even add you own curve, but this requires coding as follows:

1. Add your curve name to CurveNames enum.
2. Add your curve implementation to Curves region.
   Add call to your curve to OnRender override.
3. In the sample I use Path elements on the custom Canvas to render the curve but in real application you would probably use some more effective approach like visual layer rendering.