Point Inside 3D Convex Polygon in Python

• Download PythonSLASProc - 5.5 MB
• Download FugroViewerSetup22 - 2.9 MB

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

• 20^th January, 2016: Initial version