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):
self.__p0 = p0
self.__p1 = p1
self.__x = p1.x - p0.x
self.__y = p1.y - p0.y
self.__z = p1.z - p0.z
@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):
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):
v = GeoVector(p0, p1)
u = GeoVector(p0, p2)
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):
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):
def __init__(self, ptList, idxList):
self.__v = []
self.__idx = []
self.__n = len(ptList)
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):
self.__v = []
self.__idx = []
self.__n = len(ptList)
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):
self.__MaxUnitMeasureError = 0.001;
self.__Set3DPolygonBoundary(polygonInst)
self.__Set3DPolygonUnitError()
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):
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):
faceVerticeIndex = []
fpOutward = []
vertices = polygon.v
n = polygon.n
for i in range(0, n):
p0 = vertices[i]
for j in range(i + 1, n):
p1 = vertices[j]
for k in range(j + 1, n):
p2 = vertices[k]
trianglePlane = GeoPlane.Create(p0, p1, p2)
onLeftCount = 0
onRightCount = 0
pointInSamePlaneIndex = []
for l in range(0, n):
if(l != i and l != j and l != k):
pt = vertices[l]
dis = trianglePlane * pt
if(abs(dis) < self.MaxDisError):
pointInSamePlaneIndex.append(l)
else:
if(dis < 0):
onLeftCount = onLeftCount + 1
else:
onRightCount = onRightCount + 1
if(onLeftCount == 0 or onRightCount == 0):
verticeIndexInOneFace = [];
verticeIndexInOneFace.append(i);
verticeIndexInOneFace.append(j);
verticeIndexInOneFace.append(k);
m = len(pointInSamePlaneIndex);
if(m > 0):
for p in range(0, m):
verticeIndexInOneFace.append(pointInSamePlaneIndex[p])
if ( not Utility.ContainsList
(faceVerticeIndex, verticeIndexInOneFace)):
faceVerticeIndex.append(verticeIndexInOneFace)
if (onRightCount == 0):
fpOutward.append(trianglePlane)
else:
if (onLeftCount == 0):
fpOutward.append(-trianglePlane)
self.__NumberOfFaces = len(faceVerticeIndex)
self.__Faces = []
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))
v = []
idx = []
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(dis > 0):
return False
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:
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)
rbFile = open('../LASInput/cube.las', 'rb')
wbFile = open('../LASOutput/cube_inside.las', 'wb')
wtFile = open('../LASOutput/cube_inside.txt', 'w')
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)
rbFile.seek(0)
buf = rbFile.read(point_offset[0])
wbFile.write(buf)
insideCount = 0
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)):
rbFile.seek(record_loc)
buf = rbFile.read(record_len[0])
wbFile.write(buf)
wtFile.write("%15.6f %15.6f %15.6f\n" % ( x, y, z) )
insideCount = insideCount + 1
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