Как интерполировать, используя в полярных координатах

Question

У меня есть полярные координаты, радиус 0.05 <= r <= 1 а также 0 ≤ θ ≤ 2π, Радиус r составляет 50 значений от 0,05 до 1, а полярный угол θ составляет 24 значения от 0 до 2π.

Как я могу интерполировать 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];

и вы хотите, чтобы взаимосвязь контрольной точки.

Когда вы выбираете некоторые данные, чтобы использовать их для интерполяции, вы должны подумать, как сделать выборку данных в соответствии с вашими требованиями. Поэтому, когда вы выбираете регулярную сетку полярных координат, эти координаты при преобразовании в прямоугольные образуют круглую форму, в которой большинство точек сосредоточено в центре крикла, а когда мы перемещаемся от центра к внешним областям, расстояние между точками выросла.

%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, чтобы увеличить точность.
С другой стороны, если вы создаете регулярную сетку в прямоугольных координатах, равное значение придается каждому региону. Так что точность интерполяции будет одинаковой во всех регионах.
Для этого сначала вы создаете регулярную сетку в полярных координатах, а затем конвертируете сетку в прямоугольные координаты, чтобы вы могли вычислить экстенты (min max) точек выборки в прямоугольных координатах. На основе этих экстентов вы можете создать регулярную сетку в прямоугольных координатах. Регулярную сетку прямоугольных координат затем преобразовать в полярные координаты, чтобы получить z для точек сетки по формуле z2.

%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, вы можете сделать следующее:

1. Интерполировать полярные координаты с interp2 (интерполяция для данных с привязкой к сетке) этот подход прост, но имеет тот недостаток, что r и тета не имеют одинакового масштаба, и это может повлиять на точность интерполяции.

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

2. преобразовать полярные координаты в прямоугольные, а затем применить метод интерполяции, предназначенный для рассеянных данных. этот подход требует большего количества вычислений, но результат этого более надежен. MATLAB имеет griddata Функция, которая задает рассеянные точки, сначала генерирует триангуляцию точек, а затем создает регулярную сетку поверх треугольников и интерполирует значения точек сетки. Так что если вы хотите интерполировать значение точки [ri ti] Затем вы должны применить вторую интерполяцию, чтобы получить значение точки из интерполированной сетки.

С помощью некоторой информации из пространственного анализа online и 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)

Вывод:

Результат интерполяции связан со структурой исходных данных и методом выборки. Если метод выборки соответствует шаблону исходных данных, результат интерполяции является более точным, поэтому в случаях, когда точки сетки полярных координат следуют шаблону данных, результат интерполяции регулярной полярной координаты может быть более надежным. Но если регулярные полярные координаты не соответствуют структуре данных или структура данных, например, нерегулярный рельеф, метод интерполяции, основанный на триангуляции, может лучше представлять данные.

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

Пожалуйста, проверьте этот пример, я использовал два для цикла, внутри для цикла я использовал оператор условия, если вы прокомментируете этот оператор условия и запустите программу, вы получите правильный ответ, после того, как вы раскомментируете этот оператор условия и запустите программу, вы я получу неправильный ответ. пожалуйста, проверь это.

% 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