Intersections of 3D polygons in python

Question

Are there any open source tools or libraries (ideally in python) that are available for performing lots of intersections with 3D geometry read from an ESRI shapefile? Most of the tests will be simple line segments vs polygons.

I've looked into OGR 1.7.1 / GEOS 3.2.0, and whilst it loads the data correctly, the resulting intersections aren't correct, and most of the other tools available seem to build on this work.

Whilst CGAL would have been an alternative, it's license isn't suitable. The Boost generic geometry library looks fantastic, but the api is huge, and doesn't seem to support wkt or wkb readers out of the box.

Answer 1

Update after nearly 10 years!

Here is some code I wrote 10 years ago for the first version of my optical ray-tracing project. The new version of the code no longer uses polygon; we do everything with meshes now.

The code could be cleared up, there is a lot of defensive copying. I don't give any guarantee of its accuracy or usefulness. I've tried to pull it almost verbatim from the GitHub link above into an isolated script.

import numpy as np

def cmp_floats(a,b, atol=1e-12):
    return abs(a-b) < atol


def magnitude(vector):
   return np.sqrt(np.dot(np.array(vector), np.array(vector)))


def norm(vector):
   return np.array(vector)/magnitude(np.array(vector))


class Ray(object):
    """A ray in the global cartesian frame."""
    def __init__(self, position, direction):
        self.position = np.array(position)
        self.direction = norm(direction) 


class Polygon(object):
    """ Polygon constructed from greater than two points.
    
        Only convex polygons are allowed! 
    
        Order of points is of course important!
    """
    
    def __init__(self, points):
        super(Polygon, self).__init__()
        self.pts = points
        #check if points are in one plane
        assert len(self.pts) >= 3, "You need at least 3 points to build a Polygon"
        if len(self.pts) > 3:
            x_0 = np.array(self.pts[0])
            for i in range(1,len(self.pts)-2):
                #the determinant of the vectors (volume) must always be 0
                x_i = np.array(self.pts[i])
                x_i1 = np.array(self.pts[i+1])
                x_i2 = np.array(self.pts[i+2])
                det = np.linalg.det([x_0-x_i, x_0-x_i1, x_0-x_i2])
                assert cmp_floats( det, 0.0 ), "Points must be in a plane to create a Polygon"
                

    def on_surface(self, point):
        """Returns True if the point is on the polygon's surface and false otherwise."""
        n = len(self.pts)
        anglesum = 0
        p = np.array(point)

        for i in range(n):
            v1 = np.array(self.pts[i]) - p
            v2 = np.array(self.pts[(i+1)%n]) - p

            m1 = magnitude(v1)
            m2 = magnitude(v2)

            if cmp_floats( m1*m2 , 0. ):
                return True #point is one of the nodes
            else:
                # angle(normal, vector)
                costheta = np.dot(v1,v2)/(m1*m2)
            anglesum = anglesum + np.arccos(costheta)
        return cmp_floats( anglesum , 2*np.pi )


    def contains(self, point):
        return False


    def surface_identifier(self, surface_point, assert_on_surface = True):
        return "polygon"


    def surface_normal(self, ray, acute=False):
        vec1 = np.array(self.pts[0])-np.array(self.pts[1])
        vec2 = np.array(self.pts[0])-np.array(self.pts[2])
        normal = norm( np.cross(vec1,vec2) )
        return normal


    def intersection(self, ray):
        """Returns a intersection point with a ray and the polygon."""
        n = self.surface_normal(ray)

        #Ray is parallel to the polygon
        if cmp_floats( np.dot( np.array(ray.direction), n ), 0. ):
            return None
 
        t = 1/(np.dot(np.array(ray.direction),n)) * ( np.dot(n,np.array(self.pts[0])) - np.dot(n,np.array(ray.position)) )
        
        #Intersection point is behind the ray
        if t < 0.0:
            return None

        #Calculate intersection point
        point = np.array(ray.position) + t*np.array(ray.direction)
        
        #Check if intersection point is really in the polygon or only on the (infinite) plane
        if self.on_surface(point):
            return [list(point)]

        return None


if __name__ == '__main__':
    points = [[0,0,0],[0,0.1,0],[0.1,0.1,-0.03],[0.1,0,-0.03]]
    polygon = Polygon(points)
    ray = Ray(position=(0,0,0), direction=(0,0,1))
    print(polygon.intersection(ray))

Answer 2

You can try my latest lib Geometry3D by

pip install Geometry3D

An example is shown below

>>> from Geometry3D import *
>>>
>>> a = Point(1,1,1)
>>> b = Point(-1,1,1)
>>> c = Point(-1,-1,1)
>>> d = Point(1,-1,1)
>>> e = Point(1,1,-1)
>>> f = Point(-1,1,-1)
>>> g = Point(-1,-1,-1)
>>> h = Point(1,-1,-1)
>>> cph0 = Parallelepiped(Point(-1,-1,-1),Vector(2,0,0),Vector(0,2,0),Vector(0,0,2))
>>> cpg12 = ConvexPolygon((e,c,h))
>>> cpg13 = ConvexPolygon((e,f,c))
>>> cpg14 = ConvexPolygon((c,f,g))
>>> cpg15 = ConvexPolygon((h,c,g))
>>> cpg16 = ConvexPolygon((h,g,f,e))
>>> cph1 = ConvexPolyhedron((cpg12,cpg13,cpg14,cpg15,cpg16))
>>> a1 = Point(1.5,1.5,1.5)
>>> b1 = Point(-0.5,1.5,1.5)
>>> c1 = Point(-0.5,-0.5,1.5)
>>> d1 = Point(1.5,-0.5,1.5)
>>> e1 = Point(1.5,1.5,-0.5)
>>> f1 = Point(-0.2,1.5,-0.5)
>>> g1 = Point(-0.2,-0.5,-0.5)
>>> h1 = Point(1.5,-0.5,-0.5)
>>>
>>> cpg6 = ConvexPolygon((a1,d1,h1,e1))
>>> cpg7 = ConvexPolygon((a1,e1,f1,b1))
>>> cpg8 = ConvexPolygon((c1,b1,f1,g1))
>>> cpg9 = ConvexPolygon((c1,g1,h1,d1))
>>> cpg10 = ConvexPolygon((a1,b1,c1,d1))
>>> cpg11 = ConvexPolygon((e1,h1,g1,f1))
>>> cph2 = ConvexPolyhedron((cpg6,cpg7,cpg8,cpg9,cpg10,cpg11))
>>> cph3 = intersection(cph0,cph2)
>>>
>>> cph4 = intersection(cph1,cph2)
>>> r = Renderer()
>>> r.add((cph0,'r',1),normal_length = 0)
>>> r.add((cph1,'r',1),normal_length = 0)
>>> r.add((cph2,'g',1),normal_length = 0)
>>> r.add((cph3,'b',3),normal_length = 0.5)
>>> r.add((cph4,'y',3),normal_length = 0.5)
>>> r.show()

You can also look at the documentation Geomrtry3D