Gym - 101158J -Cover the Polygon with Your Disk（模拟退火+多边形与圆面积的交）

（一）题面：

A convex polygon is drawn on a flat paper sheet. You are trying to place a disk in your hands to cover as large area of the polygon as possible. In other words, the intersection area of the polygon and the disk should be maximized.

Input

The input consists of a single test case, formatted as follows. All input items are integers.
n r
x1 y1
. . .
xn yn
n is the number of vertices of the polygon (3 ≤ n ≤ 10). r is the radius of the disk (1 ≤ r ≤ 100). xi and yi give the coordinate values of the i-th vertex of the polygon (1 ≤ i ≤ n). Coordinate values satisfy 0 ≤ xi ≤ 100 and 0 ≤ yi ≤ 100.

The vertices are given in counterclockwise order. As stated above, the given polygon is convex. In other words, interior angles at all of its vertices are less than 180◦ . Note that the border of a convex polygon never crosses or touches itself.

Output

Output the largest possible intersection area of the polygon and the disk. The answer should not have an error greater than 0.0001 (10−4 ).

Sample Input 1

4 4
0 0
6 0
6 6
0 6

Sample Output 1

35.759506 24 

 

（二）题意：

在二维平面中有一个多边形以及一个圆，现给定多边形各个顶点的坐标以及圆的半径，让你任意平移圆，使得圆与多边形的交的面积最大，输出最大的重叠面积。

 

（三）题解：

对于多边形与圆的面积的交的求解有模板，这里不再做过多赘述。

那么我们对于给定的任意一个圆心，我们都可以求出多边形与圆的面积的交，问题是怎么得到最大的面积的交。

问题可以转化为：对于二维平面中的每一个点，我们都可以得到一个值，并且点与点之间可以通过不断地位移得到，我们的目标是找到一个目标值最优的点。到这里我觉得如果会模拟退火的话就可以想到这么写了，不会的话也没什么办法（实在想不到怎么向模拟退火上靠了（＞人＜；））。

（模拟退火不懂的可以先看看模拟退火简介）

 

（四）代码：

#include<iostream>
#include<cstring>
#include<string>
#include<cstdio>
#include<vector>
#include<cmath>
#include<ctime>
#include<algorithm>
#define eps 1e-14
#define PI acos(-1)
#define LL long long
#define RD T*(rand()*2-RAND_MAX)
#define Drand (long double)rand()/RAND_MAX
using namespace std;
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 d){
    return Vector(A.x*d,A.y*d);
}
Vector operator / (Vector A,double d){
    return Vector(A.x/d,A.y/d);
}
bool operator < (const Point &A,const Point &B){
    return A.x<B.x||(A.x==B.x&&A.y<B.y);
}
int dcmp(double x){
    if(fabs(x)<eps)return 0;
    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;
}
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);
}
struct Circle{
    Point c;
    double r;
    Circle(Point _c,double _r){
        c=_c;r=_r;
    }
    Point point(double rad){
        return Point(c.x+cos(rad)*r,c.y+sin(rad)*r);
    }
};
double TriAngleCircleInsection(Circle C, Point A, Point B){ //三角形的有向面积
    Vector OA = A-C.c, OB = B-C.c;
    Vector BA = A-B, BC = C.c-B;
    Vector AB = B-A, AC = C.c-A;
    double DOA = Length(OA), DOB = Length(OB),DAB = Length(AB), r = C.r;
    if(dcmp(Cross(OA,OB)) == 0) return 0;
    if(dcmp(DOA-C.r) < 0 && dcmp(DOB-C.r) < 0) return Cross(OA,OB)*0.5;
    else if(DOB < r && DOA >= r) {
        double x = (Dot(BA,BC) + sqrt(r*r*DAB*DAB-Cross(BA,BC)*Cross(BA,BC)))/DAB;
        double TS = Cross(OA,OB)*0.5;
        return asin(TS*(1-x/DAB)*2/r/DOA)*r*r*0.5+TS*x/DAB;
    }
    else if(DOB >= r && DOA < r) {
        double y = (Dot(AB,AC)+sqrt(r*r*DAB*DAB-Cross(AB,AC)*Cross(AB,AC)))/DAB;
        double TS = Cross(OA,OB)*0.5;
        return asin(TS*(1-y/DAB)*2/r/DOB)*r*r*0.5+TS*y/DAB;
    }
    else if(fabs(Cross(OA,OB)) >= r*DAB || Dot(AB,AC) <= 0 || Dot(BA,BC) <= 0) {
        if(Dot(OA,OB) < 0){
            if(Cross(OA,OB) < 0) return (-acos(-1.0)-asin(Cross(OA,OB)/DOA/DOB))*r*r*0.5;
            else return ( acos(-1.0)-asin(Cross(OA,OB)/DOA/DOB))*r*r*0.5;
        }
        else return asin(Cross(OA,OB)/DOA/DOB)*r*r*0.5;
    }
    else {
        double x = (Dot(BA,BC)+sqrt(r*r*DAB*DAB-Cross(BA,BC)*Cross(BA,BC)))/DAB;
        double y = (Dot(AB,AC)+sqrt(r*r*DAB*DAB-Cross(AB,AC)*Cross(AB,AC)))/DAB;
        double TS = Cross(OA,OB)*0.5;
        return (asin(TS*(1-x/DAB)*2/r/DOA)+asin(TS*(1-y/DAB)*2/r/DOB))*r*r*0.5 + TS*((x+y)/DAB-1);
    }
}
double AreaOfPolygonAreaToCircle(Circle c,Point *Polygon,int n){   //三角形面积和
    double Area=0;
    for(int i=0;i<n;i++)
        Area+=TriAngleCircleInsection(c,Polygon[i],Polygon[(i+1)%n]);
    return fabs(Area);
}
Point polygon[15];
int dir[4][2]={1,0,-1,0,0,1,0,-1};
int main(){
    #ifdef DanDan
    freopen("in.txt","r",stdin);
    #endif // DanDan

    int n;Point sp,np;
    srand(time(0));
    double r,T=200,dt=0.999;
    scanf("%d%lf",&n,&r);
    for(int i=0;i<n;i++){
        scanf("%lf%lf",&polygon[i].x,&polygon[i].y);
        sp.x+=polygon[i].x;sp.y+=polygon[i].y;
    }sp.x/=n;sp.y/=n;
    double ans=AreaOfPolygonAreaToCircle(Circle(sp,r),polygon,n),best=ans;
    while(T>eps){
        np.x=sp.x+RD;
        np.y=sp.y+RD;
        double res=AreaOfPolygonAreaToCircle(Circle(np,r),polygon,n);
        if(best<ans)best=ans;
        if(ans<res)ans=res,sp=np;
        T*=dt;
    }
    printf("%.12f\n",best);
    return 0;
}

 

（五）总结：

模拟退火练手题

对于可以用模拟退火做的题目，特征还是比较明显的

这个题的退火速率不能太快（貌似dt=0.998都不行？），可能是因为坐标值比较小

但实际上最开始我的dt=0.98，但是退火时是这么写的：

    while(T>eps){
        for(int i=0;i<15;i++){
            np.x=sp.x+RD;
            np.y=sp.y+RD;
            double res=AreaOfPolygonAreaToCircle(Circle(np,r),polygon,n);
            if(res>ans)ans=res,sp=np;
        }
        T*=dt;
    }

也就是对于每一个温度，都使其“稳定一段时间”，效果也比较好，但是为了统一写法，最后还是改了退火速率。