Python Open3D geometry foundation chapter point cloud operation

Python Open3D geometry Fundamentals (I) point cloud operation

Pages referenced in this article:

This article mainly introduces the basic usage of point cloud in Open3D:

Visual point cloud( Visualize point cloud)

How to read and visualize the point cloud:

read_point_cloud reads the point cloud from the file. It will parse the point cloud file according to the suffix. The specific type needs to be referred to File IO
draw_ Geometry visualizes the point cloud, which can be viewed from different angles with the mouse
The visual API can view this page

This looks like a dense surface, but in fact, the point cloud is rendered as surface primitives. This GUI interface supports different key functions

  • Press the "-" key to reduce the size of surface elements
  • Press the "H" key to print out the list of all keyboard commands of the GUI
  • User defined visualization can refer to Custom visualization page
  • In macOS, the GUI may not receive keyboard events. You can try python instead of python

Example code:

print("Load a ply point cloud, print it, and render it") # Load the PLY point cloud, print it and display it
ply_point_cloud = o3d.data.PLYPointCloud() # Get the point cloud data. If the download is unsuccessful, you can delete this line of code
pcd = o3d.io.read_point_cloud(ply_point_cloud.path) # ply_point_cloud.path, which can be changed to the address of your own point cloud file
print(pcd)
print(np.asarray(pcd.points))
o3d.visualization.draw_geometries([pcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024]
""" 
Load a ply point cloud, print it, and render it
[Open3D INFO] Downloading https://github.com/isl-org/open3d_downloads/releases/download/20220201-data/fragment.ply
[Open3D INFO] Downloaded to /home/runner/open3d_data/download/PLYPointCloud/fragment.ply
PointCloud with 196133 points.
[[0.65234375 0.84686458 2.37890625]
 [0.65234375 0.83984375 2.38430572]
 [0.66737998 0.83984375 2.37890625]
 ...
 [2.00839925 2.39453125 1.88671875]
 [2.00390625 2.39488506 1.88671875]
 [2.00390625 2.39453125 1.88793314]]
"""                                 

Voxel down sampling( Voxel downsampling)

Voxel downsampling uses a regular voxel mesh to create a uniform downsampling point cloud from the input point cloud. It is often used as a preprocessing step for many point cloud processing tasks. The algorithm is divided into two steps:

  1. Points are bucket into voxels
  2. Each occupied voxel generates exactly one point by averaging all points inside
print("Downsample the point cloud with a voxel of 0.05")
downpcd = pcd.voxel_down_sample(voxel_size=0.05)
o3d.visualization.draw_geometries([downpcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

Estimate vertex normal( Vertex normal estimation)

Another basic point cloud operation is the normal estimation of points. Press the "N" key on the operation interface to see the normal vector of the point. "-" and "+" can control the length of the normal vector

print("Recompute the normal of the downsampled point cloud")
downpcd.estimate_normals(
    search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30))
o3d.visualization.draw_geometries([downpcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024],
                                  point_show_normal=True)

  • estimate_normals can calculate the normal vector of each point. This function can find nearby points and calculate the Principal Axis of nearby points using Covariance Analysis.
  • This function takes an instance of the KDTreeSearchParamHybrid class as a parameter. Two key parameters radius = 0.1 and max_nn = 30 specifies the search radius and the maximum nearest neighbor. It has a search radius of 10 cm and only considers 30 neighbors at most to save computing time.

The analysis of covariance algorithm generates two opposite directions as candidate normal directions.
Without knowing the overall structure of the geometry, both directions may be right. This is the so-called normal direction problem.
Open3D attempts to orient the normals to align them with the original normals, if any. Otherwise, Open3D will make random guesses.
If you want to consider the direction, you can call the advanced direction function, such as orient_normals_to_align_with_direction and orientation_ normals_ towards_ camera_ location

View the estimated vertex normal( Access estimated vertex normal)

  • The estimated normal vector can be retrieved through the member variable normals of downpcd

    print("Print a normal vector of the 0th point")
    print(downpcd.normals[0])
    
    # Print a normal vector of the 0th point
    # [-0.27566603 -0.89197839 -0.35830543]
    
    
  • If you want to view other related variables, you can use help(downpcd)

  • The normal vector can be passed through NP Convert asarray to Numpy array

    • Sample code:
      print("Print the normal vectors of the first 10 points")
      print(np.asarray(downpcd.normals)[:10, :])
      
      """
      Print the normal vectors of the first 10 points
      [[-0.27566603 -0.89197839 -0.35830543]
       [-0.04230441 -0.99410664 -0.09981149]
       [-0.00399871 -0.99965423 -0.02598917]
       [-0.93768261 -0.07378998  0.3395679 ]
       [-0.43476205 -0.62438493 -0.64894177]
       [-0.09739809 -0.9928602  -0.06886388]
       [-0.27498453 -0.67317361 -0.68645524]
       [-0.11728718 -0.95516445 -0.27185399]
       [-0.00816546 -0.99965616 -0.02491762]
       [-0.11067463 -0.99205156 -0.05987351]]
      """
      
      
  • Examples of functions related to Numpy can be seen here link

Crop point cloud( Crop point cloud)

print("Load a polygon volume and use it to crop the original point cloud")
demo_crop_data = o3d.data.DemoCropPointCloud()
pcd = o3d.io.read_point_cloud(demo_crop_data.point_cloud_path)
vol = o3d.visualization.read_selection_polygon_volume(demo_crop_data.cropped_json_path)
chair = vol.crop_point_cloud(pcd)
o3d.visualization.draw_geometries([chair],
                                  zoom=0.7,
                                  front=[0.5439, -0.2333, -0.8060],
                                  lookat=[2.4615, 2.1331, 1.338],
                                  up=[-0.1781, -0.9708, 0.1608])

read_selection_polygon_volume reads the json file of the specified polygon selection area.
vol.crop_point_cloud(pcd) filters other points, leaving only chairs.

Point cloud coloring( Paint point cloud)

print("Paint chair")
chair.paint_uniform_color([1, 0.706, 0])
o3d.visualization.draw_geometries([chair],
                                  zoom=0.7,
                                  front=[0.5439, -0.2333, -0.8060],
                                  lookat=[2.4615, 2.1331, 1.338],
                                  up=[-0.1781, -0.9708, 0.1608])

paint_uniform_color uses a uniform color to color all points. The color on is in RGB space, ranging from 0 to 1.

Point cloud distance( Point cloud distance)

  • Open3D provides compute_ point_ cloud_ The distance method is used to calculate the distance between the target point cloud and the source point cloud. That is, it can calculate the distance between each point of the source point cloud and the nearest point of the target point cloud.
  • In the following example, we use this function to calculate the difference between two point clouds. It should be noted that this method can also calculate the Chamfer distance between two point clouds
# Load data
demo_crop_data = o3d.data.DemoCropPointCloud()
pcd = o3d.io.read_point_cloud(demo_crop_data.point_cloud_path)
vol = o3d.visualization.read_selection_polygon_volume(demo_crop_data.cropped_json_path)
chair = vol.crop_point_cloud(pcd)

dists = pcd.compute_point_cloud_distance(chair)
dists = np.asarray(dists)
ind = np.where(dists > 0.01)[0]
pcd_without_chair = pcd.select_by_index(ind)
o3d.visualization.draw_geometries([pcd_without_chair],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

Bounding box( Bounding volumes)

PointCloud geometry types have bounding boxes like all other geometry types in Open3D.
At present, the Open3D implementation of AxisAlignedBoundingBox and OrientedBoundingBox can also be used to clip geometry.

aabb = chair.get_axis_aligned_bounding_box()
aabb.color = (1, 0, 0)
obb = chair.get_oriented_bounding_box()
obb.color = (0, 1, 0)
o3d.visualization.draw_geometries([chair, aabb, obb],
                                  zoom=0.7,
                                  front=[0.5439, -0.2333, -0.8060],
                                  lookat=[2.4615, 2.1331, 1.338],
                                  up=[-0.1781, -0.9708, 0.1608])

Convex hull( Convex hull)

  • The convex hull of a point cloud is the smallest convex set containing all points. There is compute in Open3D_ convex_ Hull method is used to calculate the convex hull of point cloud Qhull Implemented.
  • In the following example code, first sample a point cloud from the mesh, and then calculate the convex hull returned as a triangular mesh. Then use the red LineSet to visualize the convex hull.
bunny = o3d.data.BunnyMesh()
mesh = o3d.io.read_triangle_mesh(bunny.path)
mesh.compute_vertex_normals()

pcl = mesh.sample_points_poisson_disk(number_of_points=2000)
hull, _ = pcl.compute_convex_hull()
hull_ls = o3d.geometry.LineSet.create_from_triangle_mesh(hull)
hull_ls.paint_uniform_color((1, 0, 0))
o3d.visualization.draw_geometries([pcl, hull_ls])

DBSCAN clustering algorithm( DBSCAN clustering)

  • Given a point cloud, such as a depth sensor, we want to group local point clouds together. To this end, we can use clustering algorithm. Open3D implements DBSCAN [Ester1996]
    This is based on density clustering algorithm.
  • This algorithm is through cluster_ Implemented by DBSCAN, this method requires two parameters eps and min_points ;
    • eps defines the distance between neighbors in the cluster,
    • min_points defines the minimum number of points needed to form a cluster.
    • This function will return labels, where tag - 1 represents noise.
  • Example code:

    ply_point_cloud = o3d.data.PLYPointCloud()
    pcd = o3d.io.read_point_cloud(ply_point_cloud.path)
    
    with o3d.utility.VerbosityContextManager(
            o3d.utility.VerbosityLevel.Debug) as cm:
        labels = np.array(
            pcd.cluster_dbscan(eps=0.02, min_points=10, print_progress=True))
    
    max_label = labels.max()
    print(f"point cloud has {max_label + 1} clusters")
    colors = plt.get_cmap("tab20")(labels / (max_label if max_label > 0 else 1))
    colors[labels < 0] = 0
    pcd.colors = o3d.utility.Vector3dVector(colors[:, :3])
    o3d.visualization.draw_geometries([pcd],
                                      zoom=0.455,
                                      front=[-0.4999, -0.1659, -0.8499],
                                      lookat=[2.1813, 2.0619, 2.0999],
                                      up=[0.1204, -0.9852, 0.1215])
    

  • This algorithm will calculate the radius in advance ε All my neighbors. If you choose a large ε It will require a lot of memory.

Plane segmentation( Plane segmentation)

  • Open3D also supports the use of RANSAC to segment geometric primitives from point clouds.
  • We can use segment_plane finds the most supported plane in the point cloud. This method has three parameters, distance_threshold ,ransac_n ,num_iterations :
    • distance_threshold defines the maximum distance from a point to an interior point in the estimation plane. (the maximum distance a point can have to an estimated plane to be considered an inlier)
    • ransac_n defines the number of points randomly sampled on an estimation plane.
    • num_iterations defines the sampling and verification frequency of the random plane.
    • The function returns the plane $(a,b,c,d) also Just yes each one individual spot That is, every point That is, every point (x,y,z) stay stay On ax+by+cz+d=0 $
    • . The function further returns a list of indices of the inlier points.
  • code

    pcd_point_cloud = o3d.data.PCDPointCloud()
    pcd = o3d.io.read_point_cloud(pcd_point_cloud.path)
    
    plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,
                                             ransac_n=3,
                                             num_iterations=1000)
    [a, b, c, d] = plane_model
    print(f"Plane equation: {a:.2f}x + {b:.2f}y + {c:.2f}z + {d:.2f} = 0")
    
    inlier_cloud = pcd.select_by_index(inliers)
    inlier_cloud.paint_uniform_color([1.0, 0, 0])
    outlier_cloud = pcd.select_by_index(inliers, invert=True)
    o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud],
                                      zoom=0.8,
                                      front=[-0.4999, -0.1659, -0.8499],
                                      lookat=[2.1813, 2.0619, 2.0999],
                                      up=[0.1204, -0.9852, 0.1215])
    # Plane equation: -0.06x + -0.10y + 0.99z + -1.06 = 0
    

Hidden point removal( Hidden point removal)

  • Suppose you want to render a point cloud from a given viewpoint, but the points in the background will leak out in the foreground because they will not be obscured by other points. Therefore, we can use the hidden point removal algorithm.
  • In Open3D [Katz2007] This algorithm is based on the visibility of the approximate point cloud of a given view, and does not need surface reconstruction or normal estimation.
  • code

    print("Convert mesh to a point cloud and estimate dimensions")
    armadillo = o3d.data.ArmadilloMesh()
    mesh = o3d.io.read_triangle_mesh(armadillo.path)
    mesh.compute_vertex_normals()
    
    pcd = mesh.sample_points_poisson_disk(5000)
    diameter = np.linalg.norm(
        np.asarray(pcd.get_max_bound()) - np.asarray(pcd.get_min_bound()))
    o3d.visualization.draw_geometries([pcd])
    

print("Define parameters used for hidden_point_removal")
camera = [0, 0, diameter]
radius = diameter * 100

print("Get all points that are visible from given view point")
_, pt_map = pcd.hidden_point_removal(camera, radius)

print("Visualize result")
pcd = pcd.select_by_index(pt_map)
o3d.visualization.draw_geometries([pcd])

Tags: Python Computer Vision open3d

Posted by N-Bomb(Nerd) on Fri, 06 May 2022 09:52:40 +0300