# Arbitrary precision spline interpolation

The current implementation of Interpolation does not allow arbitrary precision spline interpolation. Yu-Sung Chang says here that "it is not hard to implement it manually for arbitrary precision using BSplineBasis."

I am new to splines and not familiar with BSplineBasis. What is the best way to implement arbitrary precision spline interpolation in Mathematica?

I am interested primarily in an analog of the Method -> "Spline" of Interpolation which I investigated in this answer (BTW, what is the name of such a spline and of this kind of parametrization?). Also, J. M. gives in this answer an interesting implementation of a spline with centripetal parametrization, but it is not clear how to use it for interpolation.

P.S. I need an implementation where the degree of a spline is arbitrary. At the moment I am more interested in even-degree splines.

-
Yu-Sung's pointer is wolfram.com/xid/0c0rpcn4bku6m-dj3q8d, where there are two examples of such a thing (I haven't tested them). What are those example lacking? –  belisarius is forth Jan 10 '14 at 12:45
That is an example of cubic B-spline interpolation. I need an implementation where the degree of spline is arbitrary. –  Alexey Popkov Jan 10 '14 at 12:51
Please see my question about B-spline Basis Function –  Shutao TANG May 28 at 8:12

The OP linked to an answer of mine for interpolating over general point sets; for constructing a single interpolating function, a slight modification of my procedure is needed. (In particular, you don't need centripetal or chord-length parametrization in this case.)

data = {{0, 0}, {1/10, 3/10}, {1/2, 3/5}, {1, -1/5}, {2, 3}, {3, -6/5}};


To review the problems with using the built-in Interpolation[], let's try building an interpolating quartic spline from data:

m = 4; (* spline order *)
sp1 = Interpolation[data, InterpolationOrder -> m, Method -> "Spline"];


Now, evaluate:

sp1[3/2]
0.26753591659167364


Herein lies the problem: even though the contents of data and the argument fed to sp1 are exact, the output is a machine precision number. Yeesh!

So, we go back to basics and build the interpolating spline ourselves from BSplineBasis[]. First, we generate the knots. The clamped condition is the most appropriate for interpolation, so we generate the appropriate knot sequence like so:

n = Length[data];
{xa, ya} = Transpose[data]
knots = Join[ConstantArray[xa[[1]], m + 1],
If[m + 2 <= n, MovingAverage[ArrayPad[xa, -1], m], {}],
ConstantArray[xa[[-1]], m + 1]];


To generate the control points for the B-spline, we construct the relevant linear system and solve it with LinearSolve[]:

cp = LinearSolve[Outer[BSplineBasis[{m, knots}, #2, #1] &,
xa, Range[0, Length[data] - 1]], ya];


Note that all the entries of cp are exact.

At this juncture, you might think that we can use BSplineFunction[] on the control points and the knots, but alas, the function also only evaluates at machine precision. We are left with no choice but to use BSplineBasis[] once more:

sp2[x_] = cp.Table[BSplineBasis[{m, knots}, j - 1, x], {j, n}];


Try it out:

sp2[3/2]
25251843/94386740

N[%]
0.2675359165916738


As a verification that we were able to reproduce the control points used internally by Interpolation[], let's do a comparison:

sp1bs = Cases[sp1, _BSplineFunction, ∞][[1]];
sp1bs["ControlPoints"]
{0., 0.777830630128766, 1.3470257035045392, -6.670028631670085, 13.100332135107115, -1.2}

N[cp]
{0., 0.777830630128766, 1.347025703504539, -6.6700286316700845, 13.100332135107113, -1.2}


Certainly, I can't end this post without at least one picture, so:

Plot[{sp1[x], sp2[x]}, {x, 0, 3}]


Plot[sp1[x] - sp2[x], {x, 0, 3}, PlotRange -> All]


I leave the task of bundling everything in a routine as an exercise. If you want to learn more, Piegl and Tiller's The NURBS Book is the canonical reference.

-
It is good to have you back. –  Mr.Wizard May 24 at 10:31
I was able to borrow a friend's computer with a Mathematica installation, so I went for it and answered some questions. It was fun while it lasted… :) –  J. M. is back. May 24 at 10:45
Did you try contacting the Stack Exchange team about the matter? I understand if you do not wish to pursue it but I think they might help. –  Mr.Wizard May 24 at 10:48
Maybe I'll ask them sometime… thanks for the concern. :) –  J. M. is back. May 24 at 10:51
@Guess who it is. May I ask you also about a generalization of the spline interpolation (at least) for the 2D case? In the linked answer it is said that "Multi-dimension is simply tensor product version." I do not see this as a simple thing but it would be very helpful to have such a generalization. May be a separate thread would be the best place for this. –  Alexey Popkov May 26 at 14:01

Here is my implementation of quadric spline arbitrary precision interpolation which is defined exactly as in Interpolation with options Method->"Spline", InterpolationOrder->2.

## Theoretical background

Quadric spline interpolation for n datapoints is defined as a Piecewise function which consists of n-2 parabolas which are spliced together in the middles of successive datapoints (with exceptions for first and last intervals) by two conditions: successive parabolas at these middle points must be equal and their first derivatives must be equal too. So we have 2(n-3) such conditions.

Let us designate jth datapoint as {x[j],y[j]} and define jth parabola as

s[j_, xx_] = y[j + 1] + b[j] (xx - x[j + 1]) + c[j] (xx - x[j + 1])^2;


This definition automatically makes jth parabola equal to y[j + 1] at x[j + 1]. We have to add also 2 boundary conditions:

s[1, x[1]] == y[1]
s[n - 2, x[n]] == y[n]


And 2(n-3) splicing conditions:

Table[{
s[j, (x[j + 1] + x[j + 2])/2] == s[j + 1, (x[j + 1] + x[j + 2])/2],
ds[j, (x[j + 1] + x[j + 2])/2] == ds[j + 1, (x[j + 1] + x[j + 2])/2]
}, {j, 1, n - 3}]


where ds is first derivative:

ds[j_, xx_] = D[s[j, xx], xx];


Now we combine everything together and convert into matrix form:

eqs = Flatten[{
s[1, x[1]] == y[1],
Table[{
s[j, (x[j + 1] + x[j + 2])/2] == s[j + 1, (x[j + 1] + x[j + 2])/2],
ds[j, (x[j + 1] + x[j + 2])/2] == ds[j + 1, (x[j + 1] + x[j + 2])/2]
}, {j, 1, n - 3}],
s[n - 2, x[n]] == y[n]}];
arrs = Simplify@
Normal@CoefficientArrays[eqs, Flatten@Array[{b[#], c[#]} &, n - 2]];
MatrixForm /@ (4 arrs)


Looking at the matrices it is easy to see that they have periodic structure and contain only elements of the forms x[j+1]-x[j] and y[j+1]-y[j] with some numerical coefficients. This periodic structures can be expressed as SparceArrays.

## Implementation

Assuming that data contains an array of {x[j],y[j]}, we define

Δx[i_] := Subtract @@ data[[{i + 1, i}, 1]];
Δy[i_] := Subtract @@ data[[{i + 1, i}, 2]];


Two matrices in the arrs can be expressed as follows:

bB = SparseArray[{1 -> -4 Δy[1], -1 -> 4 Δy[n - 1], i_ /; OddQ[i] :> 0, i_ /; EvenQ[i] :> 4 Δy[i/2 + 1]}, 2 (n - 2)];
mM = SparseArray[
Join[{{1, 1} -> -4 Δx[1], {1, 2} -> 4 Δx[1]^2, {-1, -1} -> 4 Δx[n - 1]^2, {-1, -2} -> 4 Δx[n - 1]},
Table[Band[{2 i, 2 i - 1}] -> {{2 Δx[i + 1], Δx[i + 1]^2, 2 Δx[i + 1], -Δx[i + 1]^2}, {4, 4 Δx[i + 1], -4, 4 Δx[i + 1]}}, {i, 1, n - 3}]], {2 (n - 2), 2 (n - 2)}];


Now for obtaining coefficients b[i], c[i] we can use LinearSolve:

bcM = LinearSolve[mM, bB];


Grouping coefficients with identical indexes together:

bcM = Partition[bcM, 2];


The intervals at which the parabolas are defined:

intervList =
Partition[
Join[{data[[1, 1]]},
MovingAverage[data[[2 ;; -2, 1]], 2], {data[[-1, 1]]}], 2, 1];


Now we compile all the spline data in one object:

splineData = Table[{data[[i + 1]], bcM[[i]], intervList[[i]]}, {i, n - 2}];


splineData contains everything that is needed to construct the spline. It contains elements of the form:

{{x[i + 1], y[i + 1]}, {b[i], c[i]}, {xmin, xmax}}


It is very handy to use splineData for exact numerical integration (see here).

To produce the explicit piecewise function we define the constructor (HornerForm is already applied and gives 27% speedup):

makeSpline[splineData_List, x_Symbol] :=
Piecewise[
Append[Table[{d[[1, 2]] - d[[1, 1]] d[[2, 1]] +
d[[1, 1]]^2 d[[2, 2]] + x (d[[2, 1]] + x d[[2, 2]] -
2 d[[1, 1]] d[[2, 2]]), #1 <= x <= #2 & @@ d[[3]]}, {d, splineData}],
{Indeterminate, True}]]


It can be used as follows:

spline[\[FormalX]_] = makeSpline[splineData, \[FormalX]];


Compiling the code for creation of splineData into one function:

ClearAll[toSplineData, makeSpline];
Options[toSplineData] = {Method -> Automatic};
toSplineData[data_, OptionsPattern[]] /; MatrixQ[data, NumberQ] :=
Module[{Δx, Δy, bB, mM, bcM, intervList,
n = Length[data], dataS = Sort[data]},
Δx[i_] :=
Subtract @@ dataS[[{i + 1, i}, 1]]; Δy[i_] :=
Subtract @@ dataS[[{i + 1, i}, 2]];
bB = SparseArray[{1 -> -4 Δy[1], -1 ->
4 Δy[n - 1], i_ /; OddQ[i] :> 0,
i_ /; EvenQ[i] :> 4 Δy[i/2 + 1]}, 2 (n - 2)];
mM = SparseArray[
Join[{{1, 1} -> -4 Δx[1], {1, 2} ->
4 Δx[1]^2, {-1, -1} ->
4 Δx[n - 1]^2, {-1, -2} ->
4 Δx[n - 1]},
Table[Band[{2 i,
2 i - 1}] -> {{2 Δx[i + 1], Δx[
i + 1]^2,
2 Δx[i + 1], -Δx[i + 1]^2}, {4,
4 Δx[i + 1], -4,
4 Δx[i + 1]}}, {i, 1, n - 3}]], {2 (n - 2),
2 (n - 2)}];
bcM = LinearSolve[mM, bB, Method -> OptionValue[Method]];
bcM = Partition[bcM, 2];
intervList =
Partition[
Join[{dataS[[1, 1]]},
MovingAverage[dataS[[2 ;; -2, 1]], 2], {dataS[[-1, 1]]}], 2, 1];
Table[{dataS[[i + 1]], bcM[[i]], intervList[[i]]}, {i, n - 2}]];


Now to construct the spline from data we can evaluate

spline[\[FormalX]_] = makeSpline[toSplineData[data, \[FormalX]]];


Any suggestions and improvements are welcome!

Quadric spline interpolation with splicing points in the middle of successive data points is locally invariant and weakly depends on boundary condition. (But the same is not true for quadric spline interpolation with splicings in data points.) It gives much lesser artifacts than "usual" cubic spline interpolation. Moreover, it seems that splines of higher degree do not give any advantages over such quadric splines. As opposed to cubic spline, data points are not special points at which polynomials are spliced, and it makes much easier exact integration in these points.

-

The following is a shameful plug of J. M.'s answer that you already linked to show how to get the parametrization by using BSplineBasis[] with arbitrary precision and packed into a function:

splineInterp[data_, order_, prec_] :=
Module[{parametrizeCurve, tvals, bas, ctrlpts, knots},
parametrizeCurve[pts_ /; MatrixQ[pts, NumericQ], a : (_?NumericQ) : 1/2] :=
FoldList[
Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]];
tvals = parametrizeCurve[data];
(*knots for interpolating B-spline*)
knots =
Join[
ConstantArray[0, order + 1],
ConstantArray[1, order + 1]];
(*basis function matrix*)
bas =
Table[
BSplineBasis[{order, knots}, j - 1, N[tvals[[i]], prec]],
{i, Length[data]}, {j, Length[data]}];
ctrlpts = LinearSolve[bas, data];
Return[
Sum[
ctrlpts[[i + 1]] BSplineBasis[{order, knots}, i, #],
{i, 0, Length[testData] - 1}] &]
]


Usage

SeedRandom[42];
testData = RandomReal[{0, 1}, {10, 2}, WorkingPrecision -> prec];
f = splineInterp[testData, 5, 50];
f[1/10]
ParametricPlot[
f[t], {t, 0, 1},
Axes -> None,
Epilog ->
{Directive[Green, PointSize[Large]], Point[testData]}, Frame -> True
]

(*
{0.06025513038208479326777055464103084790562316, \
0.5927955887343273319037301352154452480202510}
*)


### Explanation

The centripetal distribution:

The function parametrizeCurve complete this algorithm

parametrizeCurve[pts_ /; MatrixQ[pts, NumericQ], a : (_?NumericQ) : 1/2] :=
FoldList[
Plus, 0, Normalize[(Norm /@ Differences[pts])^a, Total]]

-
I do not see in your answer how can I use this as interpolation in a regular way: interpolation[x] where x is not a parametrization variable t but the actual independent variable on which the dependent variable y depends by y[x]. –  Alexey Popkov Jan 10 '14 at 16:28
@AlexeyPopkov, Please see my explanation –  Shutao TANG May 29 at 13:04
@ShutaoTang Thank you for the answer to my comment, now I see that the second argument of parametrizeCurve may be arbitrary but must be lesser or equal to $1$. BTW it would be better for others if you would post this explanation under your own answer. :^) –  Alexey Popkov May 29 at 13:22
@AlexeyPopkov, My pleasue. I used the chord-length method to generate the parameters in my answer, so I think it is suitable to place this explanation in belisarius's answer:-) –  Shutao TANG May 29 at 13:34

A bug

### Implementation

Search the index of span

In generally, when $u_0 \in \left[u _i ,u _{i+1} \right)$, the B-spline function basis $N_{i-p,p},\cdots ,N_{i,p}$ are not equal to 0.

searchSpan[knots_, u0_] :=
With[{max = Max[knots]},
If[u0 == max,
Position[knots, max][[1, 1]] - 2,
Ordering[UnitStep[u0 - knots], 1][[1]] - 2]
]


Coefficient Matrix

Options[BSplineInterpolation] =
{SplineDegree -> Automatic, Method -> Automatic, SplineColor -> Black};

BSplineInterpolation[pts : {{_, _} ..},opts : OptionsPattern[]] /;
MatrixQ[pts, NumericQ] :=
Module[{n, md, sd, preKnots, knots, coeffMat, ctrlpts},
n = Length@pts - 1;
md = OptionValue[Method];
sd = OptionValue[SplineDegree] /. Automatic -> 3 /.
deg_ :> n /; deg > n;
preKnots =
Which[
md === Automatic || md === "ChordLength",
FoldList[Plus, 0, Normalize[(Norm /@ Differences[pts]), Total]] // N,
md === "Centripetal",
FoldList[Plus, 0,
Normalize[(Norm /@ Differences[pts])^(1/2), Total]] // N,
md === "EqualSpaced",
Range[0, 1, 1/n]
];
(*calculate the knots*)
knots =
Join[
ConstantArray[0, sd + 1],
1/sd (Plus @@ (preKnots[[# + 1 ;; # + sd]]) & /@ Range[1, n - sd]),
ConstantArray[1, sd + 1]] // N;
(*calculate the coefficients of matrix*)
coeffMat =
Function[{u0},
With[{i = searchSpan[knots, u0]},
Join[
ConstantArray[0, i - sd],
BSplineBasis[{sd, knots}, #, u0] & /@ Range[i - sd, i],
ConstantArray[0, n - i]]
]
] /@ preKnots;
(*solve the control points of B-Spline curve*)
ctrlpts = LinearSolve[coeffMat, pts];
(*visualize the result*)
Graphics[
{OptionValue[SplineColor],
BSplineCurve[ctrlpts, SplineDegree -> sd, SplineKnots -> knots],
Red, PointSize[Medium], Point[pts]}
]
]


### Test

pts =
{{-1.5, -2}, {-3, -1}, {-2.7, 0.5}, {-1.75, 1.3}, {0.8, -1.5},
{1.5, 0.4}, {0, 2}, {3, 2}};

BSplineInterpolation[pts, Method -> #] & /@
{"ChordLength", "EqualSpaced", "Centripetal"}


Show[
BSplineInterpolation[pts, SplineColor -> #2, Method -> #1] & @@@
{{"ChordLength", Green}, {"EqualSpaced", Black}, {"Centripetal", Pink}}]


BSplineInterpolation[pts, SplineDegree -> #] & /@ {2, 3, 4, 5}


-
The parametrizeCurve[] routine in the procedure belisarius linked to is in fact a generalization; the default setting of the second argument yields centripetal parametrization, while setting it to $1$ gives the chord-length parametrization. –  J. M. is back. May 29 at 9:33
@Guesswhoitis., Yes, got it:) –  Shutao TANG May 29 at 9:40
@AlexeyPopkov, BTW, It is the first step to implement the B-spline interpolation. I'd like to know that is it possible to speed up this implementation by optimizing some code or with the help of Compile. –  Shutao TANG May 29 at 13:26
@ShutaoTang With Compile you are limited to MachinePrecision. The limiting step is LinearSolve which cannot be compiled but (in principle) can be optimized by optimizing the structure of the matrix and applying SparseArray (see an example in my answer). –  Alexey Popkov May 29 at 13:31
@Alexey, as you might have already surmised, Lee intended that parameter of my parametrizeCurve[] as a way to smoothly adjust between the equispaced, centripetal, and chord length parametrizations; thus, classically, it is usually in $[0, 1]$. I will edit my other answer later to explain this. –  J. M. is back. May 29 at 15:38