Разбросанные данные состоят из набора точек X и соответствующие значения V, где точки не имеют структуры или порядка между их относительными расположениями. Существуют различные подходы к интерполяции разрозненных данных. В одном широко используемом подходе используется триангуляция точек Делоне.
В этом примере показано, как построить интерполированную поверхность путем триангуляции точек и подъема вершин на величину V в размер, ортогональный X.
Существуют различные варианты применения этого подхода. В этом примере интерполяция разбивается на отдельные этапы; обычно общий процесс интерполяции выполняется с помощью одного вызова функции.
Создайте рассеянный набор данных на поверхности параболоида.
X = [-1.5 3.2; 1.8 3.3; -3.7 1.5; -1.5 1.3; ... 0.8 1.2; 3.3 1.5; -4.0 -1.0; -2.3 -0.7; 0 -0.5; 2.0 -1.5; 3.7 -0.8; -3.5 -2.9; ... -0.9 -3.9; 2.0 -3.5; 3.5 -2.25]; V = X(:,1).^2 + X(:,2).^2; hold on plot3(X(:,1),X(:,2),zeros(15,1), '*r') axis([-4, 4, -4, 4, 0, 25]); grid stem3(X(:,1),X(:,2),V,'^','fill') hold off view(322.5, 30);
Создание триангуляции Делоне, подъем вершин и вычисление интерполятора в точке запроса Xq.
figure('Color', 'white') t = delaunay(X(:,1),X(:,2)); hold on trimesh(t,X(:,1),X(:,2), zeros(15,1), ... 'EdgeColor','r', 'FaceColor','none') defaultFaceColor = [0.6875 0.8750 0.8984]; trisurf(t,X(:,1),X(:,2), V, 'FaceColor', ... defaultFaceColor, 'FaceAlpha',0.9); plot3(X(:,1),X(:,2),zeros(15,1), '*r') axis([-4, 4, -4, 4, 0, 25]); grid plot3(-2.6,-2.6,0,'*b','LineWidth', 1.6) plot3([-2.6 -2.6]',[-2.6 -2.6]',[0 13.52]','-b','LineWidth',1.6) hold off view(322.5, 30); text(-2.0, -2.6, 'Xq', 'FontWeight', 'bold', ... 'HorizontalAlignment','center', 'BackgroundColor', 'none');
Этот шаг обычно включает в себя обход структуры данных триангуляции, чтобы найти треугольник, который заключает в себе точку запроса. После нахождения точки последующие шаги по вычислению значения зависят от метода интерполяции. Можно вычислить ближайшую точку в окрестности и использовать значение в этой точке (метод интерполяции ближайшего соседа). Также можно вычислить взвешенную сумму значений трех вершин окружающего треугольника (метод линейной интерполяции). Эти методы и их варианты рассматриваются в текстах и ссылках на интерполяцию разрозненных данных.
Хотя иллюстрация подсвечивает 2-D интерполяцию, этот метод можно применить к более высоким размерам. В более общих чертах, учитывая набор баллов X и соответствующие значения V, можно создать интерполятор формы V = F(X). Интерполант можно вычислить в точке запроса Xq, дать Vq = F(Xq). Это однозначная функция; для любой точки запроса Xq внутри выпуклого корпуса X, он будет производить уникальное значение Vq. Предполагается, что данные выборки соответствуют этому свойству для получения удовлетворительной интерполяции.
MATLAB ® предоставляет два способа выполнения интерполяции рассеянных данных на основе триангуляции:
scatteredInterpolant класс
griddata функция поддерживает интерполяцию 2-D разрозненных данных. griddatan функция поддерживает интерполяцию разрозненных данных в N-D; однако это нецелесообразно в измерениях выше 6-D для средних или больших наборов точек из-за экспоненциального роста памяти, необходимого для лежащей в основе триангуляции.
scatteredInterpolant класс поддерживает интерполяцию разрозненных данных в пространстве 2-D и 3-D. Использование этого класса поощряется, поскольку он является более эффективным и легко адаптируется к более широкому кругу проблем интерполяции.
griddata и griddatan функции берут набор точек выборки, X, соответствующие значения, Vи точки запроса, Xqи возвращает интерполированные значения, Vq. Синтаксис вызова аналогичен для каждой функции; основным отличием является 2-D/ 3-D griddata функция позволяет определить точки в терминах X, Y / X, Y, Z координаты. Эти две функции интерполируют разрозненные данные в предопределенных местоположениях точек сетки; цель состоит в создании данных с привязкой к сетке, отсюда и название. Интерполяция на практике носит более общий характер. Может потребоваться выполнить запрос в произвольных местах в выпуклом корпусе точек.
В этом примере показано, как griddata функция интерполирует разрозненные данные в наборе точек сетки и использует эти данные сетки для создания контурного графика.
Постройте график seamount набор данных (подводная гора - подводная гора). Набор данных состоит из набора долготы (x) и широты (y) местоположения и соответствующие seamount отметки (z), измеренный в этих координатах.
load seamount plot3(x,y,z,'.','markersize',12) xlabel('Longitude') ylabel('Latitude') zlabel('Depth in Feet') grid on
Использовать meshgrid для создания набора 2-D точек сетки в плоскости долготы-широты и последующего использования griddata для интерполяции соответствующей глубины в этих точках.
figure [xi,yi] = meshgrid(210.8:0.01:211.8, -48.5:0.01:-47.9); zi = griddata(x,y,z,xi,yi); surf(xi,yi,zi); xlabel('Longitude') ylabel('Latitude') zlabel('Depth in Feet')
Теперь, когда данные находятся в формате с сеткой, вычислите и постройте контуры.
figure [c,h] = contour(xi,yi,zi); clabel(c,h); xlabel('Longitude') ylabel('Latitude')
Также можно использовать griddata для интерполяции в произвольных местах в выпуклом корпусе набора данных. Например, глубина в координатах (211.3, -48.2) задается следующим образом:
zi = griddata(x,y,z, 211.3, -48.2);
Базовая триангуляция вычисляется каждый раз, когда griddata вызывается функция. Это может повлиять на производительность, если один и тот же набор данных многократно интерполируется с различными точками запроса. scatteredInterpolant класс, описанный в разделе Интерполяция рассеянных данных с использованием класса интерполяции рассеянных данных, является более эффективным в этом отношении.
Программное обеспечение MATLAB также обеспечивает griddatan для поддержки интерполяции в более высоких размерах. Синтаксис вызова аналогичен griddata.
griddata функция полезна при необходимости интерполяции для поиска значений в наборе предопределенных расположений точек сетки. На практике проблемы интерполяции часто являются более общими, и scatteredInterpolant класс обеспечивает большую гибкость. Класс имеет следующие преимущества:
Он создает интерполяционную функцию, которую можно эффективно запросить. То есть лежащая в основе триангуляция создается один раз и повторно используется для последующих запросов.
Метод интерполяции может быть изменен независимо от триангуляции.
Значения в точках данных могут быть изменены независимо от триангуляции.
Точки данных могут быть постепенно добавлены к существующему интерполятору без запуска полного повторного вычисления. Точки данных также можно эффективно удалять и перемещать при условии, что количество редактируемых точек мало по сравнению с общим количеством точек образца.
Она обеспечивает функции экстраполяции для аппроксимации значений в точках, которые выходят за пределы выпуклого корпуса. Дополнительные сведения см. в разделе Экстраполяция разбросанных данных.
scatteredInterpolant обеспечивает следующие методы интерполяции:
'nearest' - интерполяция ближайшего соседа, где интерполирующая поверхность является прерывистой.
'linear' - линейная интерполяция (по умолчанию), где интерполирующая поверхность C0 является непрерывной.
'natural' - интерполяция по естественному соседству, где интерполирующая поверхность является C1 непрерывной, за исключением точек выборки.
scatteredInterpolant класс поддерживает интерполяцию разрозненных данных в пространстве 2-D и 3-D. Интерполант можно создать путем вызова scatteredInterpolant и пропускают местоположения точек и соответствующие значения, и необязательно способы интерполяции и экстраполяции. См. раздел scatteredInterpolant для получения дополнительной информации о синтаксисах, которые можно использовать для создания и оценки scatteredInterpolant.
В этом примере показано, как использовать scatteredInterpolant для интерполяции рассеянной выборки peaks функция.
Создайте набор разрозненных данных.
rng default;
X = -3 + 6.*rand([250 2]);
V = peaks(X(:,1),X(:,2));Создайте интерполятор.
F = scatteredInterpolant(X,V)
F =
scatteredInterpolant with properties:
Points: [250x2 double]
Values: [250x1 double]
Method: 'linear'
ExtrapolationMethod: 'linear'
Points представляет координаты точек данных, и Values представляет связанные значения. Method свойство представляет метод интерполяции, выполняющий интерполяцию. ExtrapolationMethod свойство представляет метод экстраполяции, используемый, когда точки запроса выходят за пределы выпуклого корпуса.
Вы можете получить доступ к свойствам F таким же образом осуществляется доступ к полям struct. Например, использовать F.Points для проверки координат точек данных.
Вычислите интерполятор.
scatteredInterpolant обеспечивает подстрочную оценку интерполятора. Его оценивают так же, как функцию. Можно вычислить интерполятор следующим образом. В этом случае значение в местоположении запроса задается Vq. Можно вычислить в одной точке запроса:
Vq = F([1.5 1.25])
Vq = 1.4838
Можно также передать отдельные координаты:
Vq = F(1.5, 1.25)
Vq = 1.4838
Можно вычислить в векторе расположения точек:
Xq = [0.5 0.25; 0.75 0.35; 1.25 0.85]; Vq = F(Xq)
Vq = 3×1
0.4057
1.2199
2.1639
Вы можете оценить F в местах расположения точек сетки и распечатать результат.
[Xq,Yq] = meshgrid(-2.5:0.125:2.5); Vq = F(Xq,Yq); surf(Xq,Yq,Vq); xlabel('X','fontweight','b'), ylabel('Y','fontweight','b'); zlabel('Value - V','fontweight','b'); title('Linear Interpolation Method','fontweight','b');
Измените метод интерполяции.
Метод интерполяции можно изменить на лету. Задайте для метода значение 'nearest'.
F.Method = 'nearest';Выполните переоценку и постройте график интерполятора, как и раньше.
Vq = F(Xq,Yq); figure surf(Xq,Yq,Vq); xlabel('X','fontweight','b'),ylabel('Y','fontweight','b') zlabel('Value - V','fontweight','b') title('Nearest neighbor Interpolation Method','fontweight','b');
Измените метод интерполяции на естественный сосед, выполните переоценку и постройте график результатов.
F.Method = 'natural'; Vq = F(Xq,Yq); figure surf(Xq,Yq,Vq); xlabel('X','fontweight','b'),ylabel('Y','fontweight','b') zlabel('Value - V','fontweight','b') title('Natural neighbor Interpolation Method','fontweight','b');
Замените значения в местоположениях данных образца.
Можно изменить значения V в местах выборки данных, Xна лету. Это полезно на практике, поскольку некоторые проблемы интерполяции могут иметь несколько наборов значений в одних и тех же местоположениях. Например, предположим, что требуется интерполяция поля скорости 3-D, которое определяется расположениями (x, y, z) и соответствующих компонентных векторов скорости (Vx, Vy, Vz). Можно интерполировать каждый компонент скорости, назначив его свойству values (V) в свою очередь. Это имеет важные преимущества с точки зрения производительности, поскольку позволяет повторно использовать один и тот же интерполятор, не влекя за собой накладные расходы при каждом вычислении нового.
Следующие шаги показывают, как изменить значения в нашем примере. Значения вычисляются с помощью выражения xe-x2-y2.
V = X(:,1).*exp(-X(:,1).^2-X(:,2).^2); F.Values = V;
Вычислите интерполятор и постройте график результата.
Vq = F(Xq,Yq); figure surf(Xq,Yq,Vq); xlabel('X','fontweight','b'), ylabel('Y','fontweight','b') zlabel('Value - V','fontweight','b') title('Natural neighbor interpolation of v = x.*exp(-x.^2-y.^2)')
Добавьте дополнительные расположения точек и значения к существующему интерполятору.
Это выполняет эффективное обновление в отличие от полного повторного вычисления с использованием набора дополненных данных.
При добавлении данных образца важно добавить как местоположения точек, так и соответствующие значения.
Продолжая пример, создайте новые выборочные точки следующим образом:
X = -1.5 + 3.*rand(100,2); V = X(:,1).*exp(-X(:,1).^2-X(:,2).^2);
Добавьте новые точки и соответствующие значения в триангуляцию.
F.Points(end+(1:100),:) = X; F.Values(end+(1:100)) = V;
Вычислите уточненный интерполятор и постройте график результата.
Vq = F(Xq,Yq); figure surf(Xq,Yq,Vq); xlabel('X','fontweight','b'), ylabel('Y','fontweight','b'); zlabel('Value - V','fontweight','b');
Удаление данных из интерполятора.
Можно постепенно удалить точки данных выборки из интерполятора. Можно также удалить точки данных и соответствующие значения из интерполятора. Это полезно для удаления ложных отклонений.
При удалении данных образца важно удалить как местоположение точки, так и соответствующее значение.
Снимите 25-ю точку.
F.Points(25,:) = []; F.Values(25) = [];
Удалите точки с 5 по 15.
F.Points(5:15,:) = []; F.Values(5:15) = [];
Сохраните 150 очков и удалите остальные.
F.Points(150:end,:) = []; F.Values(150:end) = [];
Это создает более грубую поверхность при вычислении и построении графика:
Vq = F(Xq,Yq); figure surf(Xq,Yq,Vq); xlabel('X','fontweight','b'), ylabel('Y','fontweight','b'); zlabel('Value - V','fontweight','b'); title('Interpolation of v = x.*exp(-x.^2-y.^2) with sample points removed')
В этом примере показано, как интерполировать разрозненные данные, когда значение в каждом месте выборки является сложным.
Создайте образец данных.
rng('default')
X = -3 + 6*rand([250 2]);
V = complex(X(:,1).*X(:,2), X(:,1).^2 + X(:,2).^2);Создайте интерполятор.
F = scatteredInterpolant(X,V);
Создайте сетку точек запроса и вычислите интерполятор в точках сетки.
[Xq,Yq] = meshgrid(-2.5:0.125:2.5); Vq = F(Xq,Yq);
Постройте график реального компонента Vq.
VqReal = real(Vq); figure surf(Xq,Yq,VqReal); xlabel('X'); ylabel('Y'); zlabel('Real Value - V'); title('Real Component of Interpolated Value');
Постройте график мнимой составляющей Vq.
VqImag = imag(Vq); figure surf(Xq,Yq,VqImag); xlabel('X'); ylabel('Y'); zlabel('Imaginary Value - V'); title('Imaginary Component of Interpolated Value');
Многие из иллюстративных примеров в предыдущих разделах касались интерполяции наборов точек, которые были отобраны на гладких поверхностях. Кроме того, точки были относительно равномерно разнесены. Например, кластеры точек не были разделены относительно большими расстояниями. Кроме того, интерполятор был оценен хорошо внутри выпуклого корпуса в местах расположения точек.
При решении реальных проблем интерполяции данные могут быть более сложными. Это может происходить от измерительного оборудования, которое может производить неточные показания или отклонения. Базовые данные не могут изменяться плавно, значения могут резко переходить из точки в точку. В этом разделе приведены рекомендации по выявлению и устранению проблем с интерполяцией разрозненных данных.
Необходимо предварительно обработать образцы данных, которые содержат NaN значения для удаления NaN значения, поскольку эти данные не могут внести вклад в интерполяцию. Если NaN удаляется, соответствующие значения данных/координаты также должны быть удалены для обеспечения непротиворечивости. Если NaN значения присутствуют в образце данных, при вызове конструктора возникнет ошибка.
Следующий пример иллюстрирует удаление NaNs.
Создать некоторые данные и заменить некоторые записи на NaN:
x = rand(25,1)*4-2; y = rand(25,1)*4-2; V = x.^2 + y.^2; x(5) = NaN; x(10) = NaN; y(12) = NaN; V(14) = NaN;
F = scatteredInterpolant(x,y,V);
NaNs и затем построить интерполятор:nan_flags = isnan(x) | isnan(y) | isnan(V); x(nan_flags) = []; y(nan_flags) = []; V(nan_flags) = []; F = scatteredInterpolant(x,y,V);
NaN значения.X = rand(25,2)*4-2; V = X(:,1).^2 + X(:,2).^2; X(5,1) = NaN; X(10,1) = NaN; X(12,2) = NaN; V(14) = NaN;
F = scatteredInterpolant(X,V);
nan_flags = isnan(X(:,1)) | isnan(X(:,2)) | isnan(V); X(nan_flags,:) = []; V(nan_flags) = []; F = scatteredInterpolant(X,V);
griddata и griddatan вернуть NaN значения при запросе точек за пределами выпуклого корпуса с помощью 'linear' или 'natural' методы. Однако можно ожидать числовых результатов при запросе тех же точек с помощью 'nearest' способ. Это происходит потому, что ближайший сосед к точке запроса существует как внутри, так и снаружи выпуклого корпуса.
Если требуется вычислить приблизительные значения за пределами выпуклого корпуса, следует использовать scatteredInterpolant. Дополнительные сведения см. в разделе Экстраполяция разбросанных данных.
Входные данные редко являются «совершенными», и приложение может обрабатывать повторяющиеся местоположения точек данных. Две или более точек данных в одном и том же местоположении в наборе данных могут иметь различные соответствующие значения. В этом сценарии scatteredInterpolant объединяет точки и вычисляет среднее из соответствующих значений. В этом примере показано, как scatteredInterpolant выполняет интерполяцию набора данных с повторяющимися точками.
Создайте образец данных, лежащий на плоской поверхности:
x = rand(100,1)*6-3; y = rand(100,1)*6-3; V = x + y;
Введите расположение повторяющейся точки, назначив координаты точки 50 точке 100:
x(50) = x(100); y(50) = y(100);
Создайте интерполятор. Обратите внимание, что F содержит 99 уникальных точек данных:
F = scatteredInterpolant(x,y,V)
Проверьте значение, связанное с 50-й точкой:
F.Values(50)
Это значение является средним значением исходного 50-го и 100-го значений, поскольку эти две точки данных имеют одинаковое местоположение:
(V(50)+V(100))/2
В некоторых проблемах интерполяции множественные наборы значений выборки могут соответствовать одним и тем же местоположениям. Например, набор значений может быть записан в одних и тех же местоположениях в разные периоды времени. Для эффективности можно интерполировать один набор показаний, а затем заменить значения для интерполяции следующего набора.
Всегда используйте согласованное управление данными при замене значений при наличии дубликатов точек. Предположим, что у вас есть два набора значений, связанных со 100 расположениями точек данных, и вы хотите по очереди интерполировать каждый набор путем замены значений.
Рассмотрим два набора значений:
V1 = x + 4*y; V2 = 3*x + 5*y
Создайте интерполятор. scatteredInterpolant объединяет повторяющиеся местоположения, и интерполятор содержит 99 уникальных точек выборки:
F = scatteredInterpolant(x,y,V1)
F.Values = V2 означает присвоение 100 значений 99 образцам. Контекст предыдущей операции слияния потерян; количество местоположений выборки не будет соответствовать количеству значений выборки. Интерполятор потребует устранения несогласованности для поддержки запросов.В этом более сложном сценарии перед созданием и редактированием интерполятора необходимо удалить дубликаты. Используйте unique для поиска индексов уникальных точек. unique может также выводить аргументы, идентифицирующие индексы повторяющихся точек.
[~, I, ~] = unique([x y],'first','rows'); I = sort(I); x = x(I); y = y(I); V1 = V1(I); V2 = V2(I); F = scatteredInterpolant(x,y,V1)
F для интерполяции первого набора данных. Затем можно заменить значения для интерполяции второго набора данных.F.Values = V2;
scatteredInterpolant позволяет редактировать свойства, представляющие значения образца (F.Values) и метод интерполяции (F.Method). Поскольку эти свойства не зависят от лежащей в основе триангуляции, правки могут быть выполнены эффективно. Однако, как и при работе с большим массивом, следует позаботиться о том, чтобы случайно не создавать ненужные копии при редактировании данных. Копии создаются, когда несколько переменных ссылаются на массив, и этот массив затем редактируется.
Копия не сделана в следующем:
A1 = magic(4) A1(4,4) = 11
A1 и A2 не может больше совместно использовать те же данные после выполнения редактирования:A1 = magic(4) A2 = A1 A1(4,4) = 32
scatteredInterpolant содержит данные и ведет себя как массив - на языке MATLAB он называется объектом значения. Язык MATLAB обеспечивает оптимальную производительность, когда приложение структурировано в функции, находящиеся в файлах. Создание прототипов в командной строке может не обеспечить такой же уровень производительности.Следующий пример демонстрирует это поведение, но следует отметить, что прирост производительности в этом примере не обобщается на другие функции в MATLAB.
Этот код не обеспечивает оптимальной производительности:
x = rand(1000000,1)*4-2; y = rand(1000000,1)*4-2; z = x.*exp(-x.^2-y.^2); tic; F = scatteredInterpolant(x,y,z); toc tic; F.Values = 2*z; toc
Когда MATLAB выполняет программу, состоящую из функций, которые находятся в файлах, он имеет полную картину выполнения кода; это позволяет MATLAB оптимизировать производительность. При вводе кода в командной строке MATLAB не может предвидеть, что вы собираетесь ввести дальше, поэтому он не может выполнять тот же уровень оптимизации. Разработка приложений путем создания многоразовых функций является общей и рекомендуемой практикой, и MATLAB оптимизирует производительность в этой настройке.
Триангуляция Делоне хорошо подходит для задач интерполяции разрозненных данных, поскольку она имеет благоприятные геометрические свойства, которые дают хорошие результаты. Эти свойства:
Отказ от щелевидных треугольников/тетраэдров в пользу более равносторонних.
Свойство пустой окружности, неявно определяющее отношение ближайшего соседа между точками.
Свойство пустой окружности гарантирует, что на интерполированные значения влияют точки выборки в окрестности расположения запроса. Несмотря на эти качества, в некоторых ситуациях распределение точек данных может привести к плохим результатам, и это обычно происходит вблизи выпуклого корпуса набора данных выборки. Когда интерполяция дает неожиданные результаты, график данных выборки и лежащая в его основе триангуляция часто может дать представление о проблеме.
В этом примере показана интерполированная поверхность, которая ухудшается вблизи границы.
Создайте образец набора данных, в котором возникнут проблемы вблизи границы.
t = 0.4*pi:0.02:0.6*pi; x1 = cos(t)'; y1 = sin(t)'-1.02; x2 = x1; y2 = y1*(-1); x3 = linspace(-0.3,0.3,16)'; y3 = zeros(16,1); x = [x1;x2;x3]; y = [y1;y2;y3];
Теперь поднимите эти точки образца на поверхность + y2 и интерполируйте поверхность.
z = x.^2 + y.^2; F = scatteredInterpolant(x,y,z); [xi,yi] = meshgrid(-0.3:.02:0.3, -0.0688:0.01:0.0688); zi = F(xi,yi); mesh(xi,yi,zi) xlabel('X','fontweight','b'), ylabel('Y','fontweight','b') zlabel('Value - V','fontweight','b') title('Interpolated Surface');
Фактическая поверхность:
zi = xi.^2 + yi.^2;
figure
mesh(xi,yi,zi)
title('Actual Surface')
Чтобы понять, почему интерполяционная поверхность ухудшается вблизи границы, полезно посмотреть на лежащую в основе триангуляцию:
dt = delaunayTriangulation(x,y); figure plot(x,y,'*r') axis equal hold on triplot(dt) plot(x1,y1,'-r') plot(x2,y2,'-r') title('Triangulation Used to Create the Interpolant') hold off
Треугольники в красных границах имеют относительно хорошую форму; они построены из точек, которые находятся в непосредственной близости, и интерполяция хорошо работает в этом регионе. За пределами красной границы треугольники являются щелевидными и соединяют точки, удаленные друг от друга. Нет достаточного отбора проб для точного захвата поверхности, поэтому неудивительно, что результаты в этих регионах плохи. В 3-D визуальный осмотр триангуляции становится немного сложнее, но взгляд на распределение точек часто может помочь проиллюстрировать потенциальные проблемы.
MATLAB ® 4griddata способ, 'v4', не основан на триангуляции и не подвержен ухудшению поверхности интерполяции вблизи границы.
[xi,yi] = meshgrid(-0.3:.02:0.3, -0.0688:0.01:0.0688); zi = griddata(x,y,z,xi,yi,'v4'); mesh(xi,yi,zi) xlabel('X','fontweight','b'), ylabel('Y','fontweight','b') zlabel('Value - V','fontweight','b') title('Interpolated surface from griddata with v4 method','fontweight','b');
Интерполированная поверхность из griddata с использованием 'v4' способ соответствует ожидаемой фактической поверхности.