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

В этом посте я хочу разобрать реализацию алгоритма построения диаграммы Вороного за 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 - которые будут ограничивать мою "бесконечную" плоскость. Это самый удобный и простой способ избавиться от "бесконечных" многоугольников и подобных неприятных вещей.

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