python中3D多边形的交点

Question

是否有任何开源工具或库（最好是 python）可用于与从 ESRI shapefile 读取的 3D 几何体执行大量相交？大多数测试都是简单的线段与多边形。

我研究了 OGR 1.7.1 / GEOS 3.2.0，虽然它正确加载了数据，但生成的交叉点不正确，并且大多数其他可用工具似乎都建立在这项工作的基础上。

虽然 CGAL 是一种替代方案，但它的许可证并不合适。 Boost 通用几何库看起来很棒，但 API 很大，并且似乎不支持开箱即用的 wkt 或 wkb 阅读器。

Answer 1

近10年后更新！

这里是我 10 年前为我的光学光线追踪项目的第一个版本编写的一些代码。新版本代码不再使用多边形；我们现在用网格来做所有事情。

代码可以清理，有很多防御性复制。我不保证其准确性或有用性。我尝试将其从上面的 GitHub 链接中几乎逐字提取到一个独立的脚本中。

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

您可以尝试我最新的 lib Geometry3D，如下

pip install Geometry3D

所示

>>> 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()

也可以查看文档几何3D