Минимальная окружность, покрывающая множество точек

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

    1. Пусть наша окружность обязательно должна проходить через заданные три точки, тогда мы строим описывающую окружность около этого заданного треугольника.
    2. Пусть наша окружность обязательно должна проходить через две заданные точки, тогда мы перебираем все точки, берём одну из точек за третью и строим окружность согласно пункту 1. На каждом шаге перебора мы храним найденную окружность. Если текущая точка не принадлежит найденной, то найденная окружность будет равна текущей окружности, которую мы построили.
    3. Пусть наша окружность обязательно должна проходить через одну заданную точку, тогда мы аналогично пункту 2 перебираем все точки, строим текущую окружность согласно пункту 2 и, если текущая точка не принадлежит найденной окружности, приравниваем найденной текущую окружность.
    4. Пусть нам надо найти минимальную окружность, покрывающую множество точек, тогда мы перебираем все точки, строим окружности согласно пункту 3, и аналогично пунктам 2 и 3 находим искомую окружность.

    Асимптотика O(N*logN).

Листинг C++

const double inf = 1e17;
const double eps = 1e-7;

#define sqr(a) ((a) * (a))
#define dist2(a, b) (sqr (a.x - b.x) + sqr (a.y - b.y))
#define dist(a, b) (sqrt ((double) dist2 (a, b)))
#define abs(a) ((a) < 0 ?  - (a) : (a))
#define equald(a, b) (abs (a - b) <= eps)
#define lessd(a, b) (!equald (a, b) && a < b)
#define min(a, b) (lessd (a, b) ? a : b)
#define max(a, b) (lessd (a, b) ? b : a)
#define lesspxy(a, b) (lessd (a.x, b.x) || (equald (a.x, b.x) && lessd (a.y, b.y)))
#define minp(a, b) (lesspxy (a, b) ? a : b)
#define maxp(a, b) (lesspxy (a, b) ? b : a)
#define area_triangle(a, b, c) (0.5 * (a.x * b.y + b.x * c.y + c.x * a.y - a.y * b.x - b.y * c.x - c.y * a.x))
#define ccw(a, b, c) (lessd (eps, area_triangle (a, b, c)))
#define on_line(a, b, c) (equald (area_triangle (a, b, c), eps))
#define part(a, b, m, n) (point ((a.x * m + b.x * n) / (double)(m + n), (a.y * m + b.y * n) / (double)(m + n)))

class point
{
public:
    double x, y;
    point (double _x = 0, double _y = 0) : x (_x), y (_y) {};
};

void min_disk3 (point a, point b, point c, point &circ, double &r)
{
    if (on_line (a, b, c))
    {
        point p1 = minp (minp (a, b), c);
        point p2 = maxp (maxp (a, b), c);

        circ = part (p1, p2, 1, 1);
        r = 0.5 * dist (p1, p2);

        return;
    }

    point mab = part (a, b, 1, 1);
    point mac = part (a, c, 1, 1);
    double a1 = a.x - b.x;
    double b1 = a.y - b.y;
    double c1 = - a1 * mab.x - b1 * mab.y;
    double a2 = a.x - c.x;
    double b2 = a.y - c.y;
    double c2 = - a2 * mac.x - b2 * mac.y;

    double x, y;
    x = (b2 * c1 - b1 * c2) / (a2 * b1 - a1 * b2);
    if (!equald (b1, eps))
        y = - (a1 * x + c1) / b1;
    else
        y = - (a2 * x + c2) / b2;

    circ = point (x, y);
    r = dist (circ, a);
}

void min_disk2 (point p[], int n, point a, point b, point &circ, double &r)
{
    circ = part (a, b, 1, 1);
    r = dist (circ, a);
    for (int i = 0; i < n; ++ i)
        if (lessd (r, dist (p[i], circ)))
            min_disk3 (p[i], a, b, circ, r);
}
void min_disk1 (point p[], int n, point a, point &circ, double &r)
{
    circ = part (p[0], a, 1, 1);
    r = dist (circ, a);
    for (int i = 1; i < n; ++ i)
        if (lessd (r, dist (p[i], circ)))
            min_disk2 (p, i, a, p[i], circ, r);
}
void min_disk (point p[], int n, point &circ, double &r)
{
    circ = part (p[0], p[1], 1, 1);
    r = dist (circ, p[0]);
    for (int i = 2; i < n; ++ i)
        if (lessd (r, dist (p[i], circ)))
            min_disk1 (p, i, p[i], circ, r);
}

16.07.2009
14:52

По всем вопросам обращаться: rumterg@gmail.com