Вычислительная геометрия в Python: от теории к применению

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

1. Вау, это звучит сложно.
2. Ах да, выпуклый корпус.

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

К чему вся эта суета?

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

Теоретически интересно…

С теоретической точки зрения вопросы вычислительной геометрии часто чрезвычайно интересны; ответы убедительные; и пути, по которым они достигаются, разнообразны. Одни только эти качества делают эту область достойной изучения, на мой взгляд.

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

Когда мы переводим это в вычислительную геометрическую нотацию, «план этажа» галереи представляет собой простой многоугольник. И с некоторым усилием локтя мы можем доказать, что n/3 камер всегда достаточно для полигона на n вершинах, каким бы грязным он ни был. Само доказательство использует двойственные графы, некоторую теорию графов, триангуляции и многое другое.

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

И важно на практике

Как я упоминал ранее, разработка игр в значительной степени зависит от применения вычислительной геометрии (например, обнаружение столкновений часто зависит от вычисления выпуклой оболочки набора объектов); как и географические информационные системы (ГИС), которые используются для хранения и выполнения вычислений по географическим данным; и робототехника тоже (например, для задач видимости и планирования).

Почему так тяжело?

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

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

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

 poly = Polygon([Point(0, 5), Point(1, 1), Point(3, 0), Point(7, 2), Point(7, 6), Point(2, 7)]) point = Point(5.5, 2.5) poly.contains(point)

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

Более абстрактно (и игнорируя необходимость представлять эти вещи в коде), проблемы, которые мы видим в этой дисциплине, очень трудно ригоризировать («сделать строгими») в алгоритме вычислительной геометрии. Как бы мы описали сценарий «точка в многоугольнике», не используя такой тавтологический язык, как «точка находится внутри многоугольника, если она находится внутри многоугольника»? Многие из этих свойств настолько фундаментальны и настолько фундаментальны, что их трудно определить конкретно.

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

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

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

Представляем CCW

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

В основе сюжета лежит обманчиво мощная примитивная операция: против часовой стрелки, или сокращенно «против часовой стрелки». (Сейчас я вас предупрежу: CCW будет всплывать снова и снова.)

CCW принимает три точки A, B и C в качестве аргументов и спрашивает: составляют ли эти три точки поворот против часовой стрелки (по сравнению с поворотом по часовой стрелке)? Другими словами, является ли угол A -> B -> C против часовой стрелки?

Например, зеленые точки — против часовой стрелки, а красные — нет:

Почему CCW имеет значение

CCW дает нам примитивную операцию, на которой мы можем строить. Это дает нам возможность начать строгие правила и решать задачи вычислительной геометрии.

Чтобы дать вам представление о его силе, давайте рассмотрим два примера.

Определение выпуклости

Первый: если задан многоугольник, можете ли вы определить, является ли он выпуклым? Выпуклость — бесценное свойство: знание того, что ваши многоугольники выпуклы, часто позволяет повысить производительность на порядки. В качестве конкретного примера: существует довольно простой алгоритм PIP, который работает за время Log(n) для выпуклых многоугольников, но не работает для многих вогнутых многоугольников.

Интуитивно этот разрыв имеет смысл: выпуклые формы «красивы», в то время как вогнутые формы могут иметь острые края, выступающие внутрь и наружу — они просто не следуют одним и тем же правилам.

Простой (но неочевидный) алгоритм вычислительной геометрии для определения выпуклости состоит в проверке того, что каждая тройка последовательных вершин является CCW. Это займет всего несколько строк кода геометрии Python (при условии, что points указаны в порядке против часовой стрелки — если points расположены по часовой стрелке, вы захотите, чтобы все триплеты были по часовой стрелке):

 class Polygon(object): ... def isConvex(self): for i in range(self.n): # Check every triplet of points A = self.points[i % self.n] B = self.points[(i + 1) % self.n] C = self.points[(i + 2) % self.n] if not ccw(A, B, C): return False return True

Попробуйте это на бумаге с несколькими примерами. Вы даже можете использовать этот результат для определения выпуклости. (Чтобы сделать вещи более интуитивными, обратите внимание, что кривая против часовой стрелки от A -> B -> C соответствует углу меньше 180, что является широко распространенным способом определения выпуклости.)

Пересечение линий

В качестве второго примера рассмотрим пересечение сегментов линии, которое также можно решить, используя только CCW:

 def intersect(a1, b1, a2, b2): """Returns True if line segments a1b1 and a2b2 intersect.""" return ccw(a1, b1, a2) != ccw(a1, b1, b2) and ccw(a2, b2, a1) != ccw(a2, b2, b1)

Почему это так? Пересечение отрезков можно также сформулировать так: для данного отрезка с конечными точками A и B лежат ли конечные точки C и D другого отрезка по одну сторону от AB? Другими словами, если повороты из A -> B -> C и A -> B -> D идут в одном направлении, сегменты не могут пересекаться. Когда мы используем этот тип языка, становится ясно, что такая проблема является хлебом с маслом CCW.

Строгое определение

Теперь, когда мы почувствовали важность CCW, давайте посмотрим, как она вычисляется. Учитывая точки A, B и C:

 def ccw(A, B, C): """Tests whether the turn formed by A, B, and C is ccw""" return (Bx - Ax) * (Cy - Ay) > (By - Ay) * (Cx - Ax)

Чтобы понять, откуда взялось это определение, рассмотрим векторы AB и BC. Если мы возьмем их векторное произведение, AB x BC, это будет вектор вдоль оси z. Но в каком направлении (т.е. +z или -z)? Как оказалось, если векторное произведение положительно, поворот идет против часовой стрелки; в противном случае по часовой стрелке.

Это определение будет казаться неинтуитивным, если вы не очень хорошо понимаете линейную алгебру, правило правой руки и т. д. Но именно поэтому у нас есть абстракция — когда вы думаете о CCW, просто подумайте о его интуитивном определении, а не о его вычислении. Значение будет сразу ясно.

Мое погружение в вычислительную геометрию и программирование с использованием Python

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

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

Алгоритм Киркпатрика

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

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

Хотя я не буду подробно обсуждать алгоритм, вы можете узнать больше здесь.

Минимальный ограничивающий треугольник

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

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

Практические советы, приложения и проблемы

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

На практике нам приходится иметь дело с целым рядом новых проблем.

Помните КВО? В качестве приятного продолжения давайте рассмотрим еще одно из его замечательных качеств: оно защищает нас от опасностей ошибок с плавающей запятой.

Ошибки с плавающей запятой: почему CCW является королем

В моем курсе вычислительной геометрии Бернар Шазель, уважаемый профессор, опубликовавший больше статей, чем я могу сосчитать, взял за правило не упоминать углы при попытке описать алгоритм или решение.

Почему? Углы грязные. Углы "грязные". Когда вам нужно вычислить угол, вам нужно разделить или использовать некоторую аппроксимацию (например, любую, связанную с числом Пи) или некоторую тригонометрическую функцию.

Когда вам нужно вычислить угол в коде , вы почти всегда будете приближаться. Вы будете отклоняться от какой-то крошечной степени точности с плавающей запятой, что важно, когда вы проверяете равенство. Вы можете найти некоторую точку на плоскости двумя разными способами и, конечно, ожидать, что p1.x == p2.x and p1.y == p2.y . Но на самом деле эта проверка часто дает сбой. Далее (и совершенно очевидно), эти точки будут потом иметь разные хэши.

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

Итак, что мы можем с этим сделать?

почти равно

Часть проблемы вычислительной геометрии Python заключается в том, что мы требуем точности в мире, где вещи редко бывают точными. Это станет проблемой чаще, чем при обработке углов. Рассмотрим следующее:

 # Define two random points p1 = RandomPoint() p2 = RandomPoint() # Take the line through them l1 = Line(p1, p2) # Shift both points up by sqrt(2) p1.y += sqrt(2) p2.y += sqrt(2) l2 = Line(p1, p2) # Slope 'should' be the same? if abs(l1.slope - l2.slope) > 0: print "Error!" # Error!

На самом деле этот код выведет «Ошибка!» примерно в 70% случаев (эмпирически). Мы можем решить эту проблему, немного смягчившись с нашим определением равенства; то есть, жертвуя степенью точности.

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

 def almostEqual(x, y, EPSILON=1e-5): return abs(x - y) < EPSILON class Point(object): ... def __eq__(self, that): return (almostEqual(self.x, that.x) and almostEqual(self.y, that.y))

На практике это очень помогает. Редко, если вообще когда-либо, вы будете вычислять две точки, которые отличаются менее чем на 1e-5, которые на самом деле должны быть разными точками. Я настоятельно рекомендую реализовать этот тип переопределения. Аналогичные методы можно использовать для линий, например:

 class Line(object): ... def __eq__(self, that): return (almostEqual(self.slope, that.slope) and almostEqual(self.intercept, that.intercept))

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

CCW - король

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

Вот конкретный пример: допустим, у вас есть некоторый выпуклый многоугольник P , вершина v и некоторая точка u вне многоугольника. Как узнать, пересекает ли прямая uv точку P выше или ниже v или вообще не пересекается за постоянное время?

Решение грубой силы (помимо линейного времени, а не постоянного) было бы проблематичным, поскольку вам нужно было бы вычислить некоторые точные точки пересечения линий.

Один подход с постоянным временем, который я видел, включает:

• Вычисление некоторых углов с помощью arctan2 .
• Преобразование этих углов в градусы путем умножения на 180/Pi.
• Изучение отношений между этими различными углами.

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

На мой взгляд, было бы лучше полностью избежать проблемы с ошибками с плавающей запятой. Если вы потратите несколько минут, чтобы посмотреть на проблему на бумаге, вы можете получить решение, полностью основанное на CCW. Интуиция: если вершины, смежные с v , находятся по одну сторону от uv , то прямая не пересекается; иначе посмотреть, находятся ли u и v по одну сторону линии между соседними вершинами, и, в зависимости от результата, сравнить их высоты.

Вот код Python для проверки пересечения выше v (пересечение ниже просто меняет направление сравнений):

 def intersectsAbove(verts, v, u): """ Returns True if uv intersects the polygon defined by 'verts' above v. Assumes v is the index of a vertex in 'verts', and u is outside of the polygon. """ n = len(verts) # Test if two adjacent vertices are on same side of line (implies # tangency) if ccw(u, verts[v], verts[(v - 1) % n]) == ccw(u, verts[v], verts[(v + 1) % n]): return False # Test if u and v are on same side of line from adjacent # vertices if ccw(verts[(v - 1) % n], verts[(v + 1) % n], u) == ccw(verts[(v - 1) % n], verts[(v + 1) % n], verts[v]): return uy > verts[v].y else: return uy < verts[v].y

Решение не сразу очевидно невооруженным глазом, но оно выражено на языке алгоритма вычислительной геометрии: «та же сторона линии» — классический элемент этого надежного алгоритма.

Сделано лучше, чем идеально

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

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

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

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

Вот код:

 class Polygon(object): ... def interiorPoint(self): """Returns a random point interior point""" min_x = min([px for p in self.points]) max_x = max([px for p in self.points]) min_y = min([py for p in self.points]) max_y = max([py for p in self.points]) def x(): return min_x + random() * (max_x - min_x) def y(): return min_y + random() * (max_y - min_y) p = Point(x(), y()) while not self.contains(p): p = Point(x(), y()) return p def contains(self, p): for i in range(self.n): p1 = self.points[i] p2 = self.points[(i + 1) % self.n] p3 = self.points[(i + 2) % self.n] if not ccw(p1, p2, p3): return False return True

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

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

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

Обратите внимание, что этот алгоритм требует от вас триангуляции вашего многоугольника, что сразу же накладывает на алгоритм другие ограничения времени выполнения, а также необходимость наличия библиотеки для триангуляции произвольных многоугольников (я использовал poly2tri с привязками Python).

 from p2t import CDT class Triangle(object): ... def area(self): return abs((Bx * Ay - Ax * By) + (Cx * By - Bx * Cy) + (Ax * Cy - Cx * Ay)) / 2 def interiorPoint(self): r1 = random() r2 = random() # From http://www.cs.princeton.edu/~funk/tog02.pdf return (1 - sqrt(r1)) * A + sqrt(r1) * (1 - r2) * B + r2 * sqrt(r1) * C class Polygon(object): ... def triangulate(self): # Triangulate poly with hole cdt = CDT(poly.points) triangles = cdt.triangulate() def convert(t): A = Point(tax, tay) B = Point(tbx, tby) C = Point(tcx, tcy) return Triangle(A, B, C) return map(convert, triangles) def interiorPoint(self): # Triangulate polygon triangles = self.triangulate() areas = [t.area() for t in triangles] total = sum(areas) # Calculate normalized areas probabilities = [area / total for area in areas] weighted_triangles = zip(triangles, probabilities) # Sample triangles according to area r = random() count = 0 for (triangle, prob) in weighted_triangles: count += prob # Take random point from chosen triangle if count > r: return triangle.interiorPoint()

Надеюсь, дополнительные усилия очевидны из кода. Помните: как говорят в Facebook, «сделано лучше, чем идеально». То же самое относится и к задачам вычислительной геометрии.

Визуальное и автоматизированное тестирование

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

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

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

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

 class TestLocator(unittest.TestCase): ... def runLocator(self, polygons): # Pre-process regions l = Locator(polygons) # Ensure correctness for polygon in polygons: # Test 100 random interior points per region for k in range(100): target = polygon.interiorPoint() target_polygon = l.locate(target) self.assertEqual(polygon, target_polygon) self.assertTrue(target_polygon.contains(target))

Затем мы могли бы использовать метод runLocator для разных наборов полигонов, что дало бы нам хорошо диверсифицированный набор тестов.

Решения с открытым исходным кодом

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

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

• poly2tri: отличная библиотека для быстрой триангуляции полигонов. Также поддерживает (и это часто критично) полигоны с дырками в них. Написанный на C++, poly2tri также имеет привязки к Python, и его довольно легко настроить и запустить. См. мой метод triangulate выше, чтобы познакомиться с вызовами функций.
• scipy.spatial: включает функции для вычисления выпуклых оболочек, триангуляции Делоне и т. д. Быстрый (как всегда), надежный и т. д. Примечание: мне показалось полезным использовать собственный тип данных Point с методом toNumpy : def np(self): return [self.x, self.y] . Затем я мог бы легко вызвать методы scipy.spatial, например: scipy.spatial.ConvexHull(np.array(map(lambda p: p.np()), points)) .
• OpenCV: библиотека компьютерного зрения с открытым исходным кодом имеет несколько хороших автономных модулей вычислительной геометрии. В частности, я какое-то время использовал функцию минимального охватывающего треугольника, прежде чем реализовывать ее самостоятельно.

Заключение

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

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

Если вы хотите узнать больше или у вас есть ко мне вопросы, со мной можно связаться по адресу [email protected].