-
如图2所示,条纹投射测量系统主要由投影机和摄像机组成。投影机用于投射正弦条纹图案至被测物体表面,摄像机从另一个角度同步采集物面变形条纹图像。为了研究投影机非线性对测量位相的影响,假设采用相移技术。由计算机生成的第k幅(k = 0, 1, …, K−1)正弦条纹图像可表示为:
$${g_k}(u,v) = \alpha + \beta \cos(2{{\pi }}fu + 2{\rm{\pi }}{k / K})$$ (1) 式中:(u, v)表示投影机像素坐标;K表示相移步数;f表示沿u轴方向的空间频率;正数α和β分别代表条纹的背景和对比度,且满足0≤α±β≤1,以保证条纹灰度值为正值且不会过饱和。
当用投影机投射所生成条纹图像至被测物体表面时,其亮度非线性会导致条纹发生变形。在不进行光度标定的前提下,投影机的非线性模型并不能准确知晓。然而,根据Weierstrass逼近定理[38]可知,任意定义在实数域闭区间的连续函数都可由一个均匀收敛的多项式来近似表示。因此,未知的投影机非线性模型可用多项式精确逼近。例如,Yatabe等[39]利用多项式表示投影机的非线性模型并求出其解析解,从而对非线性引起的位相位差进行补偿。文中用一个N阶多项式来近似表征投影机的输出亮度值,即
$\lambda = \mathop \sum \nolimits_{n = 0}^N {c_n}{g^n}$ ,其中cn为多项式系数,则摄像机采集的变形条纹图像可表示为:$$ {I_k}(x,y) \!=\! R(x,y)\sum\limits_{n = 0}^N {{c_n}{{\{ \alpha \!+\! \beta \cos [\varPhi (x,y) + {{2{\rm{\pi }}k} / K}]\} }^n}} + D(x,y)$$ (2) 式中:(x, y)代表摄像机成像平面像素坐标;R(x, y)表示比例因子,其大小取决于被测物体表面的反射率和法向量;D(x, y)代表与环境光照相关的附加光强。将公式(2)展开并重新整理后可得:
$${I_k}(x,y) = A(x,y) + B(x,y)\sum\limits_{n = 0}^N {{d_n}{{\cos }^n}[\varPhi (x,y) + {{2{\rm{\pi }}k} / K}]} $$ (3) 式中:A(x, y)和B(x, y)分别代表条纹的背景和调制度;dn代表系数,且有d1=1。依据同步探测算法[30]可以从变形条纹图中求出包裹位相为:
$$\psi (x,y) = - \arctan \left[ {\frac{{\displaystyle\sum\nolimits_{k = 0}^{K - 1} {{I_k}(x,y)\sin (2{\rm{\pi }}{k / K})} }}{{\displaystyle\sum\nolimits_{k = 0}^{K - 1} {{I_k}(x,y)\cos (2{\rm{\pi }}{k / K})} }}} \right]$$ (4) 对公式(4)所求得包裹位相进行时域解包裹[40]可得到绝对位相Ψ(x, y)。由于投影机非线性的影响,计算所得Ψ(x, y)包含位相误差。该误差可表示为[41]:
$$\begin{split} \varepsilon (x,y) =\; & \varPsi (x,y) - \varPhi (x,y)\approx \\ & \sum\limits_{j = 1}^{{\rm{int}} [{{(N + 1)} / K}]} {{\xi _j}} \sin [jK\varPhi (x,y)] \end{split} $$ (5) 其中,第j项的系数是:
$${\xi _j} = \sum\limits_{n = 2}^N {{\rho _{j,n}}} $$ (6) 且有
$${\rho _{j,n}} = \left\{ {\begin{aligned} & {\frac{{ - {d_n}jKC_{n + 1}^{(n + 1 - jK)/2}}}{{(n + 1){2^{n - 1}}}}} \quad \begin{aligned} & \{ [1 \le j \le (n + 1)/K]{\text{且}}\\ & (K{\text{和}}j - n{\text{均是奇数}})]\} \\ & {\text{或}}\\ & \{ [1 \le j \le (n + 1)/K]{\text{且}}\\ & (K{\text{是偶数}}, \; j{\text{是奇数}})]\} \end{aligned}\\ & 0 \;\;\quad\quad\quad\quad\quad\quad\quad {{\text{否则}}} \end{aligned}} \right.$$ (7) 从公式(5)可以看出,投影机非线性引起的位相误差是正弦波的线性组合,其大小取决于相移步数K和多项式阶数N。这些误差在位相图上呈波纹状分布,其最低频率是条纹频率的K倍。
为了更加清晰的呈现出投影机非线性对计算位相的影响,图3(a)模拟了一个一维连续位相图,图3(b)为添加载波后的位相图,图3(c)模拟了投影机输入输出的非线性响应曲线。条纹背景α和调制度β的值分别设为0.5和0.4;考虑到实际采集的条纹图像可能包含不均匀的背景和调制度,故设公式(2)中的比例因子R(x, y)为高斯分布且在两端边缘处衰减为中心的50%;由环境光引入的附加光强D(x, y)的值设为0.1。
根据公式(2),最终可以生成三幅模拟的一维条纹强度曲线,其间的相对相移量为2π/3 rad。图3(d)所示为第一幅条纹强度曲线图。依据公式(2)可计算包裹位相,对其进行解包裹运算并移除载波后得到的位相曲线如图3(e)所示,从中可见投影机非线性导致的波纹状位相误差。将图3(e)中的位相值减去图3(a)中的预设位相值,可得位相误差如图3(f)所示,从中可观察到位相误差的频率是变形条纹频率的3倍,这是因为采用了3步相移算法。位相误差的振幅取决于投影机非线性的程度,在本次模拟结果中最大误差约为0.2 rad。
图 3 投影机非线性影响。(a)模拟位相;(b)载波位相;(c)投影机非线性曲线;(d)相移条纹强度曲线;(e)去载波测量位相;(f)位相误差
Figure 3. Effect of projector nonlinearity on the phase measuring results.(a) Simulated phase curve; (b) Simulated phases with a carrier added; (c) Nonlinearity curve of projector; (d) A phase-shifting intensity curve; (e) Measured phases without carrier;(f) Phase errors
-
上一节推导了由投影机非线性引起的位相误差公式,即公式(5)。如果可以求解出该公式中未知的误差系数,就可以从测量位相中剔除该误差,以获得更加精确的位相值。本节将介绍三种不同的自校正方法。第一种方法是通过拟合归一化的条纹强度值与滤波后位相的余弦值之间的依赖关系得到的投影机非线性曲线,并进行迭代运算校正位相误差。第二种和第三种方法都是根据推导的误差公式,分别从单幅位相图和两幅不同频率的位相图中来估计误差系数,以此校正测量位相。
-
本节将介绍一种基于迭代强度拟合的自校正算法[41]。该方法通过从条纹图中直接估计出投影机非线性曲线,进而对位相误差进行校正。
利用公式(4)及位相解包裹技术可求解位相,同时,可以估计出条纹的背景强度:
$$\hat A(x,y){\rm{ = }}\frac{1}{K}\sum\limits_{k = 0}^{K - 1} {{I_k}(x,y)} $$ (8) 以及调制度
$$\begin{split} \hat M(x,y) = & \frac{2}{K}\left\{ {{{\left[ {\sum\limits_{k = 0}^{K - 1} {{I_k}(x,y)\sin \left( {\frac{{2\pi k}}{K}} \right)} } \right]}^2}} +\right. \\ &{\left. {{{\left[ {\sum\limits_{k = 0}^{K - 1} {{I_k}(x,y)\left( {\frac{{2\pi k}}{K}} \right)} } \right]}^2}} \right\}^{\frac{1}{2}}} \end{split} $$ (9) 由于投影机非线性的影响,上述计算结果均存在误差。图4(a)分别用虚线和实线绘出计算所得背景和调制度曲线。从中可见,计算所得调制度也包含有波纹状的误差。
$$\begin{split} {W_k}(x,y) = & \sum\limits_{n = 0}^N {{d_n}} {\cos ^n}\left[ {{{\varPhi (x,y) + 2k\pi }/ K}} \right] \approx \\ &{{\left[ {{I_k}(x,y) - \hat A(x,y)} \right]} / {\hat M(x,y)}} \end{split} $$ (10) 利用公式(10)消除图3(d)中条纹的背景强度和调制度,可得归一化光照强度曲线如图4(b)所示。从中可见,由投影机非线性引起的谐波仍会导致归一化的强度曲线为非正弦曲线。由公式(10)可知,此时的归一化强度曲线是由cos[Φ(x, y)+2πk/K]组成的多项式,如果能得到位相的真值Φ(x, y),那么系数{dn}可以通过公式(10)计算得到。在上一节中,已经获得了一个计算位相值Ψ(x, y),但由于投影机非线性的影响,得到的位相值包含了波纹状的位相误差。若对其作平滑低通滤波处理,则可以得到较为精确的位相值,即
$$\varPhi (x,y) \approx \bar \varPsi (x,y){\rm{ = }}\varPsi (x,y) * G(x,y)$$ (11) 式中:*代表卷积计算;G(x, y)代表空域低通滤波器,例如均值平滑滤波器或者高斯平滑滤波器。选择滤波器的窗口大小应考虑到能覆盖一个波纹周期,即1/K条纹节距(K代表相移的步数)。
对图3(e)中的位相进行高斯低通滤波处理后的结果如图5(a)所示。从图中看出,波纹状的位相误差基本被剔除了,但同时位相曲线的边缘也出现了较大的误差。图5(b)给出了滤波前后的位相差,即图3(e)和图5(a)之间的差值。图5(c)显示了滤波位相(包含载波)的余弦波形。相较于初始计算的位相Ψ(x,y),滤波位相
$\bar \varPsi \left( {x,y} \right)$ 在除边缘区域外的值更接近精确位相值Φ(x, y)。基于这一事实,如果能找出归一化强度信号与滤波后的余弦曲线之间的关联性,则可以估计出公式(10)中的系数{dn}。为了使两者间的关系更清晰,图5(d)绘制了两条曲线所有像素点之间的强度关系,其中,横轴代表图5(c)中的余弦值,纵轴代表图4(b)中的归一化条纹信号值。从图中可以看出,大量的点聚集在一起形成一条曲线,而该曲线实际上就是投影机的非线性曲线。拟合该曲线可表征投影机的非线性规律。需注意,图5(d)中还存在一些偏离曲线的点,这些点的存在可能会导致多项式系数估计误差。因此,拟合前需识别并剔除这些偏离点。在实际测量中,这些误差点的来源主要分成两类:第一类是边界点,源于空域处理的边界效应;第二类是无效点,例如阴影区域或物面梯度非常大的区域。剔除偏离点后,将滤波位相代入到公式(10)中可得线性方程为:
$$\sum\limits_{n = 0}^N {{d_n}} {\cos ^n}\left[ {\bar \psi (x,y) + \frac{{2\pi k}}{K}} \right] = \frac{{{I_k}(x,y) - \hat A(x,y)}}{{\hat M(x,y)}}$$ (12) 公式(12)构成线性方程组,利用最小二乘法可从中求解多项式系数{dn},从而从条纹图中估计出投影机的非线性曲线。文中使用一个4阶多项式对图5(d)中的有效点进行拟合,其结果如图5(e)显示。比较图5(e)和图3(c)可发现,估计得到的非线性曲线与预设曲线在形状上有很大不同,这是由于相移算法对一定条纹谐波不敏感。例如3步相移算法对3的整数倍阶谐波不敏感[32],这些谐波在非线性多项式中所对应项的系数也无法被估计出来。同时,这些项的存在与大小也不会对测量结果产生任何影响。
图 5 基于光强曲线迭代拟合的自校正算法。(a)滤波位相;(b)滤波前后位相差;(c)滤波位相余弦波形;(d)归一化强度与(滤波位相余弦的关系;(e)表示(d)中有效点的多项式拟合结果;(f)校正后位相
Figure 5. Self-correcting method based on iterative intensity curve fitting. (a) Smoothed phase curve ; (b) Phase difference before and after filtering; (c) Cosine of the smoothed phases; (d) Dependence between the normalized intensities and the cosine of smoothed phases; (e) Polynomial fitting to the clustering points in (d); (f) Corrected phases
拟合得到投影机非线性曲线的多项式表示,目的是为了对位相误差进行校正。分别定义背景强度、调制度以及位相的误差为:
$$\left\{ {\begin{aligned} & {\varDelta A(x,y) = \hat A(x,y) - A(x,y)} \\ & {\varDelta M(x,y) = \hat M(x,y) - M(x,y)} \\ & {\varDelta \varPhi (x,y) = \varPsi (x,y) - \varPhi (x,y)} \end{aligned} } \right.$$ (13) 则公式(3)的一阶泰勒展开式为:
$$\begin{split} {I_k}(x,y) = & \hat A(x,y) + \hat M(x,y)\sum\limits_{n = 0}^N {{d_n}} {\cos ^n}\left[ {\varPsi (x,y){\rm{ + }}\frac{{2\pi k}}{K}} \right] - \\ &\frac{{\partial {I_k}(x,y)}}{{\partial \hat A(x,y)}}\varDelta A(x,y) - \frac{{\partial {I_k}(x,y)}}{{\partial \hat M(x,y)}}\varDelta M(s,t) - \\ & \frac{{\partial {I_k}(x,y)}}{{\partial \varPsi (x,y)}}\varDelta \varPhi (x,y)\\[-15pt] \end{split} $$ (14) 其中
$$\left\{\begin{aligned} & \frac{{\partial {I_k}(x,y)}}{{\partial \hat A(x,y)}}= 1 \\ & {\frac{{\partial {I_k}(x,y)}}{{\partial \hat M(x,y)}}= \sum\limits_{n = 0}^N {{d_n}} {{\cos }^n}[\varPsi (x,y) + 2\pi k/K]} \\ & \frac{{\partial {I_k}(x,y)}}{{\partial \varPsi (x,y)}} = - \hat M(x,y)\sin \left[ {\varPsi (x,y)+\frac{{2\pi k}}{K}} \right] \\ & \times \sum\limits_{n = 1}^N {n{d_n}} {\cos ^{n - 1}}\left[ {\varPsi (x,y)+ \frac{{2\pi k}}{K}} \right] \end{aligned} \right.$$ (15) 对于每一像素,k从0变到K−1,共有K个线性方程,而每个线性方程的表达式为:
$$\begin{split} & \frac{{\partial {I_k}(x,y)}}{{\partial \hat A(x,y)}}\varDelta A(x,y) + \frac{{\partial {I_k}(x,y)}}{{\partial \hat M(x,y)}}\varDelta M(x,y) + \frac{{\partial {I_k}(x,y)}}{{\partial \varPsi (x,y)}}\varDelta \varPhi (x,y) \\ & = \hat A(x,y) + \hat M(x,y)\sum\limits_{n = 0}^N {{d_n}} {\cos ^n}\left[ {\varPsi (x,y){\rm{ + }}\frac{{2\pi k}}{K}} \right] - {I_k}(x,y) \end{split} $$ (16) 从中求解出ΔA(x,y)、ΔM(x,y)和ΔΦ(x,y),可分别对背景、调制度及位相进行校正,即有
$$\left\{ {\begin{aligned} & {A(x,y) = \hat A(x,y) - \varDelta A(x,y)} \\ & {M(x,y) = \hat M(x,y) - \varDelta M(x,y)} \\ & {\varPhi (x,y) = \hat \varPsi (x,y) - \varDelta \varPhi (x,y)} \end{aligned} } \right.$$ (17) 根据上述校正过程,并参照图5(e)中的投影机非线性曲线即可对图3(f)中的位相误差进行校正。但考虑到公式(13)中的位相、背景强度以及调制度都包含了误差,而根据这些参数估计出来的多项式系数值是不精确的,因此,需要用迭代的方式来逐渐逼近这些系数的精确值,过程如下:
步骤1:确定迭代过程的初始值。将平滑滤波位相
$ \bar \varPsi \left( {x,y} \right) $ 及通过公式(8)~(9)计算出来的背景$\hat A\left( {x,y} \right)$ 和调制度$\hat M\left( {x,y} \right)$ 分别作为迭代初值Φ(0)(x,y)、A(0) (x,y)、和M(0) (x,y);步骤2:估计投影机非线性的多项式系数。当迭代i次后,位相、背景和调制度估计值分别为Φ(i)(x,y)、A(i) (x,y)、M(i) (x,y),将其代入到公式(12)中可得到一组线性方程组,利用最小二乘法求解得到第(i+1)次的多项式系数估计值dn(i+1);
步骤3:校正位相误差。针对每一像素(x, y),将Φ(i)(x,y)、A(i) (x,y)、M(i) (x,y)以及dn(i+1)代入到公式(16)中,求解出ΔA(x,y)、ΔM(x,y)和ΔΦ(x,y),再根据公式(17)可计算校正后的A(i+1) (x,y)、M(i+1) (x,y)和Φ(i+1)(x,y);
步骤4:重复执行步骤2和步骤3直到算法收敛。
通过上述迭代过程,图3(f)中的位相误差得以校正。经过10次迭代后,校正后的位相如图5(f)所示,其中,波纹状位相误差已经被消除。
-
当仅有一幅位相图时,想要直接从中识别并剔除位相误差较为困难。然而,如上节所述,若采用低通滤波器对计算位相进行平滑滤波处理,并剔除边界效应等产生的粗大误差,就有可能从中估计投影机非线性引起的位相误差函数[42]。
图5(a)的滤波位相在去掉载波位相前进行包裹操作,结果如图6(a)所示。由于图5(b)中的位相差主要是由投影机的非线性引起的,且根据公式(5)可知,它们是一系列正弦波的线性组合,所以结合图5(b)和图6(a),可以绘制位相误差依赖于位相的关系曲线如图6(b)所示,其中纵坐标为图5(b)中的位相差值,横坐标为图6(a)中包裹的滤波位相值。在图6(b)中,大部分的点聚集形成一条中心曲线,且该曲线在2π弧度范围内有3个循环周期,对该中心曲线进行拟合就可以得到位相误差函数曲线。对于图6(b)中偏离中心曲线的粗大误差点,可用异值点检测算法[43]来剔除。当采样点足够多时,采用3-sigma准则能够准确地去除离群点,并保留有效数据。为了保证精度,这里采用迭代的方式去执行这一过程。当获得精确的拟合曲线后,意味着位相误差系数即被确定。需要说明的是,公式(5)中的位相误差公式包含很多项,但随阶数升高,对应振幅快速衰减。因此,在实际操作中,合理保留前J项用于计算误差系数。迭代拟合过程可表示为:
图 6 基于单幅位相图的投影机非线性误差识别与校正。(a)滤波位相包裹结果; (b)位相误差与位相的依赖关系;(c)误差函数迭代拟合曲线;(d)校正位相
Figure 6. Recognizing and removing the projector nonlinearity from a single phase map. (a) Wrapped smoothed phases; (b) Dependence of phase errors on phases; (c) Iterative fitting result of phase error function; (d) Corrected phases
$$\begin{split} & \sum\nolimits_{j = 1}^J {\xi _j^{(i{\rm{ + 1}})}\sin [jK\bar \varPsi (x,y)]} = \varPsi (x,y) - \bar \varPsi (x,y)\\ & (x,y){\text{满足}}\left| {{\varDelta ^{(i)}}(x,y)} \right| \le 3\sigma _\varDelta ^{(i)} \end{split}$$ (18) 其中,括号中的上标代表迭代次数;
${\varDelta ^{\left( i \right)}}\left( {x,y} \right)$ 的计算式为:$$\begin{split} & {\varDelta ^{(i)}}(x,y) = \varPsi (x,y) - \bar \varPsi (x,y) \\ & - \sum\nolimits_{j = 1}^J {\xi _j^{(i{\rm{ + 1}})}\sin [jK\bar \varPsi (x,y)]} \end{split} $$ (19) $\sigma _\varDelta ^{\left( i \right)}$ 代表${\varDelta ^{\left( i \right)}}\left( {x,y} \right)$ 的标准差,即$$\sigma _\varDelta ^{(i)} = \sqrt {\sum\nolimits_{(x,y)} {{{{{[{\varDelta ^{(i)}}(x,y)]}^2}}/ M}} } $$ (20) 式中:
${\varDelta ^{\left( i \right)}}\left( {x,y} \right)$ 为零均值变量;M代表点的数量。第1次迭代时,将Ψ(x, y)和
$\bar \varPsi \left( {x,y} \right)$ 代入公式(18),通过最小二乘算法计算初始误差系数$\xi _j^{\left( 1 \right)}$ 。结合公式(19)和公式(20)可分别求得${\varDelta ^{\left( 1 \right)}}\left( {x,y} \right)$ 和$\sigma _\varDelta ^{\left( 1 \right)}$ ,再根据3-sigma准则剔除异值点。重复上述过程,不断更新${\varDelta ^{\left( i \right)}}\left( {x,y} \right)$ 和$\sigma _\varDelta ^{\left( i \right)}$ 直到算法收敛。图6(c)为迭代10次后得到的拟合曲线,即位相误差函数曲线。其中J=5, 且误差函数的系数为{0.200, –0.020 0, 0.002 6, –0.000 3, –0.000 04},由此表明投影机非线性引起的误差可以通过单幅位相图识别出来。依据误差系数确定位相误差,并进一步校正其影响。图6(d)给出了校正后的位相曲线,从图中观察到由投影机非线性引起的波纹状误差已经几乎被完全消除。
-
条纹投射技术常采用多频条纹,目的是用于计算条纹级次得到绝对位相,但同时还可利用其校正投影机非线性的影响。为此,我们模拟了两组不同频率的条纹图,对应的相移步数均为3步,且相对相移步距为2π/3rad。图7(a)与图7(c)为两幅不同频率条纹曲线图,其测量位相在去除载波后如图7(b)和图7(d)所示。图7(e)中,实线与虚线分别代表高、低频条纹的位相误差。从图7(e)中不难发现,在投影机存在非线性的情况下,不同频率的条纹图位相误差具有相同的振幅。与公式(5)所揭示的规律一致,位相误差振幅只与投影机的非线性、生成条纹的背景和对比度、以及相移步数有关,而与条纹频率无关。
图 7 基于双频条纹图的投影机非线性自校正方法。(a)低频条纹强度曲线;(b)低频测量位相;(c)高频条纹强度曲线;(d)高频测量位相;(e) 低频(虚线)与高频(实线)位相误差;(f)校正位相
Figure 7. Correction of the projector nonlinearity from two-frequency phase maps. (a) Low frequency fringe; (b) Measured low frequency fringe phases; (c) High frequency fringe; (d) Measured high frequency fringe phases; (e) Phase errors in Fig.(b) (dotted) and Fig.(d) (solid); (f) Corrected phases
迭代最小二乘自校正算法[44]可利用双频条纹位相确定位相误差系数。假设图7(a)和图7(c)中的两组相移条纹频率分别为fH和fL,它们对应的解包裹位相分别为ΨH和ΨL,根据公式(5),可得:
$$\left\{ {\begin{aligned} & {{\varPsi _H} - {\varPhi _H} = \sum\limits_{j{\rm{ = }}1}^{{\rm{int}} [{{(N + 1)} / K}]} {{\xi _j}\sin (jK{\varPhi _H})} } \\ & {{\varPsi _L} - {\varPhi _L} = \sum\limits_{j{\rm{ = }}1}^{{\rm{int}} [{{(N + 1)} / K}]} {{\xi _j}\sin (jK{\varPhi _L})} } \end{aligned}} \right.$$ (21) 式中:ΦH和ΦL是位相精确值,分别对应条纹频率fH和fL。根据ΦH和ΦL之间的关系
${{{\varPhi _H}} / {{\varPhi _L}}}{\rm{ = }}{{{f_H}} / {{f_L}}}$ ,公式(21)进一步写为:$$\left\{ {\begin{aligned} & {{\varPsi _H} - {\varPhi _H} = \sum\limits_{j = 1}^{{\rm{int}} [{{(N + 1)} / K}]} {{\xi _j}\sin (jK{\varPhi _H})} } \\ & {{\varPsi _L} - \frac{{{f_L}}}{{{f_H}}}{\varPhi _H} = \sum\limits_{j = 1}^{{\rm{int}} [{{(N + 1)} / K}]} {{\xi _j}\sin \left( {\frac{{{f_L}}}{{{f_H}}}jK{\varPhi _H}} \right)} } \end{aligned}} \right.$$ (22) 高阶项的振幅会快速衰减,因此保留前J项较显著的误差项。由此,公式(22)形成非线性方程组,其未知数为ξj和ΦH。设迭代初值为
$\varPhi _H^{\left( 0 \right)} = {\varPsi _H}$ 。该方程组可通过下列不动点迭代计算:$$ \begin{split} \varPhi _H^{(i + 1)} = & \frac{1}{{1 + \dfrac{{{f_L}}}{{{f_H}}}}}\left\{ {\left[ {{\varPsi _H} - \sum\limits_{j = 1}^J {\xi _j^{(i)}\sin (jK\varPhi _H^{(i)})} } \right]}+ \right. \\ & \left. {\left[ {{\varPsi _L} - \sum\limits_{j =1}^J {\xi _j^{(i)}\sin \left( {\frac{{{f_L}}}{{{f_H}}}jK\varPhi _H^{(i)}} \right)} } \right]} \right\} \end{split} $$ (23) 继续图7的数值模拟过程,取J=5,经过10次迭代,得到的误差系数为{–0.199 6, 0.020 28, –0.003 0, 0.000 9, –0.000 2},校正后的位相如图7(f)所示。从图中看出,由投影机引起的波纹状位相误差已经基本被消除。
-
在双频条纹投射技术中,在对精度要求不高的前提下,还可采用基于统计学的算法计算误差函数[45]。仅保留公式(21)第一项,可得:
$$\left\{ {\begin{aligned} & {{\varPsi _H}(x,y) - {\varPhi _H}(x,y) = {\xi _1}\sin (K{\varPhi _H}(x,y))} \\ & {{\varPsi _L}(x,y) - {\varPhi _L}(x,y) = {\xi _1}\sin (K{\varPhi _L}(x,y))} \end{aligned}} \right.$$ (24) 据此可推导
$$\begin{split} & {\varPsi _H}\left( {x,y} \right) - \frac{{{f_H}}}{{{f_L}}}{\varPsi _L}\left( {x,y} \right) =\\ & {\xi _1}\sin \left( {K{\varPhi _H}\left( {x,y} \right)} \right) - {\xi _1}\frac{{{f_H}}}{{{f_L}}}\sin \left( {K{\varPhi _L}\left( {x,y} \right)} \right) \end{split} $$ (25) 求解上式两侧均方根即可直接计算误差系数ξ1的绝对值为:
$$\begin{split} \left| {{\xi _1}} \right| = & \frac{{\sqrt {\frac{{\rm{1}}}{Q}\sum\limits_{(x,y)} {{{\left[ {{\varPsi _H}(x,y) - \frac{{{f_H}}}{{{f_L}}}{\varPsi _L}(x,y)} \right]}^2}} } }}{{\sqrt {\frac{{\rm{1}}}{Q}\sum\limits_{(x,y)} {{{\left[ {\sin \left( {K{\varPhi _H}(x,y)} \right) - \frac{{{f_H}}}{{{f_L}}}\sin \left( {K{\varPhi _L}(x,y)} \right)} \right]}^2}} } }}\approx \\ & \frac{{\sqrt {2\sum\limits_{(x,y)} {{{\left[ {{\varPsi _H}(x,y) - \frac{{{f_H}}}{{{f_L}}}{\varPsi _L}(x,y)} \right]}^2}} } }}{{\sqrt {Q\left[ {{{\left( {\frac{{{f_H}}}{{{f_L}}}} \right)}^2} + 1} \right]} }}\\[-32pt] \end{split} $$ (26) 式中:Q代表像素点的数量。为确定ξ1的代数符号,定义误差公式为:
$$\begin{split} e(x,y) = & {\varPsi _H}(x,y) - {\xi _1}\sin \left( {K{\varPsi _H}(x,y)} \right) -\\ & \left[ {\frac{{{f_H}}}{{{f_L}}}{\varPsi _L} - \frac{{{f_H}}}{{{f_L}}}{\xi _1}\sin \left( {K{\varPsi _L}(x,y)} \right)} \right] \end{split} $$ (27) 分别尝试ξ1取正负符号代入公式(27),正确的符号可使
$\sum\nolimits_{({\rm{x}},{\rm{y}})} {\left| {e(x,y)} \right|} $ 取得较小值。$\left| {{\xi _1}} \right|$ 的估计精度与图像中的条纹数量有关。条纹越多,$\left| {{\xi _1}} \right|$ 的精度越高。在估计ξ1的值并确定其代数符号之后,可借助公式(23)所示的不动点迭代校正位相误差,其中M=1。 -
本节通过实验来比较上述不同自校正算法的有效性。测量系统如图2所示,包括了一个DLP投影机(PHILIPS PPX4010, 854 pixel×480 pixel),和一个CCD摄像机(AVT Stingray F-125B)。采用3步相移算法计算位相,利用多频条纹投射方法以时域位相解包裹算法求解绝对位相。图8(a)为一幅条纹图,大小为512 pixel×512 pixel,其测量位相在去除载波成分后如图8(b)所示。从中可见波纹状的位相误差,且频率为对应条纹的3倍,是典型的投影机非线性误差。
分别使用了上述的自校正方法,对位相图进行非线性误差校正,并将其结果转化为深度图。对应的深度图如图9所示。其中,图9(a)显示的是未经非线性校正的深度图,图9(b)显示的是采用光度标定主动补偿投影机非线性后的测量结果,图9(c)显示第2.1节所述强度非线性曲线迭代拟合算法所得结果,图9(d)为第2.2节所述基于单幅位相图位相误差估计的自校正方法测量结果,图9(e)为利用第2.3.2节所述基于统计学的自校正方法从双频位相图中重建的深度图,图8(f)则是利用第2.3.1节所述最小二乘迭代自校正方法从双频位相图中重建的深度图。由图中可见,主动补偿方法及各种自校正方法均可有效地抑制投影机非线性的影响。相比较而言,主动补偿方式和统计学自校正方法计算所得深度图中尚残留较多谐波误差。
图 9 不同投影机校正算法计算所得深度图。(a)未校正;(b)光度标定方法;(c)基于光强曲线迭代拟合的自校正算法[41];(d)基于单幅位相图位相误差估计的自校正方法[42];(e)基于统计学的自校正方法[45];(f)基于双频位相的最小二乘迭代拟合自校正方法[44]
Figure 9. Depth maps of different correcting methods for projector nonlinearities. (a) Projector nonlinearity not corrected; (b) Photometric calibration;(c) Self-correcting methods based on iterative fitting of intensity curve[41]; (d) Phase error estimation from a single phase map[42]; (e) Self-correcting method based on statistics[45]; (f) Iterative least squares fitting method based on two-frequency phase maps[44]
为了评估前述自校正方法的精度与效率,图10比较了各种方法的均方根(RMS)误差及计算时间。图中的均方根误差由被测物体背景平板上的竖直截面深度测量数值求得。从图10中可见,各类校正方法在不同程度上均可有效地抑制投影机非线性的影响。其中,采用光度标定方法,因其无需在数据处理阶段校正位相误差,耗时较少。但是,在测量之前,标定投影机非线性曲线却需要耗费大量时间。在各类自标定方法中,基于迭代光强拟合的自校正算法的数据处理效率相对较低,这是因为该方法涉及迭代计算条纹的位相、背景与调制度,耗时较多。最小二乘迭代拟合的双频位相图自校正方法具有相对较高的精度,但其数据处理仍涉及迭代拟合操作,也需要一定的计算时间。基于统计学的双频位相图自校正方法计算效率很高,但其只能估计并抑制最低一项谐波的影响,精度较其他自校正方法为低。表1中以传统的基于标定的非线性补偿方法和相移算法为参照,对文中所介绍的各种自校正方法的主要性能进行了比较。其中,关于计算复杂性,各种方法的计算时间严重依赖于计算机系统的性能。随着计算硬件水平的提升,计算时间已逐渐不再是一个制约自校正算法应用的主要因素。
图 10 不同投影机非线性误差校正算法误差及计算时间比较。(a)~(f)所用方法对应于图9(a)~(f)
Figure 10. Comparison among different compensating methods with the methods in (a)-(f) corresponding to those in Fig. 9(a)-(f)
表 1 投影仪非线性补偿方法性能比较
Table 1. Performance comparisons of projector nonlinearity correcting methods
Methods Prior calibration required Number of the required fringe patterns Computational complexity Insensitivity to time-variance of the projector nonlinearity Expected accuracy Calibration-based methods (e.g. LUT and
phase-error function)Yes Small Low High Middle Increasing the number of phase shifts with phase-shifting technique No Large Low Low High Iterative least squares fitting to the intensity curve based on single-frequency fringe patterns[41] No Small High Low High Phase error estimation from a calculated
phase map[42]No Small Middle Low High Phase error estimation from two-frequency phase maps using iterative least squares fitting method[44] No Small Middle Low High Phase error estimation from two-frequency phase maps using fringe statistics[45] No Small Low Low Low
Progress in self-correcting methods of projector nonlinearity for fringe projection profilometry
-
摘要: 在条纹投射技术中,投影机光强非线性是影响测量精度的关键因素之一。投影机非线性会在条纹信号中引入高阶谐波,从而导致位相测量结果中出现波纹状误差。投影机非线性的自适应校正方法,也即自校正方法,可以从测量数据中直接估计投影机输入输出光强曲线或位相误差函数,从而避免了繁琐的前标定过程,并因此在应用中具备很强的适应性。文中拟对投影机非线性自校正算法的研究进展作一个系统性的概述。这些方法中,其一,是从采集的条纹图像中利用迭代拟合算法直接估计投影机的非线性曲线,并依据其校正位相误差;其二,是从单幅测量位相图中识别并移除由非线性引起的位相误差;其三,是利用两幅不同频率的测量位相图估计误差系数,并补偿其影响。实际测量结果表明:上述自校正方法,在无需标定数据条件下,可以有效地解决投影机非线性误差问题,有助于提高条纹投射技术的测量精度。Abstract: In fringe projection profilometry, the luminance nonlinearity of the projector has been recognized as one of the most crucial factors decreasing the measurement accuracy. It induces the ripple-like artifacts on the measured phase map. The self-adaptive correcting algorithms, i.e., self-correcting algorithms, allow us to suppress the effect of the projector nonlinearity without a prior calibration for the projector intensities or phase errors. This paper introduces the research progress in the self-correcting algorithms. Among them, the first algorithm is to determine a nonlinear curve representing the projector nonlinearity, directly from the captured fringe patterns, thus correcting the phase errors using this curve. The second one is to recognize and remove the nonlinearity-induced errors, directly from a calculated phase map. With the last one, error function coefficients are estimated from a couple of phase maps having different frequencies. Measurement results demonstrate these self-correcting algorithms to be effective in suppressing influences of the projector nonlinearity in the absence of any calibration information.
-
Key words:
- projector nonlinearity /
- fringe projection /
- phase measurement
-
图 3 投影机非线性影响。(a)模拟位相;(b)载波位相;(c)投影机非线性曲线;(d)相移条纹强度曲线;(e)去载波测量位相;(f)位相误差
Figure 3. Effect of projector nonlinearity on the phase measuring results.(a) Simulated phase curve; (b) Simulated phases with a carrier added; (c) Nonlinearity curve of projector; (d) A phase-shifting intensity curve; (e) Measured phases without carrier;(f) Phase errors
图 5 基于光强曲线迭代拟合的自校正算法。(a)滤波位相;(b)滤波前后位相差;(c)滤波位相余弦波形;(d)归一化强度与(滤波位相余弦的关系;(e)表示(d)中有效点的多项式拟合结果;(f)校正后位相
Figure 5. Self-correcting method based on iterative intensity curve fitting. (a) Smoothed phase curve ; (b) Phase difference before and after filtering; (c) Cosine of the smoothed phases; (d) Dependence between the normalized intensities and the cosine of smoothed phases; (e) Polynomial fitting to the clustering points in (d); (f) Corrected phases
图 6 基于单幅位相图的投影机非线性误差识别与校正。(a)滤波位相包裹结果; (b)位相误差与位相的依赖关系;(c)误差函数迭代拟合曲线;(d)校正位相
Figure 6. Recognizing and removing the projector nonlinearity from a single phase map. (a) Wrapped smoothed phases; (b) Dependence of phase errors on phases; (c) Iterative fitting result of phase error function; (d) Corrected phases
图 7 基于双频条纹图的投影机非线性自校正方法。(a)低频条纹强度曲线;(b)低频测量位相;(c)高频条纹强度曲线;(d)高频测量位相;(e) 低频(虚线)与高频(实线)位相误差;(f)校正位相
Figure 7. Correction of the projector nonlinearity from two-frequency phase maps. (a) Low frequency fringe; (b) Measured low frequency fringe phases; (c) High frequency fringe; (d) Measured high frequency fringe phases; (e) Phase errors in Fig.(b) (dotted) and Fig.(d) (solid); (f) Corrected phases
图 9 不同投影机校正算法计算所得深度图。(a)未校正;(b)光度标定方法;(c)基于光强曲线迭代拟合的自校正算法[41];(d)基于单幅位相图位相误差估计的自校正方法[42];(e)基于统计学的自校正方法[45];(f)基于双频位相的最小二乘迭代拟合自校正方法[44]
Figure 9. Depth maps of different correcting methods for projector nonlinearities. (a) Projector nonlinearity not corrected; (b) Photometric calibration;(c) Self-correcting methods based on iterative fitting of intensity curve[41]; (d) Phase error estimation from a single phase map[42]; (e) Self-correcting method based on statistics[45]; (f) Iterative least squares fitting method based on two-frequency phase maps[44]
图 10 不同投影机非线性误差校正算法误差及计算时间比较。(a)~(f)所用方法对应于图9(a)~(f)
Figure 10. Comparison among different compensating methods with the methods in (a)-(f) corresponding to those in Fig. 9(a)-(f)
表 1 投影仪非线性补偿方法性能比较
Table 1. Performance comparisons of projector nonlinearity correcting methods
Methods Prior calibration required Number of the required fringe patterns Computational complexity Insensitivity to time-variance of the projector nonlinearity Expected accuracy Calibration-based methods (e.g. LUT and
phase-error function)Yes Small Low High Middle Increasing the number of phase shifts with phase-shifting technique No Large Low Low High Iterative least squares fitting to the intensity curve based on single-frequency fringe patterns[41] No Small High Low High Phase error estimation from a calculated
phase map[42]No Small Middle Low High Phase error estimation from two-frequency phase maps using iterative least squares fitting method[44] No Small Middle Low High Phase error estimation from two-frequency phase maps using fringe statistics[45] No Small Low Low Low -
[1] Srinivasan V, Liu H C, Halioua M. Automated phase-measuring profilometry of 3-D diffuse object [J]. Applied Optics, 1984, 23(18): 3105−3108. doi: 10.1364/AO.23.003105 [2] Gorthi S S, Rastogi P. Fringe projection techniques: whither we are? [J]. Optics and Lasers in Engineering, 2010, 48(2): 133−140. doi: 10.1016/j.optlaseng.2009.09.001 [3] Zhang S. Recent progresses on real-time 3d shape measurement using digital fringe projection techniques [J]. Optics and Lasers in Engineering, 2010, 48(2): 149−158. doi: 10.1016/j.optlaseng.2009.03.008 [4] Anna T, Dubey S K, Shakher C, et al. Sinusoidal fringe projection system based on compat and non-mechanical scanning low-coherence Michelson interferometer for three-dimensional shape measurement [J]. Optics Communications, 2009, 282(7): 1237−1242. doi: 10.1016/j.optcom.2008.11.080 [5] Berberova N, Stoykova E, Kang H, et al. SLM-based sinusoidal fringe projection under coherent illumination [J]. Optics Communications, 2013, 304: 116−122. doi: 10.1016/j.optcom.2013.04.034 [6] Lei S, Zhang S. Flexible 3-D shape measurement using projector defocusing [J]. Optics Letters, 2009, 34(20): 3080−3082. doi: 10.1364/OL.34.003080 [7] Ayubi G A, Ayubi J A, Maritino J M, et al. Pulse-width modulation in defocused three-dimensional fringe projection [J]. Optics Letters, 2010, 35(21): 3682−3684. doi: 10.1364/OL.35.003682 [8] Wang Y, Zhang S. Optimal pulse width modulation for sinusoidal fringe generation with projector defocusing [J]. Optics Letters, 2010, 35(24): 4121−4123. doi: 10.1364/OL.35.004121 [9] Ayubi G A, Di Martino J M, Alonso J R, et al. Three-dimensional profiling with binary fringes using phase-shifting interferometry algorithms [J]. Applied Optics, 2011, 50(2): 147−154. doi: 10.1364/AO.50.000147 [10] Ayubi G A, Di Martino J M, Alonso J R, et al. Color encoding of binary fringes for gamma correction in 3-D profiling [J]. Optics Letters, 2012, 37(8): 1325−1327. doi: 10.1364/OL.37.001325 [11] Zhu J, Zhou P, Su X, et al. Accurate and fast 3D surface measurement with temporal-spatial binary encoding structured illumination [J]. Optics Express, 2016, 24(25): 28549−28560. doi: 10.1364/OE.24.028549 [12] Coggrave C R, Huntley J M. High-speed surface profilometer based on a spatial light modulator and pipeline image processor [J]. Optical Engineering, 1999, 38(9): 1573−1581. doi: 10.1117/1.602209 [13] Kakunai S, Sakamoto T, Iwata K. Profile measurement taken with liquid-crystal grating [J]. Applied Optics, 1999, 38(13): 2824−2828. doi: 10.1364/AO.38.002824 [14] Dai M, Yang F, He X. Single-shot color fringe projection for three-dimensional shape measurement of objects with discontinuities [J]. Applied Optics, 2012, 51(12): 2062−2069. doi: 10.1364/AO.51.002062 [15] Wang Z, Nguyen D A, Barnes J C. Some practical considerations in fringe projection profilometry [J]. Optics and Lasers in Engineering, 2010, 48(2): 218−225. doi: 10.1016/j.optlaseng.2009.06.005 [16] Hoang T, Pan B, Nguyen D, et al. Generic gamma correction for accuracy enhancement in fringe-projection profilometry [J]. Optics Letters, 2010, 35(12): 1992−1994. doi: 10.1364/OL.35.001992 [17] Zhang X, Zhu L, Li Y, et al. Generic nonsinusoidal fringe model and gamma calibration in phase measuring profilometry [J]. Journal of Optical Society of America A, 2012, 29(6): 1047−1058. doi: 10.1364/JOSAA.29.001047 [18] Liu K, Wang Y, Lau D L, et al. Gamma model and its analysis for phase measuring profilometry [J]. Journal of Optical Society of America A, 2010, 27(3): 553−562. doi: 10.1364/JOSAA.27.000553 [19] Babaei A, Saadatseresht M, Kofman J. Exponential fringe pattern projection approach to gamma-independent phase computation without calibration for gamma nonlinearity in 3D optical metrology [J]. Optics Express, 2017, 25(21): 24927−24938. doi: 10.1364/OE.25.024927 [20] Gai S, Da F. A novel fringe adaptation method for digital projector [J]. Optics and Lasers in Engineering, 2011, 49(4): 547−552. doi: 10.1016/j.optlaseng.2010.12.004 [21] Ye X, Cheng H-B, Wu H-Y, et al. Gamma correction for three-dimensional object measurement by phase measuring profilometry [J]. Optik, 2015, 126(24): 5534−5538. doi: 10.1016/j.ijleo.2015.09.028 [22] Li Z, Li Y. Gamma-distorted fringe image modeling and accurate gamma correction for fast phase measuring profilometry [J]. Optics Letters, 2011, 36(2): 154−156. doi: 10.1364/OL.36.000154 [23] Ma S, Quan C, Zhu R, et al. A fast and accurate gamma correction based on Fourier spectrum analysis for digital fringe projection profilometry [J]. Optics Communications, 2012, 285(5): 533−538. doi: 10.1016/j.optcom.2011.11.041 [24] Zhang S. Yau S-T. Generic nonsinusoidal phase error correction for three-dimensional shape measurement using a digital video projector [J]. Applied Optics, 2007, 46(1): 36−43. doi: 10.1364/AO.46.000036 [25] Li Z, Shi Y, Wang C. Accurate calibration method for a structured light system [J]. Optical Engineering, 2008, 47(5): 053604. doi: 10.1117/1.2931517 [26] Li Z-W, Shi Y-S, Wang C-J, et al. Complex object 3D measurement based on phase-shifting and a neural network [J]. Optics Communications, 2009, 282(14): 2699−2706. doi: 10.1016/j.optcom.2009.04.055 [27] Pan B, Kemao Q, Huang L, et al. Phase error analysis and compensation for nonsinusoidal waveforms in phase-shifting digital fringe projection profilometry [J]. Optics Letters, 2009, 34(4): 416−418. doi: 10.1364/OL.34.000416 [28] Hu Y, Xi J, Chicharo J F, et al. Inverse Function Analysis Method for Fringe Pattern Profilometry [J]. IEEE Transactions on Instrumentation and Measurement, 2009, 58(9): 3305−3314. doi: 10.1109/TIM.2009.2022382 [29] Zuo C, Feng S, Huang L, et al. Phase shifting algorithms for fringe projection profilometry: A review [J]. Optics and Lasers in Engineering, 2018, 109: 23−59. doi: 10.1016/j.optlaseng.2018.04.019 [30] Bruning J H, Herriott D R, Gallagher J E, et al. Digital wavefront measuring interferometer for testing optical surfaces and lenses [J]. Applied Optics, 1974, 13(11): 2693−2073. doi: 10.1364/AO.13.002693 [31] Stetson K A, Brohinsky W R. Electrooptic holography and its application to hologram interferometry [J]. Applied Optics, 1985, 24(21): 3631−3637. doi: 10.1364/AO.24.003631 [32] Guo H, Chen M. Fourier analysis of the sampling characterizing of the phase-shifting algorithm [C]//SPIE, 2003, 5180: 437−444. [33] Liu Z. Zibley P C, Zhang S. Motion-induced error compensation for phase shifting profilometry [J]. Optics Express, 2018, 26(10): 12632−12637. doi: 10.1364/OE.26.012632 [34] Lu Y, Zhang R, Guo H. Correction of illumination fluctuations in phase-shifting technique by use of fringe histograms [J]. Applied Optics, 2016, 55(1): 184−197. doi: 10.1364/AO.55.000184 [35] Ma S, Zhu R, Quan C, et al. Blind phase error suppression for color-encoded digital fringe projection profilometry [J]. Optics Communications, 2012, 285(7): 1662−1668. doi: 10.1016/j.optcom.2011.12.027 [36] Guo H, He H, Chen M. Gamma correction for digital fringe projection profilometry [J]. Applied Optics, 2004, 43(14): 2906−2914. doi: 10.1364/AO.43.002906 [37] Guo H, Zhao Z. Nonlinearity correction in digital fringe projection profilometry by using histogram matching technique [C]//SPIE, 2007, 6616: 66162I. [38] Banaschewski B. On the Weierstrass-Stone approximation theorem [J]. Fundamenta Mathematicae, 1957, 44(3): 249−252. doi: 10.4064/fm-44-3-249-252 [39] Yatabe K, Ishikawa K, Oikawa Y. Compensation of fringe distortion for phase-shifting 3D shape measurement by inverse map estimation [J]. Applied Optics, 2016, 55(22): 6017−6024. doi: 10.1364/AO.55.006017 [40] Zuo C, Huang L, Zhang M, et al. Temporal phase unwrapping algorithms for fringe projection profilometry: A comparative review [J]. Optics and Lasers in Engineering, 2016, (54)85: 84−103. doi: 10.1016/j.optlaseng.2016.04.022 [41] Lü F, Xing S, Guo H. Self-correction of projector nonlinearity in phase-shifting fringe projection profilometry [J]. Applied Optics, 2017, 56(25): 7204−7216. doi: 10.1364/AO.56.007204 [42] Xing S, Guo H. Directly recognizing and removing the projector nonliearity errors from a phase map in phase-shifting fringe projection profilometry [J]. Optics Communications, 2019, 435(11): 212−220. doi: 10.1016/j.optcom.2018.11.045 [43] Hodge V J, Austin J. A survey of outlier detection methodologies [J]. Artificial Intelligence Review, 2004, 22(2): 85−126. doi: 10.1023/B:AIRE.0000045502.10941.a9 [44] Xing S, Guo H. Correction of projector nonlinearity in multi-frequency phase-shifting fringe projection profilometry [J]. Optics Express, 2018, 26(13): 16277−16291. doi: 10.1364/OE.26.016277 [45] Xing S, Guo H. A statistic method of project nonlinearity correction in multifrequency phase-shifting fringe projection profilometry [C]//SPIE, 2018, 10827: 1082711.
计量
- 文章访问数: 1400
- HTML全文浏览量: 1105
- 被引次数: 0