最小矩形覆盖问题
解题思路:
本题要求用一个面积最小的矩形将所有点覆盖住,求的是一个最小覆盖矩形。
我们先求一下所有点的凸包,则最小覆盖矩形的边上是一定存在凸包上的点的。如果某一条边上不存在凸包上的点,就意味着这条边还可以向内平移得到一个更小的矩形,这与最小覆盖矩形矛盾了,因此最小覆盖矩形的每条边上都一定存在凸包上的点。
如果我们此时只确定某一个点在矩形的边上,那么此时可能存在的矩形还是有无穷多种,但是如果我们能确定某一条边在矩形的边上,那么此时的最小覆盖矩形就被唯一确定了,我们用另外三条边从凸包外面向内平移,直到每条边都被凸包上的某个点卡住,此时四条边构成的矩形就是此时的最小覆盖矩形。
因此我们可以尝试去猜想,最小覆盖矩形上的某条边是否和凸包上的某条边共线。换句话来说,就是最小覆盖矩形上的某条边中至少存在凸包上的两个点。
由于这里需要在图上做一些数学推导,这里就不详细证明了,简单说一下证明的思路,我们假设最小覆盖矩形上每条边上只有一个凸包上的点,那么此时我们在保证这四个点在边上的同时,去旋转这个矩形,这里如果用一个函数来表示矩形的面积,然后我们通过函数可以推导出矩形往顺时针或逆时针中的一个方向旋转,一定能让矩形变小,如果顺时针能变小就顺时针转,如果逆时针能变小就逆时针转,而由于里面矩形内部还有很多凸包的点,因此矩形不可能一直旋转下去,必然有一个条边会被一个其他的凸包上的点卡住,此时这条边上就有两个点了,并且矩形面积变小,与最小覆盖矩形矛盾。反证得出最小覆盖矩形的某一条边上一定有两个凸包上的点(因为定义矩形面积的函数需要配合图来理解,我大多情况都只写文字题解,主要是懒,所以只提供一个思路,可以自行证明或参考 y 总的视频讲解)
因此我们只需要枚举所有和凸包的某一条边共线的矩形即可。我们可以直接枚举凸包上的边,对于每条边 $ab$,我们将它作为矩形的底边界,对于矩形的上边界,对于某个点 $p$,如果 $p$ 是上边界上的点,等价于 $\vec{ab} \times \vec{ap}$ 最大。对于矩形的右边界,如果 $p$ 是右边界上的点,等价于 $\vec{ap}$ 在 $\vec{ab}$ 上的投影长度最长,这可以用点积来求,同理对于矩形的左边界,如果 $p$ 是左边界上的点,等价于 $\vec{ap}$ 在 $\vec{ab}$ 上的投影长度最短,也可以用点积来求。这样我们对于每条边都能求出另外三条边的位置,构成一个最小覆盖矩形,我们从每条边作为底边得出的最小覆盖矩形中取一个面积最小的就是答案。
另外可以发现,当我们逆时针枚举每条边时,假设当前边的三个边界上的点都找到,则当枚举到下一条边时。三个边界上的点也只会在逆时针方向上单调的向后走,所以我们就可以用一个四指针算法线性求出每条边对应的最小覆盖矩形,而这就是一个旋转卡壳的思想。
而当我们求出每条边对应的三个边界上的点后,我们还需要求出这个最小覆盖矩形的面积和四个顶点的坐标,假设底边为 $ab$,上边界上的点为 $p$,左边界上的点为 $u$,右边界上的点为 $v$。则此时矩形的高就是可以用三角形的面积来求,$\frac{\vec{ab} \times p}{2}$ 是三角形 $abp$ 的面积,$ab$ 的长度已知,我们就能求出 $abp$ 的高,也就是矩形的高。而矩形的宽就是 $\vec{uv}$ 在 $\vec{ab}$ 上的投影的长度。这样就能计算出矩形的面积。
然后还要求出四个顶点的坐标,我们先求出右下角的坐标,$a$ 的横坐标加上 $\vec{av}$ 在 $\vec{ab}$ 上投影的长度,就是右下角的坐标,然后我们知道矩形的长和宽,就可以从右下角推出剩下三个角的坐标。
综上所述,就是本题全部思路,计算几何题是需要基于图来分析的,这里建议画图来配合理解
C++ 代码
#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second
using namespace std;
typedef pair<double, double> PDD;
const int N = 50010;
const double eps = 1e-12;
const double INF = 1e20;
const double PI = acos(-1);
int n;
PDD q[N]; //存储每个点的坐标
int stk[N], top; //栈,存储凸包
bool st[N]; //记录每个点是否在凸包上
PDD res[5]; //逆时针存储最小覆盖矩形的四个顶点
double min_area = INF; //记录最小覆盖矩形的面积
int sign(double x) //符号函数
{
if(fabs(x) < eps) return 0;
if(x < 0) return -1;
return 1;
}
int dcmp(double x, double y) //比较函数
{
if(fabs(x - y) < eps) return 0;
if(x < y) return -1;
return 1;
}
PDD operator+(PDD a, PDD b) //重载向量加法
{
return {a.x + b.x, a.y + b.y};
}
PDD operator-(PDD a, PDD b) //重载向量减法
{
return {a.x - b.x, a.y - b.y};
}
PDD operator*(PDD a, double t) //重载数乘
{
return {a.x * t, a.y * t};
}
PDD operator/(PDD a, double t) //重载数除
{
return {a.x / t, a.y / t};
}
double operator*(PDD a, PDD b) //重载叉乘
{
return a.x * b.y - a.y * b.x;
}
double operator&(PDD a, PDD b) //重载点乘
{
return a.x * b.x + a.y * b.y;
}
double area(PDD a, PDD b, PDD c) //计算 ab 和 ac 构成的平行四边形的有向面积
{
return (b - a) * (c - a);
}
double get_len(PDD a) //计算向量 a 的模长
{
return sqrt(a & a);
}
double project(PDD a, PDD b, PDD c) //计算 ac 在 ab 上的投影的长度
{
return ((b - a) & (c - a)) / get_len(b - a);
}
PDD norm(PDD a) //求 a 的单位向量
{
return a / get_len(a);
}
PDD rotate(PDD a, double b) //将 a 顺时针旋转 b 弧度
{
return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}
void andrew()
{
sort(q, q + n); //将所有点按照横、纵坐标从小到大排序
//求下凸包
for(int i = 0; i < n; i++)
{
while(top >= 2 && sign(area(q[stk[top - 2]], q[stk[top - 1]], q[i])) <= 0)
if(area(q[stk[top - 2]], q[stk[top - 1]], q[i]) < 0) st[stk[--top]] = false;
else top--; //如果点在直线上,则不能重置标记
stk[top++] = i;
st[i] = true;
}
//求上凸包
st[0] = false;
for(int i = n - 1; i >= 0; i--)
{
if(st[i]) continue; //如果当前点已经在凸包上,直接跳过
while(top >= 2 && sign(area(q[stk[top - 2]], q[stk[top - 1]], q[i])) <= 0) top--;
stk[top++] = i;
}
top--; //起点加入了两次,删掉一次
}
void rotating_calipers() //旋转卡壳
{
for(int i = 0, a = 2, b = 1, c = 2; i < top; i++) //a 表示上边界的点,b 表示右边界的点,c 表示左边界的点
{
PDD d = q[stk[i]], e = q[stk[i + 1]];
while(dcmp(area(d, e, q[stk[a]]), area(d, e, q[stk[a + 1]])) < 0) a = (a + 1) % top; //更新 a
while(dcmp(project(d, e, q[stk[b]]), project(d, e, q[stk[b + 1]])) < 0) b = (b + 1) % top; //更新 b
if(!i) c = a; //最开始 c 要从 a 开始往右走,才能保证走到左边界
while(dcmp(project(d, e, q[stk[c]]), project(d, e, q[stk[c + 1]])) > 0) c = (c + 1) % top; //更新 c
PDD x = q[stk[a]], y = q[stk[b]], z = q[stk[c]];
double h = area(d, e, x) / get_len(e - d); //求矩形的高
double w = ((y - z) & (e - d)) / get_len(e - d); //求矩形的宽
if(h * w < min_area) //更新最小覆盖矩形
{
min_area = h * w; //更新最小覆盖矩形的面积
//更新最小覆盖矩形的四个顶点
res[0] = d + norm(e - d) * project(d, e, y); //右下角
res[3] = d + norm(e - d) * project(d, e, z); //左下角
PDD u = norm(rotate(e - d, -PI / 2));
res[1] = res[0] + u * h; //右上角
res[2] = res[3] + u * h; //左上角
}
}
}
int main()
{
scanf("%d", &n);
for(int i = 0; i < n; i++) scanf("%lf%lf", &q[i].x, &q[i].y);
andrew(); //求二维凸包
rotating_calipers(); //旋转卡壳
//找出最小的顶点坐标
int k = 0;
for(int i = 1; i < 4; i++)
if(dcmp(res[i].y, res[k].y) < 0 || (!dcmp(res[i].y, res[k].y) && dcmp(res[i].x, res[k].x) < 0))
k = i;
printf("%.5lf\n", min_area);
for(int i = 0; i < 4; i++, k++)
{
double x = res[k % 4].x, y = res[k % 4].y;
//C++ 中 -0.00...001 四舍五入后会得到 -0.00...000,因此需要特判
if(!sign(x)) x = 0;
if(!sign(y)) y = 0;
printf("%.5lf %.5lf\n", x, y);
}
return 0;
}