Category Archives: вычислительная геометрия

Алгоритм построения траектории при наличии фазовых ограничений II

Трубка внутри множества достижимости в обход препятствий

Трубка внутри множества достижимости в обход препятствий

Подчищенная и сокращённая версия (март 2015) предыдущего варианта.

Презентация (PDF)

К тому моменту я сделал тестовую реализацию алгоритма в двумерном пространстве (с визуализацией в 3D). К сожалению, запустить её на конференции не получилось, так как оказалось, что на Windows XP SP3 exe, скомпилированные Visual Studio 2013, не запускаются. На случай подобных проблем у меня, естественно, были заготовлены скриншоты:

Иллюстрации к слайдам 10–14

Иллюстрации к слайдам 20–22

Алгоритм построения траектории при наличии фазовых ограничений

Презентация по разработанному в этом году геометрическому алгоритму построения приближённой траектории линейной управляемой системы с выпуклым множеством управления и фазовыми ограничениями, заданными в виде наборов выпуклых многогранников. Алгоритм основан на двоичных разбиениях пространства (BSP).

Ссылка на презентацию

Ссылка на презентацию

Граница Парето в 3D за O(NlogN) время

Продолжая тему поста «Граница Парето в 2D за O(NlogN) или O(N) время«, я решил реализовать и опубликовать аналог в трёхмерном пространстве (обобщение на три показателя качества). Соответствующий алгоритм также описан в книге Препараты и Шеймоса «Вычислительная геометрия: Введение» вслед за двумерным алгоритмом. Алгоритм основан на принципе заметания объёма плоскостью.

  1. Отсортировать точки по координате z.
  2. Проходя точки в порядке сортировки от максимального значения координаты z до минимального, добавлять их в статус заметающей плоскости.

Статус плоскости хранит упорядочение (например, по координате x) подмножества доминирующих точек и представляет собой сечение плоскостью z = const области пространства, доминируемой перебранными на данный момент точками («лесенку», см. рисунок ниже). При добавлении новой точки P в статус:

  1. Ищем в статусе точку R, ближайшую справа к P.
    • Если она лежит выше P, то отбрасываем P и выходим.
    • Если её нет, либо она лежит ниже P, то вставляем P в статус и переходим к п.2.
  2. Ищем в статусе точку L, ближайшую слева к P и находящуюся выше P.
  3. Удаляем из статуса все точки в интервале (L, P). Если L не найдена, то удаляются все точки слева от P.
  4. Сообщаем P как недоминируемый элемент.

граница Парето 3D статус

Статус может быть реализован в виде сбалансированного двоичного дерева (желательно, прошитого, возможно, уложенного в массив), тогда, как нетрудно видеть, алгоритм найдёт все недоминируемые точки за O(NlogN) время. Была написана реализация на C++, позволяющая реализовать статус в виде политики, передаваемой параметром шаблона, либо воспользоваться std::set или любым совместимым аналогом.

Граница Парето в 2D за O(NlogN) или O(N) время

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

Предположим, что есть набор из N возможных решений {pi} некоторой задачи, и заданы показатели качества f1(p) → max и f2(p) → max, таким образом, каждому pi соответствует некоторая точка xi = (f1(pi), f2(pi)) на плоскости. Предположим, что заранее мы не можем решить, какой из показателей важнее, или как свести два показателя к одному общему, но мы можем заранее убрать из рассмотрения заведомо проигрышные решения pj, т.е. такие, для которых мы можем подобрать некоторое pk: f1(pj) ≤ f1(pk) и f2(pj) ≤ f2(pk). Тогда соответствующий набор оставшихся xk составит «границу Парето» исходного набора точек.

Очевидно, что мы можем сделать это за O(N2)-время, сравнив все значения показателей для решений попарно. Несколько лет назад мне понадобилось решать эту задачу в качестве подзадачи. Тогда я разработал для её решения алгоритм, работающий за O(NlogN)-время с помощью сортировки. Недавно мне пришлось вспомнить тот результат, кроме того, я решил обратиться к литературе по вычислительной геометрии и нашёл тот самый алгоритм в классической книге Препараты и Шеймоса «Вычислительная геометрия: Введение» (М.: Мир, 1989, с.194-195), где данная задача названа задачей «о поиске максимумов множества точек». Позволю себе высказать одно замечание. Для этого мне необходимо процитировать вводную часть из книги (с.192).

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

Далее в книге приведена формальная постановка задачи. Что интересно, указанная неформальная постановка относится к совсем другой (хотя и похожей внешне) задаче — задаче о выпуклой оболочке. Это станет очевидным, если её формализовать корректно: пусть vi — это количество дензнаков i-й валюты, ci — курс i-й валюты в некоторых условных единицах (в конце концов, всегда можно взять некоторую валюту за единицу, у неё ci = 1), тогда стоимость пакета v есть скалярное произведение векторов (vi) и (ci). Это будет максимизируемый показатель качества (см. опорная функция). Ниже показана иллюстрация с тремя «пакетами», где все три пакета являются элементами границы Парето, но «средний» из них не будет решением задачи о потенциальных королях Эревона (показано три возможных вектора).

(иллюстрация)

Итак, пусть заданы бинарные предикаты Sl, l=1..n, каждый из которых удовлетворяет аксиомам нестрогого порядка. Частным случаем таких предикатов могут быть операции сравнения точек по значению фиксированной координаты на «меньше или равно».

Определение. Элемент p доминирует элемент q, если Sl(q, p) для всех l.

Задача. Найти все p из набора pi, i=1..N, такие, что для любого j из 1..N : pj ≠ p ⇒ pj не доминирует p (для простоты формулировки задачи положим, что набор pi не содержит дубликатов, представленный ниже алгоритм отбрасывает дубликаты).

Алгоритм

Решение указанной задачи для случая n=2 (в частности, для точек на плоскости) основывается на том наблюдении, что если некоторый p находится «левее» (S1) некоторого q, то из недоминирования q над p следует, что p должен находиться «выше» (S2) q. Аналогично, некоторый r, находящийся «левее» p, должен быть «выше» p. Пусть процедура report(p) сообщает p в качестве одного из элементов решения.

  1. Переномеруем (отсортируем) pi так, чтобы для любого i=1..(N-1) и любого j=(i+1)..N было верно S1(pi, pj).
  2. Положим max=pN; report(max).
  3. Для i=(N-1)..1: если S1(max, pi), то max=pi; report(max).

граница Парето алгоритм 2D

Данный довольно простой алгоритм демонстрирует один из базовых приёмов вычислительной геометрии — «заметание» (плоскости прямой) в его простейшей форме.

Обобщённая реализация данного алгоритма на C++ здесь.

Простой анализ алгоритма даёт O(NlogN)-время для произвольного предиката S1 (сортировка с предикатом) и O(N)-память для хранения результата сортировки для последующего заметания. В случае, если pi можно отобразить в точки xi=(f1(pi), f2(pi)) на плоскости, заданные с ограниченной точностью, и сравнивать покоординатно, то сортировку можно выполнить за O(N)-время (например, с помощью блочной или поразрядной сортировок).

Дифференциальные игры 2

Доложился на семинаре Отдела динамических систем Института математики и механики УрО РАН, выкладываю презентацию и связанную с ней статью. Тема: численное построение решений (нэшевские решения) в линейных неантагонистических позиционных дифференциальных играх двух лиц с терминальными показателями качества в фазовом пространстве размерности больше двух (да, длинно, но могло быть ещё длиннее).

презентация

Работа в области дифференциальных игр

Пока работа по серии ParallelSort заморожена (и, возможно, будет перезапуск с корректировкой уже представленных результатов), я решил выложить презентации, представляющие мою текущую работу в должности научного сотрудника Института математики и механики УрО РАН в области построения численных решений в неантагонистических позиционных дифференциальных играх.

Первая презентация – с защиты моей кандидатской диссертации.

Автореферат кандидатской диссертации.

Вторая презентация была представлена на конференции «Математическая теория управления и математическое моделирование» в ИжГТУ (2012), с которой я только что вернулся.

Надо заметить, что ParallelSort имеет пересечение с этой работой – это отработка некоторых приёмов распараллеливания, которые можно применить и к описанным в презентациях алгоритмам (что немаловажно ввиду потенциально большого объёма обрабатываемых данных).

В дальнейшем я надеюсь оформить наработки в виде opensource проекта.

Библиотека вычислительной геометрии

Пожалуй, среди библиотек алгоритмов вычислительной геометрии на языке C++ наиболее широко известна достаточно развитая библиотека CGAL. В последнее время у меня возникло желание несколько «развернуть» логику построения библиотеки вычислительной геометрии для упрощения формулировки алгоритмов, которые мне приходилось реализовывать.

Базовый блок — вектор.

FixedVector: n-мерный вектор, размерность жёстко фиксирована;
BoundedVector: n-мерный вектор, размерность жёстко ограничена сверху.

Оба варианта не требуют управления памятью в духе std::vector, создавать их разумно в рамках пула (например, boost::pool_allocator) или в стеке.

CGAL различает точку и вектор в рамках аффинной структуры. Это заманчиво (проверка компилятором корректности выражений), но лично мне, скорее, мешало, поэтому пока что ограничимся только понятием вектора и работой с множествами векторов (слово «точка» ниже считается синонимом слова «вектор»).

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

  • преобразовать в (минимальную по количеству векторов) ортонормированную форму (для нормирования нужно отдельно задавать способ вычисления нормы или нормирования, непосредственно вектора можно сравнивать только на равенство, они «не умеют» считать свою длину);
  • определить разложимость по векторам базиса некоего вектора или другого базиса (принадлежность подпространству);
  • определить ортогональность базису вектора или другого базиса;
  • сложить с другим базисом (объединить вектора и упростить);
  • построить кобазис (вектора, недостающие базису подпространства до базиса полного пространства);
  • пересечение базисов = кобазис суммы кобазисов;
  • разность базисов = кобазис суммы кобазиса «уменьшаемого» с базисом «вычитаемого».

Система координат — базис и offset-вектор (можно считать точкой начала координат). Система координат задаёт линейное многообразие. Если базис заменить на кобазис (описывает левую часть системы линейных уравнений), а offset-вектор использовать как правую часть СЛУ, то это также можно использовать как описание линейного многообразия. Что удобнее, зависит от ситуации, в частности, гиперплоскость, скорее, удобнее описывать кобазисом, т.к. он состоит из одного уравнения / вектора нормали.

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

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

Алгоритмы: O(N log N + K) пересечение isobox‘ов (построение графа инцидентности), O(N log N) удаление доминируемых элементов.

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

GLU-1: Триангуляция многоугольников средствами GLU

Что такое GLU

Библиотека OpenGL Utility [html] [pdf] (сокращённо «GLU») является стандартным компонентом, поставляемым вместе с реализацией 3D API OpenGL. Как правило, доступна в составе ОС, поддерживающих OpenGL. Существуют открытые реализации (например, от SGI).

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

Программа-пример

Проект MSVC2010 под названием GLUPoly на основе FreeGLUT содержит реализацию описанного ниже метода работы с GLU для отображения загружаемых из файлов многоугольников.

Управление. Стрелки перемещают камеру, +/- увеличивает/уменьшает изображение, home возвращает в начальную позицию, enter запрашивает с консоли имя файла на отображение. На старте автоматически открывается input.txt из текущей директории.

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

Многоугольники

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

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

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

Триангуляция

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

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

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

Пересечение отрезков

Для корректного выполнения триангуляции необходимо найти все точки пересечения контуров. Эту задачу можно решать по-разному, например, простым попарным перебором звеньев всех контуров, что даёт O(N^2) сложность (где N — количество всех звеньев). Современные библиотеки используют более эффективные алгоритмы, имеющие сложность O(N log N + K), где K — количество точек пересечения. Например, алгоритмы, основанные на принципе заметания плоскости прямой, возникшие как улучшенные аналоги O( (N + K) log N) алгоритма Бентли-Оттманна (сюда можно отнести алгоритм Чазелле-Эдельсбруннера [acm], основанный на трапецеидальном разбиении плоскости). Для разрешения вырожденных случаев (вроде пересечения в одной точке сразу нескольких контуров) могут применяться «возмущения» (perturbations; небольшие изменения) координат точек.

Метод витков

Метод витков определяет, какие из треугольников, полученных в процессе триангуляции, включить в вывод алгоритма. Для каждого треугольника определено «число витков» — сумма полных оборотов (оборот против часовой стрелки добавляет единицу, оборот по часовой стрелке вычитает единицу), которые сделает вокруг внутренности треугольника точка при обхождении всех контуров. Например, можно брать все треугольники, число витков которых не равно нулю или выбирать треугольники, число витков которых нечётно (таково было поведение GLU по умолчанию в версиях 1.0 и 1.1, реализация GLU 1.2 от Microsoft также по умолчанию использует это правило, во многих случаях оно позволяет не заботиться о правильности ориентации исходных контуров). Задавая разные критерии, можно выполнять булевские теоретико-множественные операции над многоугольниками (в замыкании).

Создание объекта тесселятора GLU

Для триангуляции многоугольника нужно создать объект GLUtesselator с помощью вызова функции gluNewTess. Когда объект станет не нужен, его следует удалить, передав функции gluDeleteTess. Можно триангулировать разные многоугольники в разных потоках исполнения, используя разные объекты GLUtesselator.

Настройка тесселятора

Тесселятор поддерживает три «свойства», которые можно установить функцией gluTessProperty.

  • GLU_TESS_WINDING_RULE — критерий выбора треугольников по числу витков, варианты значений:
    • GLU_TESS_WINDING_ODD нечётное число,
    • GLU_TESS_WINDING_NONZERO число витков не равно нулю,
    • GLU_TESS_WINDING_POSITIVE число витков положительно,
    • GLU_TESS_WINDING_NEGATIVE число витков отрицательно,
    • GLU_TESS_WINDING_ABS_GEQ_TWO модуль числа витков больше единицы.
  • GLU_TESS_BOUNDARY_ONLY — выбор режима: получать результат в виде наборов треугольников (GL_FALSE) или в виде контуров, описывающих выбранные в соответствии с критерием области (GL_TRUE).
  • GLU_TESS_TOLERANCE — допустимая относительная погрешность, в случае ненулевого значения GLU может удалять слишком близко друг от друга отстоящие точки в рамках указанной погрешности.

Получение результирующих примитивов

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

  • GLU_TESS_BEGIN (_DATA) вызывается, когда начинается передача нового примитива, аналог glBegin;
  • GLU_TESS_VERTEX (_DATA) вызывается для передачи очередной вершины в виде указателя «vertexData», переданного в свою очередь тесселятору, аналог glVertexv;
  • GLU_TESS_END (_DATA) вызывается, когда примитив закончен, аналог glEnd;
  • GLU_TESS_COMBINE (_DATA) вызывается, когда необходимо породить новую вершину (найдена точка пересечения или удалены точки в соответствии с установленной допустимой погрешностью).

DATA-варианты позволяют передать колбеку дополнительный указатель «polygonData», и используются для привязки к конкретному объекту-получателю примитивов (что и сделано в примере).

Функция COMBINE нужна для того, чтобы дать пользователю возможность корректно «смешать» данные, прикреплённые к вершинам для порождения новых вершин (например, цвета и текстурные координаты) на основе четырёх исходных вершин. Минимальный пример этой функции дан в GLUPoly, более подробный пример здесь.

Процесс передачи контуров тесселятору

Для каждого отдельного многоугольника:
gluTessBeginPolygon (тесселятор, polygonData);

  для каждого контура:
  gluTessBeginContour (тесселятор);
    для каждой вершины контура:
    gluTessVertex (тесселятор, double координаты[3], vertexData);
  gluTessEndContour (тесселятор);
gluTessEndPolygon (тесселятор);

Третью координату вершины можно устанавливать в 0, триангуляция строится в проекции по направлению нормали, устанавливаемой функцией gluTessNormal. По умолчанию нормаль равна (0, 0, 1).

Структура программы-примера

Color.h содержит класс Color — представление RGBA цвета с float-компонентами («128-битный цвет»);

Vec2.h содержит класс Vec2 — двухкомпонентный вектор двойной точности;

Bbox2.h содержит класс Bbox2 — прямоугольник со сторонами параллельными осям координат, предназначенный для описания прямоугольных областей, занимаемых некими объектами, используется также для представления отображаемой области плоскости;

Poly.h + Poly.cpp определяют контур на основе стандартного vector;

Prim.h + Prim.cpp определяют класс-обёртку примитива OpenGL и класс Mesh2 (представляющий триангуляцию многоугольника), который обращается к GLU для триангуляции переданных контуров (см. шаблон конструктора, принимающий пару итераторов);

GLUPoly.cpp содержит собственно программу.