const
maxn
=
50000
;
type
point
=
record
x,y:longint;
end
;
var
n,i,x,m,ans,j:longint;
ch,p:
array
[
1
..maxn
+
1
]
of
point;
s:
array
[
1
..maxn]
of
longint;
function
cross(a,b,c:point):longint; inline;
begin
cross:
=
(b.x
-
a.x)
*
(c.y
-
a.y)
-
(b.y
-
a.y)
*
(c.x
-
a.x);
end
;
function
dist(a,b:point):longint; inline;
begin
dist:
=
sqr(a.x
-
b.x)
+
sqr(a.y
-
b.y);
end
;
function
cmp(a,b:point):boolean; inline;
begin
cmp:
=
(a.x
<
b.x)
or
(a.x
=
b.x)
and
(a.y
<
b.y);
end
;
function
max(a,b:longint):longint;
begin
if
a
>
b
then
max:
=
a
else
max:
=
b;
end
;
procedure
swap(a,b:longint); inline;
var
x:point;
begin
x:
=
p[a]; p[a]:
=
p[b]; p[b]:
=
x;
end
;
procedure
hull(l,r:longint; a,b:point);
var
x,i,j,k:longint;
y:point;
begin
x:
=
l; y:
=
p[l];
for
k:
=
l
to
r
do
if
(s[x]
<
s[k])
or
(s[x]
=
s[k])
and
(cmp(y,p[k]))
then
begin
x:
=
k; y:
=
p[k];
end
;
i:
=
l
-
1
; j:
=
r
+
1
;
for
k:
=
l
to
r
do
begin
inc(i); s[i]:
=
cross(p[k],a,y);
if
s[i]
>
0
then
swap(i,k)
else
dec(i);
end
;
for
k:
=
r
downto
l
do
begin
dec(j); s[j]:
=
cross(p[k],y,b);
if
s[j]
>
0
then
swap(j,k)
else
inc(j);
end
;
if
l
<=
i
then
hull(l,i,a,y);
inc(m); ch[m]:
=
y;
if
j
<=
r
then
hull(j,r,y,b);
end
;
begin
assign(input,
'
Maxd.in
'
); reset(input);
assign(output,
'
Maxd.out
'
); rewrite(output);
readln(n);
for
i:
=
1
to
n
do
begin
readln(p[i].x,p[i].y);
if
(x
=
0
)
or
cmp(p[i],p[x])
then
x:
=
i;
end
;
swap(
1
,x);
m:
=
1
; ch[
1
]:
=
p[
1
]; hull(
2
,n,p[
1
],p[
1
]);
ch[m
+
1
]:
=
ch[
1
]; j:
=
2
; ans:
=
0
;
for
i:
=
1
to
m
do
begin
while
cross(ch[i],ch[j],ch[i
+
1
])
<
cross(ch[i],ch[j
+
1
],ch[i
+
1
])
do
begin
inc(j);
if
j
>
m
then
j:
=
1
;
end
;
ans:
=
max(ans,dist(ch[i],ch[j]));
end
;
writeln(ans);
close(input); close(output);
end
.
转载于 http://www.cnblogs.com/Booble/archive/2011/04/03/2004865.html
一.简单枚举算法的不足
上一次介绍了一个基本的求平面最远点对的算法
即先求点集的凸包 然后枚举凸包上的点来求最远点集
这是利用了凸包上的点相比 点集中的点 一般是很少的 平均情况很好 并且我们也能AC这个问题
但是这是有局限性的 当凸包上的点达到O(N)的级别时 凸包的优化作用就不存在了
不过我们还要考虑到 凸包还起了对凸包上点集排序的作用
凸包有很多的优美的性质 我们可以加以利用 以得到更加高效的算法
旋转卡壳算法就是利用凸包特性的一类解决问题的方法
==============================
二.旋转卡壳算法
旋转卡(qiǎ)壳算法(Rotating Calipers Algorithm):
是解决一些与凸包有关问题的有效算法 就像一对卡壳卡住凸包旋转而得名
Every time one blade of the caliper lies flat against an edge of the polygon, it forms an antipodal pair with the point or edge touching the opposite blade. It turns out that the complete "rotation" of the caliper around the polygon detects all antipodal pairs and may be carried out in O(n) time.
http://en.wikipedia.org/wiki/Rotating_calipers
(图片来自:http://cgm.cs.mcgill.ca/~orm/rotcal.html)
被一对卡壳正好卡住的对应点对称为对踵点(Antipodal point)
http://en.wikipedia.org/wiki/Antipodal_point
可以证明对踵点的个数不超过3N/2个 也就是说对踵点的个数是O(N)的
对踵点的个数也是我们下面解决问题时间复杂度的保证
上第一个图是卡壳的一般情况 卡住两点 图二是卡住一条边和一个点
由于实现中 卡住两点的情况不好处理 我们通常关注第二种情况
在第二种情况中 我们可以看到 一个对踵点和对应边之间的距离比其他点要大
也就是一个对踵点和对应边所形成的三角形是最大的 下面我们会据此得到对踵点的简化求法
看一下官方的伪代码:
当时我看完了 就一个字 长... 我最讨厌冗长的程序了...
begin p0: = pn; q: = NEXT[p]; while (Area(p,NEXT[p],NEXT[q]) > Area(p,NEXT[p],q)) do q: = NEXT[q]; q0: = q; while (q ! = p0) do begin p: = NEXT[p]; Print(p,q); while (Area(p,NEXT[p],NEXT[q]) > Area(p,NEXT[p],q) do begin q: = NEXT[q]; if ((p,q) ! = (q0,p0)) then Print(p,q) else return end ; if (Area(p,NEXT[p],NEXT[q]) = Area(p,NEXT[p],q)) then if ((p,q) ! = (q0,p0)) then Print(p,NEXT[q]) else Print(NEXT[p],q) end end .
几经折腾 终于找到了一个不错的实现:http://www.cnblogs.com/DreamUp/archive/2010/09/16/1828131.html
不过不是很好理解 这里作一下说明
1 ch[m + 1 ]: = ch[ 1 ]; j: = 2 ; 2 for i: = 1 to m do 3 begin 4 while cross(ch[i],ch[j],ch[i + 1 ]) < cross(ch[i],ch[j + 1 ],ch[i + 1 ]) do 5 begin inc(j); if j > m then j: = 1 ; end ; 6 writeln(ch[i].x,' ',ch[i].y, ' ' ,ch[j].x,' ',ch[j].y); 7 end ;
上面就是旋转卡壳寻找对踵点的过程
其中叉积函数Cross(A,B,C:Point):Real 返回AB到AC的二维定义下的叉积
这里主要用到了叉积求三角形面积的功能
我们对于一条对应边<CH i,CH Next[i]>求出距离这条边最远的点CHj
则由上面第二种情况可知 CH i 和 CH j 为一对对踵点 这样让 CH i 绕行凸包一周即可得到所有的对踵点
下面面这个图 由于本人的gif图制作水平拙劣 所以不好看
需要的可以下载几何画板察看原版GSP文件 点击这里下载GSP文件
接下来考虑 如何得到距离每条对应边的的最远点呢?
稍加分析 我们可以发现 凸包上的点依次与对应边产生的距离成单峰函数
具体证明可以从凸包定义入手 用反证法解决
这样我们再找到一个点 使下一个点的距离小于当前的点时就可以停止了
而且随着对应边的旋转 最远点也只会顺着这个方向旋转 我们可以从上一次的对踵点开始继续寻找这一次的
由于内层while循环的执行次数取决于j增加次数 j最多增加O(N)次
所以求出所有对踵点的时间复杂度为O(N)
还有有两点需要注意:
1.上面这段代码及代码的分析都是需要凸包上没有三点共线的
2.Next[i] 不需要手动求 在原代码中有很好的处理
最后指出网上很多文章的一个错误 一个点的对踵点并不是离这个点最远的点!
这样子的点对是根本不满足对踵点的性质的 即最为重要的单峰分布性质
下图是一个反例:
==============================
三.旋转卡壳算法的简单应用
至此我们终于可以更高效的解决平面最远点对问题了
有一个很重要的结论是 最远点对必然属于对踵点对集合
那么我们先求出凸包 然后求出对踵点对集合 然后选出距离最大的即可
用这个算法可以47ms AC这个问题 算上凸包的时间 总复杂度为O(Nlog2N)
[Poj 2187]代码如下:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define _DEBUG 1
#define MAX_N 50010
typedef struct{
int x,y;
}Point;
Point pts[MAX_N];//postion of farms
Point chs[MAX_N];//凸包点
int n;//个数
long dist(Point p0,Point p1){
return ((p1.x-p0.x)*(p1.x-p0.x)+(p1.y-p0.y)*(p1.y-p0.y));
}
long crossMultiple(Point p0,Point p1,Point p2){//计算（p1-p0)×（p2-p0)
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}
void exchange(Point &p1,Point &p2){
Point tmp = p1;
p1 = p2;
p2 = tmp;
}
int compare(const void *p1,const void *p2){//可能次序反了
Point p3 = *(Point *)p1;
Point p4 = *(Point *)p2;
long cm = crossMultiple(pts[0],p3,p4);
if(cm < 0)//方向注意
return 1;
else if(cm == 0){
if(dist(p3,pts[0]) > dist(p4,pts[0]) )
return 1;
}
return -1;
}
int convexHull(){
int i;
int temp = 0;
for(i=1;i<n;++i){//找到最下，最左边的点
if(pts[i].y < pts[temp].y ||
(pts[i].y == pts[temp].y && pts[i].x < pts[temp].x))
temp = i;
}
exchange(pts[0],pts[temp]);//交换temp 与第0个点
qsort(pts+1,n-1,sizeof(Point),compare);//点按逆时针排序
chs[0] = pts[0];
chs[1] = pts[1];
int top = 1;
for(i=2;i<n;++i){
while(crossMultiple(chs[top-1],chs[top],pts[i]) <= 0 ){//出栈,这边可能有问题，边界判断
if(top <1)break;
top--;
}
chs[++top] = pts[i];
}
return top+1;
}
int main(){
int i,j;
#if _DEBUG == 1
freopen("POJ2187.in","r",stdin);
#endif
scanf("%d",&n);
for(i=0;i<n;++i){
scanf("%d %d",&pts[i].x,&pts[i].y);
}
int m=convexHull();//计算凸包,并返回凸包中点个数
j = 1;
long distance = 0;
long d;
chs[m] = chs[0];//注意
for(i=0;i<m;++i){//通过旋转卡壳算法得到最大距离,要考虑的问题：假如凸包上的点在一条直线上
while(crossMultiple(chs[(j+1) % m],chs[i],chs[i+1]) > crossMultiple(chs[j],chs[i],chs[i+1])){
j = (j+1) % m;
}
if((d=dist(chs[j],chs[i])) > distance)//可能有问题
distance = d;
if((d=dist(chs[j],chs[i+1])) > distance)
distance = d;
}
printf("%ld\n",distance);
return 0;
}
1361
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
