如何在极坐标中使用插值
How to interpolate using in polar coordinate
我有极坐标、半径0.05 <= r <= 1和0 ≤ θ ≤ 2π。半径r是0.05到1之间的50个值，极角θ是0到2π之间的24个值。
如何插入 r = 0.075 和 theta = pi/8？
根据评论您有以下信息
%the test point
ri=0.53224;
ti = pi/8;
%formula fo generation of Z
g=9.81
z0=@(r)0.01*(g^2)*((2*pi)^-4)*(r.^-5).*exp(-1.25*(r/0.3).^-4);
D=@(t)(2/pi)*cos(t).^2;
z2=@(r,t)z0(r).*D(t) ;
%range of vlaues of r and theta
r=[0.05,0.071175,0.10132,0.14422,0.2053, 0.29225,0.41602,0.5922,0.84299,1.2];
t=[0,0.62832,1.2566,1.885, 2.5133,3.1416,3.7699,4.3982,5.0265,5.6549,6.2832];
并且您想要测试点的插值。
当您对一些数据进行采样以将其用于插值时，您应该考虑如何根据您的要求对数据进行采样。
因此，当您对极坐标的规则网格进行采样时，这些坐标在转换为矩形时将形成一个圆形
大多数点都集中在 cricle 的中心，当我们从中心移动到外部区域时，点之间的距离增加。
%regular grid generated for r and t
[THETA R] = meshgrid(t ,r);
% Z for polar grid
Z=z2(R,THETA);
%convert coordinate from polar to cartesian(rectangular):
[X, Y] = pol2cart (THETA, R);
%plot points
plot(X, Y, 'k.');
axis equal
因此，当您使用这些点进行插值时，插值精度在中心较高，而在点之间距离增加的外部区域较低。
换句话说，使用这种采样方法，您会更加重视与外部区域相关的中心区域。
要提高精度，应增加网格点（r 和 theta）的密度，因此如果 r 和 theta 的长度为 11，则可以创建大小为 20 的 r 和 theta 以提高精度。
另一方面，如果您在 直角坐标 中创建规则网格，则每个区域都具有同等重要性。因此所有区域的插值精度都相同。
为此，您首先在极坐标中创建一个规则网格，然后将网格转换为直角坐标，以便您可以计算直角坐标中采样点的范围（最小最大）。基于此范围，您可以在直角坐标中创建规则网格
然后将直角坐标的规则网格转换为极坐标，使用 z2 公式得到网格点的 z。
%get the extent of points
extentX = [min(X(:)) max(X(:))];
extentY = [min(Y(:)) max(Y(:))];
%sample 100 points(or more or less) inside a region specified be the extents
X_samples = linspace(extentX(1),extentX(2),100);
Y_samples = linspace(extentY(1),extentY(2),100);
%create regular grid in rectangular coordinates
[XX YY] = meshgrid(X_samples, Y_samples);
[TT RR] = cart2pol(XX,YY);
Z_rect = z2(RR,TT);
对于测试点的插值说 [ri ti] 首先它转换为矩形然后使用 XX ,YY 值 z 被插值
[xi yi] = pol2cart (ti, ri);
z=interp2(XX,YY,Z_rect,xi,yi);
如果您没有选择更改数据采样方式并且只有极点网格（如与@RodyOldenhuis 讨论的那样），您可以执行以下操作：
用interp2插值极坐标（网格化数据的插值）
这种方法很简单，但缺点是 r 和 theta 的尺度不同，这可能会影响插值的准确性。
z = interp2(THETA, R, Z, ti, ri)
将极坐标转换为直角坐标，然后应用针对散点数据的插值方法。
这种方法需要更多的计算，但它的结果更可靠。
MATLAB 有 griddata 函数，给定的散点首先生成点的三角剖分，然后在三角形顶部创建规则网格并对网格点的值进行插值。
因此，如果您想对点 [ri ti] 的值进行插值，则应应用第二个插值以从插值网格中获取该点的值。
借助 spatialanalysisonline and Wikipedia 基于三角测量的线性插值的一些信息，以这种方式计算（在 Octave 中测试。在较新版本的 MATLAB 中，推荐使用 triangulation 和 pointLocation delaunay 和 tsearch 的）：
ri=0.53224;
ti = pi/8;
[THETA R] = meshgrid(t ,r);
[X, Y] = pol2cart (THETA, R);
[xi yi] = pol2cart (ti, ri);
%generate triangulation
tri = delaunay (X, Y);
%find the triangle that contains the test point
idx = tsearch (X, Y, tri, xi, yi);
pts= tri(idx,:);
%create a matrix that repesents equation of a plane (triangle) given its 3 points
m=[X(pts);Y(pts);Z(pts);ones(1,3)].';
%calculate z based on det(m)=0;
z= (-xi*det(m(:,2:end)) + yi*det([m(:,1) m(:,3:end)]) + det(m(:,1:end-1)))/det([m(:,1:2) m(:,end)]);
进一步细化：
由于已知搜索点被 4 个点包围，我们可以仅使用这些点进行三角测量。这些点形成一个梯形。梯形的每条对角线形成两个三角形，因此使用梯形的顶点我们可以形成 4 个三角形，梯形内的一个点也可以位于至少 2 个三角形中。
之前基于三角测量的方法仅使用来自一个三角形的信息，但是这里测试点的 z 可以从两个三角形的数据中插值两次，并且可以对计算出的 z 值进行平均以获得更好的近似值。
%find 4 points surrounding the test point
ft= find(t<=ti,1,'last');
fr= find(cos(abs(diff(t(ft+(0:1))))/2) .* r < ri,1,'last');
[T4 R4] = meshgrid(t(ft+(0:1)), r(fr+(0:1)));
[X4, Y4] = pol2cart (T4, R4);
Z4 = Z(fr+(0:1),ft+(0:1));
%form 4 triangles
tri2= nchoosek(1:4,3);
%empty vector of z values that will be interpolated from 4 triangles
zv = NaN(4,1);
for h = 1:4
pts = tri2(h,:);
% test if the point lies in the triangle
if ~isnan(tsearch(X4(:),Y4(:),pts,xi,yi))
m=[X4(pts) ;Y4(pts) ;Z4(pts); [1 1 1]].';
zv(h)= (-xi*det(m(:,2:end)) + yi*det([m(:,1) m(:,3:end)]) + det(m(:,1:end-1)))/det([m(:,1:2) m(:,end)]);
end
end
z= mean(zv(~isnan(zv)))
结果：
True z:
(0.0069246)
Linear Interpolation of (Gridded) Polar Coordinates :
(0.0085741)
Linear Interpolation with Triangulation of Rectangular Coordinates:
(0.0073774 or 0.0060992) based on triangulation
Linear Interpolation with Triangulation of Rectangular Coordinates(average):
(0.0067383)
结论：
与原始数据结构和采样方法有关的插值结果。如果采样方法与原始数据模式匹配插值结果更准确，那么在极坐标网格点遵循数据模式的情况下，规则极坐标插值结果会更可靠。但如果规则极坐标与数据结构不匹配或数据结构如不规则地形，基于三角剖分的插值方法可以更好地表示数据。
我不知道您尝试了什么，但是 interp2 在极坐标数据上的效果与在笛卡尔坐标上的效果一样好。这里有一些证据：
% Coordinates
r = linspace(0.05, 1, 50);
t = linspace(0, 2*pi, 24);
% Some synthetic data
z = sort(rand(50, 24));
% Values of interest
ri = 0.075;
ti = pi/8;
% Manually interpolate
rp = find(ri <= r, 1, 'first');
rm = find(ri >= r, 1, 'last');
tp = find(ti <= t, 1, 'first');
tm = find(ti >= t, 1, 'last');
drdt = (r(rp) - r(rm)) * (t(tp) - t(tm));
dr = [r(rp)-ri ri-r(rm)];
dt = [t(tp)-ti ti-t(tm)];
fZ = [z(rm, tm) z(rm, tp)
z(rp, tm) z(rp, tp)];
ZI_manual = (dr * fZ * dt.') / drdt
% Interpolate with MATLAB
ZI_MATLAB = interp2(r, t, z', ri, ti, 'linear')
结果：
ZI_manual =
2.737907208525297e-002
ZI_MATLAB =
2.737907208525298e-002
请检查这个例子，我用了两个for循环，在for循环里面我使用了条件语句，如果你评论这个条件语句和运行这个程序，你会得到正确的答案，你取消注释后条件语句和 运行 程序，你会得到错误的答案。请检查一下。
% Coordinates
r = linspace(0.05, 1, 10);
t = linspace(0, 2*pi, 8);
% Some synthetic data
%z = sort(rand(50, 24));
z=zeros();
for i=1:10
for j=1:8
if r(i)<0.5||r(i)>1
z(i,j)=0;
else
z(i,j) = r(i).^3'*cos(t(j)/2);
end
end
end
% Values of interest
ri = 0.55;
ti = pi/8;
% Manually interpolate
rp = find(ri <= r, 1, 'first');
rm = find(ri >= r, 1, 'last');
tp = find(ti <= t, 1, 'first');
tm = find(ti >= t, 1, 'last');
drdt = (r(rp) - r(rm)) * (t(tp) - t(tm));
dr = [r(rp)-ri ri-r(rm)];
dt = [t(tp)-ti ti-t(tm)];
fZ = [z(rm, tm) z(rm, tp)
z(rp, tm) z(rp, tp)];
ZI_manual = (dr * fZ * dt.') / drdt
% Interpolate with MATLAB
ZI_MATLAB = interp2(r, t, z', ri, ti, 'linear')
结果：
z1 =
0.1632
ZI_manual =
0.1543
ZI_MATLAB =
0.1582
我有极坐标、半径0.05 <= r <= 1和0 ≤ θ ≤ 2π。半径r是0.05到1之间的50个值，极角θ是0到2π之间的24个值。
如何插入 r = 0.075 和 theta = pi/8？
根据评论您有以下信息
%the test point
ri=0.53224;
ti = pi/8;
%formula fo generation of Z
g=9.81
z0=@(r)0.01*(g^2)*((2*pi)^-4)*(r.^-5).*exp(-1.25*(r/0.3).^-4);
D=@(t)(2/pi)*cos(t).^2;
z2=@(r,t)z0(r).*D(t) ;
%range of vlaues of r and theta
r=[0.05,0.071175,0.10132,0.14422,0.2053, 0.29225,0.41602,0.5922,0.84299,1.2];
t=[0,0.62832,1.2566,1.885, 2.5133,3.1416,3.7699,4.3982,5.0265,5.6549,6.2832];
并且您想要测试点的插值。
当您对一些数据进行采样以将其用于插值时，您应该考虑如何根据您的要求对数据进行采样。 因此，当您对极坐标的规则网格进行采样时，这些坐标在转换为矩形时将形成一个圆形 大多数点都集中在 cricle 的中心，当我们从中心移动到外部区域时，点之间的距离增加。
%regular grid generated for r and t
[THETA R] = meshgrid(t ,r);
% Z for polar grid
Z=z2(R,THETA);
%convert coordinate from polar to cartesian(rectangular):
[X, Y] = pol2cart (THETA, R);
%plot points
plot(X, Y, 'k.');
axis equal
因此，当您使用这些点进行插值时，插值精度在中心较高，而在点之间距离增加的外部区域较低。
换句话说，使用这种采样方法，您会更加重视与外部区域相关的中心区域。
要提高精度，应增加网格点（r 和 theta）的密度，因此如果 r 和 theta 的长度为 11，则可以创建大小为 20 的 r 和 theta 以提高精度。
另一方面，如果您在 直角坐标 中创建规则网格，则每个区域都具有同等重要性。因此所有区域的插值精度都相同。
为此，您首先在极坐标中创建一个规则网格，然后将网格转换为直角坐标，以便您可以计算直角坐标中采样点的范围（最小最大）。基于此范围，您可以在直角坐标中创建规则网格
然后将直角坐标的规则网格转换为极坐标，使用 z2 公式得到网格点的 z。
%get the extent of points
extentX = [min(X(:)) max(X(:))];
extentY = [min(Y(:)) max(Y(:))];
%sample 100 points(or more or less) inside a region specified be the extents
X_samples = linspace(extentX(1),extentX(2),100);
Y_samples = linspace(extentY(1),extentY(2),100);
%create regular grid in rectangular coordinates
[XX YY] = meshgrid(X_samples, Y_samples);
[TT RR] = cart2pol(XX,YY);
Z_rect = z2(RR,TT);
对于测试点的插值说 [ri ti] 首先它转换为矩形然后使用 XX ,YY 值 z 被插值
[xi yi] = pol2cart (ti, ri);
z=interp2(XX,YY,Z_rect,xi,yi);
如果您没有选择更改数据采样方式并且只有极点网格（如与@RodyOldenhuis 讨论的那样），您可以执行以下操作：
用
interp2插值极坐标（网格化数据的插值） 这种方法很简单，但缺点是r和 theta 的尺度不同，这可能会影响插值的准确性。z = interp2(THETA, R, Z, ti, ri)将极坐标转换为直角坐标，然后应用针对散点数据的插值方法。 这种方法需要更多的计算，但它的结果更可靠。 MATLAB 有
griddata函数，给定的散点首先生成点的三角剖分，然后在三角形顶部创建规则网格并对网格点的值进行插值。 因此，如果您想对点[ri ti]的值进行插值，则应应用第二个插值以从插值网格中获取该点的值。
借助 spatialanalysisonline and Wikipedia 基于三角测量的线性插值的一些信息，以这种方式计算（在 Octave 中测试。在较新版本的 MATLAB 中，推荐使用 triangulation 和 pointLocation delaunay 和 tsearch 的）：
ri=0.53224;
ti = pi/8;
[THETA R] = meshgrid(t ,r);
[X, Y] = pol2cart (THETA, R);
[xi yi] = pol2cart (ti, ri);
%generate triangulation
tri = delaunay (X, Y);
%find the triangle that contains the test point
idx = tsearch (X, Y, tri, xi, yi);
pts= tri(idx,:);
%create a matrix that repesents equation of a plane (triangle) given its 3 points
m=[X(pts);Y(pts);Z(pts);ones(1,3)].';
%calculate z based on det(m)=0;
z= (-xi*det(m(:,2:end)) + yi*det([m(:,1) m(:,3:end)]) + det(m(:,1:end-1)))/det([m(:,1:2) m(:,end)]);
进一步细化： 由于已知搜索点被 4 个点包围，我们可以仅使用这些点进行三角测量。这些点形成一个梯形。梯形的每条对角线形成两个三角形，因此使用梯形的顶点我们可以形成 4 个三角形，梯形内的一个点也可以位于至少 2 个三角形中。 之前基于三角测量的方法仅使用来自一个三角形的信息，但是这里测试点的 z 可以从两个三角形的数据中插值两次，并且可以对计算出的 z 值进行平均以获得更好的近似值。
%find 4 points surrounding the test point
ft= find(t<=ti,1,'last');
fr= find(cos(abs(diff(t(ft+(0:1))))/2) .* r < ri,1,'last');
[T4 R4] = meshgrid(t(ft+(0:1)), r(fr+(0:1)));
[X4, Y4] = pol2cart (T4, R4);
Z4 = Z(fr+(0:1),ft+(0:1));
%form 4 triangles
tri2= nchoosek(1:4,3);
%empty vector of z values that will be interpolated from 4 triangles
zv = NaN(4,1);
for h = 1:4
pts = tri2(h,:);
% test if the point lies in the triangle
if ~isnan(tsearch(X4(:),Y4(:),pts,xi,yi))
m=[X4(pts) ;Y4(pts) ;Z4(pts); [1 1 1]].';
zv(h)= (-xi*det(m(:,2:end)) + yi*det([m(:,1) m(:,3:end)]) + det(m(:,1:end-1)))/det([m(:,1:2) m(:,end)]);
end
end
z= mean(zv(~isnan(zv)))
结果：
True z:
(0.0069246)
Linear Interpolation of (Gridded) Polar Coordinates :
(0.0085741)
Linear Interpolation with Triangulation of Rectangular Coordinates:
(0.0073774 or 0.0060992) based on triangulation
Linear Interpolation with Triangulation of Rectangular Coordinates(average):
(0.0067383)
结论：
与原始数据结构和采样方法有关的插值结果。如果采样方法与原始数据模式匹配插值结果更准确，那么在极坐标网格点遵循数据模式的情况下，规则极坐标插值结果会更可靠。但如果规则极坐标与数据结构不匹配或数据结构如不规则地形，基于三角剖分的插值方法可以更好地表示数据。
我不知道您尝试了什么，但是 interp2 在极坐标数据上的效果与在笛卡尔坐标上的效果一样好。这里有一些证据：
% Coordinates
r = linspace(0.05, 1, 50);
t = linspace(0, 2*pi, 24);
% Some synthetic data
z = sort(rand(50, 24));
% Values of interest
ri = 0.075;
ti = pi/8;
% Manually interpolate
rp = find(ri <= r, 1, 'first');
rm = find(ri >= r, 1, 'last');
tp = find(ti <= t, 1, 'first');
tm = find(ti >= t, 1, 'last');
drdt = (r(rp) - r(rm)) * (t(tp) - t(tm));
dr = [r(rp)-ri ri-r(rm)];
dt = [t(tp)-ti ti-t(tm)];
fZ = [z(rm, tm) z(rm, tp)
z(rp, tm) z(rp, tp)];
ZI_manual = (dr * fZ * dt.') / drdt
% Interpolate with MATLAB
ZI_MATLAB = interp2(r, t, z', ri, ti, 'linear')
结果：
ZI_manual =
2.737907208525297e-002
ZI_MATLAB =
2.737907208525298e-002
请检查这个例子，我用了两个for循环，在for循环里面我使用了条件语句，如果你评论这个条件语句和运行这个程序，你会得到正确的答案，你取消注释后条件语句和 运行 程序，你会得到错误的答案。请检查一下。
% Coordinates
r = linspace(0.05, 1, 10);
t = linspace(0, 2*pi, 8);
% Some synthetic data
%z = sort(rand(50, 24));
z=zeros();
for i=1:10
for j=1:8
if r(i)<0.5||r(i)>1
z(i,j)=0;
else
z(i,j) = r(i).^3'*cos(t(j)/2);
end
end
end
% Values of interest
ri = 0.55;
ti = pi/8;
% Manually interpolate
rp = find(ri <= r, 1, 'first');
rm = find(ri >= r, 1, 'last');
tp = find(ti <= t, 1, 'first');
tm = find(ti >= t, 1, 'last');
drdt = (r(rp) - r(rm)) * (t(tp) - t(tm));
dr = [r(rp)-ri ri-r(rm)];
dt = [t(tp)-ti ti-t(tm)];
fZ = [z(rm, tm) z(rm, tp)
z(rp, tm) z(rp, tp)];
ZI_manual = (dr * fZ * dt.') / drdt
% Interpolate with MATLAB
ZI_MATLAB = interp2(r, t, z', ri, ti, 'linear')
结果： z1 = 0.1632 ZI_manual = 0.1543 ZI_MATLAB = 0.1582