CN102628676B - 一种光学三维测量中的自适应窗口傅里叶相位提取法 - Google Patents
一种光学三维测量中的自适应窗口傅里叶相位提取法 Download PDFInfo
- Publication number
- CN102628676B CN102628676B CN201210016688.3A CN201210016688A CN102628676B CN 102628676 B CN102628676 B CN 102628676B CN 201210016688 A CN201210016688 A CN 201210016688A CN 102628676 B CN102628676 B CN 102628676B
- Authority
- CN
- China
- Prior art keywords
- centerdot
- imf
- square formation
- signal
- value
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B11/00—Measuring arrangements characterised by the use of optical techniques
- G01B11/24—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
- G01B11/25—Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures by projecting a pattern, e.g. one or more lines, moiré fringes on the object
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Databases & Information Systems (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Image Analysis (AREA)
- Length Measuring Devices By Optical Means (AREA)
Abstract
一种光学三维测量中的自适应窗口傅里叶相位提取法,是对灰度正弦光栅条纹图进行处理:(1)将一行条纹信号进行经验模态分解得到一系列本征模函数。(2)计算每个本征模函数的瞬时频率及边际谱,通过分析边际谱来选择能够准确描述条纹信号变化的瞬时频率。(3)确定条纹信号的背景分量。(4)用所确定的瞬时频率自适应定位局部平稳区域并据此计算高斯窗口的尺度因子。(5)将原信号去除背景分量后,对其进行自适应窗口傅里叶变换得到包裹相位分布。逐行进行上述步骤可得条纹图的全局包裹相位,展开后得到绝对相位并转换为高度信息,从而实现对物体的三维测量。本发明自适应性强,对面形复杂或含有陡峭边缘物体的测量精度有较大提高。
Description
技术领域
本发明属于三维信息重构的领域。基于灰度正弦光栅投影,采用希尔伯特-黄变换法分析变形条纹图像并求得条纹的瞬时频率及背景分量,用所求的瞬时频率自适应计算窗口傅里叶变换的窗口尺度因子,通过对去除背景的条纹信号直接做傅里叶变换或自适应窗口傅里叶变换从而得到精确的包裹相位,这种方法有利于提高对形貌复杂物体或携带陡峭边缘物体的测量精度。
背景技术
光学三维测量技术能够准确获得物体的三维面形数据,可用于三维模型重建、物体表面轮廓测量、工业环境中的尺寸和形位参数的检测等,因此它在虚拟现实、影视特技、医学整形和美容、工业产品的外观设计、艺术雕塑和文物保护等领域都有广阔的应用前景。
光栅投影法是一种重要的三维测量技术,通过向物体表面投射正弦光栅,将物体的高度信息以相位的形式隐含在光栅中,利用CCD获得物体表面的光栅条纹图像,并使用特定算法对条纹图像进行处理,提取其中的相位,从而建立物体的三维信息。
常用的求解条纹图像相位的方法有基于时间域的方法和基于变换域的方法。基于变换域的方法由于只需采集一幅变形条纹图即可完成相位的测量,有利于实现物体的动态测量,因此而被广泛研究和应用,其中包括傅里叶变换法、小波变换法、S变换法、窗口傅里叶变换法等。傅里叶变换是一种全局的信号分析工具,不能提取局部信号的特征,只适用于平稳信号。近年来,各种时频分析技术被广泛研究以期能够准确获取条纹图像的细节相位信息。连续小波变换由于具有多分辨率分析的特性被引入光学三维测量领域,通过检测小波变换最大脊处的相位,得到变形条纹图的包裹相位。但是这种方法有个前提条件,即被检测的相位必须是近似线性变化且变化缓慢,否则该方法的理论推导是不成立的。同时,小波母函数及相关参数的选取没有成熟的理论依据,也为小波变换轮廓术的广泛应用带来了限制。S变换与小波变换求解条纹图包裹相位的原理高度相似,故它也不可避免的受到与小波变换同样的条件限制。窗口傅里叶变换也具有较好的时频分析性能,但窗口尺寸即窗口尺度因子的自适应选取一直是研究的热点。现有的选取窗口尺度因子的方法大都是先通过连续小波变换或S变换检测最大脊处的尺度因子,将该尺度因子直接作为窗口傅里叶的尺度因子,或将该尺度因子取倒数作为条纹信号的瞬时频率,再通过相关算法计算窗口的尺度因子。这些方法需要事先设定基函数及经验值,因此自适应性不好,同时也同样受到被检测相位前提条件的限制,使求解的尺度因子或瞬时频率过平滑化,而不能很好的描述条纹信号的局部特征,从而导致这些方法只能测量表面变化相对平滑的物体,而对表面变化复杂或含有陡峭边缘物体的测量会受到较大限制。
发明内容
技术问题:针对光学三维测量中窗口傅里叶变换求取条纹图包裹相位的精确性和自适应性等问题,本发明在自适应进行窗口傅里叶变换及提高包裹相位的求解精度等方面给出了解决办法。本方法借助希尔伯特-黄变换,精确且自适应地求解能准确描述条纹图变化情况的瞬时频率,从而自适应计算窗口傅里叶变换的窗口尺度,同时,在不需要额外计算的情况下可以有效去除条纹图的背景分量从而大大减少窗口傅里叶变换处理过程中零频频谱的干扰。本方法能大大提高对复杂面形或含有陡峭边缘物体的测量精度。
技术方案:一种光学三维测量中的自适应窗口傅里叶相位提取法,步骤如下:
步骤1:将灰度正弦条纹图投影到被测物体表面,用CCD对被测物体表面进行拍摄,得到一幅宽度为c、高度为l的变形条纹图像g(y,x):
g(y,x)=A(y,x)+B(y,x)cos[2πf0x+φ(y,x)],
其中,A(y,x)是背景分量,B(y,x)是物体表面反射率,f0是正弦条纹基频,φ(y,x)是待求的包裹相位分布,(y,x)表示变形条纹图像中各像素点的二维坐标,取值范围分别为:1≤y≤l,1≤x≤c,B(y,x)cos[2πf0x+φ(y,x)]为基频分量,此处将每个像素点看作一个信号;
步骤2:令y=1,n=1,用经验模式分解法即EMD对g(y,x)进行分解,方法如下:
步骤2.1:将g(y,x)记为一行信号g(x),其中x仍然满足1≤x≤c,该行信号的相位为φ(x),寻找g(x)的极大值点和极小值点,分别采用公知的三次样条插值法对这些极大值点和极小值点进行插值,然后将这些值连接得到极大值包络线gmax(x)和极小值包络线gmin(x);
步骤2.2:用原始信号g(x)减去极大值包络线gmax(x)与极小值包络线gmin(x)的平均值,得到h(x):
步骤2.3:分别计算h(x)的平均幅度mean_amp(x)和包络幅度env_amp(x):
步骤2.4:如果h(x)同时满足以下三个条件,则得到一个本征模函数IMF,记为sn(x),且sn(x)=h(x),同时将sn(x)从原信号g(x)中分离,得到新的信号g′(x)=g(x)-sn(x),否则,直接将h(x)从原信号g(x)中分离,得到新的信号g′(x)=g(x)-h(x),所述的三个条件是:
条件1:mean_amp(x)/env_amp(x)<0.5,
条件2:满足不等式mean_amp(x)/env_amp(x)<0.5的像素个数占同一整行像素总数的比例小于0.05,
条件3:极大值和极小值的个数之和与h(x)过零点的个数相等或至多相差1个,
步骤2.5:若步骤2.4得到的g′(x)也同时满足以上三个条件,令g(x)=g′(x),n=n+1,返回步骤2.1;否则分解停止,并将g′(x)记为res(x),得到最后的分解结果为:
其中n为一个IMF即s(x)的序数,1≤n≤N,N为IMF的总个数;
步骤3:通过希尔伯特-黄变换,确定能够准确描述一行条纹信号变化规律的瞬时频率,具体过程如下:
步骤3.1:对第n个IMF即sn(x)做希尔伯特变换,得到:
其中“*”为卷积运算符,τ为积分变量,sn′(x)为希尔伯特变换的结果,分别构造sn(x)的解析信号zn(x):
步骤3.2:对第n个IMF求其瞬时频率fn(x):
对第n个IMF求其边际谱Mn(f):
步骤3.3:确定含基频分量最多的IMF的序数K:
选取fK(x)作为能够准确描述该行条纹信号变化规律的瞬时频率;
步骤4:确定一行条纹信号的背景分量,具体过程如下:根据步骤3.2已得的每个IMF的瞬时频率,求出每个IMF的瞬时频率均值,并从中找到最小的瞬时频率均值,然后确定最小的瞬时频率均值所对应的IMF序数,记为Kb,则以步骤2.5所得的第Kb+1个IMF到第N个IMF以及残余分量之和即为该行条纹信号的背景分量组合;
步骤5:根据一个局部平稳区域内的瞬时频率最大值小于2倍最小值这一属性,用步骤3.3确定的瞬时频率fK(x)自适应定位一行条纹信号的局部平稳区域,具体过程如下:
步骤5.3:对方阵F建立坐标系,坐标原点为F左上角第一个元素,其坐标值记为(1,1),横坐标的坐标方向为方阵的行方向,横坐标取值范围为1~c,纵坐标的坐标方向为方阵的列方向,纵坐标的取值范围仍为1~c,
在方阵F中寻找出以方阵F的零对角线的一部分或全部为对角线的最大全零子方阵,从方阵F的坐标原点开始,沿零对角线方向,依次将这些最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标记录下来,依次用(p1,p1),(p2,p2),……,(pr,pr),(c,c)表示,其中(p1,p1)为第一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标,依次排列,(pr,pr)为倒数第二个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标,(c,c)为最后一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标,
最后,该行条纹信号的平稳区域划分情况为:(1,……,p1),(p1+1,……,p2),……,(pr+1,……,c);
步骤6:将原条纹信号g(x)减去步骤4中确定的背景分量,得到步骤1中去除A(y,x)后的基频分量g0(x),即为g0(x)=B(x,y)cos[2πf0x+φ(x,y)];
步骤7:判断步骤5确定的局部平稳区域个数是否为1,若“是”,则直接对g0(x)做快速傅里叶变换得到傅里叶频谱,记作Ff,若“否”,则对g0(x)进行自适应高斯窗口傅里叶变换,具体过程为:
步骤7.1:对g0(x)做自适应高斯窗口傅里叶变换:
其中b为平移因子,b依次取值为1、2、3、…、c,a为每一个局部平稳区域内每个像素点对应的高斯窗口的尺度因子,a=L/4,其中L为相应的局部平稳区域的长度值,单位为像素, 为高斯窗函数,
对一行条纹信号所有的信号依次进行自适应高斯窗口傅里叶变换后,得到由构成的一个二维复数方阵,大小为c×c,二维复数方阵每一行的元素为每一个窗口内信号的频谱,二维复数方阵一共有c行,代表b从1取值到c即一共有c个窗口的信号的频谱;
步骤7.2:将步骤7.1所得的二维复数方阵按列叠加,得到总频谱,记作Ff;
步骤8:根据Ff确定基频频谱的范围,并提取出来,记为Ff0;对Ff0求傅里叶逆变换,并根据对Ff0求傅里叶逆变换所得的结果,求得相位角即包裹在[-π,π]之间的相位φ(x),令y=y+1,返回步骤2,若y>l,则进入步骤9;
步骤9:对包裹相位进行展开得到绝对相位,根据经典光栅投影的相位到高度的转换公式,最终求得测量物体的三维信息。
有益效果:
相比相移法,本方法不会受到设备关于正弦性、相移的精度和速度等条件限制,且仅需要一幅变形条纹图即可完成包裹相位的提取,有利于实现动态测量。与其他基于变换域的现有技术相比,本发明具有如下优点:
首先,本发明采用希尔伯特-黄变换法求解条纹信号的瞬时频率,相比于现有常用的小波变换最大脊法或者S变换最大脊法,避免了被测相位必须线性逼近且变化缓慢的条件限制,得到的瞬时频率更加真实的描述了变形条纹信号的变化,通过本方法所得的瞬时频率确定的窗口尺度因子能真实的根据信号的局部特征进行自适应变化,从而达到窗口傅里叶精确提取局部细节相位的目的;
其次,本发明提出的自适应定位局部平稳区域的方法,不需要事先确定任何经验值,也不需要迭代计算,效率高、速度快,大大提高了整个自适应窗口傅里叶相位提取的处理效率;
最后,本发明在求解条纹图的瞬时频率时,可以很方便的“顺便”去除条纹图的背景分量,大大减少了提取基频分量时背景分量所带来的零频频谱干扰,当确定的窗口尺度因子较小时,这一干扰将会非常严重,因此本发明这一处理会有很大的优势。
综上,本发明仅用一幅变形条纹图,便能够准确地获取被测物体的包裹相位,且克服了传统方法对具有复杂面形或陡峭边缘的物体测量不高的缺点。
附图说明
图1是本发明整个过程的流程图。
图2是CCD采集到的以塑料泡沫为被测物体的变形条纹图像。
图3是步骤2中对一行条纹信号进行经验模式分解的具体过程流程图。
图4是图2中直线所代表的一行条纹信号及其EMD分解结果。
图5是图4中相应IMF的瞬时频率。
图6是图4中相应IMF的边际谱。
图7是步骤5的示意图。
图8是对图2中直线所在行条纹信号所定位的局部平稳区域。
图9是去掉背景分量的变形条纹图。
图10是所求的尺度因子分布图。
图11是所求的包裹相位图。
图12是最终恢复的相位图。
具体实施方式
下面结合附图对本发明的具体实施方式做进一步说明。在Windows操作系统下选用VC++6.0作为编程工具,对CCD采集到的变形条纹图像进行处理。该实例采用塑料泡沫作为被测物体,最终得到比较精确的含有三维信息的全场包裹相位分布和。
图1是本发明的整个过程的流程图。
光学三维测量中条纹图进行包裹相位提取时,针对现有的时频分析方法对表面复杂或含有陡峭边缘的物体的测量精度较低等问题,本发明提出基于希尔伯特-黄变换的多尺度窗口傅里叶相位提取法。首先采用经验模态分解法逐行分解采集到的变形条纹图像并进行希尔伯特-黄变换处理,逐行分析确定能够描述变形条纹变化情况的瞬时频率,同时确定条纹图的背景分量。将原条纹图去除背景分量。根据确定的瞬时频率逐行定位每一行信号的局部平稳区域,若任意一行信号的局部平稳区域个数为1,则对该行信号直接进行傅里叶变换然后提取基频分量的频谱。若任意一行信号的局部平稳区域个数不为1,根据所确定的局部平稳区域的长度,计算每一个局部平稳区域内每个信号相应的高斯窗函数的尺度因子。由不同的尺度因子引导不同尺寸的高斯窗,对相应的信号逐点进行由高斯窗函数引导的窗口傅里叶变换,每个像素点的变换结果都是一个一维复数数组即傅里叶频谱,将基频频谱提取出来。对一行信号逐个信号进行相同的处理,所有的处理结果依次按行排列,得到一个二维复数矩阵。根据高斯窗函数的时频统计特性,将所有的提取出的基频频谱按列相加,得到最终的一行信号的基频频谱。对相加后的一行信后的基频频谱进行傅里叶反变换,得到一行信号的基频分量的相位。对条纹图逐行进行处理后,最终可得整幅条纹图的全场包裹相位分布。由于所采用的方法不会使所求的瞬时频率被局部平滑化,故根据其所确定的窗口尺寸能够较好的反映条纹信号的局部信息,从而能够精确地测量条纹图细节处的相位。同时,由于背景分量被去除,故在基频提取的过程中极大地减少了零频频谱的干扰,有利于基频频谱的精确提取。
本发明的具体实施步骤如下:
步骤1:将灰度正弦条纹图投影到被测物体表面,用CCD对被测物体表面进行拍摄,得到一幅宽度为c、高度为l的变形条纹图像g(y,x):
g(y,x)=A(y,x)+B(y,x)cos[2πf0x+φ(y,x)],
其中,A(y,x)是背景分量,B(y,x)是物体表面反射率,f0是正弦条纹基频,φ(y,x)是待求的包裹相位分布,(y,x)表示变形条纹图像中各像素点的二维坐标,取值范围分别为:1≤y≤l,1≤x≤c,B(y,x)cos[2πf0x+φ(y,x)]为基频分量,此处将每个像素点看作一个信号。
图2为所得的变形条纹图像,大小为1020×1020,图中直线为y=170时的一行示例条纹信号g(x)。
步骤2:令y=1,n=1,用经验模式分解法即EMD对g(y,x)进行分解,方法如下:
步骤2.1:将g(y,x)记为一行信号g(x),其中x仍然满足1≤x≤c,该行信号的相位为φ(x),寻找g(x)的极大值点和极小值点,分别采用公知的三次样条插值法对这些极大值点和极小值点进行插值,然后将这些值连接得到极大值包络线gmax(x)和极小值包络线gmin(x)。此处也可选择其他插值法,具体的选取办法可参考文献“N.E.Huang,Z.Shen,and S.R.Long et al,“The empirical mode decompositionand the Hilbert spectrum for nonliner and non-stationary time series analysis,”J.Proc.Lond.A.Great Britain,454,903-995(1998)”,及“G.Rilling,P.Flandrin and P.Goncalves,“Onempirical mode decomposition and its algorithm,”in IEEEEURASIP Workshop onNonlinear Signal andImage Processing,NSIP-03,GRADO(I)(IEEE,2003)”;
步骤2.2:用原始信号g(x)减去极大值包络线gmax(x)与极小值包络线gmin(x)的平均值,得到h(x):
步骤2.3:分别计算h(x)的平均幅度mean_amp(x)和包络幅度env_amp(x):
步骤2.4:如果h(x)同时满足以下三个条件,则得到一个本征模函数IMF,记为sn(x),且sn(x)=h(x),同时将sn(x)从原信号g(x)中分离,得到新的信号g′(x)=g(x)-sn(x),否则,直接将h(x)从原信号g(x)中分离,得到新的信号g′(x)=g(x)-h(x),所述的三个条件是:
条件1:mean_amp(x)/env_amp(x)<0.5;
条件2:满足不等式mean_amp(x)/env_amp(x)<0.5的像素个数占同一整行像素总数的比例小于0.05;
条件3:极大值和极小值的个数之和与h(x)过零点的个数相等或至多相差1个;
步骤2.5:若步骤2.4得到的g′(x)也同时满足以上三个条件,令g(x)=g′(x),n=n+1,返回步骤2.1;否则分解停止,并将g′(x)记为res(x),得到最后的分解结果为:
其中n为一个IMF即s(x)的序数,1≤n≤N,N为IMF的总个数。
图3为一行信号分解过程的流程图。这里以变形条纹图的一行变形信号为例,即对如图2中直线所代表信号g(x)进行处理。在处理过程中,三个条件的阈值0.5、0.05和0.05为建议阈值,实际应用过程中可根据需要适当微调整,调整的前提是使分解彻底,即剩余分量不能再分离出IMF,同时又不能过分解。图4为原信号及其被分解后的结果,从中可看见,IMFs随着尺度的变化由高到底依次排列,第一个IMF为高频噪声分量,res(x)的尺度最大,趋近于0,它能够描述一行信号的整体趋势。
步骤3:通过希尔伯特-黄变换,确定能够准确描述一行条纹信号变化规律的瞬时频率,具体过程如下:
步骤3.1:对第n个IMF即sn(x)做希尔伯特变换,得到:
其中“*”为卷积运算符,τ为积分变量,sn′(x)为希尔伯特变换的结果,分别构造sn(x)的解析信号zn(x):
步骤3.2:对第n个IMF求其瞬时频率fn(x):
对第n个IMF求其边际谱Mn(f):
步骤3.3:确定含基频分量最多的IMF的序数K:
选取fK(x)作为能够准确描述该行条纹信号变化规律的瞬时频率。
如图5为图4中每个IMF所得的瞬时频率,图6为相应的边际谱。图6中,由于条纹的基频是0.05,可见IMF2的边际谱最大值所对应的频率值与基频0.05的差值最小,故IMF2为含有基频分量最多的分量,IMF-f2则被选定为能够准确描述该行条纹信号变化规律的瞬时频率。
步骤4:确定一行条纹信号的背景分量,具体过程如下:根据步骤3.2已得的每个IMF的瞬时频率,求出每个IMF的瞬时频率均值,并从中找到最小的瞬时频率均值,然后确定最小的瞬时频率均值所对应的IMF序数,记为Kb,则以步骤2.5所得的第Kb+1个IMF到第N个IMF以及残余分量之和即为该行条纹信号的背景分量组合。
步骤5:根据一个局部平稳区域内的瞬时频率最大值小于2倍最小值这一属性,用步骤3.3确定的瞬时频率fK(x)自适应定位一行条纹信号的局部平稳区域,具体过程如下:
令方阵内所有负值为0,并找到方阵中全部元素为0的零对角线,得到:
步骤5.3:对方阵F建立坐标系,坐标原点为F左上角第一个元素,其坐标值记为(1,1),横坐标的坐标方向为方阵的行方向,横坐标取值范围为1~c,纵坐标的坐标方向为方阵的列方向,纵坐标的取值范围仍为1~c,
在方阵F中寻找出以方阵F的零对角线的一部分或全部为对角线的最大全零子方阵,从方阵F的坐标原点开始,沿零对角线方向,依次将这些最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标记录下来,依次用(p1,p1),(p2,p2),……,(pr,pr),(c,c)表示,其中(p1,p1)为第一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标,依次排列,(pr,pr)为倒数第二个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标,(c,c)为最后一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标,
最后,该行条纹信号的平稳区域划分情况为:(1,……,p1),(p1+1,……,p2),……,(pr+1,……,c)。
图7为自适应定位一行条纹信号局部平稳区域的直观描述,图8为图3中对直线代表信号局部平稳区域的划分。
步骤6:将原条纹信号g(x)减去步骤4中确定的背景分量,得到步骤1中去除A(y,x)后的基频分量g0(x),即为g0(x)=B(x,y)cos[2πf0x+φ(x,y)]。
图9为去除背景分量后的条纹图。
步骤7:判断步骤5确定的局部平稳区域个数是否为1,若“是”,则直接对g0(x)做快速傅里叶变换得到傅里叶频谱,记作Ff,若“否”,则对g0(x)进行自适应高斯窗口傅里叶变换,具体过程为:
步骤7.1:对g0(x)做自适应高斯窗口傅里叶变换:
其中b为平移因子,b依次取值为1、2、3、…、c,a为每一个局部平稳区域内每个像素点对应的高斯窗口的尺度因子,a=L/4,其中L为相应的局部平稳区域的长度值,单位为像素, 为高斯窗函数,
对一行条纹信号所有的信号依次进行自适应高斯窗口傅里叶变换后,得到由构成的一个二维复数方阵,大小为c×c,二维复数方阵每一行的元素为每一个窗口内信号的频谱,二维复数方阵一共有c行,代表b从1取值到c即一共有c个窗口的信号的频谱;
步骤7.2:将步骤7.1所得的二维复数方阵按列叠加,得到总频谱,记作Ff。
步骤8:根据Ff确定基频频谱的范围,并提取出来,记为Ff0;对Ff0求傅里叶逆变换,并根据对Ff0求傅里叶逆变换所得的结果,求得相位角即包裹在[-π,π]之间的相位φ(x),令y=y+1,返回步骤2,若y>l,则进入步骤9。
步骤9:对包裹相位进行展开得到绝对相位,此处相位展开的方法采用质量图相位展开法。最后,根据经典光栅投影的相位到高度的转换公式求得测量物体的三维信息。高度转换公式如下:
其中,l、d是测量系统的几何参数,l是投影仪到测量平面的距离,d是CCD摄像头到投影仪的距离,表示相位变化量,为展开相位结果,为初始相位结果,由测量参考面决定,ω0为投影光栅的角频率,由系统标定可得。
如图10为尺度因子的全场分布图,可见用本方法能够较好地描述物体的变化规律。图11为对一整幅变形条纹图进行本发明方法处理后所得的全场包裹相位分布图。图12为所测物体的相位变化分布图,即先将全场包裹相位用质量图法进行相位展开,然后用同样方法得到没有被物体调制的参考平面条纹图的绝对相位图,将包含物体三维信息的变形条纹图的展开相位与参考平面条纹图的展开相位相减即可得到相位变化分布图。用本方法所测的包裹相位,在洪水法进行相位展开时,由于物体边缘相位测量的精度较高,基本没有误差传递,而其他方法会有大片的误差传递现象。故本方法对复杂面形或含有陡峭边缘物体的测量时精度会有具有较大的提高。
Claims (1)
1.一种光学三维测量中的自适应窗口傅里叶相位提取法,具体步骤如下:
步骤1:将灰度正弦条纹图投影到被测物体表面,用CCD对被测物体表面进行拍摄,得到一幅宽度为c、高度为l的变形条纹图像g(y,x):
g(y,x)=A(y,x)+B(y,x)cos[2πf0x+φ(y,x)],
其中,A(y,x)是背景分量,B(y,x)是物体表面反射率,f0是正弦条纹基频,φ(y,x)是待求的包裹相位分布,(y,x)表示变形条纹图像中各像素点的二维坐标,取值范围分别为:1≤y≤L1≤x≤c,B(y,x)cos[2πf0x+φ(y,x)]为基频分量,此处将每个像素点看作一个信号;
其特征在于,
步骤2:令y=1,n=1,用经验模式分解法即EMD对g(y,x)进行分解,方法如下:
步骤2.1:将g(y,x)记为一行信号g(x),其中x仍然满足1≤x≤c,该行信号的相位为φ(x),寻找g(x)的极大值点和极小值点,分别采用公知的三次样条插值法对这些极大值点和极小值点进行插值,然后将这些值连接得到极大值包络线gmax(x)和极小值包络线gmin(x);
步骤2.2:用原始信号g(x)减去极大值包络线gmax(x)与极小值包络线gmin(x)的平均值,得到h(x):
步骤2.3:分别计算h(x)的平均幅度mean_amp(x)和包络幅度env_amp(x):
步骤2.4:如果h(x)同时满足以下三个条件,则得到一个本征模函数IMF,记为sn(x),且sn(x)=h(x),同时将sn(x)从原信号g(x)中分离,得到新的信号g′(x)=g(x)-sn(x),否则,直接将h(x)从原信号g(x)中分离,得到新的信号g′(x)=g(x)-h(x),所述的三个条件是:
条件1:mean_amp(x)/env_amp(x)<0.5,
条件2:满足不等式mean_amp(x)/env_amp(x)<0.5的像素个数占同一整行像素总数的比例小于0.05,
条件3:极大值和极小值的个数之和与h(x)过零点的个数相等或至多相差1个,
步骤2.5:若步骤2.4得到的g′(x)也同时满足以上三个条件,令g(x)=g′(x),n=n+1,返回步骤2.1;否则分解停止,并将g′(x)记为res(x),得到最后的分解结果为:
其中n为一个IMF即s(x)的序数,1≤n≤N,N为IMF的总个数;
步骤3:通过希尔伯特-黄变换,确定能够准确描述一行条纹信号变化规律的瞬时频率,具体过程如下:
步骤3.1:对第n个IMF即sn(x)做希尔伯特变换,得到:
其中“*”为卷积运算符,τ为积分变量,sn′(x)为希尔伯特变换的结果,
分别构造sn(x)的解析信号zn(x):
步骤3.2:对第n个IMF求其瞬时频率fn(x):
对第n个IMF求其边际谱Mn(f):
步骤3.3:确定含基频分量最多的IMF的序数K:
选取fK(x)作为能够准确描述该行条纹信号变化规律的瞬时频率:
步骤4:确定一行条纹信号的背景分量,具体过程如下:根据步骤3.2已得的每个IMF的瞬时频率,求出每个IMF的瞬时频率均值,并从中找到最小的瞬时频率均值,然后确定最小的瞬时频率均值所对应的IMF序数,记为Kb,则以步骤2.5所得的第Kb+1个IMF到第N个IMF以及残余分量之和即为该行条纹信号的背景分量组合;
步骤5:根据一个局部平稳区域内的瞬时频率最大值小于2倍最小值这一属性,用步骤3.3确定的瞬时频率fK(x)自适应定位一行条纹信号的局部平稳区域,具体过程如下:
方阵F:
令方阵内所有负值为0,并找到方阵中全部元素为0的零对角线,得到:
步骤5.3:对方阵F建立坐标系,坐标原点为F左上角第一个元素,其坐标值记为(1,1),横坐标的坐标方向为方阵的行方向,横坐标取值范围为1~c,纵坐标的坐标方向为方阵的列方向,纵坐标的取值范围仍为1~c,
在方阵F中寻找出以方阵F的零对角线的一部分或全部为对角线的最大全零子方阵,从方阵F的坐标原点开始,沿零对角线方向,依次将这些最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标记录下来,依次用(p1,p1),(p2,p2),……,(pr,pr),(c,c)表示,其中(p1,p1)为第一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标,依次排列,(pr,pr)为倒数第二个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标,(c,c)为最后一个最大全0子方阵的零对角线上最后一个元素相对于方阵F的坐标,
最后,该行条纹信号的平稳区域划分情况为:(1,……,p1),(p1+1,……,p2),……,(pr+1,……,c);
步骤6:将原条纹信号g(x)减去步骤4中确定的背景分量,得到步骤1中去除A(y,x)后的基频分量g0(x),即为g0(x)=B(x,y)cos[2πf0x+φ(x,y)];
步骤7:判断步骤5确定的局部平稳区域个数是否为1,若“是”,则直接对g0(x)做快速傅里叶变换得到傅里叶频谱,记作Ff,若“否”,则对g0(x)进行自适应高斯窗口傅里叶变换,具体过程为:
步骤7.1:对g0(x)做自适应高斯窗口傅里叶变换:
其中b为平移因子,b依次取值为1、2、3、…、c,a为每一个局部平稳区域内每个像素点对应的高斯窗口的尺度因子,a=L/4,其中L为相应的局部平稳区域的长度值,单位为像素, 为高斯窗函数,
对一行条纹信号所有的信号依次进行自适应高斯窗口傅里叶变换后,得到由构成的一个二维复数方阵,大小为c×c,二维复数方阵每一行的元素为每一个窗口内信号的频谱,二维复数方阵一共有c行,代表b从1取值到c即一共有c个窗口的信号的频谱;
步骤7.2:将步骤7.1所得的二维复数方阵按列叠加,得到总频谱,记作Ff;
步骤8:根据Ff确定基频频谱的范围,并提取出来,记为Ff0;对Ff0求傅里叶逆变换,并根据对Ff0求傅里叶逆变换所得的结果,求得相位角即包裹在[-π,π]之间的相位φ(x),令y=y+1,返回步骤2,若y>l,则进入步骤9;
步骤9:对包裹相位进行展开得到绝对相位,根据经典光栅投影的相位到高度的转换公式,最终求得测量物体的三维信息。
Priority Applications (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210016688.3A CN102628676B (zh) | 2012-01-19 | 2012-01-19 | 一种光学三维测量中的自适应窗口傅里叶相位提取法 |
KR1020137005460A KR101475382B1 (ko) | 2012-01-19 | 2012-02-28 | 광학적 3차원 측량의 자기 적응 윈도우 푸리에 위상추출방법 |
PCT/CN2012/071731 WO2013107076A1 (zh) | 2012-01-19 | 2012-02-28 | 一种光学三维测量中的自适应窗口傅里叶相位提取法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201210016688.3A CN102628676B (zh) | 2012-01-19 | 2012-01-19 | 一种光学三维测量中的自适应窗口傅里叶相位提取法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN102628676A CN102628676A (zh) | 2012-08-08 |
CN102628676B true CN102628676B (zh) | 2014-05-07 |
Family
ID=46586983
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201210016688.3A Expired - Fee Related CN102628676B (zh) | 2012-01-19 | 2012-01-19 | 一种光学三维测量中的自适应窗口傅里叶相位提取法 |
Country Status (3)
Country | Link |
---|---|
KR (1) | KR101475382B1 (zh) |
CN (1) | CN102628676B (zh) |
WO (1) | WO2013107076A1 (zh) |
Families Citing this family (37)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104034285B (zh) * | 2014-06-25 | 2017-01-04 | 西北工业大学 | 整数线性规划搜索法的双频正弦光栅绝对相位解包裹方法 |
CN104200474B (zh) * | 2014-09-04 | 2017-03-01 | 华中科技大学 | 一种获取物体变形量的数字图像分析方法 |
CN104299211B (zh) * | 2014-09-25 | 2020-12-29 | 苏州笛卡测试技术有限公司 | 一种自由移动式三维扫描方法 |
CN104949621B (zh) * | 2015-06-04 | 2017-08-29 | 广东工业大学 | 一种光栅尺条纹的边界定位方法 |
CN105091750A (zh) * | 2015-07-30 | 2015-11-25 | 河北工业大学 | 一种基于双四步相移的投影仪标定方法 |
CN106568394A (zh) * | 2015-10-09 | 2017-04-19 | 西安知象光电科技有限公司 | 一种手持式三维实时扫描方法 |
KR20170047780A (ko) | 2015-10-23 | 2017-05-08 | 한국전자통신연구원 | 적응적 윈도우 마스크를 이용하는 로우 코스트 계산장치 및 그 방법 |
CN106802137B (zh) * | 2017-01-16 | 2019-04-02 | 四川大学 | 一种相位展开方法及系统 |
CN106815840B (zh) * | 2017-01-22 | 2020-06-05 | 飞依诺科技(苏州)有限公司 | 一种肝部扫查图像的处理方法及装置 |
CN106901697A (zh) * | 2017-03-03 | 2017-06-30 | 哈尔滨理工大学 | 一种用于测试三维傅里叶变换胸腹表面测量手段的方法 |
CN107180640B (zh) * | 2017-04-13 | 2020-06-12 | 广东工业大学 | 一种相位相关的高密度叠窗频谱计算方法 |
CN109087348B (zh) * | 2017-06-14 | 2022-04-29 | 北京航空航天大学 | 一种基于自适应区域投射的单像素成像方法 |
CN107392867A (zh) * | 2017-07-14 | 2017-11-24 | 中国科学院遥感与数字地球研究所 | 基于emd的sar图像自适应辐射均衡方法 |
CN107786816A (zh) * | 2017-09-14 | 2018-03-09 | 天津大学 | 基于曝光补偿的自适应投影方法 |
CN108253907B (zh) * | 2018-02-01 | 2020-07-21 | 深圳市易尚展示股份有限公司 | 基于希尔伯特变换相位误差校正的三维测量方法和装置 |
CN109299431B (zh) * | 2018-08-27 | 2022-11-04 | 杭州电子科技大学 | 基于希尔伯特边际谱特征的光伏组件性能一致性评价方法 |
CN109633268B (zh) * | 2018-12-20 | 2020-03-10 | 北京航空航天大学 | 一种基于b样条和直方图的方波基频辨识方法 |
CN110037724B (zh) * | 2019-03-14 | 2023-09-12 | 杭州惜尔信息技术有限公司 | 一种基于st变换的ct成像方法 |
CN110175541B (zh) * | 2019-05-13 | 2022-06-14 | 武汉大学 | 一种海平面变化非线性趋势提取的方法 |
KR102218015B1 (ko) * | 2019-08-06 | 2021-02-19 | 한국표준과학연구원 | Fpm에서 켈리브레이션 그레이팅을 이용한 회절격자 위상 측정시스템 및 측정방법 |
CN110503025B (zh) * | 2019-08-19 | 2023-04-18 | 重庆大学 | 一种基于半监督协同训练的模拟电路早期故障诊断方法 |
CN111028334A (zh) * | 2019-11-26 | 2020-04-17 | 北京理工大学珠海学院 | 基于枝切法的三维立体重构方法、系统和存储介质 |
CN110991564B (zh) * | 2019-12-24 | 2023-05-26 | 安徽工业大学 | 基于多尺度分散熵偏均值与非线性模式分解的变工况轴承故障诊断方法 |
CN111259776B (zh) * | 2020-01-13 | 2023-04-18 | 浙江大学 | 一种基于同步平均主成分时频分析的确定性信号提取方法 |
CN111521112B (zh) * | 2020-04-23 | 2021-04-27 | 西安工业大学 | 一种傅里叶及窗口傅里叶变换的联合相位重构算法 |
CN112070842B (zh) * | 2020-07-28 | 2023-03-21 | 安徽农业大学 | 一种基于正交编码条纹的多摄像机全局标定方法 |
CN111985412B (zh) * | 2020-08-21 | 2022-08-05 | 西安交通大学 | 一种高压直流输电线路雷击干扰识别方法 |
CN112434708A (zh) * | 2020-11-18 | 2021-03-02 | 西安理工大学 | 一种极坐标二维s变换图像局部谱识别方法 |
CN113281809B (zh) * | 2021-04-29 | 2023-07-14 | 西安建筑科技大学 | 一种地震信号的谱分析方法 |
CN113029042B (zh) * | 2021-05-25 | 2021-08-03 | 四川大学 | 一种高温熔融态金属表面形貌动态测量装置及方法 |
CN115248549B (zh) * | 2022-01-12 | 2024-05-24 | 浙江理工大学 | 一种打散消除杂散频谱噪声的数字全息三维重建方法 |
CN114428282B (zh) * | 2022-01-26 | 2023-04-25 | 长安大学 | 一种基于去尺度s变换的地震信号时频变换方法 |
CN114608480A (zh) * | 2022-03-16 | 2022-06-10 | 合肥因赛途科技有限公司 | 基于移相条纹投影的位相-高度关系校准方法 |
CN115200797B (zh) * | 2022-09-19 | 2022-12-16 | 山东超华环保智能装备有限公司 | 一种用于零泄露阀的泄露检测系统 |
CN116330667B (zh) * | 2023-03-28 | 2023-10-24 | 云阳县优多科技有限公司 | 一种玩具的3d打印模型设计方法及系统 |
CN116402819B (zh) * | 2023-06-08 | 2023-09-12 | 季华实验室 | 类镜面缺陷检测方法、装置、电子设备及存储介质 |
CN116485690B (zh) * | 2023-06-20 | 2023-09-05 | 中国科学院苏州生物医学工程技术研究所 | 一种x射线光栅成像莫尔条纹漂移校准方法及装置 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP4398277B2 (ja) * | 2004-01-15 | 2010-01-13 | 株式会社山武 | 3次元計測装置、3次元計測方法及び3次元計測プログラム |
CN102032877A (zh) * | 2010-11-30 | 2011-04-27 | 东南大学 | 基于小波变换的三维测量方法 |
CN102169476A (zh) * | 2011-04-14 | 2011-08-31 | 哈尔滨工业大学 | 一种基于灰色理论的希尔伯特-黄变换端点效应抑制方法 |
CN102322822A (zh) * | 2011-08-08 | 2012-01-18 | 西安交通大学 | 一种三频彩色条纹投影三维测量方法 |
Family Cites Families (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP5194963B2 (ja) * | 2008-04-03 | 2013-05-08 | 株式会社ニコン | 波形解析装置、波形解析プログラム、干渉計装置、パターン投影形状測定装置、及び波形解析方法 |
TWI409502B (zh) * | 2009-01-23 | 2013-09-21 | Univ Nat Taipei Technology | 相位資訊擷取方法及其三維形貌量測系統 |
US8861833B2 (en) * | 2009-02-18 | 2014-10-14 | International Press Of Boston, Inc. | Simultaneous three-dimensional geometry and color texture acquisition using single color camera |
JP5375201B2 (ja) * | 2009-03-02 | 2013-12-25 | 株式会社豊田中央研究所 | 三次元形状測定方法及び三次元形状測定装置 |
JP5743419B2 (ja) * | 2010-04-18 | 2015-07-01 | 国立大学法人宇都宮大学 | 形状測定方法及び装置並びに歪み測定方法及び装置 |
-
2012
- 2012-01-19 CN CN201210016688.3A patent/CN102628676B/zh not_active Expired - Fee Related
- 2012-02-28 KR KR1020137005460A patent/KR101475382B1/ko not_active IP Right Cessation
- 2012-02-28 WO PCT/CN2012/071731 patent/WO2013107076A1/zh active Application Filing
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP4398277B2 (ja) * | 2004-01-15 | 2010-01-13 | 株式会社山武 | 3次元計測装置、3次元計測方法及び3次元計測プログラム |
CN102032877A (zh) * | 2010-11-30 | 2011-04-27 | 东南大学 | 基于小波变换的三维测量方法 |
CN102169476A (zh) * | 2011-04-14 | 2011-08-31 | 哈尔滨工业大学 | 一种基于灰色理论的希尔伯特-黄变换端点效应抑制方法 |
CN102322822A (zh) * | 2011-08-08 | 2012-01-18 | 西安交通大学 | 一种三频彩色条纹投影三维测量方法 |
Non-Patent Citations (6)
Title |
---|
JP特开2011-226871A 2011.11.10 |
JP特许4398277B2 2009.10.30 |
基于三色光投影的傅里叶变换三维测量法;黄昊等;《2010年全国模式识别学术会议论文集》;20110718;第613-618页 * |
朱林等.相位测量轮廓术中一种有效补偿相移误差的新算法.《东南大学学报》.2010,第40卷全文. |
相位测量轮廓术中一种有效补偿相移误差的新算法;朱林等;《东南大学学报》;20100930;第40卷;第302-308页 * |
黄昊等.基于三色光投影的傅里叶变换三维测量法.《2010年全国模式识别学术会议论文集》.2011,全文. |
Also Published As
Publication number | Publication date |
---|---|
CN102628676A (zh) | 2012-08-08 |
WO2013107076A1 (zh) | 2013-07-25 |
KR20130119910A (ko) | 2013-11-01 |
KR101475382B1 (ko) | 2014-12-22 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102628676B (zh) | 一种光学三维测量中的自适应窗口傅里叶相位提取法 | |
CN101608908B (zh) | 数字散斑投影和相位测量轮廓术相结合的三维数字成像方法 | |
KR101214702B1 (ko) | 웨이블릿 변환에 기반한 3차원 측정 방법 | |
CN102322822B (zh) | 一种三频彩色条纹投影三维测量方法 | |
CN103267496B (zh) | 一种基于小波变换的改进窗口傅里叶三维测量法 | |
CN102620685B (zh) | 一种基于史托克维尔变换的改进窗口傅里叶三维测量法 | |
CN107392875A (zh) | 一种基于k近邻域划分的点云数据去噪方法 | |
CN103454276B (zh) | 一种基于动态序列图像的织物形态风格评价方法 | |
CN103456015A (zh) | 基于最优分数域Gabor谱特征的SAR目标检测方法 | |
CN113379818B (zh) | 一种基于多尺度注意力机制网络的相位解析方法 | |
CN103454636A (zh) | 基于多像素协方差矩阵的差分干涉相位估计方法 | |
CN103325105A (zh) | 一种高精度合成孔径雷达图像自动配准方法及设备 | |
CN101540049A (zh) | 一种高光谱图像的端元提取方法 | |
CN103713287A (zh) | 一种基于互质多基线的高程重建方法及装置 | |
CN105043301A (zh) | 一种用于三维测量的光栅条纹相位求解方法 | |
CN102012466A (zh) | 数字x射线成像系统的噪声测量方法 | |
CN104730521A (zh) | 一种基于非线性优化策略的SBAS-DInSAR方法 | |
CN105353374B (zh) | 一种用于自旋目标的单频雷达成像方法 | |
CN111521112A (zh) | 一种傅里叶及窗口傅里叶变换的联合相位重构算法 | |
CN102819840B (zh) | 一种纹理图像的分割方法 | |
CN106931905A (zh) | 一种基于非线性优化的数字莫尔条纹相位提取方法 | |
CN105303538A (zh) | 一种基于nsct和pca的高斯噪声方差估计方法 | |
CN103310456B (zh) | 基于Gaussian-Hermite矩的多时相/多模态遥感图像配准方法 | |
CN113074668B (zh) | 一种基于新型2d复小波的三维面形测量方法 | |
CN103323844B (zh) | 一种多通道干涉合成孔径雷达高程重建方法及装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20140507 Termination date: 20170119 |
|
CF01 | Termination of patent right due to non-payment of annual fee |