Метки

, , , , , ,

Постановка задачи

Нахождение пересечения сложных поверхностей часто требуется при моделировании различных процессов.  Одним из неаналитических способов представления поверхности является её аппроксимация набором треугольников различных размеров. Такой формат (faces & vertices) использует Matlab.

Рассмотрим пример, где необходимо найти пересечение двух изоповерхностей S1 и S2, соответствующих значениям S1Surf и S2Surf:
[x y z] = meshgrid(1:1:N1, 1:1:N2, 1:1:N3);
[f1, v1] = isosurface(x, y, z, S1, S1Surf);
[f2, v2] = isosurface(x, y, z, S2, S2Surf);


Здесь пары массивов (f1, f2) и (v1, v2) соответствуют так называемым faces и vertices. Эти форматы данных подробно рассмотрены в справке матлаба, поэтому лишь напомню, что v1 и v2 хранят наборы точек, соответствующих поверхностям, а f1 и f2 — тройки точек, образующих треугольники которыми аппроксимируются данные в примере изоповерхности.

Алгоритм Моллера

Для определения кривой пересечения S1 и S2 строго говоря, нужно перебрать все пары треугольников из различных поверхностей и проверить на вопрос их пересечения. Красивый алгоритм, отвечающий на вопрос, пересекаются ли два треугольника T1 (U^1U^2U^3) и T2 (V^1V^2V^3) предложил Томас Моллер [1], который и рассмотрим.

Вначале определим, не лежат ли точки одного из треугольников лишь по одну сторону от плоскости, которую задает второй (тут точка О соответствует началу координат):

N_1 = U_1U_2 \times U_1U_3,\quad d_1 = - N_1 \cdot OU_1.
d_{v2}^i = N_1 \cdot OV^i.

Если все d_{v2}^i одного знака, значит все точки первого треугольника лежат по одну сторону от плоскости второго и они не пересекаются.

Иначе ищем пересечение. Для этого вычисляем направляющий вектор прямой L = P_0 + tD пересечения плоскостей треугольников D = N_1 \times N_2 , где N_2 вычисляется аналогично N_1 , рассмотренному выше, а P_0 — некоторая точка, принадлежащая этой прямой. Моллером было предложено искать проекции вершин треугольников на прямую пересечения. Для этого вычисляем воспомогательные величины p_{v1,2} для каждого из треугольников по правилу
p_{v1} = U^i_j,\quad p_{v2} = V^i_j,\quad D_j = \max\{ |D_x|, |D_y|, |D_z| \}.
Тут j = x,y,z — ось, соответствующая максимальной по модулю проекции вектора D. Теперь пусть по разные стороны от пересечения плоскостей треугольников лежат точки с индексами 1 и 2, а 1 и 3 — по одну. Тогда можно вычислить значение линейной координаты t :
t_{1} = p_{v1} + d_{v2}^1 \dfrac{p_{v2} - p_{v1}}{d_{v2}^1 - d_{v1}^2}.
Аналогично вычисляем линейную координату для пары точек 3 и 2, и повторяем процедуру для второго треугольника.

Если отрезки на прямой L , задаваемые вычисленными координатами t пересекаются, значит пересекаются и треугольники.

Зная какие треугольники пересекаются, легко получить точки их пересечения, например пересечение сторон одного из треугольников со вторым [4].

Код алгоритма Моллера в Matlab

% Возвращает истину, если два треугольника пересекаются.
% Все точки задаются векторами-строками.
% Отрезок задается переменной line, line(1, : ) --- начало, line(2,: ) ---
% конец отрезка.
% Треугольники задаются как T1(1, : ) --- первая вершина первого
% треугольника, T1(2, : ) --- вторая вершина первого треугольника и т. д.
% Функция реализует алгоритм Tomas Moller.
%
% Входные параметры:
%   T1, T2 : Первый и второй треугольники в указанном выше формате.
%
% Выходные параметры:
%     flag : 1 если треугольники пересекаются и 0 --- если нет.
%        D : направляющий вектор прямой пересечения, нуль-вектор, если
%            такового нет.
%   points : массив точек пересечения сторон треугольников друг с другом.
%
function [flag, D, points] = getTTIntersection(T1, T2)
flag = false; D = zeros(3, 1); points = zeros(1, 3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Первейшая проверка на возможность перечения
%
N1 = cross(T1(2,: )-T1(1,: ), T1(3,: )-T1(1,: ));
d1 = - dot(N1, T1(1, : ));
dv2 = [dot(N1, T2(1,: )), dot(N1, T2(2,: )), dot(N1, T2(3,: ))] + d1;
if dv2 > 0
    return;
end
if dv2 < 0
    return;
end
N2 = cross(T2(2,: )-T2(1,: ), T2(3,: )-T2(1,: ));
d2 = - dot(N2, T2(1,: ));
% dv1 = [dot(N2, T1(1,: )), dot(N2, T1(2,: )), dot(N2, T1(3,: ))] + d2;
% if dv1 > 0
%     return;
% end
% if dv1 < 0 
%     return;
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Получение проекции
%
D = cross(N1, N2);
[mx, ind] = max(abs(D));
pv1 = [T1(1, ind) T1(2, ind) T1(3, ind)];
pv2 = [T2(1, ind) T2(2, ind) T2(3, ind)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Сравниваем интервалы. Тут нужны точки по разные стороны от отрезка
% пересечения треугольников, или что то же самое, по разные стороны от
% плоскости второго треугольника. Какая точка с какой стороны показывают
% знаки в массивах dv1, dv2.
%
if dv1(1)*dv1(2) > 0
    i1 = 1; i2 = 3; i3 = 2;
elseif dv1(1)*dv1(3) > 0
    i1 = 1; i2 = 2; i3 = 3;
elseif dv1(2)*dv1(3) > 0
    i1 = 2; i2 = 1; i3 = 3;
end
t11t = pv1(i1) + (pv1(i2) - pv1(i1))*dv1(i1)/(dv1(i1) - dv1(i2));
t21t = pv1(i3) + (pv1(i2) - pv1(i3))*dv1(i3)/(dv1(i3) - dv1(i2));
if dv2(1)*dv2(2) > 0
    i1 = 1; i2 = 3; i3 = 2;
elseif dv2(1)*dv2(3) > 0
    i1 = 1; i2 = 2; i3 = 3;
elseif dv2(2)*dv2(3) > 0
    i1 = 2; i2 = 1; i3 = 3;
end
t12t = pv2(i1) + (pv2(i2) - pv2(i1))*dv2(i1)/(dv2(i1) - dv2(i2));
t22t = pv2(i3) + (pv2(i2) - pv2(i3))*dv2(i3)/(dv2(i3) - dv2(i2));
t11 = min([t11t, t21t]); t21 = max([t11t, t21t]);
t12 = min([t12t, t22t]); t22 = max([t12t, t22t]);
if ((t11 >= t12) && (t11 <= t22)) || ((t21 >= t12) && (t21 <= t22)) || ...
        ((t12 >= t11) && (t12 <= t21)) || ((t22 >= t11) && (t22 <= t21))
    flag = true;
    % готовим точки пересечения
    ind = 1;
    p = getPoint(T1(1,: ), T1(2,: ), T2, N2);
    if size(p, 2) == 3
        points(ind, : ) = p;
        ind = ind + 1;
    end
    p = getPoint(T1(1,: ), T1(3,: ), T2, N2 );
    if size(p, 2) == 3
        points(ind, : ) = p;
        ind = ind + 1;
    end
    p = getPoint(T1(2,: ), T1(3,: ), T2, N2);
    if size(p, 2) == 3
        points(ind, : ) = p;
        ind = ind + 1;
    end
    p = getPoint(T2(1,: ), T2(2,: ), T1, N1);
    if size(p, 2) == 3
        points(ind, : ) = p;
        ind = ind + 1;
    end
    p = getPoint(T2(1,: ), T2(3,: ), T1, N1);
    if size(p, 2) == 3
        points(ind, : ) = p;
        ind = ind + 1;
    end
    p = getPoint(T2(2,: ), T2(3,: ), T1, N1);
    if size(p, 2) == 3
        points(ind, : ) = p;
    end
end
return;

Любые улучшения и конструктивная критика приветствуется 🙂

Полезные ссылки

  1. Tomas Moller «A Fast Triangle-Triangle Intersection Test».
  2. Набор различных алгоритмов на многие случаи жизни: http://progs-maker.narod.ru/algor.html.
  3. Координаты точки пересечения отрезка и плоскости: http://dubanov.exponenta.ru/russian/rus_head_1/rus_head_1.htm.
  4. Пересечение отрезка и треугольника: http://program.rin.ru/razdel/html/653.html.
  5. Изложение алгоритма Моллера на С++: http://softsurfer.com/Archive/algorithm_0105/algorithm_0105.htm.
Реклама