poj 3384 Feng Shui 半平面交
给定一个多边形，在多边形内放有两个相同的圆，使两个圆尽可能多的覆盖多边形。输出最终两个圆心的位置。
最优的放置方法必定是圆内切于两条边，那么把所有的边向内推移半径R的距离，得到新的多边形（也有可能是一个点或一条直线），求出新多边形相距最远的两个顶点，这两个顶点就是圆心的位置了。
1 #include<stdio.h> 2 #include<stdlib.h> 3 #include<iostream> 4 #include<cmath> 5 using namespace std; 6 const double pi=acos(-1.0); 7 const int N=105; 8 const double eps=1e-9; 9 int n,m,s; 10 struct point 11 { 12 double x,y; 13 }p[N],q[N],input[N]; 14 int sig(double k) 15 { 16 return (k<-eps)?-1:(k>eps); 17 } 18 double dis(point a,point b) 19 { 20 return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); 21 } 22 double xmule(point p1,point p2,point p3) 23 { 24 return ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y)); 25 } 26 void interect(point x,point y,double a,double b,double c)//求2直线交点(xy与ax0+by0+c=0的交点) 27 { 28 double u=fabs(a*x.x+b*x.y+c); 29 double v=fabs(a*y.x+b*y.y+c); 30 q[++s].x=(x.x*v+y.x*u)/(u+v); 31 q[s].y=(x.y*v+y.y*u)/(u+v); 32 } 33 void cut(double a,double b,double c)//半平面切割 34 { 35 s=0; 36 int i; 37 for(i=1;i<=m;i++)//枚举所有点 38 { 39 if(sig(a*p[i].x+b*p[i].y+c)<=0) 40 { 41 q[++s]=p[i];//符合 42 } 43 else 44 { 45 if(sig(a*p[i-1].x+b*p[i-1].y+c)<0) 46 interect(p[i-1],p[i],a,b,c);//p[i]p[i-1]与直线的交点符合 47 if(sig(a*p[i+1].x+b*p[i+1].y+c)<0) 48 interect(p[i],p[i+1],a,b,c);//p[i]p[i+1]与直线的交点符合 49 } 50 } 51 for(i=1;i<=s;i++) 52 { 53 p[i]=q[i];//更新 54 } 55 p[0]=p[s]; 56 p[s+1]=p[1]; 57 m=s; 58 } 59 double dir(point a,point b) 60 { 61 return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)); 62 } 63 int main() 64 { 65 int i,j,ans1,ans2,r; 66 double mid,max,len; 67 while(scanf("%d%d",&n,&r)!=EOF) 68 { 69 for(i=1;i<=n;i++) 70 { 71 scanf("%lf%lf",&input[i].x,&input[i].y); 72 p[i]=input[i]; 73 } 74 input[n+1]=input[1]; 75 p[0]=p[n]; 76 p[n+1]=p[1]; 77 double a,b,c; 78 point tt,ta,tb; 79 for(i=1;i<=n;i++) 80 { 81 p[i]=input[i]; 82 } 83 p[0]=p[n]; 84 p[n+1]=p[1]; 85 m=n; 86 for(i=1;i<=n;i++) 87 { 88 tt.x=input[i+1].y-input[i].y;//顺时针，注意顺序，每条边一起向“内”推进R 89 tt.y=input[i].x-input[i+1].x; 90 double k=r*1.0/sqrt(tt.x*tt.x+tt.y*tt.y); 91 tt.x=tt.x*k; 92 tt.y=tt.y*k; 93 ta.x=input[i].x+tt.x; 94 ta.y=input[i].y+tt.y; 95 tb.x=input[i+1].x+tt.x; 96 tb.y=input[i+1].y+tt.y; 97 98 a=ta.y-tb.y;//由2点求出直线ax+by+c=0 99 b=tb.x-ta.x; 100 c=ta.x*tb.y-ta.y*tb.x; 101 cut(a,b,c); 102 } 103 max=-1; 104 ans1=m; 105 ans2=m; 106 for(i=1;i<m;i++) 107 for(j=i+1;j<=m;j++) 108 { 109 len=dir(p[i],p[j]); 110 if(len-max>eps) 111 { 112 max=len; 113 ans1=i; 114 ans2=j; 115 } 116 } 117 printf("%.4f %.4f %.4f %.4f\n",p[ans1].x,p[ans1].y,p[ans2].x,p[ans2].y); 118 } 119 return 0; 120 }
