4

ESRI シェープファイルから読み取った 3D ジオメトリと多くの交差を実行するために使用できるオープン ソース ツールまたはライブラリ (理想的には Python) はありますか? テストのほとんどは、単純な線分とポリゴンです。

私は OGR 1.7.1 / GEOS 3.2.0 を調べました。データは正しくロードされますが、結果の交差は正しくなく、利用可能な他のツールのほとんどはこの作業に基づいているようです。

CGAL は代替手段でしたが、そのライセンスは適切ではありません。Boost ジェネリック ジオメトリ ライブラリは素晴らしく見えますが、API は巨大であり、そのままでは wkt または wkb リーダーをサポートしていないようです。

4

2 に答える 2

7

約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))
于 2011-01-10T22:17:14.253 に答える