基于matlab的sfm三维重建

基于matlab的sfm三维重建

想学习三维重建，但是不会c++语言，python调用opencv弄得我头大，正好matlab上也有三维重建的代码，于是编跟着案例库的案例倒弄了一阵，大致弄明白了，分享一个流程给新手做参考。
ps：几个周以后再看自己做的这个，都是啥玩意，我都不好意思看了，事实上这个算不上稠密重建，纹理映射没有，由于手机是双摄标定出来也不准确，导致三角测量时相机位置偏的太多，在之后点云拼接的时候是一团乱麻，根本拼不上，毕竟新手起步，等我把整个完整流程弄完了再来分享。

1，加载图像和相机标定

相机参数需要先标定相机，matlab带有一个相机标定工具箱，很简单。一开始想做无标定三维重建的，但是方法我没找到。
本方法重建效果主要取决于三点：图像质量，相机校正好坏，特征匹配方法。

imageDir = fullfile('D:','picture','box');  % 这里输入路径，我重建的对象是一个盒子
imds = imageDatastore(imageDir);
load('calibrationSession');  % 加载相机参数，这个是matlab标定工具箱标定好后保存的mat文件
cameraParams = calibrationSession.CameraParameters;  %我们需要的是CameraParameters

将图像都转化为灰度图，特征检测只能检测灰度图。

images = cell(1, numel(imds.Files));
for i = 1:numel(imds.Files) % 图片的数量
    I = readimage(imds, i);
    images{i} = rgb2gray(I); 
end

实验室的盒子，一共拍摄了16张图片

2、创建图像追踪和配置颜色模板

创建图像追踪，原理就是先检测第一张图片所有的特征点，然后设置一个追踪器tracker，和图像管理器vSet，在剩下图片中追踪该特征点，然后储存在vSet中，同时vSet还储存相机位置等。（vSet貌似不能储存颜色信息，所以我才把颜色单独弄出来）

以第一二张图为基准读取点云信息，用于后面的颜色匹配。

I1 = readimage(imds, 1);% 读取第1张图片
I2 = images{2};  % 这里读取的不是灰度图
[currPoints, validIdx] = step(tracker, I2);
matchedPoints = PrevPoints(validIdx, :);

特征检测，最小特征值检测方法，速度稍慢，特征多，surf特征检测方法，速度快，检测特征少，KAZE很慢，列出来的原因是因为matlab没自带sift。

I = images{1};  % 选择第一张灰度图像
prevPoints = detectMinEigenFeatures(I, 'MinQuality', 0.001,'FilterSize', 5);  
% prevPoints = detectSURFFeatures(I,'MetricThreshold',5,'NumScaleLevels',6);  
% prevPoints = detectKAZEFeatures(I,'Diffusion','region','Threshold',...
% 0.00001,'NumOctaves',4,'NumScaleLevels',5);  

创建追踪器来追踪特征点

tracker = vision.PointTracker('MaxBidirectionalError', 1, 'NumPyramidLevels',...
6,'BlockSize',[31,31]);
% 初始化点跟踪器。
PrevPoints = prevPoints.Location;
initialize(tracker, PrevPoints, I);

将检测出来的点存放在vSet中

vSet = viewSet;  % matlab自带的视图管理器
viewId = 1;
vSet = addView(vSet, viewId, 'Points', prevPoints, 'Orientation', ...
    eye(3, 'like', PrevPoints), 'Location', ...
    zeros(1, 3, 'like', PrevPoints));

MinEigen最小特征法检测的稠密点匹配（8w多个点）

kaze特征检测的稠密点（4w多个点）。

3、重建主流程

这里不分稀疏重建和稠密重建，根据特征检测时把对应参数调高或者调低，来控制稀疏还是稠密。
注意，我这里只重建第一张图所包含的内容！！！
在这里追踪器开始产生作用

for i = 2:numel(images)
    I = images{i};
    
    [currPoints, validIdx] = step(tracker, I); %追踪第i张图片的特征点，返回特征点和点的对应关系
    
    matchedPoints1 = PrevPoints(validIdx, :);  % 筛选出来的匹配点
    matchedPoints2 = currPoints(validIdx, :);
    
    [E, epipolarInliers] = estimateEssentialMatrix(matchedPoints1,...
    matchedPoints2,cameraParams);
    %警告：已达到最大审判次数。考虑增加最大距离或降低期望的信心。

    inlierPoints1 = matchedPoints1(epipolarInliers, :);
    inlierPoints2 = matchedPoints2(epipolarInliers, :);
    
    [orient, loc] = relativeCameraPose(E, cameraParams,... 
    inlierPoints1,inlierPoints2);
    
    % 将当前视图匹配添加到视图集
    vSet = addView(vSet, i, 'Points', currPoints);
    % 存储上一个视图和当前视图之间的点匹配。
    
    matches = repmat((1:size(PrevPoints, 1))', [1, 2]);
    matches = matches(validIdx, :); 
    vSet = addConnection(vSet, i-1, i, 'Matches',matches);
    % 因为储存的点为PrevPoints和currPoints，要变成 matchedPoints1，2的对应关系。
    % 获取包含前一个相机姿势的表。
    prevPose = poses(vSet, i-1);
    prevOrientation = prevPose.Orientation{1};
    prevLocation = prevPose.Location{1};
    % 计算当前相机在全局坐标系中相对于第一个视图的姿态。
    orientation = orient * prevOrientation;
    location    = prevLocation + loc * prevOrientation;
    vSet = updateView(vSet, i, 'Orientation', orientation, 'Location', location);
    
    tracks = findTracks(vSet);   % 在所有视图中查找轨迹
    camPoses = poses(vSet);    % 得到所有视角的相机姿态表
    
    % 确定三维世界点的初始位置。
    xyzPoints = triangulateMultiview(tracks, camPoses, cameraParams);
     
    % 改进三维世界点和相机的姿势.(光束平差法)
    [xyzPoints, camPoses, reprojectionErrors] = bundleAdjustment(xyzPoints, ...
        tracks, camPoses, cameraParams, 'FixedViewId', 1, ...
        'PointsUndistorted', true);
    % 保存精准的相机姿势
    vSet = updateView(vSet, camPoses);
    PrevPoints   = currPoints;  
end  

4、添加颜色

原理很简单，就是找三维点在像素平面中对应位置的颜色信息，然后添加上去。

numPixels = size(I1, 1) * size(I1, 2);   %读取图片尺寸
allColors = reshape(I1, [numPixels, 3]);   % 改变数组的维数。
colorIdx = sub2ind([size(I1, 1), size(I1, 2)], round(matchedPoints(:,2)), ...
round(matchedPoints(:, 1)));
color = allColors(colorIdx, :);
% 创建点云
goodIdx = (reprojectionErrors < 15);  % 删除误差大的点
xyzPoints = xyzPoints(goodIdx, :);
Color = color(goodIdx, :);
ptCloud = pointCloud(xyzPoints, 'Color', Color);   % 将点云添加上颜色

5、显示点云图

% 显示相机姿态
% camPoses = poses(vSet);
% figure;
% plotCamera(camPoses, 'Size', 0.2);
% hold on

% 显示三维点.
pcshow(ptCloud, 'VerticalAxis', 'y', 'VerticalAxisDir', 'down', ...
    'MarkerSize', 45);
grid on
hold off
% 查看指定卷，免得图像乱跑
loc1 = camPoses.Location{1};
xlim([loc1(1)-5, loc1(1)+4]);
ylim([loc1(2)-5, loc1(2)+4]);
zlim([loc1(3)-1, loc1(3)+20]);
camorbit(0, -30);
% 储存点云
a = cell(1,i);
a{1} = ptCloud;
save('PtCloud.mat','a');

6、重建结果

minEigen特征

surf特征

kaze特征

目前遇到问题有两个，一个是颜色匹配有点误差，还有一个是太占内存了，如果电脑内存不够可能会出现错误。
接下来研究一下如何重建一个完整的盒子，考虑再用一个循环，然后点云拼接，本人新手上路，有大神的话可以指点一下。
附带我的程序和图片：sfm三维重建