#include "Geometry/convex.hpp"
#include <bits/stdc++.h>
using namespace std;
namespace convex_algo
{
template <class T, typename = enable_if_t<is_arithmetic_v<T>>>
pair<Point<T>, Point<T>> find_furthest_pair(const Convex<T> &convex)
{
const vector<Point<T>> &ccw_points = convex.get_ccw_points();
assert(2 < ccw_points.size());
int left = distance(ccw_points.begin(), min_element(ccw_points.begin(), ccw_points.end()));
int right = distance(ccw_points.begin(), max_element(ccw_points.begin(), ccw_points.end()));
pair<Point<T>, Point<T>> result;
T max_dist = 0;
int i = left, j = right;
do
{
T dist = (ccw_points[i] - ccw_points[j]).norm();
if (max_dist < dist)
{
max_dist = dist;
result = pair(ccw_points[i], ccw_points[j]);
}
// i の次の点と j の次の点を見比べ、i と j が反時計回りに回るように更新
int next_i = (i + 1) % ccw_points.size(), next_j = (j + 1) % ccw_points.size();
if (Point<T>::cross(ccw_points[next_i] - ccw_points[i], ccw_points[next_j] - ccw_points[j]) < 0)
{
i = next_i;
}
else
{
j = next_j;
}
} while (i != left || j != right);
return result;
}
}
#line 1 "Geometry/point.hpp"
#include <bits/stdc++.h>
using namespace std;
template <class T, typename = enable_if_t<is_arithmetic_v<T>>>
class Point
{
complex<T> coordinate;
static constexpr long double EPS = 1e-8;
public:
Point() {}
Point(T x, T y) : coordinate(x, y) {}
Point(complex<T> coordinate) : coordinate(coordinate) {}
T real() const { return coordinate.real(); }
T imag() const { return coordinate.imag(); }
T norm() const { return std::norm(coordinate); }
static T dot(const Point &a, const Point &b)
{
return a.real() * b.real() + a.imag() * b.imag();
}
static T cross(const Point &a, const Point &b)
{
return a.real() * b.imag() - a.imag() * b.real();
}
Point operator+(const Point &other) const
{
return Point(this->coordinate + other.coordinate);
}
Point operator-(const Point &other) const
{
return Point(this->coordinate - other.coordinate);
}
Point operator*(const Point &other) const
{
return Point(this->coordinate * other.coordinate);
}
Point operator/(const Point &other) const
{
return Point(this->coordinate / other.coordinate);
}
bool operator<(const Point &other) const
{
if (this->real() != other.real())
return this->real() < other.real();
return this->imag() < other.imag();
}
enum class CCWResult
{
COUNTER_CLOCKWISE = 2,
CLOCKWISE = -2,
ON_LINE_BACK = 1,
ON_LINE_FRONT = -1,
ON_SEGMENT = 0
};
static CCWResult ccw(const Point &a, Point b, Point c)
{
b.coordinate = b.coordinate - a.coordinate, c.coordinate = c.coordinate - a.coordinate;
if (cross(b, c) > Point::EPS)
return CCWResult::COUNTER_CLOCKWISE;
if (cross(b, c) < -Point::EPS)
return CCWResult::CLOCKWISE;
if (dot(b, c) < 0)
return CCWResult::ON_LINE_BACK;
if (b.norm() < c.norm())
return CCWResult::ON_LINE_FRONT;
return CCWResult::ON_SEGMENT;
}
};
#line 3 "Geometry/convex.hpp"
using namespace std;
template <class T, typename = enable_if_t<is_arithmetic_v<T>>>
class Convex
{
using CCWResult = typename Point<T>::CCWResult;
vector<Point<T>> ccw_points;
set<Point<T>> upper_hull;
set<Point<T>> lower_hull;
bool allow_on_line;
T twice_area = 0;
Point<T> bounding_box_min, bounding_box_max;
T calculate_twice_area()
{
if (ccw_points.size() < 3)
{
return 0;
}
T result = 0;
for (int i = 0; i < ccw_points.size(); i++)
{
const Point<T> &p = ccw_points[i];
const Point<T> &next_p = ccw_points[(i + 1) % ccw_points.size()];
result += Point<T>::cross(p, next_p);
}
return result;
}
list<Point<T>> build_half_hull(const set<Point<T>> &points, bool upper)
{
assert(3 <= points.size());
CCWResult ccw_direction = (upper ? CCWResult::CLOCKWISE : CCWResult::COUNTER_CLOCKWISE);
list<Point<T>> half_hull;
for (auto &p3 : points)
{
if (half_hull.size() < 2)
{
half_hull.push_back(p3);
continue;
}
while (half_hull.size() >= 2)
{
auto p2 = half_hull.back();
half_hull.pop_back();
auto p1 = half_hull.back();
CCWResult ccw_result = Point<T>::ccw(p1, p2, p3);
bool is_valid_direction = (ccw_result == ccw_direction) || (allow_on_line && ccw_result == CCWResult::ON_LINE_FRONT);
if (is_valid_direction)
{
half_hull.push_back(p2);
break;
}
}
half_hull.push_back(p3);
}
return half_hull;
}
vector<Point<T>> merge_hulls(const list<Point<T>> &upper_points, const list<Point<T>> &lower_points)
{
assert(2 <= upper_points.size() && 2 <= lower_points.size());
vector<Point<T>> merged_hull(lower_points.begin(), lower_points.end());
// lower_pointsの最後とupper_pointsの最初が重複するので取り除く
merged_hull.pop_back();
merged_hull.insert(merged_hull.end(), upper_points.rbegin(), upper_points.rend());
// upper_pointsの最後とlower_pointsの最初が重複するので取り除く
merged_hull.pop_back();
return merged_hull;
}
void set_bounding_box()
{
const T MAX_INF = numeric_limits<T>::max();
const T MIN_INF = numeric_limits<T>::lowest();
T x_min = MAX_INF, y_min = MAX_INF, x_max = MIN_INF, y_max = MIN_INF;
for (const auto &point : ccw_points)
{
x_min = min(x_min, point.real());
y_min = min(y_min, point.imag());
x_max = max(x_max, point.real());
y_max = max(y_max, point.imag());
}
bounding_box_min = Point<T>(x_min, y_min);
bounding_box_max = Point<T>(x_max, y_max);
}
public:
Convex(vector<Point<T>> points, bool allow_on_line = false) : allow_on_line(allow_on_line)
{
// 重複点の削除 + ソート
set<Point<T>> point_st(points.begin(), points.end());
if (point_st.size() <= 2)
{
ccw_points = vector<Point<T>>(point_st.begin(), point_st.end());
upper_hull = point_st;
lower_hull = point_st;
twice_area = 0;
set_bounding_box();
return;
}
list<Point<T>> upper_hull = build_half_hull(point_st, true);
list<Point<T>> lower_hull = build_half_hull(point_st, false);
this->ccw_points = merge_hulls(upper_hull, lower_hull);
this->upper_hull = set<Point<T>>(upper_hull.begin(), upper_hull.end());
this->lower_hull = set<Point<T>>(lower_hull.begin(), lower_hull.end());
twice_area = calculate_twice_area();
set_bounding_box();
};
bool has_point(const Point<T> &point) const
{
return upper_hull.contains(point) || lower_hull.contains(point);
}
bool is_inside_bounding_box(const Point<T> &point) const
{
return bounding_box_min.real() <= point.real() && point.real() <= bounding_box_max.real() &&
bounding_box_min.imag() <= point.imag() && point.imag() <= bounding_box_max.imag();
}
const vector<Point<T>> &get_ccw_points() const
{
return ccw_points;
}
const T &get_twice_area() const
{
return twice_area;
}
const set<Point<T>> &get_upper_hull() const
{
return upper_hull;
}
const set<Point<T>> &get_lower_hull() const
{
return lower_hull;
}
};
#line 3 "Geometry/convex_algo/find_furthest_pair.hpp"
using namespace std;
namespace convex_algo
{
template <class T, typename = enable_if_t<is_arithmetic_v<T>>>
pair<Point<T>, Point<T>> find_furthest_pair(const Convex<T> &convex)
{
const vector<Point<T>> &ccw_points = convex.get_ccw_points();
assert(2 < ccw_points.size());
int left = distance(ccw_points.begin(), min_element(ccw_points.begin(), ccw_points.end()));
int right = distance(ccw_points.begin(), max_element(ccw_points.begin(), ccw_points.end()));
pair<Point<T>, Point<T>> result;
T max_dist = 0;
int i = left, j = right;
do
{
T dist = (ccw_points[i] - ccw_points[j]).norm();
if (max_dist < dist)
{
max_dist = dist;
result = pair(ccw_points[i], ccw_points[j]);
}
// i の次の点と j の次の点を見比べ、i と j が反時計回りに回るように更新
int next_i = (i + 1) % ccw_points.size(), next_j = (j + 1) % ccw_points.size();
if (Point<T>::cross(ccw_points[next_i] - ccw_points[i], ccw_points[next_j] - ccw_points[j]) < 0)
{
i = next_i;
}
else
{
j = next_j;
}
} while (i != left || j != right);
return result;
}
}
Geometry/convex_algo/find_furthest_pair.hpp