如何在极坐标中进行插值

3 回答

• 2

  根据评论，您有以下信息

  %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和θ，以提高精度 .
  另一方面，如果在直角坐标中创建规则网格，则每个区域都具有相同的重要性 . 因此，插值的准确性在所有区域都是相同的 .
  首先，在极坐标中创建一个规则网格，然后将网格转换为直角坐标，这样就可以计算出直角坐标中采样点的范围（最小最大值） . 基于此范围，您可以在直角坐标中创建规则网格，然后使用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);
  

  对于测试点的插值，首先将它转换为矩形，然后使用 XX ,YY 的 z 值进行插值

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

  如果您无法更改数据采样方式，并且只有@RodyOldenhuis讨论的极点网格，则可以执行以下操作：

  • 使用 interp2 （网格化数据的插值）插入极坐标这种方法很简单，但缺点是 r 和θ不是相同的比例，这可能会影响插值的精度 .

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

  • 将极坐标转换为矩形，然后应用散乱数据的插值方法 . 这种方法需要更多的计算，但结果更可靠 . MATLAB具有 griddata 函数，给定散乱点首先生成点的三角剖分，然后在三角形顶部创建规则网格并插入网格点的值 . 因此，如果要插入点 [ri ti] 的值，则应应用第二个插值以从插值网格中获取点的值 .

  在spatialanalysisonline和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)]);
  

  More refinement: 由于已知搜索点被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)))
  

  Result:

  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)
  

  Conclusion ：

  与原始数据结构和采样方法相关的插值结果 . 如果采样方法匹配原始数据的模式，插值结果更准确，那么在极坐标网格点遵循数据模式的情况下，常规极坐标插值结果可以更可靠 . 但是如果常规极坐标不匹配数据结构或数据结构如不规则地形，基于三角测量的插值方法可以更好地表示数据 .

  回复于 2024-04-27T09:39:32+08:00
• 0

  我什么都不知道你试过了，但是 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
  

  回复于 2024-04-27T09:39:32+08:00
• 0

  请检查这个例子，我使用了两个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

  回复于 2024-04-27T09:39:32+08:00