（intermediate) 圆与凸多边形相交 UVA 11177 - Fighting Against a Polygonal Monster

Problem C
Fighting against a Polygonal Monster
Input: Standard Input

Output: Standard Output

 

You're fighting against a monster. You have a special weapon that can shoot a laser beam. The laser beam can be seen as a cylinder along the z-axis, So when look from above the XY plane, the laser beam is actually a circle centered at (0,0) (yes, you cannot change the position of center). The monster is a convex polygon (well, you may think of it as a very thin prism) with n vertices surrounding the origin (i.e. (0,0) is strictly inside the monster, not on his boundary).

 

You know that, when the common area of the laser beam and the monster is at least R, the monster dies. Since larger laser beam consumes more power, you're interested in the minimum radius of the laser beam.

 

Write a program to find the minimum radius. It is guaranteed that (0 £ R £A), whereA is the area of the monster.

 

Input

The input consists at most 10 cases. Each case starts with a single integer n (3 £ n £ 50) and a floating-point number R followed by n lines of two real numbers: the coordinates of the monster. The points are arranged in counter-clockwise order or clockwise order. The last test case is followed by a single zero, which should not be processed. The meaning of N and R are given in the problem statement.

 

Output

For each test case, print the case number and the minimum radius, to 2 decimal places. Inputs will be such that small precision errors will not change the visible output if you use double-precision floating-point numbers.

 

Sample Input                             Output for Sample Input

3 1.60

-1 -1

1 -1

0 1

0

 

Case 1: 0.93

---

Problem setter: Rujia Liu

Alternative Solutions: Wenbin Tang, Yiming Li

Originally used in the Sichuan Team Selection Contest for NOI, 2006.

 

 

题意：有一个凸多边形，原点被严格包围在里面，你要以原点为圆心，画一个半径尽可能小的圆使得圆和凸多边形的公共部分的面积不少于要求的面积。

思路：二分半径我想大部分人都能想到。这个题目面积讨论起来其实很麻烦的。。。我的讨论方法是拿出凸多边形的一条边，看看这条边的两个端点在圆的什么地方，如果都在圆外，显然面积是扇形，如果两个端点都在圆内，显然是三角形，如果一个在外面，一个在里面，显然是一个扇形+三角形，如果一个在圆上，一个在圆外，是扇形，如果一个在圆上，一个在圆内，是三角形，如果两个都在外面，有两种情况，交点不在线段上面的话是扇形，如果都在线段上面是两个扇形+一个三角形。不可能一个点在线段上，一个不在的。。。但是不知道为什么我的代码会出现这种情况。。。。不过我改了个写法就过了，改成用OnSegment不会有问题了。慢慢写，把情况分清楚，你能做出来的。

代码：

#include<iostream>#include<cstdio>#include<string.h>#include<math.h>#include<string>#include<cstring>#include<map>#include<algorithm>#include<vector>using namespace std;#define eps 1e-10const double inf = 1e16;const double PI = acos(-1.0);struct Point{Point(const Point&p) { x = p.x , y = p.y; }Point (double xx=0,double yy=0) : x(xx) , y(yy) { }double x;double y;};typedef Point Vector;Vector operator+(Vector  v1,Vector  v2) { return Vector(v1.x+v2.x,v1.y+v2.y); }Vector operator-(Vector  v1,Vector  v2) { return Vector(v1.x-v2.x,v1.y-v2.y); }Vector operator*(Vector  v, double p) { return Vector(v.x*p,v.y*p); }Vector operator/(Vector  v,double p) { return Vector(v.x/p,v.y/p); }bool operator < (Point  a,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;}inline double toRad(double x) { return x * PI/180; }inline double toDegreed(double rad) { return rad*180/PI; }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 Angle(Vector A,Vector B) { return acos(Dot(A,B)/Length(A)/Length(B)); }double Cross(Vector A,Vector B) { return A.x*B.y-A.y*B.x; }double Area2(Point a,Point b,Point c) {  return fabs(Cross(b-a,c-a)); }//旋转Vector Rotate(Vector A,double rad) {return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));}//单位法线Vector Normal(Vector A) { double L = Length(A); return Vector(-A.y/L,A.x/L); }//点和直线Point GetLineIntersection(Point P,Vector v,Point Q,Vector w){Vector u = P-Q;double t = Cross(w,u) / Cross(v,w);return P+v*t;}double DistanceToLine(Point P,Point A,Point B) {Vector v1 = B-A , v2 = P-A;return fabs(Cross(v1,v2))/Length(v1);}double DistanceToSegment(Point P,Point A,Point B){if (A==B) return Length(P-A);Vector v1 = B-A , v2 = P-A , v3 = P-B;if (dcmp(Dot(v1,v2)) < 0) return Length(v2);else if (dcmp(Dot(v1,v3)) > 0) return Length(v3);else return fabs(Cross(v1,v2))/Length(v1);}//点在直线上的投影Point GetLineProjection(Point P,Point A,Point B) {Vector v = B-A;return A+v*(Dot(v,P-A)/Dot(v,v));}//线段相交(不包括端点）bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2) {double c1 = Cross(a2-a1,b1-a1) , c2 = Cross(a2-a1,b2-a1) ,     c3 = Cross(b2-b1,a1-b1) , c4 = Cross(b2-b1,a2-b1);return dcmp(c1)*dcmp(c2) < 0 && dcmp(c3)*dcmp(c4)<0;}//点在线段上（不包括端点）bool OnSegment(Point p,Point a,Point b){return dcmp(Cross(a-p,b-p))==0 && dcmp(Dot(a-p,b-p)) < 0;}//--------------------------------------------------------------------------------------------//直线和直线struct Line{Point P;//直线上任意一点Vector v;// 方向向量。它的左边就是对应的半平面double ang;//极角Line() { }Line(Point P,Vector v) { this->P = P ; this->v = v; ang = atan2(v.y,v.x); }bool operator < (const Line& L) const { return ang < L.ang; } //排序用的比较运算符Point point(double t) { return v*t+P; }};//点p在有向直线L的左边（线上不算）bool OnLeft(Line L , Point p) { return Cross(L.v,p-L.P) > 0; }//二直线交点。假定交点唯一存在Point GetIntersection(Line a,Line 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(Line* L,int n,Point* poly){sort(L,L+n);     //按极角排序int first , last;//双端队列的第一个元素和最后一个元素的下表Point *p = new Point[n];//p[i]为q[i]和q[i+1]的交点Line *q = new Line[n];//双端队列q[first=last=0] = L[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]);//计算首尾两个半平面的交点//从deque复制到输出中int m = 0;for (int i = first ; i <= last ; ++i) poly[m++] = p[i];return m ;}//--------------------------------------------//与圆相关struct Circle {Circle() { }Point c;double r;Circle(Point c, double r) : c(c) , r(r) { }Point point (double a) { return Point(c.x+cos(a)*r,c.y+sin(a)*r); }};int getLineCircleIntersection(Line L,Circle C,double &t1,double &t2,vector<Point>& sol){double a = L.v.x , b = L.P.x-C.c.x , c= L.v.y, d = L.P.y-C.c.y;double e = a*a+c*c , f = 2*(a*b+c*d) , g = b*b+d*d-C.r*C.r;double delta = f*f-4*e*g;//判别式if (dcmp(delta) < 0) return 0;//相离if (dcmp(delta)==0) {                   //相切t1 = t2 = -f/(2*e); sol.push_back(L.point(t1));return 1;}//相交t1 = (-f-sqrt(delta)) / (2*e); sol.push_back(L.point(t1));t2 = (-f+sqrt(delta)) / (2*e); sol.push_back(L.point(t2));return 2;}double angle(Vector v) { return atan2(v.y,v.x); }int getCircleCircleIntersection(Circle C1,Circle C2,vector<Point>& sol){double d = Length(C1.c-C2.c);if (dcmp(d)==0) {if (dcmp(C1.r-C2.r)==0) return -1;//两圆重合return 0;}if (dcmp(C1.r+C2.r-d) < 0) return 0;if (dcmp(fabs(C1.r-C2.r)-d) > 0) return 0;double a = angle(C2.c-C1.c);double da = acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));//向量C1C2的极角//C1C2到C1P1的角Point p1 = C1.point(a-da) , p2 = C1.point(a+da);sol.push_back(p1);if (p1==p2) return 1;sol.push_back(p2);return 2;}//国电p到圆C的切线。v[i]是第i条切线的向量。返回切线条数int getTangents(Point p,Circle C,Vector* v){Vector u= C.c-p;double dist = Length(u);if (dist < C.r) return 0;else if (dcmp(dist-C.r)==0) {v[0] = Rotate(u,PI/2);return 1;} else {double ang = asin(C.r/dist);v[0] = Rotate(u,-ang);v[1] = Rotate(u,+ang);return 2;}}int getTangents(Circle A,Circle B,Point* a, Point* b){int cnt = 0;if (A.r < B.r) { swap(A,B); swap(a,b); }double d2 = Dot(A.c-B.c,A.c-B.c);double rdiff = A.r-B.r;double rsum = A.r+B.r;if (dcmp(d2-rdiff*rdiff) < 0) return 0;//内含double base = atan2(B.c.y-A.c.y,B.c.x-A.c.x);if (d2==0 && dcmp(A.r-B.r)==0) return -1;//无限多条切线if (dcmp(d2-rdiff*rdiff)==0) {                       //内切，1条切线a[cnt] = A.point(base);b[cnt] = B.point(base);++cnt;return 1;}//有外切共线double ang = acos((A.r-B.r)/sqrt(d2));a[cnt] = A.point(base+ang); b[cnt] = B.point(base+ang); ++cnt;a[cnt] = A.point(base-ang); b[cnt] = B.point(base-ang); ++cnt;if (dcmp(d2-rsum*rsum)==0) {                  //一条内公切线a[cnt] = A.point(base);b[cnt] = B.point(PI+base);++cnt;} else if (dcmp(d2-rsum*rsum)>0) {           //两条公切线double ang = acos((A.r+B.r)/sqrt(d2));a[cnt] = A.point(base+ang); b[cnt] = B.point(PI+base+ang); ++cnt;a[cnt] = A.point(base-ang); b[cnt] = B.point(PI+base-ang); ++cnt;}return cnt;}//点p和圆的关系: 0:在圆上 1：在圆外 -1：在圆内int PointCircleRelation(Point p,Circle c) {return dcmp(Dot(p-c.c,p-c.c)-c.r*c.r);}//A在B内bool InCircle(Circle A,Circle B){if (dcmp(A.r-B.r)>0) return false;double d2 = Dot(A.c-B.c,A.c-B.c);double rdiff = A.r-B.r;double rsum = A.r+B.r;if (dcmp(d2-rdiff*rdiff) <= 0) return true;//内含或内切或重合return false;}//----------------------------------------------------------------------------//与球相关的转化//角度转换为弧度double torad(double deg) { return deg/180*PI; }//经纬度（角度）转化为空间坐标void get_coord(double R,double lat,double lng,double& x,double& y,double& z){lat = torad(lat);lng = torad(lng);x = R*cos(lat)*cos(lng);y = R*cos(lat)*sin(lng);z = R*sin(lat);}//-----------------------------------------------------------------------const int maxn = 50+5;int n;double S;Point p[maxn];const Point center = Point(0,0);bool Ok(double R){Circle c(center,R);double T = 0 , A = 0;for (int i = 0 ; i < n ; ++i) {Vector v1 = p[i] , v2 = p[(i+1)%n];int k1 = PointCircleRelation(v1,c) , k2 = PointCircleRelation(v2,c);if (k1==0 && k2==0) { T += Area2(center,v1,v2); continue; }if (k1+k2==1) { A += Angle(v1,v2); continue; }if (k1+k2<=-1) { T += Area2(center,v1,v2); continue; }vector<Point> inter;double t1 , t2;getLineCircleIntersection(Line(v1,v2-v1),c,t1,t2,inter);if (k1==-1 && k2==1) {if (OnSegment(inter[0],v1,v2)) {A += Angle(v2,inter[0]);T += Area2(center,v1,inter[0]);} else if (OnSegment(inter[1],v1,v2)) {A += Angle(v2,inter[1]);T += Area2(center,v1,inter[1]);};continue;}if (k1==1 && k2==-1) {if (OnSegment(inter[0],v1,v2)) {A += Angle(v1,inter[0]);T += Area2(center,v2,inter[0]);} else if (OnSegment(inter[1],v1,v2)) {A += Angle(v1,inter[1]);T += Area2(center,v2,inter[1]);}continue;}if (inter.size() < 2) A += Angle(v1,v2);else if (inter.size()==2) {if (dcmp(t1)>0 && dcmp(t1-1)<0) {if (dcmp(t2)>0 && dcmp(t2-1)<0) {A += Angle(v1,v2)-Angle(inter[0],inter[1]);T += Area2(center,inter[0],inter[1]);} else if (k1==1 && k2==1) while (true);}else A += Angle(v1,v2);}}return dcmp((T+A*R*R)/2-S) >= 0;}int main(){int k = 0;while (scanf("%d",&n)==1,n) {++k;scanf("%lf",&S);double left = 0 , right = 0;for (int i = 0 ; i < n ; ++i) {scanf("%lf%lf",&p[i].x,&p[i].y);right = max(right,Length(p[i]));}double mid = left+(right-left)*0.5;double ans = inf;while (dcmp(left-right)<0) {if (Ok(mid)) {ans = min(ans,mid);right = mid;} else left = mid;mid = left+(right-left)*0.5;}printf("Case %d: %.2lf\n",k,ans+eps);}}