Reconstruct multivariate spline from csapi

Using the most Excellent Curve Fitting toolbox, I encounter difficulties in reconstructing a bivariate spline that I created from fitting 2D grid data using sp = csapi({x,y}, data).

My problem here is that I want to use the coefficient matrix I get (by calling squeeze(sp.coefs)) in a C header file, and then compute the spline and all of its first and second derivatives inside my C code.

Unfortunately, the elsewise well-documented Curve Fitting Toolbox exhibits some "don't break your head too much over internal mechanisms" manner in this context, rather than giving the math for a multivariate pp-form spline in the context of the coefficient matrix.

Can anybody help me out here? Instead of being a 3D tensor with "depth" 4 (for the coefficients of the bi-cubic spline per quadrant), the coefficient matrix is a 2D matrix, and there is no documentation on which entry of this matrix goes where in the bi-cubic spline formula that I can only expect to be $\sum_{i=0}^3 \sum_{j=0}^3 x^i y^j a_{ij}$ with $a_{ij}$ being the "proper" place in the coefficient matrix.

Answer 1

[Wolfgang]

Okay, now I finally figured out what's going on with the ppform for multivariate tensor product splines. No thanks at all to the crappy documentation, though. Since some other posts are related to this question, I'll give a code example here that might be helpful if included in the next documentation as well.

Using the same bivariate example as in the toolbox' documentation, we write the following as initialization:

x = -2:.5:2; y = -1:.25:1; [xx, yy] = ndgrid(x,y);
z = exp(-(xx.^2+yy.^2));
sp = csapi({x,y}, z);

For ease of access, we do a little reshaping of the coefficients

coefs = reshape(sp.coefs, [1, 8, 4, 8, 4]);

Let's say we want to evaluate the spline at

x = [0.75; -0.2]

with

fnval(sp,[0.75 -0.2]')

which gives us the result

0.5491

Now for the manual evaluation. We can determine the proper indices of the polynomial piece to use as

i = [ find(x(1) >= sp.breaks{1}, 1, 'last')
      find(x(2) >= sp.breaks{2}, 1, 'last') ];

The value of the bicubic spline is then

y = 0;
for r1 = 1 : 4
  for r2 = 1 : 4
    r = [r1; r2];
    y = y + coefs(1, i(1), r(1), i(2), r(2)) * ...
            (x(1) - sp.breaks{1}(i(1)))^(sp.order(1) - r(1)) * ...
            (x(2) - sp.breaks{2}(i(2)))^(sp.order(2) - r(2));
  end;
end;

and stored in the variable y as

0.5491

Et voilá! Here it is. Not even really tough to figure out, if seen from my new perspective, but the way the documentation is written it was like re-inventing the whole thing. As I said, no thanks at all to the toolbox documentation!

[dominik]

Thanks! 2h of manual werent helping, but your post did.

Answer 2

[Wolfgang]

Okay, so after wasting 22 sleepless hours with the curve fitting toolbox' manual, the following are clear:

• the manual is full of typos
• explanations of both the pp-form and the B-splines are incomplete
• the solution I seek is hidden somewhere in the manual pages of ppmak, but the crucial paragraph ("Then, coefs is interpreted as an ..." up to "This is, in fact, the internal interpretation"), even after figuring out the nowhere documented notation, is non-comprehensible (e.g. $i(\mu)$ for a $l$ polynomial pieces, with $l > m$, appears to only use the first $m$ entries of dimension $i(\mu)$ in the "equivalent reading" of coefs)
• Prof. de Boor seems to have written this manual only to boost sales of his book

Well, I won't buy his book, not after this desastrous escapade. The handbook is utter crap, completely unusuable for anything but the most basic needs, and de Boor's smug comments about his readers' intelligence ("Fortunately, the toolbox is set up in such a way that there is usually no reason for you to concern yourself with these details of either form.") don't help at all.

Honestly, the price my university pays for the campus licenses of Matlab is way too high. Even though the tool per se is "okay", the documentation sucks from here to Natick, Massachusetts.