相机扭曲实验

二维图像扭曲旋转三维图像


    1.匹配两个图像之间的稀疏点集。有两种方法可以找到两个图像之间的点对应关系。
该示例使用检测min特征特征函数检测第一图像中的角点，并使用视觉点跟踪器将
它们跟踪到第二图像中。或者，您可以使用提取特征，后跟匹配特征。
    2.提取计算估值
    3.利用摄像机的姿态函数计算摄像机的运动
    4.匹配两个图像之间密集的点集。通过检测最小特征值来重新检测点，减少最小质量
以获得更多的点。然后利用视觉跟踪点将稠密点跟踪到第二幅图像。
    5.确定匹配点的三维位置，利用三角。
    6.已知物体的大小。在这个场景中有一个球体，其半径设定定值10cm。使用pcfitsphere在点云找到地球。
    7.恢复实际的尺度，从而进行度量重建。


步骤：




%读两张图片
imageDir = fullfile(toolboxdir('vision'), 'visiondata','upToScaleReconstructionImages');
images = imageSet(imageDir);
I1 = read(images, 1);
I2 = read(images, 2);
figure
imshowpair(I1, I2, 'montage');
title('Original Images');


负荷的相机参数


    本例使用摄像机校准程序计算的摄像机参数。参数存储在相机参数对象，包括相机内部
和镜头畸变系数。


% Load precomputed camera parameters负荷计算摄像机参数
load upToScaleReconstructionCameraParameters.mat


%删除镜头畸变
%透镜畸变会影响最终重建的精度。从每个失真图像功能消除了失真。这个过程就直接由镜头的径向转变弯曲的线变直。


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


%查找图像之间的对应关系
%%检测好的功能来跟踪。 减少“MinQuality”以检测更少的点，这将在整个图像中更均匀地分布。
 如果相机的运动不是很大，那么使用KLT算法进行跟踪是建立点对应的好方法。


% Detect feature points 检测特征点
imagePoints1 = detectMinEigenFeatures(rgb2gray(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');




%基本矩阵的估计
%使用估计基本矩阵函数来计算基本矩阵和找到内点满足极线约束。


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




重建匹配点的三维位置
使用较低的“MinQuality”重新检测第一个图像中的点，以获得更多的点。 跟踪新的点到第二个图像。
使用实现直接线性变换（DLT）算法的三角函数来估计与匹配点对应的三维位置。 将原点放置在与第
一个图像相对应的相机的光学中心处。


% Detect dense feature points 检测密集的特征点
imagePoints1 = detectMinEigenFeatures(rgb2gray(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.
计算相机每个位置的相机矩阵第一台摄像机位于X轴的原点。 因此，它的
旋转矩阵是一致的，其平移向量为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);


显示三维点云
使用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');


%%将一个球体装配到点云上以找到全球
通过使用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);
t = t * scaleFactor;


% Visualize the point cloud in centimeters 点云的单位cm
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');