凸包问题(分治法)
题目简述
P2742 [USACO5.1]圈奶牛Fencing the Cows
农夫约翰想要建造一个围栏用来围住他的奶牛,可是他资金匮乏。他建造的围栏必须包括他的奶牛喜欢吃草的所有地点。对于给出的这些地点的坐标,计算最短的能够围住这些点的围栏的长度。
输入数据的第一行是一个整数。表示农夫约翰想要围住的放牧点的数目 n。
第 2 到第 (n + 1) 行,每行两个实数,第 (i + 1) 行的实数 x, y分别代表第 i 个放牧点的横纵坐标。
算法思想
原题解释:
在二维坐标系上有 n 个点,现知它们的坐标,试求出用线段把所有的点围住需要的最短长度,即二维凸包的周长。
凸包问题的分治法:
本题主要突出采用分治法(快包)的策略,求出其凸包中的点。(如下所示)
分治法
把一个大问题分成几个结构相同的子问题,把子问题再分成几个更小的子问题……。然后我们就能用递归的方法,分别求这些子问题的解。最后把每个子问题的解“组装”成原来大问题的解。
在此题中我们很容易可以想出这种分治法的思路,首先我们能很清晰地发现,位于坐标系的最左边以及最右边的必定是解集中的点,即上图中的 P12 与P3 必定是解。
当我们连接P12(最左边的点)和P3(最右边的点)时,整个问题就自然地被划分成了上下两个结构相同的子问题了,采用递归的形式再分别求解上下两个子问题,最终合并子问题的解即可求出原题的解了。
子问题
又上述分析可知,原问题分分割成了如下的一个个子问题。
此时我们又能很清晰地发现,在子问题当中,距离分割线最远的一个点必然为解集中的点,即下图所示。
求解子问题数学知识
我们可以采用解析几何的知识,如果A(x1, y1),B(x2, y2),C(x3, y3)是平面上的任意三点时,那么三角形▲ABC的面积等于上面行列式绝对值的二分之一。
当且仅当C位于直线 的左侧时,该表达式的符号为正。则使用这个公式,我们可以轻易地检查出一个点是否位于两个点的左侧,并且可以求得这个点到这条直线的距离。
实现过程(附c++代码)
#include <cstdio>
#include <iostream>
#include <climits>
#include <vector>
#include <cmath>
#include <algorithm>
using namespace std;
const int MAXN = 1e5+10;
typedef struct node{
double x, y;
} point;
int n;
point p[MAXN];
vector<point> ans;
double cal(point a, point b, point c) {
return (a.x - c.x) * (b.y - c.y) - (b.x - c.x) * (a.y - c.y);
}
double dis(point p1, point p2) {
return sqrt((p2.y-p1.y)*(p2.y-p1.y)+(p2.x-p1.x)*(p2.x-p1.x));
}
double check(point a1,point a2,point b1,point b2)//检查叉积是否大于0,如果是a就逆时针转到b
{
return (a2.x-a1.x)*(b2.y-b1.y)-(b2.x-b1.x)*(a2.y-a1.y);
}
bool cmp(point p1, point p2) {
double tmp=check(p[1],p1,p[1],p2);
if(tmp>0)
return 1;
if(tmp==0&&dis(p[1],p1)<dis(p[1],p2))
return 1;
return 0;
}
void init() {
//cout << "please input the n: ";
cin >> n;
for(int i = 1; i <= n; i++) {
//cout << "please input the " << i << "'s value" << endl;
cin >> p[i].x >> p[i].y;
}
}
void divide(point a, point b, vector<point> vec) {
//cout << "a:" << a.x << a.y << endl;
//cout << "b:" << b.x << b.y << endl;
if(vec.empty()) {
ans.push_back(a);
//cout << " ans : " << a.x << a.y << endl;
return ;
}
vector<point> l, r;
double area = 0;
int index;
for(int i = 0; i < vec.size(); i++) {
double areat = cal(a,b,vec[i]);
if(areat > area) {
area = areat;
index = i;
}
}
for(int i = 0; i < vec.size(); i++) {
if(cal(a, vec[index], vec[i]) > 0) l.push_back(vec[i]);
if(cal(vec[index], b, vec[i]) > 0) r.push_back(vec[i]);
}
divide(a,vec[index],l);
divide(vec[index],b,r);
}
void solve() {
vector<point> l, r;
int lx, rx;
double xmin = INT_MAX, xmax = INT_MIN;//初始化
for(int i = 1; i <= n; i++) {
if(p[i].x < xmin) {
xmin = p[i].x;
lx = i;
}
if(p[i].x > xmax) {
xmax = p[i].x;
rx = i;
}
}
//cout << lx << " " << rx;
for(int i = 1; i <= n; i++) {
if(cal(p[lx],p[rx],p[i]) > 0) l.push_back(p[i]);
if(cal(p[lx],p[rx],p[i]) < 0) r.push_back(p[i]);
}
divide(p[lx],p[rx],l);
divide(p[rx],p[lx],r);
}
void print() {
sort(ans.begin(),ans.end(),cmp);
cout << "ans: " << endl;
for(int i = 0; i < ans.size(); i++) {
cout << ans[i].x << "," << ans[i].y << endl;
}
}
void print_len() {
sort(ans.begin(),ans.end(),cmp);
double len = 0;
for(int i = 0; i < ans.size(); i++) {
len += dis(ans[i],ans[(i+1)%ans.size()]);
}
printf("%.2lf",len);
}
int main() {
init();//初始化操作
solve();
print_len();
return 0;
}