python - python 中 3D 多边形的交点

Question

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

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

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

Answer 1

近10年后的更新!

这里是 some code I wrote 10 years ago对于我的光线追踪项目的第一个版本。新版代码不再使用多边形；我们现在用网格做所有事情。

代码可以清理，有很多防御性复制。我不对其准确性或实用性提供任何保证。我试图将它从上面的 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))