CN104200474A - 一种获取物体变形量的数字图像分析方法 - Google Patents

Abstract

本发明公开了一种获取物体变形量的数字图像分析方法，步骤S1获得变形前和变形后的物体图像，步骤S2在两幅图像中同一物体相同部分处分别选取坐标相同的分析区域；步骤S3分别对两个区域内的像素点的灰度值进行快速傅里叶变换和数学处理，得到关于分析区域内像素点的灰度值、变形量以及关于数学处理中变量的函数；步骤S4为分别计算所述数学处理中变量取四个不同定值时获得的四个不同函数；步骤S5为对步骤S4中四个不同函数进行数学处理，从而获得第五函数；步骤S6对所述第五函数进行傅里叶变换，得到包含物体变形量的脉冲函数，计算以获得该脉冲函数中最高脉冲点的位置，该位置的坐标即为物体变形造成在X轴方向和Y轴方向上变形量。

Description

一种获取物体变形量的数字图像分析方法

技术领域

本发明属于工程测量领域，更具体地，涉及一种获取物体变形量的数字图像分析方法。

背景技术

测量材料的结构变形从而得到其力学性能一直是广大力学工作者所关心的重要课题。在材料的力学行为实验中通常将材料做成标准试样，再借助引伸计获取试样变形，根据材料变形可计算出材料的力学性能。早期使用机械杠杆式引伸计，现在通常用应变式引伸计，应变式引伸计上的敏感变形元件是由弹性材料制成的悬臂梁，该悬臂梁的自由端为一刀口，工作时其与被测件紧密固定，梁上粘有测量变形的应变片。

以上采用应变片的接触式测量方法存在如下问题：在一些织物，塑料等柔性材料上无法采用引伸计的测量方式；在高温高压等恶劣环境下，常规的应变片根本无法适应恶劣的使用环境，且高温会使构件发生蠕变，即使采用高温电阻应变计也会因受到高温影响可无法准确测量；引伸计上的悬臂梁刀口会因长久使用变钝，从而导致接触力下降而造成测试误差；对于测试橡胶这类大变形材料或金属断裂疲劳试验中，引伸计悬臂梁刀口会因为材料断裂时震动和冲击而严重损坏；对于一些微小试件时，没有合适标距的引伸计而无法测得变形。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种获取物体变形量的数字图像分析方法，其目的在于获取物体由于自身变形造成的位移量，无需接触构件，不对构件造成破坏，采用图像处理技术分析物体变形，可快速获取结果，由此解决目前方法中需要接触物体表面或者处理步骤耗费时间长的技术问题。

为实现上述目的，本发明提供了一种获取物体变形量的数字图像分析方法，用于获取物体由于自身变形造成的位移量，表面位移量即为物体由于变形或者移动造成的质点间的相对移动距离，其特征在于，包含以下步骤：

S1：利用拍照设备在同一位置拍摄同一物体变形前后的两幅图像，且该两幅图像至少包括该物体的至少一个相同部分；

S2：在所述两幅图像中同一物体相同部分的区域内分别随机选取一个分析区域，该两个分析区域的形状和包含的像素数量均相同，并且该两个区域在图像中像素的起始坐标和像素的终点坐标也分别相同；

S3：分别对所述步骤S2中两个区域内的所述像素点的灰度值进行快速傅里叶变换，分别得到表示物体变形前的第一变换结果和表示物体变形后的第二变换结果，由于物体变形后图像相对于变形前图像改变了变形造成的位移量，所以第二变换结果就是第一变化结果与位移量的函数；将第二变换结果乘以一个角度变量，进一步得到第三变换结果，将所述第一变换结果和所述第三变换结果相加，以得到换算结果，所述换算结果为矩阵，对该矩阵进行共轭相乘处理，得到关于变形前图像的所述分析区域内像素点的灰度值、变形造成的位移量以及所述角度变量的函数；

S4：分别计算所述角度变量取第一定值时得到的第一函数、角度变量取第二定值时得到的第二函数、角度变量取第三定值时得到的第三函数以及角度变量取第四定值时得到的第四函数；

S5：对所述步骤S4中的第一函数，第二函数，第三函数以及第四函数进行数学处理，以抵消掉物体变形前图像的灰度值在频域中共轭相乘所引起的一个系数，从而得到第五函数；

S6：对所述第五函数进行傅里叶变换，得到包含变形造成的位移量的脉冲函数，计算获得该脉冲函数中最高脉冲点的位置，该位置的坐标即为物体变形造成的位移量。

进一步的，所述步骤S4中第一定值、第二定值、第三定值以及第四定值分别为0、π/2、π以及3π/2。

进一步的，所述步骤S5中获得的第五函数I₅如下所示：

$\begin{array}{c}{I}_5(u,v)=\frac{I_1-I_3}{\sqrt{{(I_1-I_3)}^2+{(I_4-I_2)}^2}}\\ =\cos 2\mathrm{\π}(\frac{\mathrm{dxu}}{M}+\frac{\mathrm{dyv}}{N})\end{array}$

其中，I₁、I₂、I₃以及I₄分别为第一函数、第二函数、第三函数以及第四函数，u,v分别是其傅里叶变换后的坐标单位，M,N分别是选取的所述分析区域的x,y方向的像数数量，dx和dy分别为变形造成的位移量，π为数学中角度，等于180°。

进一步的，所述步骤S6中，所述脉冲函数及其换算过程为：

$\begin{array}{c}W(x,y)=\frac{1}{\mathrm{MN}}\underset{u=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{N-1}{\mathrm{\Σ}}}{I}_5(u,v)e^{-j2\mathrm{\π}(\frac{\mathrm{ux}}{M}+\frac{\mathrm{vy}}{N})}\\ =\frac{1}{\mathrm{MN}}\underset{u=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{N-1}{\mathrm{\Σ}}}\cos \left[2\mathrm{\π}(\frac{\mathrm{udx}}{M}+\frac{\mathrm{vdy}}{N})\right]e^{-j2\mathrm{\π}(\frac{\mathrm{ux}}{M}+\frac{\mathrm{vy}}{N})}\\ =\frac{1}{2}[\mathrm{\δ}(x+\mathrm{dx},y+\mathrm{dy})+\mathrm{\δ}(x-\mathrm{dx},y-\mathrm{dy})]\end{array}$

式中，M、N分别是选取的所述分析区域的x,y方向的像数数量，x和y分别为由图像平面位移和像素点的灰度值组成的时空区域内的坐标单位，dx和dy分别变形造成的x,y方向的表面位移量，u,v分别是傅里叶变换后的坐标单位，j为虚数单位，δ(x)为狄拉克函数，e为自然对数的底数。

本发明方法认为物体由于变形造成了表面位移，该表面位移量仅仅是由于物体变形造成的，因此该位移量较小。在同一位置拍摄物体变形前后表面的图像，该表面图像因为物体变形而发生了变化，其对应的像素点也发生了移动，根据灰度值相同，求解出具有相同灰度值的像素点发生的位移量，即为由于物体变形造成的表面位移。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

1、本发明采集原始数据的方法即拍摄变形物体前后的照片，方式简单，且对各种测试环境均可拍摄，不受测试环境影响，适用范围广，且拍摄照片无需接触构件，具有不破坏构件的优点。

2、本发明方法基于图像处理技术分析物体变形，采用数学中常见的傅里叶变换或者变速傅里叶变换进行处理，不需要进行搜寻过程，可以直接得到物体变形造成的位移量，可快速获取结果，节省时间。

附图说明

图1(a)和图1(b)分别是本发明实施例中表面喷有散斑物体相同部分变形前和变形后的图片；

图2是采用本发明实施例中方法获得的第五函数的条纹图；

图3是采用本发明实施例方法获得的步骤S6中的脉冲函数的频谱图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

实施例包含以下步骤：

S1：利用相机在同一位置拍摄表面具有散斑的物体变形前后的图像，如图1所示，图1(a)是变形前散斑部分的照片，图1(b)是变形后散斑部分的照片，该两幅图像包括该相同部分，即散斑部分。

S2：在两幅图像中散斑附近的区域内分别随机选取分析区域，该两个分析区域的形状相同，均为正方形，且两个区域包含的像素数量也相同，本实施例中，在变形前和变形后的的图像中选取的像素点起始坐标均为(350，150)，结束像素点坐标均为(450，250)，分析区域均为尺寸101×101。

S3：分别对两个区域内的像素点的灰度值进行快速傅里叶变换，分别得到表示物体变形前的第一变换结果和表示物体变形后的第二变换结果，第一变换结果即为指数函数F₀(u,v)，第二变换结果即为指数函数F₁(u,v)，两个指数函数分别如下所示：

$F_0(u,v)=\frac{1}{\mathrm{MN}}\underset{x=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{f}_0(x,y)e^{-j2\mathrm{\π}(\frac{\mathrm{xu}}{M}+\frac{\mathrm{yv}}{N})}$

$\begin{array}{c}{F}_1(u,v)=\frac{1}{\mathrm{MN}}\underset{x=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{y=0}{\overset{N-1}{\mathrm{\Σ}}}{f}_0(x+\mathrm{dx},y+\mathrm{dy})e^{-2\mathrm{\πj}(\frac{\mathrm{xu}}{M}+\frac{\mathrm{yv}}{N})}\\ =F_0(u,v)e^{2\mathrm{\πj}(\frac{\mathrm{udx}}{M}+\frac{\mathrm{vdy}}{N})}\end{array}$

式中，M、N分别是选取的分析区域的x,y方向的像数数量，x和y分别为由图像平面位移和像素点的灰度值组成的时空区域内的坐标，dx和dy分别是由于变形造成像素点其在x,y方向的位移，u,v分别是傅里叶变换后的坐标单位，j为虚数单位，e为自然对数的底数，π为数学中角度，等于180°，f₀(x,y)表示变形前图像灰度值与坐标间关系的函数，f₀(x+dx,y+dy)则表示变形后图像的灰度值与坐标间的关系函数。由于可认为变形前后位移是平移关系，所以，变形前后图像像素的灰度值存在X方向和Y方向的位移值，即分别为dx和dy，也即物体变形前后产生的表面位移量。

将变形后图像灰度值函数的快速傅里叶变换结果，也即第二变换结果F₁(u,v)乘上一个角度变量得到第三变换结果F₃(u,v)，如下式所示：

式中，M、N分别是选取的所述分析区域的x,y方向的像数数量，x和y分别为由于图像平面位移和像素点的灰度值组成的时空区域内的坐标，dx和dy分别是在x,y方向的位移，u,v分别是傅里叶变换后的坐标单位，j为虚数单位，为需要增加的位移变量。

由于变形造成的位移量是一个相对运动造成的结果，其依赖于变形前的图像的灰度值，之所以增加就是为了在后期的处理中利用这个增加的位移量将变形造成的位移量独立于变形前图像的灰度值，从而为求解变形造成的位移量提供可能。

F₃(u,v)是一个与F₀(u,v)，表面位移量dx和dy以及角度变量有关的函数，然后将这个函数与变形前图像的快速傅里叶变换结果F₀(u,v)相加，得到一个函数A_i(u,v)，该函数就是换算结果，具体过程如下式所示：

式中，M、N分别是选取的所述分析区域的x,y方向的像数数量，x和y分别为由于图像平面位移和像素点的灰度值组成的时空区域内的坐标，dx和dy分别是在x,y方向的位移，u,v分别是傅里叶变换后的坐标单位，j为虚数单位，为需要增加的位移变量。

将代表该换算结果的函数A_i(u,v)共轭相乘，得到关于变形前图像的分析区域内像素点的灰度值、变形造成的位移以及角度变量的函数I_i(u,v)，具体即得到一个与F₀(u,v)，表面位移量dx和dy以及变量中有关的函数，如下式所示：

其中，B₀(u,v)＝2F₀(u,v)F₀ ^*(u,v)，式中，M、N分别是选取的所述分析区域的x,y方向的像数数量，x和y分别为由于图像平面位移和像素点的灰度值组成的时空区域内的坐标，dx和dy分别是在x,y方向的位移，u,v分别是傅里叶变换后的坐标单位，j为虚数单位，为需要增加的位移变量，*表示其为所示函数的转置。

S4：分别计算所述变量取第一定值时得到的第一函数、变量取第二定值时得到的第二函数、变量取第三定值时得到的第三函数以及变量取第四定值时得到的第四函数，本实施例中第一定值、第二定值、第三定值以及第四定值分别为π/2,π和3π/2，得到第一函数、第二函数、第三函数以及第四函数，分别对应为I₁、I₂、I₃以及I₄，具体如下：

$I_1(u,v)=B_0(u,v)(1+\cos 2\mathrm{\π}(\frac{\mathrm{dxu}}{M}+\frac{\mathrm{dyv}}{N}))$

$I_2(u,v)=B_0(u,v)(1-\sin 2\mathrm{\π}(\frac{\mathrm{dxu}}{M}+\frac{\mathrm{dyv}}{N}))$

$I_3(u,v)=B_0(u,v)(1-\cos 2\mathrm{\π}(\frac{\mathrm{dxu}}{M}+\frac{\mathrm{dyv}}{N}))$

$I_4(u,v)=B_0(u,v)(1+\sin 2\mathrm{\π}(\frac{\mathrm{dxu}}{M}+\frac{\mathrm{dyv}}{N}))$

式中，式中，M、N分别是选取的所述分析区域的x,y方向的像数数量，x和y分别为由于图像平面位移和像素点的灰度值组成的时空区域内的坐标，dx和dy分别是在x,y方向的位移，u,v分别是傅里叶变换后的坐标单位，j为虚数单位，B₀(u,v)＝2F₀(u,v)F₀ ^*(u,v)。

S5：对步骤S4中的第一函数、第二函数、第三函数以及第四函数进行数学处理以抵消掉所述四个函数中的像素点的灰度值，就可以得到关于表面位移量dx和dy的第五函数，该第五函数为如下的三角函数：

$\begin{array}{c}{I}_5(u,v)=\frac{I_1-I_3}{\sqrt{{(I_1-I_3)}^2+{(I_4-I_2)}^2}}\\ =\cos 2\mathrm{\π}(\frac{\mathrm{dxu}}{M}+\frac{\mathrm{dyv}}{N})\end{array}$

其中，I₁、I₂、I₃以及I₄分别为第一函数、第二函数、第三函数以及第四函数，u,v分别是其傅里叶变换后的坐标单位，M,N分别是选取的所述分析区域的x,y方向的像数数量，dx和dy分别为变形造成的位移量，π为数学中角度，等于180°。

用软件matlab对三角函数I₅(u,v)绘图，得到关于表面位移量dx和dy的条纹图，如图2所示。

S6：最后，再对I₅(u,v)进行傅里叶变换，可得到去掉零频分量的脉冲函数W(x,y)，其最高脉冲对应的点的位置就是所求的表面位移量dx和dy，具体如下所示：

$\begin{array}{c}W(x,y)\\ =\frac{1}{\mathrm{MN}}\underset{u=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{N-1}{\mathrm{\Σ}}}{I}_5(u,v)e^{-j2\mathrm{\π}(\frac{\mathrm{ux}}{M}+\frac{\mathrm{vy}}{N})}\\ =\frac{1}{\mathrm{MN}}\underset{u=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{N-1}{\mathrm{\Σ}}}\cos \left[2\mathrm{\π}(\frac{\mathrm{udx}}{M}+\frac{\mathrm{vdy}}{N})\right]e^{-j2\mathrm{\π}(\frac{\mathrm{ux}}{M}+\frac{\mathrm{vy}}{N})}\\ =\frac{1}{2}[\mathrm{\δ}(x+\mathrm{dx},y+\mathrm{dy})+\mathrm{\δ}(x-\mathrm{dx},y-\mathrm{dy})]\end{array}$

式中，M、N分别是选取的所述分析区域的x,y方向的像数数量，x和y分别为由图像平面位移和像素点的灰度值组成的时空区域内的坐标单位，dx和dy分别变形造成的x,y方向的表面位移量，u,v分别是傅里叶变换后的坐标单位，j为虚数单位，δ(x)为狄拉克函数，e为自然对数的底数。

利用MAYLAB或其他软件编写程序，计算峰值对应的坐标点位置，就是所求的表面位移量dx和dy的具体数值。在本实施例中为dx和dy分别为9个像素，和7个像素，使用软件绘制出该脉冲函数，如图3所示。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种获取物体变形量的数字图像分析方法，用于获取物体由于自身变形造成的位移量，表面位移量即为物体由于变形或者移动造成的表面质点间的相对移动距离，其特征在于，包含以下步骤：

S1：利用拍照设备在同一位置拍摄同一物体变形前后的两幅图像；

S2：在所述两幅图像中同一物体相同部分的区域内分别随机选取一个分析区域，该两个分析区域的形状和包含的像素数量均相同，并且该两个区域在图像中像素的起始坐标和像素的终点坐标也分别相同；

S3：分别对所述步骤S2中两个区域内的所述像素点的灰度值进行快速傅里叶变换，分别得到表示物体变形前的第一变换结果和表示物体变形后的第二变换结果，所述第二变换结果为第一变换结果与表面位移量的函数，将第二变换结果乘以一个角度变量，进一步得到第三变换结果，将所述第一变换结果和所述第三变换结果相加，以得到换算结果，所述换算结果为矩阵，对该矩阵进行共轭相乘处理，得到关于变形前图像的所述分析区域内像素点的灰度值、表面位移量以及所述角度变量的函数；

S4：分别计算所述角度变量取第一定值时得到的第一函数、角度变量取第二定值时得到的第二函数、角度变量取第三定值时得到的第三函数以及角度变量取第四定值时得到的第四函数；

S5：对所述步骤S4中的第一函数，第二函数，第三函数以及第四函数进行数学处理，以抵消掉物体变形前图像的灰度值在频域中共轭相乘所引起的系数，以得到第五函数；

S6：对所述第五函数进行傅里叶变换，得到包含变形造成的位移量的脉冲函数，计算获得该脉冲函数中最高脉冲点的位置，该位置的坐标即为物体变形造成的位移量。

2.如权利要求1所述的一种获取物体变形量的数字图像分析方法，其特征在于，所述步骤S4中第一定值、第二定值、第三定值以及第四定值分别为0、π/2、π以及3π/2。

3.如权利要求1或2所述的一种获取物体变形量的方法，其特征在于，所述步骤S5中获得的第五函数I₅如下所示：

$\begin{array}{c}{I}_5(u,v)=\frac{I_1-I_3}{\sqrt{{(I_1-I_3)}^2+{(I_4-I_2)}^2}}\\ =\cos 2\mathrm{\π}(\frac{\mathrm{dxu}}{M}+\frac{\mathrm{dyv}}{N})\end{array}$

其中，I₁、I₂、I₃以及I₄分别为第一函数、第二函数、第三函数以及第四函数，u,v分别是其傅里叶变换后的坐标单位，M,N分别是选取的所述分析区域的x,y方向的像数数量，dx和dy分别为变形造成的x轴,y轴方向的位移量。

4.如权利要求1或2所述的一种获取物体变形量的方法，其特征在于，所述步骤S6中，所述脉冲函数及其换算过程为：

$\begin{array}{c}W(x,y)=\frac{1}{\mathrm{MN}}\underset{u=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{N-1}{\mathrm{\Σ}}}{I}_5(u,v)e^{-j2\mathrm{\π}(\frac{\mathrm{ux}}{M}+\frac{\mathrm{vy}}{N})}\\ =\frac{1}{\mathrm{MN}}\underset{u=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{v=0}{\overset{N-1}{\mathrm{\Σ}}}\cos \left[2\mathrm{\π}(\frac{\mathrm{udx}}{M}+\frac{\mathrm{vdy}}{N})\right]e^{-j2\mathrm{\π}(\frac{\mathrm{ux}}{M}+\frac{\mathrm{vy}}{N})}\\ =\frac{1}{2}[\mathrm{\δ}(x+\mathrm{dx},y+\mathrm{dy})+\mathrm{\δ}(x-\mathrm{dx},y-\mathrm{dy})]\end{array}$

式中，M、N分别是选取的所述分析区域的x轴,y轴方向的像数数量，x和y分别为由图像平面位移和像素点的灰度值组成的时空区域内的坐标单位，dx和dy分别变形造成的x轴,y轴方向的表面位移量，u,v分别是傅里叶变换后的坐标单位，j为虚数单位，δ(x)为狄拉克函数，e为自然对数的底数。