结构光逆相机法重建详解+代码

代码地址：在公众号「3D视觉工坊」，后台回复「逆相机法」，即可直接下载。

注：本文的理论主要来自于参考文献1、2，代码来源于南京理工大学左超老师课题组发表的参考文献2中，笔者在其基础上稍作修改，便于大家更好理解。

逆相机法，也称为三角立体模型，其将投影仪看做“逆相机”，投影结构化光，主动标记视场内的“同名点”，利用类似双目视差原理（不完全相同）进行重建。

01 理论

1.1 单目标定

传感器存在误差，这部分请参考：​​一文图解单目相机标定算法​​。

此外，Matlab$标定程序还会计算它的投影矩阵 P，这也是我们需要的：​

1.3 重建原理

我们将投影仪看做一台“逆向”的相机，假设经过系统标定后，我们已经得到了 投影矩阵（这里换个记号）：

如我们之前所说，从世界坐标系 -> 像素坐标系，有以下关系：

那么世界坐标系、相机/投影仪的像素坐标系之间有如下关系：

需要说明的是，这个公式是对每个像素点单独计算的，这在代码里可以体现。

再来分析公式的精度跟哪些因素有关：

这种方法的优点是比较简单，并且在标定范围外进行测量，精度也比较高，缺点是：将投影仪当做逆相机，两个模型明显会不一样，显然不可能做到非常高精度。

02 实践

2.1 标定

标定系统参数：

1. 相机、投影仪的成像误差
2. 相机、投影仪的相对位置关系

(1) 相机标定

代码中，第一步，标定相机：

• 
  

num_x = 11; %number of circles in the x direction of the calibration board
num_y = 9; %number of circles in the y direction of the calibration board
dist_circ = 25; % 标定板


disp('开始相机标定...');

% 角点的像素坐标系
load('camera_imagePoints.mat'); % load the camera image points of the centers

% 角点的世界坐标系
worldPoints = generateCheckerboardPoints([num_y+1,num_x+1], dist_circ);

% 标定相机
[cameraParams,imagesUsed,estimationErrors] = estimateCameraParameters(imagePoints,worldPoints, 'EstimateTangentialDistortion', true);

figure;showReprojectionErrors(cameraParams);title('相机标定的重投影误差');
figure;showExtrinsics(cameraParams);title('相机标定的外参');

重投影误差在 0.08 pixel，精度还是相当高的，这是因为使用高精度标定板的原因，它长这样：

注：无论是棋盘格、环形标定板，标定的原理都是一样的，都是提取到角点/中心点的相机像素坐标下坐标，然后利用张正友标定法进行标定。

---

我们以第一张标定板左上角的圆心所处的点作为世界坐标系的原点，并且保存之后重建需要的参数，代码里：

• 
  

% save Rc_1, Tc_1, KK, inv_KK,
% 相机相对于世界坐标系原点的位姿（旋转+平移=外参）
Rc_1 = cameraParams.RotationMatrices(:,:,1);
Rc_1 = Rc_1';    % 转置
Tc_1 = cameraParams.TranslationVectors(1, :);
Tc_1 = Tc_1';
 % 相机内参
KK = cameraParams.IntrinsicMatrix';
save('CamCalibResult.mat', 'Rc_1', 'Tc_1', 'KK');

(2) 投影仪标定

投影仪标定会复杂一点，因为它无法直接看到圆形标定板中心的像素，这个其实也很简单，需要在相机的帮助下，我们投影相移图案，利用相移法获得。

具体地：我们知道，如果要进行单目标定，两块信息是必须已知的，从投影仪角度看：

1. 圆心的世界坐标，已知，同相机
2. 圆心的像素坐标，未知，因为无法直接看到标定板

方法是：

1. 通过投影仪项标定板投射相移图案，相机拍摄，解码获取相位图
2. ==取每个圆心的相位，跟理想的相位图（投影仪）做匹配，得到圆心在投影仪下的像素坐标==

获得这两块信息，直接进行标定即可。这里的代码假设我们已经获取了==圆心在投影仪下的像素坐标==（以后补充），后续操作同相机标定：

• 
  

%% step2: 标定投影仪
disp('开始投影仪标定...');
% 加载投影仪下圆心的像素坐标
load('projector_imagePoints.mat'); % load the projector image points of the centers
% 圆心的世界坐标
worldPoints = generateCheckerboardPoints([num_y+1,num_x+1], dist_circ);
[prjParams,imagesUsed,estimationErrors] = estimateCameraParameters(prjPoints,worldPoints, 'EstimateTangentialDistortion', true);
figure; showReprojectionErrors(prjParams); title('投影仪的重投影误差');
figure; showExtrinsics(prjParams);title('投影仪的外参');

% 保存参数
Rc_1 = prjParams.RotationMatrices(:,:,1);
Rc_1 = Rc_1';
Tc_1 = prjParams.TranslationVectors(1, :);
Tc_1 = Tc_1';
KK = prjParams.IntrinsicMatrix';
save('PrjCalibResult.mat', 'Rc_1', 'Tc_1', 'KK');

2.2 重建

测试对象：

先设置参数：

• 
  

%% step1: input parameters
width = 640; % camera width
height = 480; % camera height
prj_width = 912; % projector width

加载标定参数：

• 
  

%camera: Projection matrix Pc
load('CamCalibResult.mat');
Kc = KK;  % 相机内参
Ac = Kc * [Rc_1, Tc_1]; % 计算相机投影矩阵

%projector: Projection matrix Ap
load('PrjCalibResult.mat');
Kp = KK;  % 投影仪内参
Ap = Kp * [Rc_1, Tc_1];  % 计算投影仪投影矩阵

这里需要说明的是：

Ac = Kc * [Rc_1, Tc_1];
Ap = Kp * [Rp_1, Rc_2];

都是在计算投影矩阵：将世界坐标系下坐标映射到像素坐标系 下。

现在读取测试图片的相位图：

• 
  

% 条纹频率64，也是间距（一个周期由64个像素组成）用于计算绝对相位，另外的频率1、8用于包裹相位展开
f = 64;             % 条纹频率（单个周期条纹的像素个数）
load('up_test_obj.mat');
up_test_obj = up_test_obj / f;  % 将相位归一化到[0, 2pi]之间
figure; imshow(up_test_obj / (2 * pi)); colorbar; title("相位图, freq=" + num2str(f));
figure; mesh(up_test_obj); colorbar; title("相位图, freq=" + num2str(f));

计算投影仪的坐标：

% 计算投影仪坐标
x_p = up_test_obj / (2 * pi) * prj_width;

不同于标准的双目视差重建，我们不需要进行立体校正，直接依据$Eq. \ref{3dr}$ 进行重建：

% 3D重建
Xws = nan(height, width);
Yws = nan(height, width);
Zws = nan(height, width);

for y = 1:height
   for x = 1:width
       if ~isnan(up_test_obj(y, x))
           uc = x - 1;
           vc = y - 1;
           up = (x_p(y, x) - 1);
            % Eq. (32) in the reference paper.
           A = [Ac(1,1) - Ac(3,1) * uc, Ac(1,2) - Ac(3,2) * uc, Ac(1,3) - Ac(3,3) * uc;
                Ac(2,1) - Ac(3,1) * vc, Ac(2,2) - Ac(3,2) * vc, Ac(2,3) - Ac(3,3) * vc;
                Ap(1,1) - Ap(3,1) * up, Ap(1,2) - Ap(3,2) * up, Ap(1,3) - Ap(3,3) * up];
           
           b = [Ac(3,4) * uc - Ac(1,4);
                Ac(3,4) * vc - Ac(2,4);
                Ap(3,4) * up - Ap(1,4)];
     
           XYZ_w = inv(A) * b;
           Xws(y, x) = XYZ_w(1); Yws(y, x) = XYZ_w(2); Zws(y, x) = XYZ_w(3);
       end
   end
end

点云显示：

% 点云显示
xyzPoints(:, 1) = Xws(:);
xyzPoints(:, 2) = Yws(:);
xyzPoints(:, 3) = Zws(:);
ptCloud = pointCloud(xyzPoints);
xlimits = [min(Xws(:)), max(Xws(:))];
ylimits = [min(Yws(:)), max(Yws(:))];
zlimits = ptCloud.ZLimits;
player = pcplayer(xlimits,ylimits,zlimits);
xlabel(player.Axes,'X (mm)');
ylabel(player.Axes,'Y (mm)');
zlabel(player.Axes,'Z (mm)');
view(player,ptCloud);

03 问题

现在还有三个问题：

1. 怎么投影图案，求解绝对相位？
2. 怎么帮助投影仪“看见”标定板？进行投影仪标定？
3. 精度如何提升，比如进行gamma校正？

这些后续再讲~欢迎关注公众号！

04 参考文献

1. GPU-assisted high-resolution, real-time 3-D shape measurement，Song Zhang，OPTICS EXPRESS，2006
2. Calibration of fringe projection profilometry: A comparative review，Feng Shiji，Chao Zuo，Optics and Lasers in Engineering，2021