Data fitting using thin-plate spline/interpolation

I have several points in 3D space. I have used the surface fitting function to create a pretty accurate representation.

Now I need to generate an equation with the corresponding coefficients for this surface, so that I can use it elsewhere without having the original data points.

1) What is the equation for the thin-plate interpolation? Like when you use polynomial, sine, fourier, etc, there is a clear equation that describes the data....

2) How do I extract the coefficients? Coeffx only returns 'thinplateinterp'.

x=[0,0,0,0.5,0.5,0.5,1,1,1];
y=[0.001,0.1,10,0.001,0.1,10,0.001,0.1,10];
z=[4.35,1.37,2.31,7.04,3.58,4.14,2.08,2.35,2.43];

ylog=log10(y);

[xData, yData, zData] = prepareSurfaceData( x, ylog, z );

ft = 'thinplateinterp';

[fitresult, gof] = fit( [xData, yData], zData, ft, 'Normalize', 'on' );
Coeffx=coeffvalues(fitresult);

Answer 1

[William Rose]

@Pelajar UM,

I too am unable to find the coefficients for the thin plate psline model.

I can evaluate the fitted model in two ways:

>> fitresult(.5,-3)

ans =

7.0400

>> evaluate(Coeffx,.5,-3)

ans =

7.0400

The reult with "evaluate(Coeffx,.5,-3)" gives the answer above only if I compute the fit with 'Normalize','off'.

If I compute the fit with 'Normlaize','On', then the output from fitresult(.5,-3) is correct, and is the value shown above, but the output from evaluate(Coeffx,.5,-3) is incorrect.

The wikipedia article for thin plate spline says that there are 2(K+3) parameters in the thin plate spline model, where K=number of fitted points. In your example, K=9, so there should be 24 parameters.

I did the fit in the CurveFitter app also, with the same results. I was also unable to find the coefficients in the curve fitter app, and when I selected the "export to workspace" option, it just exported variables like fitresult, so this did not help identify the parameters.

One can also fit a thin plate spline to the same data by:

p=1; st=tpaps([xData',yData'],zData',p)

The resulting structure, st, includes st.coefs, which is a vector of 12 coefficients. By setting p=1, tpaps() returns an exact-fitting thin plate spline. As p decreases toward 0, the fit becomes increasingly approximate. See the help for tpaps().

Good luck with your fitting.

[Pelajar UM]

Thanks. I tried to use tpaps as you suggested but I get this error:

Error using tpaps

TPAPS(X,...) expects the data sites to form the COLUMNS of X. With that interpretation, you seem to be providing data sites in 1-dimensional space, while, at present, TPAPS can only handle data sites in the plane.

[William Rose]

@Pelajar UM,

For tpaps, do as follows:

x=[0,0,0,0.5,0.5,0.5,1,1,1];
y=[0.001,0.1,10,0.001,0.1,10,0.001,0.1,10];
z=[4.35,1.37,2.31,7.04,3.58,4.14,2.08,2.35,2.43];
ylog=log10(y);
p=1; 
st=tpaps([x;ylog],z,p);
fnval(st,[.5;-3])  %evaluate the thin plate spline function

ans = 7.0400

Note that the point I chose for evaluation is one of the orignal points, and the value returned is exactly the corresponding z value, as we expect it to be.

Good luck.

[Pelajar UM]

Thanks a lot. This works now.

Two questions though:

1. How can I extend the surface beyond the data points. For example, using the curve fitter app, when I zoom out, I can see the extended surface beyond ylog=1, but with tpaps, the surface is only between ylog = -3 and 1. I am using fnplt(st) because plot(st), surf(st), plot3(st), etc don't seem to work.

UPDATE: you can do this by

st.interv = {[0 1] [-3 3]};

And now you see a massive difference between the surface generated by tpaps

vs the one from Curve Fitter using thin-plate spline:

To me, the surface from curve fitter is most certainly a better representation of the data.

2. Now that I have access to the coefficients, how do I use them to replicate this surface elsewhere? Basically, what is the point of having these coefficients if you cannot plug them into something?

UPDATE: answer to question 2 is here: https://www.mathworks.com/help/curvefit/stmak.html

+page 36: http://cda.psych.uiuc.edu/matlab_pdf/splines.pdf

So my only problem now is that the tpaps results are not as good as the one from curve fitter. I have played around with p to no avail...

[William Rose]

@Pelajar UM, Congratulations on tracking down the Matlab document that contains details of the formulas for the thin plate spline. I don;t know why the surfce generated with tpaps() appears to differ from the surface generated by the curve fitter app. There are multiple options in the curve fitter app. Perhaps the settings on some of the options in CurveFitter differ from the default settings and assumptions for tpaps().

[Pelajar UM]

Thanks. So when Normalize is off, then the curve fitter app gives exactly the same result as tpaps.

Thin-plate spline interpolant:
f(x,y) = thin-plate spline computed from p
where x is normalized by mean 0.5 and std 0.4564
and where y is normalized by mean -1 and std 1.528

I wonder if there is a way to implement the normalize function for tpaps too...

Also instead of evaluate(Coeffx,0.5,-3), use fitresult(0.5,-3). It gives the correct answer with normalize = on.

[William Rose]

@Pelajar UM,

Good job figuring out the options in CurveFitter that make it consistent with tpaps().

As for fitresult() versus evaluate(), when normalize is on or off, I did attempt to point that out in my comment yesterday:

Three ways to evaluate the spline:

fitresult(.5,-3)
evaluate(Coeffx,.5,-3)  % works if the fit is done with normalize off
fnval(st,[.5,-3])       % st is obtained with tpaps()

[Pelajar UM]

Right. Sorry. I missed that.

Also UPDATE: this would give exactly the same surface as the one from curve fitter.

x=normalize(x, 'center', 'mean', 'scale', 'std');
ylog=normalize(ylog, 'center', 'mean', 'scale', 'std');

But for some weird reason, it also shifts the x and y axes (unlike curve fitter) which is something I need to avoid. For instance, it centers the x axis at 0, even though the mean is 0.5....

[William Rose]

@Pelajar UM, that is interesting and strange.

I didn't know about the normalize command. Thank you for bringing it to my attention. It will be useful. This is one of the reasons why it is nice to visit this site.

[Pelajar UM]

Last update:

Long story short, it is better that you use fitrsult to generate new data points in the desired range and then just feed these into a high-degree polynomial. It's not as accurate (also depends on how many points you are willing to generate) but in my case, I got R^2 = 0.935, by almost doubling the number of my data points. However one must note that the resulting polynomial gives absolutely nonsensical predictions outside the data range .

There is a way to implement the spline (and I think it involves denoramlizing the data after fitting), but it is just too much headache for a small increase in accuracy, especially if your end goal is to replicate this outside MATLAB.

[William Rose]

@Pelajar UM, thanks for the update.

Answer 2

[John D'Errico]

This is a common problem people have with splines in any form. They hope to see some nice simple function they can write down. Maybe put the function into a paper, or use in some way they would use other functional forms.

Sorry. That is pretty much never the case. Splines are great things for some purposes. But if you want to write down that nice looking function, it won't happen.

You can evaluate the function at any point. There are methods provided for that purpose. But if you really need the pretty function you can write down, you will need to use other tools.

x=[0,0,0,0.5,0.5,0.5,1,1,1];
y=[0.001,0.1,10,0.001,0.1,10,0.001,0.1,10];
z=[4.35,1.37,2.31,7.04,3.58,4.14,2.08,2.35,2.43];
ylog=log10(y);

surf(reshape(x,3,3),reshape(ylog,3,3),reshape(z,3,3))

Anyway, using a thin plate spline for such minimal data is probably wild overkill. At best, a simple polynomial model may be sufficient, and about as good as you can do with such limited data.

mdl = fit([x',ylog'],z','poly22')

     Linear model Poly22:
     mdl(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2
     Coefficients (with 95% confidence bounds):
       p00 =       1.457  (-0.9711, 3.886)
       p10 =       9.961  (-0.05327, 19.97)
       p01 =    -0.03542  (-1.641, 1.57)
       p20 =      -9.753  (-19.24, -0.2679)
       p11 =      0.5975  (-1.079, 2.274)
       p02 =      0.3229  (-0.2699, 0.9158)

This gives you a nice model you can write down. Not a very good model, but your data is highly limited.

[mdl(x,ylog);z]

ans = 2×9

4.4697 1.8156 1.7447 6.1156 4.0589 4.5856 2.8847 1.4256 2.5497 4.3500 1.3700 2.3100 7.0400 3.5800 4.1400 2.0800 2.3500 2.4300

As you can see, the model does predict the data reasonably well.

[xint,yint] = meshgrid(linspace(0,1),linspace(-3,1));
surf(xint,yint,mdl(xint,yint))
hold on
plot3(x,ylog,z,'ro')

You can even plot it as a nice smooth, well behaved surface.

[Pelajar UM]

Thanks for your help. I understand that my data is limited.

But can I generate more data from the spline fit and then use that data to come up with a more accurate polynomial fit?

Also the original question remains. Where are these coeffients stored and how would one use them? I am not necessarily looking for a simple equation, but the idea is to be able to generate this surface and extract data from it without having the original data points.

[William Rose]

@Pelajar UM, your question is a reasonable one, and I don't know the answer to where the coefficients can be found. I demonstrated with tpaps() how to find 12 coefficients. This is only half the number of coefficients expected, according to the wikipedia article on thin plate splines. I also demonstrated three ways you can use a spline fit to predict values at other points:

fitresult(.5,-3)
evaluate(Coeffx,.5,-3)  %works if the fit is done with normalize off
fnval(st,[.5,-3])       %st is obtained with tpaps()