CN102620685B - 一种基于史托克维尔变换的改进窗口傅里叶三维测量法 - 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)]

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

步骤2：对变形条纹图逐行做史托克维尔变换，得到变形条纹图像每一点的近似尺度因子，具体过程如下：

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

$S_{y_1}(b,f)={\∫}_{-\∞}^{+\∞}g(x,y_1)\mathrm{\ω}(b-x,f)\exp (-i2\mathrm{\πfx})\mathrm{dx}$

其中y₁表示某一行的行号，i为复数单位，f是尺度因子，f依次取值为0.001赫兹、0.002赫兹、0.003赫兹、……、1赫兹，b是平移因子，b依次取值为1、2、3、……、r，单位为像素，r为图像宽度，获得的是一个1000行r列的二维复数矩阵，ω(b-x，f)是史托克维尔变换窗口函数，窗口大小由尺度因子f控制，其表达式为：

$\mathrm{\ω}(b-x,f)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}\exp [-\frac{f^2{(b-x)}^2}{2}]$

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

求出的对应的模矩阵搜索模矩阵第x列中值最大的元素，并将模矩阵的第x列中值最大元素的行标号赋值给a_max，则a_rx＝0.001+0.001×a_max，a_rx为条纹图像的近似尺度因子分布图f_a1(x，y₁)在坐标(x，y₁)处数值，

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

步骤3：对变形条纹图逐行做史托克维尔变换并去除史托克维尔变换脊误差，得到变形条纹图像每一点的准确尺度因子，具体过程如下：

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

$S_{y_1}(b,f)={\∫}_{-\∞}^{+\∞}g(x,y_1)\mathrm{\ω}(b-x,f)\exp (-i2\mathrm{\πfx})\mathrm{dx}$

其中y₁表示某一行的行号，i为复数单位，f是尺度因子，f依次取值为0.001赫兹、0.002赫兹、0.003赫兹、……、1赫兹，b是平移因子，b依次取值为1、2、3、……、r，单位为像素，r为图像宽度，获得的是一个1000行r列的二维复数矩阵，ω(b-x，f)是史托克维尔变换窗口函数，窗口大小由尺度因子f控制，其表达式为：

$\mathrm{\ω}(b-x,f)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}\exp [-\frac{f^2{(b-x)}^2}{2}]$

步骤3.2：去除史托克维尔变换后的变形条纹图每一个像素点的脊误差，处理过程为：

S(x，y₁)＝S₀+S₁+S₂+ε₀

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

其中表示相位二阶导，对近似尺度因子逐行求导求出f′_al(x，y₁)，求出最后求出对应尺度因子f的ε₀，得到ε₀为1000行的一维复数数组在坐标(x，y₁)的像素点上的误差数组ε₀(x，y₁)，S₀，S₁，S₂分别为史托克维尔变换计算过程中的简化表达式，形式如下：

S₀(b，f)＝A(x，y)exp(-2π²)exp(-i2πfb)

S(x，y₁)中逐点减去误差ε₀(x，y₁)，求出精确史托克维尔变换脊的值S_ε(x，y₁)＝S₀+S₁+S₂，求出每一行精确史托克维尔变换脊矩阵S_1000×r，矩阵S_1000×r为1000行r列的复数矩阵，

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

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

遍历条纹图像所有坐标点，求得条纹图像的近似尺度因子分布图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}$

其频域表达式为：

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

其中，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}]$

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^2\{{\mathrm{\φ}}_A(x+1,y)-{\mathrm{\φ}}_A(x,y)\}+W^2\{{\mathrm{\φ}}_A(x,y+1)-{\mathrm{\φ}}_A(x,y)\}}$

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

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

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

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

步骤6：读取最终的展开相位结果根据经典光栅投影的从展开相位结果到物体高度h(x，y)的转换公式，最终求得测量物体的三维信息，所述的转换公式为：

其中，l、d是测量系统的几何参数，l是投影仪到测量平面的距离，d是CCD摄像头到投影仪的距离，表示相位变化量，为展开相位结果，为初始相位结果，ω₀为投影光栅的角频率。

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

附图说明

图1是本发明的整个过程的流程图。

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

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

图4是变形条纹图中某一点史托克维尔变换后的值取模值的曲线图，横坐标代表取样频率f，纵坐标代表史托克维尔变换后的模值。

图5是变形条纹图中某一点史托克维尔变换后的值去除脊误差的模值曲线图，横坐标代表取样频率f，纵坐标代表史托克维尔变换后去误差后的模值。

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

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

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

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

具体实施方式

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

图1是本发明的整个过程的流程图。

图2是本发明在步骤5中展开相位的过程具体流程图。

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

本发明基于改进史托克维尔变换三维测量方法，具体实现步骤如下：

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

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

其中，A(x，y)是背景光强分布，B(x，y)是物体表面反射率，实际处理中可以将A(x，y)和B(x，y)作为常数处理，f₀是正弦条纹频率，条纹图中固定条纹频率为20赫兹(HZ)，φ(x，y)是待求的相对相位分布图，(x，y)表示变形条纹图像的二维坐标，图3为变形条纹图像。

步骤2：对变形条纹图逐行做史托克维尔变换，得到变形条纹图像每一点的近似尺度因子，具体过程如下：

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

$S_{y_1}(b,f)={\∫}_{-\∞}^{+\∞}g(x,y_1)\mathrm{\ω}(b-x,f)\exp (-i2\mathrm{\πfx})\mathrm{dx}$

其中y₁表示某一行的行号，i为复数单位，f是尺度因子，f依次取值为0.001赫兹、0.002赫兹、0.003赫兹、……、1赫兹，【即：每隔0.001赫兹取一个值，】b是平移因子，b依次取值为1、2、3、……、r，【即：每隔1个像素取一个值，】单位为像素，r为图像宽度，获得的是一个1000行r列的二维复数矩阵，ω(b-x，f)是史托克维尔变换窗口函数，窗口大小由尺度因子f控制，其表达式为：

$\mathrm{\ω}(b-x,f)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}\exp [-\frac{f^2{(b-x)}^2}{2}]$

步骤2.2：求取条纹图像的近似尺度因子分布图f_a1(x，y₁)：求出的对应的模矩阵其形式为：

$C_{y_1}(b,f)=\sqrt{{\mathrm{imag}}^2\left[S_{y_1}(b,f)\right]+{\mathrm{real}}^2\left[S_{y_1}(b,f)\right]}$

单个点实验效果如图4，搜索模矩阵第x列中值最大的元素，并将模矩阵的第x列中值最大元素的行标号赋值给a_max，则a_rx＝0.001+0.001×a_max，a_rx为条纹图像的近似尺度因子分布图f_a1(x，y₁)在坐标(x，y₁)处数值，

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

步骤3：对变形条纹图逐行做史托克维尔变换并去除史托克维尔变换脊误差，得到变形条纹图像每一点的准确尺度因子，具体过程如下：

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

$S_{y_1}(b,f)={\∫}_{-\∞}^{+\∞}g(x,y_1)\mathrm{\ω}(b-x,f)\exp (-i2\mathrm{\πfx})\mathrm{dx}$

其中y₁表示某一行的行号，i为复数单位，f是尺度因子，f依次取值为0.001赫兹、0.002赫兹、0.003赫兹、……、1赫兹，【即：每隔0.001赫兹取一个值，】b是平移因子，b依次取值为1、2、3、……、r，【即：每隔1个像素取一个值，】单位为像素，r为图像宽度，获得的是一个1000行r列的二维复数矩阵，ω(b-x，f)是史托克维尔变换窗口函数，窗口大小由尺度因子f控制，其表达式为：

$\mathrm{\ω}(b-x,f)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}\exp [-\frac{f^2{(b-x)}^2}{2}]$

步骤3.2：去除史托克维尔变换后的变形条纹图每一个像素点的脊误差，由于S脊法只选取之前的值，而在实际使用中很容易出现相位跳变剧烈的区域，也不再是趋于零而可以消除的量，通常在变化较大区域将趋于平缓，此时对脊的求取精度影响很低，于是可以假定在跳变剧烈区域为0。因此我们可以将相位表示为

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

S(x，y₁)＝S₀+S₁+S₂+ε₀

S₀，S₁，S₂分别为史托克维尔变换计算过程中的简化表达式，形式如下：

S₀(b，f)＝A(x，y)exp(-2π²)exp(-i2πfb)

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

其中t为中间计算变量，表示相位二阶导，表示相位一阶导，其中包含变量f，可见较大时，将偏离脊频。由于史托克维尔变换是对每一点不断地使用不同频率的变换来求取脊频(其中采样频率以0.001HZ的幅度从0.001变化到1)，因此史托克维尔变换S(x，y₁)＝S₀+S₁+S₂+ε₀中每一个点逐频去除ε₀，再求取S_ε(x，y₁)＝S₀+S₁+S₂的脊便可得到较为精确的ε₀中只有一个未知量，利用以下方法近似求取

考虑后，f_b存在一个与相关的偏差设为

对f_b求导得

对近似尺度因子逐行求导求出f′_a1(x，y₁)，求出最后求出对应尺度因子f的ε₀，得到ε₀为1000行的一维复数数组在坐标(x，y₁)的误差数组ε₀(x，y₁)。S(x，y₁)中逐点减去误差ε₀(x，y₁)，求出精确史托克维尔变换脊的值S_ε(x，y₁)＝S₀+S₁+S₂，求出每一行精确史托克维尔变换脊矩阵S_1000×r，矩阵S_1000×r为1000行r列的复数矩阵，

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

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

遍历条纹图像所有坐标点，求得条纹图像的近似尺度因子分布图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}$

其频域表达式为：

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

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

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

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

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

对P_n(f-nf₀，y)进行滤波并提取出一阶频谱P₁(f-f₀，y)，再对P₁(f-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^2\{{\mathrm{\φ}}_A(x+1,y)-{\mathrm{\φ}}_A(x,y)\}+W^2\{{\mathrm{\φ}}_A(x,y+1)-{\mathrm{\φ}}_A(x,y)\}}$

其中W{}是包裹函数，当值大于2π或小于-2π时将其减去或加上2π，这样可以避免在真正的周期跳变的地方计算质量值错误。

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

步骤5.3：判断栈是否为空，如果是，说明相对相位分布图像中所有点已经得到处理，则相位展开过程结束；如果否，继续。弹出栈顶的点，展开该点的四邻点中没有处理的像素点，并将这些点入栈；

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

步骤6：读取最终的展开相位结果根据经典光栅投影的从展开相位结果到物体高度h(x，y)的转换公式，最终求得测量物体的三维信息，所述的转换公式为：

其中，l、d是测量系统的几何参数，l是投影仪到测量平面的距离，d是CCD摄像头到投影仪的距离，表示相位变化量，为展开相位结果，为初始相位结果，ω₀为投影光栅的角频率。图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)]

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

步骤2：对变形条纹图逐行做史托克维尔变换，得到变形条纹图像每一点的近似尺度因子，具体过程如下：

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

$S_{y_1}(b,f)={\∫}_{-\∞}^{+\∞}g(x,y_1)\mathrm{\ω}(b-x,f)\exp (-i2\mathrm{\πfx})\mathrm{dx}$

其中y₁表示某一行的行号，i为复数单位，f是尺度因子，f依次取值为0.001赫兹、0.002赫兹、0.003赫兹、……、1赫兹，b是平移因子，b依次取值为1、2、3、……、r，单位为像素，r为图像宽度，获得的是一个1000行r列的二维复数矩阵，ω(b-x,f)是史托克维尔变换窗口函数，窗口大小由尺度因子f控制，其表达式为：

$\mathrm{\ω}(b-x,f)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}\exp [-\frac{f^2{(b-x)}^2}{2}]$

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

求出的对应的模矩阵搜索模矩阵第x列中值最大的元素，并将模矩阵的第x列中值最大元素的行标号赋值给a_max，则a_rx＝0.001+0.001×a_max，a_rx为条纹图像的近似尺度因子分布图f_a1(x,y₁)在坐标(x,y₁)处数值，

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

步骤3：对变形条纹图逐行做史托克维尔变换并去除史托克维尔变换脊误差，得到变形条纹图像每一点的准确尺度因子，具体过程如下：

步骤3.1：将y视为常数y₁，采用一维史托克维尔变换对变形条纹图像g(x,y)的第y₁行进行处理，处理过程为：

$S_{y_1}(b,f)={\∫}_{-\∞}^{+\∞}g(x,y_1)\mathrm{\ω}(b-x,f)\exp (-i2\mathrm{\πfx})\mathrm{dx}$

其中y₁表示某一行的行号，i为复数单位，f是尺度因子，f依次取值为0.001赫兹、0.002赫兹、0.003赫兹、……、1赫兹，b是平移因子，b依次取值为1、2、3、……、r，单位为像素，r为图像宽度，获得的是一个1000行r列的二维复数矩阵，ω(b-x,f)是史托克维尔变换窗口函数，窗口大小由尺度因子f控制，其表达式为：

$\mathrm{\ω}(b-x,f)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}\exp [-\frac{f^2{(b-x)}^2}{2}]$

步骤3.2：去除史托克维尔变换后的变形条纹图每一个像素点的脊误差，处理过程为：

S(x,y₁)＝S₀+S₁+S₂+ε₀

其中S(x,y₁)为变形条纹图中任意一像素点的史托克维尔变换简化形式，其中ε₀表示由变形条纹图中每一点的相位二阶导所带来的误差，其表达式为：

其中表示相位二阶导，对近似尺度因子逐行求导求出f′_a1(x,y₁)，求出最后求出对应尺度因子f的ε₀，得到ε₀为1000行的一维复数数组在坐标(x,y₁)的像素点上的误差数组ε₀(x,y₁)，S₀，S₁，S₂分别为史托克维尔变换计算过程中的简化表达式，形式如下：

S₀(b,f)＝A(x,y)exp(-2π²)exp(-i2πfb)

其中和分别表示相位和相位的一阶导数，S(x,y₁)中逐点减去误差ε₀(x,y₁)，求出精确史托克维尔变换脊的值S_ε(x,y₁)＝S₀+S₁+S₂，求出每一行精确史托克维尔变换脊矩阵S_1000×r，矩阵S_1000×r为1000行r列的复数矩阵，

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

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

遍历条纹图像所有坐标点，求得条纹图像的近似尺度因子分布图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}$

其频域表达式为：

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

其中，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}]$

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^2\{{\mathrm{\φ}}_A(x+1,y)-{\mathrm{\φ}}_A(x,y)\}+W^2\{{\mathrm{\φ}}_A(x,y+1)-{\mathrm{\φ}}_A(x,y)\}}$

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

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

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

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

步骤6：读取最终的展开相位结果根据经典光栅投影的从展开相位结果到物体高度h(x,y)的转换公式，最终求得测量物体的三维信息，所述的转换公式为：

其中，l、d是测量系统的几何参数，l是投影仪到测量平面的距离，d是CCD摄像头到投影仪的距离，表示相位变化量，为展开相位结果，为初始相位结果，ω₀为投影光栅的角频率。