poj 1279 半平面交求多边形的核，求核的面积

http://acm.pku.edu.cn/JudgeOnline/problem?id=1279

这个问题是所谓多边形的核，在一个多边形（不一定是凸的），要求找到一个区域，使得在这个区域可见多边形的所有边界。

#include <cmath>
#include <stdio.h>
#include <string.h>
#include <algorithm>
#include <iostream>
using namespace std;

typedef double TYPE;
#define MaxPoint 1550
#define Epsilon  1e-10     /*验证*///精度的范围 ,根据不同的情况调整精度值 
#define Abs(x)   (((x)>0)?(x):(-(x)))   /*验证*/
//空间中的点，可以用来作为二维点来用 
struct POINT {/*验证*/
    TYPE x; TYPE y; TYPE z;
    POINT() : x(0), y(0), z(0) {};
    POINT(TYPE _x_, TYPE _y_, TYPE _z_ = 0) 
    : x(_x_), y(_y_), z(_z_) {};
    //要用 G++ 提交 ，可以不用这个 
    POINT operator =(const POINT &A){  
          x = A.x; 
          y = A.y;
          z = A.z;
    }
};

// 多边形 ,逆时针或顺时针给出x,y 
struct POLY {/*验证*/
       //n个点 
       int n; 
       //x,y为点的指针，首尾必须重合
       TYPE * x;   
       TYPE * y; 
       POLY() : n(0), x(NULL), y(NULL) {};
       POLY(int _n_, const TYPE * _x_, const TYPE * _y_) {
                   n = _n_;
                x = new TYPE[n + 1];
                memcpy(x, _x_, n*sizeof(TYPE)); 
                x[n] = _x_[0];
                y = new TYPE[n + 1];
                memcpy(y, _y_, n*sizeof(TYPE)); 
                y[n] = _y_[0];
       }
}; 

//判断 x 是正数还是负数 
inline int Sign(TYPE x){/*验证*/
    return x<-Epsilon?-1:x>Epsilon;
}

void Interect(POINT x,POINT y,TYPE a,TYPE b,TYPE c,int &s,POINT q[]){
    TYPE u=fabs(a*x.x+b*x.y+c);
    TYPE v=fabs(a*y.x+b*y.y+c);
    q[++s].x=(x.x*v+y.x*u)/(u+v);
    q[s].y=(x.y*v+y.y*u)/(u+v);
}
//利用半平面切割
void Cut(TYPE a,TYPE b,TYPE c,int &KarnalPoint,POINT p[])
{
    int s=0;
    int i;
    POINT q[MaxPoint];
    for(i=1; i<=KarnalPoint; i++){//遍历所有顶点是否能观察到该边
    
           if(Sign(a* p[i].x+ b*p[i].y+ c)>= 0){//因为线段是顺时针给出的，如果是逆时针就是<=0
           
                q[++s]= p[i]; //若是则存储
           }
           else{
           
                if(Sign(a* p[i-1].x+ b* p[i-1].y+ c)> 0)//逆时针就是<0
                     Interect(p[i-1], p[i], a, b, c, s, q);
                     
                if(Sign(a* p[i+1].x+ b* p[i+1].y+ c)> 0)//逆时针就是<0
                     Interect(p[i+1], p[i], a, b, c, s, q);
           }
    }
    //最后的p数组存放半平面的点集合
    for(i=1;i<=s;i++) 
           p[i]=q[i];
    p[s+1]=p[1],p[0]=p[s];
    KarnalPoint=s;
}

POLY PolygonKernal(int n, POINT point[])
{    
    int KarnalPoint= n;
    POINT p[MaxPoint];//p 的大小和 tr 的大小一样
       for(int i= 0; i< n; i++){
            p[i+1]= point[i]; //初始化边界
       }
       point[n]= point[0];
       p[n+1]= p[1];
       p[0]  = p[n];
       
       TYPE a,b,c;
       for(int i=0;i<n;i++){
            a=point[i+1].y- point[i].y  ; //计算出相邻两点所在直线ax+by+c=0
            b=point[i].x  - point[i+1].x;
            c=point[i+1].x* point[i].y- point[i].x* point[i+1].y;
            Cut(a, b, c, KarnalPoint, p);
       }
       
       TYPE X[MaxPoint],Y[MaxPoint];
       for(int i= 0; i< KarnalPoint; i++){
           X[i]= p[i].x;
           Y[i]= p[i].y;
       }
       POLY poly(KarnalPoint, X, Y);
 
      return poly;    
}
//求多边形面积  拍好序的点 (返回的有可能是负数，Abs 一下)
TYPE Area(const POLY & poly) { /*验证*/
    if ( poly.n < 3) 
         return TYPE(0); 
    double  s = poly.y[0] * (poly.x[poly.n - 1] - poly.x[1]); 
    for (int i = 1; i < poly.n; i++) { 
        s += poly.y[i] * (poly.x[i - 1] - poly.x[(i + 1) % poly.n]);    
    } 
    return s/2;
}
int main(){
    int n,t;
    POINT point[MaxPoint];
    cin>>t;
    while(t--)
    {
        cin>>n;
           for(int i=0;i<n;i++){
                scanf("%lf%lf",&point[i].x,&point[i].y);
           }
           POLY poly= PolygonKernal(n, point);
           
           TYPE s= Area(poly);
           //之前用自定义的 Abs 宏定义，出现错误，改成下面的就 AC 了
           if(s< 0) printf("%.2f\n",-s);
           else     printf("%.2f\n", s);
                   
           
    }
    return 0;
}