注册 登录  
 加关注

网易博客网站关停、迁移的公告:

将从2018年11月30日00:00起正式停止网易博客运营
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

zxj015

 
 
 

日志

 
 

poj 3525 Most Distant Point from the Sea 多边形内切圆  

2011-05-07 22:24:44|  分类: 计算几何 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
求多边形最大内接圆半径。
半平面交+二分查找:把多边形的每条边向内推R的距离,若变幻之后的仍是一个多边形,说明距离 R太小了,若无边说明R过大。
这边有两个纯几何上的问题,就是如果求直线和平移直线。
一、已知两点求直线ax+by+c=0:
        已知两点为A(x[i],y[i]),B(x[i+1],y[i+1]),取直线上的另外一点C(x,y),则有:
        (x-x[i])/(y-y[i])=(x[i+1]-x[i])/(y[i+1]-y[i]);
         化解得:
        x(y[i+1]-y[i])+y(x[i]-x[i+1])+x[i+1]y[i]-x[i]y[i+1]=0;
       即a=y[i+1]-y[i];
          b=x[i]-x[i+1];
          c=x[i+1]y[i]-x[i]y[i+1];
二、如何平移直线:
就拿最简单的过原点的直线y=(a/b)x的直线来讲吧(有常数项的直线和它类似)。
该直线的倾斜角tan∠T=a/b;
                                cos∠T=b/sqrt(a*a+b*b);
                               sin∠T=a/sqrt(a*a+b*b);
         设直线平移了mid的距离,直线上的点(x,y)变幻为(x-x1,y-y1);
        x1,y1,mid三边可组成一个以mid为斜边的直角三角形;
        且x1=sin∠T*mid=a/sqrt(a*a+b*b)*mid;
           y1=cos∠T*mdi=b/sqrt(a*a+b*b)*mid;
这样取原直线上的任意两点做如上的变幻,以新得到的这两点再做直线即可。
        
#include<iostream>
#include<cmath>
using namespace std;
const double pi=acos(-1.0);
const int N=500;
const double eps=1e-9;
int n,m,s;
struct point
{
 double x,y;
}p[N],q[N],input[N];
int sig(double k)
{
 return (k<-eps)?-1:(k>eps);
}
double dis(point a,point b)
{
 return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
double xmule(point p1,point p2,point p3)
{
 return ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y));
}
void interect(point x,point y,double a,double b,double c)//求2直线交点(xy与ax0+by0+c=0的交点)
{
 double u=fabs(a*x.x+b*x.y+c);
 double v=fabs(a*y.x+b*y.y+c);
 q[++s].x=(x.x*v+y.x*u)/(u+v);
 q[s].y=(x.y*v+y.y*u)/(u+v);
}
void cut(double a,double b,double c)//半平面切割
{
 s=0;
 int i;
 for(i=1;i<=m;i++)//枚举所有点
 {
  if(sig(a*p[i].x+b*p[i].y+c)>=0)//逆时针时为>=0
  {
   q[++s]=p[i];//符合
  }
  else
  {
   if(sig(a*p[i-1].x+b*p[i-1].y+c)>0)
    interect(p[i-1],p[i],a,b,c);//p[i]p[i-1]与直线的交点符合
   if(sig(a*p[i+1].x+b*p[i+1].y+c)>0)
    interect(p[i],p[i+1],a,b,c);//p[i]p[i+1]与直线的交点符合
  }
 }
 for(i=1;i<=s;i++)
 {
  p[i]=q[i];//更新
 }
 p[0]=p[s];
 p[s+1]=p[1];
 m=s;
}
int main()
{
 //freopen("in.txt","r",stdin);
 int i,j;
 double mid;
  while(scanf("%d",&n)!=EOF&&n)
 {
  for(i=1;i<=n;i++)
  {
   scanf("%lf%lf",&input[i].x,&input[i].y);
   p[i]=input[i];
  }
  input[n+1]=input[1];
  p[0]=p[n];
  p[n+1]=p[1];
  double max=-1;
  for(i=1;i<=n;i++)
  {
    for(j=i+1;j<=n;j++)
   {
    double temp=dis(p[i],p[j]);
    if(max<temp)//找2点距离最大的点
    {
     max=temp;
    }
   }
  }
  double l=0, r=max/2;
  while(r-l>eps)//二分半径
  {
   mid=(l+r)/2;
   double a,b,c;
   point tt,ta,tb;
   for(i=1;i<=n;i++)
   {
    p[i]=input[i];
   }
      p[0]=p[n];
      p[n+1]=p[1];
   m=n;
   for(i=1;i<=n;i++)
   {
    tt.x=input[i].y-input[i+1].y;//逆时针,注意顺序,每条边一起向“内”推进R
    tt.y=input[i+1].x-input[i].x;
    double k=mid/sqrt(tt.x*tt.x+tt.y*tt.y);
    tt.x=tt.x*k;
    tt.y=tt.y*k;
    ta.x=input[i].x+tt.x;
    ta.y=input[i].y+tt.y;
    tb.x=input[i+1].x+tt.x;
    tb.y=input[i+1].y+tt.y;

    a=ta.y-tb.y;//由2点求出直线ax+by+c=0
    b=tb.x-ta.x;
    c=ta.x*tb.y-ta.y*tb.x;
    cut(a,b,c);
   }
   if(m<=2) r=mid;
  // else if(fabs(xmule(p[1],p[2],p[3]))<eps) r=mid;
   else l=mid;
  }
  printf("%.6lf\n",l);
 }
 return 0;
}

  评论这张
 
阅读(332)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018