初學算法 - 求凸包的Garham's Scan算法的C++實現

    所謂凸包，就是一個計算幾何（圖形學）中的概念。用不嚴謹的話來說，給定二維平面上的點集，凸包就是將最外層的點鏈接起來構成的凸多邊型，它能包含點集中全部的點。維基百科對集合X的凸包（Convex Hull）有四個定義，分別爲：

• The (unique) minimal convex set containing X            ---  包含集合X的最小凸集合

• The intersection of all convex sets containing X          --- 全部包含集合X的凸集合的交集

• The set of all convex combinations of points in X.       --- 集合X中全部點的凸組合的集合

• The union of all simplices with vertices in X.                --- 集合X中全部單一頂點的集合

    對於二維凸包，不如咱們把平面上的一些點想象爲「釘子」，而你正將一個橡皮筋撐的足夠大，以致於全部「釘子」都在你的橡皮筋包圍的區域裏。如今咱們鬆開它。「啪」的一聲，橡皮筋會盡量的收縮到極致，而這時撐起橡皮筋的這些「釘子」構成的集合， 也就是凸包。

                                            

    

    經過觀察，咱們能夠知道「最左」和「最右」的兩個點必定在構成凸包的集合裏。而Garham's Scan算法也正是注意到了這點。另外，若是咱們按照順時針方向觀察凸包，如P->Q->R，在每個點上凸包都是「右拐」的（固然，也可能構成一條直線）。

                                        

    使用兩個鏈表Lupper和Llower分別表示凸包的上半部分（Upper Hull）和下半部分（Lower Hull），Garham的算法能夠經過以下僞代碼描述：

Algorithm CONVEXHULL(P)
Input. A set P of points in the plane.
Output. A list containing the vertices of CH(P) in clockwise order.
 Sort the points by x-coordinate, resulting in a sequence p1,..., pn.
 Put the points p1 and p2 in a list Lupper, with p1 as the first point.
 for i ← 3 to n
     do Append pi to Lupper.
         while Lupper contains more than two points and the last three points in Lupper do not make a right turn
             do Delete the middle of the last three points from Lupper.
 Put the points pn and pn−1 in a list Llower 6 , with pn as the first point.
 for i ← n−2 downto 1
     do Append pi to Llower.
         while Llower contains more than 2 points and the last three points in Llower do not make a right turn
             do Delete the middle of the last three points from Llower.
 Remove the first and the last point from Llower to avoid duplication of the points where the upper and lower hull meet.
 Append Llower to Lupper, and call the resulting list L.
 return L

    

    如今問題轉到了已知三點A,B,C的座標，如何肯定C點是否在向量AB的左邊了。

    對於這個問題，咱們可使用幾何學的一個結論：有向面積。三角形的有向面積能夠經過一個行列式求得：

            (不知道這個行列式怎麼計算的童鞋能夠去補補線代了~)

    若是C在AB方向的左邊，那麼這個行列式求得的值是正的，反之爲負。在設計程序時，咱們只須要知道這個行列式值的正負，而並不關心它具體的值，所以就不須要再作一個除法了。

/**
 * Calculate the doubled directed area of pqs
 */
long long Area2( Point p, Point q, Point s)
{
    return p.x*q.y - p.y*q.x
          +q.x*s.y - q.y*s.x
          +s.x*p.y - p.x*s.y;
}

//Determinate whether s is on the left side of a directed line p->q
bool ToLeft( Point p, Point q, Point s)
{
    return Area2(p,q,s)>0;
}

    整個算法能夠寫成以下的代碼：

/**
 Created on 2015-10-12 11:24:04

 Description : Generate the Convex Hull via the Graham's Scan Algorithm
               Time Cost : O(nlogn)
               Thanks to Tsinghua Univ. MOOC Computational Geometry 2015
               As the answer of CG2015 PA1-1 Convex Hull
               (http://dsa.cs.tsinghua.edu.cn/oj/course.shtml?courseid=39)

 Author: ChenZheng / Arc001
 Email : chenzheng17@yahoo.com
 Copyright 2015 Xi'an University of Posts & Telecommunications
*/

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <list>

using namespace std;

enum Point_Status {isExtreme, notExtreme};

struct Point{
    long long x,y;
    int index;
    Point_Status status;

    bool operator < (const Point & rhs) const
    {
        if(x == rhs.x)
        {
            return y < rhs.y;
        }
        return x < rhs.x;
    }
};

long long Area2( Point p, Point q, Point s); //這裏再也不重複
bool ToLeft( Point p, Point q, Point s);

Point S[100005];
list<Point> upper_hull, lower_hull;

int main()
{
    int n;
    long long ans = 1;
    scanf("%d",&n);
    for(int i=1;i<=n;i++){
        scanf("%lld%lld",&S[i].x,&S[i].y);
        S[i].index = i;
        S[i].status = notExtreme;
    }

    sort(&S[1],&S[n+1]);

    upper_hull.push_back(S[1]);
    upper_hull.push_back(S[2]);

    list<Point>::iterator it,it2,it3;
    for(int i=3; i<=n; i++)
    {
        upper_hull.push_back(S[i]);
        it = it2 = it3 = upper_hull.end();
        --it; --it; --it;
        --it2; --it2;
        --it3;
        //cout<<it3->index<<" pushed!"<<endl;

        while(upper_hull.size()>2 && ToLeft(*it,*it2,*it3))
        {
            upper_hull.erase(it2);
            //cout<<it2->index<<" deleted! "<<it3->index<<" "<<Area2(*it,*it2,*it3)<<endl;
            it = it2 = it3 = upper_hull.end();
            --it; --it; --it;
            --it2; --it2;
            --it3;
        }
        //cout<<upper_hull.size()<<endl<<endl;
    }
    upper_hull.pop_back();

    lower_hull.push_back(S[n]);
    lower_hull.push_back(S[n-1]);

    for(int i=n-2; i>=1; i--)
    {
        lower_hull.push_back(S[i]);

        it = it2 = it3 = lower_hull.end();
        --it; --it; --it;
        --it2; --it2;
        --it3;
        while(lower_hull.size()>2 && ToLeft(*(it),*(it2),*(it3)))
        {
            lower_hull.erase(it2);
            it = it2 = it3 = lower_hull.end();
            --it; --it; --it;
            --it2; --it2;
            --it3;
        }
    }
    lower_hull.pop_back();

    //cout<<lower_hull.size()<<endl;
    int count_vertex = 0;
    for(it = upper_hull.begin(); it != upper_hull.end(); it++)
    {
        S[it->index].status = isExtreme;
        //cout<<it->index<<endl;
        count_vertex++;
    }

    for(it = lower_hull.begin(); it != lower_hull.end(); it++)
    {
        S[it->index].status = isExtreme;
        //cout<<it->index<<endl;
        count_vertex++;
    }

    for(int i=1; i<=n; i++){        //這裏在計算這道題的答案，與算法關係不大
        if(S[i].status == isExtreme)
        {
            ans *= i;
            ans %= (n+1);
            //cout<<i<<endl;
        }
    }

    ans *= count_vertex;
    ans %= (n+1);

    printf("%lld\n",ans);
    return 0;
}