Normalizing and Feature Extraction of 3D Point Clouds

Point Cloud Normalization and Feature Extraction

Point clouds are 3D representations of objects or environments, consisting of numerous points in space. Normalizing these point clouds and extracting meaningful features are crucial steps in various applications, including object recognition, scene understanding, and 3D modeling. This presentation will guide you through the process of normalizing point clouds and extracting useful features using Python.

import numpy as np
import open3d as o3d

# Load a sample point cloud
pcd = o3d.io.read_point_cloud("sample_pointcloud.ply")
points = np.asarray(pcd.points)

# Visualize the original point cloud
o3d.visualization.draw_geometries([pcd])

Centering the Point Cloud

The first step in normalization is centering the point cloud. This involves subtracting the mean position from all points, effectively moving the center of the point cloud to the origin (0, 0, 0).

# Calculate the centroid (mean position) of the point cloud
centroid = np.mean(points, axis=0)

# Center the point cloud by subtracting the centroid
centered_points = points - centroid

# Create a new point cloud with centered points
centered_pcd = o3d.geometry.PointCloud()
centered_pcd.points = o3d.utility.Vector3dVector(centered_points)

# Visualize the centered point cloud
o3d.visualization.draw_geometries([centered_pcd])

Scaling the Point Cloud

After centering, we scale the point cloud to fit within a unit sphere. This step ensures that point clouds of different sizes are comparable and helps in maintaining consistent feature extraction.

# Calculate the maximum distance from the origin
max_distance = np.max(np.linalg.norm(centered_points, axis=1))

# Scale the point cloud to fit within a unit sphere
scaled_points = centered_points / max_distance

# Create a new point cloud with scaled points
scaled_pcd = o3d.geometry.PointCloud()
scaled_pcd.points = o3d.utility.Vector3dVector(scaled_points)

# Visualize the scaled point cloud
o3d.visualization.draw_geometries([scaled_pcd])

Voxel Downsampling

Voxel downsampling is a technique used to reduce the number of points in a point cloud while preserving its overall structure. This process helps in reducing computational complexity and memory usage for subsequent operations.

# Perform voxel downsampling
voxel_size = 0.05  # Adjust this value based on your needs
downsampled_pcd = scaled_pcd.voxel_down_sample(voxel_size)

# Visualize the downsampled point cloud
o3d.visualization.draw_geometries([downsampled_pcd])

# Print the number of points before and after downsampling
print(f"Original point cloud: {len(scaled_pcd.points)} points")
print(f"Downsampled point cloud: {len(downsampled_pcd.points)} points")

Normal Estimation

Normal vectors are important features that describe the local surface orientation at each point. Estimating normals is crucial for many point cloud processing tasks and feature extraction methods.

# Estimate normals
downsampled_pcd.estimate_normals(
    search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30)
)

# Visualize the point cloud with normals
o3d.visualization.draw_geometries([downsampled_pcd], point_show_normal=True)

# Access the computed normals
normals = np.asarray(downsampled_pcd.normals)

Feature Histogram Computation

Feature histograms are powerful descriptors that capture the distribution of geometric properties in a point cloud. One common histogram is the Fast Point Feature Histogram (FPFH), which encodes the relationships between points and their neighbors.

# Compute FPFH features
radius = 0.2
max_nn = 100
fpfh = o3d.pipelines.registration.compute_fpfh_feature(
    downsampled_pcd,
    o3d.geometry.KDTreeSearchParamHybrid(radius=radius, max_nn=max_nn)
)

# Access the computed FPFH features
fpfh_features = np.asarray(fpfh.data).T

print(f"FPFH feature shape: {fpfh_features.shape}")

Principal Component Analysis (PCA)

PCA is a technique used to reduce the dimensionality of data while preserving its most important characteristics. In point cloud processing, PCA can be used to find the principal axes of the point cloud and extract global features.

# Perform PCA on the point cloud
covariance_matrix = np.cov(scaled_points.T)
eigenvalues, eigenvectors = np.lign.eig(covariance_matrix)

# Sort eigenvalues and eigenvectors in descending order
sort_indices = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[sort_indices]
eigenvectors = eigenvectors[:, sort_indices]

# Print the eigenvalues and eigenvectors
print("Eigenvalues:", eigenvalues)
print("Eigenvectors:")
print(eigenvectors)

# Visualize the principal axes
pca_pcd = o3d.geometry.PointCloud()
pca_pcd.points = o3d.utility.Vector3dVector(scaled_points)
pca_pcd.paint_uniform_color([0.8, 0.8, 0.8])

for i in range(3):
    line_set = o3d.geometry.LineSet()
    line_set.points = o3d.utility.Vector3dVector([np.array([0, 0, 0]), eigenvectors[:, i]])
    line_set.lines = o3d.utility.Vector2iVector([[0, 1]])
    line_set.colors = o3d.utility.Vector3dVector([[1, 0, 0]])
    o3d.visualization.draw_geometries([pca_pcd, line_set])

Real-Life Example: 3D Object Recognition

In 3D object recognition, normalized point clouds and extracted features are used to identify and classify objects in a scene. For example, in a robotics application, a robot might need to recognize and grasp objects on a table.

import numpy as np
import open3d as o3d
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split

# Load and preprocess multiple object point clouds
def load_and_preprocess(file_path):
    pcd = o3d.io.read_point_cloud(file_path)
    points = np.asarray(pcd.points)
    centered_points = points - np.mean(points, axis=0)
    scaled_points = centered_points / np.max(np.linalg.norm(centered_points, axis=1))
    return scaled_points

# Load example objects (you would have multiple objects in a real scenario)
mug = load_and_preprocess("mug.ply")
plate = load_and_preprocess("plate.ply")
fork = load_and_preprocess("fork.ply")

# Extract features (using FPFH as an example)
def extract_features(pcd):
    pcd_o3d = o3d.geometry.PointCloud()
    pcd_o3d.points = o3d.utility.Vector3dVector(pcd)
    pcd_o3d.estimate_normals()
    fpfh = o3d.pipelines.registration.compute_fpfh_feature(
        pcd_o3d,
        o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30)
    )
    return np.asarray(fpfh.data).T.mean(axis=0)  # Use mean FPFH as a simple global descriptor

# Create a dataset
X = np.vstack([extract_features(mug), extract_features(plate), extract_features(fork)])
y = np.array(['mug', 'plate', 'fork'])

# Train a simple classifier
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
clf = RandomForestClassifier(n_estimators=100, random_state=42)
clf.fit(X_train, y_train)

# Predict on test data
accuracy = clf.score(X_test, y_test)
print(f"Model accuracy: {accuracy:.2f}")

Real-Life Example: 3D Scene Registration

Point cloud registration is the process of aligning multiple point clouds to create a complete 3D model. This is commonly used in 3D scanning, robotics, and augmented reality applications.

import numpy as np
import open3d as o3d

# Load two partial scans of the same object
source = o3d.io.read_point_cloud("partial_scan_1.ply")
target = o3d.io.read_point_cloud("partial_scan_2.ply")

# Preprocess point clouds
def preprocess_point_cloud(pcd, voxel_size):
    pcd_down = pcd.voxel_down_sample(voxel_size)
    pcd_down.estimate_normals(
        o3d.geometry.KDTreeSearchParamHybrid(radius=voxel_size * 2, max_nn=30))
    pcd_fpfh = o3d.pipelines.registration.compute_fpfh_feature(
        pcd_down,
        o3d.geometry.KDTreeSearchParamHybrid(radius=voxel_size * 5, max_nn=100))
    return pcd_down, pcd_fpfh

# Preprocess both point clouds
voxel_size = 0.05
source_down, source_fpfh = preprocess_point_cloud(source, voxel_size)
target_down, target_fpfh = preprocess_point_cloud(target, voxel_size)

# Perform global registration
result_ransac = o3d.pipelines.registration.registration_ransac_based_on_feature_matching(
    source_down, target_down, source_fpfh, target_fpfh, True,
    0.075,  # maximum_correspondence_distance
    o3d.pipelines.registration.TransformationEstimationPointToPoint(False),
    3,  # ransac_n
    [o3d.pipelines.registration.CorrespondenceCheckerBasedOnEdgeLength(0.9),
     o3d.pipelines.registration.CorrespondenceCheckerBasedOnDistance(0.075)],
    o3d.pipelines.registration.RANSACConvergenceCriteria(4000000, 500))

# Refine the alignment using ICP
result_icp = o3d.pipelines.registration.registration_icp(
    source_down, target_down, 0.02, result_ransac.transformation,
    o3d.pipelines.registration.TransformationEstimationPointToPoint())

# Transform the source point cloud
source.transform(result_icp.transformation)

# Visualize the result
o3d.visualization.draw_geometries([source, target])

Intensity and Color Features

In addition to geometric features, point clouds often include intensity or color information. These attributes can be valuable for feature extraction and object recognition tasks.

import numpy as np
import open3d as o3d
import matplotlib.pyplot as plt

# Load a colored point cloud
pcd = o3d.io.read_point_cloud("colored_pointcloud.ply")

# Extract color information
colors = np.asarray(pcd.colors)

# Compute color histograms
def compute_color_histogram(colors, bins=32):
    hist_r, _ = np.histogram(colors[:, 0], bins=bins, range=(0, 1))
    hist_g, _ = np.histogram(colors[:, 1], bins=bins, range=(0, 1))
    hist_b, _ = np.histogram(colors[:, 2], bins=bins, range=(0, 1))
    return hist_r, hist_g, hist_b

hist_r, hist_g, hist_b = compute_color_histogram(colors)

# Visualize color histograms
plt.figure(figsize=(12, 4))
plt.subplot(131)
plt.bar(range(len(hist_r)), hist_r, color='red', alpha=0.7)
plt.title('Red Channel Histogram')
plt.subplot(132)
plt.bar(range(len(hist_g)), hist_g, color='green', alpha=0.7)
plt.title('Green Channel Histogram')
plt.subplot(133)
plt.bar(range(len(hist_b)), hist_b, color='blue', alpha=0.7)
plt.title('Blue Channel Histogram')
plt.tight_layout()
plt.show()

# Visualize the colored point cloud
o3d.visualization.draw_geometries([pcd])

Local Surface Features

Local surface features capture the geometric properties of a point's neighborhood. These features are useful for describing the shape and curvature of the surface around each point.

import numpy as np
import open3d as o3d

# Load a point cloud
pcd = o3d.io.read_point_cloud("sample_pointcloud.ply")

# Estimate normals
pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30))

# Compute local surface features
def compute_local_features(pcd, radius, max_nn):
    points = np.asarray(pcd.points)
    normals = np.asarray(pcd.normals)
    tree = o3d.geometry.KDTreeFlann(pcd)
    
    features = []
    for i in range(len(points)):
        [k, idx, _] = tree.search_hybrid_vector_3d(points[i], radius, max_nn)
        if k < 3:
            features.append([0, 0, 0])  # Not enough neighbors
            continue
        
        neighbor_points = points[idx[1:]]  # Exclude the point itself
        neighbor_normals = normals[idx[1:]]
        
        # Compute PCA of the neighborhood
        covariance_matrix = np.cov(neighbor_points.T)
        eigenvalues, _ = np.linalg.eig(covariance_matrix)
        eigenvalues = np.sort(eigenvalues)[::-1]
        
        # Compute shape features
        linearity = (eigenvalues[0] - eigenvalues[1]) / eigenvalues[0]
        planarity = (eigenvalues[1] - eigenvalues[2]) / eigenvalues[0]
        sphericity = eigenvalues[2] / eigenvalues[0]
        
        features.append([linearity, planarity, sphericity])
    
    return np.array(features)

local_features = compute_local_features(pcd, radius=0.1, max_nn=30)

# Visualize the point cloud colored by a feature (e.g., planarity)
planarity = local_features[:, 1]
planarity_normalized = (planarity - np.min(planarity)) / (np.max(planarity) - np.min(planarity))
colors = np.zeros((len(planarity), 3))
colors[:, 0] = planarity_normalized  # Red channel

pcd.colors = o3d.utility.Vector3dVector(colors)
o3d.visualization.draw_geometries([pcd])

Global Descriptors

Global descriptors provide a compact representation of the entire point cloud, useful for tasks such as object classification and retrieval. One such descriptor is the Viewpoint Feature Histogram (VFH), which captures both local geometry and global viewpoint information.

import numpy as np
import open3d as o3d

def compute_vfh(pcd):
    # Estimate normals if not already present
    if not pcd.has_normals():
        pcd.estimate_normals()
    
    points = np.asarray(pcd.points)
    normals = np.asarray(pcd.normals)
    
    # Compute centroid and central direction
    centroid = np.mean(points, axis=0)
    central_normal = np.mean(normals, axis=0)
    central_normal /= np.linalg.norm(central_normal)
    
    # Initialize histogram bins
    num_bins = 45
    vfh_signature = np.zeros(3 * num_bins + 128)
    
    for i in range(len(points)):
        point = points[i]
        normal = normals[i]
        
        # Compute angle between central direction and normal
        cos_angle = np.dot(central_normal, normal)
        angle = np.arccos(np.clip(cos_angle, -1.0, 1.0))
        
        # Compute distances and angles
        distance = np.linalg.norm(point - centroid)
        azimuth = np.arctan2(point[1] - centroid[1], point[0] - centroid[0])
        elevation = np.arcsin((point[2] - centroid[2]) / distance)
        
        # Update histogram bins
        bin_distance = int((distance / np.max(distance)) * (num_bins - 1))
        bin_azimuth = int((azimuth + np.pi) / (2 * np.pi) * (num_bins - 1))
        bin_elevation = int((elevation + np.pi/2) / np.pi * (num_bins - 1))
        
        vfh_signature[bin_distance] += 1
        vfh_signature[num_bins + bin_azimuth] += 1
        vfh_signature[2 * num_bins + bin_elevation] += 1
    
    # Normalize the histogram
    vfh_signature /= len(points)
    
    return vfh_signature

# Load and compute VFH for a sample point cloud
pcd = o3d.io.read_point_cloud("sample_pointcloud.ply")
vfh = compute_vfh(pcd)

print("VFH descriptor shape:", vfh.shape)

Feature Matching and Correspondence

Feature matching is crucial for tasks like object recognition and point cloud registration. This slide demonstrates how to match features between two point clouds using the Fast Library for Approximate Nearest Neighbors (FLANN).

import numpy as np
import open3d as o3d
from sklearn.neighbors import NearestNeighbors

def match_features(features1, features2, ratio_threshold=0.8):
    # Use FLANN for fast approximate nearest neighbor search
    nn = NearestNeighbors(n_neighbors=2, algorithm='kd_tree').fit(features2)
    distances, indices = nn.kneighbors(features1)
    
    # Apply ratio test
    good_matches = []
    for i, (d1, d2) in enumerate(distances):
        if d1 < ratio_threshold * d2:
            good_matches.append((i, indices[i][0]))
    
    return good_matches

# Load two point clouds
source = o3d.io.read_point_cloud("source.ply")
target = o3d.io.read_point_cloud("target.ply")

# Compute FPFH features for both point clouds
radius = 0.1
source_fpfh = o3d.pipelines.registration.compute_fpfh_feature(
    source, o3d.geometry.KDTreeSearchParamHybrid(radius=radius, max_nn=100))
target_fpfh = o3d.pipelines.registration.compute_fpfh_feature(
    target, o3d.geometry.KDTreeSearchParamHybrid(radius=radius, max_nn=100))

# Match features
matches = match_features(source_fpfh.data.T, target_fpfh.data.T)

print(f"Number of matches: {len(matches)}")

# Visualize matched points
source_points = np.asarray(source.points)
target_points = np.asarray(target.points)
matched_source = o3d.geometry.PointCloud()
matched_target = o3d.geometry.PointCloud()

for source_idx, target_idx in matches:
    matched_source.points.append(source_points[source_idx])
    matched_target.points.append(target_points[target_idx])

matched_source.paint_uniform_color([1, 0, 0])  # Red for source
matched_target.paint_uniform_color([0, 1, 0])  # Green for target

o3d.visualization.draw_geometries([matched_source, matched_target])

Point Cloud Segmentation

Segmentation is the process of dividing a point cloud into meaningful segments or clusters. This is useful for identifying distinct objects or regions within a scene.

import numpy as np
import open3d as o3d

def segment_plane(pcd, distance_threshold=0.01, ransac_n=3, num_iterations=1000):
    plane_model, inliers = pcd.segment_plane(distance_threshold, ransac_n, num_iterations)
    return plane_model, inliers

def cluster_points(pcd, eps=0.05, min_points=10):
    labels = np.array(pcd.cluster_dbscan(eps=eps, min_points=min_points))
    max_label = labels.max()
    return labels, max_label

# Load a point cloud
pcd = o3d.io.read_point_cloud("scene.ply")

# Remove the dominant plane (e.g., floor or table surface)
plane_model, inliers = segment_plane(pcd)
non_plane_cloud = pcd.select_by_index(inliers, invert=True)

# Cluster the remaining points
labels, max_label = cluster_points(non_plane_cloud)

print(f"Number of clusters: {max_label + 1}")

# Visualize the clusters
colors = plt.get_cmap("tab20")(labels / (max_label if max_label > 0 else 1))
non_plane_cloud.colors = o3d.utility.Vector3dVector(colors[:, :3])
o3d.visualization.draw_geometries([non_plane_cloud])

Additional Resources

For further exploration of point cloud processing and feature extraction techniques, consider the following resources:

1. ArXiv paper: "PointNet: Deep Learning on Point Sets for 3D Classification and Segmentation" by Qi et al. (2017) URL: https://arxiv.org/abs/1612.00593
2. ArXiv paper: "PointNet++: Deep Hierarchical Feature Learning on Point Sets in a Metric Space" by Qi et al. (2017) URL: https://arxiv.org/abs/1706.02413
3. ArXiv paper: "FPFH: Fast Point Feature Histograms (FPFH) for 3D Registration" by Rusu et al. (2009) URL: https://arxiv.org/abs/0905.2693

These papers provide in-depth discussions on advanced techniques for point cloud processing and feature extraction, building upon the concepts covered in this presentation.