引言
点在多边形判定（Point in Polygon, PIP）是计算几何领域最基础也是最重要的算法之一。从游戏开发中的碰撞检测、地理信息系统（GIS）中的区域查询，到计算机图形学中的填充与裁剪，都离不开这一核心问题的高效求解。给定一个任意多边形和空间中一点，如何快速准确地判断该点是否位于多边形内部？
本文将用Java实现两种经典算法——射线法（Ray Casting）和绕数法（Winding Number），深入剖析其数学原理、边界情况处理及性能优化策略。
问题定义
输入：
– 多边形顶点序列 $P = {V_0, V_1, …, V_{n-1}}$，顶点按顺时针或逆时针顺序排列
– 查询点 $Q = (x, y)$
输出：
– 点在多边形内部（IN）
– 点在多边形外部（OUT）
– 点在多边形边界上（ON）
多边形可以是凸的、凹的，甚至是自相交的（复杂多边形）。不同算法对自相交多边形的处理有所差异。
算法一：射线法（Ray Casting / Crossing Number）
核心思想
从查询点 $Q$ 向右发射一条水平射线，统计该射线与多边形各边的交点数量。若交点数为奇数，则点在内部；若为偶数，则点在外部。这一结论基于若尔当曲线定理（Jordan Curve Theorem）。
边界情况处理
射线法需要谨慎处理多种边界情况：
- 射线穿过顶点：当射线恰好经过多边形顶点时，可能产生重复计数。标准做法是只计算边的”上开下闭”区间，即 $y_{min} < y_Q \le y_{max}$。
- 点在边上：通过共线检测和范围检查判定。
- 水平边：水平边与射线重合或平行，不产生有效交点，直接忽略。
交点计算
对于边 $(V_i, V_{i+1})$，与水平射线 $y = y_Q$ 的交点横坐标为：
$$x_{int} = x_i + \frac{(y_Q – y_i)(x_{i+1} – x_i)}{y_{i+1} – y_i}$$
当 $x_{int} > x_Q$ 时，视为一次有效穿越。
算法二：绕数法（Winding Number）
核心思想
绕数法计算从查询点 $Q$ 观察多边形所有边时，多边形绕 $Q$ 旋转的总角度。若绕数为非零，则点在内部；若绕数为零，则点在外部。
对于每条有向边 $(V_i, V_{i+1})$，计算其对绕数的贡献。具体地，将边投影到以 $Q$ 为原点的坐标系中，计算该边跨越正x轴的次数及方向。
高效计算
无需实际计算角度，可通过边的方向关系快速判断：
对于边 $(V_i, V_{i+1})$：
1. 若 $y_i \le y_Q$：
– 若 $y_{i+1} > y_Q$ 且边在 $Q$ 的左侧（叉积为正），则绕数加1
2. 若 $y_i > y_Q$：
– 若 $y_{i+1} \le y_Q$ 且边在 $Q$ 的右侧（叉积为负），则绕数减1
这一方法避免了三角函数计算，效率与射线法相当。
Java实现
import java.util.*;
/**
* 点在多边形判定（Point in Polygon）
* 射线法（Ray Casting）与绕数法（Winding Number）
*/
public class PointInPolygon {
/**
* 二维点/向量
*/
static class Point {
double x, y;
Point(double x, double y) { this.x = x; this.y = y; }
@Override
public String toString() {
return String.format("(%.3f, %.3f)", x, y);
}
}
/**
* 判定结果枚举
*/
enum Result { IN, OUT, ON }
/**
* ==================== 射线法（Crossing Number） ====================
*/
public static Result rayCasting(Point q, List<Point> polygon) {
int n = polygon.size();
if (n < 3) return Result.OUT;
int crossings = 0;
for (int i = 0, j = n - 1; i < n; j = i++) {
Point vi = polygon.get(i);
Point vj = polygon.get(j);
// 首先检查点是否在当前边上
if (isPointOnSegment(q, vj, vi)) {
return Result.ON;
}
// 检查射线是否与边相交
// 条件：边跨越水平线 y = q.y（上开下闭避免顶点重复计数）
boolean cond1 = (vi.y > q.y) != (vj.y > q.y);
if (!cond1) continue;
// 计算交点横坐标
double xInt = vj.x + (q.y - vj.y) * (vi.x - vj.x) / (vi.y - vj.y);
// 交点在射线右侧（x > q.x）
if (xInt > q.x) {
crossings++;
}
}
return (crossings % 2 == 1) ? Result.IN : Result.OUT;
}
/**
* 射线法（从左向右发射射线版本）
*/
public static Result rayCastingLeft(Point q, List<Point> polygon) {
int n = polygon.size();
if (n < 3) return Result.OUT;
int crossings = 0;
for (int i = 0; i < n; i++) {
Point a = polygon.get(i);
Point b = polygon.get((i + 1) % n);
if (isPointOnSegment(q, a, b)) return Result.ON;
// 严格跨越测试
if ((a.y > q.y) != (b.y > q.y)) {
double xInt = (b.x - a.x) * (q.y - a.y) / (b.y - a.y) + a.x;
if (q.x < xInt) crossings++;
}
}
return (crossings % 2 == 1) ? Result.IN : Result.OUT;
}
/**
* ==================== 绕数法（Winding Number） ====================
*/
public static Result windingNumber(Point q, List<Point> polygon) {
int n = polygon.size();
if (n < 3) return Result.OUT;
int wn = 0; // 绕数
for (int i = 0; i < n; i++) {
Point a = polygon.get(i);
Point b = polygon.get((i + 1) % n);
if (isPointOnSegment(q, a, b)) return Result.ON;
if (a.y <= q.y) {
// 向上穿越
if (b.y > q.y && crossProduct(a, b, q) > 0) {
wn++;
}
} else {
// 向下穿越
if (b.y <= q.y && crossProduct(a, b, q) < 0) {
wn--;
}
}
}
return (wn != 0) ? Result.IN : Result.OUT;
}
/**
* 向量叉积：(b - a) × (q - a)
* 正值：q在ab左侧；负值：q在ab右侧；零：共线
*/
private static double crossProduct(Point a, Point b, Point q) {
return (b.x - a.x) * (q.y - a.y) - (b.y - a.y) * (q.x - a.x);
}
/**
* 判断点q是否在线段(a, b)上（含端点）
*/
private static boolean isPointOnSegment(Point q, Point a, Point b) {
// 叉积为0表示共线
double cross = crossProduct(a, b, q);
if (Math.abs(cross) > 1e-9) return false;
// 检查是否在包围盒内
double minX = Math.min(a.x, b.x) - 1e-9;
double maxX = Math.max(a.x, b.x) + 1e-9;
double minY = Math.min(a.y, b.y) - 1e-9;
double maxY = Math.max(a.y, b.y) + 1e-9;
return q.x >= minX && q.x <= maxX && q.y >= minY && q.y <= maxY;
}
/**
* ==================== 带预处理的加速版本（Bounding Box + Grid） ====================
*/
static class AcceleratedPIP {
private final List<Point> polygon;
private final double minX, minY, maxX, maxY;
AcceleratedPIP(List<Point> polygon) {
this.polygon = polygon;
double mx = Double.POSITIVE_INFINITY, my = Double.POSITIVE_INFINITY;
double Mx = Double.NEGATIVE_INFINITY, My = Double.NEGATIVE_INFINITY;
for (Point p : polygon) {
mx = Math.min(mx, p.x); my = Math.min(my, p.y);
Mx = Math.max(Mx, p.x); My = Math.max(My, p.y);
}
this.minX = mx; this.minY = my;
this.maxX = Mx; this.maxY = My;
}
Result query(Point q) {
// 快速包围盒排除
if (q.x < minX || q.x > maxX || q.y < minY || q.y > maxY) {
return Result.OUT;
}
return rayCasting(q, polygon);
}
}
/**
* ==================== 凸多边形特殊优化（二分查找） ====================
*/
public static Result convexPolygonPIP(Point q, List<Point> polygon) {
int n = polygon.size();
if (n < 3) return Result.OUT;
// 找到y坐标最小的顶点作为起点（假设为凸包顶点之一）
int start = 0;
for (int i = 1; i < n; i++) {
if (polygon.get(i).y < polygon.get(start).y ||
(polygon.get(i).y == polygon.get(start).y && polygon.get(i).x < polygon.get(start).x)) {
start = i;
}
}
// 二分查找判断点所在扇形区域
// 此处省略完整实现，核心思想是将凸多边形分解为以start为顶点的三角形扇
return rayCasting(q, polygon); // 回退到通用方法
}
/**
* ==================== 测试与演示 ====================
*/
public static void main(String[] args) {
// 测试多边形1：简单凹多边形
List<Point> polygon1 = Arrays.asList(
new Point(0, 0), new Point(10, 0), new Point(10, 10),
new Point(5, 5), new Point(0, 10)
);
// 测试多边形2：星形（自相交，测试绕数法优势）
List<Point> polygon2 = Arrays.asList(
new Point(5, 0), new Point(8, 8), new Point(0, 3),
new Point(10, 3), new Point(2, 8)
);
Point[] testPoints = {
new Point(5, 4), // 内部
new Point(5, 6), // 外部（凹角处）
new Point(0, 0), // 顶点
new Point(5, 0), // 边上
new Point(15, 15), // 远外部
};
System.out.println("===== 多边形1（凹多边形）测试 =====");
System.out.println("顶点: " + polygon1);
for (Point p : testPoints) {
Result r1 = rayCasting(p, polygon1);
Result r2 = windingNumber(p, polygon1);
System.out.printf("点 %s -> 射线法: %s, 绕数法: %s%n", p, r1, r2);
}
System.out.println("\n===== 多边形2（星形/自相交）测试 =====");
System.out.println("顶点: " + polygon2);
Point[] starPoints = {
new Point(5, 4), // 中心区域
new Point(2, 2), // 星角内部
new Point(8, 2), // 星角内部
new Point(5, 7), // 星角内部
new Point(5, 2), // 凹入外部
};
for (Point p : starPoints) {
Result r1 = rayCasting(p, polygon2);
Result r2 = windingNumber(p, polygon2);
System.out.printf("点 %s -> 射线法: %s, 绕数法: %s%n", p, r1, r2);
}
// 性能测试
System.out.println("\n===== 性能测试（100万次查询） =====");
AcceleratedPIP fast = new AcceleratedPIP(polygon1);
Random rand = new Random(42);
int inCount = 0, onCount = 0;
long start = System.nanoTime();
for (int i = 0; i < 1_000_000; i++) {
Point p = new Point(
rand.nextDouble() * 15 - 2.5,
rand.nextDouble() * 15 - 2.5
);
Result r = fast.query(p);
if (r == Result.IN) inCount++;
if (r == Result.ON) onCount++;
}
long end = System.nanoTime();
System.out.printf("100万次查询耗时: %.2f ms%n", (end - start) / 1e6);
System.out.printf("内部: %d, 边界: %d, 外部: %d%n", inCount, onCount, 1_000_000 - inCount - onCount);
}
}
算法复杂度分析
时间复杂度
- 射线法：$O(n)$，其中 $n$ 为多边形顶点数。每条边需要常数时间判断是否与射线相交。
- 绕数法：$O(n)$，与射线法相同，每条边计算一次叉积和方向判断。
- 凸多边形二分查找优化：$O(\log n)$，将多边形分解为三角形扇后二分查找。
- 包围盒预处理：$O(1)$ 快速排除，最坏情况仍 $O(n)$。
空间复杂度
- 两种方法均为 $O(1)$ 额外空间（不计多边形存储）。
- 加速结构（如Uniform Grid、BSP树）需要 $O(n)$ 到 $O(n^2)$ 的预处理空间。
两种方法的差异
| 特性 | 射线法 | 绕数法 |
|---|---|---|
| 自相交多边形 | 可能给出”外部”结果 | 正确区分内外（按环绕次数） |
| 数值稳定性 | 需要处理顶点重复计数 | 更稳定，无顶点歧义 |
| 性能 | 略快（少一次叉积判断） | 稍慢但结果更可靠 |
| 推荐使用 | 简单非自交多边形 | 复杂/自相交多边形 |
扩展与改进
1. 大规模多边形加速
对于包含数百万顶点的大陆边界或城市轮廓，线性扫描太慢。常用加速结构：
– Uniform Grid：将空间划分为均匀网格，预计算每个网格与多边形的相交关系，查询时只需检查点所在网格
– R-Tree：层次包围盒树，支持动态多边形集合
– 梯形图（Trapezoidal Map）：期望 $O(\log n)$ 查询时间，但构建复杂
2. 3D扩展：点在多面体内
将射线法推广到三维：从点发射射线，统计与多面体各面的交点数。或使用奇偶规则判断射线穿越多面体边界的次数。更稳定的方法是计算该点对所有多面体表面的立体角，非零则在内部。
3. GPU并行批量查询
在图形渲染、粒子系统等场景中，需同时判定数百万个点。可将多边形边数据上传为纹理/SSBO，在GPU着色器中并行执行射线法，利用SIMD架构获得数量级加速。
4. 浮点精度处理
实际工程中需注意：
– 使用 epsilon 容差判断共线和相等（如 $10^{-9}$）
– 避免在顶点正上方发射射线，通过微小扰动消除退化情况
– 对于地理坐标（经纬度），需考虑球面几何而非平面几何
总结
点在多边形判定虽看似简单，实则蕴含丰富的几何与算法思想。射线法凭借简洁直观的奇偶规则成为最常用方案；绕数法则通过角度累积提供了对自相交多边形的精确语义。两种算法时间复杂度均为 $O(n)$，配合包围盒预筛选可应对绝大多数实际场景。在Java实现中，精确处理边界情况和浮点精度是确保算法鲁棒性的关键。这一基础算法也是更复杂的空间索引（R-Tree、Geohash）和计算几何算法的基石。
