Marc M
Marc M

Reputation: 51

open3d compute distance between mesh and point cloud

I have a CAD file (.stl) and several point clouds created by a laser scanner. now I want to calculate the difference between the CAD file and each point cloud.

first I started with Cloud Compare which helps a lot to get a basic understanding. (reduction of points, remove duplicates, create a mesh and compare distances)

In python, I was able to import the files and do some basic calculations. However, I am not able to calculate the distance.

here is my code:

import numpy as np 
import open3d as o3d

#read point cloud
dataname_pcd= "pcd.xyz"
point_cloud = np.loadtxt(input_path+dataname_pcd,skiprows=1)
#read mesh 
dataname_mesh = "cad.stl"
mesh = o3d.io.read_triangle_mesh(input_path+dataname_mesh)
print (mesh)

#calulate the distance
mD = o3d.geometry.PointCloud.compute_point_cloud_distance([point_cloud],[mesh])

#calculate the distance gives me this error: "TypeError: compute_point_cloud_distance(): incompatible function arguments. The following argument types are supported: 1. (self: open3d.cpu.pybind.geometry.PointCloud, target: open3d.cpu.pybind.geometry.PointCloud) -> open3d.cpu.pybind.utility.DoubleVector"

Questions: what pre transformations for mesh and point clouds are needed to calculate their distances? is there a recommended way to display the differences?

so far I just used the visualization line below

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])

Upvotes: 3

Views: 10072

Answers (2)

Rubens Benevides
Rubens Benevides

Reputation: 153

You need 2 point clouds for the function "compute point cloud distance()", but one of your geometries is a mesh, which is made of polygons and vertices. Just convert it to a point cloud:

pcd = o3d.geometry.PointCloud() # create a empty geometry
pcd.points = mesh.vertices      # take the vertices of your mesh

I'll illustrate how you can visualize the distances between 2 clouds, both captured on a moving robot (a Velodyne LIDAR) separeted by 1 meter in average. Consider 2 cloud before and after the registration, the distances between them should decrease, right? Here is some code:

import copy
import pandas as pd
import numpy as np
import open3d as o3d
from matplotlib import pyplot as plt

# Import 2 clouds, paint and show both
pc_1 = o3d.io.read_point_cloud("scan_0.pcd") # 18,421 points
pc_2 = o3d.io.read_point_cloud("scan_1.pcd") # 19,051 points
pc_1.paint_uniform_color([0,0,1])
pc_2.paint_uniform_color([0.5,0.5,0])
o3d.visualization.draw_geometries([pc_1,pc_2])

2 clouds before registration

# Calculate distances of pc_1 to pc_2. 
dist_pc1_pc2 = pc_1.compute_point_cloud_distance(pc_2)

# dist_pc1_pc2 is an Open3d object, we need to convert it to a numpy array to 
# acess the data
dist_pc1_pc2 = np.asarray(dist_pc1_pc2)

# We have 18,421 distances in dist_pc1_pc2, because cloud pc_1 has 18,421 pts.
# Let's make a boxplot, histogram and serie to visualize it.
# We'll use matplotlib + pandas. 
 
df = pd.DataFrame({"distances": dist_pc1_pc2}) # transform to a dataframe
# Some graphs
ax1 = df.boxplot(return_type="axes") # BOXPLOT
ax2 = df.plot(kind="hist", alpha=0.5, bins = 1000) # HISTOGRAM
ax3 = df.plot(kind="line") # SERIE
plt.show()

graphs before

# Load a previos transformation to register pc_2 on pc_1 
# I finded it with the Fast Global Registration algorithm, in Open3D 
T = np.array([[ 0.997, -0.062 ,  0.038,  1.161],
              [ 0.062,  0.9980,  0.002,  0.031],
              [-0.038,  0.001,  0.999,  0.077],
              [ 0.0,    0.0  ,  0.0   , 1.0  ]])
# Make a copy of pc_2 to preserv the original cloud
pc_2_copy = copy.deepcopy(pc_2)
# Aply the transformation T on pc_2_copy
pc_2_copy.transform(T)
o3d.visualization.draw_geometries([pc_1,pc_2_copy]) # show again

enter image description here

# Calculate distances
dist_pc1_pc2_transformed = pc_1.compute_point_cloud_distance(pc_2_copy)
dist_pc1_pc2_transformed = np.asarray(dist_pc1_pc2_transformed)
# Do as before to show diferences
df_2 = pd.DataFrame({"distances": dist_pc1_pc2_transformed})
# Some graphs (after registration)
ax1 = df_2.boxplot(return_type="axes") # BOXPLOT
ax2 = df_2.plot(kind="hist", alpha=0.5, bins = 1000) # HISTOGRAM
ax3 = df_2.plot(kind="line") # SERIE
plt.show()

As you can see, the distance between the point clouds has decreased.

graph after registration

Upvotes: 5

Mark Loyman
Mark Loyman

Reputation: 2180

compute_point_cloud_distance is a point to point distance. If that's what one wants, then Rubens Benevides's answer covers it, as well as the visualization part.

But for those who are interested in a point to plane distance (which is implied by the title of this question: "distance between mesh and point cloud":

Open3d offers point to mesh distance with it's "tensor" based api (open3d.t) with the sdf computation method (implemented with ray-casting). This is covered in the distance_queries tutorial.

def mesh_to_cloud_signed_distances(o3d_mesh: open3d.t.geometry.TriangleMesh, cloud: open3d.t.geometry.PointCloud) -> np.ndarray:
    scene = open3d.t.geometry.RaycastingScene()
    _ = scene.add_triangles(o3d_mesh)
    sdf = scene.compute_signed_distance(cloud.point.positions)
    return sdf.numpy()

Notes:

  1. The above is the signed distance, so you might want to abs it based on your needs.
  2. To convert o3d.geometry.TriangleMesh to o3d.t.geometry.TriangleMesh:
mesh = o3d.t.geometry.TriangleMesh.from_legacy(mesh)

Upvotes: 2

Related Questions