使用光线投影算法（Ray Casting Algorithm，RCA）解决PiP问题（判断点是否在多边形内部）

高阶近似

于 2023-11-17 20:11:05 发布

阅读量294

点赞数 1

分类专栏：算法文章标签：算法 python

本文链接：https://blog.csdn.net/weixin_45937291/article/details/134470113

版权

算法专栏收录该内容

7 篇文章 2 订阅

订阅专栏

文章目录

背景
求解算法
代码实现

背景

PiP（point-in-polygon）问题是指判断点是否在多边形内部，这是地理信息系统 (Geographic Information System，GIS) 的基本操作。如下图所示，灰色区域代表多边形，红点位于多边形外部，绿点位于多边形内部，蓝点位于边界上。

请添加图片描述

点是否在多边形内部，肉眼上很容易判断，但是使用算法来解决却不是很简单。本文介绍了一种解决PiP问题的算法，并给出了其Python代码实现。算法主要分为以下两步：

首先判断该点是否在多边形的最小外接矩形（Minimum Bounding Rectangle， MBR）内
如果该点在MBR外部，则该点一定在多边形外部；否则，使用PiP算法（如Ray Casting Algorithm (RCA)）进一步判断该点在多边形内部、边界还是外部

求解算法

MBR算法

直观上，该算法首先使用一个矩阵框住整个多边形，然后判断该点是否在矩形内，如下图所示

请添加图片描述

该矩形的位置可以通过取多边形所有顶点 $x$ 和 $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_2,y_2)$ 是多边形一个边的两个端点，设该边所在直线的方程为 $A x + B y + C = 0$ ，则有
$\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}$
则点 $x_0,y_0)$ 在边上，当前仅当下面3个条件同时满足：

$\min(x_1,x_2)\leq x_0\leq\max(x_1,x_2)$
$\min(y_1,y_2)\leq y_0\leq \max(y_1,y_2)$
$A\cdot x_0+B\cdot y_0+C=0$

判断光线是否水平向右穿过边

设 $x_1, y_1)$ 和 $x_2,y_2)$ 是多边形一个边的两个端点，设该边所在直线的方程为 $A x + B y + C = 0$ ，测试点坐标为 $x_0,y_0)$

如果 $y_1$ 等于 $y_2$ ，则光线水平向右穿过边，当且仅当 $x_0 \leq \max(x_1,x_2)$ 并且 $y_0$ 等于 $y_1$
如果 $y_1$ 不等于 $y_2$ ，可以算出光线与边所在直线的交点坐标为 $(\frac{-C - B\cdot y_0}{A}, y_0)$ ，则光线水平向右穿过边，当且仅当 $\frac{-C - B\cdot y_0}{A} \ge x_0$ 并且交点在该边上。

完整代码

输入文件

有两个输入文件，首先是 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()

运行结果

请添加图片描述

高阶近似

关注

1
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
使用光线投影算法（Ray Casting Algorithm，RCA）解决PiP问题（判断点是否在多边形内部）

PiP（point-in-polygon）问题是指判断点是否在多边形内部，这是地理信息系统 (Geographic Information System，GIS) 的基本操作。点是否在多边形内部，肉眼上很容易判断，但是使用算法来解决却不是很简单。本文介绍了一种解决PiP问题的算法，并给出了其Python代码实现。
复制链接

扫一扫