MATLAB.Exponenta
–Û·Ë͇ Matlab&Toolboxes

Spline Toolbox

Обзор средств MATLAB и ToolBox'ов для приближения данных

Работа с разбросанными данными

В MATLAB входят функции для работы с разбросанными данными, т.е. например в двумерном случае это данные zi, заданные в точках (xi, yi)i=1,2,...,k. Эти функции позволяют решать следующие задачи.

  • Интерполировать разбросанные данные различными способами. Для данных, определенных в точках на плоскости, возможно кусочно-постоянное приближение по ближайшему соседу, кусочно-линейное приближение и более гладкое приближение с непрерывными производными (функция griddata). Данные, определенные в точках трехмерного и n-мерного пространства могут быть приближены по ближайшему соседу и при помощи кусочно-линейных функций (функции griddata3 и griddatan). Значение интерполянта возвращается на ортогональной сетке, накрывающей точки, в которых заданы разбросанные данные.
  • Строить выпуклую оболочку множества точек на плоскости, в трехмерном и n-мерном пространстве (функции convhull, convhulln).
  • Строить триангуляцию Делоне для множества точек на плоскости, в трехмерном и n-мерном пространстве ( функции delaunay, delaunay3, delaunayn).
  • Строить диаграмму Вороного для множества точек на плоскости, в трехмерном и n-мерном пространстве (функции voronoi, voronoin).
  • Визуализировать выпуклую оболочку, триангуляцию Делоне и диаграмму Вороного (функции triplot, trimesh, trisurf, tetramesh).
  • Для заданной точки определять ближайшую точку к ней из заданного множества точек на плоскости (функция dsearch).
  • Для заданной точки определить содержащий ее треугольник из триангуляции Делоне (функция tsearch), или сиплекс в n-мерном случае (функция tsearchn).
  • Решать такие задачи, как:
    • нахождение площади многоугольника (функция polyarea);
    • площади пересечения четырехугольников (функция rectint).

Нахождение выпуклой оболочки множества точек и ее визуализация

Выпуклая оболочка conv(P) множества точек, принадлежащих n-мерному пространству Rn, определяется как наименьшее выпуклое множество, содержащее все эти точки. Выпуклое множество, порожденное точками (pi)i=1,2,...,k, определяется как множество точек:

p = a1p1 + a2p2 + ... + akpk

где все ai >/= 0 и a1 + a2 + ... + ak = 1 (здесь точки отождествляются с векторами, выходящими из начала координат).

В двумерном случае есть простая интуитивно понятная аналогия. Пусть на плоскости лежат точки (вбиты гвозди), мы натягиваем вокруг них веревку.

То, что получится, и будет выпуклой оболочкой (на правом рисунке ее граница изображена красным цветом). Все точки из множества находятся либо внутри выпуклой оболочки, либо на ее границе.

Выпуклая оболочка является пересечением конечного числа полупространств (полупространство - часть пространства, лежащая по одну сторону от гиперплоскости). В двумерном случае гиперплоскости это линии, а в трехмерном случае - плоскости.

Для построения выпуклой оболочки в MATLAB служит функция convhull. Ее входными аргументами являются вектора со координатами абсцисс и ординат точек, а выходными аргументами являются номера точек, образующих вершины выпуклой оболочки. В следующем примере генерируется набор из 100 точек с координатами, изменяющимися случайным образом от 0 до 1, затем при помощи функции convhull строится выпуклая оболочка. Номера точек, являющихся вершинами выпуклой оболочки, заносятся в массив kvert. Затем вершины выпуклой оболочки выделяются круглыми маркерами и рисуется сама выпуклая оболочка

x=rand(1,100);
y=rand(1,100);
figure
plot(x,y,'k.')
hold on
axis([-0.5 1.5 -0.5 1.5])
kvert=convhull(x,y)
pause(0.3)
plot(x(kvert(1)),y(kvert(1)),'ro')
for k=2:length(kvert)
    plot([x(kvert(k-1)) x(kvert(k))],[y(kvert(k-1)) y(kvert(k))],'r')
    plot(x(kvert(k)),y(kvert(k)),'ro')
    pause(0.3)
end

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

>> [kvert, s] = convhull(x, y);
>> s
s =
    0.9164 

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

В качестве примера приведем построение выпуклой оболочки точек, которые случайным образом распределены в кубе [0,1]x[0,1]x[0,1]. Координаты точек генерируются при помощи функции rand и записываются в вектор-столбцы. Из этих столбцов создается матрица, которая указывается в качестве входного аргумента функции convhulln. Функция convhulln возвращает матрицу из трех столбцов, каждый из которых содержит номера точек исходного множества, образующих некоторую грань выпуклой оболочки. Для визуализации выпуклой оболочки применяется функция trimesh, которая строит полигональный объект (patch) и возвращает указатель на него. Вид этого полигонального объекта можно изменять по своему усмотрению, включая прозрачность и цвет граней, маркеры в вершинах и т.д. Когда выпуклая оболочка найдена, то все исходное множество точек разбивается на два подмножества - точки, которые являются вершинами выпуклой оболочки, и внутренние точки. В приводимом ниже примере рядом с точками этих типов размещаются подписи разного цвета с номерами точек. Для этого используется функция text. При этом возникает вопрос, как определить номера точек, попавших внутрь выпуклой оболочки исходного множества, и тех, которые принадлежат вершинам выпуклой оболочки. Как уже говорилось, номера вершин записаны в выходном аргументе функции convhulln (двумерном массиве с тремя столбцами). Однако, они встречаются в разных строках этого массива и не по одному разу. Поэтому для выделения номеров точек, попавших в вершины выпуклой оболочки, можно поступить так: последовательно объединить содержимое столбцов матрицы с номерами вершин, воспользовавшись функцией union. Поскольку данная функция производит объединение в смысле объединения множеств, то повторы элементов исключаются. Аналогично, для нахождения номеров точек, попавших внутрь выпуклой оболочки, можно применить функцию setdiff, которая находит разность множеств. В нашем случае из множества номеров всех точек следует вычесть номера точек, образующих вершины выпуклой оболочки.

k=20; задаем двадцать точек
rand('state',0); генератору случайных чисел задаем начальное состояние
% генерируем координаты точек множества и записываем их в вектор-столбцы x, y, и z
x=rand(k,1);
y=rand(k,1);
z=rand(k,1);
% образуем двумерный массив с координатами точек множества 
X=[x y z];
% находим номера точек, образующих треугольные грани выпуклой оболочки множества
kvert=convhulln(X)
% визуализируем выпуклую оболочку при помощи функции trisurf
figure('Color','w')
h=trisurf(kvert,X(:,1),X(:,2),X(:,3))
hold on
% устанавливаем прозрачность граней, красные ребра 
% и круглые красные маркеры в вершинах выпуклой оболочки
set(h,'FaceAlpha',0.1,'FaceColor','g',...
    'EdgeColor','r','Marker','o','MarkerSize',10)
% находим номера точек, являющихся вершинами выпуклой оболочки, и записываем их в kv
kv=union(kvert(:,1),kvert(:,2))
kv=union(kv,kvert(:,3))
% подписываем номера вершин выпуклой оболочки
ht=text(x(kv), y(kv), z(kv),num2str(kv))
set(ht,'BackgroundColor','r','VerticalAlignment','Top',...
    'FontWeight','bold')
% находим номера точек, попавших внутрь оболочки, и записываем их в kinner
kinner=setdiff(1:k,kv)
% рисуем точек, попавшие внутрь оболочки, круглыми синими маркерами
plot3(x(kinner),y(kinner),z(kinner),'Color','b','LineStyle','none',...
    'Marker','o','MarkerSize',10)  
% помещаем рядом сточками, попавшими внутрь оболочки, их номера
ht=text(x(kinner), y(kinner), z(kinner),num2str(kinner'))
set(ht,'BackgroundColor','y','VerticalAlignment','Top',...
    'FontWeight','bold','Color','b')

В результате выполнения этих команд на оси графического окна выводится выпуклая оболочка, ее пронумерованные вершины и пронумерованные точки, попавшие внутрь выпуклой оболочки, в соответствии с их номерами в исходном множестве:

Вместо функции trisurf можно использовать функцию trimesh, которая в отличие от trimesh не заливает грани цветами.

При вызове функции convhulln(X) со вторым дополнительным выходным аргументом, в нем возвращается объем выпуклой оболочки

>> [kvert,V]=convhulln(X);
>> V
V =
    0.2582

В MATLAB для построения выпуклой оболочки используется функция qconvex пакета Qhull (http://www.qhull.org). Эти функция поддерживает ряд опций для задания некоторых параметров алгоритма построения выпуклой оболочки, которые можно задавать и в функциях convhull, convhulln.

Алгоритмы построения выпуклой оболочки, в том числе и алгоритм Quickhull, описаны, например, на http://algolist.manual.ru/maths/geom/convhull/. Там же приведено несколько ссылок на статьи с описаниями алгоритмов построения выпуклой оболочки множества точек.

Для визуализации выпуклой оболочки точек трехмерного пространства не обязательно использовать функцию trisurf. Вместо нее можно самому сконструировать полигональный объект с заданными способами заливки цветом ребер и граней и их прозрачностью.

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

  • Vertices - матрица из трех столбцов, в строки которой записаны координаты вершин полигонального объекта;
  • Faces - матрица, каждая строка которой содержит номера вершин, образующих различные грани полигонального объекта
  • FaceVertexCData - цвета вершин, в примере они задаются в формате RGB, цвет устанавливается зеленый, а его оттенок зависит от высоты точки над плоскостью xy.
  • FaceColor - способ заливки граней полигонального объекта цветом, в примере используется значение 'interp', т.е. плавное изменение цвета от вершины к вершине.

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

% Задание множества точек
k=30;
rand('state',0)
x=rand(k,1);
y=rand(k,1);
z=rand(k,1);
X=[x y z];
% Построение выпуклой оболочки
kvert=convhulln(X);
figure('Color','w')
% Задание цвета вершин
VColor=[zeros(size(z)) z zeros(size(z))];
% Создание полигонального объекта
patch('Vertices',X,'Faces',kvert,'FaceVertexCData',VColor,'FaceColor','interp')
% Установка трехмерного вида осей.
view(3)

Триангуляция Делоне

Если задан набор точек на плоскости, то задача триангуляции такого набора состоит в соединении всех точек непересекающимися отрезками так, чтобы новых отрезков уже нельзя было добавить без пересечения с имеющимися, например

В этом примере видно, что триангуляция может быть неединственная, например

также является триангуляцией.

Важный в приложениях класс триангуляций составляют триангуляции Делоне. Триангуляция называется триангуляцией Делоне, если внутрь окружности, описанной вокруг каждого треугольника, не попадают точки множества. В предыдущих примерах триангуляция, построенная синими линиями, является триангуляцией Делоне, т.к. она удовлетворяет этому определению (см. левый рисунок, на котором окружности нарисованы не для всех треугольников только для того, чтобы не загромождать рисунок). Триангуляция, построенная зелеными линиями, не удовлетворяет условию Делоне (см. правый рисунок).

Триангуляция Делоне может быть неединственной (в случае, когда четыре точки множества лежат на одной окружности), например обе триангуляции одного и того же набора точек, приведенные ниже, являются триангуляциями Делоне.

Можно показать, что триангуляция Делоне:

  • обладает максимальной суммой минимальных углов треугольников среди всех триангуляций множества точек;
  • минимальной суммой радиусов описанных вокруг всех треугольников окружностей по сравнению со всевозможными другими триангуляциями.

В Интернет о триангуляции Делоне можно почитать, например, на сайте http://algolist.manual.ru (в разделе геометрия), где приведены ссылки на несколько доступных для скачивания источников с описаниями алгоритмов триангуляции.

Для построения триангуляции Делоне множества точек на плоскости в MATLAB используется функция delaunay, входными аргументами которой являются вектора с координатами точек, а выходным аргументом - матрица из трех столбцов, каждая строка которой содержит номера точек, образующие вершины треугольников. Для визуализации триангуляции удобно использовать функцию triplot, входными аргументами которой являются матрица, возвращаемая функцией delaunay, и два вектора с координатами точек. Функция triplot возвращает указатель на создаваемый ею полигональный объект, который в дальнейшем можно изменять по своему усмотрению.

Построим триангуляцию Делоне следующего множества точек

x = [4.4 8.9 8.1 4.5 0.1];
y = [6.9 6.1 0.7 1.9 0.2];

Отображаем их маркерами в графическом окне

figure('Color','w')
plot(x,y,'Marker','.','MarkerEdgeColor','r','MarkerSize',35,'LineStyle','none')
axis equal
axis off
hold on

Вызываем функцию delaunay для построения триангуляции

T = delaunay(x,y)

Рисуем триангуляцию при помощи функции triplot и устанавливаем линиям толщину 2пт.

h=triplot(T,x,y)
set(h,'LineWidth',2)

Легко убедиться, что построенная триангуляция действительно подходит под определение триангуляции Делоне. Построим для каждого треугольника описанную вокруг него окружность. Для этого обойдем треугольники (строки матрицы T) в цикле и по известным формулам вычислим координаты центра и радиус описанной вокруг него окружности. Сами окружности нарисуем при помощи функции rectangle, которая позволяет скруглять углы для получения эллипсов и кругов:

[m,n]=size(T); % определяем размер матрицы T
% Проходим по треугольникам в цикле
for k=1:m
    % Записываем в p1, p2, p3 номера точек, являющихся вершинами текущего треугольника
    p1=T(k,1); 
    p2=T(k,2); 
    p3=T(k,3); 
    % Находим координаты центра (xc, yc) описанной окружности
    a=(y(p2)-y(p1))/(x(p2)-x(p1));
    b=(y(p3)-y(p2))/(x(p3)-x(p2));
    xc=(a*b*(y(p1)-y(p3))+b*(x(p1)+x(p2))-a*(x(p2)+x(p3)))*0.5/(b-a)
    yc=-1/a*(xc-(x(p1)+x(p2))*0.5)+(y(p1)+y(p2))*0.5;
    % Находим радиус
    r=sqrt((xc-x(p1))^2+(yc-y(p1))^2);
    % строим описанную окружность 
    rectangle('Curvature',[1 1],'Position',[xc-r, yc-r, 2*r, 2*r])
end

В результате убеждаемся, что триангуляция является триангуляцией Делоне.

Триангуляцию удобно использовать для визуализации функций, заданных на непрямоугольной области определения. Сначала функция вычисляется в точках этой области, затем строится триангуляция Делоне множества точек на плоскости. После этого применяется функция trimesh, либо функция trisurf для визуализации:

% генерируем точки на плоскости
rand('state',0);
r=rand(1,500);
t=rand(1,500);
x = 2*r.*sin(2*pi*t)-1;
y = 2*r.*cos(2*pi*t)-1;
% вычисляем значение функции
z = exp(-x.^2-y.^2).*sin(pi*x).*y;
% находим триангуляцию Делоне
T = delaunay(x,y)
% строим поверхность двумя способами
figure('Color','w')
subplot(1,2,1)
trimesh(T,x,y,z)
subplot(1,2,2)
trisurf(T,x,y,z)

В трехмерном случае триангуляция Делоне определяется аналогично, только строятся тетраэдры, а сферы, описанные вокруг каждого тетраэдра, не должны содержать внутри себя точек из триангулируемого множества. Для построения триангуляции точек трехмерного пространства используется функция delaunay3, во входных аргументах которой задаются векторы с координатами точек триангулируемого множества. Выходным аргументом функции delaunay3 является матрица из четырех столбцов, каждый из которых содержит номера вершин симплекса (тетраэдра) триангуляции. Для визуализации триангуляции в трехмерном пространстве удобно применять функцию tetramesh, первым входным аргументом которой является матрица, возвращаемая delaunay3, а вторым входным аргументом - матрица, каждая строка которой содержит координаты точек триангулируемого множества. Функция tetramesh возвращает указатель на создаваемый ею полигональный объект, свойства которого можно изменять.

% генерируем набор из девяти точек в трехмерном пространстве 
rand('state',0);
x = rand(9,1);
y = rand(9,1);
z = rand(9,1);
% Строим триангуляцию Делоне
T = delaunay3(x,y,z);
% Визуализируем тетраэдры
figure('Color','w')
h=tetramesh(T,[x y z])
% Задаем маркеры в вершинах и прозрачность граней
set(h,'Marker','.','MarkerEdgeColor','r','MarkerSize',30,'FaceAlpha',0.6)

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

Функции delaunay, delaunay3 и delaunayn используют функцию qdelaunay пакета Qhull (http://www.qhull.org), которая позволяют задавать ряд опций для управления процессом триангуляции.

Построение диаграммы Вороного

Если на плоскости задано множество различных точек (xi,yi)i=1,2,...,k, то диаграмма Вороного определяется как деление плоскости на k ячеек (по числу точек), каждая ячейка содержит одну точку из исходного набора и для любой ячейки выполняется следующее свойство. Ячейка, которая содержит точку (xi,yi), содержит также все точки плоскости, которые ближе к (xi,yi), чем к любой другой точке из исходного набора. Для примера приведем диаграмму Вороного, на которой точки множества заданы красным цветом, а сама диаграмма синими линиями

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

Каждая ячейка диаграммы, содержащая (xi,yi), называется многоугольником Вороного, соответствующим точке (xi,yi), ребра и вершины, соответственно, ребра диаграммы Вороного и вершины диаграммы Вороного.

Приведем некоторые свойства диаграммы Вороного, которые доказываются довольно просто в предположении, что никакие четыре точки исходного множества не лежат на одной окружности (см. например, классическую монографию, переведенную на русский язык: Ф.Препарата, М.Шеймос. Вычислительная геометрия: Введение. М. "Мир". 1989).

  1. Ребро диаграммы Вороного является отрезком прямой, перпендикулярной отрезку, соединяющему некоторые две точки из исходного набора, и проходящей ровно через середину этого отрезка. Каждое ребро принадлежит ровно двум многоугольникам.
  2. Каждый ближайший сосед некоторой точки из исходного множества определяет ребро в многоугольнике Вороного.
  3. Вершины диаграммы Вороного являются точками пересечения ровно трех ребер диаграммы Вороного.
  4. Если провести окружность с центром в вершине диаграммы Вороного, проходящую через три точки из исходного множества, которые принадлежат смежным с вершиной многоугольникам, то внутри нее не будут содержаться никакие точки из исходного множества.
  5. Очевидно, часть многоугольников диаграммы Вороного неограниченна. Многоугольник является неограниченным тогда и только тогда, когда соответствующая ему точка из исходного множества лежит на границе выпуклой оболочки всех точек множества.
  6. Диаграмма Вороного для множества из k точек содержит не более 2k-5 вершин и 3k-6 ребер.
  7. Между диаграммой Вороного и триангуляцией Делоне есть прямая связь (см. следующий раздел Как связаны выпуклая оболочка, триангуляция Делоне и диаграмма Вороного).

Для построения диаграммы Вороного заданного множества точек плоскости в MATLAB используется функция voronoi, которая допускает несколько вариантов вызова.

  • voronoi(x,y) - отображает в графическом окне множество точек, координаты которых заданы в массивах x и y, и диаграмму Вороного этого множества точек. Ребра, соответствующие бесконечным ячейкам, не строятся.
  • voronoi(x,y,T) - делает то же самое, что и voronoi(x,y), но не вычисляет предварительно триангуляцию Делоне, а использует триангуляцию, заданную в массиве T. Для получения массива T следует предварительно вызвать функцию delaunay так: T = delaunay(x,y) (см. раздел Триангуляция Делоне).
  • voronoi(hax,...) - работает аналогично предыдущим вариантам, только строит диаграмму Вороного на осях с указателем hax.
  • voronoi(...,'LineSpec') - позволяет задать цвет и стиль линий, используемых для построения диаграммы Вороного, а так же маркеры в вершинах. Например voronoi(x, y, 'ko:') рисует ребра черными пунктирными линиями и помещает круглые маркеры в вершинах диаграммы.
  • h = voronoi(...) - возвращает вектор указателей на линии, которыми прорисовываются ребра диаграммы.
  • [vx, vy] = voronoi(...) - возвращает координаты начальной и конечной точки каждого ребра (т.е. вершин диаграммы) в массивах vx и vy, которые являются матрицами из двух строк. Соответствующие координаты записаны в столбцах матриц, причем каждый столбец соответствует ребру диаграммы. При таком способе обращения к функции voronoi диаграмма Вороного не строится. Ее можно построить так: plot(vx,vy), а затем добавить исходные точки отобразив их некоторыми маркерами: hold on, plot(x,y,'.'). Для построения одновременно и диаграммы Вороного и множества точек все можно объединить в одну команду: plot(vx, vy, x, y, '.').
  • Для получения диаграммы в виде закрашенных многоугольников следует использовать функцию voronoin, которая предназначена для построения диаграммы Вороного для точек из трехмерного пространства и пространства более высоких размерностей (об этом чуть ниже).

Для построения диаграммы Вороного в функции voronoi предварительно строится триангуляция Делоне множества точек (см. раздел Триангуляция Делоне), для чего вызывается функция delaunay пакета MATLAB. Она в свою очередь, основана на функции qdelaunay пакета Qhull (http://www.qhull.org), которая реализует алгоритм Quickhull и позволяют задавать ряд опций для управления процессом триангуляции. Поэтому в функции voronoi, допустимо задание тех же самых опций, что и в delaunay. Они указываются в массиве ячеек строк последним входным аргументом в voronoi.

Приведем пример, в котором демонстрируется построение диаграммы Вороного для случайно сгенерированного набора из восьми точек и то обстоятельство, что точки, соответствующие бесконечным многоугольникам, лежат в вершинах выпуклой оболочки исходного множества точек. Точки рисуются красным цветом, грани диаграммы Вороного зеленым цветом, а грани выпуклой оболочки черным цветом. Для построения выпуклой оболочки применяется функция convhull (см. Нахождение выпуклой оболочки множества точек и ее визуализация).

 
% генерируем случайный набор из восьми точек
rand('state',0);
x = rand(8,1);
y = rand(8,1);
% строим диаграмму Вороного
figure('Color','w')
h=voronoi(x,y);
% устанавливаем зеленый цвет линий многоугольников и толщину 2пт
set(h,'Color','g','LineWidth',2)
hold on
% строим набор точек красными маркерами
plot(x,y,'Marker','.','MarkerEdgeColor','r','MarkerSize',30, 'LineStyle', 'none')
% ищем выпуклую оболочку заданного набора точек
k = convhull(x,y);
% рисуем выпуклую оболочку черными линиями
plot(x(k),y(k),'k')
% ищем пределы по осям x и y, которые содержат координаты заданных точек 
% и вершины диаграммы Вороного
X=get(h,'XData');
Y=get(h,'YData');
m=length(X)-1;
XX=cell2mat(X(1:m-1));
YY=cell2mat(Y(1:m-1));
xmin=min([XX(:);x]); xmax=max([XX(:);x]); 
ymin=min([YY(:);y]); ymax=max([YY(:);y]); 
% задаем пределы осей и делаем их невидимыми
axis([xmin xmax ymin ymax])
axis off

Для построения диаграммы Вороного в трехмерном пространстве и в пространствах большей размерности предназначена функция voronoin. Ее интерфейс оказывается полезным и при визуализации диаграммы Вороного для множества точек плоскости. Если координаты каждой точки множества записаны в строках массива X (т.е. число его столбцов совпадает с размерностью пространства), то

[V,C] = voronoin(X)

возвращает следующие массивы

  • Двумерный числовой массив V, содержащий координаты вершин диаграммы Вороного. Координаты каждой вершины записываются в строки массива V, причем первая строка всегда содержит Inf (т.е. бесконечность) и соответствует "бесконечно удаленной вершине", что нужно для отличия бесконечных ячеек диаграммы Вороного от конечных.
  • Вектор ячеек C, каждая ячейка которого содержит числовой вектор с номерами вершин, образующих ячейки диаграммы Вороного. Если один из элементов этого равняется единице (т.е. номеру бесконечно удаленной вершины, добавленной первой в массив V), то соответствующая ячейка диаграммы Вороного является неограниченной.

Используя массивы V и C, возвращаемые функцией voronoin, можно самостоятельно выбрать один из способов визуализации конечных ячеек диаграммы Вороного. Единственное, для каждой конечной ячейки придется построить выпуклую оболочку по ее вершинам (ячейки диаграммы Вороного являются выпуклыми политопами). Для этого можно использовать функцию convhulln (см. раздел Нахождение выпуклой оболочки множества точек и ее визуализация).

В следующем примере поясняется построение конечных ячеек диаграммы Вороного для точек трехмерного пространства.

% генерируем восемь случайных точек
k=8
rand('state',0)
x=rand(k,1);
y=rand(k,1);
z=rand(k,1);
% отображаем точки на графике
figure('Color','w')
plot3(x,y,z,'Marker','.','MarkerEdgeColor','r','MarkerSize',10, 'LineStyle', 'none')
axis off
% записываем координаты каждой точки в строку матрицы
X=[x y z];
% строим диаграмму Вороного
[V,C]=voronoin(X);
% отображаем содержимое массива V с координатами вершин диаграммы Вороного
V
% отображаем содержимое массива С с номерами вершин для каждой ячейки
for k=1:length(C)
    disp(C{k})
end
%проходим по ячейкам диаграммы
for k=1:length(C)
    if all(C{k}~=1) 
        % выбираем только конечные ячейки  (в множество их вершин не входит первая вершина, 
        % соответствующая бесконечно удаленной точке
        VertCell = V(C{k},:);
        % строим выпуклую оболочку вершин конечной ячейки
        KVert = convhulln(VertCell);
        % отображаем ее при помощи полигонального объекта с полупрозрачными гранями
        patch('Vertices',VertCell,'Faces',KVert,'FaceColor','g','FaceAlpha',0.5)
        pause(1)
    end
end

В результате получаем следующий массив V с координатами вершин диаграммы Вороного, в котором первая строка соответствует бесконечно удаленной точке

V =
       Inf       Inf       Inf
    1.6786   -1.3247   -0.8356
   -0.0413    1.4857    0.2195
    0.7505    0.4797    0.7848
    0.6998    0.4013    1.4033
    0.2062    0.3564    0.4473
    0.0184    0.8389    0.4706
    0.5013    0.5008    0.7148
    0.7016    1.5308    0.5868
    0.7426    0.7456    0.6857
    0.0122    1.4615    0.2417
    4.4104   -2.1563   -0.0837
    1.1532    0.2007    0.6310
    1.7639   -1.1568   -0.6208

Массив C с номерами вершин ячеек получается такой:

     1     4     5     9    10    12    13
     1     5     6     7     8
     2     3     4     6     7     8    10    11    13    14
     1     3     4     5     7     8     9    10    11
     1     2     3     9    11    12    14
     9    10    11    12    13    14
     1     2     4     5     6     8    12    13    14
     1     2     3     6     7

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

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

% генерируем случайные точки
rand('state',0);
x = rand(22,1);
y = rand(22,1);
% отображаем их красными маркерами
figure('Color','w')
plot(x,y,'Marker','.','MarkerEdgeColor','r','MarkerSize',10, 'LineStyle', 'none')
axis off
% находим вершины и многоугольники диаграммы Вороного
[V,C]=voronoin([x y]);
% проходим в цикле по многоугольникам 
for k=1:length(C)    
    if all(C{k}~=1)
        % конечные многоугольники рисуем полигональными объектами
        VertCell = V(C{k},:)
        patch('XData',VertCell(:,1),'YData',VertCell(:,2),'FaceColor','g','FaceAlpha',0.5)
        pause(0.1)
    end
end

Как связаны выпуклая оболочка, триангуляция Делоне и диаграмма Вороного

Мы рассмотрим этот вопрос на примере множества точек на плоскости (xi,yi)i=1,2,...,k.

Определим для каждой точки высоту zi = xi² + yi². Тогда точки с координатами (xi,yi,zi)i=1,2,...,k лежат на параболоиде
z(x,y) = x² + y² (на рисунке он обозначен желтым цветом).

Построим для множества точек с координатами (xi,yi,zi)i=1,2,...,k выпуклую оболочку

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

Команды, для получения этих картинок приведены ниже:

% задаем координаты точек на плоскости
x =[-0.3 -0.04 0.37 0.35 0.08 -0.34 -0.36 -0.13 0.15 0.02 -0.21]';
y =[0.29 0.38 0.25 -0.24 -0.39 -0.32 -0.02 0.08 0.05 -0.2 -0.21]';
% рисуем плоскость при помощи полигонального объекта
figure('Color','w')
patch('XData',[-0.6 0.6 0.6 -0.6], 'YData',[-0.6 -0.6 0.6 0.6], 'FaceColor','g','FaceAlpha',0.1)
% делаем оси невидимыми и устанавливаем нужный угол обзора
axis off
view(-8,12)
hold on
% рисуем точки маркерами
plot(x,y,'r.','MarkerSize',20)

	% строим параболоид
[X,Y]=meshgrid(-0.6:0.05:0.6);
Z=X.^2+Y.^2;
hs=surf(X,Y,Z,'EdgeColor','none','FaceColor','y','FaceAlpha',0.5)

% в цикле для каждой точки рисуем линию от исходной плоскости до параболоида и отмечаем
% на нем новые точки
for k=1:length(x)
line('XData',[x(k) x(k)],'YData',[y(k) y(k)], 'ZData', [0 x(k)^2+y(k)^2],...
       'Color','k','LineWidth',2,'LineStyle',':')
    plot3(x(k),y(k),x(k)^2+y(k)^2,'r.','MarkerSize',20)
end
	
	% строим выпуклую оболочку точек 
kvert=convhulln([x,y,x.^2+y.^2])
% визуализируем выпуклую оболочку
h=trimesh(kvert,x,y,x.^2+y.^2)
set(h,'Marker','.','MarkerSize',20,'MarkerEdgeColor','r')
% устанавливаем нужный угол обзора (наблюдатель под плоскостью)
view(0,-90)
% удаляем параболоид
delete(hs)

Построенный выше параболоид z(x,y) = x² + y² позволяет получить диаграмму Вороного для множества точек (xi,yi)i=1,2,...,k. Для этого проведем в каждой точке (xi,yi,zi)i=1,2,...,k касательную плоскость к параболоиду, на рисунке показана касательная плоскость в одной точке

Для каждой точки (xi,yi,zi)i=1,2,...,k уравнение касательной к параболоиду плоскости будет иметь следующий вид:

Несложно показать, что если некоторая точка (x,y) равноудалена от точек (xi,yi) и (xj,yj), то тогда zi(x,y) = zj(x,y), а если точка (x,y) дальше от (xi,yi), чем от некоторой другой точки (xj,yj), то тогда zi(x,y) > zj(x,y). Обратное тоже верно, т.е. если zi(x,y) > zj(x,y), то (x,y) дальше от (xi,yi), чем от (xj,yj) и из равенства zi(x,y) = zj(x,y) следует, что расстояния от (x,y) до (xi,yi) и (xj,yj) совпадают. Каждая плоскость zi(x,y) делит все пространство на верхнее и нижнее полупространства. Построив объединение верхних подпространств мы получим выпуклый многогранник. Проекция его граней на исходную плоскость и даст диаграмму Вороного.

Если известна диаграмма Вороного для точек (xi,yi,zi)i=1,2,...,k, то легко получить триангуляцию Делоне. Мы берем точки, находящиеся в смежных ячейках (не забывая про бесконечные ячейки) и соединяем их отрезками.

Синими линиями нарисована диаграмма Вороного для красных точек, а зелеными - триангуляция Делоне.

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

Нахождение площади многоугольника и пересечения четырехугольников

Площадь многоугольника с вершинами (xk,yk), где k = 1,2,...,n, находится по формуле:

в которой формально положено

(x0,y0) = (xn,yn).

При этом важен порядок, в котором задаются вершины.

В MATLAB для вычисления площади многоугольника имеется функция polyarea, обращение к которой выглядит следующим образом:

S = polyarea(x, y)

где x, y - вектора, содержащие координаты вершин в порядке обхода против часовой стрелки, например команды:

figure('Color','w')
x=[0 1 2 1 1 0];
y=[1 0 1 1 2 1];
patch(x,y,'g')
S=polyarea(x,y)
text(0.8, 0.5, ['S=' num2str(S)], 'FontSize', 20) text(0.8, 0.5, ['S=' num2str(S)], 'FontSize', 20)

приводят к отображению многоугольника при помощи полигонального объекта и выводят значение площади:

При задании координат вершин в векторах x и y следует придерживаться определенного направления обхода, скажем, против часовой стрелки, или по часовой стрелке. Если в предыдущем примере сменить порядок обхода на противоположный

S=polyarea(x(end:-1:1),y(end:-1:1))

то получим тот же самый результат.

Если задать многоугольник с пересекающимися ребрами, то значение площади многоугольника будет вычислено неверно, что демонстрируется в следующем примере, где вместо площади двух треугольников, образующих многоугольник, находится их разность и в результате получается площадь равная нулю:

figure('Color','w')
x=[0 2 0 2];
y=[0 0 2 2];
patch(x,y,'r')
S=polyarea(x,y)
text(0.8, 0.5, ['S=' num2str(S)], 'FontSize', 20)


Площадь равна нулю при неправильном задании вершин многоугольник

Точнее, выводится модуль разности площадей многоугольников (в данном случае это верхний и нижний треугольники), которые образованы вершинами, пронумерованными по часовой и против часовой стрелки. В этом примере верным является следующее задание многоугольной области

figure('Color','w')
x=[0 2 1 2 0 1];
y=[0 0 1 2 2 1];
patch(x,y,'g')
S=polyarea(x,y)
text(0.8, 0.5, ['S=' num2str(S)], 'FontSize', 20)


Если вершины многоугольника задать в правильном порядке, то площадь вычисляется верно

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

figure('Color','w')
X=[0 3 3 
   2 4 4
   2 4 4
   0 3 3];
Y=[0 0 3
   0 0 3
   2 1 4
   2 1 4];
patch(X,Y,'g')
S=polyarea(X,Y)
text(0.5, 1, ['S_1=' num2str(S(1))], 'FontSize', 20)
text(3.2, 0.5, ['S_2=' num2str(S(2))], 'FontSize', 20)
text(3.2, 3.5, ['S_3=' num2str(S(3))], 'FontSize', 20)


Вычисление площадей нескольких многоугольников

Если координаты вершин находятся в строках матрицы, то используется третий входной аргумент, который задает соответствующее измерение массивов: polyarea(X,Y,2). Впрочем, аналогичный результат получится и при следующем вызове функции polyarea(X',Y'), где в качестве входных аргументов polyarea передаются транспонированные матрицы (для транспонирования векторов и матриц в MATLAB предназначен апостроф).

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

[x    y    ширина    высота]

где x - абсцисса левого нижнего угла прямоугольника, а y - ордината левого нижнего угла прямоугольника. Для построения прямоугольников можно использовать низкоуровневую функцию rectangle.

Пример.
figure('Color','w')
a=[0 0 3 2];
b=[1 1 3 4];
rectangle('Position',a,'EdgeColor','r','LineWidth',3)
rectangle('Position',b,'EdgeColor','g','LineWidth',3)
S=rectint(a,b)
text(1, 1.5, ['area of intersection=' num2str(S)], 'FontSize', 16)


Площадь пересечения прямоугольников

Входные аргументы функции rectint могут быть также двумерными массивами A и B с четырьмя столбцами и различным числом строк, каждая строка которых содержит четыре элемента (x, y, ширина, высота) и задает один из прямоугольников. Тогда функция rectint вычисляет матрицу, компоненты которой содержат площади пересечений каждого прямоугольника, определяемого массивом A, с каждым прямоугольником, определяемым массивом B. Т.е., если размер массива A составляет N на 4, а размер массива B - M на 4, то

S = rectint(A,B)

возвращает прямоугольную матрицу S размера N на M и в элемент S(i,j) записана площадь пересечения i-го прямоугольника, определяемого i-ой строкой матрицы A, с j-ым прямоугольником, определяемым j-ой строкой матрицы B.

Нахождение ближайшей точки и симплекса, содержащего точку

Функция dsearch MATLAB позволяет для заданной точки плоскости с координатами определять ближайшую к ней точку из заданного множества точек на плоскости . Для этого предварительно требуется построить триангуляцию Делоне по набору точек при помощи функции delaunay (см. раздел Триангуляция Делоне).

Для примера при помощи функции rand сгенерируем набор из 10 точек, абсциссы и ординаты которых равномерно распределены на отрезке [0,1], запишем их значения в массивы x и y, соответственно, отобразим точки синими маркерами и построим триангуляцию Делоне

x=rand(1,10);
y=rand(1,10);
figure('Color','w')
plot(x,y,'Marker','.','MarkerEdgeColor','b','MarkerSize',35,'LineStyle','none')
axis equal
hold on
T = delaunay(x,y);
h=triplot(T,x,y);
set(h,'LineWidth',2)
Теперь зададим точку, для которой мы хотим определить ближайшую из заданного набора точек, и отобразим ее красным маркером
xx=0.5;
yy=0.5;
plot(xx,yy,'Marker','.','MarkerEdgeColor','r','MarkerSize',35,'LineStyle','none')
Далее вызовем функцию dsearch, которая возвращает номер ближайшей точки из набора к заданной точке с координатами (xx,yy)
k = dsearch(x,y,T,xx,yy)
и построим найденную точку зеленым маркером
xnear=x(k);
ynear=y(k);
plot(xnear,ynear,'Marker','.','MarkerEdgeColor','g','MarkerSize',35,'LineStyle','none')
В результате получим


Поиск ближайшей точки. Задана красная точка, найдена зеленая.

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

Часто возникает следующая задача. Триангуляция Делоне построена и задана еще некоторая точка. Требуется определить симплекс, которому принадлежит заданная точка плоскости. На плоскости эту задачу решает функция MATLAB, которая называется tsearch. Для пространств, размерности большей двух, предназначена функция tsearchn. Для построения триангуляции Делоне по заданным точкам плоскости используется функция delaunay, для точек трехмерного пространства - функция delaunay3 и для пространства произвольной размерности - функция delaunayn (см. раздел Триангуляция Делоне). Функции tsearch и tsearchn возвращают номер симплекса, содержащего заданную точку. Если заданная точка лежит вне выпуклой оболочки множества точек, то, разумеется, нужного симплекса не существует и функции tsearch и tsearchn возвращают NaN (не число - not a number).

Напомним, что выпуклая оболочка множества точек, принадлежащих n-мерному пространству, определяется как наименьшее выпуклое множество, содержащее все эти точки. Выпуклое множество, порожденное точками (Pi)i=1,2,...,k, определяется как множество точек:

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

Приведем пример, нахождения треугольника триангуляции, содержащего заданную точку, и закраски его цветом. Воспользуемся функцией rand для создания набора из 10 точек, абсциссы и ординаты которых равномерно распределены на отрезке [0,1], запишем их значения в массивы x и y, соответственно, и построим триангуляцию Делоне

x=rand(1,10);
y=rand(1,10);
figure('Color','w')
axis equal
hold on
T = delaunay(x,y)
h=triplot(T,x,y)
set(h,'LineWidth',2)

Зададим точку, как линейную комбинацию 1-ой, 5-ой и 7-ой точек с положительными коэффициентами, в сумме равными единице (для того, чтобы эта точка всегда оказывалась внутри выпуклой оболочки заданного множества точек)

xx=0.3*x(1)+0.5*x(5)+0.2*x(7);
yy=0.3*y(1)+0.5*y(5)+0.2*y(7);
и построим ее красным маркером
plot(xx,yy,'Marker','.','MarkerEdgeColor','r','MarkerSize',35)
Вызовем функцию tsearch для поиска нужного треугольника
k = tsearch(x,y,T,xx,yy)

запишем абсциссы и ординаты его вершин в массивы X и Y, соответственно, (каждая строка матрицы T содержит номера точек - вершин соответствующего треугольника) и построим частично прозрачный треугольник желтого цвета при помощи низкоуровневой функции patch, предназначенной для рисования полигональных объектов

X=[x(T(k,1)) x(T(k,2)) x(T(k,3))];
Y=[y(T(k,1)) y(T(k,2)) y(T(k,3))];
patch(X,Y,'y','FaceAlpha',0.3)
В результате получаем


Заданная точка в двумерном пространстве и содержащий ее симплекс

Если заданная точка лежит вне выпуклой оболочки множества точек,

xx=5;
yy=5;
то функция tsearch возвращает NaN
k = tsearch(x,y,T,xx,yy)
k =
   NaN

В качестве входных аргументов xx и yy функции tsearch могут быть и вектора, содержащие, соответственно, абсциссы и ординаты каждой из точек. Тогда выходным аргументом k функции tsearch будет вектор того же самого размера, что и xx, причем каждая компонента вектора k которого содержит или соответствующий номер треугольника, или NaN в том случае, если точка не принадлежит выпуклой оболочке.

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

k = tsearchn(X, T, XX)

где строки массива X содержат координаты точек исходного набора, T является матрицей с информацией о триангуляции исходного набора точек, а XX есть вектор с координатами точки, для которой требуется найти содержащий ее симплекс. Функция tsearchn позволяет также найти барицентрические координаты заданной точки и может принести пользу и в двумерном случае. Сначала мы рассмотрим пример о нахождении симплекса, содержащего заданную точку, а затем пример нахождения барицентрических координат заданной точки.

Сгенерируем набор из девяти точек в трехмерном пространстве, абсциссы, ординаты и аппликаты которых равномерно распределены на отрезке [0,1] и запишем координаты точек набора в строки массива X

rand('state',0);
X = rand(9,3);
Построим триангуляцию набора точек, используя функцию delaunay3 и визуализируем ее при помощи функции tetramesh
T = delaunay3(X(:,1), X(:,2), X(:,3));
figure('Color','w')
h=tetramesh(T,X);
set(h,'FaceColor','none','EdgeColor','b','LineWidth',2)
hold on
Возьмем точку, являющуюся линейной комбинацией 1-ой, 5-ой и 7-ой точек с положительными коэффициентами, в сумме равными единице (для того, чтобы эта точка всегда оказывалась внутри выпуклой оболочки заданного множества точек) и отобразим ее ее красным маркером
XX=0.3*X(1,:)+0.5*X(3,:)+0.2*X(7,:);
plot3(XX(1),XX(2),XX(3),'Marker','.','MarkerEdgeColor','r','MarkerSize',35)
Найдем номер содержащего эту точку симплекса, обратившись к функции tsearchn, и запишем его в переменную k
k = tsearchn(X,T,XX)
Рисуем ребра симплекса красными линиями и делаем его грани желтыми полупрозрачными
H=tetramesh(T(k,:),X);
set(H,'FaceColor','none','EdgeColor','r','FaceColor','y',...
    'FaceAlpha',0.3,'LineWidth',2)
Получаем


Заданная точка в трехмерном пространстве и содержащий ее симплекс

Как уже говорилось выше, функция tsearchn позволяет получить барицентрические координаты точки в симплексе (эти координаты еще называют естественными, или L-координатами). Барицентрические координаты были введены Августом Фердинандом Мёбиусом (1790 - 1868), немецким ученым, который занимался математикой и теоретической астрономией. Разберем, что такое барицентрические координаты, на двумерном примере, когда симплексами являются треугольники. Положение любой точки в треугольнике можно задавать тремя величинами, каждая из которых Lk является отношением расстояния Sk от точки до одной из его сторон к высоте Hk, опущенной на эту сторону из противоположной вершины, см. рисунок ниже для 1-ой вершины, т.е. для L1


Высота и расстояние от точки до стороны для определения барицентрических координат

Таким образом, барицентрические координаты принимают постоянные значения на линиях, параллельных соответствующей стороне, например, координата L1 принимает постоянные значения от 0 до 1 на следующих линиях (они нарисованы зеленым цветом):


Линии, на которых координата L1 принимает постоянные значения

Несложно показать, что барицентрические координаты каждой точки треугольника не являются независимыми (они и не могут быть независимыми, поскольку положение точки двумерного пространства всегда определяется ровно двумя координатами). Действительно, например координата L1 равна отношению площади заштрихованного треугольника к площади исходного треугольника (см. рисунок ниже)


Треугольник (заштрихованный), который связан с координатой L1

Сумма площадей таких треугольников для каждой из вершин равна площади исходного треугольника, поэтому

L1 + L2 + L3 = 1

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

  • матрицу в строках которой содержатся координаты точек исходного набора (в нашем случае у этой матрицы будет два столбца);
  • матрицу с информацией о триангуляции;
  • вектор с координатами точки, для которой требуется найти барицентрические координаты.

Снова сгенерируем 10 точек при помощи функции rand с координатами, равномерно распределенными на отрезке [0,1]

rand('state', 0)
x=rand(1,10);
y=rand(1,10);
figure('Color','w')
axis equal
hold on
построим триангуляцию, используя функцию delaunay
T = delaunay(x,y)
h=triplot(T,x,y)
set(h,'LineWidth',2)
Зададим теперь точку, которая принадлежит выпуклой оболочке заданного множества точек (для того, чтобы она принадлежала какому-либо треугольнику триангуляции), и отобразим ее на графике маркером. Пусть эта точка является линейной комбинацией первой, пятой и седьмой точек набора с коэффициентами, соответственно 0.2, 0.5 и 0.3
xx=0.2*x(1)+0.5*x(5)+0.3*x(7);
yy=0.2*y(1)+0.5*y(5)+0.3*y(7);
plot(xx,yy,'Marker','.','MarkerEdgeColor','r','MarkerSize',35)
и вызовем функцию tsearch с двумя выходными аргументами. Первый, как и в предыдущих примерах, содержит номер треугольника, которому принадлеит заданная точка, а второй является вектором барицентрических координат точки.
[k,p] = tsearchn([x' y'],T,[xx yy])
В результате получаем, что точка принадлежит второму треугольнику, т.к.
k =
     9
и ее барицентрические координаты имеют следующие значения.
p =

    0.1162    0.6572    0.2266
Непосредственная проверка позволяет убедиться, что сумма барицентрических координат действительно равна единице
sum(p)
ans =
     1

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

Для понимания того, что такое барицентрические координаты в случае большего числа измерений, достаточно вспомнить, что в n-мерном пространстве симплекс определяется как выпуклая линейная оболочка n+1-ой точек, не лежащих в одной гиперплоскости, т.е. для заданных точек x(1), x(2),..., x(n+1), симплексом является множество

Тогда для любой точки x из n-мерного пространства можно записать следующее представление в виде линейной комбинации

где

Числа Pk являются барицентрическими координатами точки заданной точки x из n-мерного пространства. Барицентрические координаты удобны при интерполировании сеточных функций, т.е. функций, значения которых определены в узлах некоторой сетки. В задаче интерполирования сеточных функций возникает необходимость построения таких базисных функций, которые 1) на каждом треугольнике являются полиномами от нескольких переменных (от двух для двумерных задач); 2) равны единице в одном из узлов треугольника и обращаются в ноль во всех остальных узлах треугольника. Линейные функции выглядят следующим образом


Линейные функции на треугольных элементах

и определяются через барицентрические координаты по формулам

Барицентрические координаты узлов таковы

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


Пример кусочно-линейной функции, используемой для интерполяции

Если задана сеточная функция, которая в i-ой точке из заданного набора принимает значение yi, то определенные выше функции Фi позволяют легко построить кусочно-линейный интерполянт Ф, который задается формулой

Для вычисления его в некоторой точке плоскости сначала определяют номер треугольника, которому принадлежит точка, затем находят ее барицентрические координаты и вычисляют значение интерполянта как сумму

в которой yi, yj, yk - значения заданной сеточной функции в соответствующих узлах треугольника. Связь между глобальной нумерацией узлов и локальной (1, 2, 3 в треугольнике) содержится в матрице, которую возвращает функция delaunay.

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

Квадратичная интерполяция требует задания узлов не только в вершинах треугольника, но и на его сторонах:


Треугольник с шестью узлами для квадратичной интерполяции

Барицентрические координаты узлов таковы:

Квадратичные функции, зависящие от барицентрических координат, определяются по следующим формулам:

Несложно проверить, что каждая равна единице в k-ом узле треугольника и обращается в ноль в остальных пяти узлах.

В случае кубической интерполяции требуется задание еще по одной точке на каждой стороне и дополнительной точки внутри треугольника, т.е. треугольник для кубической интерполяции содержит 10 узлов и имеет вид:


Треугольник с десятью узлами для кубической интерполяции

Для интерполяции применяются следующие кубические функции, зависящие от барицентрических координат

То, что каждая функция равна единице в k-ом узле треугольника и обращается в ноль в остальных девяти узлах проверяется несложной подстановкой барицентрических координат узлов в .

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

выражаются по следующим формулам:

Приближение разбросанных данных


Поиск по сайту:

Система Orphus

Яндекс.Метрика