2 replaced http://mathematica.stackexchange.com/ with https://mathematica.stackexchange.com/

If you want further confirmation, you can try recovering the control points and knots from ipf[], using a procedure similar to the one in my previous answermy previous answer, and compare them with the control points and knots I generated here.

If you want further confirmation, you can try recovering the control points and knots from ipf[], using a procedure similar to the one in my previous answer, and compare them with the control points and knots I generated here.

If you want further confirmation, you can try recovering the control points and knots from ipf[], using a procedure similar to the one in my previous answer, and compare them with the control points and knots I generated here.

1

I still haven't figured out how to write a routine for arbitrary dimension, but I'm posting my (incomplete!) solution in case people might have ideas on extending what I have.

## Bivariate interpolant

Here is a random bivariate polynomial, which we'll use for generating test data:

f[x_, y_] := -2 + 4 x^2 + 4 x^3 - 3 y - 5 x^2 y + 5 y^2 + 5 x y^2 + y^3


Here's some test data from f[x, y], sampled at non-equispaced points:

da = Flatten[Table[{{x, y}, f[x, y]},
{x, {-2, -4/3, -1/5, 1/9, 3/4, 1}},
{y, {0, 1/6, 3/8, 9/5, 2}}], 1];


Here's the reference interpolating function:

{p, q} = {2, 3}; (* spline degrees in the two variables *)
ipf = Interpolation[da, InterpolationOrder -> {p, q}, Method -> "Spline"];


Some preliminary processing to separate out independent and dependent variables:

{pts, vals} = Transpose[SplitBy[SortBy[da, First], #[[1, 1]] &], {3, 2, 1}];


Make the knot sequence for each independent variable:

makeKnots[list_?VectorQ, deg_Integer?NonNegative] :=
With[{n = Length[list]},
Join[ConstantArray[list[[1]], deg + 1],
If[deg + 2 <= n, MovingAverage[ArrayPad[list, -1], deg], {}],
ConstantArray[list[[-1]], deg + 1]]]

{u, v} = {pts[[1, All, 1]], pts[[All, 1, 2]]};
{uk, vk} = MapThread[makeKnots, {{u, v}, {p, q}}];


Build the control points:

{m, n} = {Length[u], Length[v]};
usol = LinearSolve[Outer[BSplineBasis[{p, uk}, #2, #1] &,
u, Range[0, m - 1], 1]];
vsol = LinearSolve[Outer[BSplineBasis[{q, vk}, #2, #1] &,
v, Range[0, n - 1], 1]];

cpts = vsol /@ Transpose[usol /@ vals];


Finally, the bivariate interpolating spline:

spf[x_, y_] = Fold[Dot, cpts, {Table[BSplineBasis[{q, vk}, k - 1, y], {k, n}],
Table[BSplineBasis[{p, uk}, k - 1, x], {k, m}]}];


Plot the two interpolants and the original function:

MapThread[Plot3D[#1[x, y], {x, -2, 1}, {y, 0, 2}, PlotLabel -> #2] &,
{{f, ipf, spf}, {"True", "InterpolatingFunction", "B-spline"}}]
// GraphicsRow


Evaluate ipf[] and spf[] at the same argument:

{ipf[-1, 1], spf[-1, 1]}
{-8.59478, -552429212/64275003}


Note that only the second function gave exact output.

The difference between ipf[] and spf[], showing good agreement:

Plot3D[ipf[x, y] - spf[x, y], {x, -2, 1}, {y, 0, 2}, PlotRange -> All]


If you want further confirmation, you can try recovering the control points and knots from ipf[], using a procedure similar to the one in my previous answer, and compare them with the control points and knots I generated here.

## Trivariate interpolant

Hopefully, you can see the similarities and differences between the previous example and this one:

(* random polynomial *)
mp[x_, y_, z_] := 3 - 7 x^2 + 2 x^3 - 2 y + 8 y^2 + 5 x y^2 + 8 y^3 + 4 x y^3 -
2 y^4 + 6 x z + 2 x^2 z - 6 x^3 z + 2 y z + 4 x y z -
2 x^2 y z - 3 y^2 z - 8 x y^2 z + 7 y^3 z - 5 z^2 + x z^2 +
2 x^2 z^2 + 6 y z^2 - 4 y^2 z^2 + 9 z^3 - 4 x z^3 -
3 y z^3 - 9 z^4
(* random data *)
da = Flatten[Table[{{x, y, z}, mp[x, y, z]}, {x, {-2, -3/2, -9/7, 7/9, 2, 3}},
{y, {1, 5/3, 5/2, 9/2, 5}}, {z, {-1, 5/7, 11/10, 2}}], 2];

{p, q, r} = {4, 3, 2}; (* B-spline degree for each variable *)
ipf = Interpolation[da, InterpolationOrder -> {p, q, r}, Method -> "Spline"];

{pts, vals} = Transpose[GatherBy[da, {#[[1, 1]] &, #[[1, 2]] &}], {4, 3, 2, 1}];

(* make knots *)
{u, v, w} = {pts[[1, 1, All, 1]], pts[[1, All, 1, 2]], pts[[All, 1, 1, 3]]};
{uk, vk, wk} = MapThread[makeKnots, {{u, v, w}, {p, q, r}}];

(* make control points *)
{l, m, n} = Length /@ {u, v, w};
usol = LinearSolve[Outer[BSplineBasis[{p, uk}, #2, #1] &, u, Range[0, l - 1], 1]];
vsol = LinearSolve[Outer[BSplineBasis[{q, vk}, #2, #1] &, v, Range[0, m - 1], 1]];
wsol = LinearSolve[Outer[BSplineBasis[{r, wk}, #2, #1] &, w, Range[0, n - 1], 1]];
cpts = Map[wsol, Transpose[Map[vsol,
Transpose[Map[usol, vals, {2}], {2, 3, 1}], {2}], {1, 3, 2}], {2}];

(* B-spline interpolant *)
spf[x_, y_, z_] = Fold[Dot, cpts, {Table[BSplineBasis[{r, wk}, k - 1, z], {k, n}],
Table[BSplineBasis[{q, vk}, k - 1, y], {k, m}],
Table[BSplineBasis[{p, uk}, k - 1, x], {k, l}]}];


Tests:

{ipf[2, 3, 1], spf[2, 3, 1]}
{383.531, 10447533501473/27240371550}

ContourPlot3D[#[x, y, z], {x, -2, 3}, {y, 1, 5}, {z, -1, 2},
BoxRatios -> Automatic, MaxRecursion -> 0] & /@
{ipf, spf} // GraphicsRow