poj3384半平面交

链接：http://poj.org/problem?id=3384

题意：一个凸多边形，用两个半径相等的圆来覆盖，要求覆盖面积最大时的圆心坐标。

思路：把凸多边形每条边向内移r的距离，用半平面交求出新的多边形，由于要使覆盖面积最大，所以两个圆心应该是多边形的顶点中距离最远的两个。要注意当这个多边形退化为一个点时，此时两个圆心是重合的。一直wa，内伤，最后发现是模板,的问题。。。《训练指南》上的半平面交模板，当半平面交退化为点和线段时返回值都是0，不适合此处，然后改了一下onleft()函数，如果最后的结果其实只有一个点，返回的m不为1，而是重复的个数，虽然这道题过了，但是算法有问题，或可用增量法求半平面交。

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
const int maxn=210;
const double eps=1e-10;
struct Point
{
    double x,y;
    Point(double x=0,double y=0):x(x),y(y) {}
};
typedef Point Vector;
Vector operator + (Vector A,Vector B)
{
    return Vector(A.x+B.x,A.y+B.y);
}
Vector operator - (Vector A,Vector B)
{
    return Vector(A.x-B.x,A.y-B.y);
}
Vector operator * (Vector A,double p)
{
    return Vector(A.x*p,A.y*p);
}
int dcmp(double x)
{
    if(fabs(x)<eps) return 0;
    else return x<0?-1:1;
}

bool operator == (const Point& a,const Point& b)
{
    return dcmp(a.x-b.x)==0 && dcmp(a.y-b.y)==0;
}
bool operator < (const Point& a,const Point& b)
{
    return a.x<b.x || (a.x==b.x && a.y<b.y);
}
double Dot(Vector A,Vector B)
{
    return A.x*B.x+A.y*B.y;
}
double Length(Vector A)
{
    return sqrt(Dot(A,A));
}
double Cross(Vector A,Vector B)
{
    return A.x*B.y-A.y*B.x;
}
Vector Normal(Vector A)
{
    double L=Length(A);
    return Vector(-A.y/L,A.x/L);
}
double dist(Point a,Point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double dist2(Point a,Point b)
{
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
struct DLine
{
    Point P;
    Vector v;
    double ang;
    DLine() {}
    DLine(Point P,Vector v):P(P),v(v)
    {
        ang=atan2(v.y,v.x);
    }
    bool operator < (const DLine& L) const
    {
        return ang<L.ang;
    }
};
bool OnLeft(DLine L,Point p)
{
    return Cross(L.v,p-L.P)>=0;
}
Point GetIntersection(DLine a,DLine b)
{
    Vector u=a.P-b.P;
    double t=Cross(b.v,u)/Cross(a.v,b.v);
    return a.P+a.v*t;
}
int HalfPlaneIntersection(DLine* L,int n,Point* poly)
{
    sort(L,L+n);
    int first,last;
    Point *p=new Point[n];
    DLine *q=new DLine[n];
    q[first=last=0]=L[0];
    for(int i=1;i<n;i++)
    {
        while(first<last && !OnLeft(L[i],p[last-1])) last--;
        while(first<last && !OnLeft(L[i],p[first])) first++;
        q[++last]=L[i];
        if(fabs(Cross(q[last].v,q[last-1].v))<eps)
        {
            last--;
            if(OnLeft(q[last],L[i].P)) q[last]=L[i];
        }
        if(first<last) p[last-1]=GetIntersection(q[last-1],q[last]);
    }
    while(first<last && !OnLeft(q[first],p[last-1])) last--;
    if(last-first <= 1) return 0;
    p[last]=GetIntersection(q[last],q[first]);
    int m=0;
    for(int i=first;i<=last;i++) poly[m++]=p[i];
    return m;
}
Point p[maxn],poly[maxn];
DLine L[maxn];
Vector v[maxn],v2[maxn];
void solve(int m)
{
    Point p1=poly[0],p2=poly[0];
    if(m==1)
    {printf("%.4f %.4f %.4f %.4f
",p1.x,p1.y,p2.x,p2.y);return;}
    double ans=0.0,d;
    for(int i=0;i<m;i++)
        for(int j=0;j<m;j++)
        {
           if(i==j) continue;
           d=dist2(poly[i],poly[j]);
           if(dcmp(d-ans)>=0)
           {
               ans=d;
               p1=poly[i],p2=poly[j];
           }
        }
    printf("%.4f %.4f %.4f %.4f
",p1.x,p1.y,p2.x,p2.y);
}

int main()
{
//    freopen("ine.cpp","r",stdin);
    int n,r;
    while(~scanf("%d%d",&n,&r))
    {
        for(int i=n-1;i>=0;i--)
            scanf("%lf%lf",&p[i].x,&p[i].y);
        for(int i=0;i<n;i++)
        {
            v[i]=p[(i+1)%n]-p[i];
            v2[i]=Normal(v[i]);
            L[i]=DLine(p[i]+v2[i]*r,v[i]);
        }
        int m=HalfPlaneIntersection(L,n,poly);
        solve(m);
    }
    return 0;
}