test/geometry/contains.test.cpp

• View this file on GitHub
• Last update: 2023-05-08 16:51:58+09:00
• Problem: https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/all/CGL_3_C

Depends on

• geometry/line.hpp
• geometry/line_segment.hpp
• point (geometry/point.hpp)
• geometry/polygon.hpp

Code

#define PROBLEM \
    "https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/all/CGL_3_C"

#include <algorithm>
#include <iomanip>
#include <iostream>
#include <vector>

#include "geometry/line.hpp"
#include "geometry/line_segment.hpp"
#include "geometry/point.hpp"
#include "geometry/polygon.hpp"

namespace ebi {

void main_() {
    int n;
    std::cin >> n;
    Polygon poly(n);
    for (auto& [x, y] : poly) {
        std::cin >> x >> y;
    }
    int q;
    std::cin >> q;
    while (q--) {
        point p;
        std::cin >> p.x >> p.y;
        std::cout << contains(poly, p) << '\n';
    }
}

}  // namespace ebi

int main() {
    std::cout << std::fixed << std::setprecision(15);
    std::cin.tie(nullptr);
    std::ios::sync_with_stdio(false);
    ebi::main_();
}

#line 1 "test/geometry/contains.test.cpp"
#define PROBLEM \
    "https://onlinejudge.u-aizu.ac.jp/courses/library/4/CGL/all/CGL_3_C"

#include <algorithm>
#include <iomanip>
#include <iostream>
#include <vector>

#line 2 "geometry/line.hpp"

#include <cassert>
#include <cmath>

#line 2 "geometry/point.hpp"

#line 6 "geometry/point.hpp"
#include <cstdint>

#line 9 "geometry/point.hpp"

namespace ebi {

constexpr long double EPS = 1e-10;

const long double PI = std::acos(-1);

namespace internal {

int sgn(long double a) {
    return (a < -EPS) ? -1 : (a > EPS) ? 1 : 0;
}

long double add(long double a, long double b) {
    if (std::abs(a + b) < EPS * (std::abs(a) + std::abs(b))) return 0;
    return a + b;
}

}  // namespace internal


long double arg_to_radian(long double arg) {
    return PI * arg / (long double)(180);
}

struct point {
    long double x, y;

    point() = default;

    point(long double x, long double y) : x(x), y(y) {}

    point &operator+=(const point rhs) noexcept {
        x = internal::add(x, rhs.x);
        y = internal::add(y, rhs.y);
        return *this;
    }

    point &operator-=(const point rhs) noexcept {
        x = internal::add(x, -rhs.x);
        y = internal::add(y, -rhs.y);
        return *this;
    }

    point &operator*=(const point rhs) noexcept {
        long double _x = internal::add(x * rhs.x, -y * rhs.y);
        long double _y = internal::add(x * rhs.y, y * rhs.x);
        x = _x;
        y = _y;
        return *this;
    }

    point &operator*=(const long double k) noexcept {
        x *= k;
        y *= k;
        return *this;
    }

    point &operator/=(const long double k) {
        assert(internal::sgn(k) != 0);
        x /= k;
        y /= k;
        return *this;
    }

    point operator+(const point &rhs) const noexcept {
        return point(*this) += rhs;
    }

    point operator-(const point &rhs) const noexcept {
        return point(*this) -= rhs;
    }

    point operator*(const point &rhs) const noexcept {
        return point(*this) *= rhs;
    }

    point operator*(const long double rhs) const noexcept {
        return point(*this) *= rhs;
    }

    point operator/(const long double rhs) const {
        return point(*this) /= rhs;
    }

    point operator-() const noexcept {
        return point(0, 0) - *this;
    }

    long double abs() const noexcept {
        return std::sqrt(internal::add(x * x, y * y));
    }

    long double dot(const point rhs) const noexcept {
        return internal::add(x * rhs.x, y * rhs.y);
    }

    long double det(const point rhs) const noexcept {
        return internal::add(x * rhs.y, -y * rhs.x);
    }

    // arctan(y/x) (単位はラジアン)

    long double arg() const {
        return std::atan2(y, x);
    }

    // x昇順, その後y昇順

    bool operator<(const point &rhs) const noexcept {
        if (internal::sgn(x - rhs.x)) return internal::sgn(x - rhs.x) < 0;
        return internal::sgn(y - rhs.y) < 0;
    }

    bool operator==(const point &rhs) const noexcept {
        if (internal::sgn(x - rhs.x) == 0 && internal::sgn(y - rhs.y) == 0)
            return true;
        else
            return false;
    }

    bool operator!=(const point &rhs) const noexcept {
        return !((*this) == rhs);
    }
};

std::ostream &operator<<(std::ostream &os, const point &a) {
    return os << a.x << " " << a.y;
}

std::istream &operator>>(std::istream &os, point &a) {
    return os >> a.x >> a.y;
}

point conj(const point &a) {
    return point(a.x, -a.y);
}

// 点a をang(ラジアン)回転する

point rot(const point &a, long double ang) {
    return point(std::cos(ang) * a.x - std::sin(ang) * a.y,
                 std::sin(ang) * a.x + std::cos(ang) * a.y);
}

point rot90(const point &a) {
    return point(-a.y, a.x);
}

long double dot(const point &a, const point &b) {
    return a.dot(b);
}

long double det(const point &a, const point &b) {
    return a.det(b);
}

long double abs(const point &a) {
    return a.abs();
}

long double norm(const point &a) {
    return internal::add(a.x * a.x, a.y * a.y);
}

int isp(const point &a, const point &b, const point &c) {
    int flag = internal::sgn(det(b - a, c - a));
    if (flag == 0) {
        if (internal::sgn(dot(b - a, c - a)) < 0) return -2;
        if (internal::sgn(dot(a - b, c - b)) < 0) return +2;
    }
    return flag;
}

// 分割統治で最近点対を求める O(N log N)

long double closest_pair(std::vector<point> p) {
    std::sort(p.begin(), p.end());
    int n = p.size();
    auto f = [&](auto &&self, int l, int r) -> long double {
        if (r - l == 1) {
            return 1e9;
        }
        int mid = (l + r) / 2;
        long double x = p[mid].x;
        long double d = std::min(self(self, l, mid), self(self, mid, r));
        std::vector<point> b;
        b.reserve(r - l);
        int j = mid;
        for (int i = l; i < mid; i++) {
            while (j < r && p[j].y <= p[i].y) {
                b.emplace_back(p[j++]);
            }
            b.emplace_back(p[i]);
        }
        while (j < r) {
            b.emplace_back(p[j++]);
        }
        for (int i = 0; i < r - l; i++) {
            p[l + i] = b[i];
        }
        b.clear();
        for (int i = l; i < r; i++) {
            if (std::abs(p[i].x - x) >= d) continue;
            for (int j = int(b.size()) - 1; j >= 0; j--) {
                if (p[i].y - b[j].y >= d) break;
                d = std::min(d, abs(p[i] - b[j]));
            }
            b.emplace_back(p[i]);
        }
        return d;
    };
    return f(f, 0, n);
}

// ∠ABCを求める(ラジアン)

long double angle(const point &A, const point &B, const point &C) {
    long double a = (B - C).abs(), b = (C - A).abs(), c = (A - B).abs();
    long double cos =
        internal::add(internal::add(a * a, c * c), -b * b) / (2.0 * c * a);
    return std::acos(cos);
}

// 符号付き

long double calc_ang(point a, point b, point c) {
    long double cos = dot((a - b), (c - b)) / (abs(a - b) * abs(c - b));
    long double sin = det((a - b), (c - b)) / (abs(a - b) * abs(c - b));
    return std::atan2(sin, cos);
}

void arg_sort(std::vector<point> &a) {
    int n = a.size();
    std::vector ps(4, std::vector<point>());
    auto idx = [](point v) -> int {
        if (v.y >= 0)
            return (v.x >= 0) ? 0 : 1;
        else
            return (v.x >= 0) ? 3 : 2;
    };
    for (auto p : a) {
        assert(!(p.x == 0 && p.y == 0));
        ps[idx(p)].emplace_back(p);
    }
    a.clear();
    a.reserve(n);
    for (int i = 0; i < 4; i++) {
        std::sort(ps[i].begin(), ps[i].end(), [](point &p1, point &p2) -> bool {
            int flag = internal::sgn(internal::add(p1.x * p2.y, -p2.x * p1.y));
            return flag == 0 ? (norm(p1) < norm(p2)) : flag > 0;
        });
        for (auto &p : ps[i]) a.emplace_back(p);
    }
    return;
}

template <class T> void arg_sort_ll(std::vector<std::pair<T, T>> &a) {
    using Point = std::pair<T, T>;
    int n = a.size();
    std::vector ps(4, std::vector<Point>());
    auto idx = [](Point v) -> int {
        if (v.second >= 0)
            return (v.first >= 0) ? 0 : 1;
        else
            return (v.first >= 0) ? 3 : 2;
    };
    for (auto p : a) {
        assert(!(p.first == 0 && p.second == 0));
        ps[idx(p)].emplace_back(p);
    }
    a.clear();
    a.reserve(n);
    for (int i = 0; i < 4; i++) {
        std::sort(ps[i].begin(), ps[i].end(), [](Point &p1, Point &p2) -> bool {
            T flag = p1.first * p2.second - p2.first * p1.second;
            return flag == 0 ? (p1.first * p1.first + p1.second * p1.second <
                                p2.first * p2.first + p2.second * p2.second)
                             : flag > 0;
        });
        for (auto &p : ps[i]) a.emplace_back(p);
    }
    return;
}

}  // namespace ebi
#line 7 "geometry/line.hpp"

namespace ebi {

struct line {
    point a, b;

    line(long double x1, long double y1, long double x2, long double y2)
        : a(x1, y1), b(x2, y2) {}

    line(const point &a, const point &b) : a(a), b(b) {}

    point proj(const point &p) const {
        return a + (b - a) * (dot(b - a, p - a) / norm(b - a));
    }

    point relf(const point &p) const {
        return proj(p) * double(2) - p;
    }

    long double distance(const point &c) const {
        return std::abs(det(c - a, b - a) / abs(b - a));
    }
};

int intersection(const line &a, const line &b) {
    if (internal::sgn(det(a.b - a.a, b.a - b.b)) != 0) {
        if (internal::sgn(dot(a.b - a.a, b.b - b.a)) == 0) {  // 垂直
            return 1;
        }
        return 0;                                                // 交差
    } else if (internal::sgn(det(a.b - a.a, b.a - a.a)) != 0) {  // 平行
        return 2;
    } else {  // 同一直線
        return 3;
    }
}

point cross_point(const point &a, const point &b, const point &c,
                  const point &d) {
    return a + (b - a) * det(c - a, d - c) / det(b - a, d - c);
}

// 交点があるか確認する！
point cross_point(const line &s, const line &t) {
    assert(intersection(s, t) < 2);
    return s.a +
           (s.b - s.a) * det(t.a - s.a, t.b - t.a) / det(s.b - s.a, t.b - t.a);
}

// 直線aと点cの距離
long double distance(const line &a, const point &c) {
    return std::abs(det(c - a.a, a.b - a.a) / abs(a.b - a.a));
}

long double distance(const line &a, const line &b) {
    if (intersection(a, b) < 2) {
        return 0;
    } else {
        return distance(a, b.a);
    }
}

}  // namespace ebi
#line 2 "geometry/line_segment.hpp"

#line 5 "geometry/line_segment.hpp"

#line 8 "geometry/line_segment.hpp"

namespace ebi {

struct line_segment {
    point a, b;

    line_segment() = default;

    line_segment(long double x1, long double y1, long double x2, long double y2)
        : a(x1, y1), b(x2, y2) {}

    line_segment(const point &a, const point &b) : a(a), b(b) {}
};

// 線分ab, cd が交わるか判定
bool intersection_line_segment(const point &a, const point &b, const point &c,
                               const point &d) {
    if (internal::sgn(isp(a, b, c) * isp(a, b, d)) <= 0 &&
        internal::sgn(isp(c, d, a) * isp(c, d, b)) <= 0) {
        return true;
    }
    return false;
}

// 線分ab, cd が交わるか判定
bool intersection(const line_segment &a, const line_segment &b) {
    return intersection_line_segment(a.a, a.b, b.a, b.b);
}

bool intersection(const line &a, const line_segment &b) {
    if (internal::sgn(det(a.b - a.a, b.a - a.a)) *
            internal::sgn(det(a.b - a.a, b.b - a.a)) <
        0) {
        return true;
    } else {
        return false;
    }
}

point cross_point(const line_segment &s, const line_segment &t) {
    assert(intersection(s, t));
    return s.a +
           (s.b - s.a) * det(t.a - s.a, t.b - t.a) / det(s.b - s.a, t.b - t.a);
}

long double distance(const line_segment &a, const point &c) {
    if (internal::sgn(dot(a.b - a.a, c - a.a)) < 0) {
        return abs(c - a.a);
    } else if (internal::sgn(dot(a.a - a.b, c - a.b)) < 0) {
        return abs(c - a.b);
    } else {
        return std::abs(det(c - a.a, a.b - a.a) / abs(a.b - a.a));
    }
}

long double distance(const line_segment &a, const line_segment &b) {
    if (intersection(a, b)) {
        return 0;
    } else {
        return std::min(std::min(distance(a, b.a), distance(a, b.b)),
                        std::min(distance(b, a.a), distance(b, a.b)));
    }
}

long double distance(const line &a, const line_segment &b) {
    if (intersection(a, b)) {
        return 0;
    } else {
        return std::min(distance(a, b.a), distance(a, b.b));
    }
}

}  // namespace ebi
#line 2 "geometry/polygon.hpp"

#line 5 "geometry/polygon.hpp"

#line 8 "geometry/polygon.hpp"

namespace ebi {

using Polygon = std::vector<point>;

long double area(const Polygon &poly) {
    long double s = 0;
    int n = poly.size();
    for (int i = 0; i < n; i++) {
        s = internal::add(s, det(poly[i], poly[(i + 1 != n) ? i + 1 : 0]));
    }
    s /= 2.0;
    return s;
}

// 凸多角形か判定. pに反時計回りで点が入っていると仮定
bool is_convex(const Polygon &poly) {
    int n = poly.size();
    for (int i = 0; i < n; i++) {
        if (isp(poly[i], poly[(i + 1 != n) ? i + 1 : 0],
                poly[(i + 2 < n) ? i + 2 : (i + 2) % n]) == -1) {
            return false;
        }
    }
    return true;
}

enum { OUT, ON, IN };

int contains(const Polygon &poly, const point &p) {
    bool in = false;
    int n = poly.size();
    for (int i = 0; i < n; i++) {
        point a = poly[i] - p;
        point b = poly[(i + 1 < n) ? i + 1 : 0] - p;
        if (a.y > b.y) std::swap(a, b);
        if (a.y <= 0 && 0 < b.y && det(a, b) < 0) in = !in;
        if (det(a, b) == 0 && dot(a, b) <= 0) return ON;
    }
    return in ? IN : OUT;
}

long double convex_diameter(const Polygon &convex_poly) {
    int n = convex_poly.size();
    int is = 0, js = 0;
    for (int i = 1; i < n; i++) {
        if (convex_poly[i].y > convex_poly[is].y) is = i;
        if (convex_poly[i].y < convex_poly[js].y) js = i;
    }
    long double max = (convex_poly[is] - convex_poly[js]).abs();
    int i, max_i, j, max_j;
    i = max_i = is;
    j = max_j = js;
    do {
        if (det(convex_poly[(i + 1 < n) ? i + 1 : 0] - convex_poly[i],
                convex_poly[(j + 1 < n) ? j + 1 : 0] - convex_poly[j]) >= 0) {
            j = (j + 1 < n) ? j + 1 : 0;
        } else {
            i = (i + 1 < n) ? i + 1 : 0;
        }
        if ((convex_poly[i] - convex_poly[j]).abs() > max) {
            max = (convex_poly[i] - convex_poly[j]).abs();
            max_i = i;
            max_j = j;
        }
    } while (i != is || j != js);
    return max;
}

Polygon convex_polygon_cut(const Polygon &poly, const line &l) {
    Polygon ret;
    int n = poly.size();
    for (int i = 0; i < n; i++) {
        const point &now = poly[i];
        const point &nxt = poly[(i + 1 < n) ? i + 1 : 0];
        long double cf = det(l.a - now, l.b - now);
        long double cs = det(l.a - nxt, l.b - nxt);
        if (internal::sgn(cf) >= 0) {
            ret.emplace_back(now);
        }
        if (internal::sgn(cf) * internal::sgn(cs) < 0) {
            ret.emplace_back(cross_point(line(now, nxt), l));
        }
    }
    return ret;
}

}  // namespace ebi
#line 13 "test/geometry/contains.test.cpp"

namespace ebi {

void main_() {
    int n;
    std::cin >> n;
    Polygon poly(n);
    for (auto& [x, y] : poly) {
        std::cin >> x >> y;
    }
    int q;
    std::cin >> q;
    while (q--) {
        point p;
        std::cin >> p.x >> p.y;
        std::cout << contains(poly, p) << '\n';
    }
}

}  // namespace ebi

int main() {
    std::cout << std::fixed << std::setprecision(15);
    std::cin.tie(nullptr);
    std::ios::sync_with_stdio(false);
    ebi::main_();
}

Back to top page