FOSSIL (상, p. 486) - 다시 풀어봄.
문제#
기하로 풀어보기#
두 볼록 다각형(Convex Hull)이 겹치는 구간도 볼록 다각형이다. 이 새로운 볼록 다각형에서 남북 방향의 최대 길이를 구하면 되는데, 남북 방향의 최대 길이는 새로운 볼록 다각형의 꼭지점에서만 생성되므로, 각 꼭지점 에서의 길이중 최대 값을 출력하면 된다.
상세 방법은 아래와 같다.
1. 교집합 볼록 다각형(Convex Hull)을 구한다.#
1) 두 볼록 다각형을 점이 아닌 선분으로 나타낸다.#
선분으로 따로 나타내면 반복문 돌리기가 쉽다.
2) 두 볼록 다각형의 교점을 구한다.#
두 선분 집합을 이중 for문으로 돌려서 교점이 있는지 확인하고, 교점이 있으면 해당 교점을 구한다. 교점은 벡터의 외적을 활용하여 구한다.
① 먼저 교차점의 존재 여부부터 파악해야 한다.
2차원에서 두 벡터 $\vec{u}=(x_1,y_1),\; \vec{v}=(x_2,y_2)$의 외적은 다음과 같은 스칼라 값으로 정의한다.
$$ \vec{u} \times \vec{v} = x_1y_2 - y_1x_2 $$그리고 이 결과값의 부호가 방향을 의미한다.
> 0: 반시계(CCW)< 0: 시계(CW)= 0: 같은 직선(colinear)
참고로 이 결과값의 크기(절대값)는 두 벡터가 이루는 평행사변형의 넓이를 의미한다.
아래 그림을 보면, $\overline{AB}$와 $\overline{CD}$는 교차한다.
그러면 $\overrightarrow{AB} \times \overrightarrow{AC} \;$와 $\overrightarrow{AB} \times \overrightarrow{AD} \;$는 부호가 다를 수밖에 없다. $\left((\overrightarrow{AB} \times \overrightarrow{AD})\cdot (\overrightarrow{AB} \times \overrightarrow{AD}) \le 0\right)$
그런데 이 것만 고려하면, 두 선분 $\overline{AB}$와 $\overline{CD}$가 만나지 않는 경우를 제대로 고려할 수 없다. 교차점 존재 여부를 파악하기 위해선 $(\overrightarrow{CD} \times \overrightarrow{CA})\cdot (\overrightarrow{CD} \times \overrightarrow{CB}) \le 0$도 만족하는지 봐야한다.
② 교차점 존재 여부를 파악했으면, 교차점($I$) 을 구한다.
$\overline{AB}$위의 모든 점은 아래와 같이 벡터와 파라미터를 사용해서 나타낼 수 있다.
$$ P(t) = A + t(B-A),\quad 0\le t \le 1 $$$\overline{CD}$도 마찬가지로 아래와 같이 나타낼 수 있다.
$$ Q(u) = C + w(D-C), \quad 0 \le u \le 1 $$그러므로 교차점 $I$는 아래의 식을 풀면 된다.
$$ A + t(B-A) = C+u(D-C) $$$$ \Rightarrow t(B-A) = (C-A) + u(D-C) $$양변에 $(D-C)$로 외적을 취한다.
$$ \Rightarrow t(B-A)\times (D-C) = (C-A)\times(D-C) $$왜냐면, $u(D-C) \times (D-C) = 0$ 이기 때문이다. (자기 자신과의 외적은 항상 0)
$$ \therefore t = \frac{(C-A)\times(D-C)}{(B-A)\times (D-C)} $$여기서 분모가 0이 된다는 것은 $\overline{AB} \parallel \overline{CD}$이거나 두 선분이 한 직선위에 있는 상황인데, 두 선분이 한 직선위에 있는 상황은 문제에서 주어지지 않는다 하였고(”한 직선 위에 있는 세 점이 주어지는 일도 없습니다.”), 위의 방식으로 교차점을 구하기 전에 교차점이 있는지 이미 확인한 상태이므로 분모가 0이 되는 상황은 발생하지 않는다.
이제 위의 t값을 P(t)에 대입하면 교차점 $I$를 구할 수 있다.
3) 한쪽 볼록 다각형이 포함되는 나머지 볼록 다각형의 정점들을 구한다.#
이제 교집합 볼록 다각형을 이루는 나머지 점들을 구해야 하는데, 이러한 점들은 볼록 다각형에 포함되는 다른 볼록 다각형의 정점들이다. 그래서 볼록 다각형 V의 정점들이 다른 볼록 다각형인 W안에 포함되는지를 확인하고, 그 반대로도 포함여부를 확인해주면 된다.
정점 $p$가 볼록 다각형 $V$에 포함되는지 확인하는 방법은 아래와 같다.
모든 $i$에 대해 점 $p$로 부터 항상 $(v[i\text{+}1]-v[i]) \times (p -v[i]) \le 0$ 이면, 정점 p가 볼록 다각형 V에 포함된다.
코드로 구현하면 아래와 같다.
void SortCvhByCcwOrder(vector<Vec>& u) {
//모든 점의 중심점을 계산한다.
double midx, midy;
Vec mid(0.0, 0.0);
for (int i = 0; i < u.size(); i++) {
mid = mid + u[i];
}
mid = mid * (1.0/u.size());
//중심점과 Convex Hull의 점점들 간의 PolarAngle값을 구하고, 그 값 순서로 정렬한다.
auto cmp = [mid](const Vec& lhs, const Vec& rhs) {
double lang = (lhs - mid).PolarAngle();
double rang = (rhs - mid).PolarAngle();
if (fabs(lang - rang) < ES) {
return (lhs - mid).Len() < (rhs - mid).Len();
}
else return lang < rang;
};
sort(u.begin(), u.end(), cmp);
}2. 교집합 볼록 다각형의 남북방향 최대 길이를 계산한다.#
1) 교집합 볼록 다각형의 정점개수가 1이하면, 남북방향 최대 길이는 0이다.#
2) 교집합 볼록 다각형을 CCW(시계 반대 방향) 방향으로 정렬한다.#
남북방향 최대 길이를 쉽게 계산하기 위해, 교집합 볼록 다각형 U를 윗껍질과 아랫껍질로 나누려고 한다. 이렇게 나누게 되면 x좌표에 해당하는 윗껍질의 y좌표와 아랫껍질의 y좌표값을 뺀 것이 좌표 x에서의 남북방향 길이가 된다.
이렇게 윗껍질과 아랫껍질로 나누려면, 볼록 다각형의 정점들이 CCW 방향으로 정렬되어 있어야 한다. 정렬 방법은 볼록 다각형 정점들의 평균이 되는 중앙점(mid)을 구하고, 이 중앙점과 볼록 다각형의 모든 정점간의 각도 순으로 정렬하면 된다. (각도는 x축부터 반시계방향으로 $[0,2\pi)$값을 갖는다.)
x=a꼴의 선분의 경우 윗껍질/아랫껍질 어느 것에도 포함시키지 않는다. 이런 경우 어차피 x=a와 맞닿아 있는 선분들을 고려하면 x=a에 해당하는 남북방향 길이가 구해지기 때문이다.
3) 교집합 볼록 다각형의 정점마다 남북방향의 길이를 측정하고 최댓값을 구한다.#
윗껍질 함수: $y = f(x)$ 아랫껍질 함수: $y = g(x)$
라고 할 때, $max(f(x)-g(x))$를 구하면 된다.
정답 코드 (기하)#
#include <bits/stdc++.h>
using namespace std;
constexpr double ES = 1e-8; //Epsilon
constexpr double NA = 105.0; //Not Available
inline double PI() {
static const double pi = acos(-1.0);
return pi;
}
inline int Sign(double x) {
return (x>-ES)-(x<ES);
}
struct Vec {
double x,y;
Vec() {}
Vec(double x, double y):x(x),y(y) {}
bool operator<(const Vec & rhs) const {
return (x!=rhs.x ? x<rhs.x : y<rhs.y);
}
bool operator==(const Vec & rhs) const {
return x==rhs.x && y==rhs.y;
}
bool operator<=(const Vec & rhs) const {
return (*this < rhs || *this == rhs);
}
Vec operator*(double k) const {
return Vec(k*x, k*y);
}
Vec operator-(const Vec & rhs) const {
return Vec(x-rhs.x,y-rhs.y);
}
Vec operator+(const Vec & rhs) const {
return Vec(x+rhs.x,y+rhs.y);
}
double PolarAngle() const {
return fmod(atan2(y,x)+2*PI(), 2*PI());
}
double Dot(const Vec & rhs) const {
return x*rhs.x + y*rhs.y;
}
double Cross(const Vec & rhs) const {
return x*rhs.y - rhs.x*y;
}
double Len() const {
return hypot(x,y);
}
Vec UnitVec() const {
return Vec(x/Len(), y/Len());
}
Vec ProjectTo(const Vec & rhs) const {
Vec uv = rhs.UnitVec();
return uv*(this->Dot(uv));
}
int Ccw(const Vec & rhs) const {
//return값은 아래와 같다.
// 1) rhs가 *this의 CCW방향 위치에 있을 때: 1
// 2) rhs가 *this의 CW방향 위치에 있을 때: -1
// 3) rhs가 *this와 방향이 동일하거나 180도 정반대 방향일 때: 0
return Sign(this->Cross(rhs));
}
void Print() const {
printf("(%.3lf, %.3lf)\n", x,y);
}
};
struct LineSeg { //선분
Vec l,r;
LineSeg() {}
LineSeg(const Vec & lhs, const Vec & rhs) {
if(rhs < lhs) l=rhs, r=lhs;
else l=lhs, r=rhs;
}
bool HasIntersection(const LineSeg & rhs, Vec & intersection) const {
Vec p = r-l;
Vec a = rhs.r-l;
Vec b = rhs.l-l;
Vec q = rhs.r-rhs.l;
Vec c = r-rhs.l;
Vec d = l-rhs.l;
if(p.Ccw(a)*p.Ccw(b)==0 && q.Ccw(c)*q.Ccw(d)==0) { //두 선분이 한 직선위에 있는 경우
if(r<rhs.l || rhs.r<l) return false;
else {
//겹치는 부분이 최소 한 점 이상인 경우, 교차점을 아무거나 return
if(rhs.l <= r) intersection = rhs.l;
else intersection = r;
return true;
}
}
else {
if(p.Ccw(a)*p.Ccw(b)<=0 && q.Ccw(c)*q.Ccw(d)<=0) {
double t = (rhs.l-l).Cross(rhs.r-rhs.l) / (r-l).Cross(rhs.r-rhs.l);
intersection = l + ((r-l)*t);
return true; //p x a(외적)와 p x b의 부호가 다르면 교차점이 존재한다.
}
else return false;
}
}
bool IsBetweenX(double x) const {
double mn = min(l.x, r.x);
double mx = max(l.x, r.x);
return mn <= x && x <= mx;
}
//*this 선분에서 x좌표값에 해당하는 y좌표(f(x)) 값을 return한다.
double F(double x) const {
assert(r.x != l.x);
if(IsBetweenX(x)) {
double y = (r.y - l.y) / (r.x - l.x) * (x - l.x) + l.y;
if (l.y <= r.y) {
if (l.y <= y && y <= r.y) return y;
}
else {
if (r.y <= y && y <= l.y) return y;
}
}
return NA;
}
};
double F(double x, vector<LineSeg> & ls) {
for(int i=0; i<ls.size(); i++) {
if(ls[i].r.x == ls[i].l.x) continue;
double y = ls[i].F(x);
if(y >= NA) continue;
return y;
}
return NA;
}
//정점p가 Ccw로 정렬된 Convex Hull 안에 포함되는지를 본다.
//<=0으로 비교했기 때문에 Convex Hull 경계에 있어도 true를 리턴한다.
bool IsIn(Vec & p, vector<Vec> & ccw_cvh) {
for(int i=1; i<ccw_cvh.size(); i++) {
if((ccw_cvh[i]-ccw_cvh[i-1]).Ccw(p-ccw_cvh[i-1]) <= 0) return false;
}
return true;
}
void SortCvhByCcwOrder(vector<Vec>& u) {
//모든 점의 중심점을 계산한다.
double midx, midy;
Vec mid(0.0, 0.0);
for (int i = 0; i < u.size(); i++) {
mid = mid + u[i];
}
mid = mid * (1.0/u.size());
//중심점과 Convex Hull의 점점들 간의 PolarAngle값을 구하고, 그 값 순서로 정렬한다.
auto cmp = [mid](const Vec& lhs, const Vec& rhs) {
double lang = (lhs - mid).PolarAngle();
double rang = (rhs - mid).PolarAngle();
if (fabs(lang - rang) < ES) {
return (lhs - mid).Len() < (rhs - mid).Len();
}
else return lang < rang;
};
sort(u.begin(), u.end(), cmp);
}
vector<Vec> v,w,u;
vector<LineSeg> ls_v, ls_w;
vector<LineSeg> lower_ls, upper_ls;
int main() {
ios::sync_with_stdio(false); cin.tie(0);
int tc; cin >> tc;
while(tc--) {
v.clear(); w.clear();
ls_v.clear(); ls_w.clear();
int n,m; cin >> n >> m;
double x,y;
for(int i=0; i<n; i++) {
cin >> x >> y;
v.emplace_back(x,y);
}
v.emplace_back(v[0]);
for(int i=0; i<m; i++) {
cin >> x >> y;
w.emplace_back(x,y);
}
w.emplace_back(w[0]);
//Convex Hull을 선분으로 다시 표현.
for(int i=0; i<n; i++) {
ls_v.emplace_back(v[i],v[i+1]);
}
for(int i=0; i<m; i++) {
ls_w.emplace_back(w[i],w[i+1]);
}
//u: 새로운 Convex Hull 정점들
u.clear();
Vec intsc;
//1) v와 w의 교점을 u에 넣는다.
for(int i=0; i<n; i++) {
for(int j=0; j<m; j++) {
if(ls_v[i].HasIntersection(ls_w[j], intsc)) {
u.emplace_back(intsc);
}
}
}
//2) v의 정점이 w Convex Hull에 속하면 u에 넣는다.
// w의 정점이 v Convex Hull에 속하면 u에 넣는다.
for(int i=0; i<n; i++) {
if(IsIn(v[i], w)) u.emplace_back(v[i]);
}
for(int j=0; j<m; j++) {
if(IsIn(w[j], v)) u.emplace_back(w[j]);
}
//u의 정점개수가 1개이하면, 문제의 답은 항상 0이다.
if(u.size()<=1) {
printf("%.8lf\n", 0.0);
continue;
}
//u를 ccw 방향 순으로 정렬한다.
SortCvhByCcwOrder(u);
u.emplace_back(u[0]);
lower_ls.clear();
upper_ls.clear();
for(int i=0; i<u.size()-1; i++) {
if(u[i+1].x > u[i].x) lower_ls.emplace_back(u[i], u[i+1]);
else if(u[i+1].x < u[i].x) upper_ls.emplace_back(u[i+1], u[i]);
}
double ans = 0.0;
for(int i=0; i<u.size()-1; i++) {
ans = max(ans, F(u[i].x, upper_ls) - F(u[i].x, lower_ls));
}
printf("%.8lf\n", ans);
}
return 0;
}주의할 점#
아래와 같은 경우도 고려되었는지 확인해야 한다.
좌측의 그림은 두 Convex Hull간의 겹치는 부분에 두 Convex Hull의 꼭지점이 하나도 포함되어 있지 않은 경우도 있다는 것을 말한다.
우측의 그림은 한 Convex Hull이 다른 Convex Hull을 완전히 포함하는 경우를 말한다.
추가 Test Case와 Output#
//Input
2
3 3
0 0 6 2 0 4
2 2 4 0 4 4
4 4
20 20 80 20 80 80 20 80
40 40 60 40 60 60 40 60
//Output
2.00000000
20.00000000Comment#
기하 문제 푸는 방식에 대해 잘 모르고 풀었다가 책을 보고 삼분 탐색으로 풀어본 뒤 나중에 계산 기하 챕터를 공부하며 다시 풀어보게 되었다.
책 풀이 (삼분 탐색)#
이 문제는 수치 해석파트에 있는 문제다. 벡터 없이도 풀 수 있고, 책에서는 삼분 탐색으로 풀고 있다. 풀이 방법은 아래와 같다.
두 볼록 다각형(Convex Hull)의 윗 껍질과 아랫 껍질을 나누어 저장하는데, 어느 볼록 다각형의 껍질인지 구분할 필요가 없이 저장하고, 두 껍질에서의 x값에 따른 y값 차이의 최소를 삼분탐색으로 최대가 되는 지점을 찾는 것이다.
정답 코드 (삼분 탐색)#
#include <bits/stdc++.h>
using namespace std;
using pdd = pair<double, double>;
using ppddpdd = pair<pdd, pdd>;
vector<pdd> v, w;
vector<ppddpdd> upper, lower;
bool IsBetween(double x1, double x2, double x) {
return (x1 <= x && x <= x2) || (x2 <= x && x <= x1);
}
double At(ppddpdd& line, double x) {
double a = (line.second.second - line.first.second) / (line.second.first - line.first.first);
return a * (x - line.first.first) + line.first.second;
}
double GetUpperY(vector<ppddpdd>& upper, double x) {
double mn_y = 100;
for (int i = 0; i < upper.size(); i++) {
if (IsBetween(upper[i].first.first, upper[i].second.first, x)) {
mn_y = min(mn_y, At(upper[i], x));
}
}
return mn_y;
}
double GetLowerY(vector<ppddpdd>& lower, double x) {
double mx_y = 0;
for (int i = 0; i < lower.size(); i++) {
if (IsBetween(lower[i].first.first, lower[i].second.first, x)) {
mx_y = max(mx_y, At(lower[i], x));
}
}
return mx_y;
}
double GetMinX(vector<pdd>& g) {
double mn = 100;
for (int i = 0; i < g.size() - 1; i++) {
mn = min(mn, g[i].first);
}
return mn;
}
double GetMaxX(vector<pdd>& g) {
double mx = 0;
for (int i = 0; i < g.size() - 1; i++) {
mx = max(mx, g[i].first);
}
return mx;
}
int main() {
ios::sync_with_stdio(false); cin.tie(0);
int tc; cin >> tc;
while (tc--) {
int n, m;
cin >> n >> m;
v.clear(); w.clear();
lower.clear(); upper.clear();
double x, y;
for (int i = 0; i < n; i++) {
cin >> x >> y;
v.emplace_back(x, y);
}
for (int i = 0; i < m; i++) {
cin >> x >> y;
w.emplace_back(x, y);
}
v.push_back(v[0]);
w.push_back(w[0]);
for (int i = 1; i < v.size(); i++) {
if (v[i].first > v[i - 1].first) lower.emplace_back(v[i - 1], v[i]);
else if (v[i].first < v[i - 1].first) upper.emplace_back(v[i], v[i - 1]);
}
for (int i = 1; i < w.size(); i++) {
if (w[i].first > w[i - 1].first) lower.emplace_back(w[i - 1], w[i]);
else if (w[i].first < w[i - 1].first) upper.emplace_back(w[i], w[i - 1]);
}
sort(lower.begin(), lower.end());
sort(upper.begin(), upper.end());
double lo = max(GetMinX(v), GetMinX(w));
double hi = min(GetMaxX(v), GetMaxX(w));
for (int z = 0; z < 100; z++) {
double lmid = (lo * 2 + hi) / 3;
double rmid = (lo + hi * 2) / 3;
double lres = GetUpperY(upper, lmid) - GetLowerY(lower, lmid);
double rres = GetUpperY(upper, rmid) - GetLowerY(lower, rmid);
if (lres < rres) lo = lmid;
else hi = rmid;
}
double ans_x = lo;
double ans_len = GetUpperY(upper, ans_x) - GetLowerY(lower, ans_x);
if (ans_len < 0) ans_len = 0.0;
printf("%.8lf\n", ans_len);
}
return 0;
}