GIS算法资料

算法概念

栅格化

    栅格化这个术语可以用于任何将[[矢量图形]]转换成[[位图|栅格图像]]的过程。这个术语用来表示在计算机上显示三维形状的流行[[渲染]]算法。栅格化目前是生成实时[[三维计算机图形]]最流行的算法。实时应用需要立即响应用户输入，并且通常需要至少每秒 24 帧的速率。
    最基础的栅格化算法将多边形表示的三维场景渲染到二维表面。多边形由三角形的集合表示，三角形由三维空间中的三个顶点表示。在最简单的实现形式中，栅格化工具将顶点数据映射到观察者显示器上对应的二维坐标点，然后对变换出的二维三角形进行合适的填充。


变换
通常使用矩阵运算进行变换，另外也可以用四元数运算。在三维顶点中添加一个齐次变量成为四维定点然后左乘一个 4 x 4 的变换矩阵，通过这种方法就可以对三维顶点进行变换。主要的变换有平移、缩放、旋转以及投射。


裁剪
一旦三角形顶点转换到正确的二维位置之后，这些位置可能位于观察窗口之外，也可能位于屏幕之内。裁剪就是对三角形进行处理以适合显示区域的过程。
最常用的技术是Sutherland-Hodgeman裁剪算法。在这种方法中，每次测试每个图像平面的四条边，对于每个边测试每个待渲染的点。如果该点位于边界之外，就剔除该点。对于与图像平的面边相交的三角形边，即边的一个顶点位于图像内部一个位于外部，那么就在交叉点插入一个点并且移除外部的点。


扫描变换
传统的栅格化过程的最后一步就是填充图像平面中的二维三角形，这个过程就是扫描变换。
第一个需要考虑的问题就是是否需要绘制给定的像素。一个需要渲染的像素必须位于三角形内部、必须未被裁掉，并且必须未被其它像素遮挡。有许多算法可以用于在三角形内进行填充，其中最流行的方法是扫描线算法。
由于很难确定栅格化引擎是否会从前到后绘制所有像素，因此必须要有一些方法来确保离观察者较近的像素不会被较远的像素所覆盖。最为常用的一种方法是深度缓存，深度缓存是一个与图像平面对应的保存每个像素深度的二维数组。每个像素进行绘制的时候都要更新深度缓存中的深度值，每个新像素在绘制之前都要检查深度缓存中的深度值，距离观察者较近的像素就会绘制，而距离较远的都被舍弃。
为了确定像素颜色，需要进行纹理或者浓淡效果计算。纹理图是用于定义三角形显示外观的位图。每个三角形顶点除了位置坐标之外都与纹理以及二维纹理坐标 (u,v) 发生关联。每次渲染三角形中的像素的时候，都必须在纹理中找到对应的纹素，这是根据在屏幕上像素与顶点的距离在与纹理坐标相关联的三角形顶点之间插值完成的。在透视投影中，插值是在根据顶点深度分开的纹理坐标上进行的，这样做就可以避免透视缩减（perspective foreshortening）问题。
在确定像素最终颜色之前，必须根据场景中的所有光源计算像素上的光照。在场景中通常有三种类型的光源。定向光是在场景中按照一个固定方向传输并且强度保持不变的光。在现实生活中，由于太阳距离遥远所以在地球上的观察者看来是平行光线并且其衰减微乎其微，所以太阳光可以看作是定向光。点光源是从空间中明确位置向所有方向发射光线的光源。在远距离的物体上的入射光线会有衰减。最后一种是聚光灯，如同现实生活中的聚光灯一样，它有一个明确的空间位置、方向以及光锥的角度。另外，经常在光照计算完成之后添加一个环境光值以补偿光栅化无法正确计算的全局照明效果。
有许多可以用于光栅化的浓淡算法。所有的浓淡处理算法都必须考虑与光源的距离以及遮蔽物体法向量与光照入射角。最快的算法让三角形中的所有像素使用同样的亮度，但是这种方法无法生成平滑效果的表面。另外也可以单独计算顶点的亮度，然后绘制内部像素的时候对顶点亮度进行插值。速度最慢也最为真实的实现方法是单独计算每点的亮度。常用的浓淡模型有 Gouraud shading 和 Phong shading。


加速技术
为了在任何栅格化引擎中获得最大的性能，只能往渲染工具中发送最少数量的多边形。人们已经开发出了一些加速技术以剔除无法看到的物体。
后向剔除


最简单的剔除多边形的方法就是剔除所有背离观察者的多边形，这就是后向剔除。由于大多数三维物体都是封闭的，所以除非观察者位于物体内部，背离观察者的多边形都会被面向观察者的多边形所遮挡。多边形的方向由它的旋绕方向（winding）或者送到渲染工具的顶点顺序所确定。一旦多边形变换到屏幕空间之后，就可以检查它是否位于相反的方向，一旦如此就丢弃这个多边形。当然，后向剔除不适合于简并的不封闭立体。
空间数据结构


许多先进的技术使用数据结构提出观察物体之外的物体或者被其它物体遮挡的物体，最为常用的数据结构有二元空间分割、八叉树以及 Cell and Portal Culling。
4进一步改进
编辑


尽管基本的栅格化过程已经出现了数十年，许多当今的应用仍然在优化、增加栅格化渲染引擎的应用范围。
纹理映射


纹理是在特定的分辨率生成的，但是由于纹理覆盖的表面与观察者之间可能是任意的距离，所以纹理也可能最后图像上有任意的尺寸。因此，屏幕上的一个像素通常并不直接对应于一个纹素，而是需要使用一些纹理滤波技术来生成任意距离的清晰图像。有许多在图像质量与计算的复杂性进行不同这种考虑的方法可以完成这项工作。
环境映射


环境映射是纹理坐标与观察点相关的纹理映射形式。例如，其中一个常用的应用程序就是用来模拟镜面反射，我们可以将整个房间内部环境映射到房间内的一个金属杯上，当观察者沿着杯子移动的时候，杯子顶点的纹理坐标也随之变化，这样就得到反射效果。
凸凹纹理映射


凸凹纹理映射是改变像素深度而不是颜色的另外一种纹理映射形式。尤其是与最新的阴影工具一起使用的时候，凸凹纹理映射使得表面显现出与光照有关的凸凹不平，从而大幅度地提高真实感。
细节层次


在许多当今的应用中，任何场景中的多边形数目都是非常大的，但是场景中的观察者只能区分近距物体的细节。细节层次算法根据物体与观察者的距离改变几何图形的复杂性。正对着观察者的物体需要进行非常复杂的渲染，而距离远的物体可以动态地简化，甚至可以完全适用二维 sprite 替代。
阴影


传统栅格化过程的光照计算没有考虑物体遮挡的因素。阴影图与阴影体是当今两种生成阴影的普通技术。
5硬件加速
编辑


从1990年代开始，用于消费型桌面电脑的硬件加速开始变得普通。图形程序员早期依赖于手工编写汇编程序以改善程序的运行速度，现在大多数的程序已经使用现有的图形 API 驱动专用的图形处理器。
最新的图形处理器带有可编程的 pixel shader，这可以大幅度地提升程序员的能力。未来的发展趋势就是完全可编程的图形流水线。

八方向栅格化算法代码

/*八方向栅格化算法*/
public List<PixelPosition_T> RasterationMap(Map_T map, Raster_T ras)//栅格化整个选中的地图数据
        {
            List<PixelPosition_T> mappixel = new List<PixelPosition_T>();
            foreach (Geometry_T geo in map.Geofeatures)
            {
                if (geo.GetType() == typeof(Path_T))
                {
                    Path_T pa = (Path_T)geo;
                    foreach (Line_T line in pa.Path)
                    {
                        List<PixelPosition_T> linepixel = RasterationLine(line, ras);
                        mappixel.AddRange(linepixel);
                    }
                }
            }
            return mappixel;
        }


    返回值是 PixelPosition_T 类型的List，输入参数是Map_T类型的map和Raster_T类型的ras；
    new一个PixelPosition_T 类型的List - mappixel；
    对于map中的Geofeatures中的每个Geometry_T类型变量geo循环
        如果geo的类型等于Path_T
            把geo转换为Path_T类型赋给pa
            对pa中的Path中的每个Line_T类型的line循环
                用line和ras调用RasterationLine，返回结果给PixelPosition_T类型的List-linepixel
                在mappixel的范围增加linepixel
    返回mappixel

GIS开发常用算法原理分析

作为计算机科学的一个分支，计算几何主要研究解决几何问题的算法。在现代工程和数学领域，计算几何在图形学、机器人技术、超大规模集成电路设计和统计等诸多领域有着十分重要的应用。在本文中，我们将对计算几何常用的基本算法做一个全面的介绍，希望对您了解并应用计算几何的知识解决问题起到帮助。


矢量的概念


如果一条线段的端点是有次序之分的，我们把这种线段成为有向线段(directed segment)。如果有向线段p1p2的起点p1在坐标原点，我们可以把它称为矢量(vector)p2。




矢量加减法


设二维矢量P = ( x1,y1 )，Q = ( x2 ,y2 )，则矢量加法定义为： P + Q = ( x1 + x2 , y1 + y2 )，同样的，矢量减法定义为： P - Q =( x1 - x2 , y1 - y2 )。显然有性质 P + Q = Q + P , P - Q = - ( Q - P )。




矢量叉积


计算矢量叉积是与直线和线段相关算法的核心部分。设矢量P =（x1,y1），Q = (x2,y2)，则矢量叉积定义为由(0,0)、p1、p2和p1+p2所组成的平行四边形的带符号的面积，即：P× Q = x1*y2 - x2*y1，其结果是一个标量。显然有性质 P× Q = - ( Q× P )和 P× ( - Q ) = - ( P× Q )。一般在不加说明的情况下，本文下述算法中所有的点都看作矢量，两点的加减法就是矢量相加减，而点的乘法则看作矢量叉积。




叉积的一个非常重要性质是可以通过它的符号判断两矢量相互之间的顺逆时针关系：




若 P× Q > 0 ,则P在Q的顺时针方向。若 P × Q < 0 ,则P在Q的逆时针方向。若 P × Q = 0 ,则P与Q共线，但可能同向也可能反向。




折线段的拐向判断


折线段的拐向判断方法可以直接由矢量叉积的性质推出。对于有公共端点的线段p0p1和p1p2，通过计算(p2 - p0)× (p1 - p0)的符号便可以确定折线段的拐向：




若(p2 - p0)× (p1 - p0) > 0,则p0p1在p1点拐向右侧后得到p1p2。




若(p2 - p0)× (p1 - p0) < 0,则p0p1在p1点拐向左侧后得到p1p2。




若(p2 - p0)× (p1 - p0) = 0,则p0、p1、p2三点共线。




判断点是否在线段上


设点为Q，线段为P1P2，判断点Q在该线段上的依据是：( Q- P1 )× ( P2 - P1 ) = 0且 Q在以 P1，P2为对角顶点的矩形内。前者保证Q点在直线P1P2上，后者是保证Q点不在线段P1P2的延长线或反向延长线上，对于这一步骤的判断可以用以下过程实现：




ON-SEGMENT(pi,pj,pk)




if min(xi,xj)<=xk<=max(xi,xj) andmin(yi,yj)<=yk<=max(yi,yj)




then return true;




else return false;




特别要注意的是，由于需要考虑水平线段和垂直线段两种特殊情况，min(xi,xj)<=xk<=max(xi,xj)和min(yi,yj)<=yk<=max(yi,yj)两个条件必须同时满足才能返回真值。




判断两线段是否相交


我们分两步确定两条线段是否相交：




(1)快速排斥试验




设以线段 P1P2为对角线的矩形为R，设以线段 Q1Q2 为对角线的矩形为T，如果R和T不相交，显然两线段不会相交。




(2)跨立试验如果两线段相交，则两线段必然相互跨立对方。若P1P2跨立Q1Q2，则矢量 ( P1 - Q1 )和( P2 - Q1 )位于矢量( Q2 - Q1 )的两侧，即( P1 - Q1 )× ( Q2 - Q1 ) * ( P2 - Q1 )× ( Q2 - Q1 ) < 0。上式可改写成( P1 -Q1 )× ( Q2 - Q1 ) * ( Q2 - Q1 )× ( P2 - Q1 ) > 0。当 ( P1 - Q1)× ( Q2 - Q1 ) = 0时，说明 ( P1 - Q1 )和 ( Q2 - Q1 )共线，但是因为已经通过快速排斥试验，所以 P1一定在线段 Q1Q2上；同理，( Q2 - Q1 )×(P2 - Q1 ) = 0说明 P2一定在线段 Q1Q2上。所以判断P1P2跨立Q1Q2的依据是：( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 )× ( P2 - Q1 ) >= 0。同理判断Q1Q2跨立P1P2的依据是：( Q1 -P1 )× ( P2 - P1 ) * ( P2 - P1 )× ( Q2 - P1 ) >= 0。具体情况如下图所示：








在相同的原理下，对此算法的具体的实现细节可能会与此有所不同，除了这种过程外，大家也可以参考《算法导论》上的实现。




判断线段和直线是否相交


有了上面的基础，这个算法就很容易了。如果线段P1P2和直线Q1Q2相交，则P1P2跨立Q1Q2，即：( P1 - Q1 ) ×( Q2 - Q1 ) * ( Q2 - Q1 )× ( P2 - Q1 ) >= 0。




判断矩形是否包含点


只要判断该点的横坐标和纵坐标是否夹在矩形的左右边和上下边之间。




判断线段、折线、多边形是否在矩形中


因为矩形是个凸集，所以只要判断所有端点是否都在矩形中就可以了。




判断矩形是否在矩形中


只要比较左右边界和上下边界就可以了。




判断圆是否在矩形中


很容易证明，圆在矩形中的充要条件是：圆心在矩形中且圆的半径小于等于圆心到矩形四边的距离的最小值。




判断点是否在多边形中


判断点P是否在多边形中是计算几何中一个非常基本但是十分重要的算法。以点P为端点，向左方作射线L，由于多边形是有界的，所以射线L的左端一定在多边形外，考虑沿着L从无穷远处开始自左向右移动，遇到和多边形的第一个交点的时候，进入到了多边形的内部，遇到第二个交点的时候，离开了多边形，……所以很容易看出当L和多边形的交点数目C是奇数的时候，P在多边形内，是偶数的话P在多边形外。




但是有些特殊情况要加以考虑。如图下图(a)(b)(c)(d)所示。在图(a)中，L和多边形的顶点相交，这时候交点只能计算一个；在图(b)中，L和多边形顶点的交点不应被计算；在图(c)和(d)中，L和多边形的一条边重合，这条边应该被忽略不计。如果L和多边形的一条边重合，这条边应该被忽略不计。








为了统一起见，我们在计算射线L和多边形的交点的时候，1。对于多边形的水平边不作考虑；2。对于多边形的顶点和L相交的情况，如果该顶点是其所属的边上纵坐标较大的顶点，则计数，否则忽略；3。对于P在多边形边上的情形，直接可判断P属于多边行。由此得出算法的伪代码如下： count ← 0;以P为端点，作从右向左的射线L; for多边形的每条边s do if P在边s上 then return true; if s不是水平的 then if s的一个端点在L上 if 该端点是s两端点中纵坐标较大的端点then count← count+1 else if s和L相交 then count← count+1; if count mod 2 = 1 then return true; else return false;




其中做射线L的方法是：设P\'的纵坐标和P相同，横坐标为正无穷大（很大的一个正数），则P和P\'就确定了射线L。




判断点是否在多边形中的这个算法的时间复杂度为O(n)。




另外还有一种算法是用带符号的三角形面积之和与多边形面积进行比较，这种算法由于使用浮点数运算所以会带来一定误差，不推荐大家使用。




判断线段是否在多边形内


线段在多边形内的一个必要条件是线段的两个端点都在多边形内，但由于多边形可能为凹，所以这不能成为判断的充分条件。如果线段和多边形的某条边内交（两线段内交是指两线段相交且交点不在两线段的端点），因为多边形的边的左右两侧分属多边形内外不同部分，所以线段一定会有一部分在多边形外(见图a)。于是我们得到线段在多边形内的第二个必要条件：线段和多边形的所有边都不内交。




线段和多边形交于线段的两端点并不会影响线段是否在多边形内；但是如果多边形的某个顶点和线段相交，还必须判断两相邻交点之间的线段是否包含于多边形内部（反例见图b)。










因此我们可以先求出所有和线段相交的多边形的顶点，然后按照X-Y坐标排序(X坐标小的排在前面，对于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的中点也在多边形内。在实际编程中，没有必要计算所有的交点，首先应判断线段和多边形的边是否内交，倘若线段和多边形的某条边内交则线段一定在多边形外；如果线段和多边形的每一条边都不内交，则线段和多边形的交点一定是线段的端点或者多边形的顶点，只要判断点是否在线段上就可以了。至此我们得出算法如下： if 线端PQ的端点不都在多边形内 then return false;点集pointSet初始化为空; for多边形的每条边s do if线段的某个端点在s上 then将该端点加入pointSet; else if s的某个端点在线段PQ上 then 将该端点加入pointSet; else if s和线段PQ相交 // 这时候已经可以肯定是内交了 then return false;将pointSet中的点按照X-Y坐标排序; for pointSet中每两个相邻点 pointSet[i] , pointSet[ i+1] do if pointSet[i] , pointSet[ i+1]的中点不在多边形中then return false; return true;这个过程中的排序因为交点数目肯定远小于多边形的顶点数目n，所以最多是常数级的复杂度，几乎可以忽略不计。因此算法的时间复杂度也是O(n)。






判断折线是否在多边形内


只要判断折线的每条线段是否都在多边形内即可。设折线有m条线段，多边形有n个顶点，则该算法的时间复杂度为O(m*n)。




判断多边形是否在多边形内




只要判断多边形的每条边是否都在多边形内即可。判断一个有m个顶点的多边形是否在一个有n个顶点的多边形内复杂度为O(m*n)。




判断矩形是否在多边形内


将矩形转化为多边形，然后再判断是否在多边形内。




判断圆是否在多边形内


只要计算圆心到多边形的每条边的最短距离，如果该距离大于等于圆半径则该圆在多边形内。计算圆心到多边形每条边最短距离的算法在后文阐述。




判断点是否在圆内


计算圆心到该点的距离，如果小于等于半径则该点在圆内。




判断线段、折线、矩形、多边形是否在圆内


因为圆是凸集，所以只要判断是否每个顶点都在圆内即可。




判断圆是否在圆内


设两圆为O1,O2，半径分别为r1, r2，要判断O2是否在O1内。先比较r1，r2的大小，如果r1<r2则O2不可能在O1内；否则如果两圆心的距离大于r1 - r2 ，则O2不在O1内；否则O2在O1内。




计算点到线段的最近点


如果该线段平行于X轴（Y轴），则过点point作该线段所在直线的垂线，垂足很容易求得，然后计算出垂足，如果垂足在线段上则返回垂足，否则返回离垂足近的端点；如果该线段不平行于X轴也不平行于Y轴，则斜率存在且不为0。设线段的两端点为pt1和pt2，斜率为：k = ( pt2.y - pt1. y ) / (pt2.x - pt1.x );该直线方程为：y = k*( x - pt1.x) + pt1.y。其垂线的斜率为 - 1 / k，垂线方程为：y = (-1/k) * (x - point.x) + point.y。




联立两直线方程解得：x = ( k^2 * pt1.x + k * (point.y - pt1.y ) + point.x ) / ( k^2 + 1)，y = k * ( x- pt1.x) + pt1.y;然后再判断垂足是否在线段上，如果在线段上则返回垂足；如果不在则计算两端点到垂足的距离，选择距离垂足较近的端点返回。




计算点到折线、矩形、多边形的最近点


只要分别计算点到每条线段的最近点，记录最近距离，取其中最近距离最小的点即可。




计算点到圆的最近距离及交点坐标


如果该点在圆心，因为圆心到圆周任一点的距离相等，返回UNDEFINED。




连接点P和圆心O，如果PO平行于X轴，则根据P在O的左边还是右边计算出最近点的横坐标为centerPoint.x - radius 或centerPoint.x +radius。如果PO平行于Y轴，则根据P在O的上边还是下边计算出最近点的纵坐标为 centerPoint.y -+radius或 centerPoint.y -radius。如果PO不平行于X轴和Y轴，则PO的斜率存在且不为0，这时直线PO斜率为k = （ P.y - O.y）/ ( P.x - O.x )。直线PO的方程为：y = k * ( x - P.x) + P.y。设圆方程为:(x - O.x ) ^2+ ( y - O.y ) ^2 = r ^2，联立两方程组可以解出直线PO和圆的交点，取其中离P点较近的交点即可。




计算两条共线的线段的交点


对于两条共线的线段，它们之间的位置关系有下图所示的几种情况。图(a)中两条线段没有交点；图 (b)和 (d)中两条线段有无穷焦点；图 (c)中两条线段有一个交点。设line1是两条线段中较长的一条，line2是较短的一条，如果line1包含了line2的两个端点，则是图(d)的情况，两线段有无穷交点；如果line1只包含line2的一个端点，那么如果line1的某个端点等于被line1包含的line2的那个端点，则是图(c)的情况，这时两线段只有一个交点，否则就是图(b)的情况，两线段也是有无穷的交点；如果line1不包含line2的任何端点，则是图(a)的情况，这时两线段没有交点。








计算线段或直线与线段的交点


设一条线段为L0 = P1P2，另一条线段或直线为L1= Q1Q2 ，要计算的就是L0和L1的交点。 1．首先判断L0和L1是否相交（方法已在前文讨论过），如果不相交则没有交点，否则说明L0和L1一定有交点，下面就将L0和L1都看作直线来考虑。




2．如果P1和P2横坐标相同，即L0平行于Y轴




a)若L1也平行于Y轴，




i.若P1的纵坐标和Q1的纵坐标相同，说明L0和L1共线，假如L1是直线的话他们有无穷的交点，假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点（该方法在前文已讨论过）； ii. 否则说明L0和L1平行，他们没有交点；




b)若L1不平行于Y轴，则交点横坐标为P1的横坐标，代入到L1的直线方程中可以计算出交点纵坐标；




3．如果P1和P2横坐标不同，但是Q1和Q2横坐标相同，即L1平行于Y轴，则交点横坐标为Q1的横坐标，代入到L0的直线方程中可以计算出交点纵坐标；




4．如果P1和P2纵坐标相同，即L0平行于X轴




a)若L1也平行于X轴，




i.若P1的横坐标和Q1的横坐标相同，说明L0和L1共线，假如L1是直线的话他们有无穷的交点，假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点（该方法在前文已讨论过）； ii. 否则说明L0和L1平行，他们没有交点；




b)若L1不平行于X轴，则交点纵坐标为P1的纵坐标，代入到L1的直线方程中可以计算出交点横坐标；




5．如果P1和P2纵坐标不同，但是Q1和Q2纵坐标相同，即L1平行于X轴，则交点纵坐标为Q1的纵坐标，代入到L0的直线方程中可以计算出交点横坐标；




6．剩下的情况就是L1和L0的斜率均存在且不为0的情况




a)计算出L0的斜率K0，L1的斜率K1；




b)如果K1 = K2




i.如果Q1在L0上，则说明L0和L1共线，假如L1是直线的话有无穷交点，假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点（该方法在前文已讨论过）； ii.如果Q1不在L0上，则说明L0和L1平行，他们没有交点。 c)联立两直线的方程组可以解出交点来这个算法并不复杂，但是要分情况讨论清楚，尤其是当两条线段共线的情况需要单独考虑，所以在前文将求两条共线线段的算法单独写出来。另外，一开始就先利用矢量叉乘判断线段与线段（或直线）是否相交，如果结果是相交，那么在后面就可以将线段全部看作直线来考虑。需要注意的是，我们可以将直线或线段方程改写为 ax+by+c=0的形式，这样一来上述过程的部分步骤可以合并，缩短了代码长度，但是由于先要求出参数，这种算法将花费更多的时间。




求线段或直线与折线、矩形、多边形的交点


分别求与每条边的交点即可。




求线段或直线与圆的交点


设圆心为O，圆半径为r，直线（或线段）L上的两点为P1,P2。




1.如果L是线段且P1，P2都包含在圆O内，则没有交点；否则进行下一步。




2.如果L平行于Y轴，




a)计算圆心到L的距离dis； b) 如果dis > r则L和圆没有交点； c)利用勾股定理，可以求出两交点坐标，但要注意考虑L和圆的相切情况。 3.如果L平行于X轴，做法与L平行于Y轴的情况类似；




4.如果L既不平行X轴也不平行Y轴，可以求出L的斜率K，然后列出L的点斜式方程，和圆方程联立即可求解出L和圆的两个交点；




5.如果L是线段，对于2，3，4中求出的交点还要分别判断是否属于该线段的范围内。




凸包的概念


点集Q的凸包(convexhull)是指一个最小凸多边形，满足Q中的点或者在多边形边上或者在其内。下图中由红色线段表示的多边形就是点集Q={p0,p1,...p12}的凸包。






凸包的求法


现在已经证明了凸包算法的时间复杂度下界是O(n*logn),但是当凸包的顶点数h也被考虑进去的话，Krikpatrick和Seidel的剪枝搜索算法可以达到O(n*logh)，在渐进意义下达到最优。最常用的凸包算法是Graham扫描法和Jarvis步进法。本文只简单介绍一下Graham扫描法，其正确性的证明和Jarvis步进法的过程大家可以参考《算法导论》。




对于一个有三个或以上点的点集Q，Graham扫描法的过程如下：




令p0为Q中Y-X坐标排序下最小的点设<p1,p2,...pm>为对其余点按以p0为中心的极角逆时针排序所得的点集（如果有多个点有相同的极角，除了距p0最远的点外全部移除压p0进栈S压p1进栈S压p2进栈S for i← 3 to m do while由S的栈顶元素的下一个元素、S的栈顶元素以及pi构成的折线段不拐向左侧对S弹栈压pi进栈S return S;




此过程执行后，栈S由底至顶的元素就是Q的凸包顶点按逆时针排列的点序列。需要注意的是，我们对点按极角逆时针排序时，并不需要真正求出极角，只需要求出任意两点的次序就可以了。而这个步骤可以用前述的矢量叉积性质实现。

 GIS算法源码集合

 
其他GIS相关代码下载索引 http://www.mygis.com.cn/codeindex10.htm




1．深度优先实现的路径分析源码 http://www.mygis.com.cn/codes/21-481-481.htm  
2．用遗传算法编写的tsp源码 http://www.mygis.com.cn/codes/21-3961-3961.htm  
3．最短路径算法源码(VB)  http://www.mygis.com.cn/codes/21-579-579.htm  
4．Delphi最短路径 http://www.mygis.com.cn/codes/21-4778-4778.htm  
5．建立公交网络拓扑关系 http://www.mygis.com.cn/codes/21-3962-3962.htm  
6．VC实现公交换乘代码 http://www.mygis.com.cn/codes/21-4671-4671.htm  
7．经纬度BL换算到高斯平面直角坐标XY（高斯投影正算）的源码及算法 http://www.mygis.com.cn/codes/21-1686-1686.htm  
8．Dijkstra算法函数(c++)  http://www.mygis.com.cn/codes/21-11595-11595.htm  
9．C画等值线代码 http://www.mygis.com.cn/codes/21-1708-1708.htm  
10．经典三角插值源代码 http://www.mygis.com.cn/codes/21-2424-2424.htm   
11．delanay三角网算法的VB程序 http://www.mygis.com.cn/codes/21-5483-5483.htm  
12．道格拉斯-普克压缩算法 http://www.mygis.com.cn/codes/21-10424-10424.htm  
13．高斯投影变换_cpp  http://www.mygis.com.cn/codes/21-1721-1721.htm  
14．WGS84墨卡托投影转换C#源码 http://www.mygis.com.cn/codes/21-5588-5588.htm  
15．读写Ini文件 http://www.mygis.com.cn/codes/21-2689-2689.htm  
16．一个VC做的gis系统源代码 http://www.mygis.com.cn/codes/21-479-479.htm  
17．鼠标截取程序源码，可在MapBasic中调用 http://www.mygis.com.cn/codes/21-3593-3593.htm  
18．c#读mif文件/txt文件代码 http://www.mygis.com.cn/codes/21-3810-3810.htm  
19．串口连接GPS实例VB（适合初学者） http://www.mygis.com.cn/codes/21-3860-3860.htm  
20．串口连接GPS实例VC  http://www.mygis.com.cn/codes/21-10945-10945.htm  
21．北京市地理信息公众查询系统mo+c#  http://www.mygis.com.cn/codes/21-8821-8821.htm  
22．PDA上GPS的一个例子 eVB  http://www.mygis.com.cn/codes/21-5639-5639.htm  
23．VB.NET+MapObject查询系统源代码 http://www.mygis.com.cn/codes/21-6299-6299.htm  
24．ArcMap连接SQL Server的三种方法 http://www.mygis.com.cn/codes/21-4677-4677.htm  
25．将点状图层中各点x,y坐标输出到文本文件中 http://www.mygis.com.cn/codes/21-4679-4679.htm  
26．AO中直接加载ArcSDE矢量数据 http://www.mygis.com.cn/codes/21-4682-4682.htm  
27．地理信息系统二次开发实例教程-C#和MapObjects实现 http://www.mygis.com.cn/codes/21-4701-4701.htm  
28．GeoDatabase xml转换为文件xml的源码 http://www.mygis.com.cn/codes/21-5460-5460.htm  
29．MapX打包程序脚本原码 http://www.mygis.com.cn/codes/21-4459-4459.htm  
30．MO的索引图和放大窗口代码_VB  http://www.mygis.com.cn/codes/21-7371-7371.htm  
31．VC版的MO View  http://www.mygis.com.cn/codes/21-6383-6383.htm  
32．C#网络编程 http://www.mygis.com.cn/codes/21-9024-9024.htm  
33．C#播放声音 http://www.mygis.com.cn/codes/21-9025-9025.htm  
34．C#数据库编程一 http://www.mygis.com.cn/codes/21-9022-9022.htm  
35．C#数据库编程二 http://www.mygis.com.cn/codes/21-9023-9023.htm  
36．C#自定义事件 http://www.mygis.com.cn/codes/21-9032-9032.htm  
37．给地图加上比例尺 http://www.mygis.com.cn/codes/21-5461-5461.htm  
38．开发OLE容器的一个代码 http://www.mygis.com.cn/codes/21-2422-2422.htm  
39．AO的一个例子 http://www.mygis.com.cn/codes/21-10943-10943.htm  
40．VC实现GIS矢量图形的绘制 http://www.mygis.com.cn/codes/21-10997-10997.htm   
41．判断点是否在多边形中 http://www.mygis.com.cn/codes/21-11474-11474.htm  
42．VC下面ADO & DAO编程 http://www.mygis.com.cn/codes/21-2519-2519.htm  
43．读写中国地球空间数据交换格式源码 http://www.mygis.com.cn/codes/21-11886-11886.htm  
44．1_1万地形图54与80坐标系换算程序(源码） http://www.mygis.com.cn/codes/21-11889-11889.htm  
45．用VC写的miniGIS  http://www.mygis.com.cn/codes/21-12195-12195.htm  
46．GIS源程序（VC++） http://www.mygis.com.cn/codes/21-12196-12196.htm  
47．获取系统进程讯息 http://www.mygis.com.cn/codes/21-12706-12706.htm  
48．短消息代码 http://www.mygis.com.cn/codes/21-13799-13799.htm