Пользователь
0,0
рейтинг
2 июня 2012 в 00:13

Разработка → Построение минимальных выпуклых оболочек tutorial


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


Понятие минимальной выпуклой оболочки


Пусть на плоскости задано конечное множество точек A. Оболочкой этого множества называется любая замкнутая линия H без самопересечений такая, что все точки из A лежат внутри этой кривой. Если кривая H является выпуклой (например, любая касательная к этой кривой не пересекает ее больше ни в одной точке), то соответствующая оболочка также называется выпуклой. Наконец, минимальной выпуклой оболочкой (далее кратко МВО) называется выпуклая оболочка минимальной длины (минимального периметра). Я не проверял (похоже, это можно доказать от противного), но кажется очевидным, что минимальная оболочка просто обязана быть выпуклой. Все введенные понятия иллюстрируются следующим рисунком.

Главной особенностью МВО множества точек A является то, что эта оболочка представляет собой выпуклый многоугольник, вершинами которого являются некоторые точки из A. Поэтому задача поиска МВО в конечном итоге сводится к отбору и упорядочиванию нужных точек из A. Упорядочивание является необходимым по той причине, что выходом алгоритма должен быть многоугольник, т.е. последовательность вершин. Наложим дополнительно условие на порядок расположения вершин — направление обхода многоугольника должно быть положительным (напомню, что положительным называется обход фигуры против часовой стрелки).

Задача построения МВО считается одной из самых простых задач вычислительной геометрии, для нее существует много различных алгоритмов. Ниже мы рассмотрим два таких алгоритма — Грэхема (Graham scan) и Джарвиса (Jarvis march). Их описание иллюстрируется кодом на Питоне. Обоим алгоритмам потребуется функция rotate, побробно описанная в предыдущем моем посте. Напомню, что эта функция определяет, с какой стороны от вектора AB находится точка C (положительное возвращаемое значение соответствует левой стороне, отрицательное — правой).
def rotate(A,B,C):
  return (B[0]-A[0])*(C[1]-B[1])-(B[1]-A[1])*(C[0]-B[0])


Алгоритм Грэхема (Graham scan)


Этот алгоритм является трехшаговым. На первом шаге ищется любая точка в A, гарантированно входящая в МВО. Нетрудно сообразить, что такой точкой будет, например, точка с наименьшей x-координатой (самая левая точка в A). Эту точку (будем называть ее стартовой) перемещаем в начало списка, вся дальнейшая работа будет производиться с оставшимися точками. По некоторым соображениям, исходный массив точек A нами меняться не будет, для всех манипуляций с точками будем использовать косвенную адресацию: заведем список P, в котором будут хранится номера точек (их позиции в массиве A). Итак, первый шаг алгоритма заключается в том, чтобы первой точкой в P оказалась точка с наименьшей x-координатой. Код:
def grahamscan(A):
  n = len(A) # число точек
  P = range(n) # список номеров точек
  for i in range(1,n):
    if A[P[i]][0]<A[P[0]][0]: # если P[i]-ая точка лежит левее P[0]-ой точки
      P[i], P[0] = P[0], P[i] # меняем местами номера этих точек 

Второй шаг в алгоритме Грэхема — сортировка всех точек (кроме P[0]-ой), по степени их левизны относительно стартовой точки R=AP[0]. Будем говорить, что B<C, если точка С находится по левую сторону от вектора RB.

Для выпонения такого упорядочивания можно применять любой алгоритм сортировки, основанный на попарном сравнении элементов, например, быструю сортировку. По некоторым причинам (главная из которых — корявость* рук), я буду использовать сортировку вставками.
*я буду очень признателен тем, кто сможет мне объяснить, как применить в данном случае встроенную питоновскую сортировку...
Итак, сортировка вставками (не забываем про косвенную адресацию и про то, что нулевая точка не сортируется):
  for i in range(2,n):
    j = i
    while j>1 and (rotate(A[P[0]],A[P[j-1]],A[P[j]])<0): 
      P[j], P[j-1] = P[j-1], P[j]
      j -= 1

Результат сортировки можно проиллюстрировать следующим рисунком.

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

Переходим к третьему действию. Все, что нам осталось сделать, так это срезать углы. Для этого нужно пройтись по всем вершинам и удалить те из них, в которых выполняется правый поворот (угол в такой вершине оказывается больше развернутого). Заводим стек S (реально список) и помещаем в него первые две вершины (они, опять же, гарантированно входят в МВО).
  S = [P[0],P[1]]

Затем просматриваем все остальные вершины, и отслеживаем направление поворота в них с точки зрения последних двух вершин в стеке S: если это направление отрицательно, то можно срезать угол удалением из стека последней вершины. Как только поворот оказывается положительным, срезание углов завершается, текущая вершина заносится в стек.
  for i in range(2,n):
    while rotate(A[S[-2]],A[S[-1]],A[P[i]])<0:
      del S[-1] # pop(S)
    S.append(P[i]) # push(S,P[i])

В итоге в стеке S (который теперь можно рассматривать, как список) оказывается искомая последовательность вершин, причем в нужной нам ориентации, определяющая МВО заданного множества точек A.
  return S


Сложность первого и последнего шагов алгоритма является линейной по n (хотя в последнем случае имеется вложенный цикл, однако, каждая вершина внутри этого цикла ровно один раз заносится в стек, и не более одного раза оттуда удаляется), следовательно, сложность всего алгоритма определяется вторым шагом — сортировкой, именно поэтому сортировка вставкой оказывается не лучшим вариантом при больших n. Если ее заменить на быструю сортировку, то получим суммарную сложность алгоритма O(nlogn). Можно ли улучшить это время? Если алгоритм основан на попарном сравнении точек (как у нас), то доказано, что данная оценка в общем случае не улучшаема. С этой точки зрения алгоритм Грэхема оптимален. Тем не менее у него имеется не очень хорошая особенность — он не является адаптивным в том смысле, что не важно, сколько вершин в итоге войдет в МВО (три, пять, десять или n), все равно время будет линейно-логарифмическим. Такой адаптивностью обладает алгоритм Джарвиса, к рассмотрению которого мы плавно и переходим.
Полный код алгоритма Грэхема
def grahamscan(A):
  n = len(A) # число точек
  P = range(n) # список номеров точек
  for i in range(1,n):
    if A[P[i]][0]<A[P[0]][0]: # если P[i]-ая точка лежит левее P[0]-ой точки
      P[i], P[0] = P[0], P[i] # меняем местами номера этих точек 
  for i in range(2,n): # сортировка вставкой
    j = i
    while j>1 and (rotate(A[P[0]],A[P[j-1]],A[P[j]])<0): 
      P[j], P[j-1] = P[j-1], P[j]
      j -= 1
  S = [P[0],P[1]] # создаем стек
  for i in range(2,n):
    while rotate(A[S[-2]],A[S[-1]],A[P[i]])<0:
      del S[-1] # pop(S)
    S.append(P[i]) # push(S,P[i])
  return S


Алгоритма Джарвиса


Алгоритм Джарвиса (другое название — алгоритм заворачивания подарков) концептуально устроен проще алгоритма Грэхема. Он является двухшаговым и не требует сортировки. Первый шаг точно такой же — нам нужна стартовая точка, которая гарантированно входит в МВО, берем самую левую точку из A.
def jarvismarch(A):
  n = len(A)
  P = range(n)
  for i in range(1,n):
    if A[P[i]][0]<A[P[0]][0]: 
      P[i], P[0] = P[0], P[i]  

На втором шаге алгоритма строится МВО. Идея: делаем стартовую вершину текущей, ищем самую правую точку в A относительно текущей вершины, делаем ее текущей и т.д. Процесс завершается, когда текущей вновь окажется стартовая вершина. Как только точка попала в МВО, больше ее можно не учитывать. Поэтому заведем еще один список H, в котором в правильном порядке будут храниться вершины МВО. В этот список сразу же заносим стартовую вершину, а в списке P эту вершину переносим в конец (где мы ее в конце концов найдем и завершим алгоритм).
  H = [P[0]]
  del P[0]
  P.append(H[0])

Теперь организуем бесконечный цикл, на каждой итерации которого ищем самую левую точку из P относительно последней вершины в H. Если эта вершина стартовая, то прерываем цикл, если нет — то переносим найденную вершину из P в H. После завершения цикла в H находится искомая оболочка, которую мы и возвращаем в качестве результата.
  while True:
    right = 0
    for i in range(1,len(P)):
      if rotate(A[H[-1]],A[P[right]],A[P[i]])<0:
        right = i
    if P[right]==H[0]: 
      break
    else:
      H.append(P[right])
      del P[right]
  return H

Хм, мне удалось рассказать об алгоритме Джарвиса, не используя картинок. Следующий рисунок иллюстрирует все!

Оценим сложность алгоритма Джарвиса. Первый шаг линеен по n. Со вторым все интереснее. У нас имеется вложенный цикл, число внешних итераций равно числу вершин h в МВО, число внутренних итераций не превышает n. Следовательно, сложность всего алгоритма равна O(hn). Необычным в этой формуле является то, что сложность определяется не только длиной входных данных, но и длиной выхода (output-sensitive algorithm). А дальше как карты точки лягут. В худшем случае все точки из A входят в МВО (т.е. A уже само по себе выпуклый многоугольник), тогда h=n и сложность подскакивает до квадратичной. В лучшем случае (при условии, что точки из A не лежат на одной прямой) h=3 и сложность становится линейной. Осталось заранее понять, какой у нас случай, что сделать не так просто (если у вас нет машины времени**), можно только исходить из характера задачи — если точек много и они равномерно заполняют некоторую область, то (возможно) Джарвис будет быстрее, если же данные собраны на границе области, то быстрее будет Грэхем, как-то так…
**Машина времени вообще полезная штука с точки зрения алгоритмов, любая задача, требующая триллиона лет вычислений, с ее помощью может быть решена практически мгновенно — запускаем программу, садимся в машину времени, «летим» в будущее, считываем результат, возвращаемся назад. Осталось придумать, как обеспечить бесперебойную работу компьютера на пару триллионов лет...

Полный код алгоритма Джарвиса
def jarvismarch(A):
  n = len(A)
  P = range(n)
  # start point
  for i in range(1,n):
    if A[P[i]][0]<A[P[0]][0]: 
      P[i], P[0] = P[0], P[i]  
  H = [P[0]]
  del P[0]
  P.append(H[0])
  while True:
    right = 0
    for i in range(1,len(P)):
      if rotate(A[H[-1]],A[P[right]],A[P[i]])<0:
        right = i
    if P[right]==H[0]: 
      break
    else:
      H.append(P[right])
      del P[right]
  return H      


Вместо заключения


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

Спасибо всем за внимание!
Николай Ершов @nickme
карма
130,0
рейтинг 0,0
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Спецпроект

Самое читаемое Разработка

Комментарии (55)

  • 0
    Недавно для диплома впервые написал Грэхема, был счастлив по самое «не хочу»: Джарвиса — самый наивный алгоритм, который только можно представить.
    А насчёт выбора для использования — вообще без вариантов. Для общего набора Грэхем c O(NlogN) уделывает весьма специфичный квадратный алгоритм.
  • –11
    кажется очевидным, что минимальная оболочка просто обязана быть выпуклой

    С чего бы это?
    • +6
      Эта оболочка не минимальна, её можно укоротить.
    • +5
      Она не минимальная. Соедините самую левую точку с самой верхней напрямую — уменьшите периметр, ибо длина одной стороны треугольника в любом случае не больше суммы длин двух других сторон.
      • +4
        И правда, не подумал.
    • +2
      В статье минимальная оболочка определена как оболочка с наименьшим периметром. На этой картинке явно периметр можно уменьшить, убрав впадины.
    • 0
      Выпуклость минимальной оболочки — прямое следствие из неравенства треугольника
  • +1
    Интересно, познавательно и понятно изложенно. Спасибо! Буду ждать новых постов :)
  • 0
    Спасибо, картинки крутые, сразу понятна суть алгоритма.

    Что-то подобное в школе писал, только никаких алгоритмов не знал, сам придумал. Кажется я брал все пары точек, строил уравнение прямой через них и подставлял в уравнение координаты оставшихся точек. Если все получившиеся выражения имели одинаковый знак, то строил отрезок.
  • +1
    Я наверное Джарвис, потомучто именно этот алгоритм сразу пришел в голову глядя на вторую картинку из трех обрисованых множеств
  • +2
    Есть алгоритм пострения выпуклой оболочки методом «разделяй и властвуй». Тоже O(n log n), момедленнее гэхема на практике, но красиво!
    moais.imag.fr/membres/denis.trystram/SupportsDeCours/convexHull.pdf
    • +4
      там на самом деле минимум 3 подхода разделяй и властвуй. а к ним еще неплохо рассказать алгоритм построения опорной прямой. а если писать все методы, то надо экла-туссена еще расписать к примеру. но это уже получится листов 10-15 формата A4, с учетом того,, как подробно расписал их автор.

      а в дополнение к этим 2 достаточно известным яростро напрашивается алгоритм Чана, который использует грехэма и джарвиса.
      • 0
        а вообще, отлично, очень клевые картинки, все смотрится. заменить питон на псевдокод и можно в простой школе как материал давать. рад бы был если автор так же что-нить более сложное и интересное расписал.
      • +1
        Добавлю, что ассимптотика у алгоритма Чана O(n log h), где h — количество точек в выпуклой оболочке.
    • 0
      Он интересен в том смысле, что похожим подходом можно трехмерную выпуклую оболочку строить.
    • 0
      Спасибо за ссылку, почитаю… На самом деле на википедии есть описания многих алгоритмов, осталось найти время их поизучать: Выпуклая оболочка, Convex_hull и Convex_hull_algorithms
  • 0
    кажется очевидным, что минимальная оболочка просто обязана быть выпуклой

    Это элементарно доказывается с помощью свойства |AB| + |BC| >= |AC| (правило треугольника) и верно, кстати говоря, для любого нормированного пространства, а не только для плоскости.
    • +1
      Извиняюсь, в нормированном линейном пространстве.
    • 0
      Ага, неравенство треугольника я и имел в виду, когда писал про «очевидность» выпуклости. Думаю, что при строгом доказательстве для произвольных кривых возникнут некоторые трудности, нет? Для ломаных линий (т.е. для многоугольников, в контексте рассматриваемой темы) это не проблема, согласен…
      • +1
        Единственная трудность — придумать, куда воткнуть интеграл или предел.

        А дальше — то же самое неравенство треугольника.
      • 0
        Произвольных множеств точек, вы хотели сказать? Если множество конечно, то это все еще ломаная, если же бесконечно, но вложено в ограниченное подмножество плоскости, то тут веселее, понадобится предельные переход. Однако, если бесконечное множество точек на плоскости вложено в некоторое ограниченное ее подмножество, то никаких проблем нет: алгоритм точно такой же: индуктивное применение неравенства треугольника. Только индукция тут трансфинитная.
        • 0
          Можно проще и без пределов. Легко доказать что для любой не выпуклой оболочки существует выпуклая оболочка меньшей длины содержащая данную. Алгоритм построения этой оболочки — соединяем ближайшие две точки в которых нарушается выпуклость (касательная пересекает оболочку) прямой. Из того факта, что минимальное расстояние между двумя точками — отрезок прямой следует, что построенная оболочка имеет длину меньше чем данная.

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

    Приятно видеть подобные посты, пишите почаще
    • +1
      Я, пока пост готовил, тоже успел придумать суперэффективный алгоритм, при ближайшем рассмотрении оказалось, что это просто-напросто усложненная версия алгоритма Джарвиса :)
  • +2
    А можно пару примеров реального практического применения? В смысле где это требуется?
    • +1
      Ну я вот сейчас пишу бота для старкрафта. Там надо иметь общее представление о своей территории и о территории противника и ограничить снаружи, то есть нарисовать огибающую линию.

      Подобный алгоритм применяется для того чтобы запихать все свои здания в некоторый регион и наречь его потом «базой». Здания — это по сути точки на плоскости. Остальное думаю понятно.
      • НЛО прилетело и опубликовало эту надпись здесь
        • 0
          Нет, читал, чтобы быть в курсе, но писать не пробовал.
  • 0
    Я хочу заметить что с машиной времени можно решать любую задачу намного легче: считываем данные, считываем из будущего ответ Х, через минуту начинаем проверять этот Х, если он не подходит (можно убедиться в теории и за миллиард лет, но особенно хорошо алгоритм работает с задачами из класса NP, где ответ проверяется быстро), то отправляем в прошлое Х+1, если же подходит отправляем Х. Так как отправленное число должно совпасть с полученным, Х гарантированно будет правильным ответом на задачу.
    • 0
      Надо эту мысль обдумать на досуге. А вообще-то, «Назад в будущее» какое-то получается! :)
  • 0
    Не могу не добавить. То же самое в пару строчек на J: dr-klm.livejournal.com/42312.html

    :)
    • +1
      Мощно! Справедливости ради, не сомневаюсь, что приведенный выше питоновский код можно существенно подсократить, просто цель моего поста — объяснение алгоритма, а не демонстрация возможностей языка программирования :) В любом случае спасибо за ссылку, записал ее в избранное…
  • 0
    Спасибо, интересно!
    Не вижу больше ни одной статьи, кроме этой, с тегом «вычислительная геометрия».
    Может, про circle packing как-нибудь напишите?
    • 0
      В топике ошибка была в теге :(
  • +2
    По поводу сортировки на Питоне.
    wiki.python.org/moin/HowTo/Sorting

    def cmp_to_key(mycmp):
        'Convert a cmp= function into a key= function'
        class K(object):
            def __init__(self, obj, *args):
                self.obj = obj
            def __lt__(self, other):
                return mycmp(self.obj, other.obj) < 0
            def __gt__(self, other):
                return mycmp(self.obj, other.obj) > 0
            def __eq__(self, other):
                return mycmp(self.obj, other.obj) == 0
            def __le__(self, other):
                return mycmp(self.obj, other.obj) <= 0
            def __ge__(self, other):
                return mycmp(self.obj, other.obj) >= 0
            def __ne__(self, other):
                return mycmp(self.obj, other.obj) != 0
        return K
    


    base = P[0]
    del P[0]
    P.sort(key = cmp_to_key(lambda x, y: rotate(A[base], A[x], A[y])))
    P.insert(0, base)
    


    Как-то так.
    PS давно на питоне не писал, извиняюсь за глупые ошибки
    • 0
      Спасибо, сейчас опробуем!
      • +1
        Все заработало, отлично! Только порядок сортировки надо поменять:
          P.sort(key = cmp_to_key(lambda x, y: rotate(A[base], A[x], A[y])), reverse=True)
        

        Либо ключ (вернее, функцию сравнения) слегка изменить:
          P.sort(key = cmp_to_key(lambda x, y: rotate(A[base], A[y], A[x])))
        

        Работают оба варианта. А функция cmp_to_key, как написано по вашей ссылке, должна входить в модуль functools.
  • 0
    помню, на областной олимпиаде по информатик была подобная задача. через тангенсы решили.
    • 0
      Да. Это, впрочем, старая задача, я ее видел на олимпиаде Челябинской области в 1991 году.
      И решал через поиск максимального угла между соседними вершинами оболочки. И, что характерно, решил :)
    • 0
      Сейчас на студенческих олимпиадах эту задачу в чистом виде не встретишь.
      Обычно с применения алгоритма Грэхема следует начинать решение некоторых задач, чтобы перейти от произвольного множества точек к выпуклому многоугольнику.
      Например, в задаче нахождения круга минимального радиуса, содержащей данное множество точек.
      Или в задаче нахождения наибольшего по площади треугольника, вершинами которого являются точки из данного множества.
      • +1
        Тестовое задание в одно из подразделений CSoft: построить прямоугольник минимальной площади вокруг набора точек. Дополнительно предлагалось обобщить для 3х-мерного случая (для прямоугольного параллелепипеда).
      • НЛО прилетело и опубликовало эту надпись здесь
        • 0
          Еcть рандомизированный алгоритм за O(n). Кажется, оптимальнее не придумаешь.
          • НЛО прилетело и опубликовало эту надпись здесь
            • 0
              Welzl algorithm. По имени автора.
          • НЛО прилетело и опубликовало эту надпись здесь
  • 0
    С интересом бы почитал книгу на подобные темы, написанную понятным языком (как эта статья). Возможно ктото может посоветовать? И по распознаванию образов тоже было бы интересно.
    Спасибо за стастью.
    • +5
      Кормен, «Алгоритмы. Построение и анализ»?

  • 0
    Отличный стиль — читается легко, несмотря на ненависть к питону. Хочется добавки, интересует нахождение мво на растре, наборе дисплейных регулярных пикселей- там стопудово грехэм быстрее, но может есть нечто оригинальное?
  • 0
    С помощью чего вы делали картинки?
    • +1
      PowerPoint :)
  • НЛО прилетело и опубликовало эту надпись здесь
    • 0
      Если это вот этот алгоритм, то «Its average case complexity is considered to be O(nlogn), whereas in the worst case it takes O(n2) (quadratic)»… Может он самый быстрый в том же смысле, что и быстрая сортировка? Ссылкой не поделитесь? Нашел оригинальную статью.
      • НЛО прилетело и опубликовало эту надпись здесь
  • 0
    Ой спасибо тебе добрый человек! Очень пост выручил!

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