您的位置:首页 > 其它

【计算几何】 Andrew凸包算法 + 旋转卡壳(以求点对最长距离为例) -- 以 POJ 2187 Beauty Contest 为例

2012-09-06 09:07 501 查看
Part One 求凸包

凸包的算法有很多种,我这里使用的是是Andrew那个凸包算法,就是先循环所有点求出下凸包,再逆向循环所有点定求出上凸包的那个算法。

时间复杂度为O(NlogN)

一个java的Applet演示算法执行过程

Monotone Chain Convex Hull 算法伪代码:
输入:一个在平面上的点集P
点集 P 按 先x后y 的递增排序
声明两个栈:U,L
U - 上凸包
L - 下凸包
// 构建上凸包
for i = 1,2,3,4, .... n:
while U 中至少有两个点 and
U 中的栈首的两个点 与 P[i] 不为逆时针时:
U 出栈一次
P[i] 进入 U 栈
// 构建下凸包
for i=n,n-1,n-2, .... 1:
while L 中至少有两个点 and
L 中的栈首的两个点 与 P[i] 不为逆时针是:
L 出栈一次
P[i] 进入 L 栈

按顺序合并U L两个栈: U+L => H 则 集合H 为 点集P 的凸包点集合, 点的方向为逆时针

原文:
Input: a list P of points in the plane.
Sort the points of P by x-coordinate (in case of a tie, sort by y-coordinate).
Initialize U and L as empty lists.
The lists will hold the vertices of upper and lower hulls respectively.

for i = 1, 2, ..., n:
while L contains at least two points and the sequence of last two points
of U and the point P[i] does not make a counter-clockwise turn:
remove the last point from U
append P[i] to U

for i = n, n-1, ..., 1:
while L contains at least two points and the sequence of last two points
of L and the point P[i] does not make a counter-clockwise turn:
remove the last point from L
append P[i] to L

Remove the last point of each list (it's the same as the first point of the other list).
Concatenate L and U to obtain the convex hull of P.
Points in the result will be listed in counter-clockwise order.


LINK 1:Monotone Chain Convex Hull 英文原文算法伪代码来自这里www.algorithmist.com

LINK 2:
【精】各种凸包的演算法笔记+模板+分析

Part Tow 旋转卡壳(以求点对最长距离为例)



贡献两个比较好的连接(具体分析就不再赘述了,详见例题POJ 2187注释):

LINK 1:
详细的介绍了旋转卡壳算法的证明和各种应用(多达十几种)

LINK 2:
详细的介绍了旋转卡壳算法在求点对最长距离中应用(中文版)

// POJ 2187 Beauty Contest -- Andrew凸包算法 + 点对距离最长(旋转卡壳)
// 大意:给你一堆点,让你求出相隔最远的两点距离并输出
// 要求:最后结果不需要开平方 应为距离的平方
//
/*test data
5
1 3
2 1
4 1
5 3
3 5
== 17

*/

#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;

const double Pzero = 0.000001;
const int MAXP = 50005;

struct POINT{
int x,y;
}p[MAXP],CHp[MAXP*2];

double unsqrtDistance(POINT a, POINT b){
return (a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y);
}

double XJ(POINT a, POINT b, POINT p){
// p 关于 a->b 的左右位置 左为正,  a, b, p 为 逆时针顺序给出 则为正
// vector a->b (b.x-a.x , b.y-a.y)
// vector a->p (p.x-a.x , p.y-a.y)
// a->b X a->p == (b.x-a.x)*(p.y-a.y) - (p.x-a.x)*(b.y-a.y)

return ( (b.x-a.x)*(p.y-a.y) - (p.x-a.x)*(b.y-a.y) );
}

bool comp(POINT a, POINT b){
return (a.x < b.x) || ( (a.x == b.x) && (a.y < b.y));
}

int AndrewConvexHull(int n){
// 所得 数组CHp[] 为凸包的逆时针序列

int nCHp = 0; // 用于CHp[]栈顶指针, 实际意义为凸包顶点数

// 升序排序 先x后y
sort(p , p + n , comp);

// 上半凸包 p从前往后循环
for (int i = 0; i < n; i++){
while(nCHp >= 2 && XJ(CHp[nCHp-2], CHp[nCHp-1], p[i]) >= 0)
nCHp--; // 栈首两个点和当前的p点 呈顺时针了 则CHp[]出栈一个, 否则入栈一个, 下同
CHp[nCHp++] = p[i];
}

// 下半凸包 p从后往前循环 (p[0]被计算过两次)
// 把上面的处理后的末端点留作这次循环的栈底,因此 m = nCHp + 1
for (int i = n - 2, m = nCHp + 1; i >= 0; i--){
while(nCHp >= m && XJ(CHp[nCHp-2], CHp[nCHp-1], p[i]) >= 0)
nCHp--;
CHp[nCHp++] = p[i];
}
// 起点被计算了两遍, 因此求下半凸包中加进去的最后一个值p[0]要去掉
return --nCHp;
}

double rotating_calipers(POINT *ch, int n){
// 旋转卡壳 求点对最长距离     效率 O(n)
// ch[] 凸包 为逆时针给出的
// 计算最长距离的时候,枚举了每一条边 |p p+1| 及 离该边最远的点q 的情况, 每一次计算点到边两端的距离,并更新一下ans(取最大值)
// 如果遇到下一个点与该线段围成的三角形面积(即叉积) 比这一次变大了,则比较下一个点 (即q++)
// 如果遇到下一个点与该线段围成的三角形面积(即叉积) 比这一次变小了,则计算两个距离(因为会出现平行的情况)并更新ans
// POJ 2187 暴力可过 肯定凸包的顶点数不多

int q = 1;double ans = 0;
ch
= ch[0]; // ch[q+1] 可能会访问到
for(int p = 0; p < n; p++){
// 此处遍历的线段为 |p->p+1| 点为 q
// 所以叉积的意义为 线段|p->p+1| 和 点q 围成的三角形的有向面积,逆时针为正,因此在这里均为正
while(XJ(ch[p+1], ch[q], ch[p]) > XJ(ch[p+1], ch[q+1], ch[p]))
q = (q + 1) % n; // 移到下一个点,因为q可能会转几圈,因此对n取余

// 这里要比较q点到线段 |p->p+1| 两个端点的距离,因为有可能在旋转的时候两个边平行
ans = max( ans, max(unsqrtDistance(ch[p], ch[q]), unsqrtDistance(ch[p+1], ch[q+1])) );
}
return ans;
}

void InitInput(int n){
for(int i=0;i<n;i++){
scanf("%d%d", &p[i].x, &p[i].y);
}
}

int main()
{
//	freopen("in.txt","r",stdin);
int n;
while( ~scanf("%d",&n) ){
InitInput(n);
int nCHp = AndrewConvexHull(n); // 返回凸包的顶点数
printf("%.0lf\n", rotating_calipers(CHp, nCHp) );
}
return 0;
}


http://cgm.cs.mcgill.ca/~orm/width.html 凸包的宽度(英文版)
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: