Point Inside 3D Convex Polygon in Python
5.00/5 (2 votes)
An algorithm to determine if a point is inside a 3D convex polygon for a given polygon vertices in Python
Introduction
This is a Python version of the original C++ algorithm which can be found here.
Background
Please refer to the original C++ algorithm here.
Using the Code
The algorithm is wrapped into a Python class library folder GeoProc. The main test program PythonLASProc reads point cloud data from a LAS file, then writes all inside points into a new LAS file. A LAS file is able to be displayed with a freeware FugroViewer.
To consume the GeoProc library, prepare a GeoPolygonProc instance first like this:
gp = [p1, p2, p3, p4, p5, p6, p7, p8];
gpInst = GeoPolygon(gp);
procInst = GeoPolygonProc(gpInst);
In the above code, gp is a GeoPoint array; gpInst is a GeoPolygon instance created from a GeoPoint array as its vertices; procInst is a GeoPolygonProc instance created from a GeoPolygon instance.
In the main test program PythonLASProc.py, the pronInst is used to check the 3D polygon boundary first to pre-filter the inside points for coarse screening, then use its main public method PointInside3DPolygon to get the actual inside points as below:
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)):
Here are all Python GeoProc classes:
The codes are almost self-evident with comments.
GeoPoint.py
Python has operator overloading, public function __add__ is a operator + overloading, which is referenced in 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
This class declares two read-write properties and three read-only properties, public function __mul__ is a multiple operator * overloading, which is referenced in 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
The class declares 4 public data members, one public static method Create, one unary operator - overloading __neg__, one multiple operator overloading __mul__, these methods are referenced in 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
The class declares two read-only array properties.
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
The class declares two read-only array properties.
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
The class delcares a static method to check if a 2d array constains a 1d array as its item. This method is used to avoid adding duplicate faces in 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
Please refer to the same section in original C++ version for the explanation.
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
The unpack and pack are used in binary file access which is similar with 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.');
Points of Interest
Python syntax is very special, i.e., indentation is used as code block boundary.
History
- 20th January, 2016: Initial version