python-cython-多边形裁剪点云（python的最快代码，比shapely快至少100倍）

龙格、库塔的仆人

已于 2024-01-09 14:52:38 修改

阅读量486

点赞数 8

文章标签： python 开发语言

于 2024-01-05 16:18:08 首次发布

本文链接：https://blog.csdn.net/loceccc/article/details/135412137

版权

本文介绍了如何使用Cython和numpy库在一个test.pyx文件中实现点是否在多边形内的快速判断，以及如何通过setup.py进行编译。实测处理1.1亿个点的速度达到6.9秒，展示了Cython在大数据处理中的优势。

摘要由CSDN通过智能技术生成

1.新建一个test.pyx文件，下面是代码（polygon为list，pts为numpy）

import numpy as np
cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)
# polygon = [[502990.4322, 2518211.1518], [502984.9641, 2518202.5316],[502992.6705, 2518199.6960], [502997.7703, 2518208.2879]]
# pts=numpy.array([[502990.4322, 2518211.1518,0],[502990.4322, 2518211.1518,0]])	
def is_point_inside_polygon_cython(polygon, pts):
    polygon=polygon.astype(np.float)
    cdef double[:,:] polygon_view=polygon
    pts = pts.astype(np.float)
    cdef double[:, :] pts_view = pts
    cdef int n = polygon.shape[0]
    cdef int m = pts.shape[0]
    cdef double xints
    cdef int inside
    cdef double p1x, p1y, p2x, p2y
    cdef double x, y,z
    cdef Py_ssize_t k
    cdef Py_ssize_t i
    data =[]

    for k in range(m):
        x = pts_view[k,0]
        y = pts_view[k,1]
        z = pts_view[k,2]
        inside = False
        p1x, p1y = polygon_view[0, 0], polygon_view[0, 1]

        for i in range(n + 1):
            p2x, p2y = polygon_view[i % n, 0], polygon_view[i % n, 1]

            if y > min(p1y, p2y) and y <= max(p1y, p2y) and x <= max(p1x, p2x):
                if p1y != p2y:
                    xints =(y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside

            p1x, p1y = p2x, p2y

        if inside:
            data.append([x, y,z])

    return data

2.新建一个setup.py

from setuptools import setup
from Cython.Build import cythonize
import numpy as np

setup(
    ext_modules=cythonize("test.pyx"),
    include_dirs=[r"C:\Users\Huangchao\AppData\Local\Programs\Python\Python37\Lib\site-packages\numpy\core\include"]
    # include_dirs为自己的numpy里include路径
)

3.setup.py和test.pyx放在一个文件夹，文件夹路径打开powershell，输入下面代码进行编译。

 python setup.py build_ext --inplace

4.运行代码，实测1.1亿个点，速度为6.9s

import time
import numpy as np
import laspy
import test


def read_las_file(file_path):  # 读取las格式数据，生成numpy的格式
    las = laspy.read(file_path)
    xyz = las.xyz
    xy = xyz[:, :2]
    return xy, xyz  # 返回点云的xy，xyz
path1 = r"0625-1右.las"	# 自己的las格式点云
point_cloud1 = read_las_file(path1)  # 读取点云

# 示例使用
polygon = [[502990.4322, 2518211.1518], [502984.9641, 2518202.5316],
           [502992.6705, 2518199.6960], [502997.7703, 2518208.2879]]

t1 = time.time()
data3 = test.is_point_inside_polygon_cython(np.array(polygon), point_cloud1[1])  # 只在范围内的点云
t2 = time.time()
print(len(data3), t2 - t1)