Search code examples
python3dmeshopen3d

Create planar TriangleMesh from points in border


I'm using Python, Scipy and Open3D.

Problem:

I have N points (say, ~2000) that roughly correspond to a noisy plane. I managed to find all the Nv points corresponding to the points at the border of the plane. How can I quickly generate a TriangleMesh using these border points?

Here is my work so far:

Since I have a model for this plane, I flatten the whole pointcloud to make it fit the model.

After that, I can compute the convex hull directly using Open3D, but it's annoying because I cannot turn off the warnings. Instead, what I did is: using the model of the plane, I calculated the rotation matrix to align it with the XY plane, then used Scipy's ConvexHull implementation on the X and Y elements of these rotated points.

Using this, I find the indexes of all the original flattened points that are in the border. So far, so good.

Question:

How can I use this to automatically create a mesh now, though?

I know that I need the points array, of course, but also a vertices array showing how to connect the points to create a full triangle mesh. Is there any kind of algorithm I can directly apply in my case?


Solution

  • Alright, saurabheights said, I can use triangulation. Here's the code I came up with on Python, for future reference:

    import numpy as np
    from scipy.spatial import ConvexHull, Delaunay
    from open3d.geometry import TriangleMesh, PointCloud
    from open3d.utility import Vector3iVector, Vector3dVector
    
    def get_mesh(plane_normal, plane_distance, points):
        
        points = np.asarray(points)
    
        # needs signed distance
        distances = points.dot(plane_normal) + plane_distance
        points -= (distances * plane_normal[..., np.newaxis]).T
        
        rot = get_rotation_from_axis(plane_normal)
        projection = (rot @ points.T).T[:, :2]
        
        chull = ConvexHull(projection)
        borders = projection[chull.vertices]
        
        triangles = Delaunay(borders).simplices
        
        # needed to make plane visible from both sides
        triangles = np.vstack([triangles, triangles[:, ::-1]]) 
        
        mesh = TriangleMesh()
        mesh.vertices = Vector3dVector(points[chull.vertices])
        mesh.triangles = Vector3iVector(triangles)
        
        return mesh
    

    Where rotation is just a convenience function that gives the rotation matrix that aligns the axis z with the plane's normal vector.