5 计算几何 计算几何是计算机科学中的一个分支，是专门研究有关几何对象问题的。它在图像分析、模式识别、计算机图形学等中应用甚广。本章主要介绍几个基本计算几何问题的简单并行算法和它们的MPI编程实现，包括包含问题、相交问题和凸壳问题等。 包 含 问 题 包含问题(Inclusion Problem)是计算几何中最基本的问题之一。简单的来说，包含问题就是判断点与多边形的位置关系，即点在多边形内，点在多边形外，和点在多边形上。我们这里讨论的仅是点与任意多边形之间的关系，不考虑与任意曲线封闭图形之间的关系。 一个能够画在平面图上而没有任何相交边的图形称为平面图（Planar Graph）。由一些相邻的多边形（称为单图）所组成的图形称为平面细图（Planar Subvision），其中各多边形中除了顶点外，无两边相交。给定一个平面细图和一个顶点p，确定哪个多边形包含了p，此即所谓的包含问题。最简单的情况是判断点p是否在一个多边形Q中。 包含问题及其串行算法 判断点在多边形中的算法的基本思想是先求点和边的交点个数，然后根据交点个数确定点和边的关系。基本步骤是：①过p点向x轴的负半轴方向做一条射线；②求此射线与多边型Q诸边的交点；③判断：交点的个数若为奇数，则p位于Q内；否则p位于Q外。这种测试，对于n边形可在O(n)步内完成。很清楚，这是最优的串行算法，因为读入n条边也至少需要Ω(n)步。 但是需要注意几种特殊的边界情况：首先，过点p的射线与多边形顶点相交，则交点只计数一次；其次，若点p在多边形边上，则也认为点p在多边形中；最后，如果射线与水平边重叠且点p不在该边上，则交点不计数。 算法17.1 单处理机上包含问题算法 输入：点p坐标(px,py)；多边形Q的n条边E1,E2, ,…,En的两个端点坐标集合S 输出：点和多边形的关系：true(p位于Q内)；false(p位于Q外) Begin (1) count=0 (2) while i<n do if p is on Ei then return true else if y=py(x≤px) intersects Ei left q and q is not Ei low-end then count=count+1 end if end if end while (3) if ((count mod 2)=1) then return true else return false end if End 算法的总运行时间为t(n)=O(n)。 包含问题并行算法 我们介绍两种简单方法来实现上面串行算法的并行化。一种方法是将串行算法中的步骤2并行化。假设系统有p个处理器，多边形有n条边，将n条边平均分配到这p个处理器中，每个处理器最多处理条边。具体的处理方法为过p点向x的负半轴方向做一条射线，判断：如果点是在边上则返回为；线与边无交点时返回0；如果有交点，那么就有2种情况，如果交点是边的下顶点，则返回为0，否则，即交点在边上，则返回为1；把一个处理器中的所有返回值加起来，然后将该值发送到主处理器(my_rank=0)，最后主处理器根据点的个数来判断点与多边形的关系。 另一种方法同样也是将串行算法中的步骤2并行化。假设有p=2α-1个处理器，为正整数。p个处理器的编号从根开始自上而下，自左而右逐级向下推进。每个处理器存储多边形Q的一条边，边由其两个端点的迪卡尔坐标表示，点p的坐标为()。开始时，根读进，然后传送给树中的其余处理器。当接收到p的坐标时，他就确定：穿过p的射线是否和Q的边相交；对于特殊的情况需要利用图形学的知识处理之。如果条件满足，则产生“1”输出；否则为“0”。将个处理器的输出相加，若和为奇数，则p位于Q内；如果为偶数，则p位于Q外；如果和为，则输出p在多边形上。 算法17.2 多处理机上包含问题算法 输入：点p坐标(px,py)；多边形Q的n条边E1,E2, ,…,En的两个端点坐标集合S 输出：点和多边形的位置关系：true(p位于Q内)；false(p位于Q外) Begin (1) count=0 (2) for all Pj where 1≤j≤p par-do resultj=0 end for (3) for all Pj where 1≤j≤p par-do for i=1 to  do if p is on Ei×p+j then return resultj= else if y=py(x≤px) intersects Ei×p+j left q and q is not Ei×p+j low-end then resultj= resultj+1 end if end if end for end for (4) for j=1 to p par-do Pj sends resultj to P0 end for for j=1 to p do count= count+ resultj end for if ((count mod 2)=1) or (count=) then return true else return false end if End 在有p个处理器的情况下，上述算法的时间复杂度为O()。 MPI源程序请参见所附光盘。 相交问题 在很多实际应用中，要求确定一组几何物体(目标)是否相交。例如，在模式分类中，必须确定代表不同类别的空间中的不同区域是否具有共同的子区域；在集成电路设计中，重要的是要避免导线的交叉和元件的重叠；在计算机图形学中，要求消去三维景象的二维表示中的隐线和隐面等等。像如上的这些问题都可归结为物体的相交问题(Intersection Problem)。 设有平面上的两个多边形(允许有边相交)R和Q，如果多边形R的一条边和Q的一条边相交，则称R和Q是相交的。所以两个多边形的相交问题可以转化为线段与多边形的相交问题。三维空间的相交问题与二维平面上的相交问题并没有实质的区别，只是在判断边的相交时比二维问题上判断边的相交要麻烦，因为三维空间上的点坐标是与3个值有关的。 下面描述的算法都是针对二维平面上的多边形相交而言的。 两多边形相交问题及其串行算法 最基本的相交问题是判断两条线段是否相交。而边与多边形相交就是判断一条边和多条边中的某条边是否相交的算法。要是一个多边形的某条边与另一个多边形的一条边相交，则就称两个多边形相交。这样两个多边形相交的问题就转化为多条边与一个多边形相交的问题。 判断线段是否相交的关键是求两条直线的交点，即判断此交点是否在线段上。这和包含问题有些相似，需要判断点与边的关系。不同点是只须判断点是否在边上。 算法17.3 单处理机上的两多边形相交问题算法 输入：多边形R的n条边E1,E2,,…,En的两个端点坐标集合S1，多边形Q的m条边 F1,F2,,…,Fm的两个端点坐标集合S2 输出：两个多边形是否相交：true(两多边形相交)；false(两多边形不相交) Begin for i=1 to n do for j=1 to m do if (Ei intersects Fj ) then return true end if end for end for return false End 显然上述算法所需时间为O(mn)。 相交问题的并行算法 下面我们给出两多边形相交问题的朴素并行算法：对于多边形R的每一条边，要确定其是否与多边形Q相交；如果R的边中有一条边与Q相交，那么就可以断定多边形R与Q是相交的。假设R有n条边，Q有m条边，总共有p个处理器P1 ,P2 ,…,Pp。对于R中的每条边依次判断是否与Q相交。 算法17.4 两多边形相交问题的并行算法 输入：多边形R的n条边E1,E2,,…,En的两个端点坐标集合S1，多边形Q的m条边 F1,F2,,…,Fm的两个端点坐标集合S2 输出：两个多边形是否相交：true(两多边形相交)；false(两多边形不相交) Begin (1) for i=1 to m do 将Fi广播到所有处理器上 end for (2) for j=1 to m do 将Ej广播到所有处理器上 end for for all Pk where 1≤k≤p par-do for i=1 to  do for j=1 to m do if (Ei×p+k intersects Fj ) then resultk=true end if end for end for end for (4) 将各个处理器上resultk返回到主处理器，如果其中有一个为真，则两多边形相交，否则两多边形不相交 End MPI源程序请参见章末附录。 凸壳问题 凸壳（Convex Hull）问题是计算几何中一个重要问题，它的应用很广泛。例如，在图象处理中可以通过构造凸壳找到数字图象中的凹面；在模式识别中，可视模式的凸壳能够作为描述模式外形的重要特征；在分类中，一组物体的凸壳就可勾画出这些物体的所属的类；在计算机图形学中，使用一组点的凸壳可以显示出点簇；在几何问题中，集合S中最远两点就是凸壳的顶点，等等。 凸壳问题及其串行算法 给定平面中的点的集合，所谓S之凸壳简记为CH(S)，就是包含S所有点的最小的凸多边形。实际中，假定平面上的n个点，用n个钉在木板上的图钉表示，用一条橡皮带缠绕着这些图钉，然后放松橡皮带，在撑紧橡皮带所构成的平面图形就是凸多边形。  图 17.1 图示求凸壳的方法 参照图 17.1，串行算法的基本思路如下：①识别极点（Extreme Point）：S中那些最大X坐标(XMAX)，最大Y坐标(YMAX)，最小X坐标(XMIN)，最小Y坐标(YMIN)的那些顶点称为极点；②识别凸边(Hull Edge)：线段是CH(S)的一条凸边，当且仅当S的其余n－2个顶点均位于穿过和的(无限长的)直线的同侧。由此定义可知，和一定是CH(S)的顶点；③识别凸点(Hull Point)：令和是CH(S)上的两连续顶点。假定取为坐标原点，那么在所有S中的点，和相对于X轴所形成的正的或负的夹角是最小。这些点称为凸点。 很明显，极点都是CH(S)的顶点，任何位于由极点所围成的多边形内的点都不是CH(S)的顶点，那些在多边形外的点可以归入由XMAX、XMIN、YMAX、YMIN所围成的四边形与由极点所围成的多边形所形成的四个三角区中。这样问题就归结成为求找这四个三角区中顶点的凸边，再把这些凸边连接起来就可求得CH(S)。 算法17.5 单处理机上求凸壳的串行算法 输入：n点集合 输出：返回包含S的凸壳的顶点表列CH(S) Begin (1) 识别极点，它们都是CH(S)的顶点 (2) 将顶点归入有极点确定的四个三角区中，不在四个三角区中的顶点不需要处理 (3) 计算顶点的极角，并排序： (3.1)依次对属于同一三角区域的点计算极角。然后将有最小极角点的序号作为自己的Nextindex值 (3.2) 然后处理器按照点的Nextindex索引输出点 End 凸壳问题并行算法 假定n个处理器排成网孔形状，处于第i行和第j列的处理器用P(i,j)表示。点i的坐标用表示。首先来介绍两个基本概念和术语：①拓扑结构的转化：给定的n个处理器，可以认为其是直线拓扑结构。现在将一维拓扑结构改变为2维的，则P(i,j)的i=1,2,3,4，而j=1,2,…,n/4。②行主处理器：行主处理器是一个人为设定的处理器，纯粹是为了处理上的方便。在本算法中取p(i,1)为行主处理器。 算法17.6 凸壳问题并行算法 输入：n点集合 输出：返回凸壳顶点表列CH(S) Begin (1) 计算极点： (1.1) 第1行的行主处理器向第2,3,4行的行主处理器发送顶点坐标 (1.2) 第1、2、3、4行的行主处理器分别计算XMAX、 YMIN、 YMAX、 XMIN； 并把它们的坐标分别存储在4个行主处理器P(1,1)、P(2,1)、P(3,1)、P(4,1)中 (1.3)确定极点后，将四条由极点组成的边存储到每一行的行主处理器上 (2) 确定S中的顶点是否在四个三角区中： (2.1)四个行主处理器同时判断顶点是否处于自身所在的区域，其中属于自己区域内 的点，存储到行主处理器中本行要处理的顶点表列 (2.2)将行主处理器上的表列传递到本行中其余的处理器上 (3) 计算顶点的极角，并排序输出： (3.1)每一行上的第i个处理器计算表列中第j个点极角，其中j是mod行处理器数 等于i的数。然后将有最小极角点的序号作为自己的Nextindex值 (3.2)处理器将计算的Nextindex值传递到每行中的主处理器 (3.3)然后每个行主处理器按照从极点开始按Nextindex索引依次输出所得到的极点 上的点 End MPI源程序请参见所附光盘。 小结 本章主要介绍了计算几何中包含问题、相交问题和凸壳等三个基本问题的并行算法及其MPI编程及其实现。这些算法可参见文献[1]。读者如欲进一步学习可参考文献[2]、[3]和[4]。 参考文献 陈国良 编著. 并行算法的设计与分析（修订版）. 高等教育出版社, 2002.11 Akl S G. Parallel Computational Geometry. Prentic-Hall, 1992 Atallah M J. Goodrich M T. Efficient Parallel Solutions to Some Geometric Problems. J. of Parallel and Distributed Computing, 1986, 3: 492-507 周培德 著. 计算几何-算法设计与分析. 清华大学出版社,广西科学技术出版社, 2000.5 附录 包含问题并行算法的MPI源程序 1. 源程序including.c #include <mpi.h> #include <stdio.h> int n; double xtemp[20],ytemp[20]; double x,y; int s,mys; int group_size,my_rank; /*判断射线与线段是否有交点*/ int cal_inter(int number,int i,double x,double y) { double x1,y1,x2,y2,temp; int result; result=0; if(number+i*group_size>=n) return result; x1=xtemp[number+i*group_size]; y1=ytemp[number+i*group_size]; x2=xtemp[(number+1+i*group_size)%n]; y2=ytemp[(number+1+i*group_size)%n]; if(y1>y2) { temp=x1; x1=x2; x2=temp; temp=y1; y1=y2; y2=temp; } /*判断竖直边的情况*/ if(x1==x2) { if((x>x1)&&(y<=y2)&&(y>y1)) result=1; else result=0; /*点在竖直边上，应该对result赋一个比 较的大的值，这里是100*/ if((x==x1)&&((y-y1)*(y2-y)>=0)) result=100; } else { /*非竖直边，非水平边*/ if (y1!=y2) { temp=x2+(y-y2)*(x2-x1)/(y2-y1); /*交点刚好在边上，且不为下顶 点*/ if((temp<x)&&(y<=y2)&&(y>y1)) result=1; else result=0; /*点在边上，应该对result赋一个 比较的大的值，这里是100*/ if((temp==x)&&((y-y2)* (y1-y)>=0)) result=100; } else { /*点在水平边上，应该对result赋 一个比较的大的值，这里是100*/ if((y==y1)&&((x1-x)*(x-x2)>=0)) result=100; } } return result; } main(int argc,char* argv[]) { int i; MPI_Init(&argc,&argv); MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); MPI_Comm_size(MPI_COMM_WORLD, &group_size); /*各处理器计数器初始化，对应于 算法17.2步骤(2)*/ mys=0; /*主处理器读入多边形顶点和要判断的点 的坐标*/ if(my_rank==0) { printf("请输入点的个数:"); scanf("%d",&n); printf("请输入各点的坐标\n"); for(i=0;i<n;i++) { printf("%d:",i); scanf("%lf",&xtemp[i]); scanf("%lf",&ytemp[i]); } printf("请输入要判断点的坐标\n"); scanf("%lf %lf",&x,&y); } MPI_Barrier(MPI_COMM_WORLD); /*把多边形的顶点数、顶点坐标与要判别的点 的坐标播送给所有进程*/ MPI_Bcast(&n,1,MPI_INT,0, MPI_COMM_WORLD); MPI_Bcast(&x,1,MPI_DOUBLE,0, MPI_COMM_WORLD); MPI_Bcast(&y,1,MPI_DOUBLE,0, MPI_COMM_WORLD); MPI_Bcast(xtemp,n,MPI_DOUBLE,0, MPI_COMM_WORLD); MPI_Bcast(ytemp,n,MPI_DOUBLE,0, MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD); /*每一个处理器处理n/group_size条边上的情 况并求和，对应于算法17.2步骤(3)*/ for(i=0;i<n/group_size+1;i++) { mys+=cal_inter(my_rank,i,x,y); } MPI_Barrier(MPI_COMM_WORLD); /*把mys的值规约到s，对应于 算法17.2步骤(4)和(5)*/ MPI_Reduce(&mys,&s,1,MPI_INT,MPI_SUM, 0,MPI_COMM_WORLD); /*根据s值确定输出结果，对应于 算法17.2步骤(6)*/ if(my_rank==0) { if(s>=100) printf("vertex p is in polygon\n"); else if(s%2==1) printf("vertex p is in polygon\n"); else printf("vertex p is out of polygon\n"); } MPI_Finalize(); } 2. 运行实例 编译：mpicc including.c –o including 运行：可以使用命令 mpirun –np SIZE including来运行该串匹配程序，其中SIZE是所使用的处理器个数。本实例中使用了SIZE=4个处理器。 mpirun –np 4 including 运行结果： 请输入点的个数:6 请输入各点的坐标 0：1 0 1：2 1 2：4 0 3：2 4 4：1 2 5：0 3 请输入要判断点的坐标 2 3 结果：vertex p is in polygon 说明：输入顶点的(x,y)坐标时中间用空格隔开。