3D reconstruction from multiple depth maps [closed]

Question

I'm working on reconstructing a 3D object from multiple depth maps taken at different angles. The object is rotating on a turntable, and the hyperspectral camera is stationary. The depth data is stored in a 570x570 2D array for each photo, and I'm using rotation matrices to align the points.

I'm running into issues with overlap and incorrect rotations when trying to merge the depth maps. The code to extract the depth information for every pixel from the camera works, but I'm struggling with properly aligning and merging the depth maps.

How can I fix the rotation and merging of depth maps to get a correct 3D model? Note: I'm working in Python. Should I use something else than rotation matrices?

This is the point where I've got now. I'm able to convert the datacubes (that the hyperspectral camera generates after taking a picture), into 2D depth maps. The next step in my code will be to merge the depth maps into a 3D-shape, but that's where I'm stuck. Code:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Settings
num_pictures = 10  # This is a choice of the user.
rotation_angle = 360 / num_pictures  # Angle difference per picture
image_size = 570  # 2D-array is 570x570
camera_distance = 162  # Depth distance from camera to rotation-axis / center of object

# Simulate a object_distance_map (dummy data for tests)
object_distance_map = np.random.uniform(100, 300, (image_size, image_size))
object_distance_map -= camera_distance  # Correctie toepassen

# Assenstelsel define
x_range = np.linspace(-image_size // 2, image_size // 2, image_size)
y_range = np.linspace(0, image_size, image_size)

X, Y = np.meshgrid(x_range, y_range)
Z = object_distance_map  # Depth information

# 3D-puntenwolk initiëren
points_3d = np.column_stack((X.flatten(), Y.flatten(), Z.flatten()))

# Function to create a rotation matrix around the y-axis
def rotation_matrix_y(theta):
    theta_rad = np.radians(theta)
    return np.array([
        [np.cos(theta_rad), 0, np.sin(theta_rad)],
        [0, 1, 0],
        [-np.sin(theta_rad), 0, np.cos(theta_rad)]
    ])

# Transform 3D points per rotation angle
all_points = []
for i in range(num_pictures):
    theta = i * rotation_angle
    R = rotation_matrix_y(theta)
    rotated_points = points_3d @ R.T  
    all_points.append(rotated_points)

# Combine all points
all_points = np.vstack(all_points)

# Visualisation
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(all_points[:, 0], all_points[:, 1], all_points[:, 2], s=0.5, alpha=0.5)
ax.set_xlabel("X-axis")
ax.set_ylabel("Y-axis") # y-axis points upwarts
ax.set_zlabel("Z-axis") # z-axis should be depth-axis
ax.set_title("3D-reconstruction of object")
plt.show()

Note: object_distance_map is the 570x570 array, which contains the depth-value between camera and a certain pixel of 1 image. If we take 10 pictures, we have 10 maps as input.

Answer 1

You can estimate the surface using the alpha shape method in open3d. The example below is for a point cloud distribution that is similar to yours, but using fewer points and normalised coordinates.

The raw 3D point cloud used for testing, coloured by camera angle:

Surface reconstructed using the alpha shape method:

Also visualised using trimesh (more responsive interactivity compared to matplotlib):

The main part of the code takes your entire point cloud all_points and uses it to instantiate an open3d PointCloud object, o3d_point_cloud:

# ---- To open3d point cloud object ----
o3d_point_cloud = open3d.geometry.PointCloud()
o3d_point_cloud.points = open3d.utility.Vector3dVector(all_points)

The surface normals are estimated and aligned before reconstructing the surface. Tune alpha= to get a tighter fit (small alpha) or a more regularised/smoother fit (larger alpha).

# Estimate normals prior to surface recon
o3d_point_cloud.estimate_normals()
o3d_point_cloud.orient_normals_consistent_tangent_plane(k=100)

# Surface reconstruction using alpha shape
mesh = open3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(
    o3d_point_cloud, alpha=0.07
)
mesh.compute_vertex_normals()

I finally smooth the resulting mesh:

mesh = mesh.filter_smooth_simple(number_of_iterations=2)

Reproducible example

Imports and data for testing:

import numpy as np

#Static pyplots
%matplotlib notebook

#Interactive pyplots
# %matplotlib widget

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import open3d


# ---- Image settings ----
num_pictures = 7

#Angle difference per picture
rotation_angle = 360 / num_pictures

#2D-array shaped (image_size, image_size)
image_size = 570 // 10

#Depth distance from camera to rotation-axis / center of object
camera_distance = 0.5

# Define axes
x_range = np.linspace(0, 1, image_size)
y_range = np.linspace(0, 1, image_size)
X, Y = np.meshgrid(x_range, y_range)


# ---- Simulate a object_distance_map for testing ----
object_distance_map = np.random.uniform(low=0.1, high=0.3, size=(image_size, image_size))

#Correction
object_distance_map -= camera_distance

# Depth information
Z = object_distance_map

# Initiate 3D dot cloud
points_3d = np.column_stack((X.flatten(), Y.flatten(), Z.flatten()))

# Function to create a rotation matrix around the y-axis
# Assuming format: R dot [3 x n] coordinates array
def rotation_matrix_y(theta_deg):
    c, s = [trig_fn(np.deg2rad(theta_deg)) for trig_fn in (np.cos, np.sin)]

    return np.array([
        [c, 0, s],
        [0, 1, 0],
        [-s, 0, c]
    ])

# Transform 3D points per rotation angle
points_per_angle = []
for i in range(num_pictures):
    theta_deg = i * rotation_angle
    R = rotation_matrix_y(theta_deg)
    rotated_points = (R @ points_3d.T).T
    points_per_angle.append(rotated_points)
all_points = np.vstack(points_per_angle)

Convert to PointCloud object for surface reconstruction:

# ---- To open3d point cloud object ----
o3d_point_cloud = open3d.geometry.PointCloud()
o3d_point_cloud.points = open3d.utility.Vector3dVector(all_points)

# Estimate normals prior to surface recon
o3d_point_cloud.estimate_normals()
o3d_point_cloud.orient_normals_consistent_tangent_plane(k=100)

# Surface reconstruction using alpha shape
mesh = open3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(
    o3d_point_cloud, alpha=0.07
)
mesh.compute_vertex_normals()

#Smooth the mesh
mesh = mesh.filter_smooth_simple(number_of_iterations=2)

Visualise point cloud and surface using matplotlib:

# ---- Visualisation ----
plt.close('all')
fig = plt.figure(figsize=(10, 10), num=0, clear=True)
ax = fig.add_subplot(projection='3d')
ax.view_init(azim=95, elev=6, roll=0)

# Raw point cloud (colouring by angle for visualisation)
for i, points in enumerate(points_per_angle):
    ax.scatter(*points.T, s=1, c=plt.get_cmap('copper')(i / num_pictures))

#Estimated surface
vertices, triangles = np.asarray(mesh.vertices), np.asarray(mesh.triangles)
points_per_tri = [vertices[tri] for tri in triangles] #[triangle 1 (3 x 3d), ...]
ax.add_collection(Poly3DCollection(points_per_tri, alpha=0.3))

#Formatting
ax.set_box_aspect(aspect=None, zoom=1.2)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("x")
ax.set_title("3D-reconstruction of object")
plt.show()

Visualise mesh using trimesh:

# ---- Interactive view mesh using trimesh ----
import trimesh

tri_mesh = trimesh.Trimesh(
    vertices=np.asarray(mesh.vertices),
    faces=np.asarray(mesh.triangles),
)
tri_mesh.show()