计算几何基础：点在多边形内判定算法详解

算法计算几何多边形GIS

引言：一个看似简单的几何问题

给定一个二维平面上的点和一个多边形，判断该点是否位于多边形内部——这是计算几何中最基础也是最重要的问题之一。

这个问题的应用场景远比想象中广泛。在 GIS（地理信息系统）领域，地理围栏需要判断一个 GPS 坐标是否落在某个行政区划或配送区域内；在游戏引擎中，碰撞检测的核心步骤之一就是判断角色是否进入了某个不规则区域；在图形渲染中，扫描线填充算法需要逐像素判断哪些点落在多边形内部以进行着色。

尽管问题描述简洁，但一旦考虑凹多边形、自交多边形、射线穿过顶点等边界情况，实现一个工程上鲁棒的算法并非易事。本文将系统讲解两种经典算法——射线法（Ray Casting） 与回转数法（Winding Number），深入分析其数学原理、边界处理策略与工程优化方法。

问题定义与数学基础

多边形的分类

在讨论算法之前，有必要明确多边形的分类，因为不同类型的多边形对算法的适用性有直接影响。

简单多边形（Simple Polygon）： 边与边之间除了在顶点处相连外不存在交叉。简单多边形又可细分为凸多边形（任意内角小于 180 度）和凹多边形（存在内角大于 180 度的顶点）。凸多边形的判定可以利用叉积在 O(log n) 时间内完成，但凹多边形无法利用这一性质。

复杂多边形（Complex Polygon）： 边与边之间存在交叉（自交），形成的区域可能存在重叠。自交多边形的"内部"定义取决于填充规则的选择，这正是射线法与回转数法在语义上产生分歧的关键场景。

Jordan 曲线定理

射线法的数学基础是 Jordan 曲线定理：平面上任何一条简单闭曲线将平面划分为恰好两个区域——一个有界的内部区域和一个无界的外部区域。 从内部任意一点出发的射线，在到达无穷远处（外部）的过程中，必须穿过曲线奇数次；从外部出发的射线则穿过偶数次（包括零次）。

这个定理直觉上显而易见，但严格证明并不简单。对于多边形这一特殊的分段线性闭曲线，我们可以在离散层面上精确计算射线与每条边的交点，从而将连续几何问题转化为可计算的算法。

射线法（Ray Casting Algorithm）

核心原理：奇偶规则

射线法的核心思路极为简洁：从测试点 P 向任意方向发射一条射线（通常选择水平向右方向，即正 x 方向），统计该射线与多边形所有边的交点数量。根据 Jordan 曲线定理：

  • 交点数为奇数 -> 点在多边形内部
  • 交点数为偶数 -> 点在多边形外部

直观理解：当测试点在多边形内部时，射线从内部出发，第一次穿越边界即"出去"，随后每一对穿越都对应一次"进入-出去"的循环。整个过程产生 1 + 2k（奇数）个交点。当测试点在外部时，每次穿越都是"进入-出去"成对出现，产生 2k（偶数）个交点。

算法流程

对于一个由 n 个顶点 V_0, V_1, ..., V_{n-1} 定义的多边形和测试点 P(x, y)：

  1. 初始化交点计数器 count = 0
  2. 遍历多边形的每条边 (V_i, V_{i+1})（其中 V_n = V_0 构成闭合）
  3. 对于每条边，判断水平射线是否与该边相交
  4. 若相交，count 加 1
  5. 若 count 为奇数，则点在多边形内

判断射线与边是否相交的具体条件为：边的两个端点必须分布在射线所在水平线的两侧（一个在上方，一个在下方或恰好在线上），且交点的 x 坐标位于测试点的右侧。

边界情况的严格处理

射线法的优雅之处在于其概念的简洁性，但工程实现的难点在于边界情况的处理。以下逐一分析。

射线穿过多边形顶点： 当射线恰好经过多边形的某个顶点时，该顶点同时是两条边的端点，如果简单地对两条边分别判断交点，可能会将同一个交叉点计数两次或零次，导致奇偶性错误。

解决方案是采用半开区间约定（half-open interval）：对于每条边，只有当边的一个端点严格位于射线下方（y < P_y），另一个端点位于射线上方或恰好在射线上（y >= P_y）时，才视为相交。这样，对于穿过顶点的情况：如果该顶点连接的两条边分别位于射线的两侧（一条边从下方来，另一条边向上方去），则恰好产生一个交点；如果两条边位于射线的同一侧（如一条边从下方来到顶点，又向下方去），则交点计数为零或二，对奇偶性无影响。

边与射线重合： 当多边形的某条边恰好与射线共线（即水平且 y 坐标等于 P_y）时，该边上的所有点都是"交点"。处理策略同样依赖半开区间约定：水平边的两个端点 y 坐标相同，均不满足"一上一下"的条件，因此整条水平边不产生任何交点。这在数学上等价于将水平边"忽略"。

自交多边形： 射线法天然适用于自交多边形，此时它实际上遵循的是奇偶填充规则（Even-Odd Rule）。自交区域被穿过偶数次的边界包围，在奇偶规则下被视为"外部"——即重叠区域通过 XOR 操作被剔除。这在某些场景下是期望的行为（如 SVG 的 fill-rule: evenodd），但在另一些场景下可能不符合预期。

算法实现

以下是射线法的 C 语言实现，采用水平向右射线与半开区间约定：

bool pointInPolygon(double x, double y,
                    double *polyX, double *polyY, int polySides) {
    int i, j = polySides - 1;
    bool oddNodes = false;

    for (i = 0; i < polySides; i++) {
        if ((polyY[i] < y && polyY[j] >= y)
        ||  (polyY[j] < y && polyY[i] >= y)) {
            if (polyX[i] + (y - polyY[i]) / (polyY[j] - polyY[i])
                * (polyX[j] - polyX[i]) < x) {
                oddNodes = !oddNodes;
            }
        }
        j = i;
    }
    return oddNodes;
}

逐行解析：

  • 外层 for 循环遍历多边形每条边 (V_i, V_j)，其中 j 始终指向前一个顶点，形成循环遍历
  • 第一个 if 判断：边的两个端点是否分布在水平射线的两侧（采用半开区间 < 和 >=）
  • 第二个 if 判断：计算射线与边的交点 x 坐标，判断交点是否在测试点的左侧（即测试点在交点右侧），若是则翻转 oddNodes 标志

交点 x 坐标的计算基于线性插值：给定边的两个端点 (x_i, y_i) 和 (x_j, y_j)，射线 y = P_y 与该边的交点 x 坐标为：

x_intersect = x_i + (P_y - y_i) / (y_j - y_i) * (x_j - x_i)

复杂度分析

  • 时间复杂度：O(n)，其中 n 为多边形的边数。算法仅需遍历所有边一次，每条边的计算为常数时间
  • 空间复杂度：O(1)，仅使用常数个辅助变量

边界点的归属

需要特别说明的是，当测试点恰好位于多边形边上时，射线法的返回值是不确定的——可能返回"内"也可能返回"外"，取决于浮点运算的精度和具体的边方位。在工程实践中，由于边在数学上是无面积的一维对象，这种不确定性在绝大多数场景下可以忽略。如果应用需要对边界点做确定性判断，需要额外增加点到边的距离判断逻辑。

回转数法（Winding Number Algorithm）

原理：绕测试点的回转数

回转数法从另一个视角定义"内部"：计算多边形边界绕测试点旋转的总角度。具体而言，回转数（Winding Number） 定义为：沿多边形边界按顶点顺序遍历一圈时，从测试点观察，边界总共绕测试点逆时针旋转了多少圈。

  • 回转数不为零 -> 点在多边形内部
  • 回转数为零 -> 点在多边形外部

对于简单多边形，回转数只可能是 +1、-1 或 0，此时回转数法与射线法的结果完全一致。两者的差异体现在自交多边形上。

与射线法的语义差异

考虑一个"8"字形自交多边形，其两个环重叠区域的回转数为 +2 或 -2。在回转数法下，重叠区域属于"内部"（回转数非零）；在射线法（奇偶规则）下，重叠区域属于"外部"（射线穿过偶数条边）。

这两种语义分别对应图形学中的两种填充规则：

  • 奇偶规则（Even-Odd Rule）： 射线法的行为，自交区域被 XOR 剔除
  • 非零规则（Non-Zero Rule）： 回转数法的行为，只要回转数非零即为内部

SVG 和 CSS 中的 fill-rule 属性正是控制这一行为的：evenodd 对应射线法语义，nonzero 对应回转数法语义。

实现思路

直观的实现是逐边累加角度变化量，但涉及反三角函数计算，效率较低。高效的实现方法同样基于边与水平线的交叉计数，但区分方向：

function windingNumber(P, polygon):
    wn = 0
    for each edge (V_i, V_{i+1}):
        if V_i.y <= P.y:
            if V_{i+1}.y > P.y:                // 向上穿越
                if isLeft(V_i, V_{i+1}, P) > 0: // P 在边的左侧
                    wn += 1
        else:
            if V_{i+1}.y <= P.y:               // 向下穿越
                if isLeft(V_i, V_{i+1}, P) < 0: // P 在边的右侧
                    wn -= 1
    return wn != 0

其中 isLeft 函数通过叉积判断点 P 相对于有向边 (V_i, V_{i+1}) 的方位：

isLeft(V_i, V_{i+1}, P) =
    (V_{i+1}.x - V_i.x) * (P.y - V_i.y) -
    (P.x - V_i.x) * (V_{i+1}.y - V_i.y)

该实现的时间复杂度同样为 O(n)，且不涉及三角函数运算，全部计算基于加减乘运算，在性能上与射线法相当。

射线法与回转数法的对比

维度 射线法 回转数法
填充规则 奇偶规则（Even-Odd） 非零规则（Non-Zero）
简单多边形 结果一致 结果一致
自交多边形 重叠区域视为外部 重叠区域视为内部
时间复杂度 O(n) O(n)
实现难度 较低 略高（需处理方向）
数值稳定性 涉及除法 仅涉及乘法（叉积）

在大多数工程场景中，如果处理的是简单多边形，两种算法可以互换使用。对于自交多边形，需要根据业务语义选择合适的填充规则。

算法优化

包围盒预检测（Bounding Box）

最直接的优化是在执行射线法之前，先判断测试点是否落在多边形的轴对齐包围盒（AABB）内。包围盒通过预计算多边形所有顶点的 x、y 坐标极值得到：

if (x < minX || x > maxX || y < minY || y > maxY)
    return false;  // 点在包围盒之外，必定不在多边形内

包围盒预计算为 O(n)，一次计算后可复用。后续每次判定的预检测为 O(1)，能快速排除大量明显在多边形外部的点。

边的早期跳过优化

在射线法的实现中，Nathan Mercer 提出了一种优化：在判断边的两端点是否跨越水平线之后，增加一个额外条件 polyX[i] <= x || polyX[j] <= x，用于跳过那些完全位于测试点右侧的边——这些边不可能与向右发射的射线产生在测试点左侧的交点。

if ((polyY[i] < y && polyY[j] >= y)
||  (polyY[j] < y && polyY[i] >= y)) {
    if (polyX[i] <= x || polyX[j] <= x) {  // 新增条件
        if (polyX[i] + (y - polyY[i]) / (polyY[j] - polyY[i])
            * (polyX[j] - polyX[i]) < x) {
            oddNodes = !oddNodes;
        }
    }
}

这一优化避免了对完全位于右侧的边执行浮点除法和乘法运算，在多边形边数较多且测试点靠近多边形中心时效果显著。

空间索引加速

当需要对同一个多边形进行大量点的判定（如 GIS 服务中的批量坐标归属查询），O(n) 的逐边遍历会成为性能瓶颈。此时可以引入空间索引进行预处理：

网格法（Grid）： 将多边形的包围盒划分为 m x m 的均匀网格。预处理阶段，计算每个网格单元与多边形的关系：完全在内部、完全在外部、或与边界相交。查询时，先定位测试点所在的网格单元，若为完全内部或完全外部则直接返回，仅对边界单元执行射线法。预处理复杂度为 O(n * m)，查询复杂度降至 O(1)（最优情况）到 O(n/m)（最坏情况）。

梯形分解（Trapezoidal Decomposition）： 将多边形分解为若干梯形区域，构建搜索结构。预处理复杂度 O(n log n)，单次查询复杂度 O(log n)。适用于对查询延迟有严格要求的场景。

R-Tree 索引： 将多边形的边组织为 R-Tree。查询时，先通过 R-Tree 筛选出与水平射线可能相交的边，再对这些边执行精确判断。当多边形边数达到数千甚至数万级别时，R-Tree 能将实际需要检查的边数降低一到两个数量级。

批量判定的预处理策略

在实际工程中，常见的场景是：多边形固定不变，需要对海量动态点进行判定。典型如地理围栏服务——行政区划多边形固定，而 GPS 坐标持续涌入。

此时最高效的策略是基于扫描线的预处理：将包围盒在 y 方向上分为若干水平带（slab），每个水平带内记录穿过该带的边集合。查询时先通过 y 坐标定位到水平带，再只对该带内的边执行射线法。这种方法将每次查询的平均边检查数量从 n 降低到 n/k（k 为水平带数量），是 GIS 工程实践中常用的优化手段。

工程应用

GIS 地理围栏

地理围栏（Geofencing）是点在多边形内判定最典型的工程应用。在物流配送中，需要判断用户地址是否在配送区域内；在出行服务中，需要判断车辆是否驶入或驶出电子围栏；在气象服务中，需要判断某个坐标属于哪个气象分区。

GIS 场景的特殊性在于坐标系。GPS 坐标是经纬度（球面坐标），而射线法假设平面直角坐标系。对于小范围区域（如城市级别），可以直接将经纬度视为平面坐标进行计算，误差可控。对于大范围区域（如跨越数千公里的国界），需要先进行墨卡托投影或其他地图投影，将球面坐标转换为平面坐标后再执行算法。

此外，GIS 中的多边形边数通常较大——一个省级行政区划的边界可能由数万个顶点组成。此时空间索引（如 R-Tree）和网格预处理几乎是必需的优化手段。PostGIS 中的 ST_Contains 函数内部正是采用了 R-Tree 索引加速的射线法实现。

地图服务：行政区划归属查询

大型地图服务需要支持"根据经纬度查询所属行政区划"的高频请求。其本质是对数百个行政区划多边形的批量点包含测试。工程实现通常分为两级：

  1. 粗筛： 使用 R-Tree 对所有行政区划的包围盒建立索引，根据查询点快速定位到候选区划（通常只有 1-3 个）
  2. 精确判定： 对候选区划执行射线法或回转数法

这种两级架构将原本 O(N * n) 的问题（N 个区划，每个区划 n 条边）降低到接近 O(log N + n') 的复杂度，其中 n' 是候选区划的平均边数。

游戏引擎碰撞检测

在游戏开发中，碰撞检测通常采用多阶段策略：

  1. 宽阶段（Broad Phase）： 使用包围盒（AABB 或 OBB）快速排除不可能碰撞的物体对
  2. 窄阶段（Narrow Phase）： 对宽阶段筛选出的候选对执行精确碰撞检测

点在多边形内判定在窄阶段发挥作用，尤其对于不规则形状的碰撞区域（如地形边界、不规则触发区域）。对于凸多边形，通常使用 GJK（Gilbert-Johnson-Keerthi）算法，但对于凹多边形区域，射线法仍然是最直接可靠的方案。

数据可视化：区域着色与热力图

在前端数据可视化中，Choropleth Map（分级统色地图）需要将数据点分配到对应的地理区域。例如，将每个用户的坐标归属到对应的省份，再根据各省份的用户密度进行着色。这本质上就是大规模的点在多边形内批量判定。

D3.js 和 ECharts 等可视化库内部均实现了射线法。在 Canvas 渲染场景中，CanvasRenderingContext2D.isPointInPath() 方法底层也是基于类似的算法实现。

总结

点在多边形内判定是计算几何中的基石问题。射线法以其概念简洁、实现高效的特点成为工程实践中的首选方案，其 O(n) 的时间复杂度和对凹多边形的天然支持使其适用于绝大多数场景。回转数法在处理自交多边形时提供了不同的语义选择（非零填充规则 vs 奇偶填充规则），是理解图形填充行为的关键。

在工程应用中，算法本身的效率往往不是瓶颈——真正的挑战在于如何通过包围盒预检测、空间索引、网格预处理等手段，将海量查询的平均复杂度从 O(n) 降低到接近 O(1)。理解算法原理是基础，而将其融入具体的工程架构中进行系统性优化，才是将理论转化为生产力的关键。

加载导航中...

评论