Python 中三维凸多边形内的点
5.00/5 (2投票s)
一个算法，用于确定给定多边形顶点时，一个点是否位于给定三维凸多边形内。
引言
这是原始 C++ 算法的 Python 版本，该算法可以在 这里找到。
背景
请参考原始 C++ 算法 这里。
Using the Code
该算法被封装成一个 Python 类库文件夹 GeoProc。主测试程序 PythonLASProc 从 LAS 文件读取点云数据，然后将所有内部点写入一个新的 LAS 文件。可以使用免费软件 FugroViewer 显示 LAS 文件。
要使用 GeoProc 库，首先准备一个 GeoPolygonProc 实例，如下所示
gp = [p1, p2, p3, p4, p5, p6, p7, p8];
gpInst = GeoPolygon(gp);
procInst = GeoPolygonProc(gpInst);
在上面的代码中，gp 是一个 GeoPoint 数组；gpInst 是从 GeoPoint 数组创建的 GeoPolygon 实例，作为其顶点；procInst 是从 GeoPolygon 实例创建的 GeoPolygonProc 实例。
在主测试程序 PythonLASProc.py 中，pronInst 首先用于检查 3D 多边形边界，以预先筛选内部点进行粗略筛选，然后使用其主要的 public 方法 PointInside3DPolygon 获取实际的内部点，如下所示
if (x > procInst.x0 and x < procInst.x1 and
y > procInst.y0 and y < procInst.y1 and
z > procInst.z0 and z < procInst.z1):
if (procInst.PointInside3DPolygon(x, y, z)):
以下是所有 Python GeoProc 类
代码几乎不言自明，并带有注释。
GeoPoint.py
Python 具有运算符重载，public 函数 __add__ 是一个 运算符 + 重载，在 GeoVector.py 中引用。
class GeoPoint(object):
def __init__(self, x, y, z):
self.x = x
self.y = y
self.z = z
def __add__(self, p):
return GeoPoint(self.x + p.x, self.y + p.y, self.z + p.z)
GeoVector.py
该类声明了两个 读写属性 和三个 只读属性，public 函数 __mul__ 是一个乘法 运算符 * 重载，在 GeoPlane.py 中引用。
from GeoPoint import *
class GeoVector(object):
def __init__(self, p0, p1): # GeoPoint p0, p1
self.__p0 = p0 # read write (get set)
self.__p1 = p1 # read write (get set)
self.__x = p1.x - p0.x # read only
self.__y = p1.y - p0.y # read only
self.__z = p1.z - p0.z # read only
@property
def p0(self):
return self.__p0
@p0.setter
def p0(self, p0):
self.__p0 = p0
self.__x = self.p1.x - p0.x
self.__y = self.p1.y - p0.y
self.__z = self.p1.z - p0.z
@property
def p1(self):
return self.__p1
@p1.setter
def p1(self, p1):
self.__p1 = p1
self.__x = p1.x - self.p0.x
self.__y = p1.y - self.p0.y
self.__z = p1.z - self.p0.z
@property
def x(self):
return self.__x
@property
def y(self):
return self.__y
@property
def z(self):
return self.__z
def __mul__(self, v): # v: GeoVector
x = self.y * v.z - self.z * v.y
y = self.z * v.x - self.x * v.z
z = self.x * v.y - self.y * v.x
p0 = v.p0
p1 = p0 + GeoPoint(x, y, z)
return GeoVector(p0, p1)
GeoPlane.py
该类声明了 4 个 公共数据成员，一个公共 静态方法 Create，一个 一元运算符 - 重载 __neg__，一个乘法运算符重载 __mul__，这些方法在 GeoPolygonProc.py 中引用。
from GeoVector import *
class GeoPlane(object):
def __init__(self, a, b, c, d):
self.a = a
self.b = b
self.c = c
self.d = d
@staticmethod
def Create(p0, p1, p2): # p0, p1, p2: GeoPoint
v = GeoVector(p0, p1)
u = GeoVector(p0, p2)
# normal vector
n = u * v;
a = n.x
b = n.y
c = n.z
d = - (a * p0.x + b * p0.y + c * p0.z)
return GeoPlane(a, b, c, d)
def __neg__(self):
return GeoPlane(-self.a, -self.b, -self.c, -self.d)
def __mul__(self, pt): # pt: GeoPoint
return (self.a * pt.x + self.b * pt.y + self.c * pt.z + self.d)
GeoFace.py
该类声明了两个 只读数组属性。
class GeoFace(object):
# ptList: vertice GeoPoint Array
# idxList: vertice index Integer Array
def __init__(self, ptList, idxList):
#alloc instance array memory
self.__v = [] # vertice point array
self.__idx = [] # vertice index array
self.__n = len(ptList) # number of vertices
for i in range(0, self.__n):
self.__v.append(ptList[i])
self.__idx.append(idxList[i])
@property
def v(self):
return self.__v
@property
def idx(self):
return self.__idx
@property
def n(self):
return self.__n
GeoPolgyon.py
该类声明了两个 只读数组属性。
class GeoPolygon(object):
def __init__(self, ptList): # ptList: vertice GeoPoint Array
#alloc instance array memory
self.__v = [] # vertice point array
self.__idx = [] # vertice index array
self.__n = len(ptList) # number of vertices
for pt in ptList:
self.__v.append(pt)
self.__idx.append(ptList.index(pt))
@property
def v(self):
return self.__v
@property
def idx(self):
return self.__idx
@property
def n(self):
return self.__n
Utility.py
该类声明了一个 static 方法，用于检查 2d 数组是否包含 1d 数组作为其项。此方法用于避免在 GeoPolygonProc.py 中添加重复的面。
class Utility(object):
@staticmethod
def ContainsList(listList, listItem):
listItem.sort()
for i in range(0, len(listList)):
temp = listList[i]
if(len(temp) == len(listItem)):
temp.sort()
itemEqual = True
for j in range(0, len(listItem)):
if(temp[j] != listItem[j]):
itemEqual = False
break
if(itemEqual):
return True
else:
return False
return False
GeoPolygonProc.py
请参考原始 C++ 版本的相同部分以获取说明。
from Utility import *
from GeoPlane import *
from GeoPoint import *
from GeoFace import *
class GeoPolygonProc(object):
def __init__(self, polygonInst): # polygonInst: GeoPolygon
self.__MaxUnitMeasureError = 0.001;
# Set boundary
self.__Set3DPolygonBoundary(polygonInst)
# Set maximum point to face plane distance error,
self.__Set3DPolygonUnitError()
# Set faces and face planes
self.__SetConvex3DFaces(polygonInst)
@property
def x0(self):
return self.__x0
@property
def x1(self):
return self.__x1
@property
def y0(self):
return self.__y0
@property
def y1(self):
return self.__y1
@property
def z0(self):
return self.__z0
@property
def z1(self):
return self.__z1
@property
def NumberOfFaces(self):
return self.__NumberOfFaces
@property
def FacePlanes(self):
return self.__FacePlanes
@property
def Faces(self):
return self.__Faces
def __Set3DPolygonBoundary(self, polygon): # polygon: GeoPolygon
# List<GeoPoint>
vertices = polygon.v;
n = polygon.n;
xmin = vertices[0].x
xmax = vertices[0].x
ymin = vertices[0].y
ymax = vertices[0].y
zmin = vertices[0].z
zmax = vertices[0].z
for i in range(1, n):
if (vertices[i].x < xmin): xmin = vertices[i].x
if (vertices[i].y < ymin): ymin = vertices[i].y
if (vertices[i].z < zmin): zmin = vertices[i].z
if (vertices[i].x > xmax): xmax = vertices[i].x
if (vertices[i].y > ymax): ymax = vertices[i].y
if (vertices[i].z > zmax): zmax = vertices[i].z
self.__x0 = xmin
self.__x1 = xmax
self.__y0 = ymin
self.__y1 = ymax
self.__z0 = zmin
self.__z1 = zmax
def __Set3DPolygonUnitError(self):
self.MaxDisError = ((abs(self.__x0) + abs(self.__x1) +
abs(self.__y0) + abs(self.__y1) +
abs(self.__z0) + abs(self.__z1)) / 6 * self.__MaxUnitMeasureError)
def __SetConvex3DFaces(self, polygon): # polygon: GeoPolygon
# vertices indexes for all faces, 2d list
faceVerticeIndex = []
# face planes for all faces, 1d list
fpOutward = []
vertices = polygon.v
n = polygon.n
for i in range(0, n):
# triangle point 1
p0 = vertices[i]
for j in range(i + 1, n):
# triangle p
p1 = vertices[j]
for k in range(j + 1, n):
# triangle point 3
p2 = vertices[k]
trianglePlane = GeoPlane.Create(p0, p1, p2)
onLeftCount = 0
onRightCount = 0
# indexes of points that is in same plane with face triangle plane
pointInSamePlaneIndex = []
for l in range(0, n):
# any point other than the 3 triangle points
if(l != i and l != j and l != k):
pt = vertices[l]
dis = trianglePlane * pt
# next point is in the triangle plane
if(abs(dis) < self.MaxDisError):
pointInSamePlaneIndex.append(l)
else:
if(dis < 0):
onLeftCount = onLeftCount + 1
else:
onRightCount = onRightCount + 1
# This is a face for a CONVEX 3d polygon.
# For a CONCAVE 3d polygon, this maybe not a face.
if(onLeftCount == 0 or onRightCount == 0):
verticeIndexInOneFace = [];
# add 3 triangle plane vertices
verticeIndexInOneFace.append(i);
verticeIndexInOneFace.append(j);
verticeIndexInOneFace.append(k);
m = len(pointInSamePlaneIndex);
# add other vertices in the same triangle plane
if(m > 0):
for p in range(0, m):
verticeIndexInOneFace.append(pointInSamePlaneIndex[p])
# if verticeIndexInOneFace is a new face
if ( not Utility.ContainsList
(faceVerticeIndex, verticeIndexInOneFace)):
#print(verticeIndexInOneFace)
# add it in the faceVerticeIndex list
faceVerticeIndex.append(verticeIndexInOneFace)
# add the trianglePlane in the face plane list fpOutward.
if (onRightCount == 0):
fpOutward.append(trianglePlane)
else:
if (onLeftCount == 0):
fpOutward.append(-trianglePlane)
#else:
# possible reasons :
# 1. the plane is not a face of a convex 3d polygon,
# it is a plane crossing the convex 3d polygon.
# 2. the plane is a face of a concave 3d polygon
# number of faces
self.__NumberOfFaces = len(faceVerticeIndex)
# face list
self.__Faces = []
# face planes list
self.__FacePlanes = []
for i in range(0, self.__NumberOfFaces):
self.__FacePlanes.append(GeoPlane(fpOutward[i].a, fpOutward[i].b,
fpOutward[i].c, fpOutward[i].d))
# face vertices
v = []
# face vertices indexes
idx = []
# number of vertices of the face
count = len(faceVerticeIndex[i])
for j in range(0, count):
idx.append(faceVerticeIndex[i][j])
v.append(GeoPoint(vertices[ idx[j] ].x,
vertices[ idx[j] ].y,
vertices[ idx[j] ].z))
self.__Faces.append(GeoFace(v, idx))
def PointInside3DPolygon(self, x, y, z):
pt = GeoPoint(x, y, z)
for i in range(0, self.__NumberOfFaces):
dis = self.__FacePlanes[i] * pt
# If the point is in the same half space
# with normal vector for any face of the cube,
# then it is outside of the 3D polygon
if(dis > 0):
return False
# If the point is in the opposite half space with normal vector for all 6 faces,
# then it is inside of the 3D polygon
return True
PythonLASProc.py
unpack 和 pack 用于二进制文件访问，类似于 PHP。
import struct
import sys
sys.path.append('../GeoProc')
from GeoPoint import *
from GeoPolygon import *
from GeoPolygonProc import *
print('Start Processing ...')
try:
# prepare GeoPolygonProc instance
p1 = GeoPoint( - 27.28046, 37.11775, - 39.03485)
p2 = GeoPoint( - 44.40014, 38.50727, - 28.78860)
p3 = GeoPoint( - 49.63065, 20.24757, - 35.05160)
p4 = GeoPoint( - 32.51096, 18.85805, - 45.29785)
p5 = GeoPoint( - 23.59142, 10.81737, - 29.30445)
p6 = GeoPoint( - 18.36091, 29.07707, - 23.04144)
p7 = GeoPoint( - 35.48060, 30.46659, - 12.79519)
p8 = GeoPoint( - 40.71110, 12.20689, - 19.05819)
gp = [p1, p2, p3, p4, p5, p6, p7, p8];
gpInst = GeoPolygon(gp)
procInst = GeoPolygonProc(gpInst)
# open files
rbFile = open('../LASInput/cube.las', 'rb')
wbFile = open('../LASOutput/cube_inside.las', 'wb')
wtFile = open('../LASOutput/cube_inside.txt', 'w')
# read LAS public header
rbFile.seek(96)
r_point_offset = rbFile.read(4)
r_variable_len_num = rbFile.read(4)
r_record_type = rbFile.read(1)
r_record_len = rbFile.read(2)
r_record_num = rbFile.read(4)
rbFile.seek(131)
r_x_scale = rbFile.read(8)
r_y_scale = rbFile.read(8)
r_z_scale = rbFile.read(8)
r_x_offset = rbFile.read(8)
r_y_offset = rbFile.read(8)
r_z_offset = rbFile.read(8)
point_offset = struct.unpack('I', r_point_offset)
record_len = struct.unpack('H', r_record_len)
record_num = struct.unpack('L', r_record_num)
x_scale = struct.unpack('d', r_x_scale)
y_scale = struct.unpack('d', r_y_scale)
z_scale = struct.unpack('d', r_z_scale)
x_offset = struct.unpack('d', r_x_offset)
y_offset = struct.unpack('d', r_y_offset)
z_offset = struct.unpack('d', r_z_offset)
# read LAS header
rbFile.seek(0)
buf = rbFile.read(point_offset[0])
# write LAS header
wbFile.write(buf)
insideCount = 0
# read points
for i in range(0, record_num[0]):
record_loc = int(point_offset[0]) + int(record_len[0]) * int(i)
rbFile.seek(record_loc)
xi = struct.unpack('l', rbFile.read(4))
yi = struct.unpack('l', rbFile.read(4))
zi = struct.unpack('l', rbFile.read(4))
x = (int(xi[0]) * x_scale[0]) + x_offset[0]
y = (int(yi[0]) * y_scale[0]) + x_offset[0]
z = (int(zi[0]) * z_scale[0]) + z_offset[0]
if (x > procInst.x0 and x < procInst.x1 and
y > procInst.y0 and y < procInst.y1 and
z > procInst.z0 and z < procInst.z1):
if (procInst.PointInside3DPolygon(x, y, z)):
# read point record
rbFile.seek(record_loc)
buf = rbFile.read(record_len[0])
# write binary point record
wbFile.write(buf)
#write text point record
wtFile.write("%15.6f %15.6f %15.6f\n" % ( x, y, z) )
insideCount = insideCount + 1
# update point number of ground LAS
wbFile.seek(107)
wbFile.write(struct.pack('i', insideCount))
finally:
wbFile.close()
wtFile.close()
rbFile.close()
print('Completed.');
关注点
Python 语法非常特殊，例如，缩进用于代码块边界。
历史
- 2016 年 1 月 20 日：初始版本