Про быстрое пересечение полуплоскостей и диаграмму Вороного

Наконец-то дошли руки до написания этой статьи. Здесь я хочу в первую очередь поговорить о нахождении пересечения полуплоскостей за O(nlogn).

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

Итак, еще раз суть алгоритма в нескольких пунктах:

• Сначала добавим к нашим полуплоскостям 4 дополнительные - эти полуплоскости будут образовывать так называем bounding-box. Данное действие избавит от лишней возни с бесконечными пересечениями и т.п.
• Дальше нужно отсортировать все полуплоскости по углу их нормального вектора.
• Теперь будем добавлять полуплоскости в таком порядке и после каждого добавления будем поддерживать дэк полуплоскостей, которые образуют стороны многоугольника пересечения.
• Добавляя новую полуплоскость мы, возможно, должны выкинуть часть полуплоскостей с обоих концов дэка. 

Остается еще пару моментов, которые не оговорены в этом списке.

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

Рассмотрим такую ситуацию.

Добавилась черная полуплоскость и пересечение стало пусто.

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

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

Полуплоскость First не содержит в себе точку А - добавлять черную полуплоскость в дэк не надо.

Обратная ситуация - добавить надо.

Вроде бы это все. И в конце этой части статьи привожу код:

//You should normalize plane normal vectors for avoiding problems with precision
Plane pl[N], res[N];
int half(const Point &v)
{
 if (Gr(v.y, 0) || (Eq(v.y, 0) && Gr(v.x, 0)))
  return 1;
 return 2;
}

bool cmpPlane(const Plane &a, const Plane &b)
{
 if (half(a.v) != half(b.v))
  return half(a.v) < half(b.v);
 if (!Eq(a.v * b.v, 0))
  return a.v * b.v > 0;
 return b.contain(a.M);
}

bool eqPlane(const Plane &a, const Plane &b)
{
 return Eq(a.v * b.v, 0);
}

int halfPlaneIntersect(Point *pts)
{
 //Add bounding-box
 Point box[4];
 box[0] = Point(-MAXC, -MAXC);
 box[1] = Point(MAXC, -MAXC);
 box[2] = Point(MAXC, MAXC);
 box[3] = Point(-MAXC, MAXC);
 for (int i = 0; i < 4; i++) pl[cntPl++] = Plane(box[i], box[(i + 1) % 4] - box[i]);
 
 //Sort plane and delete planes with same normal vectors
 sort(pl, pl + cntPl, cmpPlane);
 cntPl = unique(pl, pl + cntPl, eqPlane) - pl;
 
 //Main part of algorithm
 int l = 0, r = 0;
 for (int i = 0; i < cntPl; i++)
 {
  while (r - l > 1 && !pl[i].contain(res[r - 1].getPoint(res[r - 2])))
   r--;
  while (r - l > 1 && !pl[i].contain(res[l].getPoint(res[l + 1])))
   l++;
  if (r - l > 0 && LsEq(res[r - 1].v * pl[i].v, 0))
   return 0;//Empty intersection
  if (r - l < 2 || res[l].contain(pl[i].getPoint(res[r - 1])))
   res[r++] = pl[i];
 }
 //Create result polygon
 int cntPts = 0;
 for (int i = l; i < r - 1; i++)
  pts[cntPts++] = res[i].getPoint(res[i + 1]);
 pts[cntPts++] = res[l].getPoint(res[r - 1]);
 return cntPts;
}

В приведенном коде опущены реализации классов Point и Plane.
Также нельзя оставлять без внимания возможные проблемы с точностью, избежать некоторые из которых можно предварительно отнормировав нормальные вектора плоскостей на единчную длину (спасибо Илье Кучумову за то, что обратил на это внимание).
Плоскость в данной реализации я задаю точкой (M), принадлежащей прямой, образующей плоскость, и направляющим вектором(v) (причем таким, что все точки, принадлежащие плоскости, удовлетворяют условию v * (P - M) >= 0; где * - векторное произведение).
Скажу также, что класс Plane содержит в себе два метода - getPoint(Plane) и contain(Point). Первый возвращает точку пересечения прямых, образующих плоскости; второй - проверяет, принадлежит ли данная точка полуплоскости (строгой полуплоскости).

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

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

Ура! Мы научились строить диаграмму Вороного за O(n²logn)!

Немного кода по диаграмме Вороного

В этом посте я хочу разобрать реализацию алгоритма построения диаграммы Вороного за O(n³).
Вспомним на словах саму суть алгоритма:

1. Для каждой точки, как обычно, строим ячейку Вороного. В данном случае мы хотим строить её за O(n²).
2. Будем поддерживать многоугольник, полученный пересечением первых i полуплоскостей (что это за полуплоскости и откуда они взялись можно узнать в предыдущем посте).
3. Добавлять новую полуплоскость будем за O(n), пользуясь тем фактом, что многоугольник пересечения не сильно меняется при добавлении новой полуплоскости.

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

Point hull[N];
Point globP, globV;

bool badPoint(const Point &A)
{
 return Ls((A - globP) % globV, 0); 
}

void addPoint(Point A, int pos, int &n)
{
 for (int i = n; i > pos; i--)
  hull[i] = hull[i - 1];
 hull[pos] = A;
 n++;
}

// A - точка, принадлежащая прямой, образующей полуплоскость
// u - нормальный вектор этой прямой
// r - количество точек в многоугольнике
void addHalfPlane(Point A, Point u, int &r) 
{             
 globP = A;
 globV = u;
 for (int i = 0; i < r; i++)
 {
  Point B = hull[i];
  Point C = hull[(i + 1) % r];
  Point P;
  if (intersectLine(B, C - B, A, u.ort(), P))
  {
   if (onSegment(B, C, P) && B != P && C != P)
    addPoint(P, i + 1, r);
  }
 }
 r = remove_if(hull, hull + r, badPoint) - hull;
}

Поясню некоторые моменты:

1. hull[] - массив, в котором хранится наш многоугольник пересечения
2. addPoint() - insert в массив
3. badPoint() - функция, которая возвращает true, если точка не принадлежит текущей полуплоскости

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

1. intersectLine(A, v, B, u, P) - возвращает true, если прямые пересекаются и сохраняет в P искомую точку. В противном случае, возвращает false.
2. onSegment(A, B, P) - проверяет, правда ли, что точка P лежит на отрезке [A;B]
3. Также, в данном коде используется переопределенный оператор % для двух точек, который обозначает скалярное произведение векторов.

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

void solve(Point A, int n)
{
 int r = 0;  
        hull[r++] = Point(0, 0);
 hull[r++] = Point(INF, 0);
 hull[r++] = Point(INF, INF);
 hull[r++] = Point(0, INF);

 for (int i = 0; i < n; i++)
 {
  if (p[i] == A)
   continue;
  Point M = (p[i] + A) / 2.;
  Point u = (A - p[i]);
  addHalfPlane(M, u, r);
 }
}

В начале я также добавил 4 точки - bounding box - которые будут ограничивать мою "бесконечную" плоскость. Это самый удобный и простой способ избавиться от "бесконечных" многоугольников и подобных неприятных вещей.

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

Немного слов про диаграмму Вороного.

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

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

Опишу несколько свойств диаграммы Вороного, которыми я скорее всего буду пользоваться дальше:

1.  Диаграмма Вороного - это планарный граф. С помощью теоремы Эйлера можно несложно доказать, что если диаграмма Вороного построена на множестве из n точек, то она содержит O(n) ребер и вершин (я знаю доказательство, которое использует тот факт, что степень каждой вершины хотя бы 3).
2. Ячейка диаграммы Вороного для точки P является пересечением полуплоскостей, образованных серединными перпендикулярами отрезков вида [A;P] (A - точка из нашего множества) и содержащих точку P.

Синим выделена ячейка Вороного для точки G
edge - это ребро диаграммы Вороного
vertex - вершина

Все приведенные ниже алгоритмы будут строить лишь множество ячеек Вороного, не объединяя их в общий планарный граф.
Соответственно, самый просто алгоритм за O(n⁴)

1. Зафиксируем точку, ячейку Вороного которой мы хотим построить - точка P.
2. Построим все серединные перпендикуляры отрезков вида [A;P].
3. Пересечем все серединные перпендикуляры между собой.
4. Для каждой точки пересечения проверим, что она принадлежит каждой полуплоскости. Точку, лежащую во всех полуплоскостях, назовем подходящей.
5. Отсортируем все подходящие точки по полярному углу.
6. Полученное множество точек - ячейка Вороного точки P.

    Видно, что для каждой точки ячейка построится за время O(n³) (для каждой из точек пересечения, которых порядка n², мы делаем проверку за n; алгоритм сортировки в таком случае не портит асимптотику).

    Из картинки видно, что ячейка Вороного может быть бесконечной. Данный случай можно определить однозначно по существованию "разрыва" - пары последовательных точек в нашем отсортированном списке, угол между которыми не меньше 180°.

угол обозначенный зеленым цветом - "разрыв"

    Далее разберем на словах алгоритм за O(n³).
В основе этого алгоритма лежит идея пересечения n полуплоскостей за квадрат. Как это делать?
Известно, что пересечение полуплоскостей - это выпуклый многоугольник (возможно бесконечный, как на рисунке). Тогда алгоритм пересечения за квадрат следующий:

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

Пересекаем многоугольник с зеленой полуплоскостью

Вуа-ля! Пересекли. Добавились точки J, K. От старого многоугольника остались точки E, D, C.

    Дальше - больше! Перед Вашими глазами сейчас предстанет алгоритм алгоритм за O(n²logn).
    Идея та же, что и в предыдущем алгоритме - пересечение полуплоскостей. Только можно сделать это оптимальнее - за nlogn. 
    Вспомним, что у нас есть. Мы для каждой точки хотим пересечь n полуплоскостей за nlogn. Заметим, что многоугольник пересечения полуплоскостей можно задавать массивом полуплоскостей, обладающим таким свойством - точка пересечения прямых, образующих пары соседних полуплоскостей - это вершина искомого многоугольника (см. рисунок для лучшего понимания).

Выделенный многоугольник можно задать последовательностью

 I, III, IV, V

    Идея алгоритма состоит в следующем:

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

Заметим, что наш массив - это стэк (т.к. мы только удаляем с конца и добавляем в конце).

Ниже приведена gif-анимация процесса построения пересечения полуплоскостей (они пронумерованы римскими числами в порядке сортировки).

В процессе добавления полуплоскости III мы должны удалить полуплоскость II, т.к. она добавляет плохую точку BADPOINT

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

1. Отсортировать все полуплоскости - O(nlogn)
2. Пробег по отсортированному массиву с добавлением и извлечением из стэка. Заметим, что каждая полуплоскость будет добавлена в стэк ровно 1 раз, соответственно удалена тоже не больше одного раза. Итого данная часть алгоритм работает за O(n).

    И наконец, последний алгоритм, который я хотел затронуть (и который я в силах осознать) - это алгоритм за O(n²).
 В данном алгоритме используется тот факт, что количество ребер в диаграмме Вороного - O(n).
Действительно, если бы мы научились находить каждое ребро одной ячейки Вороного за время T(n), то мы автоматически придумываем алгоритм построения диаграммы Вороного за O(n T(n)). Круто =)
    Давайте начнем придумывать. Для начала хочется от чего-то оттолкнуться. Действительно, найдем такую точку, которая по любому есть в ячейке Вороного для данной точки. Эта точка несложно ищется - давайте среди всех середин отрезков M = (середина [P;A], где P - выбранная точка, A - любая другая точка из исходного множества) выберем ближайшую. Утверждается, что эта точка точно присутствуем в ячейке для точки P (т.е. она лежит на границе многоугольника, который мы хотим найти). Если таких точек несколько - возьмем любую.

K - ближайшая точка для G

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

    Теперь поймем, как по прямой, содержащей ребро ответа и точке на ребре узнать следующее ребро. Ну это не сильно сложно. Давайте пересечем все остальные серединные перпендикуляры с текущей прямой. Возьмем 2 (или 1) точки пересечения - самый ближайшие точки с каждой из сторон к текущей точке на ребре (см. рисунок срочно!):

Разберем немного рисунок. Точка K - опорная точка. Мы пересекли четыре серединных перпендикуляра. P, Q - точки пересечения с прямой f, причем самые ближайшие к точке K. Понятно - что эти точки являются вершинами конечного многоугольника, а серединные перпендикуляры, их содержащие - следующими ребрами. Поэтому достаточно запустить рекурсивную функцию от изменившихся параметров - точки P и прямой n и точки Q и прямой j. 

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

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

    Таким образом, мы каждое ребро мы научились строить за O(n). Т.е. итоговая асимптотика такая, какую я заявил в начале рассказа про этот алгоритм =)

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