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

Пусть дано множество из [math]n[/math] отрезков и требуется найти все точки их пересечения. Очевидно, что задачу можно решить полным перебором за [math]O(n^2)[/math]; ясно также, что любой алгоритм будет в худшем случае работать за [math]\Theta(n^2)[/math] (нетрудно привести пример, когда количество пересечений квадратично, а алгоритм обязан сообщить о каждом пересечении). Однако существуют алгоритмы, которые оптимальнее, если количество точек пересечения отрезков невелико. Так алгоритм Бентли-Оттмана (англ. Bentley-Ottmann) позволяет решить задачу о пересечении отрезков, используя [math]O((n+I)\log{n})[/math] времени и [math]O(n)[/math] памяти, где [math]I[/math] — количество пересечений.

Формальная постановка задачи

Дано: [math]S[/math] — множество замкнутых отрезков на плоскости. Найти все точки их взаимного пересечения, охарактеризовав каждую точку пересечения множеством отрезков, которые участвуют в этом пересечении.

Описание алгоритма

Заметающая прямая и события

Воспользуемся методом заметающей прямой (sweep line), расположенной горизонтально и двигающейся вниз (в сторону уменьшения y-координаты). Нас будут интересовать события (events, event points) трёх типов:

• верхний конец отрезка,
• нижний конец отрезка,
• точка пересечения пары отрезков;

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

Приоритет событий

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

Нам потребуется две структуры данных.

Во-первых, мы будем хранить очередь событий [math]Q[/math] в виде сбалансированного бинарного дерева поиска, что позволит извлекать следующее событие и вставлять новое событие в очередь за [math]O(\log{m})[/math], где [math]m[/math] — количество элементов в очереди. Дерево будет отсортировано согласно порядку, введённому выше. Причём вместе с каждым событием мы будем хранить список отрезков, верхней точкой которых он является.

Во вторых, будем хранить статус [math]T[/math] заметающей прямой: множество отрезков, пересекающих заметающую прямую в данный момент времени, упорядоченных слева направо. От статуса нам потребуется оптимальные вставка и удаление отрезков, поэтому по-прежнему удобно воспользоваться бинарным деревом поиска.

Статус заметающей прямой

Соседние отрезки в статусе

Главная идея заключается в том, что мы будем проверять, пересекаются ли два отрезка, если они являются соседними в статусе. Это означает, что каждый отрезок мы будем проверять на пересечение с его левым и правым соседями в статусе. Далее, по ходу выполнения алгоритма, у отрезка могут измениться соседи; когда это происходит мы снова проверяем, не пересекает ли отрезок его новых соседей в статусе?

Далее приведен псевдокод алгоритма, а ниже подробно расписана обработка события определенного типа.

findIntersections(S)
   Инициализировать Q и T
   for s in S
      вставить концы s в Q (вместе с верхним концом вставить сам s)
   while not Q.empty()
      p = Q.pop()
      handleEventPoint(p)

Рассмотрим обработку событий.

• Верхний конец отрезка. В этом случае мы вставим отрезок в статус и проверим, не пересекает ли он соседние отрезки в статусе. Нас будут интересовать только пересечения ниже заметающей прямой. Естественно, если пересечения будут обнаружены, то они должны быть вставлены в [math]Q[/math].

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

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

Обработка верхних концов и пересечений

Обработка нижних концов

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

handleEventPoint(p)
   U(p) = множество отрезков, верхний конец которых есть p // Напомним, что мы храним такие отрезки в Q вместе c p
   // Далее найдём в T все отрезки, которые содержат p (они будут соседними в дереве)
   L(p) = множество отрезков, нижний конец которых есть p
   C(p) = множество отрезков, содержащих внутри себя p
   if [math] \vert L(p) \cup C(p) \cup U(p) \vert \gt  1[/math]
      report(p, [math]L(p) \cup C(p) \cup U(p)[/math]) // сообщить о p как о точки пересечения отрезков [math]L(p) \cup C(p) \cup U(p)[/math]
   T.remove([math]L(p) \cup C(p)[/math])
   T.insert([math]C(p) \cup U(p)[/math])
   // отрезки должны быть вставлены в порядке,
   // в котором они пересекают горизонтальную линию,
   // лежащую немного ниже заметающей прямой
   // (при удалении-вставке отрезков из C(p) - те поменяют порядок на обратный
   if [math]C(p) \cup U(p) = \emptyset[/math]
      s_l = левый сосед p в T
      s_r = правый сосед p в T
      findIntersection(s_l, s_r, p)
   else
      s1 = самый левый отрезок из [math]C(p) \cup U(p)[/math] в T
      s_l = левый сосед s1 в T
      findIntersection(s_l, s1, p)
      s2 = самый правый отрезок из [math]C(p) \cup U(p)[/math] в T
      s_r = правый сосед s2 в T
      findIntersection(s2, s_r, p)

findIntersection(s_l, s_r, p)
   if not intersects(s_l, s_r)
      return
   x = точка пересечения s_l и s_r
   if x лежит ниже заметающей прямой или на ней справа от p 
      Q.insert(x) // должно работать корректно, если x уже есть в Q

Поддержание статуса при обработке события

Доказательство корректности

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

Лемма:

Пусть [math]p[/math] — точка пересечения нескольких отрезков, причем [math]p[/math] не совпадает ни с одним из концов отрезков, участвующих в пересечении. Отсортируем отезки по углу вокруг [math]p[/math] и пусть отрезки [math]s_i[/math] и [math]s_j[/math] окажутся соседними в сортировке. Тогда существует событие [math]q[/math] с приоритетом выше, чем у [math]p[/math], такое что при обработке этого события [math]s_i[/math] и [math]s_j[/math] будут соседними в статусе.

Доказательство:

[math]\triangleright[/math]

Рассмотрим событие [math]q[/math] которое будет обработано непосредственно перед [math]p[/math]. Предположим, что при обработке этого события [math]s_i[/math] и [math]s_j[/math] не будут соседними в статусе. Это возможно только если между ними есть третий отрезок [math]s[/math]. Но тогда или [math]s[/math] пересекает какой-то из отрезков [math]s_i, s_j[/math] в точке с приоритетом выше, чем у [math]p[/math], но ниже, чем у [math]q[/math]; или [math]s[/math] пересекает их оба в точке [math]p[/math], что противоречит тому, что они соседние в сортировке по углу. Следовательно при обработке [math]q[/math] [math]s_i[/math] и [math]s_j[/math] будут соседними в статусе.

[math]\triangleleft[/math]

Теорема:

Алгоритм сообщает о всех точках пересечения

Доказательство:

[math]\triangleright[/math]

Воспользуемся индукцией по событиям, отсортированным в порядке, введённом выше ([math]p \lt q \Leftrightarrow p_y \lt q_y \lor p_y = q_y \land p_x \lt q_x[/math]). Пусть [math]p[/math] — точка пересечения. Предположим что все события [math]q, q \lt p[/math] были обработаны корректно. Тогда [math]p[/math] будет обнаружена. Действительно, если [math]p[/math] — конец некоторого отрезка, то [math]p[/math] добавлена в очередь событий в начале работы. Все содержащие её отрезки из статуса, который будет текущим при обработке [math]p[/math] будут найдены. Для остальных отрезков, содержащих [math]p[/math], верно, что [math]p[/math] — их верхний конец, и они будут найдены, т.к. мы храним их вместе с [math]p[/math] в очереди. Если же [math]p[/math] не является концом ни одного из отрезков, то по лемме найдётся событие [math]q[/math] с приоритетом выше, чем у [math]p[/math], такое что при обработке этого события пара отрезков, пересекающихся в [math]p[/math] будут соседними в статусе. Следовательно в этом случае мы также обнаружим [math]p[/math].

[math]\triangleleft[/math]

Оценка времени работы

Теорема:

Время работы алгоритма — [math]O((n+I)\log{n})[/math], где [math]I[/math] — количество пересечений.

Доказательство:

[math]\triangleright[/math]

Инициализация [math]Q[/math] может быть выполнена за [math]O(n\log{n})[/math], инициализация [math]T[/math] — за [math]O(1)[/math]. Далее, по ходу алгоритма мы обрабатываем каждое событие. Обработка включает в себя удаление события ([math]O(\log{n})[/math]), вызов функции findIntersection до двух раз, что может вызвать вставку новых событий [math]O(\log{n})[/math]. Также при обработке события [math]p[/math] мы [math]m(p) = \vert L(p) \cup C(p) \cup U(p) \vert[/math] раз выполняем операции вставки, удаления, поиска над [math]T[/math]. Каждая из этих операций требует [math]O(\log{n})[/math] времени. Пусть [math]m = \sum_{p} {m(p)}[/math]. Тогда время работы алгоритма — [math]O(m\log{n})[/math] Покажем, что [math]m = O(n + I)[/math], где [math]I[/math] — количество пересечений. Для этого рассмотрим планарный граф, вершинами которого являются концы отрезков, а также их точки пересечения, а ребрами — части отрезков, их соединяющие. Рассмотрим событие [math]p[/math]. Ему соответствует вершина графа со степенью [math]O(m(p))[/math], следовательно [math]m = O(\sum_{p} {deg(p)}) = O(E) = O(V) = O(2n + I)[/math], где [math]deg(p)[/math] — степень [math]p[/math], [math]E[/math] — число ребер в графе, [math]V[/math] — число вершин. (Предпоследний переход следует из формулы Эйлера.) Итак, время работы алгоритма: [math]O((n+I)\log{n})[/math].

[math]\triangleleft[/math]

Объём памяти

Очевидно, что статус в худшем случае занимает [math]O(n)[/math] памяти, однако очередь событий может занимать [math]O(n + I)[/math] памяти. Если модифицировать алгоритм так, что в очереди будут храниться только точки пересечения отрезков, соседних в статусе, то объём используемой памяти сократится до [math]O(n)[/math]. Для этого нужно удалять точки пересечения отрезков, которые перестали быть соседними. Перед тем, как мы дойдём до точки удаленной точки, обязательно найдётся соседняя в статусе пара отрезков, которые пересекаются в этой точке, и она снова будет вставлена в очередь. Ясно, что такая оптимизация не влияет на время работы алгоритма.

Источники

• Mark de Berg, Otfried Cheong, Marc van Kreveld, Mark Overmars (2008), Computational Geometry: Algorithms and Applications (3rd edition), Springer-Verlag, ISBN 978-3-540-77973-5 Chapter 2: Line Segment Intersection: pp.20–29.