判断点是否在多边形内的5种方法.docx
《判断点是否在多边形内的5种方法.docx》由会员分享,可在线阅读,更多相关《判断点是否在多边形内的5种方法.docx(11页珍藏版)》请在冰豆网上搜索。
判断点是否在多边形内的5种方法
判断点是否在多边形内的5种方法
判断点是否在多边形内(凸包和任意多边形分类讨论)
/*
POJ1548:
判断是否为凸包,判断点(圆是否在凸包内),其中判定点是否在多边形内是主要部分
SampleInput
51.51.52.0
1.01.0
2.02.0
1.752.0
1.03.0
0.02.0
51.51.52.0
1.01.0
2.02.0
1.752.5
1.03.0
0.02.0
1
SampleOutput
HOLEISILL-FORMED
PEGWILLNOTFIT
*/
//法1、2:
叉积判定、面积法判定(适用于凸包)。
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
usingnamespacestd;
#definemaxn10005
doublex_multi(pointp1,pointp2,pointp3)
{
return(p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y);
}
boolpoint_is_inside()//叉积判断点在凸包内部!
只针对于凸多边形。
圆心连接每一条边的端点得到的叉积必须同向。
以此可以延伸出面积法判定点是否在凸包内部。
这两种方法都局限于在凸多边形
{
pointp1;
p1.x=pegx,p1.y=pegy;
inti,flag=1;
doubletmp1=0.0,tmp2;
for(i=0;i{
tmp2=Fabs(x_multi(p1,p[i],p[(i+1)%n]));
if(tmp1*tmp2<-eps)
{
flag=0;
break;
}
tmp1=tmp2;
}
returnflag;
}
doubleLen_ab(pointp1,pointp2)
{
returnsqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}
boolcircle_is_inside()//判断圆是否在凸包内
{
if(pegr==0.0)
returntrue;
inti;
doubleans;
pointpeg,t;
peg.x=pegx,peg.y=pegy;
for(i=0;i{
t=peg;
t.x+=p[i].y-p[(i+1)%n].y;
t.y+=p[(i+1)%n].x-p[i].x;
if(Fabs(x_multi(peg,t,p[i])*x_multi(peg,t,p[(i+1)%n]))==1)//如果垂足不在线段上则选择到线段两个端点距离较小的
ans=Fabs(Len_ab(peg,p[i])-Len_ab(peg,p[(i+1)%n]))==-1?
Len_ab(peg,p[i]):
Len_ab(peg,p[(i+1)%n]);
else
ans=fabs(x_multi(peg,p[(i+1)%n],p[i]))/Len_ab(p[i],p[(i+1)%n]);//否则利用面积/底边得到最小距离
if(ans-pegr<-eps)
returnfalse;
}
returntrue;
}
intmain()
{
inti,j;
while(scanf("%d",&n)&&n>=3)
{
scanf("%lf%lf%lf",&pegr,&pegx,&pegy);
for(i=0;iscanf("%lf%lf",&p[i].x,&p[i].y);
doubletmp1=0.0,tmp2;
boolflag=true;
for(i=0;i{
tmp2=Fabs(x_multi(p[i],p[(i+1)%n],p[(i+2)%n]));//精度控制,否则一直wrong
if(tmp1*tmp2<-eps)
{
flag=false;
break;
}
tmp1=tmp2;
}
if(!
flag)
{
puts("HOLEISILL-FORMED");
continue;
}
if(!
point_is_inside()||!
circle_is_inside())//判断圆是否在凸多边形内
{
puts("PEGWILLNOTFIT");
continue;
}
puts("PEGWILLFIT");
}
return0;
}
//法3:
射线法判定点是否在多边形内部(适用于任意多边形):
做一条水平射线计算与多边形的交点个数num,如果num&1则表示在多边形内,否则在多边形外。
其中射线正好交与多边形端点或者与多边形的边平行需要特判。
。
boolOnsegment(pointp1,pointp2,pointp3)
{
doublemin_x=min(p1.x,p2.x);
doublemin_y=min(p1.y,p2.y);
doublemax_x=max(p1.x,p2.x);
doublemax_y=max(p1.y,p2.y);
if(p3.x>=min_x&&p3.x<=max_x&&p3.y>=min_y&&p3.y<=max_y)
returntrue;
returnfalse;
}
boolIs_intersected(pointp1,pointp2,pointp3,pointp4)//线段相交
{
doubled1=x_multi(p1,p2,p3);
doubled2=x_multi(p1,p2,p4);
doubled3=x_multi(p3,p4,p1);
doubled4=x_multi(p3,p4,p2);
if(d1*d2<0.0&&d3*d4<0.0)
returntrue;
//if(d1==0.0&&Onsegment(p1,p2,p3))//由于前面的特判,低处的交点不作为计算
//returntrue;
if(d2==0.0&&Onsegment(p1,p2,p4))
returntrue;
if(d3==0.0&&Onsegment(p3,p4,p1))
returntrue;
if(d4==0.0&&Onsegment(p3,p4,p2))
returntrue;
returnfalse;
}
doubleDot(pointp1,pointp2,pointp3)//点积
{
return(p2.x-p1.x)*(p3.x-p1.x)+(p2.y-p1.y)*(p3.y-p1.y);
}
intpointonsegment(pointp0,pointp1,pointp2)//判断点是否在线段上
{
returnFabs(x_multi(p0,p1,p2))==0&&Fabs(Dot(p0,p1,p2))<=0;
}
boolpoint_is_inside()
{
inti,num=0;
pointp1,peg,p2,p3;
p1.x=999999999.0,p1.y=pegy,peg.x=pegx,peg.y=pegy;//p1坐为在射线极远处的一个点,可以将射线看做线段
for(i=0;i{
if(p[i].y==p[(i+1)%n].y)//如果和多边形的边平行,则判断起点是否在多边形的该边上,避免了和边重合算作无数多个点
{
if(pointonsegment(peg,p[i],p[(i+1)%n]))//判断点是否在线段上
returntrue;
}
else
{
p2=p[i],p3=p[(i+1)%n];
if(p2.y>p3.y)//画图知在与多边形端点相交的时候直接计算交的次数都无法直接判定是否在多边形的内外。
一条线段的端点有高低之分,此时规定高点的交点为有效交点
swap(p2,p3);
if(Is_intersected(peg,p1,p2,p3))
num++;
}
}
returnnum%2==1;
}
//法4:
角度和判定法,适用于任意多边形。
如果点在多边形内则点连接多边形每条边的到的角度*叉积和为360.在边上的话为180.
doublex_multi(pointp1,pointp2,pointp3)
{
return(p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y);
}
doubleLen_ab(pointp1,pointp2)
{
returnsqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}
doubleDot(pointp1,pointp2,pointp3)
{
return(p2.x-p1.x)*(p3.x-p1.x)+(p2.y-p1.y)*(p3.y-p1.y);
}
intpointonsegment(pointp1,pointp2,pointp3)
{
returnFabs(x_multi(p1,p2,p3))==0&&Fabs(Dot(p1,p2,p3))<=0;
}
boolpoint_is_inside()
{
inti;
doublesum=0.0;
for(i=0;i{
if(p0==p[i])//点在多边形端点上
returntrue;
if(p[i]==p[(i+1)%n])//去重点
continue;
if(pointonsegment(p0,p[i],p[(i+1)%n]))//点在多边形边上
returntrue;
doublea=Len_ab(p0,p[i]);
doubleb=Len_ab(p0,p[(i+1)%n]);
doublec=Len_ab(p[i],p[(i+1)%n]);
sum+=Fabs(x_multi(p0,p[i],p[(i+1)%n]))*acos((a*a+b*b-c*c)/(2.0*a*b));//计算角度和,叉积大于0则加上,小于0则减去
}
sum=fabs(sum);
if(Fabs(sum-2.0*pi)==0)
returntrue;
returnfalse;
}
/法5:
改进弧长法(适用于任意多边形)——权威算法!
精度很高!
以下对于该算法的解析摘自
该算法只需O
(1)的附加空间,时间复杂度为O(n),系数很小;最大的优点是具有很高的精度,只需做乘法和减法,若针对整数坐标则完全没有精度问题.而且实现起来也非常简单,比转角法和射线法都要好写且不易出错.
有关"弧长法"的介绍:
"弧长法要求多边形是有向多边形,一般规定沿多边形的正向,边的左侧为多边形的内侧域.以被测点为圆心作单位圆,将全部有向边向单位圆作径向投影,并计算其中单位圆上弧长的代数和,若代数和为0,则点在多边形外部;若代数和为2π,则点在多边形内部;若代数和为π,则点在多边形上."
根据上面的介绍,其实弧长法就是转角法,但它的改进方法比较比较厉害:
将坐标原点平移到被测点P,这个新坐标系将平面划分为4个象限,对每个多边形顶点P,只考虑其所在的象限,然后按邻接顺序访问多边形的各个顶点P,分析P[i]和P[i+1],有下列三种情况:
(1)P[i+1]在P[i]的下一象限,此时弧长和加π/2;
(2)P[i+1]在P[i]的上一象限,此时弧长和减π/2;
(3)P[i+1]在P[i]的相对象限,首先计算f=p[i+1].y*p[i].x-p[i+1].x*p[i].y(叉积),若f=0,则点在多边形上;若f<0,弧长和减π;若f>0,弧长和加π.最后对算出的代数和和上述的情况一样判断即可.
实现的时候还有两点要注意:
1>若P的某个坐标为0时,一律当正号处理;
2>若被测点和多边形的顶点重合时要特殊处理.
还有一个问题那就是当多边形的某条边在坐标轴上而且两个顶点分别在原点的两侧时会出错,如边(3,0)……>(-3,0),按以上的处理,象限分别是第一和第二,这样会使代数和加π/2,有可能导致最后结果是被测点在多边形外,而实际上被测点是在多边形上(该边穿过该点).
对于这点,处理办法是:
每次计算P[i]和P[i+1]时,就计算叉积和点积,判断该点是否在该边上,是则判断结束,否则继续上述过程,这样牺牲了时间,但保证了正确性.
具体实现的时候,由于只需知道当前点和上一点的象限位置,所以附加空间只需O
(1).实现的时候可以把上述的"π/2"改成1,"π"改成2,这样便可以完全使用整数进行计算,不必考虑顶点的顺序,逆时针和顺时针都可以处理,只是最后的代数和符号不同而已
intget_tmp(pointp0)
{
returnp0.x>=0?
(p0.y>=0?
0:
3):
(p0.y>=0?
1:
2);
}
boolpoint_is_inside()
{
inttmp1,tmp2,sum=0,i;
pointp0,p1;
p0.x=pegx,p0.y=pegy;
p1.x=p[0].x-p0.x,p1.y=p[0].y-p0.y;
tmp1=get_tmp(p1);
for(i=0;i{
if(p[i]==p0)
break;
intt0=Fabs(x_multi(p0,p[i],p[(i+1)%n]));
intt1=Fabs((p[i].x-p0.x)*(p[(i+1)%n].x-p0.x));
intt2=Fabs((p[i].y-p0.y)*(p[(i+1)%n].y-p0.y));
if(!
t0&&t1<=0&&t2<=0)//被测点在多边形边上
break;
p1.x=p[(i+1)%n].x-p0.x,p1.y=p[(i+1)%n].y-p0.y;
tmp2=get_tmp(p1);//计算象限
switch((tmp2-tmp1+4)%4)
{
case1:
{sum++;break;}
case2:
{
if(t0>0)sum+=2;
elsesum-=2;
break;
}
case3:
{sum--;break;}
}
tmp1=tmp2;
}
if(ireturntrue;
returnfalse;
}