Dense 3-D Reconstruction of Asteroid Surface from Image Sequence

R2026b

Since R2026a

This example shows how to reconstruct a dense 3-D model of an asteroid from a sequence of images captured by NASA’s OSIRIS-REx spacecraft during its approach to the asteroid Bennu [1]. The workflow demonstrates a complete dense reconstruction pipeline using only RGB imagery. Starting with the camera intrinsics and the image sequence, the example first performs incremental structure-from-motion (SfM) to estimate camera poses, a view graph, and a sparse 3-D reconstruction. Dense optical flow from the RAFT deep learning model is then used to compute depth maps between overlapping views, which are filtered and combined to form a dense colored point cloud of the asteroid surface. Finally, Poisson surface reconstruction is applied to generate a smooth, watertight 3-D mesh suitable for visualization and further analysis.

Download Images

This example uses a subset of the OSIRIS-REx Camera Suite (OCAMS) Bundle [1] containing 37 images of the surface of the asteroid Bennu. You can download the data to a temporary directory using a web browser or by running the following code:

% Location of the compressed data set
url = "https://ssd.mathworks.com/supportfiles/3DReconstruction/Bennu_preliminary_survey.zip";

% Store the data set in a temporary folder
downloadFolder = tempdir;
filename = fullfile(downloadFolder, "Bennu_preliminary_survey.zip");

% Uncompressed data set
imageFolder = fullfile(downloadFolder, "Bennu_preliminary_survey");

if ~exist(imageFolder, "dir") % download only once
    disp('Downloading Asteroid Dataset (15 MB)...');
    websave(filename, url);
    unzip(filename, downloadFolder);
end

Create an image datastore to load the images of the asteroid.

imds = imageDatastore(imageFolder);

The images are captured by the MapCam Camera in the OSIRIS-REX Camera Suite. MapCam has an effective focal length of 125 mm [2]. Convert this focal length into pixel units and create a cameraIntrinsicsToOpenCV object to store the intrinsic parameters of the camera.

pixelSize = 8.5e-3; % in millimeters
focalLength = 125/pixelSize;
intrinsics = cameraIntrinsics([focalLength, focalLength], [512, 512], [1024, 1024]);

Display sample images.

figure
montage(imds, Size=[3 3])
title("Sample Images")
truesize;

Sparse Reconstruction Using Structure from Motion

Reconstruct a sparse structure of the asteroid from the images using SfM. See the Structure from Motion from Multiple Views example for more in-depth exploration of the pipeline. The entire pipeline can take 5-10 minutes to complete.

Create a view graph connecting images with overlapping views of the asteroid. NumSimilarImages controls how many visually similar images are connected in the view graph. Larger values increase connectivity but also computation. For small datasets, 10–20 is typical. For large datasets, you may need 30–50.

sfmObj = sfm(imds,intrinsics);
sfmObjViewGraph = connectImagePairs(sfmObj,NumSimilarImages=20);
sfmObjVerified = verifyImagePairs(sfmObjViewGraph,MinNumInliers=30,MaxDistance=4,Verbose=true);

82 out of 274 edges passed geometric verification.
Average 2D inliers count: 1006.926829
Minimum 2D inliers count: 68
Maximum 2D inliers count: 2337

Initialize the 3-D reconstruction using a pair of views. Decrease MinMedianAngle and increase MaxTriangulationError to relax the conditions for the initial view selection if the default values result in an unsuccessful triangulation.

sfmObjInit = triangulateInitialViews(sfmObjVerified,MinMedianAngle=10,MaxTriangulationError=8,Verbose=true);

79 of 82 edges have more than 100 feature matches.
69 of 79 edges are selected after applying the grid threshold.
Views 1 and 2 are selected for initial triangulation.
Median triangulation angle: 10.313281

Incrementally reconstruct the scene and estimate camera positions of the remaining views of the asteroid. Decrease MinNumInliers in challenging cases where the reconstruction would otherwise remain incomplete.

sfmObjFinal = reconstruct(sfmObjInit,MinNumInliers=10,Verbose=true);

Retrieve the camera positions estimated using structure-from-motion.

camPoses = poses(sfmObjFinal);
sparsePoints = pointCloud(sfmObjFinal);

Plot the cameras and sparse 3-D points. The cameras appear to follow a "spiral trajectory" because the camera was approaching the asteroid while it was spinning.

figure
pcshow(sparsePoints)
xlabel("X")
ylabel("Y")
zlabel("Z")
hold on
plotCamera(camPoses, Color="green", Size=0.1, Opacity=0.1)

View Selection for Dense Depth Estimation

Instantiate an opticalFlowRAFT model for computing dense pixel correspondences between image pairs.

raft = opticalFlowRAFT;

For each view in the image collection, select its best neighboring view for dense reconstruction. Dense reconstruction depends on finding pixel-wise matches between two images and then triangulating to obtain 3D points in world-coordinates for each of the matched pixels. A numerically stable triangulation computation requires the two views to have sufficient viewpoint variation between the two camera positions, along with a minimum number of mutually visible image features.

The minOverlap threshold rejects view pairs that are too far apart and do not share common features. The ideal value is dataset dependent and has to be tuned to obtain a reasonable number of pairs covering multiple view points for subsequent dense reconstruction. Reduce this threshold to consider more image pairs as valid for dense depth estimation. More fine-grained control of the view selection process can be performed by modifying the default values specified in the included supporting function helperSelectDepthViewPairs.

minOverlap = 0.3;
validViewPairs = helperSelectDepthViewPairs(sfmObjFinal.ViewGraph, sfmObjFinal.WorldPoints, minOverlap);

numValidViewPairs = size(validViewPairs,1);
numViews = sfmObjFinal.NumImages;
numImages = length(imds.Files);
numConn = sfmObjFinal.ViewGraph.NumConnections;
disp(numValidViewPairs + " valid views out of " + numViews+".");

33 valid views out of 37.

Visualize Views Selected for Dense Reconstruction

Visualize the selected views, along with the sparse 3D points from SfM, to ensure that the object of interest is being captured from multiple view points.

figure
pcshow(sparsePoints);
xlabel("X")
ylabel("Y")
zlabel("Z")
hold on
for i = 1:numValidViewPairs
    viewId = validViewPairs(i,1);
    plotCamera(AbsolutePose=camPoses(viewId), Color="green", Size=0.1, Opacity=0.1);
end
hold off

Estimate Depth Map from Image Pair

Use the slider to select a particular image pair index, pairIndex, and compute a depth map of the first image using dense optical flow and triangulation from relative camera pose.

pairIndex = 10;

maxDepthThreshold = 7;     % reject depth at points that are too far away to be estimated correctly

viewPair = validViewPairs(pairIndex,:);

Extract 3D World Points seen in view1 from the SfM results.

[pointIndices, featureIndices] = findWorldPointsInView(sfmObjFinal.WorldPoints, viewPair(1));
xyzPoints = sfmObjFinal.WorldPoints.WorldPoints(pointIndices,:);

Extract 2D keypoints detected in view1 from the SfM results.

view1 = findView(sfmObjFinal.ViewGraph, viewPair(1));
view2 = findView(sfmObjFinal.ViewGraph, viewPair(2));
featurePoints = view1.Points{:}(featureIndices);

Load the two images.

I1 = readimage(imds, viewPair(1));
I2 = readimage(imds, viewPair(2));

Compute relative pose between camera positions of the two views.

pose1 = view1.AbsolutePose;
pose2 = view2.AbsolutePose;
p1Inv = invert(pose1);
relPose = p1Inv.A * pose2.A;
relPose = rigidtform3d(relPose);

Compute depth map using dense pixel correspondences from optical flow and the relative camera pose between two views.

[depthMap, flow, isValid] = helperComputeDepthMap(...
    raft, I1, I2, relPose, intrinsics, featurePoints, maxDepthThreshold);

Scale the estimated depth map values to the same coordinates as the sparse SfM keypoints.

scaleFactor = helperComputeDepthMapScale(featurePoints,pose1,xyzPoints,depthMap);
depthMap = scaleFactor * depthMap;

Use graythresh and imbinarize to compute a binarization mask on I1 and apply to the depth map. Since the background is black and visible parts of the asteroid are strongly illuminated in this image sequence, thresholding the image removes most of the background areas.

I1Gray = im2gray(I1);
otsuT = graythresh(I1Gray);
BW = imbinarize(I1Gray,otsuT);
depthMap(BW==0) = nan;

Visualize the two images and the estimated depth map. Invalid depth map pixels are marked using NaN values.

figure
subplot(1,3,1)
imshow(I1)
title("Image " + viewPair(1))
subplot(1,3,2)
imshow(I2)
title("Image " + viewPair(2))
subplot(1,3,3)
imagesc(depthMap)
axis image
colorbar
xticks([])
yticks([])
title("Depth Image " + viewPair(1))
truesize;

Dense Reconstruction on Entire Image Sequence

Run on all the valid image pairs in the image sequence and reconstruct a dense point cloud. The code below implements a simple and straightforward approach that accumulates all the two-view reconstructions from all the selected views, along with their color or intensity information, into a single point cloud. Depth estimates at the image boundaries are often inconsistent and are pruned out. Depth values for extremely distant points are also error prone and are discarded. The maximum threshold on the acceptable depth values is data dependent and can be chosen based on visual inspection of the visualized depth map.

This step takes about 8 minutes on a Linux machine with an NVIDIA GeForce RTX 3090 GPU.

% Specify parameters
offset = 50;                         % Ignore image points at the image boundary
depthRange = [0 maxDepthThreshold];  % Ignore depth values that are too far away

ptCloudsFull = {};
imageSize = size(I1);
[X, Y] = meshgrid(offset:2:imageSize(2)-offset, offset:2:imageSize(1)-offset);
numNan = 0;
depthFactor = 1;

for i = 1:numValidViewPairs
    viewPair = validViewPairs(i,:);
    
    fprintf(i + "/" + numValidViewPairs + " ");

    I1 = readimage(imds, viewPair(1));
    I2 = readimage(imds, viewPair(2));

    view1 = findView(sfmObjFinal.ViewGraph, viewPair(1));
    view2 = findView(sfmObjFinal.ViewGraph, viewPair(2));
    
    % 3D world points visible from view1
    [pointIndices, featureIndices] = findWorldPointsInView(sfmObjFinal.WorldPoints, viewPair(1));
    xyzPoints = sfmObjFinal.WorldPoints.WorldPoints(pointIndices,:);
    
    % 2D keypoints for visible world points in view1
    featurePoints = view1.Points{:}(featureIndices);

    % Camera poses of the two views
    pose1 = view1.AbsolutePose;
    pose2 = view2.AbsolutePose;
    
    % Compute relative pose between cameras
    p1Inv = invert(pose1);
    relPose = p1Inv.A * pose2.A;
    relPose = rigidtform3d(relPose);
    
    % Compute depth map using dense pixel correspondences from optical flow
    % and relative camera pose between two views from SfM.
    [depthMap, flow, isValid] = helperComputeDepthMap(...
        raft, I1, I2, relPose, intrinsics, featurePoints, maxDepthThreshold);

    if ~isValid
        disp("Numerically unstable depth estimation.")
        numNan = numNan + 1;
        continue;
    end
    
    % Scale the estimated depth map values to the same coordinates as the
    % sparse SfM keypoints.
    scaleFactor = helperComputeDepthMapScale(featurePoints,pose1,xyzPoints,depthMap);
    depthMap = scaleFactor * depthMap;

    % Foreground selection- compute a binarization mask on I1 and apply to
    % the depth map.
    I1Gray = im2gray(I1);
    otsuT = graythresh(I1Gray);
    BW = imbinarize(I1Gray,otsuT);
    depthMap(BW==0) = nan;
    
    [xyzPoints, validIndex] = helperReconstructFromRGBD([X(:), Y(:)], ...
        depthMap, intrinsics, pose1, depthFactor, depthRange);

    colors = zeros(numel(X), 1, 'like', I1);
    for j = 1:numel(X)
        colors(j, 1:3) = I1(Y(j), X(j), :);
    end
    ptCloudsFull{end+1} = pointCloud(xyzPoints, Color=colors(validIndex, :));

end

Aggregate Reconstructed Point Clouds

Merge the point clouds into a single point cloud by using the pccat (Point Cloud Toolbox) function.

ptCloudsMerged = pccat(cat(1,ptCloudsFull{:}));

Apply a denoising step on the point cloud by using the pcdenoise (Point Cloud Toolbox) function.

ptCloudsMerged = pcdenoise(ptCloudsMerged);

Keep only the largest cluster in the point cloud and remove isolated components that likely result from depth map merging errors. Set clusterDistance to define the maximum distance between points in the same cluster. You can determine an appropriate clusterDistance value by visually inspecting the cleaned point cloud.

clusterDistance = 0.01;
[labels, numClusters] = pcsegdist(ptCloudsMerged, clusterDistance);
clusterSizes = zeros(numClusters, 1);
for i = 1:numClusters
    clusterSizes(i) = sum(labels == i);
end
[~, largestClusterLabel] = max(clusterSizes);
largestClusterIndices = find(labels == largestClusterLabel);
ptCloudClean = select(ptCloudsMerged, largestClusterIndices);

Visualize the point cloud using pcviewer (Point Cloud Toolbox).

figure
pcviewer(ptCloudClean)

Compute Surface Mesh from Dense Point Cloud

After generating a dense colored point cloud from the input images, you can convert it into a continuous 3-D surface mesh using Poisson surface reconstruction by using the pc2surfacemesh (Point Cloud Toolbox) function. The resulting mesh provides a smooth, topologically consistent representation of the asteroid, suitable for visualization, measurement, and further geometric analysis. This step can take several minutes to complete due to the large size of the point cloud being meshed.

mesh = pc2surfacemesh(ptCloudClean,"poisson");

Display the reconstructed surface mesh.

viewer = viewer3d(BackgroundColor="black",BackgroundGradient="off",RenderingQuality="high");

surfaceMeshShow(mesh, Parent=viewer);

References

[1] OSIRIS-REx Camera Suite (OCAMS) Bundle https://arcnav.psi.edu/urn:nasa:pds:orex.ocams

[2] ORX OCAMS Instrument Kernel https://naif.jpl.nasa.gov/pub/naif/ORX/kernels/ik/orx_ocams_v04.ti