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

Дан выпуклый многоугольник с 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");
	}

}