快包_分治法求解凸包问题

凸包问题(分治法)

题目简述

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;
}