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

Наконец-то дошли руки до написания этой статьи. Здесь я хочу в первую очередь поговорить о нахождении пересечения полуплоскостей за 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)!