Persistence of tensor su nuvole di punti

Questo post testa il metodo presentato nell'articolo Line Structure Extraction from LiDAR Point Cloud Based on the Persistence of Tensor Feature Wang, X.; Lyu, H.; He, W.; Chen, Q. Line Structure Extraction from LiDAR Point Cloud Based on the Persistence of Tensor Feature. Appl. Sci. 2022, 12, 9190. https://doi.org/10.3390/app12189190

Il metodo permette di estrarre lineazioni, piani da una nuvola di punti (in pratica una generalizzazione dei RANSAC ma in questo caso non si riesce a discriminare tra piani a differente orientazione)

L'algoritmo e' stato testato sulla nuvola di punti derivante da questo modello 

L'unico parametro del metodo e' il radius a riga 55

[k, idx, _] = pcd_tree.search_knn_vector_3d(pcd.points[i], 100)

Esempio su dati di una scala (dataset compreso nei file di esempio)

import open3d as o3d

import numpy as np

def extract_tensor_features(pcd_path, neighborhood_radius=0.05):

"""

Extracts linear, planar, and spherical features from a point cloud using

the Tensor Feature method.

Args:

pcd_path (str): The path to the input.ply file.

neighborhood_radius (float): The radius for neighborhood search.

"""

# Step 1: Load the point cloud

try:

pcd = o3d.io.read_point_cloud(pcd_path)

if len(pcd.points) == 0:

print(f"Error: The point cloud file at {pcd_path} is empty.")

return

print(f"Loaded point cloud with {len(pcd.points)} points from {pcd_path}")

except Exception as e:

print(f"Error loading point cloud: {e}")

return

# Step 2: Estimate normals, which is necessary for tensor encoding

pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamRadius(radius=neighborhood_radius))

# CRITICAL CHECK: Ensure normals were successfully created

if not pcd.has_normals():

print("\nNormal estimation failed. This could be due to:")

print("- The point cloud is empty or contains only a few points.")

print("- The 'neighborhood_radius' is too small for the point cloud's density.")

print("Please check your input file or try increasing the 'neighborhood_radius' parameter.")

return

pcd.orient_normals_consistent_tangent_plane(k=10)

print("Normals estimated and oriented.")

# Create a KDTree for efficient neighborhood searches

pcd_tree = o3d.geometry.KDTreeFlann(pcd)

# Prepare arrays to store classification results

points = np.asarray(pcd.points)

num_points = len(points)

classification_labels = np.zeros(num_points, dtype=int)

colors = np.zeros((num_points, 3))

# Define colors for visualization

LINEAR_COLOR = np.array((1.0, 0.0, 0.0)) # Red for linear features (e.g., edges, lines)

PLANAR_COLOR = np.array((0.0, 1.0, 0.0)) # Green for planar features (e.g., surfaces, walls)

SPHERICAL_COLOR = np.array((0.0, 0.0, 1.0)) # Blue for spherical features (e.g., noise, isolated points)

# Step 3 & 4: Compute tensor, decompose, and classify

for i in range(num_points):

# Find neighbors within the specified radius

[k, idx, _] = pcd_tree.search_knn_vector_3d(pcd.points[i], 100)

# If not enough neighbors, skip the point

if len(idx) < 10:

continue

neighbors = points[idx, :]

# Compute the covariance tensor

center = np.mean(neighbors, axis=0)

cov_tensor = np.cov((neighbors - center).T)

# Eigen-decomposition

eigenvalues, _ = np.linalg.eig(cov_tensor)

eigenvalues = np.sort(eigenvalues)[::-1] # Sort descending

lambda1, lambda2, lambda3 = eigenvalues

# Calculate shape features based on eigenvalues

# Linear feature: one large eigenvalue

linearity = (lambda1 - lambda2) / lambda1 if lambda1 > 0 else 0

# Planar feature: two large eigenvalues

planarity = (lambda2 - lambda3) / lambda1 if lambda1 > 0 else 0

# Spherical feature: all eigenvalues are similar

sphericity = lambda3 / lambda1 if lambda1 > 0 else 0

# Classify the point based on the dominant feature

if linearity > 0.8: # Heuristic threshold for linearity

classification_labels[i] = 1

colors[i] = LINEAR_COLOR

elif planarity > 0.8: # Heuristic threshold for planarity

classification_labels[i] = 2

colors[i] = PLANAR_COLOR

else:

classification_labels[i] = 3

colors[i] = SPHERICAL_COLOR

# Step 5: Visualize the results

pcd.colors = o3d.utility.Vector3dVector(colors)

o3d.visualization.draw_geometries([pcd], window_name="Point Cloud Feature Extraction")

# Print a summary of the results

linear_points = np.sum(classification_labels == 1)

planar_points = np.sum(classification_labels == 2)

spherical_points = np.sum(classification_labels == 3)

print("\n--- Classification Summary ---")

print(f"Linear Features: {linear_points} points ({linear_points/num_points*100:.2f}%)")

print(f"Planar Features: {planar_points} points ({planar_points/num_points*100:.2f}%)")

print(f"Spherical Features: {spherical_points} points ({spherical_points/num_points*100:.2f}%)")

if __name__ == "__main__":

file_path = "./output_depth/pointcloud.ply"

extract_tensor_features(file_path)