CN103267496A - 一种基于小波变换的改进窗口傅里叶三维测量法 - Google Patents

Abstract

本发明公开了一种基于小波变换的改进窗口傅里叶三维测量法，其主要目的在于求解出精确的条纹图像的相位分布，并由相位分布得到物体的三维形貌信息。其实现步骤为：将一幅黑白正弦调制条纹图像投影到被测物体上并用CCD进行采集；对采集到的变形条纹图像逐行进行小波变换，提取小波变换脊，然后计算并去除求取脊过程中相位二阶导所带来的误差，最后得到精确的窗口大小矩阵。将精确窗口矩阵代入窗口傅里叶变换，经过滤波等步骤求得变形条纹图的相对相位信息。建立条纹图质量图，然后采用洪水算法相位展开，得到条纹图像的绝对相位分布。根据相位高度转换公式由绝对相位分布获取被测物体的三维信息。

Description

一种基于小波变换的改进窗口傅里叶三维测量法

技术领域

本发明属于三维信息重构的技术领域，基于调制光栅投影，结合小波变换和窗口傅里叶变换法，将条纹投影轮廓术用于物体的三维测量。

背景技术

基于光学投影的三维测量，广泛应用于产品检测和加工控制、医疗领域、文物保护领域、航空航天领域、文化领域等。在众多的三维测量技术中，光学三维测量技术以其非接触测量、实时性较好等特点成为占主导地位的三维测量技术。

光学三维测量技术是以现代光学为基础，融光电子学，计算机图像处理，图形学，信号处理等科学为一体的现代测量技术。它将光学图像作为信息检测和存储的实体，从图像中获得三维测量所需信息。相对于其它三维信息测量方法，光学测量技术具有高速度、高精度的特点，从上世纪60年代以来，得到了广泛的研究和应用。基于编码光栅投影的三维测量技术最具代表性，应用最为广泛。

基于光栅投影的三维测量就是将光栅图样投影到被测物表面，由摄像机获取变形的光栅图像，并由相位与高度的关系来确定出被测物体相对参考平面的高度信息。当光栅投影到物体表面上时，周期性光栅的相位就受到物体表面高度轮廓的调制，形成变形光栅，变形光栅中带有物体的三维信息。准确获取变形光栅图像的相位分布在整个三维测量过程中起着关键作用。

通常基于光栅投影的三位测量法分为多帧处理相移法和单帧处理变换测量法。傅里叶变换轮廓术(FTP)是一种主动式单帧测量法，最早由Takeda等提出，众多学者对此作了深入的研究。由此开启了傅里叶变换法在三维测量领域的新的应用。但是由于傅里叶变换(FT)是全局变换，不具有局部分析的能力，故局部的错误无法单独分析解决。基于此，具有局部分析优势的窗口傅里叶变换(WFT)被引入三维测量领域。局部信号分析使得信号处理结果更为精确，而其中关键就是前期的窗口选择。窗口越大，频率分辨率越好而空间分辨率越差；窗口越小，频率分辨率越差而空间分辨率越好。因此窗口的选择对于三维测量的精度有着重大影响。根据海森伯格测不准原理，我们不能给出一个完全准确的窗口大小，但是可以不断的优化。因此如何选取自适应窗口成为一个难点，于是小波脊法等被引入三维测量，对于窗口大小的选择给出了详尽的分析，对窗口傅里叶测量法的精度和速度都有提高。但是由于窗口傅里叶处理高分辨率的图像比较缓慢，因此在实时性有要求的情况下，有必要尽量减少窗口大小的计算时间。

由于傅里叶变换法得到的相位值是介于0-2π的，所以需要进行相位展开得到绝对相位值。但是因为窗口傅里叶变换在物体高度跳变和陡峭区域一般脊法求取其窗口大小带有很大误差，导致相位求解的过程中会出现误差，此误差在相位展开时会影响接下来待处理点的相位展开精度，所以会导致相位展开时产生连续错误，即拉线现象。直接展开相位的扫描线方法速度较快，但鲁棒性较弱，极易产生拉线现象。双频法通过投射含有两种频率分量的正弦光栅，经傅里叶变换滤波，得到两种频率下的相位值，用低频下得到的相位值修正高频下得到的相位值，从而修正高频下得到的相位不连续处的错误值。但由于双频法也是用傅里叶变换求解初始相位，所以也存在着两种频率分量之间的频谱混叠问题。

发明内容

发明目的：针对小波变换在求取窗口傅里叶窗口大小时速度缓慢以及高度跳变和陡峭区域求取窗口大小误差大的问题，本发明的目的在于不影响窗口大小求取准确性的前提下，解决跳变区域求取相位不准确的问题，提供一种基于小波变换的改进窗口傅里叶三维测量法。

技术方案：一种基于小波变换的改进窗口傅里叶三维测量法，包括如下步骤：

步骤1)，将黑白条纹图像投影到被测物体表面，使用CCD摄像头对被测物体表面进行拍摄，得到一幅高度为c、宽度为r的变形条纹图像g(x，y)：

g(x，y)＝A(x，y)+B(x，y)cos [2πf₀x+φ(x，y)]    (1)式

其中，A(x，y)是背景光强分布，B(x，y)是物体表面反射率，f₀是正弦条纹频率，φ(x，y)是待求的相对相位分布图，(x，y)表示变形条纹图像的二维坐标；

步骤2)，对所述变形条纹图g(x，y)逐行做小波变换，得到变形条纹图像每一点的近似尺度因子，具体过程如下：

步骤2.1)，将y视为常数y₁，采用一维小波变换对所述变形条纹图像g(x，y)的第y₁行进行处理，处理过程为： $W_{y_1}(a_1,b)={\∫}_{-\∞}^{+\∞}g(x,y_1){M^{*}}_{a,b}\left(x\right)\mathrm{dx}$  (2)式

其中，y₁表示某一行的行号；a是尺度因子，取值范围为10到90，每隔0.2取一个值；b是平移因子，取值范围为1到条纹的宽度r，每隔1取一个值，单位为像素；经一维小波变换获得的

是一个400行r列的二维复数矩阵，称第y₁行小波变换矩阵，a₁是矩阵

中元素的行标号；其中，

M^* _a，b(x)为M_a，b(x)的共轭函数；根据式(3)的小波函数M(x)得到M^* _a，b(x)表达式如式(4)：

$M\left(x\right)=\frac{1}{{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp \left(2\mathrm{\πi}{f}_{c}x\right)\exp \left(\frac{{-x}^2}{{2f}_b^2}\right)$   (3)式

${M^{*}}_{a,b}\left(x\right)=\frac{1}{\sqrt{a}{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp (-2\mathrm{\πi}{f}_c\frac{x-b}{a})\exp [-\frac{{(x-b)}^2}{2f_b^2{a}^2}];$   (4)式

其中，f_b是小波函数的带宽，f_c是小波函数的中心频率，i为单位复数，x为自变量；

步骤2.2)，求取条纹图像g(x，y)的近似尺度因子分布图f_a1(x，y₁)：

计算二维复数矩阵

的对应的模矩阵

搜索模矩阵第x列中值最大的元素，并将模矩阵

的第x列中值最大元素的行标号赋值给a_max，则条纹图像的近似尺度因子分布图f_a1(x，y₁)在坐标(x，y₁)处数值a_rx为a_rx＝10+0.2×a_max；

遍历条纹图像所有坐标点，求得条纹图像的近似尺度因子分布图f_a1(x，y)；

步骤3)，对经所述步骤2.1)小波变换后的变形条纹图去除小波变换脊误差，得到变形条纹图像每一点的准确尺度因子，具体过程如下：

步骤3.1)，去除变形条纹图像每一个坐标点的脊误差，处理过程为：

W(x，y₁)＝W₀+W₁+W₂+ε₀    (5)式

其中W(x，y₁)为变形条纹图中任意一坐标点的小波变换简化形式，其中ε₀表示由变形条纹图中每一点的相位二阶导所带来的误差，其表达式为：

  (6)式

其中，t为中间计算变量，

表示相位二阶导，对近似尺度因子逐行求导求出

求出

最后求出对应尺度因子a的ε₀，得到ε₀为400行的一维复数数组在坐标(x，y₁)的坐标点上的误差数组ε₀(x，y₁)；W₀，W₁，W₂分别为小波变换计算过程中的简化表达式，形式如下：

$W_0(a,b)=A(x,y)\frac{a\sqrt{\mathrm{\π}}}{{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp \left({-2f}_b^2{f}_c^2{\mathrm{\π}}^2\right)$   (7)式

  (8)式

  (9)式

其中，

为相位、

表示相位导数，W(x，y₁)中逐点减去误差ε₀(x，y₁)，求出精确小波变换脊的值W_ε(x，y₁)＝W₀+W₁+W₂，求出每一行精确小波变换脊矩阵W_400×r，矩阵W_400×r为400行r列的复数矩阵；

步骤3.2)，求取条纹图像的精确尺度因子分布图f_a2(x，y)：

求出W_400×r的对应的模矩阵C_400×r(a₁，b)，搜索模矩阵C_400×r(a₁，b)第x列中值最大的元素，并将模矩阵C_400×r(a₁，b)的第x列中值最大元素的行标号赋值给a_amax，则条纹图像的近似尺度因子分布图f_a2(x，y₁)在坐标(x，y₁)处数值a_acr为：a_acr＝10+0.2×a_amax；

遍历条纹图像所有坐标点，求得条纹图像的精确尺度因子分布图f_a2(x，y)；

步骤4)，对变形条纹图逐行做窗口傅里叶变换，求出变形条纹图相对相位分布图，具体过程如下：

步骤4.1)，将y视为常数，采用一维窗口傅立叶变换对变形条纹图像g(x，y)的每行进行处理，一维窗口傅立叶变换过程为：

$\mathrm{WF}(b,\mathrm{\ξ})=\underset{-\∞}{\overset{+\∞}{\∫}}g\left(x\right)W_{\mathrm{\δ}}(x-b)\exp (-\mathrm{j\ξx})\mathrm{dx}$   (10)式

其频域表达式为：

$\mathrm{WF}(f_s,y)=\underset{n=0}{\overset{n=+\∞}{\mathrm{\Σ}}}{P}_n(f_s-{\mathrm{nf}}_0,y)$   (11)式

其中，g(x)为一行条纹图像；WF(b，ξ)表示一维窗口傅立叶变换，ξ表示频域计算因子，δ表示一维窗口傅立叶的窗口尺度因子，δ取值为一维窗口的位置(x-b，y)对应的精确尺度因子分布图f_a2(x，y)相应位置上的点的值，W_δ(x-b)表示窗口函数，表达式为：

$W_{\mathrm{\δ}}(x-b)=\frac{{\mathrm{\δ}}^{-1}}{\sqrt{2\mathrm{\π}}}\exp [-\frac{{(b-x)}^2}{{2\mathrm{\δ}}^2}]$   (12)式

其中，n表示阶次，依次取值为0，1，2，……，无穷；P_n(f_s-nf₀，y)表示任一点傅里叶变换后对应的n阶频谱；f_s表示频域的变量；

步骤4.2)，对傅里叶变换后的频谱滤波并提取相位信息，具体过程如下：

对P_n(f_s-nf₀，y)进行滤波并提取出一阶频谱P₁(f_s-f₀，y)，再对P₁(f_s-f₀，y)进行逆傅里叶变换，得到含有相位信息的B(x，y)exp[iφ(x，y)]，计算B(x，y)exp[iφ(x，y)]的角度值即可得到含有物体高度信息的变形条纹图相对相位值φ(x，y)，得到的相位值是介于0-2π之间；

遍历条纹图像所有坐标点，得到变形条纹图的相对相位分布图φ_A(x，y)；

步骤5)，建立条纹图像质量图，采用洪水算法展开相对相位分布图φ_A(x，y)，得到实际相位图

具体过程如下：

步骤5.1)，利用相对相位分布图中的相位梯度建立质量图，质量图可以按照以下公式计算：

$\mathrm{\Δ}(x,y)=\sqrt{W_{\mathrm{rap}}^2\{{\mathrm{\φ}}_A(x+1,y)-{\mathrm{\φ}}_A(x,y)\}+W_{\mathrm{rap}}^2\{{\mathrm{\φ}}_A(x,y+1)-{\mathrm{\φ}}_A(x,y)\}}$   (13)式

其中W_rap{}是包裹函数，当值大于2π或小于-2π时将其减去或加上2π；

步骤5.2)，在相对相位分布图中央找到质量值最高的坐标点，选择该点作为相位展开的起始点，将该点放入一个空的栈中；

步骤5.3)，判断栈是否为空；如果是，则相位展开过程结束，得到展开相位结果

并进入步骤6)；如果否，弹出栈顶的点，展开该点的四邻点中没有处理的坐标点，并将这些未处理的点入栈；

步骤5.4)，按照质量图中的质量值将栈中所有点排序，质量最高的点放在栈顶，转到步骤5.3)继续处理；

步骤6)，读取最终的展开相位结果根据经典光栅投影的从展开相位结果

到物体高度h(x，y)的转换公式，最终求得测量物体的三维信息；其中转换公式h(x，y)为：

  (14)式

其中，l、d是测量***的几何参数，l是投影仪到测量平面的距离，d是CCD摄像头到投影仪的距离；

表示相位变化量，

为展开相位结果，

为初始相位结果，ω₀为投影光栅的角频率。

有益效果：与现有技术相比，本发明具有以下优点：首先，本发明相比于较为常用的基于小波变换的三维测量法，利用小波变换的算法对频率敏感的特殊性提高求取窗口大小的速度并且降低跳变剧烈区域窗口求取误差。通过求取并去除相位二阶导对小波变换的脊频求取误差，能够解决在跳变剧烈区域解相位误差偏大的问题，可以克服物体具有较大跳变和陡峭斜面带来的解相位误差，得到更加精确的相位，并且提高窗口傅里叶的抗噪性，达到了较好的测量精度和较强的鲁棒性；相比于傅里叶测量法，避免了全局变换带来的无法局部相位计算的问题，使得计算的相位信息更加准确；相比于相移法，本发明只需投影一副黑白调制条纹图像，可以实现动态测量；相比于基于斯尔维托克变换的的改进窗口傅里叶三维测量法，由于小波变换的多分辨率特性，能够将信号分解进行处理，能够提高相位的准确率。其次，本发明在相对相位求取中得到较为精确的相位信息，使得使用普通的全局洪水解包裹法就能得到准确的实际相位信息。最后，本发明在处理过程中表现出较强的鲁棒性，有着很好的实际应用的能力。综上，本发明能够准确地得到被测物体的三维高度信息，具有较好的鲁棒性。

附图说明

图1是本发明的流程图；

图2是步骤5中采用洪水算法展开相位的具体过程流程图；

图3是CCD采集到的以摩托车护板为被测物体的变形条纹图像；

图4是变形条纹图中某一点小波变换后的值取模值的曲线图，横坐标代表尺度因子的倒数，纵坐标代表小波变换后的模值；

图5是变形条纹图中某一点小波变换后的值去除脊误差的模值曲线图，横坐标代表尺度因子的倒数，纵坐标代表小波变换后去误差后的模值；

图6是获得的精确窗口大小分布图；

图7是窗口傅立叶变换后提取出的相位分布图；

图8是利用洪水解包裹法展开相对相位所获得的绝对相位分布图；

图9是通过相位分布图得到的实际高度点云图。

具体实施方式

下面结合附图对本发明做更进一步的解释。

针对基于小波变换的三维测量法在高度跳变剧烈区域解相位不精确的问题，本发明采用小波变换法求取窗口尺度因子，并且利用其对频率因子的敏感性去除小波脊求取误差，改善跳变剧烈区域窗口大小求取精度的问题，从而提高这些区域解相位精度。在此基础上，利用小波变换可以使用快速傅里叶变换的特性，改善计算窗口大小时所占用的时间，使得整个过程速度大大加快。本发明不仅解决了跳变剧烈区域相位求取精度的问题，而且提高了处理效率。在得到精确的绝对相位分布后，根据经典光栅投影的相位到高度转换公式，最终得到测量物体的三维信息。

在windows操作***下使用VC++6.0作为编程工具，对CCD采集到的黑白变形条纹图像进行处理。该实例采用塑料护板作为被测物体，最终得到比较精确的含有护板三维信息的全场绝对相位分布。

如图1所示为本发明方法的流程图，具体实现步骤如下：

步骤1)，将黑白条纹图像投影到被测物体表面，使用CCD摄像头对被测物体表面进行拍摄，得到一幅高度为c、宽度为r的变形条纹图像g(x，y)，本实施例中c为1020，r为1020：

g(x，y)＝A(x，y)+B(x，y)cos [2πf₀x+φ(x，y)]    (1)式

其中，A(x，y)是背景光强分布，B(x，y)是物体表面反射率，f₀是正弦条纹频率，本实施例中f₀为0.05，φ(x，y)是待求的相对相位分布图，(x，y)表示变形条纹图像的二维坐标；图3为护板的变形条纹图像。

步骤2)，对变形条纹图g(x，y)逐行做小波变换，得到变形条纹图像每一点的近似尺度因子，具体过程如下：

步骤2.1)，将y视为常数y₁，采用一维小波变换对变形条纹图像g(x，y)的第y₁行进行处理，处理过程为：

${W_y}_1(a_1,b)={\∫}_{-\∞}^{+\∞}g(x,y_1){M^{*}}_{a,b}\left(x\right)\mathrm{dx}$   (2)式

其中，y₁表示某一行的行号；a是尺度因子，取值范围为10到90，每隔0.2取一个值；b是平移因子，取值范围为1到条纹的宽度r，每隔1取一个值，单位为像素；经一维小波变换获得的

是一个400行r列的二维复数矩阵，称第y₁行小波变换矩阵，a₁是矩阵

中元素的行标号；其中，

M^* _a，b(x)为M_a，b(x)的共轭函数；根据式(3)的小波函数M(x)得到M^* _a，b(x)表达式如式(4)：

$M\left(x\right)=\frac{1}{{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp \left(2\mathrm{\π}{\mathrm{if}}_{c}x\right)\exp \left(\frac{{-x}^2}{{2f}_b^2}\right)$   (3)式

${M^{*}}_{a,b}\left(x\right)=\frac{1}{\sqrt{a}{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp \left({-2\mathrm{\πif}}_c\frac{x-b}{a}\right)\exp [-\frac{{(x-b)}^2}{{2f}_b^2{a}^2}];$   (4)式

其中，f_b是小波函数的带宽，本实施例中取值为0.8，f_c是小波函数的中心频率，本实施例中取值为1，i为单位复数，x为自变量。

步骤2.2)，求取条纹图像g(x，y)的近似尺度因子分布图f_a1(x，y₁)：

求出二维复数矩阵的对应的模矩阵

其形式为：

  (5)式

单个点实验效果如图4所示，搜索模矩阵

第x列中值最大的元素，并将模矩阵的第x列中值最大元素的行标号赋值给a_max，则条纹图像的近似尺度因子分布图f_a1(x，y₁)在坐标(x，y₁)处数值a_rx为a_rx＝10+0.2×a_max；遍历条纹图像所有坐标点，求得条纹图像的近似尺度因子分布图f_a1(x，y)如图4所示。

步骤3)，对经步骤2.1)小波变换后的变形条纹图去除小波变换脊误差，得到变形条纹图像每一点的准确尺度因子，具体过程如下：

步骤3.1)，去除经步骤3.1)小波变换后的变形条纹图像每一个坐标点的脊误差，由于小波脊法只选取

之前的值，而在实际使用中很容易出现相位跳变剧烈的区域，

也不再是趋于零而可以消除的量，通常在

变化较大区域将趋于平缓，此时

对脊的求取精度影响很低，于是可以假定在跳变剧烈区域

为0。因此我们可以将相位表示为：

  (6)式

这样即使在处理过程中出现剧烈跳变区域，由于考虑了带来的误差，将大大减少计算误差，于是小波变换处理过程中简化过程为：

W(x，y₁)＝W₀+W₁+W₂+ε₀    (7)式

W₀，W₁，W₂分别为小波变换计算过程中的简化表达式，形式如下：

$W_0(a,b)=A(x,y)\frac{a\sqrt{\mathrm{\π}}}{{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp (-2f_b^2{f}_c^2{\mathrm{\π}}^2)$   (8)式

  (9)式

  (10)式

其中W(x，y₁)为变形条纹图中任意一像素点的小波变换简化形式，其中ε₀表示由于变形条纹图中每一点的相位二阶导所带来的误差，其表达式为：

  (11)式

其中t为中间计算变量，表示相位二阶导，表示相位一阶导，其中包含变量a，可见

较大时，

将偏离脊频。由于小波变换是对每一点不断地使用不同的尺度因子a来求取脊频(其中尺度因子a以0.2的幅度从10变化到90)，因此小波变换W(x，y₁)＝W₀+W₁+W₂+ε₀中每一个点逐频去除ε₀，再求取W_ε(x，y₁)＝W₀+W₁+W₂的脊，便可得到较为精确的

ε₀中只有

一个未知量，利用以下方法近似求取

考虑

后，f存在一个与

相关的偏差设为

  (12)式

对f求导得

  (13)式

对近似尺度因子逐行求导求出

求出最后求出对应尺度因子a的ε₀，得到ε₀为400行的一维复数数组在坐标(x，y₁)的误差数组ε₀(x，y₁)。W(x，y₁)中逐点减去误差ε₀(x，y₁)，求出精确小波变换脊的值W_ε(x，y₁)＝W₀+W₁+W₂，求出每一行精确小波变换脊矩阵W_400×r，矩阵W_400×r为400行r列的复数矩阵。

步骤3.2)，求取条纹图像的精确尺度因子分布图f_a2(x，y)：

求出W_400×r的对应的模矩阵C_400×r(a₁，b)，搜索模矩阵C₄₀₀×r(a₁，b)第x列中值最大的元素，并将模矩阵C_400×r(a₁，b)的第x列中值最大元素的行标号赋值给a_amax，则条纹图像的近似尺度因子分布图f_a2(x，y₁)在坐标(x，y₁)处数值a_acr为：a_acr＝10+0.2×a_amax；

遍历条纹图像所有坐标点，求得条纹图像的精确尺度因子分布图f_a2(x，y)如图6所示。

步骤4)，对变形条纹图逐行做窗口傅里叶变换，求出变形条纹图相对相位分布图，具体过程如下：

步骤4.1)，将y视为常数，采用一维窗口傅立叶变换对变形条纹图像g(x，y)的每行进行处理，一维窗口傅立叶变换过程为：

$\mathrm{WF}(b,\mathrm{\ξ})=\underset{-\∞}{\overset{+\∞}{\∫}}g\left(x\right)W_{\mathrm{\δ}}(x-b)\exp (-\mathrm{j\ξx})\mathrm{dx}$   (14)式

其频域表达式为：

$\mathrm{WF}(f_s,y)=\underset{n=0}{\overset{n=+\∞}{\mathrm{\Σ}}}{P}_n(f_s-{\mathrm{nf}}_0,y)$   (15)式

其中，g(x)为一行条纹图像；WF(b，ξ)表示一维窗口傅立叶变换，ξ表示频域计算因子，δ表示一维窗口傅立叶的窗口尺度因子，δ取值为一维窗口的位置(x-b，y)对应的精确尺度因子分布图f_a2(x，y)相应位置上的点的值，W_δ(x-b)表示窗口函数，表达式为：

$W_{\mathrm{\δ}}(x-b)=\frac{{\mathrm{\δ}}^{-1}}{\sqrt{2\mathrm{\π}}}\exp [-\frac{{(b-x)}^2}{{2\mathrm{\δ}}^2}]$   (16)式

其中，n表示阶次，依次取值为0，1，2，……，无穷；P_n(f_s-nf₀，y)表示任一点傅里叶变换后对应的n阶频谱；f_s表示频域的变量；

步骤4.2)，对傅里叶变换后的频谱滤波并提取相位信息，具体过程如下：

对P_n(f_s-nf₀，y)进行滤波并提取出一阶频谱P₁(f_s-f₀，y)，再对P₁(f_s-f₀，y)进行逆傅里叶变换，得到含有相位信息的B(x，y)exp[iφ(x，y)]，计算B(x，y)exp[iφ(x，y)]的角度值即可得到含有物体高度信息的变形条纹图相对相位值φ(x，y)，得到的相位值是介于0-2π之间；

遍历条纹图像所有坐标点，得到变形条纹图的相对相位分布图φ_A(x，y)，包裹相位实验图如图7所示。

步骤5)，建立条纹图像质量图，采用洪水算法展开相对相位分布图φ_A(x，y)，得到实际相位图

    具体过程如下：

步骤5.1)，利用相对相位图中的相位梯度建立质量图，相位梯度也可以指示该处相位求解的可靠性，相位梯度越大，说明该区域可能存在物体高度不连续或者求解错误等问题，可靠性较低，相位梯度越小，说明该区域应该属于平缓区域，可靠性较高；质量图可以按照以下公式计算：

$\mathrm{\Δ}(x,y)=\sqrt{W_{\mathrm{rap}}^2\{{\mathrm{\φ}}_A(x+1,y)-{\mathrm{\φ}}_A(x,y)\}+W_{\mathrm{rap}}^2\{{\mathrm{\φ}}_A(x,y+1)-{\mathrm{\φ}}_A(x,y)\}}$   (17)式

其中W_rap{}是包裹函数，当值大于2π或小于-2π时将其减去或加上2π。

步骤5.2)，在相对相位分布图中央找到质量值最高的坐标点，选择该点作为相位展开的起始点，将该点放入一个空的栈中；

步骤5.3)，判断栈是否为空；如果是，则相位展开过程结束，得到展开相位结果并进入步骤6)；如果否，弹出栈顶的点，展开该点的四邻点中没有处理的坐标点，并将这些未处理的点入栈。

步骤5.4)，按照质量图中的质量值将栈中所有点排序，质量最高的点放在栈顶，转到步骤5.3)继续处理，最终相位展开实验图如图8所示。

步骤6)，读取最终的展开相位结果

根据经典光栅投影的从展开相位结果

到物体高度h(x，y)的转换公式，最终求得测量物体的三维信息；其中转换公式h(x，y)为：

  (18)式

其中，l、d是测量***的几何参数，l是投影仪到测量平面的距离，本实施例中l为164厘米，d是CCD摄像头到投影仪的距离，本实施例中d为58厘米；

表示相位变化量，

为展开相位结果，

为初始相位结果，ω₀为投影光栅的角频率。图9为测量物体三维信息。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (1)

1.一种基于小波变换的改进窗口傅里叶三维测量法，其特征在于，包括如下步骤：

步骤1)，将黑白条纹图像投影到被测物体表面，使用CCD摄像头对被测物体表面进行拍摄，得到一幅高度为c、宽度为r的变形条纹图像g(x，y)：

g(x，y)＝A(x，y)+B(x，y)cos[2πf₀x+φ(x，y)]    (1)式

其中，A(x，y)是背景光强分布，B(x，y)是物体表面反射率，f₀是正弦条纹频率，φ(x，y)是待求的相对相位分布图，(x，y)表示变形条纹图像的二维坐标；

步骤2)，对所述变形条纹图g(x，y)逐行做小波变换，得到变形条纹图像每一点的近似尺度因子，具体过程如下：

步骤2.1)，将y视为常数y₁，采用一维小波变换对所述变形条纹图像g(x，y)的第y₁行进行处理，处理过程为：

${W_y}_1(a_1,b)={\∫}_{-\∞}^{+\∞}g(x,y_1){M^{*}}_{a,b}\left(x\right)\mathrm{dx}$   (2)式

其中，y₁表示某一行的行号；a是尺度因子，取值范围为10到90，每隔0.2取一个值；b是平移因子，取值范围为1到条纹的宽度r，每隔1取一个值，单位为像素；经一维小波变换获得的

是一个400行r列的二维复数矩阵，称第y₁行小波变换矩阵，a₁是矩阵

中元素的行标号；其中，

M^* _a，b(x)为M_a，b(x)的共轭函数；根据式(3)的小波函数M(x)得到M^* _a，b(x)表达式如式(4)：

$M\left(x\right)=\frac{1}{{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp \left(2{\mathrm{\πif}}_{c}x\right)\exp \left(\frac{{-x}^2}{{2f}_b^2}\right)$   (3)式

${M^{*}}_{a,b}\left(x\right)=\frac{1}{\sqrt{a}{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp (-{2\mathrm{\πif}}_c\frac{x-b}{a})\exp [-\frac{{(x-b)}^2}{{2f}_b^2{a}^2}];$   (4)式

其中，f_b是小波函数的带宽，f_c是小波函数的中心频率，i为单位复数，x为自变量；

步骤2.2)，求取条纹图像g(x，y)的近似尺度因子分布图f_a1(x，y₁)：

计算二维复数矩阵

的对应的模矩阵

搜索模矩阵

第x列中值最大的元素，并将模矩阵

的第x列中值最大元素的行标号赋值给a_max，则条纹图像的近似尺度因子分布图f_a1(x，y₁)在坐标(x，y₁)处数值a_rx为a_rx＝10+0.2×a_max；

遍历条纹图像所有坐标点，求得条纹图像的近似尺度因子分布图f_a1(x，y)；

步骤3)，对经所述步骤2.1)小波变换后的变形条纹图去除小波变换脊误差，得到变形条纹图像每一点的准确尺度因子，具体过程如下：

步骤3.1)，去除变形条纹图像每一个坐标点的脊误差，处理过程为：

W(x，y₁)＝W₀+W₁+W₂+ε₀    (5)式

其中W(x，y₁)为变形条纹图中任意一坐标点的小波变换简化形式，其中ε₀表示由变形条纹图中每一点的相位二阶导所带来的误差，其表达式为：

  (6)式

其中，t为中间计算变量，

表示相位二阶导，对近似尺度因子逐行求导求出

求出

最后求出对应尺度因子a的ε₀，得到ε₀为400行的一维复数数组在坐标(x，y₁)的坐标点上的误差数组ε₀(x，y₁)；W₀，W₁，W₂分别为小波变换计算过程中的简化表达式，形式如下：

$W_0(a,b)=A(x,y)\frac{a\sqrt{\mathrm{\π}}}{{\left(f_b^2\mathrm{\π}\right)}^{1/4}}\exp \left({-2f}_b^2{f}_c^2{\mathrm{\π}}^2\right)$   (7)式

  (8)式

  (9)式

其中，

为相位、

表示相位导数，W(x，y₁)中逐点减去误差ε₀(x，y₁)，求出精确小波变换脊的值W_ε(x，y₁)＝W₀+W₁+W₂，求出每一行精确小波变换脊矩阵W_400×r，矩阵W_400×r为400行r列的复数矩阵；

步骤3.2)，求取条纹图像的精确尺度因子分布图f_a2(x，y)：

求出W_400×r的对应的模矩阵C_400×r(a₁，b)，搜索模矩阵C_400×r(a₁，b)第x列中值最大的元素，并将模矩阵C_400×r(a₁，b)的第x列中值最大元素的行标号赋值给a_amax，则条纹图像的近似尺度因子分布图f_a2(x，y₁)在坐标(x，y₁)处数值a_acr为：a_acr＝10+0.2×a_amax；

遍历条纹图像所有坐标点，求得条纹图像的精确尺度因子分布图f_a2(x，y)；

步骤4)，对变形条纹图逐行做窗口傅里叶变换，求出变形条纹图相对相位分布图，具体过程如下：

步骤4.1)，将y视为常数，采用一维窗口傅立叶变换对变形条纹图像g(x，y)的每行进行处理，一维窗口傅立叶变换过程为：

$\mathrm{WF}(b,\mathrm{\ξ})=\underset{-\∞}{\overset{+\∞}{\∫}}g\left(x\right)W_{\mathrm{\δ}}(x-b)\exp (-\mathrm{j\ξx})\mathrm{dx}$   (10)式

其频域表达式为：

$\mathrm{WF}(f_s,y)=\underset{n=0}{\overset{n=+\∞}{\mathrm{\Σ}}}{P}_n(f_s-{\mathrm{nf}}_0,y)$   (11)式

其中，g(x)为一行条纹图像；WF(b，ξ)表示一维窗口傅立叶变换，ξ表示频域计算因子，δ表示一维窗口傅立叶的窗口尺度因子，δ取值为一维窗口的位置(x-b，y)对应的精确尺度因子分布图f_a2(x，y)相应位置上的点的值，W_δ(x-b)表示窗口函数，表达式为：

$W_{\mathrm{\δ}}(x-b)=\frac{{\mathrm{\δ}}^{-1}}{\sqrt{2\mathrm{\π}}}\exp [-\frac{{(b-x)}^2}{{2\mathrm{\δ}}^2}]$   (12)式

其中，n表示阶次，依次取值为0，1，2，……，无穷；P_n(f_s-nf₀，y)表示任一点傅里叶变换后对应的n阶频谱；f_s表示频域的变量；

步骤4.2)，对傅里叶变换后的频谱滤波并提取相位信息，具体过程如下：

对P_n(f_s-nf₀，y)进行滤波并提取出一阶频谱P₁(f_s-f₀，y)，再对P₁(f_s-f₀，y)进行逆傅里叶变换，得到含有相位信息的B(x，y)exp[iφ(x，y)]，计算B(x，y)exp[iφ(x，y)]的角度值即可得到含有物体高度信息的变形条纹图相对相位值φ(x，y)，得到的相位值是介于0-2π之间；

遍历条纹图像所有坐标点，得到变形条纹图的相对相位分布图φ_A(x，y)；

步骤5)，建立条纹图像质量图，采用洪水算法展开相对相位分布图φ_A(x，y)，得到实际相位图

具体过程如下：

步骤5.1)，利用相对相位分布图中的相位梯度建立质量图，质量图可以按照以下公式计算：

$\mathrm{\Δ}(x,y)=\sqrt{W_{\mathrm{rap}}^2\{{\mathrm{\φ}}_A(x+1,y)-{\mathrm{\φ}}_A(x,y)\}+W_{\mathrm{rap}}^2\{{\mathrm{\φ}}_A(x,y+1)-{\mathrm{\φ}}_A(x,y)\}}$   (13)式

其中W_rap{}是包裹函数，当值大于2π或小于-2π时将其减去或加上2π；

步骤5.2)，在相对相位分布图中央找到质量值最高的坐标点，选择该点作为相位展开的起始点，将该点放入一个空的栈中；

步骤5.3)，判断栈是否为空；如果是，则相位展开过程结束，得到展开相位结果并进入步骤6)；如果否，弹出栈顶的点，展开该点的四邻点中没有处理的坐标点，并将这些未处理的点入栈；

步骤5.4)，按照质量图中的质量值将栈中所有点排序，质量最高的点放在栈顶，转到步骤5.3)继续处理；

步骤6)，读取最终的展开相位结果

根据经典光栅投影的从展开相位结果

到物体高度h(x，y)的转换公式，最终求得测量物体的三维信息；其中转换公式h(x，y)为：

  (14)式

其中，l、d是测量***的几何参数，l是投影仪到测量平面的距离，d是CCD摄像头到投影仪的距离；

表示相位变化量，

为展开相位结果，

为初始相位结果，ω₀为投影光栅的角频率。

Images