稀疏 3D 重建 MATLAB 示例

Question

我有一个立体相机系统，我正在使用我自己的图像和相机校准文件尝试这个 MATLAB 的计算机视觉工具箱示例 (http://www.mathworks.com/help/vision/ug/sparse-3-d-reconstruction-from-multiple-views.html)。我使用了加州理工学院的相机校准工具箱 (http://www.vision.caltech.edu/bouguetj/calib_doc/)。

首先，我根据第一个示例分别尝试了每个相机，并找到了每个相机的固有相机校准矩阵并保存了它们。我还使用 Caltech 工具箱对左右图像进行了不失真处理。因此，我从 MATLAB 示例中注释掉了该代码。

这是内在的相机矩阵：

K1=[1050 0 630;0 1048 460;0 0 1];

K2=[1048 0 662;0 1047 468;0 0 1];

顺便说一句，这些是 Bumblebee XB3 相机的右侧和中间镜头。

问题：它们不应该是一样的吗？

然后我根据第五个示例进行了立体校准。我从中保存了旋转矩阵（R）和平移矩阵（T）。因此，我将 MATLAB 示例中的代码注释掉了。

这里是旋转和平移矩阵：

R=[0.9999 -0.0080 -0.0086;0.0080 1 0.0048;0.0086 -0.0049 1];

T=[120.14 0.55 1.04];

然后我将所有这些图像、校准文件和相机矩阵输入到 MATLAB 示例中，并试图找到 3-D 点云，但结果并不乐观。我在这里附上代码。我认为这里有两个问题：

1- 我的极线约束值太大！（16次方）

2- 我不确定相机矩阵以及我如何从 R 和 Caltech 工具箱中的 T 计算它们！

附：就特征提取而言，这很好。

如果有人能帮忙就好了。

clear
close all
clc

files = {'Left1.tif';'Right1.tif'};
for i = 1:numel(files)
files{i}=fullfile('...\sparse_matlab', files{i});
images(i).image = imread(files{i});
end
figure;
montage(files); title('Pair of Original Images')

% Intrinsic camera parameters
load('Calib_Results_Left.mat')
K1 = KK;
load('Calib_Results_Right.mat')
K2 = KK;

%Extrinsics using stereo calibration
load('Calib_Results_stereo.mat')
Rotation=R;
Translation=T';
images(1).CameraMatrix=[Rotation; Translation] * K1;
images(2).CameraMatrix=[Rotation; Translation] * K2;

% Detect feature points and extract SURF descriptors in images
for i = 1:numel(images)
%detect SURF feature points
images(i).points = detectSURFFeatures(rgb2gray(images(i).image),...
    'MetricThreshold',600);
%extract SURF descriptors
[images(i).featureVectors,images(i).points] = ...
    extractFeatures(rgb2gray(images(i).image),images(i).points);
end

% Visualize several extracted SURF features from the Left image
figure; imshow(images(1).image);
title('1500 Strongest Feature Points from Globe01');
hold on;
plot(images(1).points.selectStrongest(1500));

indexPairs = ...
matchFeatures(images(1).featureVectors, images(2).featureVectors,...
'Prenormalized', true,'MaxRatio',0.4) ;
matchedPoints1 = images(1).points(indexPairs(:, 1));
matchedPoints2 = images(2).points(indexPairs(:, 2));
figure;

% Visualize correspondences
showMatchedFeatures(images(1).image,images(2).image,matchedPoints1,matchedPoints2,'montage'    );
title('Original Matched Features from Globe01 and Globe02');

% Set a value near zero, It will be used to eliminate matches that
% correspond to points that do not lie on an epipolar line.
epipolarThreshold = .05;
for k = 1:length(matchedPoints1)
% Compute the fundamental matrix using the example helper function
% Evaluate the epipolar constraint
epipolarConstraint =[matchedPoints1.Location(k,:),1]...
    *helperCameraMatricesToFMatrix(images(1).CameraMatrix,images(2).CameraMatrix)...
    *[matchedPoints2.Location(k,:),1]';

%%%% here my epipolarConstraint results are bad %%%%%%%%%%%%%

% Only consider feature matches where the absolute value of the
% constraint expression is less than the threshold.
valid(k) = abs(epipolarConstraint) < epipolarThreshold;
end

validpts1 = images(1).points(indexPairs(valid, 1));
validpts2 = images(2).points(indexPairs(valid, 2));
figure;
showMatchedFeatures(images(1).image,images(2).image,validpts1,validpts2,'montage');
title('Matched Features After Applying Epipolar Constraint');

% convert image to double format for plotting
doubleimage = im2double(images(1).image);
points3D = ones(length(validpts1),4); % store homogeneous world coordinates
color = ones(length(validpts1),3);    % store color information

% For all point correspondences
for i = 1:length(validpts1)
% For all image locations from a list of correspondences build an A
pointInImage1 = validpts1(i).Location;
pointInImage2 = validpts2(i).Location;
P1 = images(1).CameraMatrix'; % Transpose to match the convention in
P2 = images(2).CameraMatrix'; % in [1]
A = [
    pointInImage1(1)*P1(3,:) - P1(1,:);...
    pointInImage1(2)*P1(3,:) - P1(2,:);...
    pointInImage2(1)*P2(3,:) - P2(1,:);...
    pointInImage2(2)*P2(3,:) - P2(2,:)];

% Compute the 3-D location using the smallest singular value from the
% singular value decomposition of the matrix A
[~,~,V]=svd(A);

X = V(:,end);
X = X/X(end);

% Store location
points3D(i,:) = X';

% Store pixel color for visualization
y = round(pointInImage1(1));
x = round(pointInImage1(2));
color(i,:) = squeeze(doubleimage(x,y,:))';
end

% add green point representing the origin
points3D(end+1,:) = [0,0,0,1];
color(end+1,:) = [0,1,0];

% show images
figure('units','normalized','outerposition',[0 0 .5 .5])
subplot(1,2,1); montage(files,'Size',[1,2]); title('Original Images')

% plot point-cloud
hAxes = subplot(1,2,2); hold on; grid on;
scatter3(points3D(:,1),points3D(:,2),points3D(:,3),50,color,'fill')
xlabel('x-axis (mm)');ylabel('y-axis (mm)');zlabel('z-axis (mm)')
view(20,24);axis equal;axis vis3d
set(hAxes,'XAxisLocation','top','YAxisLocation','left',...
'ZDir','reverse','Ydir','reverse');
grid on
title('Reconstructed Point Cloud');

Answer 1

首先，计算机视觉系统工具箱现在包括用于校准单个相机的Camera Calibrator App，以及support for programmatic stereo camera calibration。您使用这些工具会更容易，因为您使用的示例和 Caltech Calibration Toolbox 使用的约定有些不同。

该示例使用前乘约定，即行向量 * 矩阵，而 Caltech 工具箱使用后乘约定（矩阵 * 列向量）。这意味着如果您确实使用加州理工学院的相机参数，您将不得不转置内在矩阵和旋转矩阵。这可能是您的问题的主要原因。

就您的两台相机之间的内在特性而言，这是完全正常的。所有相机都略有不同。

这也有助于查看您用于三角测量的匹配特征。鉴于您正在重建一个细长的物体，看到重建的点在 3D 中形成一条线似乎并不令人惊讶......

您也可以尝试校正图像并进行密集重建，如example I've linked to above。