§ 15. Введение в вычислительную геометрию

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

Построение трехмерных многогранников

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

Трехмерный многогранник представляет собой множество треугольных граней, каждая из которых является смежной хотя бы одной другой грани. Вместе эти грани образуют поверхность или замкнутый объем. В программе такой многогранник хранится в двух массивах: массив вершин и массив граней. Массив точек содержит координаты всех точек, из которых образуются грани. Точки в этом массиве не повторяются, чтобы сэкономить память. Массив индексов содерижт индексы точек из первого массива, которые образуют грани. Для каждой грани используется по три индекса, а для экономии памяти можно использовать целочисленные типы данных по 1, 2, 4, 8 байт в зависимости от количества индексов. Для того чтобы массивы были совместимы с вызовами графических библиотек (OpenGL и пр.), координаты и индексы должны быть расположены в памяти непрерывно.

В дополнение к массивам вершин и граней используют массивы нормалей к поверхности многогранника. В зависимости от решаемой задачи нормали могут быть определены либо в центре каждой грани, либо в каждой вершине. Нормаль к поверхности — это вектор, который направлен перпендикулярно поверхности и длина которого равна единице. Нормаль к поверхности грани вычисляется как $\vec{n} = \frac{\vec{a} \times \vec{b}}{|\vec{a} \times \vec{b}|}$, где векторы $\vec{a}$ и $\vec{b}$ определяются из координат вершин треугольной грани как $\vec{a} = \vec{v}_2 - \vec{v}_1, \vec{b} = \vec{v}_3 - \vec{v}_1$, то есть это два вектора, исходящие из одной вершины.

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

Итого, при загрузке многогоранника из файла необходимо

• зафиксировать порядок обхода вершин граней,
• сгенерировать массивы нормалей,
• зафиксировать направление нормалей по отношению к объекту.

В следующем фрагменте кода вычисляются нормали в центре каждой грани и в каждой вершине.

const auto& vertices = this->_vertices;
const auto& faces = this->_faces;
const auto nfaces = faces.size();
auto& vnormals = this->_vertex_normals;
vnormals.assign(vertices.size(), T{0});
auto& fnormals = this->_face_normals;
fnormals.assign(faces.size(), T{0});
for (size_type i=0; i<nfaces; ++i) {
    const auto& face = this->_faces[i];
    const auto i0 = face[0];
    const auto i1 = face[1];
    const auto i2 = face[2];
    const auto& v0 = vertices[i0];
    const auto& v1 = vertices[i1];
    const auto& v2 = vertices[i2];
    vertex_type a = v1-v0, b = v2-v0;
    vertex_type n = cross(a, b);
    vnormals[i0] += n, vnormals[i1] += n, vnormals[i2] += n;
    fnormals[i] = unit(n);
}
for (auto& n : vnormals) { n = unit(n); }

Нормаль в вершине вычисляется как сумма нормалей к граням, в которые входит эта вершина.

Массовые свойства трехмерных многогранников

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

Объем вычисляется как объем тетраэдра, основание которого является гранью, а вершину можно выбрать произвольно, но у всех граней она должна быть одна и та же: $V = \frac{1}{6} (\vec{e}_1 \times \vec{e}_2) \cdot \vec{e}_3$, где $\vec{e}_1 = \vec{v}_2-\vec{v}_1, \vec{e}_2 = \vec{v}_3-\vec{v}_1, \vec{e}_3 = \vec{v}_4 - \vec{v}_1$, а $\vec{v}_1$ — вершина тетраэдра и $\vec{v}_{2,3,4}$ — основание. Объем, вычисленный таким образом, имеет знак, что позволяет точно вычислять объем невыпуклых объектов: объемы пересекающихся тетраэдров вычитаются друг из друга из-за различных направлений нормалей к граням по отношению к вершине тетраэдра. Центр масс вычисляется как взвешенная по объему сумма центров граней. Тензор инерции вычисляется аналогично центру масс, но суммируются не отдельные координаты $x,y,z$, а произведения их комбинаций: $x^2, y^2, z^2, xy, xz, yz$.

Телесный угол

Телесный угол, стянутый поверхностью $S$, относительно начала координат равен по величине площади сферы с единичным радиусом, покрытой проекцией поверхности $S$ на эту сферу. В общем случае телесный угол вычисляется как

$\Omega = \iint\limits_S \frac{\vec{n} \cdot dS}{r^2},$

где $dS$ — площадь бесконечно малой площадки поверхности, $\vec{n}$ — нормаль к площадке поверхности, $r$ — расстояние от начала координат до центра площадки.

Телесный угол, стянутый треугольником, относительно начала координат вычисляется как

$\Omega = 2 \, \text{atan} \frac{\vec{a}\cdot(\vec{b}\times\vec{c})} {|\vec{a}||\vec{b}||\vec{c}| + |\vec{a}|(\vec{b}\cdot\vec{c}) + |\vec{b}|(\vec{a}\cdot\vec{c}) + |\vec{c}|(\vec{a}\cdot\vec{b})}.$

Используя эту формулу, телесный угол трехмерного многогранника вычисляется как сумма телесных углов всех его граней. Сумма будет равна $4\pi$, если начало координат находится внутри многогранника, $2\pi$, если начало координат находится на одной из граней и нулю, если начало координат находится снаружи многогранника.

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

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

Триангуляция многоугольников

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

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

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

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

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

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

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

Метод Рупперта строится на следующем факте, известном из геометрии: если треугольник $ABC$ имеет угол $\angle BCA = \theta$, то угол $\angle BPA = 2\theta$, где $P$ является центром описанной окружности.

Метод отрезания ушей

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

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

• вторая по порядку вершина является выпуклой: внутренний по отношению к многоугольнику угол при этой вершине по величине меньше $\pi$,
• ребро, соединяющее первую и последнюю вершины полностью находится внутри многоугольника и
• внутри треугольника нет остальных вершин многоугольника.

Порядок точки относительно кривой

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

Для того чтобы вычислить порядок точки относительно многогранника, необходимо пройтись по всем ребрам и для каждого из них определить положение относительно точки. Положение ребра, образованного точками $\vec{a}$ и $\vec{b}$, относительно точки $\vec{p}$ определяется как $t = (\vec{a}-\vec{p})\times(\vec{b}-\vec{p})$. В двухмерном случае это скаляр. Если точка $\vec{p}$ находится не ниже точки $\vec{a}$ и ниже точки $\vec{b}$ и $t > 0$, то мы увеличиваем порядок на единицу. Если же точка $\vec{p}$ находится ниже точки $\vec{a}$ и не ниже точки $\vec{b}$ и $t < 0$, то мы уменьшаем порядок на единицу.

🌊🌊🌊

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

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

Задания

Порядок точки относительно кривой1 балл

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

Отрезание ушей2 балла

Напишите программу, которая триангулирует двухмерный многоугольник методом отрезания ушей.

Видео

2022

Доска

2022