So, I got pretty interested in this topic b/c I was looking to do something similar.

So, it turns out that Interpolating the x and y components separately actually made the calculation simpler because I could choose which values of T to use.

So given a data set of at least four points (x1,y1,) (x2,y2) (x3,y3) (x4,y4), we separate it

(0,x1) (1,x2) (2,x3) (3,x4) and (0,y1) (0,y2) (0,y3) and (0,y4)

The function created with four points can always be represented as F(x)=ax^3+bx^2+cx+d, so we can store the function as a set of four values (a,b,c,d).

To get a value at T=t

double Px(Tuple<double, double, double, double> coefficents, double t)
{
double ret = coefficents.Item1 * Math.Pow(t, 3);
ret += coefficents.Item2 * Math.Pow(t, 2);
ret += coefficents.Item3 * t;
ret += coefficents.Item4;
return ret;
}

so the full polynomail is

y1*1/((x1-x2)(x1-x3)(x1-x4))*(x-x2)(x-x3)(x-x4)

+ y2*1/((x2-x1)(x2-x3)(x2-x4))*(x-x1)(x-x3)(x-x4)

+ y3*1/((x3-x1)(x3-x2)(x3-x4))*(x-x1)(x-x2)(x-x4)

+ y4*1/((x4-x1)(x4-x2)(x4-x3))*(x-x1)(x-x2)(x-x3)

But, we already choose the X(T) values (0,1,2,3) so we can calculate that part and put it as a constant

y1*-1/6*(x-1)(x-2)(x-3)

+ y2*1/2*(x-0)(x-2)(x-3)

+ y3*-1/2*(x-0)(x-1)(x-3)

+ y4*1/6*(x-0)(x-1)(x-2)

And we can expand to get

y1* -1/6 *( x^3-6x^2+11x-6)

+y2* 1/2*(x^3-5x+6x)

+y3* -1/2*(x^3-4x^2+3x

+y4* 1/6*(x^3-3x^2+2x)

Then distribute to get

-(1/6)y1x^3+y1x^2 -(11/6)y1x+1y1

+ (y2/2)x^3-(5/2)y2*x^2+3y2x

+ -(y3/2)x^3+2y3x^2-(3/2)y3x

+ (y4/6)x^3-(y4/2)x^2+(y4/3)x

Now we can collect the terms and it leaves us with the following (v1,v2,v3, and v4 are the points to parameterize by t):

Tuple<double, double, double, double> GetCoeffInT(double v1, double v2, double v3, double v4)
{
double a = (-v1 / 6.0) + (v2 / 2.0) + (-v3 / 2.0) + (v4 / 6.0);
double b = v1 - (5.0 / 2.0) * v2 + 2 * v3 - v4 / 2.0;
double c = (-11.0 / 6.0) * v1 + 3.0 * v2 - (3.0 / 2.0) * v3 + v4 / 3.0;
double d = v1;
return new Tuple<double, double, double, double>(a, b, c, d);
}

Now we take our four X/Y cords and convert them to two sets of co-efficents

Tuple<Tuple<double, double, double, double>,Tuple<double, double, double, double>> GetCoefficents(Point p1, Point p2, Point p3, Point p4)
{
Tuple<double, double, double, double> xCoeff = GetCoeffInT(p1.X, p2.X, p3.X, p4.X);
Tuple<double,double,double,double> yCoeff=GetCoeffInT(p1.Y,p2.Y,p3.Y,p4.Y);
return new Tuple<Tuple<double, double, double, double>, Tuple<double, double, double, double>>(xCoeff, yCoeff);
}

The last part is to draw a line of n>4 data points

double step = .01;
for (int i = 3;i < DataPoints.Count; i++)
{
int ii = i;
var v = GetCoefficents(DataPoints[ii-3], DataPoints[ii-2], DataPoints[ii-1], DataPoints[ii]);
//On the first segment draw from T=0 to T=2
//On the last segment draw from T=1 to T=3
//On every other segment draw from T=1 to T=2
for (double t = (i==3?0:1); t <= (i==DataPoints.Count-1?3.0:2.0); t += step)
{
Point p1 = new Point((int)Px(v.Item1, t - step), (int)Px(v.Item2, t - step));
Point p2 = new Point((int)Px(v.Item1, t), (int)Px(v.Item2, t));
g.DrawLine(Pens.Green, p1, p2);
}
}

Even more interesting is that is it trivial to expand this to any number of dimentions.