2D interpolation problem with scattered data

Question

I have a function F which allows me to calculate 'y' based on x and z values. i.e. y = F(x,z). This function is very computationally intensive; it takes around 10 mins per single evaluation.

What I really need is to find a value of "z" given a pair of (x,y) values, i.e. I need z = G(x,y). Solving F is not possible.

To numerically implement G, I've gone through and calculated many y values using F(x,z) using grid in x and z.

Attached is a plot of the z values in a contour plot. The red dots represents the (x,y) points for which I have the z values for. It's hard to see so zoomed-out, by the spacing in y of these points is irregular. On the right is also a surf plot showing that G(x,y) is smoothed and well-behaved.

figure(1); clf; subplot(1,2,1); hold all; box on; grid on;
scatter3(X_mesh(:),Y_mesh(:),Z_mesh(:),5,'filled','ro')
contourf(X_mesh,Y_mesh,Z_mesh,20)
xlabel('x'); ylabel('y'); zlabel('z'); hcolorbar = colorbar; ylabel(hcolorbar,'z')
subplot(1,2,2); hold all; box on; grid on;
surf(X_mesh,Y_mesh,Z_mesh);
xlabel('x'); ylabel('y'); zlabel('z');

However, when I try to interpolate this surface, the smoothness of the surface gets messed up and I get lots of noise. It happens if I use griddata or scatteredInterpolant. I've tried different methods: 'linear','natural','cubic', etc. and it's all the same. (Apologies for some streaks on the figures. I'm having graphics cards issues right now.)

figure(2); clf; subplot(1,2,1); hold all; box on; grid on;
[Xout_mesh, Yout_mesh] = meshgrid(linspace(5,60,100),linspace(0.4,300,100));
Zout_mesh = griddata(X_mesh(:),Y_mesh(:),Z_mesh(:),Xout_mesh,Yout_mesh,'linear');
contourf(X_mesh,Y_mesh,Z_mesh);
xlabel('x'); ylabel('y'); zlabel('z'); hcolorbar = colorbar; ylabel(hcolorbar,'z')
subplot(1,2,2); hold all; box on; grid on;
surf(X_mesh,Y_mesh,Z_mesh);
xlabel('x'); ylabel('y'); zlabel('z');

Does anyone have a clue what's going on? I can't use any regular spaced data interpolation because my y points are not quite regularly spaced. Below I should a plot where I checked the triangular meshing and it looks fine.

figure(3); clf; 
subplot(1,2,1); hold all; box on; grid on;
dt = delaunayTriangulation(X_mesh(:),Y_mesh(:));
xlabel('x'); ylabel('y');
triplot(dt)
subplot(1,2,2); hold all; box on; grid on;
triplot(dt)
xlabel('x'); ylabel('y');

For completeness, X_mesh, Y_mesh and Z_mesh outputted to console are below. All three matrices are 19x51 doubles.

X_mesh:

  Columns 1 through 20
 1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1
 2     2     2     2     2     2     2     2     2     2     2     2     2     2     2     2     2     2     2     2
 3     3     3     3     3     3     3     3     3     3     3     3     3     3     3     3     3     3     3     3
 4     4     4     4     4     4     4     4     4     4     4     4     4     4     4     4     4     4     4     4
 5     5     5     5     5     5     5     5     5     5     5     5     5     5     5     5     5     5     5     5
 6     6     6     6     6     6     6     6     6     6     6     6     6     6     6     6     6     6     6     6
 7     7     7     7     7     7     7     7     7     7     7     7     7     7     7     7     7     7     7     7
 8     8     8     8     8     8     8     8     8     8     8     8     8     8     8     8     8     8     8     8
10    10    10    10    10    10    10    10    10    10    10    10    10    10    10    10    10    10    10    10
12    12    12    12    12    12    12    12    12    12    12    12    12    12    12    12    12    12    12    12
14    14    14    14    14    14    14    14    14    14    14    14    14    14    14    14    14    14    14    14
16    16    16    16    16    16    16    16    16    16    16    16    16    16    16    16    16    16    16    16
18    18    18    18    18    18    18    18    18    18    18    18    18    18    18    18    18    18    18    18
20    20    20    20    20    20    20    20    20    20    20    20    20    20    20    20    20    20    20    20
25    25    25    25    25    25    25    25    25    25    25    25    25    25    25    25    25    25    25    25
30    30    30    30    30    30    30    30    30    30    30    30    30    30    30    30    30    30    30    30
40    40    40    40    40    40    40    40    40    40    40    40    40    40    40    40    40    40    40    40
50    50    50    50    50    50    50    50    50    50    50    50    50    50    50    50    50    50    50    50
60    60    60    60    60    60    60    60    60    60    60    60    60    60    60    60    60    60    60    60

Y_mesh:

  Columns 1 through 12
0.0242    0.0539    0.0764    0.0977    0.1112    0.1184    0.1203    0.1179    0.1118    0.1028    0.0914    0.0782
0.0521    0.1240    0.1890    0.2671    0.3370    0.3995    0.4552    0.5046    0.5478    0.5854    0.6175    0.6445
0.0801    0.1949    0.3036    0.4412    0.5717    0.6958    0.8138    0.9257    1.0318    1.1321    1.2265    1.3151
0.1081    0.2659    0.4186    0.6163    0.8087    0.9961    1.1788    1.3567    1.5298    1.6981    1.8613    2.0193
0.1361    0.3369    0.5337    0.7919    1.0465    1.2980    1.5465    1.7919    2.0341    2.2728    2.5079    2.7391
0.1641    0.4080    0.6489    0.9677    1.2849    1.6009    1.9158    2.2295    2.5419    2.8525    3.1613    3.4678
0.1921    0.4791    0.7642    1.1436    1.5235    1.9044    2.2861    2.6687    3.0519    3.4354    3.8189    4.2020
0.2202    0.5502    0.8795    1.3197    1.7625    2.2083    2.6571    3.1089    3.5633    4.0201    4.4790    4.9396
0.2762    0.6925    1.1102    1.6719    2.2406    2.8167    3.4001    3.9908    4.5886    5.1931    5.8039    6.4208
0.3322    0.8347    1.3409    2.0243    2.7191    3.4256    4.1439    4.8741    5.6157    6.3686    7.1322    7.9063
0.3883    0.9770    1.5716    2.3768    3.1976    4.0347    4.8882    5.7579    6.6438    7.5454    8.4624    9.3944
0.4443    1.1193    1.8024    2.7292    3.6762    4.6440    5.6327    6.6422    7.6724    8.7231    9.7937   10.8840
0.5004    1.2616    2.0332    3.0817    4.1549    5.2533    6.3773    7.5267    8.7015    9.9013   11.1258   12.3745
0.5564    1.4038    2.2639    3.4343    4.6336    5.8628    7.1221    8.4114    9.7308   11.0799   12.4583   13.8657
0.6966    1.7595    2.8409    4.3156    5.8305    7.3866    8.9843   10.6237   12.3048   14.0273   15.7910   17.5954
0.8367    2.1152    3.4178    5.1970    7.0275    8.9106   10.8468   12.8364   14.8793   16.9756   19.1248   21.3267
1.1169    2.8267    4.5718    6.9598    9.4216   11.9588   14.5722   17.2623   20.0294   22.8735   25.7942   28.7914
1.3972    3.5381    5.7257    8.7227   11.8158   15.0072   18.2979   21.6888   25.1801   28.7722   32.4648   36.2576
1.6774    4.2495    6.8797   10.4856   14.2101   18.0556   22.0238   26.1154   30.3312   34.6713   39.1359   43.7246

Z_mesh:

  Columns 1 through 12
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000
1.0000    2.5000    4.0000    6.0000    8.0000   10.0000   12.0000   14.0000   16.0000   18.0000   20.0000   22.0000

Update: I've found a clue to the problem. Below is a plot of the original (uninterpolated) data with shading interp turned on using "surf" and "trisurf" plotting. Surf produces a pretty smooth surface, whereas with trisurf streaks start appearing. These, I believe, are the same streaks as seen with griddata or scatteredInterpolant, which uses a triangular mesh.

Does anyone know how I can access the interpolated values from the surface object?

Can you try rho_out_mesh = interpn(X_mesh(:),Y_mesh(:),Z_mesh(:),Xout_mesh,Yout_mesh) and share the results? griddata tries to interpolate on a surface and given non-uniformity of y points, it may be estimating an uneven surface. — Abhineet Gupta
Thanks for the suggestion. I can't use the interpn function because the Y_mesh spacing isn't regular. (It's close but not quite). I get the following error: Error using griddedInterpolant The grid vectors are not strictly monotonic increasing. Error in interpn (line 147) F = griddedInterpolant(X, V, method,extrap); I tried sorting X_mesh, Y_mesh and Z_mesh first and get the same error. — Kent
Ahh sorry. All 3 matrices are 19x51, so Matlab couldn't show all 51 columns on the screen. For X_mesh, since Matlab outputted them as integers, it showed 20 columns, whereas Y_mesh and Z_mesh 12 columns are displayed. — Kent
I think the problem is that you have many different values of z near x=0 and y=0. It looks to me that you should get better result if you look at your data as a function y=f(x,z), rather than z=f(x,y). The alternative is to detrend (or rather rotate in 3D space), interpolate, then add the trend again. — Cris Luengo

Dev-iL Dev-iL · Accepted Answer · 2019-07-04T08:16:36

You appear not to be using the right tool for the job. A gridded or linear interpolant is not what this problem requires. The right approach here is mentioned in Cris' comment - that is, to rotate the dataset, or simply reorder the dimensions.

Playing around with the data you posted (first 12 columns) in cftool I can tell you that y=f(x,z) can be expressed almost perfectly using a polynomial surface of degree 1-3 both parameters. Here are some results for example:

poly32 model (9 coefficients)

Linear model Poly32:
     f(x,y) = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2 + p30*x^3 + p21*x^2*y 
                    + p12*x*y^2
Coefficients (with 95% confidence bounds):
       p00 =     0.00897  (-0.0102, 0.02813)
       p10 =   -0.002627  (-0.004814, -0.0004397)
       p01 =    -0.00157  (-0.004805, 0.001665)
       p20 =   0.0002045  (0.0001266, 0.0002824)
       p11 =     0.02715  (0.02698, 0.02732)
       p02 =   -0.001751  (-0.001885, -0.001617)
       p30 =  -2.844e-06  (-3.702e-06, -1.985e-06)
       p21 =   8.535e-06  (6.588e-06, 1.048e-05)
       p12 =   0.0002796  (0.000274, 0.0002852)

Goodness of fit:
  SSE: 0.1613
  R-square: 1
  Adjusted R-square: 1
  RMSE: 0.02714

poly12 model (5 coefficients)

Linear model Poly12:
     f(x,y) = p00 + p10*x + p01*y + p11*x*y + p02*y^2
Coefficients (with 95% confidence bounds):
       p00 =      0.4112  (0.3266, 0.4958)
       p10 =    -0.02307  (-0.02588, -0.02026)
       p01 =     -0.1155  (-0.1304, -0.1006)
       p11 =     0.03397  (0.03375, 0.03418)
       p02 =    0.003119  (0.002503, 0.003735)

Goodness of fit:
  SSE: 7.398
  R-square: 0.9995
  Adjusted R-square: 0.9995
  RMSE: 0.1821

Or programmatically,

x = [1:1:8, 10:2:20 25 30:10:60].';
z = [1:1.5:4, 6:2:22];
[X,Z] = ndgrid(x,z);

Y = [
0.0242    0.0539    0.0764    0.0977    0.1112    0.1184    0.1203    0.1179    0.1118    0.1028    0.0914    0.0782
0.0521    0.1240    0.1890    0.2671    0.3370    0.3995    0.4552    0.5046    0.5478    0.5854    0.6175    0.6445
0.0801    0.1949    0.3036    0.4412    0.5717    0.6958    0.8138    0.9257    1.0318    1.1321    1.2265    1.3151
0.1081    0.2659    0.4186    0.6163    0.8087    0.9961    1.1788    1.3567    1.5298    1.6981    1.8613    2.0193
0.1361    0.3369    0.5337    0.7919    1.0465    1.2980    1.5465    1.7919    2.0341    2.2728    2.5079    2.7391
0.1641    0.4080    0.6489    0.9677    1.2849    1.6009    1.9158    2.2295    2.5419    2.8525    3.1613    3.4678
0.1921    0.4791    0.7642    1.1436    1.5235    1.9044    2.2861    2.6687    3.0519    3.4354    3.8189    4.2020
0.2202    0.5502    0.8795    1.3197    1.7625    2.2083    2.6571    3.1089    3.5633    4.0201    4.4790    4.9396
0.2762    0.6925    1.1102    1.6719    2.2406    2.8167    3.4001    3.9908    4.5886    5.1931    5.8039    6.4208
0.3322    0.8347    1.3409    2.0243    2.7191    3.4256    4.1439    4.8741    5.6157    6.3686    7.1322    7.9063
0.3883    0.9770    1.5716    2.3768    3.1976    4.0347    4.8882    5.7579    6.6438    7.5454    8.4624    9.3944
0.4443    1.1193    1.8024    2.7292    3.6762    4.6440    5.6327    6.6422    7.6724    8.7231    9.7937   10.8840
0.5004    1.2616    2.0332    3.0817    4.1549    5.2533    6.3773    7.5267    8.7015    9.9013   11.1258   12.3745
0.5564    1.4038    2.2639    3.4343    4.6336    5.8628    7.1221    8.4114    9.7308   11.0799   12.4583   13.8657
0.6966    1.7595    2.8409    4.3156    5.8305    7.3866    8.9843   10.6237   12.3048   14.0273   15.7910   17.5954
0.8367    2.1152    3.4178    5.1970    7.0275    8.9106   10.8468   12.8364   14.8793   16.9756   19.1248   21.3267
1.1169    2.8267    4.5718    6.9598    9.4216   11.9588   14.5722   17.2623   20.0294   22.8735   25.7942   28.7914
1.3972    3.5381    5.7257    8.7227   11.8158   15.0072   18.2979   21.6888   25.1801   28.7722   32.4648   36.2576
1.6774    4.2495    6.8797   10.4856   14.2101   18.0556   22.0238   26.1154   30.3312   34.6713   39.1359   43.7246];


[xData, yData, zData] = prepareSurfaceData( X, Z, Y );

% Set up fittype and options.
ft = fittype( 'poly12' ); % << play around with this until an acceptable error is reached.

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

As for your concern regarding finding z=f(x,y): the solution depends on how many times you'd have to evaluate this expression, how often, and with which accuracy.

If done a small number of times, you can treat this like an equation-solving problem. Specifically, assuming we know y = f(x,z), x0 and y0 we'll end up with an expression like in the poly12 case:
```
y0 = p00 + p10*x0 + p01*z + p11*x0*z + p02*z^2
```
which, after rearranging, yields:
```
0 = (p02)*z^2 + (p01 + p1*x0)*z + (p00 + p10*x0 - y0)
```
Which is a simple zero-finding problem with one unknown. In the second degree polynomial case, the above has a closed solution, but in general, fzero should be able to handle it given the range of z ([0 100]).
If you need to do this many times, or for many points in one go, you could create fine grids for X and Z, evaluate Y using the polynomial surface equation, then perform 2d interpolation on your now-resampled grid, using either interp2 (or qinterp2) etc.

2D interpolation problem with scattered data

1 Answers