geometry/convex_polygon_diameter.hpp

• View this file on GitHub
• Last update: 2023-06-01 23:47:10+09:00
• Include: #include "geometry/convex_polygon_diameter.hpp"
• Link: http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_B
• Link: https://en.wikipedia.org/wiki/Rotating_calipers)

Depends on

• geometry/ccw.hpp
• geometry/point.hpp
• geometry/polygon.hpp
• geometry/polygon_is_convex.hpp

Required by

• geometry/all.hpp
• 最遠点対 (geometry/farthest_pair.hpp)

Verified with

• verify/aoj_cgl/aoj_cgl_4_b.test.cpp

Code

#pragma once

#include "geometry/polygon.hpp"
#include "geometry/polygon_is_convex.hpp"

// convex polygon diameter
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_B
// return {index1, index2, diameter}
// using the method of rotating calipers (https://en.wikipedia.org/wiki/Rotating_calipers)
// complexity: O(n)
template <typename T> std::tuple<int, int, T> convex_polygon_diameter(const Polygon<T> &p) {
    assert(polygon_is_convex(p));
    int n = (int)p.size();
    assert(n >= 2);
    if (n == 2) {
        return {0, 1, abs(p[0] - p[1])};
    }
    auto [it_min, it_max] = std::minmax_element(p.begin(), p.end(), compare_x<T>);
    int idx_min = it_min - p.begin();
    int idx_max = it_max - p.begin();

    T maxdis = norm(p[idx_max] - p[idx_min]);
    int maxi = idx_min, i = idx_min, maxj = idx_max, j = idx_max;
    do {
        int ni = (i + 1 == n ? 0 : i + 1), nj = (j + 1 == n ? 0 : j + 1);
        if (sign(cross(p[ni] - p[i], p[nj] - p[j])) < 0) {
            i = ni;
        } else {
            j = nj;
        }
        if (norm(p[i] - p[j]) > maxdis) {
            maxdis = norm(p[i] - p[j]);
            maxi = i;
            maxj = j;
        }
    } while (i != idx_min || j != idx_max);
    return {maxi, maxj, abs(p[maxi] - p[maxj])};
}

#line 2 "geometry/convex_polygon_diameter.hpp"

#line 2 "geometry/polygon.hpp"

#line 2 "geometry/point.hpp"

// point
template <typename T> struct Point {
    static T EPS;
    static constexpr T PI = 3.1415926535'8979323846'2643383279L;
    static void set_eps(const T &e) { EPS = e; }
    T x, y;
    Point(const T x = T(0), const T y = T(0)) : x(x), y(y) {}
    Point &operator+=(const Point &p) {
        x += p.x;
        y += p.y;
        return *this;
    }
    Point &operator-=(const Point &p) {
        x -= p.x;
        y -= p.y;
        return *this;
    }
    Point &operator*=(const Point &p) { return *this = Point(x * p.x - y * p.y, x * p.y + y * p.x); }
    Point &operator*=(const T &k) {
        x *= k;
        y *= k;
        return *this;
    }
    Point &operator/=(const Point &p) { return *this = Point(x * p.x + y * p.y, -x * p.y + y * p.x) / (p.x * p.x + p.y * p.y); }
    Point &operator/=(const T &k) {
        x /= k;
        y /= k;
        return *this;
    }

    Point operator+() const { return *this; }
    Point operator-() const { return Point(-x, -y); }

    friend Point operator+(const Point &a, const Point &b) { return Point(a) += b; }
    friend Point operator-(const Point &a, const Point &b) { return Point(a) -= b; }
    friend Point operator*(const Point &a, const Point &b) { return Point(a) *= b; }
    friend Point operator*(const Point &p, const T &k) { return Point(p) *= k; }
    friend Point operator/(const Point &a, const Point &b) { return Point(a) /= b; }
    friend Point operator/(const Point &p, const T &k) { return Point(p) /= k; }
    // for std::set, std::map, compare_arg, ...
    friend bool operator<(const Point &a, const Point &b) { return a.x == b.x ? a.y < b.y : a.x < b.x; }
    friend bool operator>(const Point &a, const Point &b) { return a.x == b.x ? a.y > b.y : a.x > b.x; }
    // I/O
    friend std::istream &operator>>(std::istream &is, Point &p) { return is >> p.x >> p.y; }
    friend std::ostream &operator<<(std::ostream &os, const Point &p) { return os << '(' << p.x << ' ' << p.y << ')'; }
};

// template
template <typename T> inline int sign(const T &x) { return x < -Point<T>::EPS ? -1 : (x > Point<T>::EPS ? 1 : 0); }
template <typename T> inline bool equal(const T &a, const T &b) { return sign(a - b) == 0; }
template <typename T> inline T radian_to_degree(const T &r) { return r * 180.0 / Point<T>::PI; }
template <typename T> inline T degree_to_radian(const T &d) { return d * Point<T>::PI / 180.0; }

// contain enum
constexpr int IN = 2;
constexpr int ON = 1;
constexpr int OUT = 0;

// equal (point and point)
template <typename T> inline bool equal(const Point<T> &a, const Point<T> &b) { return equal(a.x, b.x) and equal(a.y, b.y); }
// inner product
template <typename T> inline T dot(const Point<T> &a, const Point<T> &b) { return a.x * b.x + a.y * b.y; }
// outer product
template <typename T> inline T cross(const Point<T> &a, const Point<T> &b) { return a.x * b.y - a.y * b.x; }
// rotate Point p counterclockwise by theta radian
template <typename T> inline Point<T> rotate(const Point<T> &p, const T &theta) { return p * Point<T>(std::cos(theta), std::sin(theta)); }
// compare (x, y)
template <typename T> inline bool compare_x(const Point<T> &a, const Point<T> &b) { return equal(a.x, b.x) ? sign(a.y - b.y) < 0 : sign(a.x - b.x) < 0; }
// compare (y, x)
template <typename T> inline bool compare_y(const Point<T> &a, const Point<T> &b) { return equal(a.y, b.y) ? sign(a.x - b.x) < 0 : sign(a.y - b.y) < 0; }
// compare by (arg(p), norm(p)) [0, 360)
template <typename T> inline bool compare_arg(const Point<T> &a, const Point<T> &b) {
    // https://ngtkana.hatenablog.com/entry/2021/11/13/202103
    assert(!equal(a, Point<T>(0, 0)));
    assert(!equal(b, Point<T>(0, 0)));
    if ((Point<T>(0, 0) < Point<T>(a.y, a.x)) == (Point<T>(0, 0) < Point<T>(b.y, b.x))) {
        return (a.x * b.y == a.y * b.x) ? norm(a) < norm(b) : a.x * b.y > a.y * b.x;
    } else {
        return Point<T>(a.y, a.x) > Point<T>(b.y, b.x);
    }
}
// |p| ^ 2
template <typename T> inline T norm(const Point<T> &p) { return p.x * p.x + p.y * p.y; }
// |p|
template <typename T> inline T abs(const Point<T> &p) { return std::sqrt(norm(p)); }
// arg
template <typename T> inline T arg(const Point<T> &p) { return std::atan2(p.y, p.x); }
// polar
template <typename T> inline Point<T> polar(const T &rho, const T &theta = T(0)) { return rotate(Point<T>(rho, 0), theta); }
// EPS
template <> double Point<double>::EPS = 1e-9;
template <> long double Point<long double>::EPS = 1e-12;
template <> long long Point<long long>::EPS = 0;
template <> __int128_t Point<__int128_t>::EPS = 0;
// change EPS
// using Double = double;
// using Pt = Point<Double>;
// Point<Double>::set_eps(new_eps);
#line 4 "geometry/polygon.hpp"

// polygon
template <typename T> using Polygon = std::vector<Point<T>>;
template <typename T> std::istream &operator>>(std::istream &is, Polygon<T> &p) {
    for (auto &&pi : p) is >> pi;
    return is;
}
template <typename T> std::ostream &operator<<(std::ostream &os, const Polygon<T> &p) {
    for (auto &&pi : p) os << pi << " -> ";
    return os;
}
#line 2 "geometry/polygon_is_convex.hpp"

#line 2 "geometry/ccw.hpp"

#line 4 "geometry/ccw.hpp"

// counter clockwise
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_1_C
constexpr int COUNTER_CLOCKWISE = 1;  // a-b-c counter clockwise
constexpr int CLOCKWISE = -1;         // a-b-c clockwise
constexpr int ONLINE_BACK = 2;        // c-a-b line
constexpr int ONLINE_FRONT = -2;      // a-b-c line
constexpr int ON_SEGMENT = 0;         // a-c-b line
template <typename T> int ccw(const Point<T> &a, Point<T> b, Point<T> c) {
    b = b - a, c = c - a;
    if (sign(cross(b, c)) == 1) return COUNTER_CLOCKWISE;
    if (sign(cross(b, c)) == -1) return CLOCKWISE;
    if (sign(dot(b, c)) == -1) return ONLINE_BACK;
    if (norm(b) < norm(c)) return ONLINE_FRONT;
    return ON_SEGMENT;
}
#line 5 "geometry/polygon_is_convex.hpp"

// check polygon is convex (not strictly, 0 <= angle <= 180 degrees)
// angle = 180 degrees -> ON_SEGMENT
// angle = 0 degrees -> ONLINE_FRONT or ONLINE_BACK
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_3_B
template <typename T> bool polygon_is_convex(const Polygon<T> &p) {
    int n = (int)p.size();
    assert(n >= 3);
    bool okccw = true, okcw = true;
    for (int i = 0; i < n - 2; i++) {
        int res = ccw(p[i], p[i + 1], p[i + 2]);
        if (res == CLOCKWISE) okccw = false;
        if (res == COUNTER_CLOCKWISE) okcw = false;
        if (!okccw and !okcw) return false;
    }
    {
        int res = ccw(p[n - 2], p[n - 1], p[0]);
        if (res == CLOCKWISE) okccw = false;
        if (res == COUNTER_CLOCKWISE) okcw = false;
        if (!okccw and !okcw) return false;
    }
    {
        int res = ccw(p[n - 1], p[0], p[1]);
        if (res == CLOCKWISE) okccw = false;
        if (res == COUNTER_CLOCKWISE) okcw = false;
        if (!okccw and !okcw) return false;
    }
    return true;
}
#line 5 "geometry/convex_polygon_diameter.hpp"

// convex polygon diameter
// http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=CGL_4_B
// return {index1, index2, diameter}
// using the method of rotating calipers (https://en.wikipedia.org/wiki/Rotating_calipers)
// complexity: O(n)
template <typename T> std::tuple<int, int, T> convex_polygon_diameter(const Polygon<T> &p) {
    assert(polygon_is_convex(p));
    int n = (int)p.size();
    assert(n >= 2);
    if (n == 2) {
        return {0, 1, abs(p[0] - p[1])};
    }
    auto [it_min, it_max] = std::minmax_element(p.begin(), p.end(), compare_x<T>);
    int idx_min = it_min - p.begin();
    int idx_max = it_max - p.begin();

    T maxdis = norm(p[idx_max] - p[idx_min]);
    int maxi = idx_min, i = idx_min, maxj = idx_max, j = idx_max;
    do {
        int ni = (i + 1 == n ? 0 : i + 1), nj = (j + 1 == n ? 0 : j + 1);
        if (sign(cross(p[ni] - p[i], p[nj] - p[j])) < 0) {
            i = ni;
        } else {
            j = nj;
        }
        if (norm(p[i] - p[j]) > maxdis) {
            maxdis = norm(p[i] - p[j]);
            maxi = i;
            maxj = j;
        }
    } while (i != idx_min || j != idx_max);
    return {maxi, maxj, abs(p[maxi] - p[maxj])};
}

Back to top page