计算几何--二维几何常用算法(多边形和凸包)

内容参考书籍——《算法竞赛入门经典训练指南》　　

　　在程序中，用顶点数组表示多边形，其中各个顶点按照逆时针顺序排列。

　　判断点是否在多边形内。采用转角法，基本思想是计算多边形相对于判定点转了多少度，具体来说，将多边形每条边的转角加起来，如果是360°，说明在多边形内；如果是0°，说明在多边形如果是180°则在多边形边界上。该方法在处理一些弧形多边形时丝毫不受影响，只需要每一段的终点到起点的转角累加起来即可。另外这个三角形甚至可以不是简单多边形（即可以自交）。

　　然而，直接计算会使用大量的反三角函数，不仅速度慢且容易产生精度误差。在算法竞赛中，我们并不会这样做，而是假想有一条向右的射线，统计多边形穿过这条射线正反多少次，把这个数记为绕数wn(Winding Number)，逆时针穿过时，wn加1，顺时针穿过时，wn减1。

　　注意程序实现时，判断是否穿过，以及穿过方向时，需要用叉积判断输入点在边的左边还是右边。

　　点在凸多边形内的判定更简单，只需要判断是否在所有边的左边(假设各顶点按照逆时针顺序排序)即可

　　凸包。

　　凸包就是把定点包围在内部的、面积最小的凸多边形。基于水平序的Andrew算法(比原始的Graham更快且更稳定)。首先把所有点按x从小到大排序(如果x相同，按照y从小到大排序)，删除重复点后得到序列p1,p2,...,然后把p1和p2放到凸包中。从p3开始，当新点在凸包“前进”方向的左边时继续，否则依次删除最近加入凸包的点，直到新点在左边。

　　如下图所示，新点P₁₈在向量P₁₀P₁₅(当前“前进方向“)的右边，因此需要从凸包上删除P₁₅和P₁₀，让P₈的下一个点为P₁₈。重复这个过程，直到碰到最右边的P_n，就求出了“下凸包”。然后反过来从P_n开始再做一次，求出“上凸包”，合并起来就是完整的凸包。

　　这个算法在排序后仅仅是从左到右和从右到左各扫描了一次，时间复杂度为O（n）。加上排序后时间复杂度也仅为O（nlogn）。

代码如下：

//判断该点与多边形关系
int isPointInPolygon(Point p, Polygon poly)
{
    int wn = 0;
    int n = v.size();
    for (int i = 0; i < n; ++i)
    {
        if (isPointOnSegment(p,poly[i], poly[(i+1)%n])) return -1;
        int k = dcmp(Cross(poly[(i+1)%n]-poly[i]),p-poly[i]);
        int d1 = dcmp(poly[i].y-p.y);
        int d2 = dcmp(poly[(i+1)%n].y-p.y);
        if (k>0 && d1<=0 && d2>0) wn++;
        if (k<0 && d2<=0 && d1>0) wn--;        
    }
    if (wn != 0)return 1;
    return 0;
}
//计算凸包，输入点数组p，个数为p，输出点数组ch。函数返回凸包顶点数。
//输入不能有重复点。函数执行完之后输入点的顺序被破坏。
//如果不希望在凸包的边上有输入点，把两个<=改成<
//在精度要求高时建议用dcmp比较
int ConvexHull(Point* p, int n, Point* ch)
{
    sort(p,p+n);
    int m = 0;
    for (int i = 0; i < n; ++i)
    {
        while(m>1 && Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]) <= 0) m--;
        ch[m++] = p[i];
    }
    int k = m;
    for (int i = n-2; i >= 0; --i)
    {
        while(m >k && Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]) <= 0) m--;
        ch[m++] = p[i];
    }
    if (n > 1) m--;
    return m;
}

 求多边形的重心

hdu 1115:给一个n多边形，求重心。

#include <bits/stdc++.h>
using namespace std;
struct Point
{
    double x,y;
    Point(double X = 0,double Y = 0){x=X,y=Y;}
    Point operator + (Point B) {return Point (x+B.x,y+B.y);}
    Point operator - (Point B) {return Point (x-B.x,y-B.y);}
    Point operator * (double k) {return Point (x*k,y*k);}
    Point operator / (double k) {return Point (x/k,y/k);}
};
typedef Point Vector;
double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}
double Polygon_area(Point *p,int n)     //求多边形有向面积，有正负不能取绝对值
{
    double area = 0;
    for (int i = 0; i < n; ++i)
    {
        area += Cross(p[i],p[(i+1)%n]);
    }
    return area/2;
}
Point Polygon_center(Point *p,int n) //求重心
{
    Point ans(0,0);
    if (Polygon_area(p,n) == 0) return ans;
    for (int i = 0; i < n; ++i)
    {
        ans = ans+(p[i]+p[(i+1)%n])*Cross(p[i],p[(i+1)%n]);
    }
    return ans/Polygon_area(p,n)/6;
}
int main(int argc, char const *argv[])
{
    int t,n,i;
    Point center;
    Point p[100000];
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d",&n);
        for (int i = 0; i < n; ++i)
        {
            scanf("%lf %lf",&p[i].x,&p[i].y);
        }
        center = Polygon_center(p,n);
        printf("%.2f %.2f\n",center.x,center.y);
    }
    return 0;
}

hdu 1392:输入n个点，求凸包的周长。

代码如下：

#include <bits/stdc++.h>
using namespace std;
const int maxn = 104;
const double eps = 1e-8;
int sgn(double x)         //判断x是否为0
{
    if (fabs(x)<eps) return 0;
    else return x<0? -1:1;
}
struct Point
{
    double x,y;
    Point(){}
    Point(double x, double y):x(x),y(y){}
    Point operator + (Point B) {return Point (x+B.x,y+B.y);}
    Point operator - (Point B) {return Point (x-B.x,y-B.y);}
    bool operator == (Point B) {return sgn(x-B.x) == 0 && sgn(y-B.y) == 0;}
    bool operator < (Point B)    //用于sort()排序
    {
        return sgn(x-B.x)<0 || (sgn(x-B.x) == 0 && sgn(y-B.y)<0);
    }
};
typedef Point Vector;
double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}//叉积
double Distance(Point A, Point B){return hypot(A.x-B.x,A.y-B.y);}
//求凸包。凸包顶点放在ch中，返回值是凸包的顶点数
int ConvexHull(Point* p, int n, Point* ch)
{
    sort(p,p+n);              //对点排序：按从大到小排序，如果x相同，按y排序
    n=unique(p,p+n)-p;        //去重复点
    int m = 0;
    //求下凸包，如果p[i]是右拐弯的，这个点不在凸包上，往回退
    for (int i = 0; i < n; ++i)
    {
        while(m>1 && Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]) <= 0) m--;
        ch[m++] = p[i];
    }
    int k = m;                //求上凸包
    for (int i = n-2; i >= 0; --i)
    {
        while(m >k && Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]) <= 0) m--;
        ch[m++] = p[i];
    }
    if (n > 1) m--;
    return m;                //返回凸包顶点数
}
int main(int argc, char const *argv[])
{
    int n;
    Point p[maxn],ch[maxn];
    while(scanf("%d",&n) && n)
    {
        for (int i = 0; i < n; ++i)
        {
            scanf("%lf%lf",&p[i].x,&p[i].y);
        }
        int v = ConvexHull(p,n,ch);//凸包顶点数
        double ans = 0;
        if (v==1) ans=0;
        else if(v==2)ans=Distance(ch[0],ch[1]);
        else
        {
            for (int i = 0; i < v; ++i) //计算凸包周长
            {
                ans+=Distance(ch[i],ch[(i+1)%v]);
            }
        }
        printf("%.2f\n",ans);
    }    
    return 0;
}