POJ1113 Wall

题目来源：http://poj.org/problem?id=1113
题目大意：

　　如图所示，给定N个顶点构成的一个多边形和一个距离值L。建立一个围墙，把这个多边形完全包含在内，且围墙距离多边形任一点的距离不超过指定的距离L。求出满足条件的围墙长度最小值。

 

输入：第一行为N和L。3 <= N <= 1000， 1 <= L <= 1000。接下来N行每行两个整数代表一个点的坐标，(-10000 <= Xi, Yi <= 10000) ，每个点都不重合，输入保证了多边形的每条边不会相交，顶点顺序为沿多边形顺时针方向遍历顺序。

输出：满足题述要求的最小长度，保留精度到整数位。PS：结果四舍五入即可。

---

Sample Input

9 100
200 400
300 400
300 300
400 300
400 400
500 400
500 200
350 200
200 200

Sample Output

1628

---

最近恰好碰到了凸包相关的问题，就翻出来这道最基础的给做了。

首先要想清楚的是题中给出的不一定是凸多边形，假设是凹多边形，当我们把围墙也设为像凹多边形的形状时可以使得围墙距离恰好满足限制，但这不是我们需要的最优解情况，因为此时的围墙长度不是最短的。围墙最短的情况对应的应该是：对所有顶点求出其凸包（直观简单的理解就是一个最小的凸多边形把所有的点都包含在其内部），然后把凸包往外扩距离L，这样得到的还不是最短的，还需要在原来凸包的各顶点处画一个半径为L的圆，然后与凸包外扩后的边与切点处相连（相当于把角给磨圆了？），这样最终组成的形状才是满足要求的最短围墙。

有了上述思路后，我们就能明白围墙的最短长度应该等于 凸包的边长 + 半径L的圆的周长。因为根据上面的分析，我们在纸上画一画就可以发现围墙到凸包之间的部分实际上是由一些矩形和一些扇形组成的，这些扇形恰好会组成一个圆。

---

然后剩下的问题就是如何求顶点集合的凸包了。 

对于凸包有两种形象的理解：

1. 把点看做是钉在板子上的钉子，用一个橡皮筋把所有顶点框在里面，橡皮筋的形状就是顶点的凸包。

2. 同样把点看做钉在板子上的钉子，拿一根绳子，找一个麻绳，一端绑在最外的一颗钉子上，然后拉着麻绳逆时针(或顺时针)一圈，把所有钉子包括进去。这样最终麻绳的形状也是顶点集的凸包。

上面的第二种理解实际上就是一种求凸包的最基础的算法卷包裹法(Gift Wrapping)的思想。但是这里面涉及到两个细节：

1. 第一个顶点的选择：选最左下的或最右上的都可。

2. 下一个顶点的选择：通过计算向量叉乘来选取。

回顾一下叉乘的意义：三维情况下，两向量的叉乘得到的是一个新的向量，方向是与原向量确定平面垂直的方向(有符号), 大小是原向量形成的平行四边形的面积。在二维情况下，数值大小的意义不变，符号表示原向量的方向关系。cross<a, b> = |a| * |b| * sin<a, b> > 0, 说明  a 到 b 的角为正向，即 b 在 a 的逆时针方向， 反之在其顺时针方向， 等于0说明平行。所以显然叉乘不满足交换律。此外，cross<(p,a), (p, b)> 的数值大小表示三角形pab的面积的两倍，所以当ab固定不变时，通过比较叉乘的数值大小可以比较点到线段的距离(这里的距离也是有符号的)。

所以利用叉乘运算可以解决卷包裹算法中的下一个顶点选择问题。一个特殊情况是：多点共线。共线时选点策略是：

卷包裹算法复杂度：O(n*h)。其中h为凸包上的点数。所以如果给出的点集，位于凸包上的点数比较少时，算法较快，但是最坏情况下可能所有点都位于凸包上，这样算法复杂度就是O(n^2)。

本题卷包裹算法的代码实现：

/////////////////////////////////////////////////////////////////
//        POJ1113 Wall    (卷包裹Gift Wrapping算法版)
//        Memory: 232K        Time: 79MS
//        Language: C++        Result : Accepted
/////////////////////////////////////////////////////////////////

#include <iostream>
#include <algorithm>
#include <cmath>

using namespace std;

#define PI 3.141592653

using namespace std;

struct Point {
    int x, y;
};

Point p[1000];
int n, L, ans[1000], cnt = 0;
bool mark[1000];

double dis(Point &p1, Point &p2) {
    double delta_x = p1.x - p2.x;
    double delta_y = p1.y - p2.y;
    return sqrt(delta_x * delta_x + delta_y * delta_y);
}

int cross(Point &ps, Point &pe, Point &qs, Point &qe) {
    //二维向量叉乘公式：
    //cross<(x1, y1), (x2, y2)> = x1 * y2 - x2 * y1 = || * |b| * sin<a, b>
    //cross<p1, p2> = |p1| * |p2| * sin<p1, p2>
    //> 0 表示 p1 在 p2 的 顺时针方向， < 0 在逆时针方向， = 0 表示共线
    //函数所求为向量 ps->pe 和 qs->qe 的叉乘
    return (pe.x - ps.x) * (qe.y - qs.y) - (qe.x - qs.x) * (pe.y - ps.y);
}

bool cmp(Point &a, Point &b, Point &v) {
    //比较到参考点极角大小，极角相等时比较距参考点的距离
    int res = cross(v, a, v, b);
    if (res == 0) {
        return dis(a, v) < dis(b, v);
    } else {
        return res > 0;
    }
}

void GiftWrapping(void) {
    int last = 0, next;
    ans[cnt++] = 0;
    //mark[0] = true;
    while (true) {
        //找到一个不在凸包上的点
        next = 0;
        for (int i = 1; i < n; ++i) {
            if (!mark[i]) {
                next = i;
                break;
            }
        }
        //与其它不在凸包上的点比较
        for (int i = next + 1; i < n; ++i) {
            if (!mark[i] && cmp(p[i], p[next], p[last])) {
                next = i;
            }
        }
        //没有不在凸包上的点了 或 应该与起点相连了
        if (next == 0 || last != 0 && cmp(p[0], p[next], p[last])) {
            ans[cnt++] = 0;
            break;
        } else {
            ans[cnt++] = next;
            mark[next] = true;
            last = next;
        }
    }
}

int main(void) {
    cin >> n >> L;
    double circle = 2 * L * PI;
    int ll = 0;
    for (int i = 0; i < n; ++i) {
        cin >> p[i].x >> p[i].y;
        //找出y坐标最小的中x坐标最小的点（左下角点，一定在凸包上）
        if (p[i].y < p[ll].y || p[i].y == p[ll].y && p[i].x < p[ll].x) {
            ll = i;
        }
    }
    //将左下角点放在p[0]处，作为参考点
    Point t = p[ll];
    p[ll] = p[0];
    p[0] = t;
    //运行核心Gift Wrapping算法
    GiftWrapping();
    //计算围墙最小周长
    double length = circle;
    while (cnt > 1) {
        length += dis(p[ans[cnt - 1]], p[ans[cnt - 2]]);
        --cnt;
    }
    //四舍五入
    cout << (int)(length + 0.5) << endl;
    return 0;
}

---

 另一种常用的求凸包算法：Graham扫描算法(Graham Scan Algorithm)。算法分为两部分：排序和扫描。

1. 排序。

首先找到最左下的点（纵坐标最小的点中横坐标最小的点），该点一定在凸包上，然后设该点为参考点。

求出其它每个点到该参考点的极角值。然后对所有点按极角大小排序。排序后的结果如下面的图中a所示。

由于参考点是最左下的点，所以其它所有点到它的极角会在180°的变化范围之间。

极角的求法仍然是利用向量叉乘。

2. 扫描。

需要借助一个栈。栈里保存的是顶点，从栈底到栈顶的顶点连接起来形成的是我们当前得到的凸包的一部分边。初始化时，前三个点（含参考点）入栈。然后按排序后的顺序依次扫描其它的点，不断地扩展凸包的边，直到围成一个凸多边形。扫描时的具体做法：（参考图片）每扫描到一个新的顶点，如果它与当前栈顶顶点连成的边与最靠近栈顶的两个顶点连成的边，从当前“半凸包”的内部看构成的是一个不大于180度的角，即新边与栈顶边呈“左转”关系，那么将扫描到的点入栈；若形成的是“右转”关系，则栈顶元素出栈，重新计算边的关系，重复此操作知道满足“左转”关系为止，然后继续扫描下一个顶点，直至所有顶点扫描完毕将顶点首尾相接得到一个完整的凸多边形。

算法的复杂度分析：排序过程O(nlg2n), 扫描过程O(n). 故最坏情况下的复杂度优于卷包裹算法。

Graham算法极角排序法的实现：

/////////////////////////////////////////////////////////////////
//        POJ1113 Wall    (Graham scan 算法, 极角排序版)
//        Memory: 236K        Time: 79MS
//        Language: C++        Result : Accepted
/////////////////////////////////////////////////////////////////

#include <iostream>
#include <algorithm>
#include <cmath>

#define PI 3.141592653

using namespace std;

struct Point {
    int x, y;
};

Point p[1000];
int n, L, stack[1000], top = 1;

double dis(Point &p1, Point &p2) {
    double delta_x = p1.x - p2.x;
    double delta_y = p1.y - p2.y;
    return sqrt(delta_x * delta_x + delta_y * delta_y);
}

int cross(Point &ps, Point &pe, Point &qs, Point &qe) {
    //二维向量叉乘公式：
    //cross<(x1, y1), (x2, y2)> = x1 * y2 - x2 * y1 = || * |b| * sin<a, b>
    //cross<p1, p2> = |p1| * |p2| * sin<p1, p2>
    //> 0 表示 p1 在 p2 的 顺时针方向， < 0 在逆时针方向， = 0 表示共线
    //函数所求为向量 ps->pe 和 qs->qe 的叉乘
    return (pe.x - ps.x) * (qe.y - qs.y) - (qe.x - qs.x) * (pe.y - ps.y);
}

bool cmp(Point &a, Point &b) {
    //比较到参考点极角大小，极角相等时比较距参考点的距离
    int res = cross(p[0], a, p[0], b);
    if (res == 0) {
        return dis(a, p[0]) < dis(b, p[0]);
    } else {
        return res > 0;
    }
}
void GrahamScan(void) {
    //顶点排序
    sort(p + 1, p + n, cmp);
    //前三个点入栈
    stack[0] = 0;
    stack[1] = 1;
    stack[2] = 2;
    top = 3;
    //扫描过程
    for (int i = 3; i < n; ++i) {
        while (cross(p[stack[top - 2]], p[stack[top - 1]], p[stack[top - 1]], p[i]) < 0) {
            --top;
        }
        stack[top++] = i;
    }
    stack[top++] = 0;
}

int main(void) {
    cin >> n >> L;
    double circle = 2 * L * PI;
    int ll = 0;
    for (int i = 0; i < n; ++i) {
        cin >> p[i].x >> p[i].y;
        //找出y坐标最小的中x坐标最小的点（左下角点，一定在凸包上）
        if (p[i].y < p[ll].y || p[i].y == p[ll].y && p[i].x < p[ll].x) {
            ll = i;
        }
    }
    //将左下角点放在p[0]处，作为参考点
    Point t = p[ll];
    p[ll] = p[0];
    p[0] = t;
    //运行核心Gramham扫描算法
    GrahamScan();
    //计算围墙最小周长
    double length = circle;
    while (top > 1) {
        length += dis(p[stack[top - 1]], p[stack[top - 2]]);
        --top;
    }
    //四舍五入
    cout << (int)(length + 0.5) << endl;
    return 0;
}

上面代码耗费236k、79ms，如果把cin, cout改为scanf和printf，可以降到164k、32ms.(有一次提交居然混到了0ms，POJ判时间随机性还真大==!)

处了上面的极角排序法，其实对点的扫描顺序其实还有另一种方法：水平序法。

水平序法即对所有点先按y升序排，y相等的按x升序排。然后对所有顶点扫描两次。第一次按排序后顺序扫描，这样可以得到凸包的右侧部分（右链），然后从接着再逆序扫描一次所有顶点，将把把凸包延长完整，即得到左侧部分（左链）。在扫描右链时从最左下的顶点开始，到最右上的顶点结束，而这两个点都一定会出现在凸包中，这样第一次扫描时将找到凸包的右侧，第二次扫描时从最右上扫描到最左下，这样就找到了左侧部分，连起来就是整个凸包了。

这样做的好处是：排序的比较函数简单，所以排序快，而且可以避免极角比较法中遇到共线问题的特殊处理。故在数据量大且共线情况多时，这种算法应该会优于极角比较法。

/////////////////////////////////////////////////////////////////
//        POJ1113 Wall    (Graham scan 算法, 水平序版)
//        Memory: 236K        Time: 79MS
//        Language: C++        Result : Accepted
/////////////////////////////////////////////////////////////////

#include <iostream>
#include <algorithm>
#include <cmath>

#define PI 3.141592653

using namespace std;

struct Point {
    int x, y;
};

Point p[1000];
int n, L, stack[1000], top = 1;

double dis(Point &p1, Point &p2) {
    double delta_x = p1.x - p2.x;
    double delta_y = p1.y - p2.y;
    return sqrt(delta_x * delta_x + delta_y * delta_y);
}

int cross(Point &ps, Point &pe, Point &qs, Point &qe) {
    //二维向量叉乘公式：
    //cross<(x1, y1), (x2, y2)> = x1 * y2 - x2 * y1 = || * |b| * sin<a, b>
    //cross<p1, p2> = |p1| * |p2| * sin<p1, p2>
    //> 0 表示 p1 在 p2 的 顺时针方向， < 0 在逆时针方向， = 0 表示共线
    //函数所求为向量 ps->pe 和 qs->qe 的叉乘
    return (pe.x - ps.x) * (qe.y - qs.y) - (qe.x - qs.x) * (pe.y - ps.y);
}

bool cmp(Point &a, Point &b) {
    //水平序，y小到大，y相等时x小到大
    if (a.y == b.y) {
        return a.x < b.x;
    }
    else {
        return a.y < b.y;
    }
}
void GrahamScan(void) {
    //顶点排序
    sort(p, p + n, cmp);
    //前两个点入栈
    stack[0] = 0;
    stack[1] = 1;
    top = 2;
    //扫描过程
    //第一次扫描形成凸包右半边（右链）
    for (int i = 2; i < n; ++i) {
        while (cross(p[stack[top - 2]], p[stack[top - 1]], p[stack[top - 1]], p[i]) < 0) {
            --top;
        }
        stack[top++] = i;
    }
    //第二次扫描形成凸包左半边（左链）
    for (int i = n - 2; i >= 0; --i) {
        while (cross(p[stack[top - 2]], p[stack[top - 1]], p[stack[top - 1]], p[i]) < 0) {
            --top;
        }
        stack[top++] = i;
    }
}

int main(void) {
    cin >> n >> L;
    double circle = 2 * L * PI;
    int ll = 0;
    for (int i = 0; i < n; ++i) {
        cin >> p[i].x >> p[i].y;
        ////找出y坐标最小的中x坐标最小的点（左下角点，一定在凸包上）
        //if (p[i].y < p[ll].y || p[i].y == p[ll].y && p[i].x < p[ll].x) {
        //    ll = i;
        //}
    }
    //将左下角点放在p[0]处，作为参考点
    /*Point t = p[ll];
    p[ll] = p[0];
    p[0] = t;*/
    //运行核心Gramham扫描算法
    GrahamScan();
    //计算围墙最小周长
    double length = circle;
    while (top > 1) {
        length += dis(p[stack[top - 1]], p[stack[top - 2]]);
        --top;
    }
    //四舍五入
    cout << (int)(length + 0.5) << endl;
    return 0;
}

---

接下来还有另一种经典的凸包求解算法：QuickHull算法。 这是一种分治算法，有些类似于QuickSort的思想。

　　假设现有一个向量，起点 A 和终点 B 都是凸包上的点，我们找到点集中在向量 AB “左侧”的所有点，并找出这些左侧点中距离 AB 最远的点 P ，（求法可参考上面向量叉乘中的介绍），那么 P 一定也是凸包上的点，现在我们就可以得到两个向量 AP 和 PB , 这样我们再递归调用QuickHull。最终就可以得到完整的凸包了。伪代码：

 void QuickHull(PointSet P, Vector S /*S.p,S.q 为端点*/){
 　　/* P 所有点都在 S 左侧*/
    找到 P 中距离 S 最远的 点 V;
    Vector A <- { S.p, V };
    Vector B <- { V, S.q };
    PointSet Q <- 在 P 中 且在 A 左侧的点;
    PointSet R <- 在 P 中 且在 B 左侧的点; /* 划分 */
    QuickHull(Q, A); /* 分治 */
    output(Point V); /* 按中序输出 保证顺序*/
    QuickHull(R, B); /* 分治 */
}

 这里有一个算法演示：Demo

 附上本题QuickHull版本的AC代码：

/////////////////////////////////////////////////////////////////
//        POJ1113 Wall    (QuickHull算法版)
//        Memory: 272K        Time: 32MS
//        Language: C++        Result : Accepted
/////////////////////////////////////////////////////////////////

#include <iostream>
#include <algorithm>
#include <cmath>

#define PI 3.141592653

using namespace std;

struct Point {
    int x, y;
};

Point p[1000], ans[1000];
int n, L, cnt = 0;
double d[1000];

//点点距
double dis(Point &p1, Point &p2) {
    double delta_x = p1.x - p2.x;
    double delta_y = p1.y - p2.y;
    return sqrt(delta_x * delta_x + delta_y * delta_y);
}

int cross(Point &ps, Point &pe, Point &qs, Point &qe) {
    return (pe.x - ps.x) * (qe.y - qs.y) - (qe.x - qs.x) * (pe.y - ps.y);
}

bool cmp(Point &a, Point &b) {
    return a.x == b.x ? a.y < b.y : a.x < b.x;
}


void QuickHull(int l, int r, Point a, Point b) {
    int k = l;
    for (int i = l; i <= r; ++i) {
        if (d[i] > d[k] || d[i] == d[k] && cmp(p[k], p[i])) {
            k = i;
        }
    }
    Point pt = p[k];    //距离最远的点
    int x = l - 1, y = r + 1;
    //找出左半段线段左边的点
    for (int k = l; k <= r; ++k) {
        d[++x] = cross(p[k], a, p[k], pt);
        if (d[x] > 0) {
            swap(p[k], p[x]);
        } else {
            --x;
        }
    }
    //找出右半段线段右边的点
    for (int k = r; k >= l; --k) {
        d[--y] = cross(p[k], pt, p[k], b);
        if (d[y] > 0) {
            swap(p[k], p[y]);
        } else {
            ++y;
        }
    }
    if (l <= x) {
        QuickHull(l, x, a, pt);
    }
    ans[cnt++] = pt;
    if (y <= r) {
        QuickHull(y, r, pt, b);
    }
}

int main(void) {
    cin >> n >> L;
    double circle = 2 * L * PI;
    cin >> p[0].x >> p[0].y;
    int minn = 0;
    for (int i = 1; i < n; ++i) {
        cin >> p[i].x >> p[i].y;
        if (cmp(p[i], p[minn])) {
            minn = i;
        }
    }
    //左下角点
    swap(p[0], p[minn]);
    //运行核心QuickHull算法
    ans[cnt++] = p[0];
    QuickHull(1, n - 1, p[0], p[0]);
    ans[cnt++] = p[0];
    //计算围墙最小周长
    double length = circle;
    while (cnt > 1) {
        length += dis(ans[cnt - 1], ans[cnt - 2]);
        --cnt;
    }
    //四舍五入
    cout << (int)(length + 0.5) << endl;
    return 0;
}

至今第一道用了三种算法写了四个版本代码的题T^T. 算是把二维情况的几种基本算法搞清楚了，可是实际要用的是求解高维凸包，还需要继续研究研究=。=.