NKoj 2118 Handy Service（计算几何）

题目链接：http://acm.nankai.edu.cn/p2118.html

 

题 意：给出一个多边形和多边形外两点A，B，求在不穿过多边形的条件下，两点A，B间的最短距离。

 

最短路问题，可以在能够直接相连的顶点间连边，然后求A点到B的最短路。

建图的关键在于判断两点的连线是否穿越了多边形。

 

有这样一个结论：如果线段穿越多边形，则1：线段与多边形中的边有内交关系

或2： 在所有与顶点相交的点中的相邻交点连线在多边形内。

 

http://bbs.pep.com.cn/thread-241056-1-1.html

[算法]判断线段是否在多边形内判断线段是否在多边形内线段在多边形内的一个必要条件是线段的两个端点都在多边形内；如果线段和多边形的某条边内交（两线段内交是指两线段相交且交点不在两线段的端点），因为多边形的边的左右两侧分属多边形内外不同部分，所以线段一定会有一部分在多边形外。于是我们得到线段在多边形内的第二个必要条件：线段和多边形的所有边都不内交；线段和多边形交于线段的两端点并不会影响线段是否在多边形内；但是如果多边形的某个顶点和线段相交，还必须判断两相邻交点之间的线段是否包含与多边形内部。因此我们可以先求出所有和线段相交的多边形的顶点，然后按照X-Y坐标排序，这样相邻的两个点就是在线段上相邻的两交点，如果任意相邻两点的中点也在多边形内，则该线段一定在多边形内。证明如下：命题1：如果线段和多边形的两相邻交点P1 ，P2的中点P' 也在多边形内，则P1, P2之间的所有点都在多边形内。证明：假设P1,P2之间含有不在多边形内的点，不妨设该点为Q，在P1, P'之间，因为多边形是闭合曲线，所以其内外部之间有界，而P1属于多边行内部，Q属于多边性外部，P'属于多边性内部，P1-Q-P'完全连续，所以P1Q和QP'一定跨越多边形的边界，因此在P1,P'之间至少还有两个该线段和多边形的交点，这和P1P2是相邻两交点矛盾，故命题成立。证毕由命题1直接可得出推论：推论2：设多边形和线段PQ的交点依次为P1,P2,……Pn，其中Pi和Pi+1是相邻两交点，线段PQ在多边形内的充要条件是：P，Q在多边形内且对于i =1, 2,……, n-1，Pi ,Pi+1的中点也在多边形内。在实际编程中，没有必要计算所有的交点，首先应判断线段和多边形的边是否内交，倘若线段和多边形的某条边内交则线段一定在多边形外；如果线段和多边形的每一条边都不内交，则线段和多边形的交点一定是线段的端点或者多边形的顶点，只要判断点是否在线段上就可以了。至此我们得出算法如下：1. if 线端PQ的端点不都在多边形内2. then return false;3. 点集pointSet初始化为空;4. for 多边形的每条边s5. do if 线段的某个端点在s上6. then 将该端点加入pointSet;7. else if s的某个端点在线段PQ上8. then 将该端点加入pointSet;9. else if s和线段PQ相交 // 这时候可以肯定是内交10. then return false;11. 将pointSet中的点按照X-Y坐标排序，X坐标小的排在前面， 对于X坐标相同的点，Y坐标小的排在前面；12. for pointSet中每两个相邻点 pointSet , pointSet[ i+1]13. do if pointSet , pointSet[ i+1] 的中点不在多边形中14. then return false;15. return true;这个算法的复杂度也是O(n)。其中的排序因为交点数目肯定远小于多边形的顶点数目n，所以最多是常数级的复杂度，几乎可以忽略不计。

另判断点是否在多边形内：

http://www.cppblog.com/w2001/archive/2008/09/23/31694.html

1. 叉乘判别法（只适用于凸多边形）想象一个凸多边形，其每一个边都将整个2D屏幕划分成为左右两边，连接每一边的第一个端点和要测试的点得到一个矢量v，将两个2维矢量扩展成3维的，然后将该边与v叉乘，判断结果3维矢量中Z分量的符号是否发生变化，进而推导出点是否处于凸多边形内外。这里要注意的是，多边形顶点究竟是左手序还是右手序，这对具体判断方式有影响。2. 面积判别法（只适用于凸多边形）第四点分别与三角形的两个点组成的面积分别设为S1,S2,S3，只要S1+S2+S3>原来的三角形面积就不在三角形范围中.可以使用海伦公式 。推广一下是否可以得到面向凸多边形的算法？（不确定）3. 角度和判别法（适用于任意多边形）double angle = 0;realPointList::iterator iter1 = points.begin();for (realPointList::iterator iter2 = (iter1 + 1); iter2 < points.end(); ++iter1, ++iter2) { double x1 = (*iter1).x - p.x; double y1 = (*iter1).y - p.y; double x2 = (*iter2).x - p.x; double y2 = (*iter2).y - p.y; angle += angle2D(x1, y1, x2, y2); }if (fabs(angle - span::PI2) < 0.01) return true;else return false;另外，可以使用bounding box来加速。if (p.x < (*iter)->boundingBox.left || p.x > (*iter)->boundingBox.right || p.y < (*iter)->boundingBox.bottom || p.y > (*iter)->boundingBox.top) 。。。。。。对于多边形来说，计算bounding box非常的简单。只需要把水平和垂直方向上的最大最小值找出来就可以了。对于三角形：第四点分别与三角形的两个点的交线组成的角度分别设为j1,j2,j3，只要j1+j2+j3>360就不在三角形范围中。4. 水平/垂直交叉点数判别法（适用于任意多边形）注意到如果从P作水平向左的射线的话，如果P在多边形内部，那么这条射线与多边形的交点必为奇数，如果P在多边形外部，则交点个数必为偶数（0也在内）。所以，我们可以顺序考虑多边形的每条边，求出交点的总个数。还有一些特殊情况要考虑。假如考虑边(P1,P2)，1)如果射线正好穿过P1或者P2,那么这个交点会被算作2次，处理办法是如果P的从坐标与P1,P2中较小的纵坐标相同，则直接忽略这种情况2)如果射线水平，则射线要么与其无交点，要么有无数个，这种情况也直接忽略。3)如果射线竖直，而P0的横坐标小于P1,P2的横坐标，则必然相交。4)再判断相交之前，先判断P是否在边(P1,P2)的上面，如果在，则直接得出结论：P再多边形内部。

 

代码：

#include<stdio.h>#include<stdlib.h>#include<string.h>#include<math.h>#define Max(x,y) (x>y?x:y)#define Min(x,y) (x>y?y:x)#define Mid(x,y) ((x+y)*0.5)#define M 32const double inf=10e12;const double eps=1e-8;const double Pi=acos(-1);typedef struct talP{double x,y;}Point;Point p[M];double dis[M][M];int n;int tmp[M];int tn;double Dis(Point p,Point q){p.x-=q.x,p.y-=q.y;return sqrt(p.x*p.x+p.y*p.y);}double XMulti(Point p,Point q,Point o){return (p.x-o.x)*(q.y-o.y)-(p.y-o.y)*(q.x-o.x);}bool PInLine(Point t,Point p,Point q){return t.x>=Min(p.x,q.x)&&t.x<=Max(p.x,q.x)&&t.y>=Min(p.y,q.y)&&t.y<=Max(p.y,q.y);}int Cmp(const void *a,const void *b){return *(int *)a-*(int *)b;}int TCmp(double x){return (x>eps)-(x<-eps);}bool InConv(Point o){int i,l,r;double x1,x2,y1,y2;double tmp,angl=0;for(i=1;i<=n;i++){l=i;r=i%n+1;x1=p[l].x-o.x;y1=p[l].y-o.y;x2=p[r].x-o.x;y2=p[r].y-o.y;if(TCmp(x1*y2-x2*y1)==0&&o.x>=Min(p[l].x,p[r].x)&&o.x<=Max(p[l].x,p[r].x)&&o.y>=Min(p[l].y,p[r].y)&&o.y<=Max(p[l].y,p[r].y))return false;}for(i=1;i<=n;i++){l=i;r=i%n+1;x1=p[l].x-o.x;y1=p[l].y-o.y;x2=p[r].x-o.x;y2=p[r].y-o.y;tmp=(x1*x2+y1*y2)/Dis(p[l],o)/Dis(p[r],o);angl+=TCmp(x1*y2-x2*y1)*acos(tmp);}return (fabs(angl)>2*Pi-eps);}bool SegCross(int f,int t){int l,r,i;Point c;double r1,r2,r3,r4;tn=0;for(i=1;i<=n;i++){l=i,r=i%n+1;r1=XMulti(p[l],p[f],p[t]);r2=XMulti(p[r],p[f],p[t]);r3=XMulti(p[f],p[r],p[l]);r4=XMulti(p[t],p[r],p[l]);if(r1*r2<0&&r3*r4<0)return true;if(r1==0&&PInLine(p[l],p[f],p[t]))tmp[tn++]=l;if(r2==0&&PInLine(p[r],p[f],p[t]))tmp[tn++]=r;}qsort(tmp,tn,sizeof(tmp[0]),Cmp);tmp[tn]=tmp[0];for(i=0;i<tn;i++){if(tmp[i]==tmp[i+1]) continue;c.x=Mid(p[tmp[i]].x,p[tmp[i+1]].x);c.y=Mid(p[tmp[i]].y,p[tmp[i+1]].y);if(InConv(c)) return true;}return false;}void Init(){int i,j;for(i=0;i<=n+1;i++){for(j=0;j<=n+1;j++)dis[i][j]=inf;}for(i=0;i<=n+1;i++) dis[i][i]=0;for(i=0;i<=n;i++){for(j=i+1;j<=n+1;j++){if(!SegCross(i,j))dis[i][j]=dis[j][i]=Dis(p[i],p[j]);}}}void Floyd(){int i,j,k;for(k=0;k<=n+1;k++){for(i=0;i<=n+1;i++){if(dis[i][k]==inf||i==k) continue;for(j=0;j<=n+1;j++){if(i==j||dis[k][j]==inf) continue;if(dis[i][j]-dis[k][j]>dis[i][k])dis[i][j]=dis[k][j]+dis[i][k];}}}}int main(){int i;scanf("%d",&n);for(i=1;i<=n;i++)scanf("%lf%lf",&p[i].x,&p[i].y);scanf("%lf%lf%lf%lf",&p[0].x,&p[0].y,&p[n+1].x,&p[n+1].y);Init();Floyd();printf("%.3f/n",dis[0][n+1]);return 0;}