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

e-maxx algorithm оригинал: C/C++ #algorithm #emaxx #geometry #math

Static

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

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

} }

C# решение

auto-draft, проверить перед отправкой

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;
    };
    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 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");
            }
    }
}

C++ решение

сопоставлено/оригинал

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");
        }
}

Java решение

auto-draft, проверить перед отправкой

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;
    };
    struct ang {
            int a, b;
    };
    boolean 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 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; // очередной запрос
                    boolean 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");
            }
    }
}

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

Вакансии для этой задачи

Показаны активные вакансии с пересечением по тегам задачи.

Все вакансии

Активных вакансий пока нет.