/********************************************* *                                          *
*        計算幾何學庫函數                  *
*                                          *
*                  copyright starfish      *
*                          2000/10/24      *
*                                          *
\*********************************************/

#include <stdlib.h>


#define infinity 1e20
#define EP 1e-10

struct TPoint {
    float x,y;
    };

struct TLineSeg {
    TPoint a,b;
    };



//求平面上兩點之間的距離
float distance(TPoint p1,TPoint p2)
{
    return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));  
}


/******************************************************** *                                                      *
*  返回(P1-P0)*(P2-P0)的叉積。                        *
*  若結果為正，則<P0,P1>在<P0,P2>的順時針方向；        *
*  若為0則<P0,P1><P0,P2>共線；                        *
*  若為負則<P0,P1>在<P0,P2>的在逆時針方向;            *
*  可以根據這個函數確定兩條線段在交點處的轉向,        *
*  比如確定p0p1和p1p2在p1處是左轉還是右轉，只要求      *
*  (p2-p0)*(p1-p0)，若<0則左轉，>0則右轉，=0則共線    *
*                                                      *
\********************************************************/

float multiply(TPoint p1,TPoint p2,TPoint p0)
{
    return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));    
}



//確定兩條線段是否相交
int intersect(TLineSeg u,TLineSeg v)
{
    return( (max(u.a.x,u.b.x)>=min(v.a.x,v.b.x))&&
            (max(v.a.x,v.b.x)>=min(u.a.x,u.b.x))&&
            (max(u.a.y,u.b.y)>=min(v.a.y,v.b.y))&&
            (max(v.a.y,v.b.y)>=min(u.a.y,u.b.y))&&
            (multiply(v.a,u.b,u.a)*multiply(u.b,v.b,u.a)>=0)&&
            (multiply(u.a,v.b,v.a)*multiply(v.b,u.b,v.a)>=0));
}

//判斷點p是否在線段l上
int online(TLineSeg l,TPoint p)
{
    return( (multiply(l.b,p,l.a)==0)&&( ((p.x-l.a.x)*(p.x-l.b.x)<0 )||( (p.y-l.a.y)*(p.y-l.b.y)<0 )) );
}


//判斷兩個點是否相等
int Euqal_Point(TPoint p1,TPoint p2)
{
return((abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP));
}

//一種線段相交判斷函數，當且僅當u,v相交并且交點不是u,v的端點時函數為true;
int intersect_A(TLineSeg u,TLineSeg v)
{
    return((intersect(u,v))&&
          (!Euqal_Point(u.a,v.a))&&
          (!Euqal_Point(u.a,v.b))&&
          (!Euqal_Point(u.b,v.a))&&
          (!Euqal_Point(u.b,v.b)));
}


/************************************************* *                                              *
*      判斷點q是否在多邊形Polygon內，          *
*      其中多邊形是任意的凸或凹多邊形，        *
*      Polygon中存放多邊形的逆時針頂點序列      *
*                                              *
\*************************************************/

int InsidePolygon(int vcount,TPoint Polygon[],TPoint q)
{
    int c=0,i,n;
    TLineSeg l1,l2;
        
    l1.a=q;
    l1.b=q;
    l1.b.x=infinity;
    n=vcount;
    for (i=0;i<vcount;i++)
    {
        l2.a=Polygon[i];
        l2.b=Polygon[(i+1)%n];
        if (  (intersect_A(l1,l2))||
              (
                (online(l1,Polygon[(i+1)%n]))&&
                (
                (!online(l1,Polygon[(i+2)%n]))&&
                (multiply(Polygon[i],Polygon[(i+1)%n],l1.a)*multiply(Polygon[(i+1)%n],Polygon[(i+2)%n],l1.a)>0)
                ||
                (online(l1,Polygon[(i+2)%n]))&&
                (multiply(Polygon[i],Polygon[(i+2)%n],l1.a)*multiply(Polygon[(i+2)%n],Polygon[(i+3)%n],l1.a)>0)        
                )
              )
          ) c++;
        }
        return(c%2!=0);
}


/******************************************** *                                          *
*      計算多邊形的面積                    *
*                                          *
*    要求按照逆時針方向輸入多邊形頂點    *
*    可以是凸多邊形或凹多邊形            *
*                                          *
\********************************************/

float area_of_polygon(int vcount,float x[],float y[])
{
  int i;
  float s;
  if (vcount<3) return 0;
  s=y[0]*(x[vcount-1]-x[1]);
  for (i=1;i<vcount;i++)
    s+=y[i]*(x[(i-1)]-x[(i+1)%vcount]);
  return s/2;
  }


/*====================================================================================
  尋找凸包 graham 掃描法

  PointSet為輸入的點集；
  ch為輸出的凸包上的點集，按照逆時針方向排列;
  n為PointSet中的點的數目
  len為輸出的凸包上的點的個數；

      設下一個掃描的點PointSet[i]=P2,當前棧頂的兩個點ch[top]=P1,ch[top-1]=P0,
    如果P1P2相對于P0P1在點P1向左旋轉(共線也不行)，則P0,P1一定是凸包上的點；
    否則P1一定不是凸包上的點，應該將其出棧。
    比如下面左圖中點1就一定是凸包上的點，右圖中點1就一定不是凸包上的點，因為
    可以連接點0,2構成凸包的邊
          2
          |
          |                _____2
          1              1
          /              /
    ____0/          ____0/

====================================================================================*/

void Graham_scan(TPoint PointSet[],TPoint ch[],int n,int &len)
{
    int i,j,k=0,top=2;
    TPoint tmp;
    
    //選取PointSet中y坐標最小的點PointSet[k]，如果這樣的點右多個，則取最左邊的一個
    for(i=1;i<n;i++)
        if ((PointSet[i].y<PointSet[k].y)||((PointSet[i].y==PointSet[k].y)&&(PointSet[i].x<PointSet[k].x)))
        k=i;
    tmp=PointSet[0];
    PointSet[0]=PointSet[k];
    PointSet[k]=tmp; //現在PointSet中y坐標最小的點在PointSet[0]
    for (i=1;i<n-1;i++) //對頂點按照相對PointSet[0]的極角從小到大進行排序，極角相同的按照距離PointSet[0]從近到遠進行排序
        {  k=i;
            for (j=i+1;j<n;j++)
                if ( (multiply(PointSet[j],PointSet[k],PointSet[0])>0)
                    ||
                    ( (multiply(PointSet[j],PointSet[k],PointSet[0])==0)
                      &&(distance(PointSet[0],PointSet[j])<distance(PointSet[0],PointSet[k])) )
                  ) k=j;
            tmp=PointSet[i];
            PointSet[i]=PointSet[k];
            PointSet[k]=tmp;
        }
    ch[0]=PointSet[0];
    ch[1]=PointSet[1];
    ch[2]=PointSet[2];  
    for (i=3;i<n;i++)
        {
            while (multiply(PointSet[i],ch[top],ch[top-1])>=0) top--;
            ch[++top]=PointSet[i];
            }
    len=top+1;
}