[SciPy-User] 3d convex hull

[Pauli Virtanen]

Wed, 24 Aug 2011 23:56:36 +0800, Ning Guo wrote:
> I also want to try Delaunay function. But I cannot get enough info from
> the documentation. I want to output the Delaunay tetrahedral in 3D and
> need the vertex indices, facet areas and normals. How can I use the
> function in scipy.spatial? Now I have all the points with id and
> position.

Suppose you have 20 points in 3-D:

    import numpy as np
    import scipy.spatial

    points = np.random.rand(20, 3)
    tri = scipy.spatial.Delaunay(points)

The indices of the vertices of tetrahedron number `j` are
in `tri.vertices[j]`. The facet areas and normals can be computed
for each tetrahedron via vector cross products:

    tetra_points = tri.points[tri.vertices] # (N, 4, 3) array

    face_normals = np.empty_like(tetra_points)
    face_normals[:,0] = np.cross(tetra_points[:,0], tetra_points[:,1])
    face_normals[:,1] = np.cross(tetra_points[:,0], tetra_points[:,2])
    face_normals[:,2] = np.cross(tetra_points[:,0], tetra_points[:,3])
    face_normals[:,3] = np.cross(tetra_points[:,1], tetra_points[:,2])

    face_normal_lengths = np.sqrt(np.sum(face_normals**2, axis=2))

    face_normals /= face_normal_lengths[:,:,np.newaxis]
    face_areas = 0.5 * face_normal_lengths

One important point I don't know at the moment is if those normals
actually point away from the center of the tetrahedra. You'd have
to check the Qhull documentation to check whether they have a winding
convention that guarantees certain ordering of the vertices of
the simplices.

(Note: the above is untested code, so check it works first :)

-- 
Pauli Virtanen