public static class Spline
{
    [System.Diagnostics.DebuggerDisplay("({X},{Y})")]
    public partial struct Vec2
    {
        public float X, Y;
        public Vec2(float x, float y) {this.X = x; this.Y = y;}
        public static implicit operator PointF(Vec2 v) { return new PointF(v.X, v.Y); }
        public static implicit operator Vec2(PointF p) { return new Vec2(p.X, p.Y); }
        public static Vec2 operator +(Vec2 v1, Vec2 v2) { return new Vec2(v1.X + v2.X, v1.Y + v2.Y); }
        public static Vec2 operator -(Vec2 v1, Vec2 v2) { return new Vec2(v1.X - v2.X, v1.Y - v2.Y); }
        public static Vec2 operator *(Vec2 v, float f) { return new Vec2(v.X * f, v.Y * f); }
        public static Vec2 operator /(Vec2 v, float f) { return new Vec2(v.X / f, v.Y / f); }
    }

    /// <summary>
    /// '贝塞尔'内插。结果不包括头尾点
    /// </summary>
    public static PointF[] InterpolateBezier(PointF p0, PointF p1, PointF p2, PointF p3, int samples)
    {
        PointF[] result = new PointF[samples];
        for (int i = 0; i < samples; i++)
        {
            float t = (i + 1) / (samples + 1.0f);
            result[i] =
                (Vec2)p0 * (1 - t) * (1 - t) * (1 - t) +
                (Vec2)p1 * (3 * (1 - t) * (1 - t) * t) +
                (Vec2)p2 * (3 * (1 - t) * t * t) +
                (Vec2)p3 * (t * t * t);
        }
        return result;
    }

    public static PointF[] InterpolateCardinalSpline(PointF p0, PointF p1, PointF p2, PointF p3, int samples)
    {
        const float tension = 0.5f;
        Vec2 u = ((Vec2)p2 - (Vec2)p0) * (tension / 3) + p1;
        Vec2 v = ((Vec2)p1 - (Vec2)p3) * (tension / 3) + p2;
        return InterpolateBezier(p1, u, v, p2, samples);
    }

    /// <summary>
    /// '基数样条'内插法。 points为通过点，samplesInSegment为两个样本点之间的内插数量。
    /// </summary>
    public static PointF[] CardinalSpline(PointF[] points, int samplesInSegment)
    {
        List<PointF> result = new List<PointF>();
        for (int i = 0; i < points.Length - 1; i++)
        {
            result.Add(points[i]);
            result.AddRange( InterpolateCardinalSpline(
                points[Math.Max(i-1, 0)],
                points[i],
                points[i+1],
                points[Math.Min(i+2, points.Length-1)],
                samplesInSegment
                ));
        }
        result.Add(points[points.Length - 1]);
        return result.ToArray();
    }
}

------解决方案--------------------
随便搜了一个文章：http://blog.csdn.net/lixiaohuprogram/article/details/4325677

使用多项式方程来插值，假设你有n个点，就可以使用一个n-1阶多项式来插值它。因为n-1阶多项式正好是n个顶点（原点）。所以当你有3个点，那么就可以使用2阶曲线来插值它，得到过这三个点的二阶曲线上每一个点的坐标。

但是假设你有51个点，得到50阶多项式，显然计算量太大。于是有的使用最小二乘法，来拟合出低阶的多项式公式。不过这种方式不适合点数很多的拟合。

另外一种方式，比如有1000个点，每2个点作为一个样条（其实是两个点再加上两个斜率做为参数），保证样条两端的倒数（斜率）与相邻的样条的斜率相同，从而顺滑。

这些计算都非常简单，只不过用基本的矩阵解线性方程组的几个方法就行了。