How to get the values of biharmonic interpolation - 3D surface fitting

Hi there,

I used the 'biharmonicinterp' method for 3D surface fitting and got the following surface based on 4 given data points (blue dots), which was more reasonable than other fitting method.

However, it is hard to get the interpolation values between these points. For example, I would like to slice the surface and figure put the z value change along a certain axis, but could not get the previous interpolated values. I try to use other cfit methods, but the cross section is not consistent with the 3D surface.

Could you help me to get the interpolated values by using 'biharmonicinterp' method? Or directly slice the fitted 3D surface to obtain the z value change? Thank you very much.

The code is attached below:

% Data
xData = [0;85;0;-85];
yData = [85;0;-85;0];
zData = [0.5695;0.1766;-0.2984;-0.0129];

% Set up fittype and options.
ft = 'biharmonicinterp';

% Fit model to data.
[fitresult, gof, output] = fit( [xData, yData], zData, ft, 'Normalize', 'on' );

% Plot fit with data.
figure(1)
h = plot( fitresult, [xData, yData], zData );
legend('off')
grid on
colormap(jet)
shading interp
colorbar

% Label axes
xlabel( 'X', 'Interpreter', 'none' );
ylabel( 'Y', 'Interpreter', 'none' );
zlabel( 'Z', 'Interpreter', 'none' );

% Tick label
xticks(-80:20:80)
yticks(-80:20:80)

% Paper position
set(gcf,'PaperPositionMode','auto');

[Hao Shi]

Thank you John for your generous help and clear explanation, your method works very well! Actually I have another question for this interpolation problem, I have tried several days but cannot get through it.

You may notice the points in the 4 corners seem unreasonable, I would like to remove these points but do not know how to make it. Could we interpolate the data in a polar cooedinates (Center: [0, 0], Radius: 85) or directly remove the points outside the circle (X^2 + y^2 > 85^2)?

Thank you very much and I will appreciate your any help.

[John D'Errico]

Yes, there are several things you could reasonably do. But before I suggest some of them, do you seriously only have 4 data points to want to define the surface?

Yes, the surface as plotted may seem strange in the corners, but do you really care? Should you be even trying to interpolate in those corners, extrapolating that far out given so little data?

My point is that you should not care what the surface predicts in the corners if you don't need to evaluate it there. So the answer is, if you want only to produce a pretty plot, then just don't plot the corners. If you are using the data to extrapolate there, WHY? That is just a bad idea, using any model.

[Hao Shi]

Dear John,

Thank you for your answer. The thing is that I only have four data points at each time, because of the limited measurement locations in my experiment, and we could show sth between these points by using the interpolation as it meets the expectations.

However, my facility is a cylindrical domain, so practically we do not need the points in the corners. In the Cartesian coordinate system, we defined these points at the beginning, can you teach me how to get rid of this or how to remove these points finally? Thank you very much.

I was thinking, maybe we could define another cylindrical domain [x,y], just as you defined the cross-section coordinates, then calculate the z based on the fit results, finally plot them out. Could you give me some suggestions?

[John D'Errico]

I'd go with the simplest solution.

You have two angles, you called them x and y. But just rotate the universe. x and y both vary from -85 to 85. I'll transform the problem as

u = (x + y)/sqrt(2)

v = (x - y)/sqrt(2)

We also can go the other way of course, using a symmetric transformation. (That sqrt(2) comes from a coordinate rotation of course, since we have sin(pi/4)=cos(pi/4)=1/sqrt(2).)

x = (u + v)/sqrt(2)

y = (u - v)/sqrt(2)

Why would we do this? Because it makes the problem now trivially simple to interpolate, where the interpolation lives only on the domain of interest.

xData = [0;85;0;-85];
yData = [85;0;-85;0];
zData = [0.5695;0.1766;-0.2984;-0.0129];

uData = (xData + yData)/sqrt(2);
vData = (xData - yData)/sqrt(2);

format short g
[uData,vData,zData]
ans =
       60.104      -60.104       0.5695
       60.104       60.104       0.1766
      -60.104       60.104      -0.2984
      -60.104      -60.104      -0.0129

plot(uData,vData,'ro')
xlabel u
ylabel v

zdata does not change of course, since the only thing that has changed is how we perceive the data. This is just an orthogonal rotation. So now you can do your fits in terms of z, as a function of u and v. The world does not change just because we look at it from a 45 degree angle.

To predict a value based on x and y, first rotate them into the (u,v) orientation. Then interpolate. as you will. Now no extrapolation need be done. As long as your point falls inside the diamond shaped region in (x,y), it will fall inside the square in (u,v). For example:

ft = 'biharmonicinterp';
% Fit model to data.
[fitresult, gof, output] = fit( [uData, vData], zData, ft, 'Normalize', 'on' );
plot(fitresult)
xlabel u
ylabel v
zlabel z

Is it overkill to fit a bi-cubic surface through 4 data points? Yes, I still think that to be true. Honestly, I think the model below is as much as you can justify:

[fitresult2, gof, output] = fit( [uData, vData], zData, 'a + b*x + c*y + d*x*y', 'Normalize', 'on' );
plot(fitresult2)
xlabel u
ylabel v
zlabel z

The latter case also goes EXACTLY through the 4 data points. It seems to be almost a planar fit, though it is not perfectly so. If you look carefully, rotating the plot around, there is some small curvature. Not a lot, but your 4 data points almost fall in a plane themselves.

Don't overfit your data. There is no reason to suppose the shape of that surface is what the biharmonic fit produces. And there is insufficient information to ask for any more than the alternate fit I chose.

[Hao Shi]

Dear John,

Thank you so much for your kind help, I really appreciate it. The orthogonal rotation is a good option to get rid of the corner points, it works good for me. Meanwhile, I also try to define the points in the circular domain, then fit the zdata again with the generated fitresult model, finally we could see a pefect circle at each time. I have leanered a lot from your suggesitions, thanks a lot!

[John D'Errico]

I had seriously thought of suggesting a polar coordinate transformation, which might have also worked reasonably, at least at first glance. I did not because it seemed to be a confusing solution. That is, convert a problem that is alreadty defined in terms of angles, into a problem where we introduce another angle due to the polar coordinate system. It just seems like something I did not want to do, not when a simpler solution was available.

However, since you brought it up, yes, I might as well discuss what happens if we do. That is, convert the region of interest into a polar coordinate problem. Now you would tranform the (x,y) plane into an (r,theta) domain. But here we see a problem suddenly arises. We have 4 data points, all of which lie at one radius.

[theta,r] = cart2pol(xData,yData);

[theta,r]
ans =
           1.5707963267949                        85
                         0                        85
          -1.5707963267949                        85
          3.14159265358979                        85

So, we have 4 data points now, ALL of which lie st the same radius from the origin. They are all on a circle. So we have a function, where you might decide to interpolate circularly, around the perimeter of that circle. But now inferring what happens for smaller radii inside that circle is far more difficult.

Yes, you could make up some value at r==0. But if you knew that point, you should have provided it in the first place. You cannot use interpolation to decide what value should happen at the center of the circle, and then fit a model to that. This is far too mathematically incestuous, far too dangerous a thing to do.

And since the simple coordinate rotation provides such a good solution, just use that. And of course, I know that I hound people about this, but getting more data, IF you could actually get the data from some experimentation would always be a very good thing, greatly improving the predictive behavior of any model.

[Hao Shi]

I see, I will try to get more data from additional positions, thank you very much John!

Answer 1

[John D'Errico]

I don't see the problem. However, I would note that you have only 4 points, so the surface is pretty difficult to create in an intelligent way. And the use of any interpolation across the center of a domain like this can produce some interesting interpolation artifacts, if you know exactly where to look. (Having worked with interpolation in multiple dimensions for many years, I can attest to that as fact.) Anyway, what you did seems to produce something that I fully expect.

For example, how would I get a slice through that region? For eample, fix x at zero, then vary y?

[xData,yData,zData]
ans =
            0           85       0.5695
           85            0       0.1766
            0          -85      -0.2984
          -85            0      -0.0129

So my expectation would have z varying from -0.2984 to 0.5695. The shape of that curve along this path will probably be something smooth, but it will have some shape. Not necesarily just a straight line.

xint = zeros(1,50);
yint = linspace(-85,85,50);
zint = fitresult(xint,yint);
plot(yint,zint,'b')

So not at all unreasonable. Consistent with the limited data at the endpoints. Now, without copying all of the code you wrote, I'll add these last two lines to your code.

hold on
H = plot3(xint,yint,zint,'ko');

I fail to see the problem. I see a line with a lot of black o's that follows the shape of the surface you created, from corner to corner of that region. Now, where do you see a problem in what I did? I don't. The curve looks eminently reasonable, and is consistent with the plotted surface. Did you expect a straight line connecting those two end points?

Answer 2

[darova]

• re-create surface manually

• Use interp2

[x1,y1] = meshgrid(-50:5:50);
z1 = fitresult(x1,y1);
xx = linspace(-50,50);
yy = 2*xx;
zz = interp2(x1,y1,z1,xx,yy);
surf(x1,y1,z1)
line(xx,yy,zz,'linew',3)