电子海图中点在多边形内的判定

1          算法背景
       在电子海图系统中，经常需要用到判断一个点(可能是一个船位，或者其它点状物标)是否在某个多边形区域内(矩形或圆也可看作特殊的多边形)，比如某个警戒区，作业区，禁航区，台风影响区以及其它用户关注的区域。一般点在多边形的边线上时也作为在多边形内的一种特例。

2          假设
2.1          经度范围
       在各种电子海图系统中，经度范围的表示方式有多种，这里我们约定经度的取值范围为-180~180。其中国际变更线既是180也是-180。

2.2      两点间走小圆弧
       由于地球是圆形的，如果只给出两个点，那么它们之间可以形成两条连线(可以有大圆弧和小圆弧两条连线)。这里我们假设我们所需判断的多边形的任何两条连线都是走小圆弧。这个假设一般情况下也是成立的，用户不会需要绘制一个大小超过180度范围的区域。

3          算法特点
       电子海图的坐标系和一般的数学坐标系不同。数学坐标系在一个方向上始终保持递增或者递减，而电子海图中，当穿过经度180度的国际变更线时，经度值会变回-180度。因此，对于区域跨国际变更线的情况要进行特殊处理。

4          算法比较
       对于点在多边形内的判断，在数学上有多种策略。常见的有面积法，角度法，射线法等；在开发工具中，还有CRgn对象的PtInRegion方法。由于面积法和角度法需要大量的计算，因此不予考虑。CRgn的PtInRegion虽然代码简单，但CRgn本身就是很消耗资源的一个对象，同时也受计算机开发工具的限制，也不在考虑之列。此处选择射线法作为解决方案的基础。

5          线段法
       线段法来源与数学方法中的射线法。

       射线法的定义是从要判定的点向任意方向形成一条射线，计算该射线与多边形所有边的相交情况，如果相交的边的数量是奇数，则该点在多边形内；反之则在多边形外。

       结合电子海图的特点，纬度方向呈单向递增/递减，并且纬度值不可能大于90度。因此产生一个新的方法—线段法(笔者的命名)，即不再是发出一条射线，而是从需要判定的点保持经度不变，向正北延伸直到90度纬线，形成一条线段，判定多边形的各条边与该线段的相交情况，以此判断出点是否在多边形内。

       线段法中的线段命名为判定线段。

 

6          判定过程
6.1      判定多边形区域是否跨越国际变更线
       由于跨越国际变更线的多边形顶点坐标不是在递增坐标系下，因此必须进行判定。如果确实跨越了国际变更线，则需要将所有经度值小于0的点加上360度，即可转换为递增序列。同时，如果需要判定的点的经度值小于0，也需要加上180度。

       是否跨越国际变更线的判定方法如下：

       逐条判断多边形的每条边，如果相邻两点的经度值乘积小于0，并且两者的绝对值之和大于180，则确认该多边形跨越国际变更线。

       因为相邻两点乘积小于0时，有两种可能，或者跨越0度线，或者跨越180度线。如果跨越的是0度线，那么两点的经度值的绝对值之和必然小于180。这里用到了我们两点间走小圆弧的假设，只有在此条件下判定才能成立。

 

6.2      顶点过滤
       如果需要判定的点就是多边形的某个顶点，则直接认为点在多边形内。判定点与顶点重合会影响后续的判断逻辑。这里需要注意两个点坐标判定相等的一个精度问题。

6.3      经纬度过滤
       如果多边形任意一条边不满足以下条件，则可不予处理：

1．边的两点的纬度值必须有一个在需要判定的点的纬度到90度之间(开区间)；

2．需要判定的点的经度值必须在两边的经度值之间(开区间)；

       条件中都为开区间，以此过滤判定线段正好经过多边形顶点的情况。这种情况下，经过一个点，等于与两条边相交，对最终判定没有影响。

6.4      交点计算
       经过以上几个步骤后，还剩下以下图示的情况，需要过滤边ab。

 

   

       图中，边ab和判定线段pp0满足6.3条件。因此需要计算p0点在边ab上的投影点才能确定两线段是否相交。由于pp0与ab的交点的经度肯定就是p0点的经度，因此只需要计算交点的纬度值。如果纬度值大于p0，则交点在pp0上，确定该边与pp0相交；如果纬度值等于p0，则表明p0就在ab上，直接认定p0点在区域内。

       所有边均判定完成后，计算与pp0相交的边的数量，如果是奇数，则认定p0在多边形内，否则不在。

       p0点在ab边上的投影点计算：

       假设a(x1,y1),b(x2,y2),p0在ab线上的投影p1(x0,y0)。那么现在未知值就是y0。根据简单的数学公式，可得以下算式：（在前几步骤中，已经排除了x1和x2相等的可能，并且确定p1点在线段ab上）

       (y1-y0)(x1-x0) = (y1-y2)/(x1-x2)

转换为计算机公式就是：y0 = y1 – (y1-y2)*(x1-x0)/(x1-x2)

       根据已知条件即可计算出y0。

 

7          总结
       本方法采用了纯数学的算法。与具体开发语言和工具没有关联。可最大限度的进行共享。而且选用了计算量最小的数学公式，保证了判定的效率。当然，由于笔者能力有限，或许还有更为简单的算法，或许本文中还有部分错误和遗漏之处，欢迎大家批评指正，互相学习。

 

附源代码

view plaincopy to clipboardprint?
//判断指定的坐标是否在水域范围内  
BOOL CStressSeaArea::IsPosIn(DOUBLE_LONG_LAT dllPos)  
{  
    int nSize = m_StressSeaAreaInfo.arLatLong.GetSize();  
    if(nSize <= 2)  
        return FALSE;  
    BOOL bCross180 = IsCrossLong180();  
    //如果穿过180度经线，负经度值需要转换为正经度值  
    if(bCross180 && dllPos.dLongitude < 0)  
    {  
        dllPos.dLongitude += 360;  
    }  
    int nCrossNum = 0;  
    for(int i=0; i<nSize; i++)  
    {  
        DOUBLE_LONG_LAT dllPos1 = m_StressSeaAreaInfo.arLatLong.GetAt(i);  
        DOUBLE_LONG_LAT dllPos2;  
        if(i!=nSize-1)//最后一点与第一点形成最后一条线段  
            dllPos2 = m_StressSeaAreaInfo.arLatLong.GetAt(i+1);  
        else 
            dllPos2 = m_StressSeaAreaInfo.arLatLong.GetAt(0);  
        if(bCross180)  
        {  
            if(dllPos1.dLongitude < 0)  
                dllPos1.dLongitude += 360;  
            if(dllPos2.dLongitude < 0)  
                dllPos2.dLongitude += 360;  
        }  
        //如果顶点就是dllPos，则认为点在水域内  
        if((fabs(dllPos1.dLatitude - dllPos.dLatitude) < 0.000001) &&   
            (fabs(dllPos1.dLongitude - dllPos.dLongitude) < 0.000001))  
            return TRUE;  
        if((fabs(dllPos2.dLatitude - dllPos.dLatitude) < 0.000001) &&   
            (fabs(dllPos2.dLongitude - dllPos.dLongitude) < 0.000001))  
            return TRUE;  
        //如果两点中有一个点的纬度在90~dllPos.dLatitude间，则满足纬度条件  
        if(((dllPos1.dLatitude < 90) && (dllPos1.dLatitude > dllPos.dLatitude)) ||  
            ((dllPos2.dLatitude < 90) && (dllPos2.dLatitude > dllPos.dLatitude)))  
        {  
            //继续判断经度是否满足条件。如果dllPos.dLongitude在两点的经度之间，则满足条件  
            if(((dllPos.dLongitude > dllPos1.dLongitude) && (dllPos.dLongitude < dllPos2.dLongitude)) ||  
                ((dllPos.dLongitude > dllPos2.dLongitude) && (dllPos.dLongitude < dllPos1.dLongitude)))  
            {  
                //计算两条线段的交点的纬度值  
                float fLat = dllPos1.dLatitude - (dllPos1.dLatitude-dllPos2.dLatitude)*(dllPos1.dLongitude-dllPos.dLongitude)/(dllPos1.dLongitude-dllPos2.dLongitude);  
                //如果纬度值就是dllPos.dLatitude，表明点在多边形边线上，认为在多边形内  
                if(fabs(fLat - dllPos.dLatitude) < 0.00001)  
                    return TRUE;  
                //如果纬度值大于dllPos.dLatitude，表明两条线段是相交线段，计数加1  
                if(fLat > dllPos.dLatitude)  
                    nCrossNum++;  
            }  
        }  
    }  
    if(nCrossNum %2 == 1)  
        return TRUE;  
    return FALSE;  
}  
 
 
BOOL CStressSeaArea::IsCrossLong180()  
{  
    int nSize = m_StressSeaAreaInfo.arLatLong.GetSize();  
    if(nSize <= 2)  
        return FALSE;  
    //判断水域是否跨180经线  
    double dLong1 = 0;  
    int i=0;  
    for(i=0; i<=nSize; i++)  
    {  
        DOUBLE_LONG_LAT pos;  
        if(i==nSize)  
            pos = m_StressSeaAreaInfo.arLatLong.GetAt(0);  
        else 
            pos = m_StressSeaAreaInfo.arLatLong.GetAt(i);  
        if(i!=0)  
        {  
            if(pos.dLongitude * dLong1 < 0)//如果相邻两点的经度值符号相反，则表明肯定跨0度或者180度经线  
            {  
                //如果绝对值之和大于180，那么就是跨180度经线(这里有个假设：两点间经度差小于180)  
                if(fabs(dLong1) + fabs(pos.dLongitude) >= 180)  
                {  
                    return TRUE;  
                }  
            }  
        }  
        dLong1 = pos.dLongitude;  
    }  
    return FALSE;  
}