二分查询答案,判断每一个新形成的向量合在一块能否形成半平面交
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <cmath> 5 #include <algorithm> 6 using namespace std; 7 #define N 110 8 #define eps 1e-7 9 10 int dcmp(double x) 11 { 12 if(fabs(x)<eps) return 0; 13 else return x<0?-1:1; 14 } 15 16 struct Point{ 17 double x , y; 18 Point(double x=0 , double y=0):x(x),y(y){} 19 }po[N] , poly[N]; 20 21 typedef Point Vector; 22 Vector vec[N]; //记录每条边上对应的法向量 23 24 Vector operator+(Vector a , Vector b) { return Vector(a.x+b.x , a.y+b.y); } 25 Vector operator-(Point a , Point b) { return Vector(a.x-b.x , a.y-b.y); } 26 Vector operator*(Vector a , double b) { return Vector(a.x*b , a.y*b); } 27 Vector operator/(Vector a , double b) { return Vector(a.x/b , a.y/b); } 28 bool operator==(const Point &a , const Point &b) { return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0; } 29 30 double Dot(Vector a , Vector b) { return a.x*b.x+a.y*b.y; } 31 double Cross(Vector a , Vector b) { return a.x*b.y - b.x*a.y; } 32 double Length(Vector a) { return sqrt(Dot(a,a)); } 33 double Angle(Vector A , Vector B) { return acos(Dot(A,B)) / Length(A) / Length(B); } 34 double Area2(Point A , Point B , Point C) { return Cross(B-A , C-A); } 35 36 Vector Rotate(Vector A , double rad) { return Vector(A.x*cos(rad)-A.y*sin(rad) , A.x*sin(rad)+A.y*cos(rad)); } 37 38 Vector Normal(Vector a) 39 { 40 double l = Length(a); 41 return Vector(-a.y/l , a.x/l); 42 } 43 44 struct Line{ 45 Point P; 46 Vector v; 47 double ang; 48 Line(){} 49 Line(Point P , Vector v):P(P),v(v){ang = atan2(v.y,v.x);} 50 bool operator<(const Line &m)const { 51 return ang<m.ang; 52 } 53 }line[N] , L[N]; 54 55 bool OnLeft(Line L , Point P) 56 { 57 return Cross(L.v , P-L.P) > 0; 58 } 59 60 Point GetIntersection(Line a , Line b) 61 { 62 Vector u = a.P-b.P; 63 double t = Cross(b.v , u) / Cross(a.v , b.v); 64 return a.P+a.v*t; 65 } 66 67 /***半平面交的主过程,返回形成半平面交点的个数,无法形成就返回0***/ 68 int HalfplaneIntersection(Line *L , int n , Point *poly) 69 { 70 sort(L , L+n); 71 int first , last; 72 Point *p = new Point[n]; 73 Line *q = new Line[n]; 74 q[first=last=0] = L[0]; 75 for(int i=1 ; i<n ; i++) 76 { 77 while(first<last && !OnLeft(L[i] , p[last-1])) last--; 78 while(first<last && !OnLeft(L[i] , p[first])) first++; 79 q[++last] = L[i]; 80 if(fabs(Cross(q[last].v , q[last-1].v)) < eps) 81 { 82 last--; 83 if(OnLeft(q[last] , L[i].P)) q[last]=L[i]; 84 } 85 if(first<last) p[last-1] = GetIntersection(q[last-1] , q[last]); 86 } 87 while(first < last && !OnLeft(q[first] , p[last-1])) last--; 88 //删除无用平面 89 if(last-first<=1) return 0; 90 p[last] = GetIntersection(q[last] , q[first]); 91 92 //从deque复制到输出中 93 int m=0; 94 for(int i=first ; i<=last ; i++) poly[m++] = p[i]; 95 return m; 96 } 97 98 int main() 99 { 100 // freopen("a.in" , "r" , stdin); 101 int n; 102 while(scanf("%d" , &n) , n) 103 { 104 for(int i=0 ; i<n ; i++) 105 scanf("%lf%lf" , &po[i].x , &po[i].y); 106 107 for(int i=0 ; i<n ; i++) vec[i] = Normal(po[(i+1)%n]-po[i]); 108 double l=0 , r=20000; 109 while(r-l>eps){ 110 double m=(l+r)/2; 111 for(int i=0 ; i<n ; i++) L[i] = Line(po[i]+vec[i]*m , po[(i+1)%n]-po[i]); 112 if(HalfplaneIntersection(L , n , poly)>0) l=m; 113 else r=m; 114 } 115 printf("%.6f\n" , l); 116 } 117 return 0; 118 }