e-maxx_algo

Центры тяжести многоугольников и многогранников

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

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

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

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

Двумерный случай: многоугольники

На самом деле, говоря о центре масс двумерной фигуры, можно иметь в виду одну из трёх следующих задач:

●Центр масс системы точек — т.е. вся масса сосредоточена только в вершинах многоугольника.

●Центр масс каркаса — т.е. масса многоугольника сосредоточена на его периметре.

●Центр масс сплошной фигуры — т.е. масса многоугольника распределена по всей его площади.

Каждая из этих задач имеет самостоятельное решение, и будет рассмотрена ниже отдельно.

Центр масс системы точек

Это самая простая из трёх задач, и её решение — известная физическая формула центра масс системы материальных точек:

— искомый радиус-вектор центра масс.

В частности, если все точки имеют одинаковую массу, то координаты центра масс есть

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

и совпадает с точкой пересечения медиан:

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

и, выражая отсюда , мы и получаем требуемую формулу.

Центр масс каркаса

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

Но тогда каждую сторону многоугольника можно заменить одной точкой — серединой этого отрезка (т.к. центр масс однородного отрезка есть середина этого отрезка), с массой, равной длине этого отрезка.

Теперь мы получили задачу о системе материальных точек, и применяя к ней решение из предыдущего пункта, мы находим:

Наша цель — показать, что вектора и коллинеарны.

Обозначим через и точки, являющиеся центрами масс треугольников №3 и №4. Тогда, очевидно, центром масс совокупности этих двух треугольников будет точка , являющаяся серединой отрезка . Более того, вектор от точки к точке совпадает с вектором .

Искомый центр масс треугольника лежит посередине отрезка, соединяющего точки и (поскольку мы разбили треугольник на две части равных площадей: №1-№2 и №3-№4):

подобен треугольнику с коэффициентом , то этот же вектор равен . Отсюда получаем уравнение:

откуда находим:

Более того, попутно мы доказали, что центроид делит каждую медиану в отношении , считая от вершины.

Случай многоугольника

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

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

Окончательная формула получается следующей:

где — центроид -го треугольника в триангуляции заданного многоугольника, — площадь -го треугольника триангуляции, — площадь всего многоугольника.

Триангуляция выпуклого многоугольника — тривиальная задача: для этого, например, можно взять треугольники , где .

Случай многоугольника: альтернативный способ

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

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

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

Случай произвольного многогранника

Перейдём теперь к общему случаю — случаю произвольного многогранника.

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

Пересечение окружности и прямой

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

Решение

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

Предположим, не теряя общности, что центр окружности находится в начале координат (если это не так, то перенесём его туда, исправив соответствующе константу C в уравнении прямой). Т.е. имеем окружность с центром в (0,0) радиуса r и прямую с уравнением Ax + By + C = 0.

Сначала найдём ближайшую к центру точку прямой - точку с некоторыми координатами (x0,y0). Вопервых, эта точка должна находиться на таком расстоянии от начала координат:

|C|

----------

sqrt(A2+B2)

Во-вторых, поскольку вектор (A,B) перпендикулярен прямой, то координаты этой точки должны быть пропорциональны координатам этого вектора. Учитывая, что расстояние от начала координат до искомой точки нам известно, нам нужно просто нормировать вектор (A,B) к этой длине, и мы получаем:

A C x0 = - -----

A2+B2

B C y0 = - -----

A2+B2

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

Зная ближайшую к центру окружности точку, мы уже можем определить, сколько точек будет содержать ответ, и даже дать ответ, если этих точек 0 или 1.

Действительно, если расстояние от (x0, y0) до начала координат (а его мы уже выразили формулой - см. выше)

больше радиуса, то ответ - ноль точек. Если это расстояние равно радиусу, то ответом будет одна точка - (x0,y0). А вот в оставшемся случае точек будет две, и их координаты нам предстоит найти.

Итак, мы знаем, что точка (x0, y0) лежит внутри круга. Искомые точки (ax,ay) и (bx,by), помимо того что должны принадлежать прямой, должны лежать на одном и том же расстоянии d от точки (x0, y0), причём это расстояние легко найти:

C2

d = sqrt ( r2 - ----- ) A2+B2

Заметим, что вектор (-B,A) коллинеарен прямой, а потому искомые точки (ax,ay) и (bx,by) можно получить, прибавив к точке (x0,y0) вектор (-B,A), нормированный к длине d (мы получим одну искомую точку), и вычтя этот же вектор

(получим вторую искомую точку). Окончательное решение такое:

d2 mult = sqrt ( ----- )

A2+B2

ax = x0 + B mult

ay = y0 - A mult bx = x0 - B mult by = y0 + A mult

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

Реализация

Как и было указано в начале описания, предполагается, что окружность расположена в начале координат. Поэтому входные параметры - это радиус окружности и коэффициенты A,B,C уравнения прямой.

double r, a, b, c; // входные данные

double x0 = -a*c/(a*a+b*b), y0 = -b*c/(a*a+b*b); if (c*c > r*r*(a*a+b*b)+EPS)

puts ("no points");

else if (abs (c*c - r*r*(a*a+b*b)) < EPS) { puts ("1 point");

cout << x0 << ' ' << y0 << '\n';

}

else {

double d = r*r - c*c/(a*a+b*b); double mult = sqrt (d / (a*a+b*b)); double ax,ay,bx,by;

ax = x0 + b * mult; bx = x0 - b * mult; ay = y0 - a * mult; by = y0 + a * mult; puts ("2 points");

cout << ax << ' ' << ay << '\n' << bx << ' ' << by << '\n';

}

Пересечение двух окружностей

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

Решение

Сведём нашу задачу к задаче о Пересечении окружности и прямой.

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

x2 + y2 = r12

(x - x2)2 + (y - y2)2 = r22

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

x2 + y2 = r12

x (-2x2) + y (-2y2) + (x22 + y22 + r12 - r22) = 0

Таким образом, мы свели задачу о пересечении двух окружностей к задаче о пересечении первой окружности

иследующей прямой:

Ax + By + C = 0, A = -2x2,

B = -2y2,

C = x22 + y22 + r12 - r22.

А решение последней задачи описано в соответствующей статье.

Единственный вырожденный случай, который надо рассмотреть отдельно - когда центры окружностей совпадают. Действительно, в этом случае вместо уравнения прямой мы получим уравнение вида 0 = С, где C - некоторое число, и этот случай будет обрабатываться некорректно. Поэтому этот случай нужно рассмотреть отдельно: если радиусы окружностей совпадают, то ответ - бесконечность, иначе - точек пересечения нет.

Построение выпуклой оболочки обходом Грэхэма

Даны N точек на плоскости. Построить их выпуклую оболочку, т.е. наименьший выпуклый многоугольник, содержащий все эти точки.

Мы рассмотрим метод Грэхэма (Graham) (предложен в 1972 г.) с улучшениями Эндрю (Andrew) (1979 г.). С его помощью можно построить выпуклую оболочку за время O (N log N) с использованием только операций сравнения, сложения и умножения. Алгоритм является асимптотически оптимальным (доказано, что не существует алгоритма с лучшей асимптотикой), хотя в некоторых задачах он неприемлем (в случае параллельной обработки или при online-обработке).

Описание

Алгоритм. Найдём самую левую и самую правую точки A и B (если таких точек несколько, то возьмём самую нижнюю среди левых, и самую верхнюю среди правых). Понятно, что и A, и B обязательно попадут в выпуклую оболочку. Далее, проведём через них прямую AB, разделив множество всех точек на верхнее и нижнее подмножества S1 и S2 (точки, лежащие на прямой, можно отнести к любому множеству - они всё равно не войдут в оболочку). Точки A и

B отнесём к обоим множествам. Теперь построим для S1 верхнюю оболочку, а для S2 - нижнюю оболочку, и объединим их, получив ответ. Чтобы получить, скажем, верхнюю оболочку, нужно отсортировать все точки по абсциссе, затем пройтись по всем точкам, рассматривая на каждом шаге кроме самой точки две предыдущие точки, вошедшие в оболочку. Если текущая тройка точек образует не правый поворот (что легко проверить с помощью Ориентированной площади), то ближайшего соседа нужно удалить из оболочки. В конце концов, останутся только

точки, входящие в выпуклую оболочку.

Итак, алгоритм заключается в сортировке всех точек по абсциссе и двух (в худшем случае) обходах всех точек, т. е. требуемая асимптотика O (N log N) достигнута.

Реализация

struct pt {

double x, y;

};

bool cmp (pt a, pt b) {

return a.x < b.x || a.x == b.x && a.y < b.y;

}

bool cw (pt a, pt b, pt c) {

return a.x*(b.y-c.y)+b.x*(c.y-a.y)+c.x*(a.y-b.y) < 0;

}

bool ccw (pt a, pt b, pt c) {

return a.x*(b.y-c.y)+b.x*(c.y-a.y)+c.x*(a.y-b.y) > 0;

}

void convex_hull (vector<pt> & a) { if (a.size() == 1) return;

sort (a.begin(), a.end(), &cmp); pt p1 = a[0], p2 = a.back(); vector<pt> up, down; up.push_back (p1); down.push_back (p1);

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

if (i==a.size()-1 || cw (p1, a[i], p2)) {

while (up.size()>=2 && !cw (up[up.size()-2], up[up.

size()-1], a[i]))

up.pop_back(); up.push_back (a[i]);

}

if (i==a.size()-1 || ccw (p1, a[i], p2)) {

while (down.size()>=2 && !ccw (down[down.size()-2], down[down.size()-1], a[i]))

down.pop_back(); down.push_back (a[i]);

}

}

a.clear();

for (size_t i=0; i<up.size(); ++i) a.push_back (up[i]);

for (size_t i=down.size()-2; i>0; --i) a.push_back (down[i]);

}