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

Original 曹博 计算机视觉工坊 6/17

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

注：本文的理论主要来自于参考文献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 boardnum_y = 9; %number of circles in the y direction of the calibration boarddist_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 parameterswidth = 640; % camera widthheight = 480; % camera heightprj_width = 912; % projector width

加载标定参数：

%camera: Projection matrix Pcload('CamCalibResult.mat');Kc = KK;  % 相机内参Ac = Kc * [Rc_1, Tc_1]; % 计算相机投影矩阵%projector: Projection matrix Apload('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   endend

点云显示：

% 点云显示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);

现在还有三个问题：

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

备注：作者也是我们「3D视觉从入门到精通」特邀嘉宾：一个超干货的3D视觉学习社区本文仅做学术分享，如有侵权，请联系删文。下载1在「计算机视觉工坊」公众号后台回复：深度学习，即可下载深度学习算法、3D深度学习、深度学习框架、目标检测、GAN等相关内容近30本pdf书籍。下载2在「计算机视觉工坊」公众号后台回复：计算机视觉，即可下载计算机视觉相关17本pdf书籍，包含计算机视觉算法、Python视觉实战、Opencv3.0学习等。下载3在「计算机视觉工坊」公众号后台回复：SLAM，即可下载独家SLAM相关视频课程，包含视觉SLAM、激光SLAM精品课程。

重磅！计算机视觉工坊-学习交流群已成立

扫码添加小助手微信，可申请加入3D视觉工坊-学术论文写作与投稿 微信交流群，旨在交流顶会、顶刊、SCI、EI等写作与投稿事宜。

同时也可申请加入我们的细分方向交流群，目前主要有ORB-SLAM系列源码学习、3D视觉、CV&深度学习、SLAM、三维重建、点云后处理、自动驾驶、CV入门、三维测量、VR/AR、3D人脸识别、医疗影像、缺陷检测、行人重识别、目标跟踪、视觉产品落地、视觉竞赛、车牌识别、硬件选型、深度估计、学术交流、求职交流等微信群，请扫描下面微信号加群，备注：”研究方向+学校/公司+昵称“，例如：”3D视觉 + 上海交大 + 静静“。请按照格式备注，否则不予通过。添加成功后会根据研究方向邀请进去相关微信群。原创投稿也请联系。

▲长按加微信群或投稿

▲长按关注公众号

3D视觉从入门到精通知识星球：针对3D视觉领域的视频课程（三维重建系列、三维点云系列、结构光系列、手眼标定、相机标定、orb-slam3等视频课程）、知识点汇总、入门进阶学习路线、最新paper分享、疑问解答五个方面进行深耕，更有各类大厂的算法工程人员进行技术指导。与此同时，星球将联合知名企业发布3D视觉相关算法开发岗位以及项目对接信息，打造成集技术与就业为一体的铁杆粉丝聚集区，近2000星球成员为创造更好的AI世界共同进步，知识星球入口：

学习3D视觉核心技术，扫描查看介绍，3天内无条件退款

 圈里有高质量教程资料、可答疑解惑、助你高效解决问题觉得有用，麻烦给个赞和在看~