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 }
 



原文地址:https://www.cnblogs.com/zxj015/p/2740210.html