如何內插使用極座標

Question

我有極座標，半徑爲0.05 <= r <= 1和0 ≤ θ ≤ 2π。半徑r是在0.05到1之間的50個值，極角θ是0到2π之間的24個值。

如何內插r = 0.075和theta = pi/8？

Answer 1

基於評論你有以下信息

%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，則可以創建r和theta，大小爲20，以提高精度。
另一方面，如果您在的直角座標中創建了一個規則的網格，則對每個區域都賦予相同的重要性。所以插值的準確性在所有地區都是一樣的。
首先在極座標中創建一個規則網格，然後將網格轉換爲直角座標，以便可以計算直角座標中採樣點的範圍（最小值）。基於此範圍，您可以在直角座標系中創建規則網格 直角座標系的常規網格然後轉換爲極座標系，從而使用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]先將其轉換爲矩形，然後使用zXX ,YY值內

[xi yi] = pol2cart (ti, ri); 
z=interp2(XX,YY,Z_rect,xi,yi); 

如果你沒有選擇，改變你的採樣數據，只有具有與@RodyOldenhuis討論的極點網格，您可以執行以下操作：

1. In使用interp2（網格數據插值）的terpolate極座標 這種方法非常簡單，但缺點是r和θ不是相同的比例尺，這可能會影響插值的準確性。

   z = interp2(THETA, R, Z, ti, ri)

2. 轉換極座標矩形，然後應用插值方法是散亂數據。這種方法需要更多的計算，但結果更可靠。 MATLAB有griddata函數，給定的散點首先生成點的三角剖分，然後在三角形頂部創建一個規則網格並插入網格點的值。 因此，如果要插入點[ri ti]的值，則應該應用第二次插值來從插值網格中獲取點的值。

隨着從spatialanalysisonline和Wikipedia線性插值根據三角測量的一些信息來計算這種方式幫助（在八度測試在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) 

結論：

與原始數據的結構和採樣方法的內插結果。如果採樣方法與原始數據匹配模式的插值結果更加準確，那麼在極座標網格點跟隨數據模式的情況下，正常極座標插值結果可以更加可靠。但是，如果常規的極座標與數據結構或數據結構不匹配，例如不規則地形，則基於三角測量的插值方法可以更好地表示數據。

Answer 2

我不知道你曾經嘗試過什麼，但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 

Answer 3

請這個例子中，我對循環使用的兩個，在for循環中，我用條件語句，如果u評論這個條件語句，並運行程序，你會得到正確的答案，在你取消註釋這個條件語句並運行程序後，你會得到錯誤的答案。請檢查一下。

% 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