e-maxx_algo

Нахождение площади объединения треугольников. Метод вертикальной декомпозиции

Даны N треугольников. Требуется найти площадь их объединения.

Решение

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

Итак, у нас имеется N треугольников, которые могут как угодно пересекаться друг с другом. Избавимся от этих пересечений с помощью вертикальной декомпозиции: найдём все точки пересечения всех отрезков

(образующих треугольники), и отсортируем найденные точки по их абсциссе. Пусть мы получили некоторый массив B. Будем двигаться по этому массиву. На i-ом шаге рассматриваем элементы B[i] и B[i+1]. Мы имеем вертикальную полосу между прямыми X = B[i] и X = B[i+1], причём, согласно самому построению массива B, внутри этой полосы отрезки никак не пересекаются друг с другом. Следовательно, внутри этой полосы треугольники обрезаются до трапеций, причём стороны этих трапеций внутри полосы не пересекаются вообще. Будем двигаться по сторонам

этих трапеций снизу вверх, и складывать площади трапеций, следя за тем, чтобы каждый кусок был учитан ровно один раз. Фактически, этот процесс очень напоминает обработку вложенных скобок. Сложив площади трапеций внутри каждой полосы, и сложив результаты для всех полос, мы и найдём ответ - площадь объединения треугольников.

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

векторе элемент будет содержать нижнюю сторону самой нижней трапеции. Теперь мы просто идём по полученному вектору. Пусть i - текущий элемент; это означает, что i-ый кусок - это нижняя сторона некоторой трапеции, некоторого блока (который может содержать несколько трапеций), площадь которого мы хотим сразу прибавить к ответу. Поэтому мы устанавливаем некий счётчик треугольников равным 1, и поднимаемся по отрезкам вверх, и увеличиваем счётчик, если мы встречаем сторону какого-то треугольника в первый раз, и уменьшаем счётчик, если мы встречаем треугольник во второй раз. Если на каком-то отрезке j счётчик стал равным нулю, то мы

нашли верхнюю границу блока - на этом мы останавливаемся, прибавляем площадь трапеции, ограниченной отрезками i и j, и i присваиваем j+1, и повторяем весь процесс заново.

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

Реализация

struct segment {

int x1, y1, x2, y2;

};

struct point { double x, y;

};

struct item {

double y1, y2; int triangle_id;

};

struct line {

int a, b, c;

};

const double EPS = 1E-7;

void intersect (segment s1, segment s2, vector<point> & res) {

line l1 = { s1.y1-s1.y2, s1.x2-s1.x1, l1.a*s1.x1+l1.b*s1.y1 }, l2 = { s2.y1-s2.y2, s2.x2-s2.x1, l2.a*s2.x1+l2.b*s2.y1 };

double det1 = l1.a * l2.b - l1.b * l2.a; if (abs (det1) < EPS) return;

point p = { (l1.c * 1.0 * l2.b - l1.b * 1.0 * l2.c) / det1,

(l1.a * 1.0 * l2.c - l1.c * 1.0 * l2.a) / det1 };

if (p.x >= s1.x1-EPS && p.x <= s1.x2+EPS && p.x >= s2.x1-EPS && p.x <= s2.x2+EPS)

res.push_back (p);

}

double segment_y (segment s, double x) {

return s.y1 + (s.y2 - s.y1) * (x - s.x1) / (s.x2 - s.x1);

}

bool eq (double a, double b) { return abs (a-b) < EPS;

}

vector<item> c;

bool cmp_y1_y2 (int i, int j) { const item & a = c[i]; const item & b = c[j];

return a.y1 < b.y1-EPS || abs (a.y1-b.y1) < EPS && a.y2 < b.y2-EPS;

}

int main() {

int n; cin >> n;

vector<segment> a (n*3); for (int i=0; i<n; ++i) {

int x1, y1, x2, y2, x3, y3;

scanf ("%d%d%d%d%d%d", &x1,&y1,&x2,&y2,&x3,&y3); segment s1 = { x1,y1,x2,y2 };

segment s2 = { x1,y1,x3,y3 }; segment s3 = { x2,y2,x3,y3 }; a[i*3] = s1;

a[i*3+1] = s2; a[i*3+2] = s3;

}

for (size_t i=0; i<a.size(); ++i) if (a[i].x1 > a[i].x2)

swap (a[i].x1, a[i].x2), swap (a[i].y1, a[i].y2);

vector<point> b; b.reserve (n*n*3);

for (size_t i=0; i<a.size(); ++i)

for (size_t j=i+1; j<a.size(); ++j) intersect (a[i], a[j], b);

vector<double> xs (b.size()); for (size_t i=0; i<b.size(); ++i)

xs[i] = b[i].x;

sort (xs.begin(), xs.end());

xs.erase (unique (xs.begin(), xs.end(), &eq), xs.end());

double res = 0; vector<char> used (n); vector<int> cc (n*3); c.resize (n*3);

for (size_t i=0; i+1<xs.size(); ++i) { double x1 = xs[i], x2 = xs[i+1]; size_t csz = 0;

for (size_t j=0; j<a.size(); ++j) if (a[j].x1 != a[j].x2)

if (a[j].x1 <= x1+EPS && a[j].x2 >= x2-EPS) { item it = { segment_y (a[j],

x1), segment_y (a[j], x2), (int)j/3 };

cc[csz] = (int)csz; c[csz++] = it;

}

sort (cc.begin(), cc.begin()+csz, &cmp_y1_y2);

double add_res = 0;

for (size_t j=0; j<csz; ) { item lower = c[cc[j++]];

used[lower.triangle_id] = true; int cnt = 1;

while (cnt && j<csz) {

char & cur = used[c[cc[j++]].triangle_id]; cur = !cur;

if (cur) ++cnt; else --cnt;

}

item upper = c[cc[j-1]];

add_res += upper.y1 - lower.y1 + upper.y2 - lower.y2;

}

res += add_res * (x2 - x1) / 2;

}

cout.precision (8); cout << fixed << res;

}

Проверка точки на принадлежность выпуклому многоугольнику

Дан выпуклый многоугольник с N вершинами, координаты всех вершин целочисленны (хотя это не меняет суть решения); вершины заданы в порядке обхода против часовой стрелки (в противном случае нужно просто отсортировать их). Поступают запросы - точки, и требуется для каждой точки определить, лежит она внутри

этого многоугольника или нет (границы многоугольника включаются). На каждый запрос будем отвечать в режиме on-line за O (log N). Предварительная обработка многоугольника будет выполняться за O (N).

Алгоритм

Решать будем бинарным поиском по углу.

Один из вариантов решения таков. Выберем точку с наименьшей координатой X (если таких несколько, то выбираем самую нижнюю, т.е. с наименьшим Y). Относительно этой точки, обозначим её Zero, все остальные вершины многоугольника лежат в правой полуплоскости. Далее, заметим, что все вершины многоугольника уже упорядочены по углу относительно точки Zero (это вытекает из того, что многоугольник выпуклый, и уже упорядочен против часовой стрелки), причём все углы находятся в промежутке (-π/2 ; π/2].

Пусть поступает очередной запрос - некоторая точка P. Рассмотрим её полярный угол относительно точки Zero. Найдём бинарным поиском две такие соседние вершины L и R многоугольника, что полярный угол P лежит между полярными углами L и R. Тем самым мы нашли тот сектор многоугольника, в котором лежит точка P, и нам остаётся только проверить, лежит ли точка P в треугольнике (Zero,L,R). Это можно сделать, например, с помощью Ориентированной площади треугольника и Предиката "По часовой стрелке", достаточно посмотреть, по

часовой стрелке или против находится тройка вершин (R,L,P).

Таким образом, мы за O (log N) находим сектор многоугольника, а затем за O (1) проверяем принадлежность точки треугольнику, и, следовательно, требуемая асимптотика достигнута. Предварительная обработка

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

Замечания по реализации

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

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

Заметим, что полярный угол точки (X,Y) относительно начала координат однозначно определяется дробью Y/X,

при условии, что точка находится в правой полуплоскости. Более того, если у одной точки полярный угол меньше, чем у другой, то и дробь Y1/X1 будет меньше Y2/X2, и обратно.

Таким образом, для сравнения полярных углов двух точек нам достаточно сравнить дроби Y1/X1 и Y2/X2, что уже можно выполнить в целочисленной арифметике.

Реализация

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

struct pt {

int x, y;

};

struct ang {

int a, b;

};

bool operator < (const ang & p, const ang & q) { if (p.b == 0 && q.b == 0)

return p.a < q.a;

return p.a * 1ll * q.b < p.b * 1ll * q.a;

}

long long sq (pt & a, pt & b, pt & c) {

return a.x*1ll*(b.y-c.y) + b.x*1ll*(c.y-a.y) + c.x*1ll*(a.y-b.y);

}

int main() {

int n; cin >> n;

vector<pt> p (n); int zero_id = 0;

for (int i=0; i<n; ++i) {

scanf ("%d%d", &p[i].x, &p[i].y);

if (p[i].x < p[zero_id].x || p[i].x == p[zero_id].x && p[i].y

< p[zero_id].y)

zero_id = i;

}

pt zero = p[zero_id];

rotate (p.begin(), p.begin()+zero_id, p.end()); p.erase (p.begin());

--n;

vector<ang> a (n);

for (int i=0; i<n; ++i) {

a[i].a = p[i].y - zero.y; a[i].b = p[i].x - zero.x; if (a[i].a == 0)

a[i].b = a[i].b < 0 ? -1 : 1;

}

for (;;) {

pt q; // очередной запрос bool in = false;

if (q.x >= zero.x)

if (q.x == zero.x && q.y == zero.y) in = true;

else {

ang my = { q.y-zero.y, q.x-zero.x }; if (my.a == 0)

my.b = my.b < 0 ? -1 : 1; vector<ang>::iterator it = upper_bound (a.

begin(), a.end(), my);

if (it == a.end() && my.a == a[n-1].a && my.

b == a[n-1].b)

it = a.end()-1;

if (it != a.end() && it != a.begin()) { int p1 = int (it - a.begin()); if (sq (p[p1], p[p1-1], q) <= 0)

in = true;

}

}

puts (in ? "INSIDE" : "OUTSIDE");

}

}

Нахождение вписанной окружности в выпуклом многоугольнике с помощью тернарного поиска

Дан выпуклый многоугольник с N вершинами. Требуется найти координаты центра и радиус наибольшей вписанной окружности.

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

Алгоритм

Определим функцию Radius (X, Y), возвращающую радиус вписанной в данный многоугольник окружности с центром в точке (X;Y). Предполагается, что точки X и Y лежат внутри (или на границе) многоугольника. Очевидно, эту функцию легко реализовать с асимптотикой O (N) - просто проходим по всем сторонам многоугольника, считаем для каждой расстояние до центра (причём расстояние можно брать как от прямой до точки, не обязательно рассматривать как отрезок), и возвращаем минимум из найденных расстояний - очевидно, он и будет наибольшим радиусом.

Итак, нам нужно максимизировать эту функцию. Заметим, что, поскольку многоугольник выпуклый, то эта функция будет пригодна для тернарного поиска по обоим аргументам: при фиксированном X0 (разумеется, таком,

что прямая X=X0 пересекает многоугольник) функция Radius(X0, Y) как функция одного аргумента Y будет сначала возрастать, затем убывать (опять же, мы рассматриваем только такие Y, что точка (X0, Y)

принадлежит многоугольнику). Более того, функция max (по Y) { Radius (X, Y) } как функция одного аргумента X будет сначала возрастать, затем убывать. Эти свойства ясны из геометрических соображений.

Таким образом, нам нужно сделать два тернарных поиска: по X и внутри него по Y, максимизируя значение функции Radius. Единственный особый момент - нужно правильно выбирать границы тернарных поисков, поскольку вычисление функции Radius за пределами многоугольника будет некорректным. Для поиска по X никаких сложностей нет, просто выбираем абсциссу самой левой и самой правой точки. Для поиска по Y находим те отрезки многоугольника, в которые попадает текущий X, и находим ординаты точек этих отрезков при абсциссе X (вертикальные отрезки не рассматриваем).

Осталось оценить асимптотику. Пусть максимальное значение, которое могут принимать координаты - это C1, а требуемая точность - порядка 10-C2, и пусть C = C1 + C2. Тогда количество шагов, которые должен будет совершить каждый тернарный поиск, есть величина O (log C), и итоговая асимптотика получается: O (N log2 C).

Реализация

Константа steps определяет количество шагов обоих тернарных поисков.

В реализации стоит отметить, что для каждой стороны сразу предпосчитываются коэффициенты в уравнении прямой,

исразу же нормализуются (делятся на sqrt(A2+B2)), чтобы избежать лишних операций внутри тернарного поиска.

const double EPS = 1E-9; int steps = 60;

struct pt {

double x, y;

};

struct line {

double a, b, c;

};

double dist (double x, double y, line & l) { return abs (x * l.a + y * l.b + l.c);

}

double radius (double x, double y, vector<line> & l) { int n = (int) l.size();

double res = INF;

for (int i=0; i<n; ++i)

res = min (res, dist (x, y, l[i])); return res;

}

double y_radius (double x, vector<pt> & a, vector<line> & l) { int n = (int) a.size();

double ly = INF, ry = -INF; for (int i=0; i<n; ++i) {

int x1 = a[i].x, x2 = a[(i+1)%n].x, y1 = a[i].y, y2 = a

[(i+1)%n].y;

if (x1 == x2) continue;

if (x1 > x2) swap (x1, x2), swap (y1, y2); if (x1 <= x+EPS && x-EPS <= x2) {

double y = y1 + (x - x1) * (y2 - y1) / (x2 - x1); ly = min (ly, y);

ry = max (ry, y);

}

}

for (int sy=0; sy<steps; ++sy) { double diff = (ry - ly) / 3;

double y1 = ly + diff, y2 = ry - diff;

double f1 = radius (x, y1, l), f2 = radius (x, y2, l); if (f1 < f2)

ly = y1;

else

ry = y2;

}

return radius (x, ly, l);

}

int main() {

int n; vector<pt> a (n);

... чтение a ...

vector<line> l (n);

for (int i=0; i<n; ++i) {

l[i].a = a[i].y - a[(i+1)%n].y; l[i].b = a[(i+1)%n].x - a[i].x;

double sq = sqrt (l[i].a*l[i].a + l[i].b*l[i].b); l[i].a /= sq, l[i].b /= sq;

l[i].c = - (l[i].a * a[i].x + l[i].b * a[i].y);

}

double lx = INF, rx = -INF; for (int i=0; i<n; ++i) {

lx = min (lx, a[i].x); rx = max (rx, a[i].x);

}

for (int sx=0; sx<stepsx; ++sx) { double diff = (rx - lx) / 3;

double x1 = lx + diff, x2 = rx - diff;

double f1 = y_radius (x1, a, l), f2 = y_radius (x2, a, l); if (f1 < f2)

lx = x1;

else

rx = x2;

}

double ans = y_radius (lx, a, l); printf ("%.7lf", ans);

}

Нахождение вписанной окружности в выпуклом многоугольнике методом

"сжатия сторон" ("shrinking sides") за

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

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

Спасибо Ивану Красильникову (mf) за описание этого красивого алгоритма.

Алгоритм

Итак, дан выпуклый многоугольник. Начнём одновременно и с одинаковой скоростью сдвигать все его стороны параллельно самим себе внутрь многоугольника:

Пусть, для удобства, это движение происходит со скоростью 1 координатная единица в секунду (т.е. время в какомто смысле равно расстоянию: спустя единицу времени каждая точка преодолеет расстояние, равное единице).

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

Итак, нам надо научиться эффективно моделировать этот процесс сжатия. Для этого научимся для каждой стороны определять время, через которое она сожмётся в точку.

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

получаем формулу:

Теперь мы умеем за определять время, через которое сторона сожмётся в точку.

Занесём эти времена для каждой стороны в некую структуру данных для извлечения минимума, например, красно-чёрное дерево ( в языке C++).

Теперь если мы извлечём сторону с наименьшим временем , то эта сторона первой сожмётся в точку — в момент времени . Если многоугольник ещё не сжался в точку/отрезок, то эту сторону надо удалить из многоугольника, и продолжить алгоритм для оставшихся сторон. При удалении стороны мы должны соединить

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

При реализации для каждой стороны придётся хранить номер её правого и левого соседа (тем самым как бы построив двусвязный список из сторон многоугольника). Это позволяет реализовать удаление стороны и связывание двух её соседей за .

Если при удалении стороны оказывается, что её стороны-соседи параллельны, то это означает, что многоугольник после этого сжатия вырождается в точку/отрезок, поэтому мы можем сразу останавливать алгоритм и возвращать в качестве ответа время исчезнования текущей стороны (так что проблем с параллельными сторонами не возникает).

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

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

Реализация

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

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

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

const double EPS = 1E-9; const double PI = ...;

struct pt {

double x, y; pt() { }

pt (double x, double y) : x(x), y(y) { } pt operator- (const pt & p) const {

return pt (x-p.x, y-p.y);

}

};

double dist (const pt & a, const pt & b) {

return sqrt ((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));

}

double get_ang (const pt & a, const pt & b) {

double ang = abs (atan2 (a.y, a.x) - atan2 (b.y, b.x)); return min (ang, 2*PI-ang);

}

struct line {

double a, b, c;

line (const pt & p, const pt & q) { a = p.y - q.y;

b = q.x - p.x;

c = - a * p.x - b * p.y; double z = sqrt (a*a + b*b); a/=z, b/=z, c/=z;

}

};

double det (double a, double b, double c, double d) { return a * d - b * c;

}

pt intersect (const line & n, const line & m) { double zn = det (n.a, n.b, m.a, m.b); return pt (

-det (n.c, n.b, m.c, m.b) / zn,

-det (n.a, n.c, m.a, m.c) / zn

);

}

bool parallel (const line & n, const line & m) { return abs (det (n.a, n.b, m.a, m.b)) < EPS;

}

double get_h (const pt & p1, const pt & p2,

const pt & l1, const pt & l2, const pt & r1, const pt & r2)

{

pt q1 = intersect (line (p1, p2), line (l1, l2)); pt q2 = intersect (line (p1, p2), line (r1, r2)); double l = dist (q1, q2);

double alpha = get_ang (l2 - l1, p2 - p1) / 2; double beta = get_ang (r2 - r1, p1 - p2) / 2;

return l * sin(alpha) * sin(beta) / sin(alpha+beta);

}

struct cmp {

bool operator() (const pair<double,int> & a, const pair<double,int> & b) const {

if (abs (a.first - b.first) > EPS) return a.first < b.first;

return a.second < b.second;

}

};

int main() { int n;

vector<pt> p;

... чтение n и p ...

vector<int> next (n), prev (n); for (int i=0; i<n; ++i) {

next[i] = (i + 1) % n; prev[i] = (i - 1 + n) % n;

}

set < pair<double,int>, cmp > q; vector<double> h (n);

for (int i=0; i<n; ++i) { h[i] = get_h (

p[i], p[next[i]], p[i], p[prev[i]],

p[next[i]], p[next[next[i]]]

);

q.insert (make_pair (h[i], i));

}

double last_time; while (q.size() > 2) {

last_time = q.begin()->first; int i = q.begin()->second; q.erase (q.begin());

next[prev[i]] = next[i]; prev[next[i]] = prev[i];

int nxt = next[i], nxt1 = (nxt+1)%n, prv = prev[i], prv1 = (prv+1)%n;

if (parallel (line (p[nxt], p[nxt1]), line (p[prv], p[prv1]))) break;

q.erase (make_pair (h[nxt], nxt)); q.erase (make_pair (h[prv], prv));

h[nxt] = get_h ( p[nxt], p[nxt1], p[prv1], p[prv],

p[next[nxt]], p[(next[nxt]+1)%n]

);

h[prv] = get_h ( p[prv], p[prv1],

p[(prev[prv]+1)%n], p[prev[prv]],