Задача 2D-интерполяции с разбросанными данными - PullRequest
Новая
       73

Задача 2D-интерполяции с разбросанными данными

2 голосов
/ 03 июля 2019

У меня есть функция F, которая позволяет мне вычислять 'y' на основе значений x и z. то есть y = F (x, z). Эта функция очень интенсивна в вычислительном отношении; на одну оценку требуется около 10 минут.

Что мне действительно нужно, так это найти значение «z» по паре значений (x, y), т. Е. Мне нужно z = G (x, y). Решение F невозможно.

Чтобы численно реализовать G, я прошел и рассчитал много значений y, используя F (x, z), используя сетку по x и z.

Прилагается график значений z на контурном графике. Красные точки представляют точки (x, y), для которых у меня есть значения z. Трудно увидеть настолько уменьшенное изображение, поскольку расстояние между этими точками нерегулярно. Справа также график сёрфинга, показывающий, что G (x, y) сглаживается и ведет себя хорошо. 

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

Однако, когда я пытаюсь интерполировать эту поверхность, гладкость поверхности нарушается, и я получаю много шума. Это происходит, если я использую griddata или scatteredInterpolant. Я пробовал разные методы: «линейный», «натуральный», «кубический» и т. Д., И все это одинаково (Извините за некоторые полосы на рисунках. У меня сейчас проблемы с видеокартами.) 

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

Кто-нибудь знает, что происходит? Я не могу использовать обычную интерполяцию данных с разнесенными точками, потому что мои точки y расположены не совсем равномерно. Ниже приведен график, где я проверил треугольную сетку, и она выглядит хорошо. 

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

Для полноты X_mesh, Y_mesh и Z_mesh, выводимые на консоль, приведены ниже. Все три матрицы 19x51 двойные.

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

Обновление: Я нашел ключ к проблеме. Ниже приведен график исходных (неинтерполированных) данных с включенной функцией затенения с использованием графиков "surf" и "trisurf". Surf производит довольно гладкую поверхность, тогда как с полосами trisurf начинают появляться. Я полагаю, что это те же полосы, что и для griddata или scatteredInterpolant, в котором используется треугольная сетка.

Кто-нибудь знает, как я могу получить доступ к интерполированным значениям из поверхностного объекта?

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

1 Ответ

1 голос
/ 04 июля 2019

Вы, кажется, не используете правильный инструмент для задания . Сетчатый или линейный интерполант не является тем, что требуется для этой проблемы. Правильный подход здесь упоминается в комментарии Cris - то есть, чтобы повернуть набор данных или просто изменить порядок измерений.

Играя с данными, которые вы разместили (первые 12 столбцов) в cftool Я могу вам сказать, что y=f(x,z) можно выразить почти идеально, используя полиномиальную поверхность степени 1-3 обоих параметров. Вот некоторые результаты, например:

  • poly32 модель (9 коэффициентов) 
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 модель (5 коэффициентов) 
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

или программно, 

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

Что касается вашего беспокойства по поводу поиска z=f(x,y): решение зависит от того, сколько раз вам придется оценивать это выражение, как часто и с какой точностью.

  • Если сделать это небольшое количество раз, вы можете рассматривать это как задачу решения уравнений. В частности, предполагая, что мы знаем y = f(x,z), x0 и y0, мы получим выражение, как в случае poly12: 

    y0 = p00 + p10*x0 + p01*z + p11*x0*z + p02*z^2

    который после перестановки дает: 

    0 = (p02)*z^2 + (p01 + p1*x0)*z + (p00 + p10*x0 - y0)

    Что представляет собой простую задачу поиска нуля с одним неизвестным. В случае полинома второй степени вышеупомянутое имеет закрытое решение, но в общем случае fzero должно быть в состоянии справиться с ним, учитывая диапазон z ([0 100]).

  • Если вам нужно сделать это много раз или для множества точек за один раз, вы можете создать точные сетки для X и Z, оценить Y с использованием уравнения для полиномиальной поверхности, а затем выполнить 2d интерполяция в вашей теперь измененной сетке, используя interp2 (или qinterp2) и т. д.

...