MAXimal

добавлено: 10 Jun 2009 21:03
редактировано: 28 Oct 2015 22:35

Нахождение пары ближайших точек

Постановка задачи

Даны n точек p_i на плоскости, заданные своими координатами (x_i,y_i). Требуется найти среди них такие две точки, расстояние между которыми минимально:

 \min_{\scriptstyle i,j=0 \ldots n-1, \atop \scrip[...]

Расстояния мы берём обычные евклидовы:

 \rho(p_i,p_j) = \sqrt{ (x_i-x_j)^2 + (y_i-y_j)^2 [...]

Тривиальный алгоритм — перебор всех пар и вычисление расстояния для каждой — работает за O(n^2). Ниже описывается алгоритм, работающий за время O(n \log n). Этот алгоритм был предложен Препаратой (Preparata) в 1975 г. Препарата и Шамос также показали, что в модели дерева решений этот алгоритм асимптотически оптимален.

Алгоритм

Построим алгоритм по общей схеме алгоритмов "разделяй-и-властвуй": алгоритм оформляем в виде рекурсивной функции, которой передаётся множество точек; эта рекурсивная функция разбивает это множество пополам, вызывает себя рекурсивно от каждой половины, а затем выполняет какие-то операции по объединению ответов. Операция объединения заключается в обнаружении случаев, когда одна точка оптимального решения попала в одну половину, а другая точка — в другую (в этом случае рекурсивные вызовы от каждой из половинок отдельно обнаружить эту пару, конечно, не смогут). Основная сложность, как всегда, заключается в эффективной реализации этой стадии объединения. Если рекурсивной функции передаётся множество из n точек, то стадия объединения должна работать не более, чем O(n), тогда асимптотика всего алгоритма T(n) будет находиться из уравнения:

 T(n) = 2 T(n/2) + O(n).

Решением этого уравнения, как известно, является T(n) = O (n \log n).

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

 p_i < p_j \Longleftrightarrow (x_i < x_j) \lor \B[...]

Тогда возьмём среднюю после сортировки точку p_m (m = \lfloor n/2 \rfloor), и все точки до неё и саму p_m отнесём к первой половине, а все точки после неё — ко второй половине:

 A_1 = \{ p_i\ |\ i = 0 \ldots m \},
 A_2 = \{ p_i\ |\ i = m+1 \ldots n-1 \}.

Теперь, вызвавшись рекурсивно от каждого из множеств A_1 и A_2, мы найдём ответы h_1 и h_2 для каждой из половинок. Возьмём лучший из них: h = \min (h_1, h_2).

Теперь нам надо произвести стадию объединения, т.е. попытаться обнаружить такие пары точек, расстояние между которыми меньше h, причём одна точка лежит в A_1, а другая — в A_2. Очевидно, что для этого достаточно рассматривать только те точки, которые отстоят от вертикальной прямой раздела на расстояние, меньшее h, т.е. множество B рассматриваемых на этой стадии точек равно:

 B = \{ p_i\ |\ | x_i - x_m | < h \}.

Для каждой точки из множества B надо попытаться найти точки, находящиеся к ней ближе, чем h. Например, достаточно рассматривать только те точки, координата y которых отличается не более чем на h. Более того, не имеет смысла рассматривать те точки, у которых y-координата больше y-координаты текущей точки. Таким образом, для каждой точки p_i определим множество рассматриваемых точек C(p_i) следующим образом:

 C(p_i) = \{ p_j\ |\ p_j \in B,\ \ y_i - h < y_j \[...]

Если мы отсортируем точки множества B по y-координате, то находить C(p_i) будет очень легко: это несколько точек подряд до точки p_i.

Итак, в новых обозначениях стадия объединения выглядит следующим образом: построить множество B, отсортировать в нём точки по y-координате, затем для каждой точки p_i \in B рассмотреть все точки p_j \in C(p_i), и каждой пары (p_i,p_j) посчитать расстояние и сравнить с текущим наилучшим расстоянием.

На первый взгляд, это по-прежнему неоптимальный алгоритм: кажется, что размеры множеств C(p_i) будут порядка n, и требуемая асимптотика никак не получится. Однако, как это ни удивительно, можно доказать, что размер каждого из множеств C(p_i) есть величина O(1), т.е. не превосходит некоторой малой константы вне зависимости от самих точек. Доказательство этого факта приведено в следующем разделе.

Наконец, обратим внимание на сортировки, которых вышеописанный алгоритм содержит сразу две: сначала сортировка по парам (x,y), а затем сортировка элементов множества B по y. На самом деле, от обеих этих сортировок внутри рекурсивной функции можно избавиться (иначе бы мы не достигли оценки O(n) для стадии объединения, и общая асимптотика алгоритма получилась бы O(n \log^2 n)). От первой сортировки избавиться легко — достаточно предварительно, до запуска рекурсии, выполнить эту сортировку: ведь внутри рекурсии сами элементы не меняются, поэтому нет никакой необходимости выполнять сортировку заново. Со второй сортировкой чуть сложнее, выполнить её предварительно не получится. Зато, вспомнив сортировку слиянием (merge sort), которая тоже работает по принципу разделяй-и-властвуй, можно просто встроить эту сортировку в нашу рекурсию. Пусть рекурсия, принимая какое-то множество точек (как мы помним, упорядоченное по парам (x,y)) возвращает это же множество, но отсортированное уже по координате y. Для этого достаточно просто выполнить слияние (за O(n)) двух результатов, возвращённых рекурсивными вызовами. Тем самым получится отсортированное по y множество.

Оценка асимптотики

Чтобы показать, что вышеописанный алгоритм действительно выполняется за O(n \log n), нам осталось доказать следующий факт: |C(p_i)| = O(1).

Итак, пусть мы рассматриваем какую-то точку p_i; напомним, что множество C(p_i) — это множество точек, y-координата которых лежит в отрезке [y_i-h; y_i], а, кроме того, по координате x и сама точка p_i, и все точки множества C(p_i) лежат в полосе шириной 2h. Иными словами, рассматриваемые нами точки p_i и C(p_i) лежат в прямоугольнике размера 2h \times h.

Наша задача — оценить максимальное количество точек, которое может лежать в этом прямоугольнике 2h \times h; тем самым мы оценим и максимальный размер множества C(p_i). При этом при оценке надо не забывать, что могут встречаться повторяющиеся точки.

Вспомним, что h получалось как минимум из двух результатов рекурсивных вызовов — от множеств A_1 и A_2, причём A_1 содержит точки слева от линии раздела и частично на ней, A_2 — оставшиеся точки линии раздела и точки справа от неё. Для любой пары точек из A_1, равно как и из A_2, расстояние не может оказаться меньше h — иначе бы это означало некорректность работы рекурсивной функции.

Для оценки максимального количества точек в прямоугольнике 2h \times h разобьём его на два квадрата h \times h, к первому квадрату отнесём все точки C(p_i) \cap A_1, а ко второму — все остальные, т.е. C(p_i) \cap A_2. Из приведённых выше соображений следует, что в каждом из этих квадратов расстояние между любыми двумя точками не менее h.

Покажем, что в каждом квадрате не более четырёх точек. Например, это можно сделать следующим образом: разобьём квадрат на 4 подквадрата со сторонами h/2. Тогда в каждом из этих подквадратов не может быть больше одной точки (т.к. даже диагональ равна h / \sqrt{2}, что меньше h). Следовательно, во всём квадрате не может быть более 4 точек.

Итак, мы доказали, что в прямоугольнике 2h \times h не может быть больше 4 \cdot 2 = 8 точек, а, следовательно, размер множества C(p_i) не может превосходить 7, что и требовалось доказать.

Реализация

Введём структуру данных для хранения точки (её координаты и некий номер) и операторы сравнения, необходимые для двух видов сортировки:

struct pt {
	int x, y, id;
};
 
inline bool cmp_x (const pt & a, const pt & b) {
	return a.x < b.x || a.x == b.x && a.y < b.y;
}
 
inline bool cmp_y (const pt & a, const pt & b) {
	return a.y < b.y;
}
 
pt a[MAXN];

Для удобной реализации рекурсии введём вспомогательную функцию upd\_ans(), которая будет вычислять расстояние между двумя точками и проверять, не лучше ли это текущего ответа:

double mindist;
int ansa, ansb;
 
inline void upd_ans (const pt & a, const pt & b) {
	double dist = sqrt ((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y) + .0);
	if (dist < mindist)
		mindist = dist,  ansa = a.id,  ansb = b.id;
}

Наконец, реализация самой рекурсии. Предполагается, что перед её вызовом массив a[] уже отсортирован по x-координате. Рекурсии передаётся просто два указателя l, r, которые указывают, что она должна искать ответ для a[l \ldots r]. Если расстояние между r и l слишком мало, то рекурсию надо остановить, и выполнить тривиальный алгоритм поиска ближайшей пары и затем отсортировать подмассив по y-координате.

Для слияния двух множеств точек, полученных от рекурсивных вызовов, в одно (упорядоченное по y-координате), мы используем стандартную функцию STL merge(), и создаём вспомогательный буфер t[] (один на все рекурсивные вызовы). (Использовать inplace\_merge() нецелесообразно, т.к. она в общем случае работает не за линейное время).

Наконец, множество B хранится в том же массиве t.

void rec (int l, int r) {
	if (r - l <= 3) {
		for (int i=l; i<=r; ++i)
			for (int j=i+1; j<=r; ++j)
				upd_ans (a[i], a[j]);
		sort (a+l, a+r+1, &cmp_y);
		return;
	}
 
	int m = (l + r) >> 1;
	int midx = a[m].x;
	rec (l, m),  rec (m+1, r);
	static pt t[MAXN];
	merge (a+l, a+m+1, a+m+1, a+r+1, t, &cmp_y);
	copy (t, t+r-l+1, a+l);
 
	int tsz = 0;
	for (int i=l; i<=r; ++i)
		if (abs (a[i].x - midx) < mindist) {
			for (int j=tsz-1; j>=0 && a[i].y - t[j].y < mindist; --j)
				upd_ans (a[i], t[j]);
			t[tsz++] = a[i];
		}
}

Кстати говоря, если все координаты целые, то на время работы рекурсии можно вообще не переходить к дробным величинам, и хранить в mindist квадрат минимального расстояния.

В основной программе вызывать рекурсию следует так:

sort (a, a+n, &cmp_x);
mindist = 1E20;
rec (0, n-1);

Обобщение: поиск треугольника с минимальным периметром

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

По сути, для решения этой задачи алгоритм остаётся прежним: мы разделяем поле на две половинки вертикальной прямой, вызываем решение рекурсивно от обеих половинок, выбираем минимум minper из найденных периметров, строим полоску толщиной minper / 2, и в ней перебираем все треугольники, способные улучшить ответ. (Отметим, что у треугольника с периметром \le minper длиннейшая сторона \le minper/2.)

Задачи в online judges

Список задач, которые сводятся к поиску двух ближайших точек: