介绍
单目散斑结构光是一种基于三角测量的非接触光学3D测量方法,具有系统构成简单、成本低、时间分辨率高、近距测量精度高等优势。市面上已经推出了很多基于单目散斑结构光技术的3D相机,知名的有Microsoft Kinect, Intel Realsense, 以及Orbbec Astra系列等,同时也有许多定制化的模组,被应用于机器人避障、人脸解锁、刷脸支付、物流体积测量等诸多场景。
由于系统装配误差、镜头畸变等因素,实际生产出来的单目散斑结构光相机与理想设计一定存在偏差。为了保证深度测量的准确性,需要对这些因素导致的系统误差进行标定。尽管在网上可以搜到很多介绍单目散斑结构光相机的测距原理、立体匹配的博客、论文,但介绍单目散斑结构光相机标定的内容缺较为稀少。原因有两点,其一是标定方案是各个厂商的核心秘密之一,其二是单目散斑结构光的标定在学术上实际也没有被很好解决,缺少类似张正友相机标定法这样简单、准确、灵活、易用的方法。本文介绍笔者提出的一种单目散斑结构光标定方法,具有模型简单、灵活易用的特点。代码和验证数据已上传至https://github.com/zxcn/MSSLSJointCalibration.git。此方法已申请专利,将此方法用于商业目的需获得笔者的授权许可。下面,废话不多说,开始进入正题。
标定模型
单目散斑结构光相机由散斑投影仪、相机、算力芯片组成。投影仪用于向测量场景中投射散斑图像。散斑图案具有全局或者视差匹配范围内的唯一性,便于后续算法匹配。相机用于捕获场景散斑图像。根据场景中物体深度的不同,散斑会发生不同程度的偏移,通过将场景散斑图像与记录的垂直于相机光轴某一距离平面的参考散斑图像进行立体匹配,可计算偏移量的大小,然后根据系统参数转换为深度。算力芯片用于运行立体匹配算法以及滤波算法,负责输出视差、深度、点云等测量结果。
本文用于表示单目散斑结构光系统的模型非常简单。使用计算机视觉中常用的“针孔+畸变”的相机模型描述相机,具体参数有焦距
f
x
,
f
y
f_x,f_y
fx,fy,主点
c
x
,
c
y
c_x,c_y
cx,cy,以及由Brown模型描述的相机畸变,包含径向畸变参数
k
1
,
k
2
,
k
3
k_1,k_2,k_3
k1,k2,k3,以及切向畸变参数
p
1
,
p
2
p_1,p_2
p1,p2。这些参数统称为相机内参。对于相机坐标系中的3D点
P
=
[
x
,
y
,
z
]
T
P=[x,y,z]^T
P=[x,y,z]T,经过透视投影、畸变映射、像素采样后,成为像素坐标中的2D点
U
=
[
u
,
v
]
T
U=[u,v]^T
U=[u,v]T,过程表示如下
x
′
=
x
/
z
y
′
=
y
/
z
r
2
=
x
′
2
+
y
′
2
x
′
′
=
x
′
(
1
+
k
1
r
2
+
k
2
r
4
+
k
3
r
6
)
+
2
p
1
x
′
y
′
+
p
2
(
r
2
+
2
x
′
2
)
y
′
′
=
y
′
(
1
+
k
1
r
2
+
k
2
r
4
+
k
3
r
6
)
+
2
p
2
x
′
y
′
+
p
1
(
r
2
+
2
y
′
2
)
u
=
f
x
x
′
′
+
c
x
v
=
f
y
y
′
′
+
c
y
x'=x/z\\ y'=y/z\\ r^2=x'^2+y'^2\\ x''=x'(1+k_1r^2+k_2r^4+k_3r^6)+2p_1x'y'+p_2(r^2+2x'^2)\\ y''=y'(1+k_1r^2+k_2r^4+k_3r^6)+2p_2x'y'+p_1(r^2+2y'^2)\\ u=f_xx''+c_x\\ v=f_yy''+c_y
x′=x/zy′=y/zr2=x′2+y′2x′′=x′(1+k1r2+k2r4+k3r6)+2p1x′y′+p2(r2+2x′2)y′′=y′(1+k1r2+k2r4+k3r6)+2p2x′y′+p1(r2+2y′2)u=fxx′′+cxv=fyy′′+cy用
U
=
cam
(
P
)
U=\text{cam}(P)
U=cam(P)表示这一过程。不同于条纹结构光常用的用于描述DLP投影仪的逆相机模型,本文将散斑投影仪看成是一个从中心投射光线的点,用其在相机坐标系中的坐标
P
c
=
[
x
c
,
y
c
,
z
c
]
T
P_c=[x_c,y_c,z_c]^T
Pc=[xc,yc,zc]T描述。而对于参考平面,用相机坐标系中的平面方程
A
x
+
B
y
+
z
=
D
Ax+By+z=D
Ax+By+z=D来描述。引入平面方程后,可以使用未知距离和倾角的参考平面采集参考散斑图像,省去了调垂直过程,极大降低了标定实现难度。
一个理想的单目散斑结构光系统如图a所示。投影仪位于相机
x
x
x轴上,参考平面垂直与相机
z
z
z轴。系统预先记录平面上的参考散斑图像,当测量实际场景时,散斑位置发生变化,通过立体匹配找到散斑在图像中位置的变化量,即视差
d
d
d,然后可以通过如下公式计算物体深度
Depth
=
D
1
+
d
⋅
D
f
x
x
c
\text{Depth}=\frac{D}{1+d\cdot\frac{D}{f_xx_c}}
Depth=1+d⋅fxxcDD而实际中,由于系统装配误差、镜头畸变、参考平面倾斜和距离不准确等因素,上述理想假设不再成立,如图b所示,单目散斑结构光系统通常偏离理想状态,需要对系统进行标定。尽管有很多标定方法已经被提出,但这些方法应用上有以下不足:
- 要求参考平面于散斑垂直,需要精细的调节。
- 需要在参考平面上添加标记点,辅助计算参考平面位姿,破坏了参考散斑图像的完整性。
- 需要投射散斑图案的仿真图,代替真实参考平面散斑图做匹配。但仿真散斑图像与真实散斑图像难免存在偏差,会导致额外的匹配误差,同时仿真散斑图案对于个人开发者是难以获得的。
下面将介绍笔者提出的方法是如何灵活地实现标定,并校正这些误差,恢复成如图c中所示理想的状态。
数据采集和处理
实现标定所需的装置主要有圆点标定板和平面物体。圆点标定板的圆需要有比背景更高的反射率,有利于实现高精度的散斑匹配。可适当调大圆的直径与圆间距的比例,让圆点囊括更多的散斑信息。根据散斑清晰成像的距离,选择合适的圆点标定板的尺寸。本文实验中使用散斑结构光系统最近工作距离为600mm,标定板在1m处,整个标定板尺寸为90x70cm,由打印店用A0纸打印,成本约20元,贴在平面物体上作为标定板使用(也可以使用相对平整的墙面)。如果你的散斑结构光相机工作于红外波段,则还需一个卤素投光灯用于提供红外波段的均匀照明(加上三脚架大概100元左右)。数据采集过程如下:
- 摆放圆点标定板为几种不同的姿态,采集这些姿态下均匀照明下的标定板图像以及打开投影仪时散斑照明下的标定板图像。
- 放置平面物体于相机前,采集参考散斑图像。
这样就完成了标定数据的采集。下图给出了实采数据的示意,图a是其中一个姿态的泛光标定板图像,图b是对应姿态的散斑标定板图像,图c是参考散斑图像。
对采集的图像做如下处理:
- 如上图a所示,使用圆点提取算法,检测泛光标定板图像,获得圆点中心坐标 U ~ b \tilde{U}_b U~b。
- 如上图b所示,在对应姿态的散斑标定板图像中,提取 U ~ b \tilde{U}_b U~b附近的散斑区域。
- 如上图c所示,在参考散斑图像中,使用数字图像相关(digital image correlation, DIC)算法,匹配得到 U ~ b \tilde{U}_b U~b的同名点 U ~ r \tilde{U}_r U~r。
DIC算法可以实现亚像素级准度的散斑匹配,为了保证DIC匹配对透视变形和镜头畸变的鲁棒性,建议使用二阶的DIC形函数。数据处理的代码实现可以参考GitHub仓库中的run_circlegrid_detect_and_dic_match.m
。
相机标定
使用泛光标定板图像标定相机。标定板上的圆点坐标为
P
g
=
[
x
g
,
y
g
,
0
]
T
P_g=[x_g,y_g,0]^T
Pg=[xg,yg,0]T,其在相机坐标系中的点坐标
P
b
P_b
Pb与标定板姿态
R
(
ϕ
)
,
t
R(\phi),t
R(ϕ),t有关,表示为
P
b
=
R
(
ϕ
)
P
g
+
t
P_b = R(\phi)P_g+t
Pb=R(ϕ)Pg+t其中,
ϕ
\phi
ϕ为轴角,
R
(
ϕ
)
R(\phi)
R(ϕ)为经罗德里格斯公式变换后的旋转矩阵。
P
b
P_b
Pb经相机成像后,像素坐标为
U
b
=
cam
(
P
b
)
U_b=\text{cam}(P_b)
Ub=cam(Pb)。使用圆点提取算法,处理泛光照明下标定板图案,得到圆点像素坐标
U
~
b
\tilde{U}_b
U~b。通过最小化如下目标函数得到相机主点、焦距、畸变、以及标定板位姿
min
Ω
b
∑
i
,
j
∣
∣
U
~
b
i
,
j
−
U
b
i
,
j
∣
∣
2
\min_{\Omega_b}\sum_{i,j}||\tilde{U}_b^{i,j}-U_b^{i,j}||^2
Ωbmini,j∑∣∣U~bi,j−Ubi,j∣∣2其中,
i
i
i为标定板姿态索引,
j
j
j为圆点索引,优化参数
Ω
b
\Omega_b
Ωb包含
f
x
,
f
y
,
c
x
,
c
y
,
k
1
,
k
2
,
k
3
,
p
1
,
p
2
,
ϕ
i
,
t
i
f_x,f_y,c_x,c_y,k_1,k_2,k_3,p_1,p_2,\phi^i,t^i
fx,fy,cx,cy,k1,k2,k3,p1,p2,ϕi,ti。此步为标准的相机标定过程,优化初值可通过张正友提出的解析方法获得。
投影仪和参考平面标定
使用散斑标定板图像和参考散斑图像标定投影仪和参考平面。如图所示,从
P
c
P_c
Pc发射,经过圆点
P
b
P_b
Pb处的光线与参考平面的交点坐标
P
r
=
[
x
r
,
y
r
,
z
r
]
T
P_r=[x_r,y_r,z_r]^T
Pr=[xr,yr,zr]T表示为
s
=
D
−
[
A
,
B
,
1
]
P
c
[
A
,
B
,
1
]
(
P
b
−
P
c
)
P
r
=
P
c
+
s
(
P
b
−
P
c
)
s=\frac{D-[A,B,1]P_c}{[A,B,1](P_b-P_c)}\\ P_r=P_c+s(P_b-P_c)
s=[A,B,1](Pb−Pc)D−[A,B,1]PcPr=Pc+s(Pb−Pc)其经相机成像后,像素坐标为
U
r
=
cam
(
P
r
)
U_r=\text{cam}(P_r)
Ur=cam(Pr)。利用圆点检测算法获得圆点坐标
U
~
b
\tilde{U}_b
U~b后,使用DIC算法在参考散斑图上匹配圆点处的散斑图案,得到参考图上同名散斑点坐标
U
~
r
\tilde{U}_r
U~r。带入已知的
Ω
b
\Omega_b
Ωb,通过最小化如下目标函数获得参考平面方程和投影仪中心坐标
min
Ω
r
∑
i
,
j
∣
∣
U
~
r
i
,
j
−
U
r
i
,
j
∣
∣
2
\min_{\Omega_r}\sum_{i,j}||\tilde{U}_r^{i,j}-U_r^{i,j}||^2
Ωrmini,j∑∣∣U~ri,j−Uri,j∣∣2其中,优化参数
Ω
r
\Omega_r
Ωr包含
A
,
B
,
D
,
x
c
,
y
c
,
z
c
A,B,D,x_c,y_c,z_c
A,B,D,xc,yc,zc。可以用理论的或者粗略测量的基线长度和参考平面距离分别作为
x
c
x_c
xc和
D
D
D的初值,其余参数初值设为零。
联合标定
最后,进行相机、投影仪、参考平面参数的联合优化,以获得最优的标定参数。目标函数如下
min
Ω
∑
i
,
j
∣
∣
U
~
b
i
,
j
−
U
b
i
,
j
∣
∣
2
+
∑
i
,
j
∣
∣
U
~
r
i
,
j
−
U
r
i
,
j
∣
∣
2
\min_{\Omega}\sum_{i,j}||\tilde{U}_b^{i,j}-U_b^{i,j}||^2+\sum_{i,j}||\tilde{U}_r^{i,j}-U_r^{i,j}||^2
Ωmini,j∑∣∣U~bi,j−Ubi,j∣∣2+i,j∑∣∣U~ri,j−Uri,j∣∣2其中,优化参数
Ω
\Omega
Ω包含
Ω
b
\Omega_b
Ωb和
Ω
r
\Omega_r
Ωr中的所有参数。笔者发现此目标函数在无精确初值的情况下可以很好收敛,实际使用中可以跳过投影仪和参考平面参数的优化过程,使用同样的初值直接进行联合优化。联合标定的代码实现可以参考GitHub仓库中的run_joint_optimization.m
。
散斑图像校正
获得标定参数后,需要对参考散斑图像和场景散斑图像进行校正以使得单目散斑结构光相机成为理想状态。这种状态下,散斑立体匹配仅需在行方向进行,可以直接应用BlockMatch或SGM算法计算视差。校正过程如下:
首先,建立虚拟的极线校正相机坐标系,其原点与原始相机坐标系原点重合,在原始相机坐标系中的旋转矩阵
R
′
R'
R′由Fusiello极线校正方法获得,表示为
R
′
=
[
e
1
,
e
2
,
e
3
]
e
1
=
[
x
c
,
y
c
,
z
c
]
T
/
x
c
2
+
y
c
2
+
z
c
2
e
2
=
[
−
y
c
,
x
c
,
0
]
T
/
x
c
2
+
y
c
2
e
3
=
e
1
×
e
2
R'=[e_1,e_2,e_3]\\e_1=[x_c,y_c,z_c]^T/\sqrt{x_c^2+y_c^2+z_c^2}\\ e_2=[-y_c,x_c,0]^T/\sqrt{x_c^2+y_c^2}\\e_3=e_1{\times}e_2
R′=[e1,e2,e3]e1=[xc,yc,zc]T/xc2+yc2+zc2e2=[−yc,xc,0]T/xc2+yc2e3=e1×e2给定虚拟相机的焦距
f
x
′
,
f
y
′
f_x',f_y'
fx′,fy′,主点
c
x
′
,
c
y
′
c_x',c_y'
cx′,cy′,可以用通用的去畸变和极线校正算法对场景散斑图像进行校正。而对于参考散斑图像校正,需给定垂直于虚拟相机
z
z
z轴虚拟参考平面的距离
D
′
D'
D′。令
[
u
′
,
v
′
]
T
[u',v']^T
[u′,v′]T表示虚拟相机像素坐标,则其观察到的虚拟参考平面上的点在原始相机坐标系中的坐标
P
b
′
P_b'
Pb′为
P
b
′
=
D
′
R
′
[
(
u
′
−
c
x
′
)
/
f
x
′
(
v
′
−
c
y
′
)
/
f
y
′
1
]
P_b'=D'R'\begin{bmatrix}(u'-c_x')/f_x'\\(v'-c_y')/f_y'\\1\end{bmatrix}
Pb′=D′R′
(u′−cx′)/fx′(v′−cy′)/fy′1
由投影仪中心
P
c
P_c
Pc向
P
b
′
P_b'
Pb′投射的光线,将于参考平面交与
P
r
′
P_r'
Pr′,然后被成像为像素坐标
U
r
′
U_r'
Ur′。将
U
r
′
U_r'
Ur′作为插值索引,对采集的参考散斑图像进行重新采样,即可获得虚拟参考散斑图像。注意,给定的
D
′
D'
D′最好接近
D
D
D以保留更多重合的FOV,当参考平面倾角过大时,插值过程应考虑加入照度的校正。场景散斑图像和参考散斑图像校正的代码实现可以参考GitHub仓库中的CorrectSceneImage.m
和CorrectReferenceImage.m
。
实验验证
用仿真和真实实验对提出的标定方法做了验证。代码和数据已上传至https://github.com/zxcn/MSSLSJointCalibration.git。
仿真实验验证
仿真数据在GitHub仓库images/Sim
文件夹下。仿真数据的圆点检测与DIC散斑匹配结果如下:
仿真值:
标定结果:
真实实验验证
使用消费级单目散斑结构光相机进行实验,真实实验数据在GitHub仓库images/Real
文件夹下。首先运行run_circlegrid_detect_and_dic_match.m
进行圆点检测和DIC匹配,然后运行run_joint_optimization.m
进行标定。
设计参数与标定结果对比如下图。标定结果与设计值接近,有效表征了系统偏差。
最后运行run_depth_and_pointcloud_calc.m
进行参考散斑图像和场景散斑图像的校正、视差计算、深度和点云的转换。下图是石膏像的极线校正效果和视差计算结果。通过放大的区域,可以看到极线校正有很高的精度,得到的视差图也十分完整。
与相机输出点云对比如下图,左图是相机输出的点云,右图的是用带入上述标定结果计算的点云。虽然立体匹配算法的不同导致点云细腻程度不同,但整体效果上是相似的,说明了标定方法的有效性和正确性。
展望
- 拓展至多个距离的参考平面标定,进一步提高标定精度,同时通过多距离参考散斑图像的融合,获得更大深度FOV。
- 用于单目条纹结构光标定,好处是无需预知投射的条纹相位,适用于不具有像素化结构的条纹投影装置。