基于两个视图的运动重建
运动重建 (SfM) 是基于一组二维图像估计场景的三维结构的过程。本示例将向您展示如何根据两个图像估计已标定相机的位姿，在未知缩放因子条件下重建场景的三维结构，然后通过检测一个已知尺寸的目标来恢复实际的缩放因子。
概述
本示例演示了如何利用一组由经过相机标定器应用标定的相机拍摄的二维图像，重建三维场景。该算法包括以下步骤：
在两幅图像之间匹配一组稀疏点。有多种方法可以找到两幅图像之间的点对应关系。此示例使用
detectSIFTFeatures函数检测第一个图像中的角点，并使用vision.PointTracker函数将这些角点追踪到第二个图像中。或者，您可以先使用extractFeatures，然后使用matchFeatures。使用
estimateEssentialMatrix估计基础矩阵。使用
estrelpose函数计算相机的运动。在两幅图像之间匹配一组稠密的点。使用
detectSIFTFeatures并缩小ContrastThreshold的范围，重新检测该点，以获取更多点。然后使用vision.PointTracker将稠密点追踪到第二个图像中。使用
triangulateMultiview确定匹配点的三维位置。检测一个已知尺寸的目标。在这个场景中有一个地球仪，已知其半径为 10 厘米。使用
pcfitsphere在点云中查找地球模型。恢复实际比例尺，从而实现度量重建。
读取一组图像
将一组图像加载到工作区中。
imageDir = fullfile(toolboxdir("vision"),"visiondata","upToScaleReconstructionImages"); images = imageDatastore(imageDir); I1 = readimage(images, 1); I2 = readimage(images, 2); figure imshowpair(I1, I2, 'montage'); title("Original Images");
加载相机参数
本示例使用了由相机标定器计算出的相机参数。这些参数存储在 cameraIntrinsics 对象中，包括相机内参和镜头畸变系数。
% Load precomputed camera intrinsics data = load("sfmCameraIntrinsics.mat"); intrinsics = data.intrinsics;
消除镜头畸变
镜头畸变可能会影响最终重建结果的准确性。您可以使用 undistortImage 函数消除每个图像中的畸变。该过程可矫正因镜头径向畸变而产生的线条弯曲。
I1 = undistortImage(I1, intrinsics); I2 = undistortImage(I2, intrinsics); figure imshowpair(I1, I2, "montage"); title("Undistorted Images");
找出两幅图像之间的点对应关系
检测适合追踪的特征。如果相机的位移不大，那么使用 KLT 算法进行跟踪是建立点对应关系的一种有效方法。
% Detect feature points. Use an ROI to eliminate spurious % features around the edges of the image. border = 50; roi = [border, border, size(I1, 2)- 2*border, size(I1, 1)- 2*border]; imagePoints1 = detectSIFTFeatures(im2gray(I1), ROI=roi, ContrastThreshold=0.015); % Visualize detected points figure imshow(I1); title("150 Strongest Corners from the First Image"); hold on plot(selectStrongest(imagePoints1, 1000), ShowScale=false, ShowOrientation=false);
% Create the point tracker tracker = vision.PointTracker(MaxBidirectionalError=1, NumPyramidLevels=5); % Initialize the point tracker imagePoints1 = imagePoints1.Location; initialize(tracker, imagePoints1, I1); % Track the points [imagePoints2, validIdx] = step(tracker, I2); matchedPoints1 = imagePoints1(validIdx, :); matchedPoints2 = imagePoints2(validIdx, :); % Visualize correspondences figure showMatchedFeatures(I1, I2, matchedPoints1, matchedPoints2); title("Tracked Features");
求算本质矩阵
使用 estimateEssentialMatrix 函数计算本质矩阵，并找出满足极线约束条件的内点。
% Estimate the fundamental matrix [E, epipolarInliers] = estimateEssentialMatrix(... matchedPoints1, matchedPoints2, intrinsics); % Find epipolar inliers inlierPoints1 = matchedPoints1(epipolarInliers, :); inlierPoints2 = matchedPoints2(epipolarInliers, :); % Display inlier matches figure showMatchedFeatures(I1, I2, inlierPoints1, inlierPoints2); title("Epipolar Inliers");
计算相机位姿
计算第二台相机相对于第一台相机的位置和方向。请注意，loc 是一个平移单位向量，因为平移只能计算到比例为止。
relPose = estrelpose(E, intrinsics, inlierPoints1, inlierPoints2);
重建匹配点的三维位置
使用较低的 ContrastThreshold 值对第一个图像中的点进行重新检测，以获取更多点。将新点跟踪到第二个图像中。使用 triangulate 函数（该函数实现了直接线性变换 (DLT) 算法 [1]）来估计与配对点对应的三维位置。将原点置于与第一幅图像相对应的相机光心处。
% Detect dense feature points. imagePoints1 = detectSIFTFeatures(im2gray(I1), ROI = roi, ContrastThreshold=0.005); % Create the point tracker tracker = vision.PointTracker(MaxBidirectionalError=1, NumPyramidLevels=5); % Initialize the point tracker imagePoints1 = imagePoints1.Location; initialize(tracker, imagePoints1, I1); % Track the points [imagePoints2, validIdx] = step(tracker, I2); matchedPoints1 = imagePoints1(validIdx, :); matchedPoints2 = imagePoints2(validIdx, :); % Compute the camera matrices for each position of the camera % The first camera is at the origin looking along the Z-axis. Thus, its % transformation is identity. camMatrix1 = cameraProjection(intrinsics, rigidtform3d); camMatrix2 = cameraProjection(intrinsics, pose2extr(relPose)); % Compute the 3-D points points3D = triangulate(matchedPoints1, matchedPoints2, camMatrix1, camMatrix2); % Get the color of each reconstructed point numPixels = size(I1, 1) * size(I1, 2); allColors = reshape(I1, [numPixels, 3]); colorIdx = sub2ind([size(I1, 1), size(I1, 2)], round(matchedPoints1(:,2)), ... round(matchedPoints1(:, 1))); color = allColors(colorIdx, :); % Create the point cloud ptCloud = pointCloud(points3D, Color=color);
显示三维点云
使用 plotCamera 函数可视化相机的位置和方向，使用 pcshow 函数可视化点云。
% Visualize the camera locations and orientations cameraSize = 0.3; figure plotCamera(Size=cameraSize, Color="r", Label="1", Opacity=0); hold on grid on plotCamera(AbsolutePose=relPose, Size=cameraSize, ... Color="b", Label="2", Opacity=0); % Visualize the point cloud pcshow(ptCloud, VerticalAxis="y", VerticalAxisDir="down", MarkerSize=45); % Rotate and zoom the plot camorbit(0, -30); camzoom(1.5); % Label the axes xlabel("x-axis"); ylabel("y-axis"); zlabel("z-axis") title("Up to Scale Reconstruction of the Scene");
将球体拟合到点云上以确定地球的形状
使用 pcfitsphere 函数将一个球体拟合到三维点云上，从而在点云中找到地球。
% Detect the globe globe = pcfitsphere(ptCloud, 0.1); % Display the surface of the globe plot(globe); title("Estimated Location and Size of the Globe"); hold off
场景的度量重建
地球的实际半径为 10 厘米。现在，您可以以厘米为单位确定三维点的坐标了。
% Determine the scale factor scaleFactor = 10 / globe.Radius; % Scale the point cloud ptCloud = pointCloud(points3D * scaleFactor, Color=color); relPose.Translation = relPose.Translation * scaleFactor; % Visualize the point cloud in centimeters cameraSize = 2; figure plotCamera(Size=cameraSize, Color="r", Label="1", Opacity=0); hold on grid on plotCamera(AbsolutePose=relPose, Size=cameraSize, ... Color="b", Label="2", Opacity=0); % Visualize the point cloud pcshow(ptCloud, VerticalAxis="y", VerticalAxisDir="down", MarkerSize=45); camorbit(0, -30); camzoom(1.5); % Label the axes xlabel("x-axis (cm)"); ylabel("y-axis (cm)"); zlabel("z-axis (cm)") title("Metric Reconstruction of the Scene");
总结
例如，本示例向您展示了如何通过两张由经过标定的相机拍摄的图像，恢复相机运动并重建场景的三维结构。
参考资料
[1] Hartley, Richard, and Andrew Zisserman.Multiple View Geometry in Computer Vision.Second Edition.Cambridge, 2000.
另请参阅
opticalFlowRAFT | detectSIFTFeatures | estimateEssentialMatrix | estrelpose
