Структура от движения из двух видов

Структура от движения (SfM) является процессом оценки 3-D структуры сцены из набора 2-D изображений. В этом примере показано, как оценить положения калиброванной камеры из двух изображений, восстановить 3-D структуру сцены до неизвестного масштабного коэффициента, а затем восстановить фактический масштабный коэффициент путем обнаружения объекта известного размера.

Обзор

В этом примере показано, как восстановить 3-D сцену из пары изображений 2-D сделанных камерой, калиброванной с помощью Camera Calibrator. Алгоритм состоит из следующих шагов:

  1. Сопоставьте разреженный набор точек между двумя изображениями. Существует несколько способов нахождения соответствия точек между двумя изображениями. Этот пример обнаруживает углы в первом изображении, используя detectMinEigenFeatures function, и отслеживает их во второе изображение используя vision.PointTracker. Также можно использовать extractFeatures далее следуют matchFeatures.

  2. Оцените фундаментальную матрицу, используя estimateFundamentalMatrix.

  3. Вычислите движение камеры с помощью cameraPose функция.

  4. Совпадайте с плотным набором точек между двумя изображениями. Повторно обнаружить точку используя detectMinEigenFeatures с уменьшенной 'MinQuality' чтобы получить больше точки. Затем отследите плотные точки во второе изображение, используя vision.PointTracker.

  5. Определите 3-D местоположения совпадающих точек используя triangulate.

  6. Обнаружение объекта известного размера. В этой сцене есть глобус, радиус которого, как известно, составляет 10 см. Использование pcfitsphere чтобы найти глобус в облаке точек.

  7. Восстановите фактическую шкалу, что приведет к метрической реконструкции.

Чтение пары изображений

Загрузите пару изображений в рабочую область.

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');

Загрузка параметров камеры

Этот пример использует параметры камеры, вычисленные приложением Camera Calibrator. Параметры хранятся в cameraParams объект, и включают в себя внутреннюю составляющую камеры и коэффициенты искажения объектива s.helpview (fullfile (docroot, 'toolbox', 'matlab', 'matlab _ prog', 'matlab _ prog.map'), 'nested _ functions')

% Load precomputed camera parameters
load upToScaleReconstructionCameraParameters.mat

Удаление искажения объектива

Искажение объектива может повлиять на точность окончательной реконструкции. Вы можете удалить искажение из каждого из изображений, используя undistortImage функция. Этот процесс выпрямляет линии, которые изогнуты радиальным искажением линзы.

I1 = undistortImage(I1, cameraParams);
I2 = undistortImage(I2, cameraParams);
figure 
imshowpair(I1, I2, 'montage');
title('Undistorted Images');

Поиск соответствия точек между изображениями

Обнаружение хороших функций для отслеживания. Уменьшите 'MinQuality' обнаружить меньше точек, которые были бы более равномерно распределены по изображению. Если движение камеры не очень большое, то отслеживать с помощью алгоритма KLT - хороший способ установить соответствия точек.

% Detect feature points
imagePoints1 = detectMinEigenFeatures(im2gray(I1), 'MinQuality', 0.1);

% Visualize detected points
figure
imshow(I1, 'InitialMagnification', 50);
title('150 Strongest Corners from the First Image');
hold on
plot(selectStrongest(imagePoints1, 150));

% 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');

Оцените Фундаментальную Матрицу

Используйте estimateFundamentalMatrix функция, чтобы вычислить фундаментальную матрицу и найти inlier точки, которые удовлетворяют эпиполярному ограничению.

% Estimate the fundamental matrix
[fMatrix, epipolarInliers] = estimateFundamentalMatrix(...
  matchedPoints1, matchedPoints2, 'Method', 'MSAC', 'NumTrials', 10000);

% Find epipolar inliers
inlierPoints1 = matchedPoints1(epipolarInliers, :);
inlierPoints2 = matchedPoints2(epipolarInliers, :);

% Display inlier matches
figure
showMatchedFeatures(I1, I2, inlierPoints1, inlierPoints2);
title('Epipolar Inliers');

Вычисление положения камеры

Вычислите поворот и перемещение между позами камеры, соответствующими двум изображениям. Обратите внимание, что t является вектором модуля, потому что перемещение может быть вычислено только до масштаба.

[R, t] = cameraPose(fMatrix, cameraParams, inlierPoints1, inlierPoints2);

Восстановите 3-D местоположения совпадающих точек

Повторно обнаружить точки в первом изображении используя более низкие 'MinQuality' чтобы получить больше точки. Отследите новые точки во втором изображении. Оцените местоположения 3-D, соответствующие совпадающим точкам, используя triangulate функция, которая реализует алгоритм прямого линейного преобразования (DLT) [1]. Поместите источник в оптический центр камеры, соответствующий первому изображению.

% Detect dense feature points
imagePoints1 = detectMinEigenFeatures(im2gray(I1), 'MinQuality', 0.001);

% 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 X-axis. Thus, its
% rotation matrix is identity, and its translation vector is 0.
camMatrix1 = cameraMatrix(cameraParams, eye(3), [0 0 0]);
camMatrix2 = cameraMatrix(cameraParams, R', -t*R');

% 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);

Отображение 3-D облака точек

Используйте 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('Location', t, 'Orientation', R, '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');

Подбор сферы к облаку точек для поиска глобуса

Найдите глобус в облаке точек, подгоняя сферу к 3-D точкам с помощью 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 см. Теперь можно определить координаты точек 3-D в сантиметрах.

% Determine the scale factor
scaleFactor = 10 / globe.Radius;

% Scale the point cloud
ptCloud = pointCloud(points3D * scaleFactor, 'Color', color);
t = t * scaleFactor;

% Visualize the point cloud in centimeters
cameraSize = 2; 
figure
plotCamera('Size', cameraSize, 'Color', 'r', 'Label', '1', 'Opacity', 0);
hold on
grid on
plotCamera('Location', t, 'Orientation', R, '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');

Сводные данные

Этот пример показал, как восстановить движение камеры и восстановить 3-D структуру сцены из двух изображений, сделанных калиброванной камерой.

Ссылки

[1] Хартли, Ричард и Эндрю Зиссерман. Несколько видов геометрии в Компьютерное Зрение. Второе издание. Кембридж, 2000.