1 引言

三维(3D)重建作为非接触式测量方法,广泛应用于增材制造、医疗、文物恢复、3D测量等领域^[1]。相位展开作为3D重建的重要环节,会影响重建模型的精度。因此,研究一种鲁棒性强、准确性和效率高的相位展开算法对物体的3D重建具有重要意义。

在条纹图分析中,相位φ(x,y)通常由包含反正切函数的表达式获得,其中,(x,y)为对应点的坐标。反正切函数的返回值在[-π,π]之间^[2-3],其包裹的相位φ(x,y)为锯齿形,通过相位展开可得到物体高度的真实映射。为了提高相位展开算法的抗噪性,人们提出了多种算法^[4-5],如全局算法、区域算法、路径跟踪算法。刘稳等^[6]提出了基于取整最小二乘的全局相位展开方法,用非迭代方法处理相位跳变部分的梯度,降低了误差点对信息传播的影响。Costantini^[7]结合路径分割法和最小范数法的优点,将相位展开的最小化问题转化为求解网络中的最小费用流问题,鲁棒性较强,但计算量大。钟何平等^[8]在进行相位展开时,先计算相位质量图,并将包裹相位分为高、低质量区域,然后利用质量引导与最小不连续相结合的相位展开方法将相位展开和区域分割结果发送至主线程,最后主线程用最小不连续算法优化对相位展开的不同分块区域边界。Ding等^[9]针对合成孔径雷达干涉(InSAR)相位展开中的相关误差,提出了一种用己知InSAR轨道模型模拟相位修正前干涉相位图的相位展开算法。区域算法是对图像进行区域分割,其鲁棒性没有全局算法强,但计算量较低;廖薇等^[10]为了快速准确地对含噪包裹位相图进行相位展开,仿真对比了几种不同的路径跟踪算法,但这些方法都存在处理效率低或对噪声处理效果不佳的问题。

针对上述问题,本文提出了用Canny边缘检测算法检测真实边缘,获取真实边缘后,用二阶差分函数计算需要展开像素处的左右、上下、对角像素,得出边缘可靠性。然后将相邻边缘可靠性值较高的像素归为一组,采用逐区域扩展方法沿可靠性高的边缘进行展开。最后将区域展开连接成全局范围内的相位展开。实验结果表明,该算法可以提高物体廓形边缘的识别能力,且计算量较小,是一种鲁棒性和准确性较高的相位展开方法。

2 三维重建原理

结构光3D视觉是在光学三角法基础上实现的,数字光学处理(DLP)投影仪作为光源将一定频率的结构光投射到物体表面,先得到同一频率的4张正弦光栅图像,其相位分别为0、π/2、π、3π/2。被投射的光栅在物体表面上会发生畸变,得到调制后的光栅图像。然后利用CCD相机采集光栅的二维畸变图像并传输到个人笔记本电脑(PC)中。最后将畸变的光栅图像与参考光栅图像的相位进行对比,解相位得到物体的三维廓形,本系统的工作原理如图1所示。

图 1. 结构光3D重建系统的原理图

Fig. 1. Schematic diagram of the structured light 3D reconstruction system

下载图片 查看所有图片

相移法根据光波干涉原理将满足测量要求的正弦分布光栅用DLP投影系统投射到被测物体表面,受物面高度分布的调制,条纹会发生畸变。用CCD相机采集的二维(2D)畸变条纹图的灰度值I_N可表示为

IN=a(x,y)+b(x,y)cos[2πfx+Δφ(x,y)],(1)

式中,a(x,y)、b(x,y)分别为相机接收到的光强和物面背景光强,Δφ(x,y)为隐含物体表面高度信息的相位函数,f为光栅的频率,N(N≥3)为畸变光栅的图像数量。

3D重建的基本原理是三角测量原理^[11],如图2所示,投影仪向物体投射正弦光栅后,可通过相机采集因物体高度变化而发生畸变的调制图像。其中,o_c、o_p分别为相机和投影仪的光轴中心,o_c-x_cy_cz_c、o_p-x_py_pz_p分别为相机坐标系和投影仪坐标系。

图 2. 三角测量法的原理

Fig. 2. Schematic diagram of the triangulation method

下载图片 查看所有图片

相邻条纹图像之间的相移量为2π/N,采集移相前后的畸变条纹图像,用N步相移算法求出物面上的相位分布^[12-13],第n张畸变光栅的相位函数可表示为

Δφ(x,y)=arctan∑n=0N-1In(x,y)sin2nπN∑n=0N-1In(x,y)cos2nπN,(2)

四步相移的相位函数为

Δφ(x,y)=arctanI2(x,y)-I4(x,y)I1(x,y)-I3(x,y)。(3)

3D重建的关键是准确高效地获取物体的深度信息,相位包裹是反三角函数将深度信息截断至-π到π内。此时包裹相位实际是按投射的正弦光栅周期分布的,因此需要对包裹相位进行相位展开,获取连续相位并通过计算相位Δφ(x,y)与高度h的映射关系,得到物体真实的高度为

h=lΔφ(x,y)2πfd+Δφ(x,y),(4)

式中,d为相机镜头轴心点到投影仪镜头轴心点的距离,l为投影仪镜头轴心点到参考平面的距离。

3 相位展开原理

相位展开是给每个包裹的像元加上2π的整数倍,得到一个连续的相位。但实际上相位不连续是由环境噪声、物体形状的不规则性、采样丢失等原因导致,相邻包裹中相邻像素间的相位差超过π,即为不连续现象。设任意像元(m,n)处的包裹相位为φ_m,n,展开相位为ϕ_m,n,则^[14]

ϕm,n=φm,n+2πkm,n,(5)

式中,k_m,n为整数,是从一点到另外一点的积分,与路径无关。

当相邻像素间的相位差大于π时,说明像元间存在不连续现象,而实际情况中环境噪声、欠采样、物体不规则都会导致相邻像素间的相位差发生剧烈变化,即相邻像素为离散数据。当相机采集图像的灰度发生剧烈变化时,相邻像素间的相位差可能会产生一个局部极值,而二阶差分对图像灰度变化的观测和边缘定位效果较好。

二阶差分算法能很好地检测相位图中可能存在的不一致,度量相位函数的凹度/凸度^[15]Q_i,j可表示为

Qi,j=[h2(i,j)+v2(i,j)]1/2,(6)

式中,h(i,j)、v(i,j)分别为水平像素值和垂直像素值。在二维包裹相位图中,假设在尺寸为3 pixel×3 pixel的窗口中,像素的可靠性由像素周围的8个相邻像素决定。根据可靠性排序算法进行相位展开,对可靠性高的区域进行展开,采用的二阶差分函数D(i,j)可表示为

D(i,j)=[H2(i,j)+V2(i,j)+D12(i,j)+D22(i,j)]1/2H(i,j)=γ{[ϕ(i-1,j)-ϕ(i,j)]-[ϕ(i,j)-ϕ(i+1,j)]}V(i,j)=γ{[ϕ(i,j-1)-ϕ(i,j)]-[ϕ(i,j)-ϕ(i,j+1)]}D1(i,j)=γ{[ϕ(i-1,j-1)-ϕ(i,j)]-[ϕ(i,j)-ϕ(i+1,j+1)]}D2(i,j)=γ{[ϕ(i-1,j+1)-ϕ(i,j)]-[ϕ(i,j)-ϕ(i+1,j-1)]},(7)

式中,H(i,j)、V(i,j)、D₁(i,j)、D₂(i,j)分别为水平、垂直、两个对角像素间的二阶差分值,γ为相位展开操作,ϕ(i,j)为像素值,令

y=γ(x),(8)γ(x)=(xmod2π)-π,(9)

式中,mod为取模运算,γ(x)为相位展开后的函数值。

坐标为(i,j)的像素可靠性函数R(i,j)为

R(i,j)=1D(i,j)。(10)

从(10)式可以看出,二阶差分函数值D(i,j)越小,相位展开的可靠性就越高。尺寸为3 pixel×3 pixel窗口内像素的可靠性值如图3(a)所示,其路径展开方式为图3(a)中相邻像素的和,如图3(b)所示。其中,矩形框内的数值为沿该路径展开得到的可靠性值,矩形框所在的边线为展开路径。

图 3. 相邻像素间的可靠性值。(a)单个像素的可靠性值;(b)展开路径的可靠性值

Fig. 3. Reliability value between adjacent pixels. (a) Reliability value of a single pixel; (b) reliability value of the expanded path

下载图片 查看所有图片

图4为路径展开的方式,图4(a)中的可靠性值最大为2.3,即方框内数值为0的区域,先对该区域的相邻像素进行展开(第一次);随后是可靠性值较大的2.2、2.1区域,分别对应方框内数值为1、2的区域,空白区域表示尺寸为3 pixel×3 pixel窗口内未展开的像素区域。图4(b)中2.1为当前较大的可靠性值,与已展开的0区域进行合并展开(第二次);图4(c)中2.0为当前最大的可靠性值,方框内为1的区域沿其中一个数值为2.0的路径展开,方框内为0的区域沿数值为2.0的路径展开(第三次);最后合并方框内为2的区域,得到全局展开图,如图4(d)所示,方框内的数值全为0(第四次)。

图 4. 可靠性排序的过程。(a)第一次; (b)第二次; (c)第三次; (d)第四次

Fig. 4. Process of reliability sequencing. (a) First time; (b) second time; (c) third time; (d) fourth time

下载图片 查看所有图片

可靠性排序算法的核心是对整张图像进行逐像素循环扫描,而二阶差分算法对边缘比较敏感,在相位展开后可得到稳定的物体廓形,为后期重建物体提供了基础。

图5分别为Laplacian、LoG和Canny算子对耳蜗模型的边缘检测结果,可以发现,图5(a)中Laplacian算子的边缘检测效果最差,无法得到耳蜗模型的真实边缘;图5(b)中的LoG算子能较好地检测出耳蜗模型的廓形;而图5(c)中的Canny算子能得到更准确的耳蜗模型边缘,因此,实验采用Canny边算子进行边缘检测。

图 5. 三种算子的边缘检测结果。(a) Laplacian算子;(b) LoG算子;(c) Canny算子

Fig. 5. Edge detection results of the three operators. (a) Laplacian operator; (b) LoG operator; (c) Canny operator

下载图片 查看所有图片

Canny算子在检测物体边缘时不会丢失重要的边缘信息,且能避免产生虚假边缘。为了提高3D重建的准确度和效率,将Canny算子与可靠性排序算法相结合,大大提高了物体廓形的识别率,避免噪声等因素在物体重构过程中产生的拉丝现象^[15]。Canny算子的检测步骤如下。

1) 用高斯滤波器进行降噪处理

3D重建过程中,环境噪声会造成相位不连续。因此,相位展开过程中需对图像进行降噪,去除孤立噪点。采用高斯滤波器对数据进行卷积运算,可以减少部分类型的噪声,提升边缘检测效果,(2k+1)×(2k+1)的滤波器方程^[16]可表示为

Hij=12πσ2×exp-i-(k+1)2+j-(k+1)22σ2,(11)

式中,1≤i,j≤(2k+1),k为正整数,σ为方差。

2) 计算像素的梯度

设像元(m,n)的大小为f(x,y),梯度为G[f(x,y)],则像素的梯度可表示为

G[f(x,y)]={[Δxf(x,y)]2+[Δyf(x,y)]2}1/2。(12)

式中,Δ_x为函数f(x,y)对x的偏导,Δ_y为函数f(x,y)对y的偏导。相机采集的图像边缘是指向各个方向的,因此,用Canny算子检测出水平、垂直和对角的边缘,返回Δ_x和Δ_y后再确定像元(m,n)的梯度。

3) 非极大值抑制和边缘连接处理

Canny算子选取的高斯标准差为0.6,边缘检测算法的低阈值为50,高阈值为150。图6(a)为实际采集的被物体调制后的图像,由于形貌复杂导致椭圆区域的信息丢失;图6(b)和图6(c)分别为Canny算子检测过程中经非极大值抑制(NMS)处理以及边缘连接后的图像。可以看出,在相位展开过程中进行边缘检测时,需要预先消除图像中的虚假边缘,检测出耳蜗模型的实际边缘,以降低后续相位展开过程的错误率、提高检测效率。图7为相位展开的具体步骤。

图 6. 边缘检测处理的过程。(a)丢失信息的耳蜗调制图;(b)经NMS处理后的图像;(c)经边缘处理后的图像

Fig. 6. Process of edge detection processing. (a) Ear modulation map with missing information; (b) image obtained after NMS processing; (c) image obtained after edge processing

下载图片 查看所有图片

图 7. 相位展开流程图

Fig. 7. Flow chart of phase unwrapping

下载图片 查看所有图片

4 仿真与实验

4.1 相位展开算法的误差分析

为了验证本算法进行相位展开的准确性,对函数f(x)=sin(πx)+cos(πx)的相位进行展开,结果如图8(a)所示。对函数f(x)添加系数为0.1的随机噪声,用本算法和传统Itoh算法根据某一点与其相邻点的相位差进行相位展开,然后对该函数相位展开后得到的一维数据进行拟合,结果如图8(b)所示。其中,实线为理想拟合曲线,符号“+”为本算法相位展开后得到的曲线,通过cftool(数据拟合工具箱)得到本算法与理想相位展开拟合曲线的均方根误差(RMSE)为0.0121;实心点为传统算法相位展开后的数据点,其与理想相位展开数据点的RMSE为0.0562,这表明相比传统算法,本算法的解包裹效果有明显改进。

图 8. 相位展开的误差分析。(a)相位函数的3D描述;(b)不同算法的相位展开误差

Fig. 8. Error analysis of phase unwrapping. (a) 3D description of the phase function; (b) phase unwrapping errors of different algorithms

下载图片 查看所有图片

为了进一步验证本算法的准确性,通过四组仿真实验验证相位展开的效果,仿真环境包括添加噪声和不连续或信息丢失的情形,结果如图9所示。可以发现,在图9(a)中的理想环境下,本算法能对原始展开相位进行包裹后精确地重展开;在图9(b)中添加系数为0.1的随机噪声后,本算法仍然可以精确地实现相位展开,这表明本算法的抗噪能力较强。在图9(c)中信息丢失和图9(d)中添加系数为0.1的随机噪声和存在不连续或信息丢失的情况下,本算法也能很好地进行包裹和相位展开,这表明本算法具有较高的鲁棒性和准确性。表1为四种仿真环境下相位展开所需的时间,实验环境:电脑为戴尔笔记本电脑,CPU为i5-8300H。

图 9. 本算法在不同环境下的相位展开结果。(a)理想情况;(b)随机噪声;(c)信息丢失;(d)随机噪声和信息丢失

Fig. 9. Phase unwrapping results of our algorithm in different environments. (a) Ideal situation; (b) random noise; (c) information loss; (d) random noise and information loss

下载图片 查看所有图片

4.2 相位展开算法的实际测试结果

实验的外部设备主要包括DLP投影仪和CCD相机,上位机为64位的Windows10笔记本电脑,搭建的系统如图10所示。首先,通过自行编写的软件产生正弦结构光栅,光栅的分辨率为1920 pixel×1080 pixel,频率为48 Hz,正弦光栅为I=128+127sin(2πf+φ)。然后,用DLP投影仪将光栅投射到耳蜗模型表面,用CCD相机采集调制光栅,结果如图11(b)所示。其中,椭圆标记区域的信息因物体自身形貌复杂而缺失。采用边缘检测及可靠性排序算法进行3D相位展开,能解决耳蜗模型在实际数据采集过程中存在环境噪声以及自身形貌遮挡导致的相位展开混乱问题。基于边缘检测及可靠性排序算法进行相位展开,得到的廓形显示效果较好,这表明本算法的鲁棒性较高。

图 10. 搭建的系统实物图

Fig. 10. Physical map of the built system

下载图片 查看所有图片

在相同环境下采集一组耳蜗模型的调制光栅,分别用枝切法、最小二乘法、以及本算法进行仿真对比。图11(a)为对相机采集的光栅进行相位展开得到的结果,图11(b)为相机采集的耳蜗模型调制图像,通过提取图像的感兴趣区域(ROI),得到大小为402 pixel×316 pixel的窗口,最后采用四步相移法验证相位展开的效果。

图 11. 光栅展开后的图像和耳蜗模型调制图。(a)光栅相位的展开图像;(b)耳蜗模型的调制光栅图

Fig. 11. Image and ear model modulation after the raster is expanded. (a) Unfolded image of grating phase; (b) modulated grating image of ear model

下载图片 查看所有图片

将采集的图像用Goldstein枝切法进行3D重建,采集时处于自然光环境,由于耳蜗模型自身的形貌复杂,未采集到图11(b)中的椭圆区域。同时,相机采集数据的过程中也会存在噪声、信息丢失等情况导致相位展开不连续,残差点会密集分布,枝切线的设置会产生错误,导致相位图中产生“孤岛”现象^[17],造成整个相位的展开错误。图12为枝切法对采集图像的仿真,可以看出,耳蜗模型的廓形没有恢复,这表明Goldstein枝切法在实际重建过程中并不适用。

图 12. Goldstein枝切法得到的结果。(a)耳蜗模型的2D图;(b)相位展开后的耳蜗模型;(c)耳蜗模型的真实高度

Fig. 12. Results obtained by the Goldstein branch cutting method. (a) 2D image of the ear model; (b) ear model after phase expansion; (c) true height of the ear model

下载图片 查看所有图片

图13分别为最小二乘法和本算法得到的包裹相位图和相位展开图,由于实际采集耳蜗模型调制图时存在许多噪声,可以看出,使用最小二乘法进行相位展开时存在耳蜗模型廓形拉丝现象,原因是最小二乘法无法限制误差在空间中的传播,而本算法得到的耳蜗模型廓形更清晰。

图 13. 两种算法的包裹相位和相位展开图。(a)~(b)最小二乘法;(c)~(d)本算法

Fig. 13. Wrapping phase and phase unfolding diagrams of the two algorithms. (a)-(b) Least square method; (c)-(d) our algorithm

下载图片 查看所有图片

用Matlab处理采集的数据,图14为最小二乘法和本算法进行相位展开后,对耳蜗模型的一面进行重建的仿真图。可以发现,图14(a)的二维图中存在拉丝现象,原因是模型实际采集过程中存在误差,且最小二乘法无法解决误差在空间传播的问题^[18],导致重建的图14(c)也出现了断层现象。而本算法采用边缘检测及可靠性排序算法,可以很好地解决或缩小误差的传播,得到的耳蜗模型三维廓形更清晰。

图 14. 两种算法的相位展开结果。(a)最小二乘法(2D);(b)本算法(2D);(c)最小二乘法(3D);(d)本算法(3D)

Fig. 14. Phase unwrapping results of the two algorithms. (a) Least squares algorithm (2D); (b) our algorithm (2D); (c) least squares algorithm (3D); (d) our algorithm (3D)

下载图片 查看所有图片

用本算法结合四步相移法求取耳蜗模型的廓形信息,其中,相机光心到参考平面的距离为L=26 cm,投影仪光心到相机光心连线的距离D=10 cm,投影光栅的频率f=1/48(具体频率根据物体的大小调整)。相机和相机间的最佳夹角在20°~30°之间。图15为本算法得到耳蜗模型一面的高度信息,其重建的映射高度为25.8 cm。

图 15. 3D重建的耳蜗模型

Fig. 15. Cochlear model obtained by 3D reconstruction

下载图片 查看所有图片

5 结论

针对耳蜗模型,首先,采用Canny边缘检测对采集图像的进行掩码输出,获取图像的真实边缘,对ROI区域内的数据进行运算。其次,通过二阶差分,给出可靠性函数评判标准,依据可靠性值的高低,使每个区域沿着可靠性值高的边缘展开,根据可靠性排序进行相位展开,最后实现全局相位展开。实验结果表明,本算法可以解决枝切法中因噪声密集产生的“孤岛”现象,同时也避免了最小二乘法不能解决误差在空间中传播的问题。且本算法对于信息缺失的相位也能正确地进行相位展开,在三维重建过程中的鲁棒性和准确性较高。