文章目录
背景
PiP(point-in-polygon)问题是指判断点是否在多边形内部,这是地理信息系统 (Geographic Information System,GIS) 的基本操作。如下图所示,灰色区域代表多边形,红点位于多边形外部,绿点位于多边形内部,蓝点位于边界上。
点是否在多边形内部,肉眼上很容易判断,但是使用算法来解决却不是很简单。本文介绍了一种解决PiP问题的算法,并给出了其Python代码实现。算法主要分为以下两步:
-
首先判断该点是否在多边形的最小外接矩形(Minimum Bounding Rectangle, MBR)内
-
如果该点在MBR外部,则该点一定在多边形外部;否则,使用PiP算法(如Ray Casting Algorithm (RCA))进一步判断该点在多边形内部、边界还是外部
求解算法
MBR算法
直观上,该算法首先使用一个矩阵框住整个多边形,然后判断该点是否在矩形内,如下图所示
该矩形的位置可以通过取多边形所有顶点 x x x 和 y y y 坐标的最小值和最大值获得。判断一个点是否在矩形内也很简单,这里不再赘述。
使用MBR算法的目的是进行初步排查,在矩阵外部的点一定在多边形外部,但是在矩阵边缘或者内部的点则需要进一步判断
光线投影算法(Ray Casting Algorithm,RCA)
RCA算法的原理是从测试点开始朝一个方向画一条射线,计算它穿过多边形边界的次数,如果该线穿过边界奇数次,则该点位于多边形内部,否则该点位于多边形之外。下图给出一个示例,所有光线从测试点水平向右发出。
特殊情况
该算法理解起来并不难,计算光线穿过多边形边界的次数,也就是计算光线和多边形的边相交的次数。但是存在几种特殊情况:
-
点可能位于多边形的边界上,在这种情况下,无论穿过边界的次数为奇数次还是偶数次,该算法都应返回"边界"
-
光线和多边形的顶点或者边重合,下图给出了一个示例。解决方法是:当光穿过多边形一个边的某个顶点时,如果该边的另外一个顶点在该光线之下,这次穿过才计数
代码实现
伪代码
这里首先给出整个算法的伪代码实现
outside_points = []
boundary_points = []
inside_points = []
# MBR算法
in_MBR_points = []
for p in test_points:
if p在MBR外部:
把p加入到列表outside_points中
else:
把p加入到列表in_MBR_points中
# RCA算法
for p in in_MBR_points:
pass_num = 0
for line in polygon.lines():
if p在line上:
把p加入到列表boundary_points中
break
else:
if 光线水平向右穿过line:
p1,p2为line的两个端点
if 光线穿过p1并且p2在光线之下:
pass_num += 1
elif 光线穿过p2并且p1在光线之下:
pass_num += 1
elif 光线穿过线段但是不经过端点:
pass_num += 1
if p不在列表boundary_points中:
if pass_num是奇数:
把p加入到列表inside_points中
else:
把p加入到列表outside_points中
关键步骤的实现原理
判断点是否在多边形的边上
设
(
x
1
,
y
1
)
(x_1, y_1)
(x1,y1) 和
(
x
2
,
y
2
)
(x_2,y_2)
(x2,y2) 是多边形一个边的两个端点,设该边所在直线的方程为
A
x
+
B
y
+
C
=
0
Ax+By+C=0
Ax+By+C=0,则有
A
=
y
2
−
y
1
B
=
x
1
−
x
2
C
=
x
2
⋅
y
1
−
x
1
⋅
y
2
\begin{align} A&=y_2-y_1\\ B&=x_1-x_2 \\ C&=x_2\cdot y_1-x_1\cdot y_2 \end{align}
ABC=y2−y1=x1−x2=x2⋅y1−x1⋅y2
则点
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0) 在边上,当前仅当下面3个条件同时满足:
- min ( x 1 , x 2 ) ≤ x 0 ≤ max ( x 1 , x 2 ) \min(x_1,x_2)\leq x_0\leq\max(x_1,x_2) min(x1,x2)≤x0≤max(x1,x2)
- min ( y 1 , y 2 ) ≤ y 0 ≤ max ( y 1 , y 2 ) \min(y_1,y_2)\leq y_0\leq \max(y_1,y_2) min(y1,y2)≤y0≤max(y1,y2)
- A ⋅ x 0 + B ⋅ y 0 + C = 0 A\cdot x_0+B\cdot y_0+C=0 A⋅x0+B⋅y0+C=0
判断光线是否水平向右穿过边
设 ( x 1 , y 1 ) (x_1, y_1) (x1,y1) 和 ( x 2 , y 2 ) (x_2,y_2) (x2,y2) 是多边形一个边的两个端点,设该边所在直线的方程为 A x + B y + C = 0 Ax+By+C=0 Ax+By+C=0,测试点坐标为 ( x 0 , y 0 ) (x_0,y_0) (x0,y0)
- 如果 y 1 y_1 y1 等于 y 2 y_2 y2,则光线水平向右穿过边,当且仅当 x 0 ≤ max ( x 1 , x 2 ) x_0 \leq \max(x_1,x_2) x0≤max(x1,x2) 并且 y 0 y_0 y0 等于 y 1 y_1 y1
- 如果 y 1 y_1 y1 不等于 y 2 y_2 y2,可以算出光线与边所在直线的交点坐标为 ( − C − B ⋅ y 0 A , y 0 ) (\frac{-C - B\cdot y_0}{A}, y_0) (A−C−B⋅y0,y0),则光线水平向右穿过边,当且仅当 − C − B ⋅ y 0 A ≥ x 0 \frac{-C - B\cdot y_0}{A} \ge x_0 A−C−B⋅y0≥x0 并且交点在该边上。
完整代码
输入文件
有两个输入文件,首先是 polygon.csv
,表示多边形的顶点,如下所示
id,x,y
1,0,1
2,0,6
3,1,7
4,3,7
5,4,6
6,4,4
7,3,4
8,3,5
9,2,6
10,1,5
11,1,2
12,2,1
13,3,2
14,2,2
15,2,3
16,4,3
17,4,1
18,3,0
19,1,0
20,0,1
然后是 input.csv
,表示需要测试的顶点,如下所示
id,x,y
1,1.0,5.5
2,-1.0,0.5
3,2.5,3.5
4,1.5,5.5
5,0.5,5.5
6,2.5,0.0
7,0.5,1.0
8,3.5,0.0
9,4.0,7.5
10,-1.0,5.5
11,2.5,6.0
12,-0.5,4.5
13,2.0,4.0
14,1.5,2.0
15,0.0,3.5
16,3.0,1.5
17,1.0,-1.0
18,3.0,0.5
19,3.0,-1.0
20,5.0,1.0
21,4.0,2.5
22,4.0,2.0
23,2.5,2.0
24,3.0,2.0
25,5.0,6.0
26,3.5,2.5
27,-0.5,6.5
28,0.0,1.0
29,4.0,6.5
30,-1.0,1.5
31,0.0,6.0
32,3.5,5.0
33,1.5,-1.0
34,-0.5,0.0
35,0.0,2.0
36,0.0,4.0
37,-0.5,1.5
38,0.5,5.0
39,2.0,8.0
40,1.5,2.5
41,2.0,4.5
42,0.5,8.0
43,4.5,-1.0
44,3.0,3.5
45,5.0,0.5
46,-0.5,3.0
47,-0.5,-0.5
48,4.0,0.0
49,4.5,2.5
50,4.5,5.5
51,5.0,4.0
52,4.0,8.0
53,2.5,3.0
54,2.0,7.5
55,5.0,7.5
56,1.5,0.0
57,1.0,7.5
58,0.0,0.5
59,2.0,0.5
60,0.0,5.0
61,5.0,5.5
62,1.5,4.0
63,3.5,-0.5
64,-0.5,2.0
65,5.0,6.5
66,1.5,1.5
67,4.0,0.5
68,1.5,5.0
69,2.5,8.0
70,0.0,-1.0
71,2.5,-0.5
72,4.0,4.0
73,0.0,0.0
74,-1.0,2.0
75,3.0,8.0
76,2.0,6.5
77,1.0,0.0
78,-1.0,7.0
79,1.5,6.0
80,4.0,6.0
81,4.5,7.5
82,0.0,7.5
83,4.5,0.0
84,2.5,-1.0
85,3.5,4.0
86,0.0,4.5
87,1.0,8.0
88,1.5,3.5
89,3.5,2.0
90,4.0,1.0
91,5.0,1.5
92,1.0,4.5
93,2.5,6.5
94,-1.0,6.0
95,2.0,1.0
96,3.0,1.0
97,5.0,4.5
98,4.0,5.5
99,-1.0,-1.0
100,1.5,4.5
geometry.py
首先是 geometry.py
文件,该文件主要是定义了一些基本的类,比如表示点、边和多边形的类
import matplotlib.pyplot as plt
import copy
class Geometry:
def __init__(self, name):
self.__name = name
def get_name(self):
return self.__name
class Point(Geometry):
def __init__(self, name, x, y):
super().__init__(name)
self._x = x
self._y = y
def get_pos(self):
return self._x, self._y
def __eq__(self, point):
return self.get_name() == point.get_name() and point.get_x() == self._x and point.get_y() == self._y
class Line(Geometry):
def __init__(self, name, point_1, point_2):
super().__init__(name)
self.__p1 = point_1
self.__p2 = point_2
def get_points(self):
return self.__p1, self.__p2
def point_in_line(self, p):
p1x, p1y = self.__p1.get_pos()
p2x, p2y = self.__p2.get_pos()
p0x, p0y = p.get_pos()
A, B, C = self.__get_line_equation()
if Line.__is_mid(p0x, p1x, p2x) and Line.__is_mid(p0y, p1y, p2y):
return A * p0x + B * p0y + C == 0
return False
def is_intersects_right_ray(self, p):
# if x0 < max([x1,x2]) and disjoint_is_right_line(line, p) and PiP.__is_mid(y0, y1, y2):
p0x, p0y = p.get_pos()
A, B, C = self.__get_line_equation()
if A == 0:
p1x, p1y = self.__p1.get_pos()
p2x, _ = self.__p2.get_pos()
return p0x <= max([p1x, p2x]) and p1y == p0y
else:
# 交点坐标
intersect_x, intersect_y = (-C - B * p0y) / A, p0y
# 交点在发射点的右侧且交点在线段上
return intersect_x >= p0x and self.point_in_line(Point('', intersect_x, intersect_y))
def __is_mid(x0, x1, x2):
'''
返回x0的值是否在x1和x2之间
'''
return min([x1, x2]) <= x0 <= max([x1, x2])
def __get_line_equation(self):
p1x, p1y = self.__p1.get_pos()
p2x, p2y = self.__p2.get_pos()
A = p2y - p1y
B = p1x - p2x
C = p2x * p1y - p1x * p2y
return A, B, C
class Polygon(Geometry):
def __init__(self, name, points):
super().__init__(name)
self.__points = copy.copy(points)
def get_points(self):
return self.__points
def lines(self):
res = []
points = self.get_points()
point_a = points[0]
for point_b in points[1:]:
res.append(Line(point_a.get_name() + '-' + point_b.get_name(), point_a, point_b))
point_a = point_b
res.append(Line(point_a.get_name() + '-' + points[0].get_name(), point_a, points[0]))
return res
def get_MBR(self):
'''
min_x, max_x, min_y, max_y
'''
min_x, min_y = self.__points[0].get_pos()
max_x, max_y = min_x, min_y
for p in self.__points[1:]:
xi, yi = p.get_pos()
if xi < min_x:
min_x = xi
if xi > max_x:
max_x = xi
if yi < min_y:
min_y = yi
if yi > max_y:
max_y = yi
return min_x, max_x, min_y, max_y
def display(self):
xs = []
ys = []
for p in self.__points:
x, y = p.get_pos()
xs.append(x)
ys.append(y)
plt.plot(xs, ys)
utils.py
接下来是 utils.py
,定义一些辅助函数,如读取csv文件。
import csv
from geometry import Point
def read_points_from_csv(file_path):
points = []
with open(file_path, 'r') as csv_file:
csv_reader = csv.reader(csv_file)
next(csv_reader) # 跳过Header
for row in csv_reader:
points.append(Point(row[0], float(row[1]), float(row[2])))
return points
PiP.py
最后是 PiP.py
文件,用于求解PiP问题
import matplotlib.pyplot as plt
import utils
import geometry
class PiP:
def __init__(self, polygon, test_points):
self.__polygon = polygon
self.__test_p = test_points
self.__boundary_p = []
self.__inside_p = []
self.__outside_p = []
self.__in_MBR_p = []
def solve(self):
self.__MBR()
self.__RCA()
self.display()
def __MBR(self):
min_x, max_x, min_y, max_y = polygon.get_MBR()
for p in self.__test_p:
px, py = p.get_pos()
if min_x <= px <= max_x and min_y <= py <= max_y:
self.__in_MBR_p.append(p)
else:
self.__outside_p.append(p)
def __RCA(self):
for p in self.__in_MBR_p:
pass_num = 0
_, y0 = p.get_pos()
for line in polygon.lines():
p1, p2 = line.get_points()
_, y1 = p1.get_pos()
_, y2 = p2.get_pos()
if line.point_in_line(p):
self.__boundary_p.append(p)
else:
if line.is_intersects_right_ray(p):
# 光线向右穿过line
if y0 == y1 and y2 < y0:
pass_num += 1
elif y0 == y2 and y1 < y0:
pass_num += 1
elif y0 != y1 and y0 != y2:
pass_num += 1
if p not in self.__boundary_p:
if pass_num % 2:
self.__inside_p.append(p)
else:
self.__outside_p.append(p)
def display(self):
self.__polygon.display()
self.__display_points(self.__boundary_p, 'bo')
self.__display_points(self.__inside_p, 'go')
self.__display_points(self.__outside_p, 'ro')
def __display_points(self, points, style='bo'):
for p in points:
x, y = p.get_pos()
plt.plot([x], [y], style)
if __name__ == '__main__':
poly_points = utils.read_points_from_csv('./polygon.csv')
test_points = utils.read_points_from_csv('./input.csv')
polygon = geometry.Polygon('G', poly_points)
pip_solver = PiP(polygon, test_points)
pip_solver.solve()
plt.savefig('./result.png', dpi=400)
plt.show()