DEVHIDE
Login Or Sign up

2D interpolation problem with scattered data

669 Views Asked by At

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');

Contour and surf plot of the input data

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');

Attempt to interpolate data

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');

Delaunay Triangular mesh

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?

Comparing surf and trisurf of the original (uninterpolated) data with shading interp turned on

matlab interpolation surface smoothing gridding
Original Q&A
1

There are 1 best solutions below

0
On

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.

Related Questions in MATLAB

Related Questions in INTERPOLATION

Related Questions in SURFACE

Related Questions in SMOOTHING

Related Questions in GRIDDING

Trending Questions

Popular # Hahtags

javascript python java c# php android html jquery c++ css ios sql mysql r reactjs node.js arrays c asp.net json python-3.x ruby-on-rails .net sql-server swift django angular objective-c pandas excel

Popular Questions

Copyright © 2021 Jogjafile Inc.