本系列文章力求以简洁易懂的文字介绍计算几何中的基本概念，使读者快速入门，故不追求难度和深度，仅起到抛砖引玉的作用。
好久没有写计算几何题目了~（悲），最近想找题目好好练一练，因此撰写这份博客🤔
这道题目便是 2023 ICPC 国际大学生程序设计竞赛亚洲区域赛（南京站）的 Problem B，当时的比赛赛场上呢是无人解出，这道题的难度属于夺冠水平~
首先我们来看一下这道题目的题意，简化后就是：
1.给定一个矩形（OBB），要求找到一个坐标轴平行矩形（AABB） ；
2.两个矩形的相交面积比上面积并就是所谓的交并比（IOU）；
3.我们需要求出这个矩形使IOU最大。
上述是对一些交并比求解的可以在草稿纸上模拟的几种情况（读者们可以自己试试），总结出来重合部分多边形的边数在4-8，我们在求解的时候只需要多边形类求交就行了~😄
我们一般怎么求交并比呢🤔？
这个就是我们之前所学习的基本计算几何内容了：
多边形面积交&多边形面积并,由于思路较为基础，基础部分已经做介绍，这里我们直接上源码可做调用：
1.我们一般采用vector存储多边形，点类 or 线类（不含极角排序）：
2.多边形求交：
namespace PolyIntersectArea{
//ConvexPolygonIntersectArea
db CPIA(Polygon a,Polygon b)
{
Polygon p,tmp;
p = b;
for(int i = 0;i < a.sze()n&& b.sze() > 2;i++)
{
int sflag = sgn((a[i + 1] - a[i]) ^ (p[0] - a[i])),eflag;
tmp = Polygon(0);
for(int j = 0;j < b.sze();j++,sflag = eflag)
{
if(sflag >= 0) tmp.add(p[j]);
eflag = sgn((a[i + 1] - a[i]) ^ (p[j + 1] - a[i]));
if((sflag ^ eflag) == -2)
tmp.add(Line(a[i],a[i + 1]).crossPoint(Line(p[j],p[j + 1])));
}
p = tmp;
b.p.resize(tmp.sze());
}
if(b.sze() < 3) return 0.0;
return p.getArea();
}
//simplePolygonIntersectArea
db SPIA(Polygon a,Polygon b)
{
Polygon t1(3),t2(3);
db res = 0,num1,num2;
t1[0] = a[0],t2[0] = b[0];
for(int i = 2;i < a.sze();i++)
{
t1[1] = a[i - 1],t1[2] = a[i];
num1 = sgn((t1[1] - t1[0]) ^ (t1[2] - t1[0]));
if(num1 < 0) swap(t1[1],t1[2]);
for(int j = 2;j < b.sze();j++)
{
t2[1] = b[j - 1],t2[2] = b[j];
num2 = sgn((t2[1] - t2[0]) ^ (t2[2] - t2[0]));
if(num2 < 0) swap(t2[1],t2[2]);
res = res + CPIA(t1,t2) * num1 * num2;
}
}
return res;
}
}3.多边形求并：
namespace PolyUnion{
Point dir(Line in)
{
return in.p - in.v;
}
db pos(Point p,Line in)
{
return ((p - in.p) * dir(in)) / dir(in).len2();
}
db gao(vector <Polygon> po)
{
int n = po.size();
db res = 0;
for(int i = 0;i < n;i++){
for(int ii = 0;ii < po[i].sze();ii++)
{
Point A = po[i][ii],B = po[i][(ii + 1) % p[i].sze()];
Line AB = Line(A,B);
vector<pair<db , int>>c;
for(int j = 0;j < n;j++)
if(i != j)
{
for(int jj = 0;jj < po[j].sze();jj++)
{
Point C = po[j][jj],D = po[j][(jj + 1) % po[j].sze()];
Line CD = Line(C,D);
int f1 = sgn((B - A) ^ (C - A));
int f2 = sgn((B - A) ^ (D - A));
if(! f1 && ! f2)
{
if(i < j && sgn((dir(AB) * dir(CD))) > 0)
{
c.push_back(make_pair(pos(C,AB),1));
c.push_back(make_pair(pos(D,AB),-1));
}
continue;
}
db s1 = (D - C) ^ (A - C);
db s2 = (D - C) ^ (B - C);
db t = s1 / (s1 - s2);
if(f1 >= 0 && f2 < 0) c.push_back(make_pair(t,1));
if(f1 < 0 && f2 >= 0) c.push_back(make_pair(t,-1));
}
}
c.push_back(make_pair(0.0,0));
c.push_back(make_pair(1.0,0));
sort(c.begin(),c.end());
db s = 0.5 * (A ^ B),z = min(max((db)c[0].second,(db)0.0),(db)1.0);
for(int j = 1,k = c[0].second;j < (int)c.size();j++)
{
db w = min(max(c[j].first,(db)0.0),(db)1.0);
if(k == 0) res = res + s * (w - z);
k = k + c[j].second;
z = w;
}
}
}
return fabs(res);
}
}以上就是多边形面积交并的分别求法，(⊙﹏⊙)其实我们可以明显发现，讲这个归讲这个，在这道题目中呢我们其实不需要这么麻烦，因为我们已经知道两个都是矩形了，到时候直接写一段大模拟求交点和多边形就OK了😄！
好我们来剖析一下这道题的类型：
这是一道函数求极值的问题~
我们一般对于这种问题会有一下考虑：
1. 是否单调/单峰/凸/ 考虑 ——> 二分？三分？
2. 是否可导/可微？=> 解方程（组）/梯度下降/牛顿迭代
3. 局部最小陷阱？=> 爬山/模拟退火
4. 有没有约束？=> 拉格朗日乘子法/KKT/线性规划...
解题步骤以及相关技巧：
1.利用对称性可以简化问题。答案显然只和 OBB 的长宽比以 及角度有关，我们对 OBB 做平移，镜像都不影响答案。
2.不妨把这个 OBB 平移，使得它的中心位于原点。
3.AABB 参数化为 (x, y, w, h)，其中心位于 (x, y), 上下左右边界分别为y + (h / 2),y - (h / 2),x - (w / 2),x + (w / 2)。
所以有 IOU ( 0 , 0 , w , h ) >= IOU ( x , y , w , h ),AABB的中心放在原点，也就是与OBB的中心重合的位置我们视为最优点。
4.设 g(w, h) = IOU(0, 0,w, h)，u(h) 是面积并关于 h 的导数， 那 g(w, h) 可看成关于 h 的函数表示如下:
5.注意到 u(h) 是关于 h 单调不递减且满足 0 ≤ u(h) ≤ w, 这 里的分子也是单调不递增的，而分母恒大于 0，故 g 关于 h是单峰的。对称地，g 关于 w 也是单峰的。
6.g(w, h) 整体是单峰函数。可以利用三分法，爬山法，梯度下降法等求解，都可以收敛到全局最优解。接下来介绍公式方法:
为了方便讨论，可以通过镜像以及等比缩放等操作（不影响答案），使得 OBB 长边为 cos ϕ，短边 sin ϕ，长边和 x-轴的夹角为 θ。(0 <= θ,ϕ <= 45°）
进行分类讨论，最优解时，交集可能是矩形（θ = 0），六边形（OBB 细长且倾斜度较大时）和八边形。我们也可以对每种情况列出函数，然后求解梯度等于 0 的方程来获得公式解：
解法雀实非常的明确且抽象,代码如下·,可以好好理解一下：
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef double LD;
typedef pair <int, int> pii;
#define cp const point &
const LD eps = 0;
int sgn (LD x) { return x > eps ? 1 : (x < -eps ? -1 : 0); }
LD sqr (LD x) { return x * x; }
struct point {
LD x, y;
point () {}
point (LD xx, LD yy) { x = xx, y = yy; }
point operator + (cp a) const { return {x + a.x, y + a.y}; }
point operator - (cp a) const { return {x - a.x, y - a.y}; }
point operator * (LD a) const { return {x * a, y * a}; }
point operator / (LD a) const { return {x / a, y / a}; }
point rot90() const { return {-y, x}; }
bool operator == (cp a) const {
return x == a.x && y == a.y;
}
void read() {
int xx, yy;
cin >> xx >> yy;
x = xx, y = yy;
}
};
LD det (cp a, cp b) { return a.x * b.y - b.x * a.y; }
LD dot (cp a, cp b) { return a.x * b.x + a.y * b.y; }
LD dis (cp a, cp b) { return sqrt (sqr (a.x - b.x) + sqr(a.y - b.y)); }
bool turn_left(cp a, cp b, cp c) {
return sgn (det (b - a, c - a)) > 0;
}
#define cl const line &
struct line {
point s, t;
line () {}
line (point ss, point tt) { s = ss, t = tt; }
bool operator == (cl a) const {
return s == a.s && t == a.t;
}
};
point line_inter (cl a, cl b) {
LD s1 = det (a.t - a.s, b.s - a.s);
LD s2 = det (a.t - a.s, b.t - a.s);
return (b.s * s2 - b.t * s1) / (s2 - s1);
}
bool turn_left (cl l, cp p) { return turn_left(l.s, l.t, p); }
line h[8];
LD hpi_nosort() {
line q[8]; int l = 0, r = -1;
point ret[8];
q[0] = h[0]; q[1] = h[1]; r = 1;
ret[1] = line_inter(q[0], q[1]);
for (int t = 2; t < 8; t++) {
auto &i = h[t];
while (l < r && !turn_left(i, ret[r]))
-- r;
while (l < r && !turn_left(i, ret[l + 1]))
++ l;
++ r; q[r] = i;
if (l != r) ret[r] = line_inter(q[r - 1], q[r]);
}
ret[l] = line_inter(q[r], q[l]);
LD area = 0;
for (int i = l; i <= r; i++) area += det(ret[i], ret[i == r ? l : i + 1]);
return area;
}
LD ans, sq;
LD check (LD x, LD y) {
h[1] = {{-1, -y}, {1, -y}};
h[3] = {{x, -1}, {x, 1}};
h[5] = {{1, y}, {-1, y}};
h[7] = {{0, 1}, {0, -1}};
LD a = hpi_nosort();
LD f = x * y * 4;
return a / (f + sq - a);
}
const LD Rp = (sqrt(5) - 1) / 2, Lp = 1 - Rp;
LD check (LD x) {
LD l = 0, r = 1, lv = -1, rv = -1, v = 0;
for (int t = 42; t; t--) {
LD b = (r - l) * Rp;
LD lmid = r - b, rmid = l + b;
if (lv == -1) lv = check(x, lmid);
if (rv == -1) rv = check(x, rmid);
if (rv < lv) rv = lv, lv = -1, r = rmid;
else lv = rv, rv = -1, l = lmid;
}
ans = max(ans, v = max(lv, rv));
return v;
}
void work() {
vector <point> p;
for (int i = 0; i < 4; i++) {
point u; u.read();
p.push_back(u);
}
if (!turn_left(p[0], p[1], p[2])) {
reverse(p.begin(), p.end());
}
if (sgn(det(p[1] - p[0], {1, 0})) == 0 || sgn(det(p[2] - p[1], {1, 0})) == 0) {
cout << "1\n";
return;
}
auto center = (p[0] + p[2]) / 2;
{
int k = 0;
for (int i = 0; i < 4; i++) {
if (p[i].x < p[k].x) k = i;
}
vector <point> q;
for (int i = 0; i < 4; i++) q.push_back(p[(i + k) % 4] - center);
p = move(q);
}
sq = abs(det(p[1] - p[0], p[2] - p[1]));
LD Lx = p[2].x;
LD Ly = p[3].y;
for (auto &[x, y] : p) {
x /= Lx;
y /= Ly;
}
sq /= Lx * Ly;
h[0] = {p[0], p[1]};
h[2] = {p[1], p[2]};
h[4] = {p[2], p[3]};
h[6] = {p[3], p[0]};
LD xl = 0, xr = 1, lv = -1, rv = -1;
ans = check(1 - 1e-9);
for (int t = 38; t; t--) {
LD b = (xr - xl) * Rp;
LD lmid = xr - b, rmid = xl + b;
if (lv == -1) lv = check(lmid);
if (rv == -1) rv = check(rmid);
if (rv < lv) rv = lv, lv = -1, xr = rmid;
else lv = rv, rv = -1, xl = lmid;
}
check((xl + xr) / 2);
cout << (double)ans << '\n';
}
int main() {
ios::sync_with_stdio(false); cin.tie(0);
int T = 1;
cin >> T;
cout << fixed << setprecision(10);
for (int ca = 1; ca <= T; ca ++) {
work();
}
}
扫一扫
1735
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
