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

Источник: e-maxx.ru/algo, страница PDF 269.

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

given

точек

на плоскости, заданные своими координатами

. it is required find среди них такие две

точки, расстояние между которыми минимально: Расстояния мы берём обычные евклидовы: Тривиальный Thuật toán — перебор всех пар и вычисление расстояния для каждой — работает за .

Ниже описывается Thuật toán, работающий за время

. Этот Thuật toán был предложен Препаратой

(Preparata) в 1975 г. Препарата и Шамос также показали, что в модели дерева решений этот Thuật toán асимптотически оптимален.

Thuật toán

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

точек,

то стадия объединения должна работать не более, чем

, тогда Asymptotic complexity всего Thuật toánа

будет находиться из уравнения:

Lời giảiм этого уравнения, как известно, является

. Итак, перейдём к построению Thuật toánа. Чтобы в будущем прийти к эффективной реализации стадии

объединения, разбивать множество точек на два будем согласно их

-координатам: фактически мы проводим

некоторую вертикальную прямую, разбивающую множество точек на два подмножества Ví dụно одинаковых размеров. Такое разбиение удобно произвести следующим образом: отсортируем точки стандартно как пары чисел, т.е.:

Тогда возьмём среднюю после сортировки точку

(

), и все точки до неё и саму

отнесём к

первой половине, а все точки после неё — ко второй половине:

Теперь, вызвавшись рекурсивно от каждого из множеств

и

, мы найдём ответы

и

для каждой из

половинок. Возьмём лучший из них: . Теперь нам надо произвести стадию объединения, т.е. попытаться обнаружить такие пары точек,

расстояние между которыми меньше

, причём одна точка лежит в

, а другая — в

. Очевидно, что для

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

расстояние, меньшее

, т.е. множество

рассматриваемых на этой стадии точек равно:

Для каждой точки из множества

надо попытаться find точки, находящиеся к ней ближе, чем

. НаVí dụ,

достаточно рассматривать только те точки, координата

которых отличается не более чем на

. Более того, не

имеет смысла рассматривать те точки, у которых

-координата больше

-координаты текущей точки. Таким образом,

для каждой точки

определим множество рассматриваемых точек

следующим образом:

Если мы отсортируем точки множества

по

-координате, то находить

будет очень легко: это несколько

точек подряд до точки

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

, отсортировать в нём точки по

-координате, затем для каждой точки

рассмотреть все точки

, и каждой пары

посчитать расстояние и сравнить с текущим наилучшим расстоянием. На первый взгляд, это по-прежнему неоптимальный Thuật toán: кажется, что размеры множеств

будут порядка

,

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

из множеств

есть величина

, т.е. не превосходит некоторой малой константы вне зависимости от самих точек. Chứng minh этого факта приведено в следующем разделе. Наконец, обратим внимание на сортировки, которых вышеописанный Thuật toán содержит сразу две: сначала сортировка

по парам (

,

), а затем сортировка elementов множества

по

. На самом деле, от обеих этих сортировок

внутри рекурсивной функции можно избавиться (иначе бы мы не достигли оценки

для стадии объединения,

и общая Asymptotic complexity Thuật toánа получилась бы

). От первой сортировки избавиться легко —

достаточно предварительно, до запуска рекурсии, выполнить эту сортировку: ведь внутри рекурсии сами elementы не меняются, поэтому нет никакой необходимости выполнять сортировку заново. Со второй сортировкой чуть сложнее, выполнить её предварительно не получится. Зато, вспомнив сортировку слиянием (merge sort), которая тоже работает по принципу разделяй-и-властвуй, можно просто встроить эту сортировку в нашу рекурсию. Пусть рекурсия, принимая какое-то множество точек (как мы помним, упорядоченное по парам ) returns это же множество, но отсортированное уже по координате

. Для этого достаточно просто

выполнить слияние (за

) двух результатов, возвращённых рекурсивными вызовами. Тем самым

получится отсортированное по

множество.

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

Чтобы показать, что вышеописанный Thuật toán действительно выполняется за

, нам осталось

доказать следующий факт: .

Итак, пусть мы рассматриваем какую-то точку

; напомним, что множество

— это множество точек,

-

координата которых не больше

, но и не меньше

, а, кроме того, по координате

и сама точка

, и все

точки множества

лежат в полосе шириной

. Иными словами, рассматриваемые нами точки

и

лежат в прямоугольнике размера

. Наша Bài toán — оценить максимальное количество точек, которое может лежать в этом прямоугольнике

;

тем самым мы оценим и maximum размер множества

(он будет на единицу меньше, т.к. в

этом прямоугольнике лежит ещё и точка

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

Вспомним, что

получалось как минимум из двух результатов рекурсивных вызовов — от множеств

и

, причём

содержит точки слева от линии раздела и частично на ней,

— оставшиеся точки линии раздела и точки справа

от неё. Для любой пары точек из

, равно как и из

, расстояние не может оказаться меньше

— иначе бы

это означало некорректность работы рекурсивной функции.

Для оценки максимального количества точек в прямоугольнике

разобьём его на два квадрата

, к

первому квадрату отнесём все точки

, а ко второму — все остальные, т.е.

. Из

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

не превосходит

. Но тогда это означает, что в каждом квадрате не более четырёх точек!

Действительно, пусть есть квадрат

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

. Докажем, что в квадрате не может быть больше 4 точек. НаVí dụ, это можно сделать следующим образом:

разобьём этот квадрат на 4 квадрата со сторонами

. Тогда в каждом из этих маленьких квадратов не может

быть больше одной точки (т.к. даже диагональ равна

, что меньше

). Следовательно, во всём квадрате никак

не может быть более 4 точек.

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

не может быть больше

точек, а, следовательно,

размер множества

не может превосходить

, что и требовалось доказать.

Cài đặt

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

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];

Для удобной реализации рекурсии введём вспомогательную функцию

, которая будет вычислять

расстояние между двумя точками и проверять, не лучше ли это текущего ответа:

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;

} Наконец, Cài đặt самой рекурсии. Предполагается, что перед её вызовом mảng

уже отсортирован по

-координате. Рекурсии передаётся просто два указателя

,

, которые указывают, что она должна искать ответ

для

. Если расстояние между

и

слишком мало, то рекурсию надо остановить, и выполнить

тривиальный Thuật toán поиска ближайшей пары и затем отсортировать подmảng по -координате. Для слияния двух множеств точек, полученных от рекурсивных вызовов, в одно (упорядоченное по

-координате),

мы используем стандартную функцию STL

, и создаём вспомогательный буфер

(один на все

рекурсивные вызовы). (Использовать

нецелесообразно, т.к. она в общем случае работает не

за линейное время).

Наконец, множество

хранится в том же mảngе

.

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];

} } Кстати говоря, если все координаты целые, то на Thời gian chạy рекурсии можно вообще не переходить к

дробным величинам, и хранить в

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

sort (a, a+n, &cmp_x);

mindist = 1E20;

rec (0, n-1);

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

минимальным периметром

Описанный выше Thuật toán интересно обобщается и на эту задачу: среди заданного множества точек выбрать три различных точки так, чтобы сумма попарных расстояний между ними была наименьшей. По сути, для решения этой задачи Thuật toán остаётся прежним: мы разделяем поле на две половинки вертикальной прямой, вызываем Lời giải рекурсивно от обеих половинок, выбираем минимум

из

найденных периметров, строим полоску толщиной

, и в ней перебираем все треугольники,

могущие улучшить ответ. (Отметим, что у треугольника с периметром

длиннейшая

сторона

.)

Задачи в online judges

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

● UVA 10245 "The Closest Pair Problem" [Complexity: низкая]

● SPOJ #8725 CLOPPAIR "Closest Point Pair" [Complexity: низкая]

● CODEFORCES Командная олимпиада школьников Саратова - 2011 "Минимальная сумма"

[Complexity: средняя]

● Google CodeJam 2009 Final "Min Perimeter" [Complexity: средняя]

● SPOJ #7029 CLOSEST "Closest Triple" [Complexity: средняя]

using System;
using System.Collections.Generic;
using System.Linq;

public static class AlgorithmDraft
{
    // Auto-generated C# draft from the original e-maxx C/C++ listing. Review before production use.
    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];
    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;
    }
    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];
                    }
    }
    sort (a, a+n, &cmp_x);
    mindist = 1E20;
    rec (0, n-1);
}

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];
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;
}
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];
                }
}
sort (a, a+n, &cmp_x);
mindist = 1E20;
rec (0, n-1);

import java.util.*;
import java.math.*;

public class AlgorithmDraft {
    // Auto-generated Java draft from the original e-maxx C/C++ listing. Review before production use.
    struct pt {
            int x, y, id;
    };
    inline boolean cmp_x (const pt & a, const pt & b) {
            return a.x < b.x || a.x == b.x && a.y < b.y;
    }
    inline boolean cmp_y (const pt & a, const pt & b) {
            return a.y < b.y;
    }
    pt a[MAXN];
    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;
    }
    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];
                    }
    }
    sort (a, a+n, &cmp_x);
    mindist = 1E20;
    rec (0, n-1);
}