scatteredInterpolant: what is linear interpolation in 2d?

Question

SA-W 2023-7-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1990648-scatteredinterpolant-what-is-linear-interpolation-in-2d

评论： Torsten 2023-7-5

%coords
supportPts = [3 3; 3.3 3; 3 3.25; 3.3 3.25; 3.6 3; 3.6 3.25; 3 3.5; 3.3 3.5; 3.6 3.5];
%values
Fval = [0 0.1121 0.064604 0.184942 0.233029 0.352622 0.121444 0.206496 0.375685]';
%interpolation object
Interp = scatteredInterpolant(supportPts(:,1),supportPts(:,2),Fval);
%evaluate at center of bottom left element
Interp(3.15, 3.125)
ans = 0.0884

See the above example with nine points that represent four axis-parrallel elements.

I would have expected that the value of the interpoland at the center of the bottom left element is the mean value of the coresponding four corner values, i.e.,

%expected value at (3.15,3.125)
0.25*(0+0.064604+0.184942+0.1121)
ans = 0.0904

but this is clearly not so.

So what interpolation scheme is scatteredInterpoland using (linear) to calculate values in the interior of the elements? At the element boundaries, I double-checked that there is just 1d interpolation involved.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2023-7-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1990648-scatteredinterpolant-what-is-linear-interpolation-in-2d#answer_1265783

编辑：Matt J 2023-7-1

No, the support points are being divided into triangular elements, not rectangles. It is puzzling that you would be using scatteredInterpolant instead of griddedInterpolant for data like this.

31 个评论
显示 29更早的评论隐藏 29更早的评论

Matt J 2023-7-2

编辑：Matt J 2023-7-2

在 MATLAB Online 中打开

Say, supportPts is my input (nx2) matrix which contains n pairs of x and y coordinates in arbitrary order that define a rectangular grid.

As Bruno says, you would sort them using sortrows() and then unique().

[P,is]=sortrows(supportPts);

x=unique(P(:,1));

y=unique(P(:,2));

Fval=reshape(Fval(is),numel(x),numel(y));

Interp=griddedInterpolant({x,y},Fval)

But why would you even be starting with data in such a form? If it's a rectangular grid, you would usually be given the sampling of x and y in separate vectors from the beginning, which is the input format that meshgrid, ndgrid, and various other grid operation commands require.

[X,Y]=ndgrid(3:0.3:3.6  , 3:0.25:3.5), Z=X.^2+Y.^2
X = 3×3
    3.0000    3.0000    3.0000
    3.3000    3.3000    3.3000
    3.6000    3.6000    3.6000
Y = 3×3
    3.0000    3.2500    3.5000
    3.0000    3.2500    3.5000
    3.0000    3.2500    3.5000
Z = 3×3
   18.0000   19.5625   21.2500
   19.8900   21.4525   23.1400
   21.9600   23.5225   25.2100

SA-W 2023-7-3

编辑：SA-W 2023-7-3

and the default is bilienar interpolation.

I expected a different output from griddedInterpolant when comparing the results to those of my visualization tool (paraview).

In my finite-element framework, I have bilinear interpolation functions

N_0 = 0.25*(1-u_1)*(1-u_2)

N_1 = 0.25*(1+u_1)*(1-u_2)

N_2 = 0.25*(1+u_1)*(1+u_2)

N_3 = 0.25*(1-u_i)*(1+u_2)

where u_1,u_2 \in [-1,1] are reference coordinates. If you think of a rectangle with corners [-1,-1], [1,-1], [1,1], [-1,1] , then the N_i are interpolation functions. For instance, N_0 is one at [-1,-1] and zero at the remaining three vertices. So if we have a point [x,y] of a rectangular element, this point is mapped into the reference coordinates [u_1,u_2] and the interpolated value is obtained as

val = N_i * E_i (i=1,...,4),

where the E_i are the local vertex values of the element where the point [x,y] lives.

This is clearly a bilinear interpolation and I expected this to be the same what griddedInterpolant does with linear method. However, if I evaluate the interpoland from my question at the point [3.2,3.2], I obtain 0.13714 in matlab and 0.13081 in paraview. At the element centers, the output is the same, also along the element "boundaries", only in the interior of the elements there is a mismatch.

Given that, is griddedInterpolant doing something different than what I described?

SA-W 2023-7-3

编辑：SA-W 2023-7-3

在 MATLAB Online 中打开

@Bruno Luong @Matt J

If you find MATLAB does not return this let us know.

I managed to create an example that illustrates the differences I see between griddedInterpolant and the formulae (sum of basis functions times nodal values):

%coordinates
supportPts = [3 3; 3.3 3; 3 3.25; 3.3 3.25; 3.6 3; 3.6 3.25; 3 3.5; 3.3 3.5; 3.6 3.5];
%nodal values
Fval = [0 0.1121 0.064604 0.184942 0.233029 0.352622 0.121444 0.206496 0.375685]';
%gridded interpolation object
[P,is] = sortrows(supportPts); 
Inv1 = unique(P(:,1));
Inv2 = unique(P(:,2));
Fval = reshape(Fval(is),numel(Inv1),numel(Inv2));
Interp = griddedInterpolant({Inv1,Inv2},Fval);
%vertices and corresponding values of bottom left element
vertices = [3 3; 3.3 3; 3.3 3.25; 3 3.25];
vertexValues = [0 0.1121 0.184942 0.064604];
%point where interpoland should be evaluated (currently center)
ptReal = [3.15 3.125];
%determine the coordinates on the reference element (-1,1)x(-1,1)
fun = @(x)myFun(x, ptReal, vertices);
ptRef = fsolve(fun, [1 1]);
Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
disp('value formulae with basis functions: ');
value formulae with basis functions: 
fEval(ptRef, vertexValues)
ans = 0.0904
disp('value griddedInterpolant: ');
value griddedInterpolant: 
Interp(ptReal(1), ptReal(2))
ans = 0.0904
%function to determine reference coordinates of ptEval
function F = myFun(x, ptEval, vertexValues)
N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));
F = N_0.*vertexValues(1,:) + N_1.*vertexValues(2,:) + ...
    N_2.*vertexValues(3,:) + N_3.*vertexValues(4,:) - ptEval;
end
%function to interpolate using local basis functions
function fval = fEval(x, values)
N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));
N = [N_0 N_1 N_2 N_3];
fval = N*values';
end

As you can see, both methods return the same value at the element center. But if you define

ptReal = [3.2 3.2];

we get

%value formulae with basis functions: 
ans =
0.056050010755214
%value griddedInterpolant: 
ans =
0.032302000000000

0.05605 is also what I see in my visualization tool. Note that et the center (ptReal = [3.15 3.125]), both methods give the same value.

So where is the "error" in my calculation?

Stephen23 2023-7-4

编辑：Stephen23 2023-7-4

在 MATLAB Online 中打开

You have accidentally transposed your data. This is because SORTROWS does not return the data in NDGRID format as GRIDDEDINTERPOLANT requires when using grid vectors as you did:

https://www.mathworks.com/help/matlab/ref/griddedinterpolant.html#bvh2f0p-5

First lets have a look at an NDGRID example, to see what order it has:

[X,Y] = ndgrid([3,3.3,3.6],[3,3.25,3.5,Inf])
X = 3×4
    3.0000    3.0000    3.0000    3.0000
    3.3000    3.3000    3.3000    3.3000
    3.6000    3.6000    3.6000    3.6000
Y = 3×4
    3.0000    3.2500    3.5000       Inf
    3.0000    3.2500    3.5000       Inf
    3.0000    3.2500    3.5000       Inf

Note the linear order of elements in X is: 3, 3.3, 3.6, 3, 3.3... etc. Now lets try SORTROWS:

supportPts = [3,3;3.3,3;3,3.25;3.3,3.25;3.6,3;3.6,3.25;3,3.5;3.3,3.5;3.6,3.5];
Fval = [0;0.1121;0.064604;0.184942;0.233029;0.352622;0.121444;0.206496;0.375685];
[P,is] = sortrows(supportPts) % is P in NDGRID order? (hint: no)
P = 9×2
    3.0000    3.0000
    3.0000    3.2500
    3.0000    3.5000
    3.3000    3.0000
    3.3000    3.2500
    3.3000    3.5000
    3.6000    3.0000
    3.6000    3.2500
    3.6000    3.5000
is = 9×1
     1
     3
     7
     2
     4
     8
     5
     6
     9

From this point on all results are rubbish, because P is clearly not in NDGRID order and therefore neither will FVAL be when you use those indices. You have effectively swapped the dim1 and dim2, transposing your data. This is easy to see when you plot the data (which you should always do!):

scatter3(supportPts(:,1),supportPts(:,2),Fval(:),99,'g','filled')

hold on

Inv1 = unique(P(:,1));

Inv2 = unique(P(:,2));

Fone = reshape(Fval(is),numel(Inv1),numel(Inv2))

Fone = 3×3

0 0.1121 0.2330 0.0646 0.1849 0.3526 0.1214 0.2065 0.3757

[M1,M2] = ndgrid(Inv1,Inv2) % expected order of GRIDDEDINTERPOLANT

M1 = 3×3

3.0000 3.0000 3.0000 3.3000 3.3000 3.3000 3.6000 3.6000 3.6000

M2 = 3×3

3.0000 3.2500 3.5000 3.0000 3.2500 3.5000 3.0000 3.2500 3.5000

scatter3(M1(:),M2(:),Fone(:),42,'r','filled') % Ooops, swapped dimensions!

hold off

You can see that the data are reflected X-Y. A simple solution is to specify the SORTROWS column order:

[P,is] = sortrows(supportPts,[2,1]);
Inv1 = unique(P(:,1));
Inv2 = unique(P(:,2));
Ftwo = reshape(Fval(is),numel(Inv1),numel(Inv2));
Interp = griddedInterpolant({Inv1,Inv2},Ftwo);

Lets plot those data:

scatter3(supportPts(:,1),supportPts(:,2),Fval(:),99,'g','filled')

hold on

scatter3(M1(:),M2(:),Ftwo(:),42,'r','filled') % Aaaah, much better!

hold off

Now try your various interpolation points:

vertices = [3 3; 3.3 3; 3.3 3.25; 3 3.25];
vertexValues = [0 0.1121 0.184942 0.064604];
pt1 = [3.15,3.125];
fun = @(x)myFun(x, pt1, vertices);
ptRef = fsolve(fun, [1 1]);
Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
fEval(ptRef, vertexValues)
ans = 0.0904
val1 = Interp(pt1(1), pt1(2))
val1 = 0.0904
pt2 = [3.2,3.2];
fun = @(x)myFun(x, pt2, vertices);
ptRef = fsolve(fun, [1 1]);
Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
fEval(ptRef, vertexValues)
ans = 0.1308
val2 = Interp(pt2(1), pt2(2))
val2 = 0.1308

Alternatively you could avoid those grid vectors entirely:

obj2 = griddedInterpolant(...
    reshape(P(:,1),3,3),... or M1 & M2
    reshape(P(:,2),3,3),Ftwo);
obj2(pt1(1),pt1(2))
ans = 0.0904
obj2(pt2(1),pt2(2))
ans = 0.1308

And finally checking against another routine:

interpn(M1,M2,Ftwo,pt1(1),pt1(2))
ans = 0.0904
interpn(M1,M2,Ftwo,pt2(1),pt2(2))
ans = 0.1308

"0.05605 is also what I see in my visualization tool."

Note that the value for [3.2,3.2] cannot be less than that for [3.15,3.125] because you move closer to larger values along both axes, so stating that the expected value is smaller simply makes no sense. Rather than blindly copying results from a tool, you should always perform some basic sanity checks yourself.

Lets plot them. Try rotating and see how they line up (unfortunately this is not (yet) interactive on this forum). I also added a contour to show where all solutions == 0.05605 lie, i.e. nowhere near [3.2,3.2].

mesh(M1,M2,Ftwo,'EdgeColor','k','FaceColor','none')

hold on

contour3(M1,M2,Ftwo,[0.05605,0.05605],'-r') % == 0.05605

scatter3(supportPts(:,1),supportPts(:,2),Fval(:),99,'g','filled')

scatter3(pt1(1),pt1(2),val1,64,'b','filled')

scatter3(pt2(1),pt2(2),val2,64,'k','filled')

So far everything seems to be working as expected. Unmodified support functions:

function F = myFun(x, ptEval, vertexValues)
N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));
F = N_0.*vertexValues(1,:) + N_1.*vertexValues(2,:) + ...
    N_2.*vertexValues(3,:) + N_3.*vertexValues(4,:) - ptEval;
end
%function to interpolate using local basis functions
function fval = fEval(x, values)
N_0 = 0.25*(1-x(1))*(1-x(2));
N_1 = 0.25*(1+x(1))*(1-x(2));
N_2 = 0.25*(1+x(1))*(1+x(2));
N_3 = 0.25*(1-x(1))*(1+x(2));
N = [N_0 N_1 N_2 N_3];
fval = N*values';
end

Lesson: always always always look at your data.

Bruno Luong 2023-7-4

编辑：Bruno Luong 2023-7-4

在 MATLAB Online 中打开

Here is the code to check the bilinear formula. No need for fsolve.

vertices = [3 3; 3.3 3; 3.3 3.25; 3 3.25];
vertexValues = [0 0.1121 0.184942 0.064604];
pt1 = [3.15,3.125];
pt2 = [3.2,3.2];
[XY,is] = sortrows(vertices, [2 1]);
x = unique(XY(:,1));
y = unique(XY(:,2));
nx = length(x);
ny = length(y);
sz = [nx ny];
X = reshape(XY(:,1),sz);
Y = reshape(XY(:,2),sz);
Z = reshape(vertexValues(is),sz);
f=griddedInterpolant(X,Y,Z);
format long
f(pt1(1),pt1(2))
ans = 
   0.090411500000000
bilinearformula(x, y, Z, pt1(1), pt1(2))
ans = 
   0.090411500000000
f(pt2(1),pt2(2))
ans = 
   0.130810133333333
bilinearformula(x, y, Z, pt2(1), pt2(2))
ans = 
   0.130810133333333
%% No check for overflowed so the formula is clearer
function zq = bilinearformula(x, y, Z, xq, yq)
i = discretize(xq,x);
j = discretize(yq,y);
dx = x(i+1)-x(i);
dy = y(j+1)-y(j);
wx2 = (xq-x(i))./dx;
wx1 = 1-wx2;
wy2 = (yq-y(j))./dy;
wy1 = 1-wy2;
B = [wx1.*wy1, wx1.*wy2, wx2.*wy2, wx2.*wy1]; % basis functions evaluate at (xq,yq)
Zij = [Z(i,j); Z(i,j+1); Z(i+1,j+1); Z(i+1,j)];
zq = B * Zij;
end

Bruno Luong 2023-7-5

编辑：Bruno Luong 2023-7-5

@SA-W I think the basis functions are not wrong, but finding the reference coordinates of a point in real space and working just with the basis functions on (-1,1)x(-1,1) is an alternative to your bilinearformula(x, y, Z, xq, yq), right?

Right. When I make the comment I overlook the fsolve part to map the rectangle to (-1,1)x(-1x1).

I would never do that. You want to know the formula for bilinear interpolation to understand how it works, and then you use a numerical function that does not provide any formula. Fortunately @Stephen23 give you a big help to figure out the problem of your code.

But more interesting, I just think you approach could be applied where the grid is topological equivalent to generic quadrilateral mesh, such as pixels of some non-linear camera projection (fisheye) and not grided data.

BTW for rectangular gridded the cross-term x*y of the mapping vanish, and essentially the solution of fsolve can be determined separately for x and y, and if you plug the formula together you find my formula, consist, compact and fast to compute.

For quatrilateral I would look for analytical inversion rather than a big fsolve hammer, if it exists.

Torsten 2023-7-5

在 MATLAB Online 中打开

syms x_left x_right y_low y_high

syms u_left_low u_right_low u_left_high u_right_high

syms x y

% First way

% Linearly interpolate between (x_left,y_low,u_left_low) and (x_right,y_low,u_right_low);

f_low(x) = u_left_low*(x-x_right)/(x_left-x_right) + u_right_low*(x-x_left)/(x_right-x_left);

% Linearly interpolate between (x_left,y_high,u_left_high) and (x_right,y_high,u_right_high);

f_high(x) = u_left_high*(x-x_right)/(x_left-x_right) + u_right_high*(x-x_left)/(x_right-x_left);

% Linearly interpolate between f_low(x) and f_high(x)

f1_inter(x,y) = f_low(x)*(y-y_high)/(y_low-y_high) + f_high(x)*(y-y_low)/(y_high-y_low)

f1_inter(x, y) =

% Second way

% Linearly interpolate between (x_left,y_low,u_left_low) and (x_left,y_high,u_left_high);

f_left(y) = u_left_low*(y-y_high)/(y_low-y_high) + u_left_high*(y-y_low)/(y_high-y_low);

% Linearly interpolate between (x_right,y_low,u_right_low) and (x_right,y_high,u_right_high);

f_right(y) = u_right_low*(y-y_high)/(y_low-y_high) + u_right_high*(y-y_low)/(y_high-y_low);

% Linearly interpolate between f_left(y) and f_right(y)

f2_inter(x,y) = f_left(y)*(x-x_right)/(x_left-x_right) + f_right(y)*(x-x_left)/(x_right-x_left)

f2_inter(x, y) =

x_left = 3;

x_right = 3.3;

y_low = 3;

y_high = 3.25;

u_left_low = 0;

u_right_low = 0.1121;

u_left_high = 0.064604;

u_right_high = 0.184942;

xq = 3.15;

yq = 3.125;

double(subs(f1_inter(xq,yq)))

ans = 0.0904

double(subs(f2_inter(xq,yq)))

ans = 0.0904

xq = 3.2;

yq = 3.2;

double(subs(f1_inter(xq,yq)))

ans = 0.1308

double(subs(f2_inter(xq,yq)))

ans = 0.1308

请先登录，再进行评论。

scatteredInterpolant: what is linear interpolation in 2d?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

31 个评论
显示 29更早的评论隐藏 29更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

scatteredInterpolant: what is linear interpolation in 2d?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

31 个评论 显示 29更早的评论隐藏 29更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

31 个评论
显示 29更早的评论隐藏 29更早的评论