二维点集的凸包及其直径(1)

关注

最新推荐文章于 2026-08-01 14:57:04 发布 861 阅读 AI 写同款· GEO 优化 ›

前言 :因为前几天做了一个有关凸包的题,并答应crackerwang写个blog解释一下我的算法.因为我比较懒的原因,一直拖到现在才写.预计一共有两篇,第一篇介绍求二维点集凸包的O(N*logN)时间复杂度的算法.第二篇介绍求凸包直径的O(N)时间复杂度的算法.
下面首先给出http://acm.tju.edu.cn/toj/showp2847.html 该题的C++代码,本文将使用Java 代码来描述.

cpp 代码

/**
* TJU 2847 的C++解法
* Written By Eastsun
*/
#include＜stdio.h＞
#include＜math.h＞
#include＜algorithm＞
#include＜functional＞
#define S(arr,a,b,c) ((arr[b].x-arr[a].x)*(arr[c].y-arr[a].y)-(arr[c].x-arr[a].x)*(arr[b].y-arr[a].y))
#define J(arr,a,b,c,d) ((arr[b].x-arr[a].x)*(arr[d].y-arr[c].y)-(arr[d].x-arr[c].x)*(arr[b].y-arr[a].y))
#define Q(x) ((x)*(x))
#define D(arr,a,b) (Q(arr[a].x-arr[b].x)+Q(arr[a].y-arr[b].y))
using namespace std;
struct Point{
double x,y;
};
struct myCmp: public binary_function＜Point,Point, bool ＞{
bool operator()( const Point& p1, const Point& p2) const {
if (p1.x==p2.x) return p1.y＜=p2.y;
else return p1.x＜p2.x;
}
};
int main(){
int n;
while (scanf( "%d" ,&n),n){
Point *ps = new Point[n];
for ( int i=0;i＜n;i++) scanf( "%lf%lf" ,&ps[i].x,&ps[i].y);
sort(ps,ps+n,myCmp());
int *v = new int [2*n];
int p,q;
p =q =n;
for ( int i=0;i＜n;i++){
v[p] =v[q] =i;
while (p-n＞=2&&S(ps,v[p],v[p-1],v[p-2])＞=0){v[p-1] =v[p]; p--;}
while (n-q＞=2&&S(ps,v[q+2],v[q+1],v[q])＞=0){v[q+1] =v[q]; q++;}
p++;
q--;
}
int len =p -q -2;
Point *pps = new Point[len+2];
for ( int i=q+1;i＜p;i++) pps[i-q] =ps[v[i]];
pps[0] =pps[len];
int i=0,j=1;
while (J(pps,i,i+1,j,j+1)＞0) j++;
double max =0;
while (i＜=len){
double det =J(pps,i,i+1,j,j+1);
if (det＜0) i++;
else if (det＞0) j =(j+1＞len?j+1-len:j+1);
else {
i++;
continue ;
}
double tmp =D(pps,i,j);
if (tmp＞max) max =tmp;
}
delete [] ps;
delete [] pps;
delete [] v;
printf( "%.2lf\n" ,sqrt(max));
}
return 0;
}

一:凸包的定义及其应用
对于二维点集,其凸包(convex hull)是指包含所有点的最小的凸多边形.直观上看,如果用一条橡皮筋将所有的点圈住,当橡皮筋拉紧后的形状就是这些点的凸包. 下面就是凸包的示意图:

convex hull

       那么,我们为什么需要凸包这个概念呢?它又能解决什么问题?
      首先,凸包上的点相对原有的点集,我们可以想象,其数量将大大减少.研究表明,对于二维情形,凸包顶点数m(k) =O(n^1/3).更一般的,对于k维球体中均匀分布的n个点,其凸包顶点数m(k) =O(n^(k-1)/(k+1)),可见凸包可大大降低平均意义下的时空复杂度.
     另一方面,凸包相对原有点集增加了一个"序",原来是一个杂乱无章的点集,而现在是一个性质优美的凸多边形,研究起来方便很多.

      关于其应用,这里我们只针对TJU2847来说,原题是求一个点集中任意两点的距离的最大值.显然,如果我们直接考虑这个点集中任意两点的距离,时间复杂度是O(N^2),下面我们可以看到,当求出其凸包后,我们可以在O(k)时间内求出这个值(这里的k是指凸包顶点个数k<=N).再加上构造凸包的时间复杂度O(N*logN),我们在O(N*logN)时间复杂度内解决了这个问题.

二:构造凸包的算法
      二维点集构造凸包的算法有:卷包裹法(Gift-Wrapping),Graham-Scan扫描法以及分治算法,增量算法等等.这里我采用的是Graham-Scan算法的一种变形--x-y排序.(至于原本的极角排序的Graham-Scan算法,有兴趣的可以参阅《Introduction to Algorithms》第VII章的Computational Geometry一节,里面有详细的说明及正确性证明):
      首先将给定点集对x的值进行排序,x值相同的按y的值排序.简单来说就是从左到右,从下到上排序.
      排序后,记这些点为ps[0,1,..,k],显然点列的第一个点与最后一个点都在凸包上(想一想那个橡皮筋),然后我们从第一点开始到最后一个点进行如下处理:
      构造一个堆栈(数组)stack[],stack[0] =ps[0],然后维持栈中点都是按右手系旋转(或者说从stack[0]到stack[p],所有的点都是按逆时针排列),对于每一个新增加的点,都首先检查栈顶的两个点与其是不是保持右手系,如果是,把这个点加入栈;否则,将栈顶点去掉,继续检查...一直到符合要求为止.

这样处理后的结果是:stack[]中得到一条连接ps[0]与ps[k]的,并且维持逆时针顺转的链.这个链就是我们要求的凸包的下半部分.用伪代码来描叙就是:

java 代码

//ps中保存了已经排好序的点
Point[] stack = new Point[ps.length];
int index = 0 , p = 0 ;
while (index
stack[p] =ps[index++];
while (p>= 2 && stack[p- 2 ],stack[p- 1 ],stack[p] 不是右手系){
stack[p- 1 ] =stack[p];
p --;
}
p ++;
}

显然,我们再从ps[k]到ps[0]做类似的处理就可以得到凸包的上半部分.当然,事实上我们可以把这两个工作一起做了:维持一个双头栈,使其头部的三个点为右手系,其尾部的三个点为左手系.这样经过一次扫描就可以得到整个凸包.伪代码如下:

java 代码

//ps中保存了已经排好序的点
Point[] stack = new Point[ 2 *ps.length];
int index = 0 , head =ps.length,tail =ps.length;
while (index
stack[head] =stack[tail] =ps[index++];
while (head-ps.length>= 2 && stack[head- 2 ],stack[head- 1 ],stack[head] 不是右手系){
stack[head- 1 ] =stack[head];
head --;
}
while (ps.length-tail>= 2 && stack[tail+ 2 ],stack[tail+ 1 ],stack[tail] 不是左手系){
stack[tail+ 1 ] =stack[tail];
tail ++;
}
head ++;
tail --;
}

可以看出,整个算法并不复杂,或者可以说很优美.哦,等等,伪代码里面还有个判断"左手系","右手系"是怎么来的?这就涉及到一个数学概念:叉积

三:叉积
    叉积,又叫外积,向量积.这里我不对叉积做太深入的探讨,只做一些概念性的描叙(且没有考虑数学上的准确性,只是为了解决问题方便),有兴趣的可以找本大学的<解析几何>之类的数学书看看.
    定义:对于二维平面的两个向量a =[xa,ya],b =[xb,yb],其叉积 a×b =xa*yb -xb*ya (其实叉积的定义是在三维空间中的..)
    另一方面,对于点X1 =[x1,y1],X2 =[x2,y2],对应的向量 X1X2 =[x2-x1,y2-y1]
    性质:
    对于平面上的三个点A,B,C所组成的三角形ABC的面积是向量AB与AC叉积的绝对值的一半,即 2*S(A,B,C) =|AB×AC|,并且当A,B,C三点按逆时间顺转时,AB×CD为正;当A,B,C三点按顺时间顺转时,AB×AC为负;A,B,C三点共线时,其叉积为0.

四:具体代码
有了上面的知识,下面给出具体的求凸包的Java代码.
首先给出用于表示Point的类:

java 代码

package convex;
public class Point implements Comparable＜Point＞{
public double x,y;
public static double distanceSq(Point p1,Point p2){
return (p2.x-p1.x)*(p2.x-p1.x)+(p2.y-p1.y)*(p2.y-p1.y);
}
public Point(){
this ( 0.0 , 0.0 );
}
public Point( double x, double y){
this .x =x;
this .y =y;
}
/**
* 实现Comparable的compareTo方法,将Point按从左到右,从下到上排序
*/
public int compareTo(Point p){
if (x ==p.x) return y==p.y? 0 :(y>p.y? 1 :- 1 );
else return x>p.x? 1 :- 1 ;
}
}

下面给出求凸包的代码:

java 代码

package convex;
import static java.util.Arrays.sort;
/**
* 用于求点集凸包以及求凸包直径的算法类
*/
public class Algorithm{
//避免生成该类实例
private Algorithm(){}
/**
* 构造Point数组ps中从fromIndex到toIndex的Point的凸包,时间复杂度O(N*logN)
* 注意: 该方法可能改变ps数组中从fromIndex到toIndex元素的顺序
* 前置条件: 一个任意的二维点集,不同点的个数>=2
*@param ps 保存二维点集的Point数组
*@param fromIndex 点集在数组中的起始位置
*@param toIndex 点集在数组中的结束位置
*@param convex 用来保存凸多边形的顶点的Point数组,其顶点将按逆时针排列
*@param offset 保存数据在convex数组中的起始位置
*@return length 凸多边形的顶点数目
*/
public static int getPointsConvexClosure(Point[] ps, int fromIndex, int toIndex,Point[] convex, int offset){
sort(ps,fromIndex,toIndex);
int len =toIndex -fromIndex;
Point[] tmp = new Point[ 2 *len];
int up =len, down =len;
for ( int index =fromIndex;index
tmp[up] =tmp[down] =ps[index];
while (len-up>= 2 &&multiply(tmp[up+ 2 ],tmp[up+ 1 ],tmp[up])>= 0 ){
tmp[up+ 1 ] =tmp[up];
up++;
}
while (down-len>= 2 &&multiply(tmp[down- 2 ],tmp[down- 1 ],tmp[down])<= 0 ){
tmp[down- 1 ] =tmp[down];
down--;
}
up --;
down ++;
}
System.arraycopy(tmp,up+ 1 ,convex,offset,down-up- 2 );
return down-up- 2 ;
}
    /**
    * 计算向量ab与ac的叉积
    */
    private static double multiply(Point a,Point b,Point c){
        return multiply(a,b,a,c);
    }
    /**
    * 计算向量ab与cd的叉积
    */
    private static double multiply(Point a,Point b,Point c,Point d){
        return (b.x-a.x)*(d.y-c.y)-(d.x-c.x)*(b.y-a.y);
    }
}