For one thing: That is not a spline interpolation, that is just a Bézier curve. Those are simply polynomial in specific bases.
A spline interpolation would consist of piecewise smooth functions so that the whole function has smooth differentiability to some degree.
So what you actually want: Find a cubic Bézier curve that describes your points optimally. But then, it does not matter that much if you have Bézier curves or just a regular monomial base.
What you can then do is simply regression:
Take some $t_1,\ldots,t_N$ and the model
$$ x(t) = a+bt+ct^2+dt^3 $$
$$ y(t) = e+ft+gt^2+ht^3 $$
or for Bézier representation
$$ x(t) = aB_{0,3}+bB_{1,3}+cB_{2,3}+dB_{3,3} $$
$$ y(t) = eB_{0,3}+fB_{1,3}+gB_{2,3}+hB_{3,3} $$
Then consider the quadratic error
$$ E(a,b,c,d,e,f,g,h,t_1,\ldots,t_N) = \sum_{j=1}^N (x(t_j)-P_x^j)^2 + (y(t_j)-P_y^j)^2 $$
Then differentiate and equate to $0$ (easy if $t_1,\ldots,t_N$ are fixed, harder if they are variable as that will then involve finding multivariate roots, which might be solvable using Groebner bases) to find the paramters that minimize this error.
You can then take the minimum and the maximum of the $t_i$ and transform $x(t),y(t)$ to $[0,1]$.