「Gym101158J」Cover the Polygon with Your Disk【模拟退火】

wzw1105

于 2019-09-27 13:11:03 发布

阅读量230

点赞数

分类专栏：计算几何文章标签：计算几何

本文链接：https://blog.csdn.net/qq_40791842/article/details/101528093

版权

计算几何专栏收录该内容

12 篇文章 0 订阅

订阅专栏

Cover the Polygon with Your Disk

Input: Standard Input Time Limit: 5 seconds

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$
$x_1\ y_1$
.
.
.
$x_n\ y_n$
$n$ is the number of vertices of the polygon $\leq n \leq 10).$ r is the radius of the disk $\leq r \leq 100).$ $x_i$ and $y_i$ give the coordinate values of the $i - t h$ vertex of the polygon $\leq i \leq n).$ Coordinate values satisfy $\leq x_i \leq 100$ and $\leq y_i \leq 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

Sample Output 1

35.759506

Sample Input 2

Sample Output 2

2.113100

Sample Input 3

Sample Output 3

0.019798

Sample Input 4

Sample Output 4

3.137569

Sample Input 5

Sample Output 5

177.728187

Sample Input 6

Sample Output 6

7181.603297

题意

就是给你一个凸多边形，以及某一个圆的半径 $r$ ，然后让你去找一个平面上的圆心，使得这个半径为 $r$ 的圆与凸多边形相交的面积最大

题解

因为是凸多边形，不难想到三分套三分去搞
然而可以模拟退火去做，别问我为啥，因为模拟退火就是找一个最优的实数点的
但是这个不用模拟退火中的那个“不是比局部最优解更优的答案对应的随机出来的点以一定概率去接受这个点”，用了答案不对？希望路过的大佬赐教，去掉就对了，感觉这不是模拟退火而变成了爬山算法？

代码

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-8;
const int maxn=20005;
#define inf 0x3f3f3f3f3f3f3f3f
#define pi acos(-1.0)
 
int sgn(double k) {return k<-eps?-1:(k<eps?0:1);}
double Acos(double a) {return a<=-1?pi:(a>=1?0:acos(a));}
double Sqrt(double a) {return sgn(a)==0?0:sqrt(a);}
 
typedef struct point{  //点结构体
    double x,y;
    point(double a=0,double b=0) {x=a;y=b;} 
    point operator+(point other) {return point(x+other.x,y+other.y);}
    point operator-(point other) {return point(x-other.x,y-other.y);}
    point operator*(double k) {return point(x*k,y*k);}
    point operator/(double k) {return point(x/k,y/k);} 
    //点乘，即数量积，内积(ABcos<A,B>)
    double operator*(point other) {return x*other.x+y*other.y;}
    //叉乘，即向量积，外积(ABsin<A,B>)
    double operator^(point other) {return x*other.y-y*other.x;}
    friend bool operator<(const point &p1,const point &p2) {
        if(sgn(p1.x-p2.x)==0) return p1.y<p2.y;
        return p1.x<p2.x;
    }
    static bool compare_x(const point &p1,const point &p2) {return p1.x<p2.x;}
    static bool compare_y(const point &p1,const point &p2) {return p1.y<p2.y;}
    friend bool operator==(const point &p1,const point &p2) {
        return sgn(p1.x-p2.x)==0&&sgn(p1.y-p2.y)==0;
    }
    friend bool operator!=(const point &p1,const point &p2) {
        return sgn(p1.x-p2.x)!=0||sgn(p1.y-p2.y)!=0;
    }
    //o-p长度
    friend double len(point p) {return Sqrt(p.x*p.x+p.y*p.y);}
    friend double len2(point p) {return p.x*p.x+p.y*p.y;}
    //op与x轴夹角，范围[-pi,pi]
    friend double angle(point p) {return atan2(p.y,p.x);}    
    //判断点和有向直线的关系: -1:直线s-e方向左侧 0:直线上 1:直线右侧
    friend int point_to_line(point p,point s,point e) {return -sgn((e-s)^(p-s));}
    //op1-op2无方向角度，即与原点连线夹脚，返回值的范围[0,pi]
    friend double angle(point p1,point p2) {return Acos((p1*p2)/len(p1)/len(p2));}
    //将向量p长度变为len返回一个新的向量，方向不变
    friend point extend(point p,double l) { 
        if(sgn(len(p))==0) return point(0,l);  //0向量
        return point(p*(l/len(p)));
    }
    //点p2绕着点p1顺时针方向旋转a角度，a为弧度,返回新的point
    friend point rotate(point p1,point p2,double a) {
        point vec=p2-p1;
        double xx=vec.x*cos(a)+vec.y*sin(a);
        double yy=vec.y*cos(a)-vec.x*sin(a);
        return point(p1.x+xx,p1.y+yy);
    }
    //三角形面积
    friend double area(point p1,point p2,point p3) {
        return fabs(((p2-p1)^(p3-p1))/2);
    }
}Vector;
 
 
struct line{
    point s,e;
    line(point a,point b) {s=a;e=b;}
    line(double a=0,double b=0,double c=0,double d=0) {s=point(a,b);e=point(c,d);}
    //半平面交所用极角排序函数
    static bool compare_half_plane(line a,line b){
        double A=angle(a.e-a.s),B=angle(b.e-b.s);
        if(sgn(A-B)!=0) return A<B;
        return sgn((a.e-a.s)^(b.e-a.s))>0;
    }
    //计算直线p1-p2是否与直线q1-q2的交点，注意不能平行
    /*需要重载point*double  */
    friend point intersect(line l1,line l2) {
        double k=((l2.e-l2.s)^(l2.s-l1.s))/((l2.e-l2.s)^(l1.e-l1.s));
        return l1.s+(l1.e-l1.s)*k;
    }
    //返回直线倾斜角  0<=angle<=pi
    friend double angle(line p) {
        double k=atan2(p.e.y-p.s.y,p.e.x-p.s.x);
        if(sgn(k)==0) k+=pi;
        if(sgn(k-pi)==0) k-=pi;
        return k;
    }
 
    //两条直线平行或重合
    friend  bool parallel(line l1,line l2) {
        return sgn((l1.e-l1.s)^(l2.e-l2.s))==0;
    }
 
    //判断点是否在线段上
    friend bool onseg(point p,line l) {
        return sgn((l.s-p)^(l.e-p))==0&&sgn((l.s-p)*(l.e-p))<=0;
    } 
    //判断点和有向直线的关系: -1:直线s-e方向左侧 0:直线上 1:直线右侧
    friend int point_to_line(point p,line l) {
        return -sgn((l.e-l.s)^(p-l.s));
    }
    friend bool onright(point p,line l) {
        return ((l.e-l.s)^(p-l.s))<=0;
    }
    //线段与线段关系: 0:平行没有交点 1:平行有交点 2:不平行没有交点 3:不平行有交点
    friend int seg_to_seg(line l1,line l2) {
        if(sgn((l1.e-l1.s)^(l2.e-l2.s))==0) {
            if(onseg(l1.s,l2)||onseg(l1.e,l2)||onseg(l2.s,l1)||onseg(l2.e,l1)) return 1;
            else return 0;
        }
        point inter=intersect(l1,l2);
        if(onseg(inter,l1)&&onseg(inter,l2)) return 3;
        return 2;
    }
 
    //两直线关系: 0:平行 1:重合 2：相交
    friend int line_to_line(line l1,line l2) {
        if(parallel(l1,l2)) return point_to_line(l1.s,l2)==0;
        return 2;
    }
    //点a到直线l的距离
    friend double dis_point_line(point p,line l) { 
        if(onseg(p,l)) return 0;
        double cos0=(p-l.s)*(l.e-l.s)/(len(p-l.s)*len(l.e-l.s));
        return len(p-l.s)*sin(Acos(cos0));
    }
    //点到线段l的距离
    friend double dis_point_seg(point p,line l) {
        if(sgn((p-l.s)*(l.e-l.s))<0||sgn((p-l.e)*(l.s-l.e))<0) return min(len(p-l.s),len(p-l.e));
        return dis_point_line(p,l);
    }
    //线段到线段距离
    friend double dis_seg_seg(line l1,line l2) {
        if(seg_to_seg(l1,l2)==1||seg_to_seg(l1,l2)==3) return 0;
        return min(dis_point_seg(l1.s,l2),min(dis_point_seg(l1.e,l2),min(dis_point_seg(l2.s,l1),dis_point_seg(l2.e,l1))));
    }
 
    //返回点a在直线l上的垂足
    friend point chuizu(point p,line l) {
        return l.s+((l.e-l.s)*((l.e-l.s)*(p-l.s)))/len2(l.e-l.s);
    } 
    //点a关于l的对称点
    friend point duichen(point p,line l) {
        point q=chuizu(p,l);
        return point(2*q.x-p.x,2*q.y-p.y);
    }
 
};
 
struct circle{
    point o;double r;
    circle(point a,double b) {o=a;r=b;}
    circle(double a=0,double b=0,double c=0) {o.x=a,o.y=b,r=c;}
    bool operator==(circle c) {return o==c.o&&sgn(r-c.r)==0;}
    double area() {return pi*r*r;}
 
    //求三角形外接圆
    circle(point a,point b,point c) {
        point p1=rotate((a+b)/2,b,3*pi/2),p2=rotate((b+c)/2,c,3*pi/2);
        o=intersect(line((a+b)/2,p1),line((b+c)/2,p2));
        r=len(o-a);
    }
 
    //求三角形内切圆，nouse无实际作用，只是用来区别上面外接圆
    circle(point a,point b,point c,int nouse){
        if(sgn((b-a)^(c-b))<0) swap(b,c);
        point p1=rotate(a,b,2*pi-(angle(c-a,b-a)/2)),p2=rotate(b,c,2*pi-(angle(a-b,c-b)/2));
        o=intersect(line(a,p1),line(b,p2));
        r=dis_point_line(o,line(b,c));
    }
 
    //点和圆的关系:0:圆内 1:圆上 2:圆外
    friend int point_to_circle(point p,circle c) {
        return 1+sgn(len(p-c.o)-c.r);
    }
    //直线和圆的关系:0:与圆相交 1:圆上 2:圆外
    friend int line_to_circle(line l,circle c) {
        return 1+sgn(dis_point_line(c.o,l)-c.r);
    }
 
    //返回直线与圆交点个数，交点坐标存在引用p1,p2中
    friend int line_circle_cross(line l,circle c,point &p1,point &p2) {
        if(line_to_circle(l,c)==2) return 0;
        double d=dis_point_line(c.o,l);
        point p=chuizu(c.o,l);
        d=Sqrt(c.r*c.r-d*d);
        if(sgn(d)==0) {
            p1=p2=p;
            return 1;
        }
        p1=p+extend(l.e-l.s,d);p2=p+extend(l.e-l.s,-d);
        return 2;
    }
 
    //返回三角形oab与圆相交的面积
    friend double circle_triangle_area(circle c,point a,point b) {
        if(sgn((a-c.o)^(b-c.o))==0) return 0.0;
        point q[5],p1,p2;int len=0;q[len++]=a;line l=line(a,b);
        if(line_circle_cross(l,c,q[1],q[2])==2) {
            if(sgn((a-q[1])*(b-q[1]))<0) q[len++]=q[1];
            if(sgn((a-q[2])*(b-q[2]))<0) q[len++]=q[2];
        }
        q[len++]=b;double res=0;
        if(len==4&&sgn((q[0]-q[1])*(q[2]-q[1]))>0) swap(q[1],q[2]);
        for(int i=0;i<len-1;i++) {
            if(point_to_circle(q[i],c)==2||point_to_circle(q[i+1],c)==2) {
                double ang=angle(q[i]-c.o,q[i+1]-c.o);
                res+=c.r*c.r*ang/2;
            }else res+=fabs((q[i]-c.o)^(q[i+1]-c.o))/2.0;
        }
        return res;
    }
 
};
 
struct polygon{
    int n;
    point p[maxn];
    line l[maxn];
 
    //求简单多边形与圆相交的面积
    friend double polygon_circle_area(polygon &pol,circle c) {
        double ans=0;
        for(int i=1;i<=pol.n;i++) {
            int j=i%pol.n+1;
            if(sgn((pol.p[j]-c.o)^(pol.p[i]-c.o))>=0) ans+=circle_triangle_area(c,pol.p[i],pol.p[j]);
            else ans-=circle_triangle_area(c,pol.p[i],pol.p[j]);
        }
        return fabs(ans);
    }
 
    void input() {
        scanf("%d",&n);
        for(int i=1;i<=n;i++) scanf("%lf %lf",&p[i].x,&p[i].y);
    }
    void output() {
        printf("n=%d\n",n);
        for(int i=1;i<=n;i++) printf("%lf %lf\n",p[i].x,p[i].y);
    }
};
 
polygon pol;
double r;
int n;
namespace simulated_annealing {
    double ansx,ansy; //全局最优解的坐标
    double ans;       // 全局最优解，温度
    const double delta=0.996; //降温系数
    inline double dis(double a,double b,double c,double d) {
        return sqrt((a-c)*(a-c)+(b-d)*(b-d));
    }
    inline double calc_energy(double a,double b) {  //计算决策点在(a,b)处时的能量
        return polygon_circle_area(pol,circle(a,b,r));
    }   
    inline void init() {     //初始化
        ansx=ansy=0;
        for(int i=1;i<=n;i++) ansx+=pol.p[i].x,ansy+=pol.p[i].y;
        ansx/=n,ansy/=n;
        ans=calc_energy(ansx,ansy);
    }
    inline void simulate() {
        double nowx=ansx,nowy=ansy; //当前搜索到的点
        double t=100000;
        while(t>1e-14) {
            double nxtx=nowx+(rand()*2-RAND_MAX)*t;
            double nxty=nowy+(rand()*2-RAND_MAX)*t;
            double nxt_energy=calc_energy(nxtx,nxty);
            double delta_energy=nxt_energy-ans;        //能量差
            if(delta_energy>0) nowx=ansx=nxtx,nowy=ansy=nxty,ans=nxt_energy;      //能量降低一定接受新点
            else if(exp(-delta_energy/t)*RAND_MAX>rand()) nowx=nxtx,nowy=nxty;   //以一定的概率接受这个点
            t*=delta;       //降温
        }
    }
    inline void solve() {
        init();
        simulate();
        simulate();
        simulate();
        simulate();
        simulate();
        simulate();
        simulate();
    }
}
using namespace simulated_annealing;
 
int main()
{
    scanf("%d %lf",&pol.n,&r);n=pol.n;
    for(int i=1;i<=pol.n;i++) scanf("%lf %lf",&pol.p[i].x,&pol.p[i].y);
    solve();
    printf("%.10lf\n",ans);
}

wzw1105

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
「Gym101158J」Cover the Polygon with Your Disk【模拟退火】

Problem J Cover the Polygon with Your Disk Input: Standard Input Time Limit: 5 secondsA convex polygon is drawn on a flat paper sheet. You are trying to place a disk in your handsto cover as larg...
复制链接

扫一扫

专栏目录