CN109801291A - 一种运动磨粒多表面三维形貌的获取方法 - Google Patents
一种运动磨粒多表面三维形貌的获取方法 Download PDFInfo
- Publication number
- CN109801291A CN109801291A CN201910185038.3A CN201910185038A CN109801291A CN 109801291 A CN109801291 A CN 109801291A CN 201910185038 A CN201910185038 A CN 201910185038A CN 109801291 A CN109801291 A CN 109801291A
- Authority
- CN
- China
- Prior art keywords
- formula
- abrasive grain
- point
- image
- abrasive
- 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.)
- Granted
Links
- 239000006061 abrasive grain Substances 0.000 title claims abstract description 129
- 238000000034 method Methods 0.000 title claims abstract description 58
- 230000033001 locomotion Effects 0.000 claims abstract description 46
- 239000002245 particle Substances 0.000 claims abstract description 34
- 239000000284 extract Substances 0.000 claims abstract description 13
- 238000013507 mapping Methods 0.000 claims abstract description 12
- 238000002156 mixing Methods 0.000 claims abstract description 4
- 239000011159 matrix material Substances 0.000 claims description 25
- 239000013598 vector Substances 0.000 claims description 15
- 238000013519 translation Methods 0.000 claims description 11
- 238000000605 extraction Methods 0.000 claims description 10
- 238000012545 processing Methods 0.000 claims description 10
- 238000001914 filtration Methods 0.000 claims description 7
- 230000003287 optical effect Effects 0.000 claims description 7
- 230000000007 visual effect Effects 0.000 claims description 7
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 6
- 238000001514 detection method Methods 0.000 claims description 6
- 230000000694 effects Effects 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 6
- 238000012216 screening Methods 0.000 claims description 6
- 230000009466 transformation Effects 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 4
- 238000005286 illumination Methods 0.000 claims description 4
- 230000000877 morphologic effect Effects 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000000227 grinding Methods 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 235000013399 edible fruits Nutrition 0.000 claims description 2
- 238000009499 grossing Methods 0.000 claims 1
- 238000004458 analytical method Methods 0.000 abstract description 13
- 238000005516 engineering process Methods 0.000 abstract description 10
- 230000003068 static effect Effects 0.000 abstract description 3
- 238000005211 surface analysis Methods 0.000 abstract 1
- 238000012544 monitoring process Methods 0.000 description 6
- 238000012423 maintenance Methods 0.000 description 3
- XEEYBQQBJWHFJM-UHFFFAOYSA-N Iron Chemical compound [Fe] XEEYBQQBJWHFJM-UHFFFAOYSA-N 0.000 description 2
- 238000005299 abrasion Methods 0.000 description 2
- 230000007812 deficiency Effects 0.000 description 2
- 230000008030 elimination Effects 0.000 description 2
- 238000003379 elimination reaction Methods 0.000 description 2
- 238000012876 topography Methods 0.000 description 2
- 230000002776 aggregation Effects 0.000 description 1
- 238000004220 aggregation Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 229910052742 iron Inorganic materials 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
Landscapes
- Length Measuring Devices By Optical Means (AREA)
Abstract
一种运动磨粒多表面三维形貌的获取方法，通过已有的运动磨粒视频采集平台，实现对运动磨粒的自动跟踪，获取同一运动磨粒在不同视角下的图像序列，之后利用SIFT算法提取磨粒表面的特征点，并在相邻帧定位特征点对应的位置，然后通过计算特征点的三维坐标实现空间点云模型的建立，同时通过SFS算法重建磨粒表面形貌，最后通过融合算法将SFS算法重建的磨粒表面形貌与空间点云模型进行融合，实现稠密重建，并通过纹理映射方法重构出具有真实感的运动磨粒多表面的三维形貌；本发明有效地将自动跟踪技术与三维重建技术相结合，解决了目前磨粒分析技术中只能从单一表面分析静态磨粒的问题，实现了磨粒的三维表面全信息精准提取。
Description
技术领域
本发明属于机器磨损状态监测领域的磨粒分析技术，特别涉及一种运动磨粒多表面三维形貌的获取方法。
背景技术
摩擦副之间的磨损是机械设备产生故障的主要原因之一，据统计，由磨损引起的故障大约占机械设备故障的75％。磨损状态监测通过分析机械设备的运行状态，为机械设备提供故障预测，维护策略，可以提高机械设备的安全性和可靠性，充分延长服役期限，使机械设备的维护成本与停机损失达到最小。随着现代化工业的高速发展，对机械设备的安全性和可靠性的要求越来越高，因此对机械设备进行磨损状态监测具有重要意义。
磨粒作为磨损过程中的产物，携带着丰富的磨损信息，其数量、大小、颜色和形貌等特征表征了磨损的微观机制，反映了机械设备的磨损状态。因此磨粒分析技术是机器磨损状态监测领域的重要技术，为机械设备的状态监测与维修决策提供有效依据。根据获取磨粒特征参数的维度，磨粒分析技术分为二维磨粒分析与三维磨粒分析。二维磨粒分析的主要应用是铁谱图像处理，该方法能够提取磨粒的平面二维特征，如浓度、等效尺寸、形状参数等，但不能获取磨粒的空间特征信息，导致其无法反映具有不规则性磨粒的二维图像。三维磨粒分析是指借助于激光扫描共聚焦显微镜(LSCM)、原子力显微镜(AFM)等仪器获取磨粒的三维信息。该方法操作复杂，仪器昂贵，使用范围有限。此外，无论是二维磨粒分析还是三维磨粒分析技术，在目前的应用中总存在着以下的不足：①只能对静止条件下的磨粒进行分析，然而在静止条件下获取的磨粒特征参数会因拍摄磨粒图像的条件和磨粒的聚集等原因产生不同的结果。②只能从单一视角提取磨粒的特征，然而磨粒的形状具有不规则性，单一的表面不能真实反映磨粒的性质。
总体而言，通过磨粒分析技术实现磨损状态表征已成为机械设备磨损状态监测的重要手段，已经在工业领域如航空、船舶、风力发电等方面获得了广泛的应用，但由于其方法存在的不足，如只能对静态磨粒分析、只能提取单一表面的磨粒特征信息，而无法全面地、准确地提取磨粒的信息。
发明内容
为了解决上述现有技术的缺陷，本发明的目的在于提供一种运动磨粒多表面三维形貌的获取方法，通过已有的运动磨粒多视角视频采集平台，使磨粒在流道内作翻滚运动，利用卡尔曼滤波法实现目标运动磨粒的自动跟踪，并基于关键帧提取模型提取磨粒有效多视角图像序列；之后为消除磨粒图像的形变，利用Matlab Camera Calibration Toolbox标定摄像机和校正磨粒图像；然后利用SIFT算法提取与匹配磨粒表面的特征点，再通过求解投影矩阵计算出各个特征点的三维坐标，建立空间点云模型；同时，通过阴影恢复形状法求解出磨粒表面的相对高度，恢复磨粒表面的形貌；最后通过融合算法将空间点云模型与SFS恢复的磨粒表面形貌进行融合，实现稠密重建，并通过纹理映射方法重构出具有真实感的运动磨粒多表面的三维形貌。
为实现上述目的，本发明采用如下技术方案：
一种运动磨粒多表面三维形貌的获取方法，包括以下步骤：
步骤1：通过多视角运动磨粒视频采集平台，获取磨粒的运动视频，利用卡尔曼滤波法对其进行分析，实现运动磨粒的自动跟踪，并基于关键帧提取模型选择磨粒有效多视角图像序列；
步骤2：对摄像机进行标定，利用Matlab相机标定工具箱Camera CalibrationToolbox校准后得到的相机内外参数对步骤1中获得的图像序列进行校正；
步骤3：利用SIFT算法提取并匹配步骤2中校正后的运动磨粒图像中的特征点，然后通过求解图像的投影矩阵，结合最小二乘法计算出每对匹配点所对应的空间坐标，重建出空间点云；
步骤4：利用SFS算法计算磨粒表面的相对高度，恢复磨粒表面的三维形貌；
步骤5：利用Power Curst算法对点云进行网格化处理，然后通过融合算法将SFS算法重建的磨粒表面形貌与网格化后的空间点云进行融合，实现稠密重建；
步骤6：纹理映射，利用OpenGL纹理映射方法将原始图像的纹理映射到步骤5中稠密重建的结果上，获得运动磨粒的多表面的三维形貌。
步骤1中所述的利用卡尔曼滤波法对其进行分析，实现运动磨粒的自动跟踪，并基于关键帧提取模型选择磨粒有效多视角图像序列，具体步骤如下：
步骤1.1：采用Matlab的视频工具箱中VideoRader函数按帧读取视频；
步骤1.2：根据磨粒的质心与形态参数能反映磨粒的实时运动的原理，选取磨粒的状态变量如公式(1)所示；
公式(1)：
式中x,y,l,w,A分别表示磨粒的质心横坐标、纵坐标、长轴、短轴、观测视角面积，vx,vy,vl,vw,vA分别表示上述值的变化速度；
步骤1.3：记当前时刻目标的状态向量为其卡尔曼预测为uk表示外部输入，将预测位置与外界输入相结合，得到当前状态如公式(2)所示；
公式(2)：
步骤1.4：记当前预测结果的不确定性为将之前的不确定性与本次的不确定性Q进行求和，如公式(3)所示；
公式(3)：
步骤1.5：记当前预测的卡尔曼增益亦即权重为Kk，由预测结果的不确定性和观测结果的不确定性R计算Kk，如公式(4)所示；
公式(4)：
步骤1.6：记当前的观测结果为zk，则当前的观测结果是预测结果和观测结果的加权平均所得到的如公式(5)所示；
公式(5)：
步骤1.7：记本次状态估计的不确定性为Pk，利用公式(6)更新Pk；
公式(6)：
公式(2)至公式(6)中A表示状态转移矩阵，B表示外部输入系数，H表示单位矩阵；
步骤1.8：利用Matlab对步骤1.3.2至步骤1.3.7进行编程，实现磨粒的自动跟踪；
步骤1.9：关键帧筛选方法，基于颜色、形状、特征点图像不变量构建颜色特征、形状特征和特征点对特征限制性因子，进而建立关键帧筛选模型，去除图像序列中冗余图像，提取出被跟踪的磨粒的有效图像序列；
所述步骤2包括如下具体步骤：
步骤2.1：首先对棋盘格图像进行拍照，拍照时摄像机位置固定，通过旋转和平移棋盘格在不同角度下不少于3幅棋盘格图像；
步骤2.2：启动Matlab相机标定工具箱Camera Calibration Toolbox，载入步骤2.1中拍摄棋盘格图像；
步骤2.3：对步骤2.1中所拍摄的棋盘格图像进行角点检测，点击相机标定工具箱上Extract grid cornder按钮，输入棋盘格单元方格的尺寸，自动提取角点；然后依次循环操作，对每一幅棋盘格图像进行角点检测；
步骤2.4：点击相机标定工具箱上的Calibration按钮进行标定，最终得到摄像机内外参数与畸变参数，完成摄像机的标定过程；
步骤2.5：利用步骤2.4中得到的标定结果，对运动磨粒多视角图像进行校正。
所述步骤3包括如下具体步骤：
步骤3.1：利用SIFT算法对磨粒特征点进行提取，首先利用高斯函数进行尺度空间极值点的提取，然后利用2×2的Hessian矩阵剔除半圆效应较强的极值点，最后未剔除的点即为磨粒的特征点；
步骤3.2：通过公式(7)和公式(8)计算每个特征点在规定邻域内的梯度和方向，并以梯度直方图的形式表示；然后根据梯度直方图的信息计算出每个特征点方向，尺度和位置信息，建立128维的SIFT特征点描述子；
公式(7)：
公式(8)：
公式(7)和公式(8)中L(x,y)表示高斯尺度空间；
步骤3.3：以欧式距离作为相似性度量准则对SIFT特征点进行匹配，利用公式(9)计算步骤3.1中提取的磨粒特征点的欧式距离：
公式(9)：
式中siftp、siftq分别表示待匹配特征点p、q的128维描述子；
规定匹配度阈值T为0.75，如果最小欧式距离min(d)与次小欧式距离secmin(d)的比值满足公式(10)，则特征点匹配成功；
公式(10)：
步骤3.4：求解投影矩阵M，获取磨粒图像时摄像机保持不动、磨粒运动，于是通过所假设磨粒不动，摄像机经过旋转与平移拍摄磨粒的另一视角的图像，进而求得匹配点的空间坐标；磨粒图像的投影矩阵定义如公式(11)所示；
公式(11)：
式中K为3×4摄像机内参矩阵，R为3×3旋转矩阵，t为3×1平移矢量；
旋转矩阵R与平移矢量t均可由奇异值分解法求解，故投影矩阵已被求解；
步骤3.5：计算匹配点对所对应的空间坐标；设M1j和M2j(j＝1,2,3)分别是M1、M2的3个行向量，用齐次坐标表示的空间三维点w＝(Xw,Yw,Zw,1)T在两幅图像上的投影分别是p1＝(u1,v1,1)T和p2＝(u2,v2,1)T，则由摄像机的理论成像模型，对每一对匹配点，可得：
公式(12)：
通过最小二乘法求得每对匹配点所对应的空间坐标；重建出空间点云。
所述步骤4包括如下具体步骤：
步骤4.1：对磨粒图像进行预处理，将磨粒图像进行颜色空间转换，即由RGB颜色空间转换为HSI颜色空间，并提取HSI颜色空间的亮度分量I，如公式(13)所示；另外，对磨粒图像进行平滑处理消除噪声；
公式(13)：
步骤4.2：选择磨粒的光照模型为朗伯模型；
步骤4.3：计算光源的方向，直接估计光源的世界坐标(以摄像机镜头的光心为原点，光轴为z轴的坐标系)，并进行单位化处理，即光源方向的单位向量s：
公式(14)：s＝(sx,sy,sz)
步骤4.4：利用Pentland方法对磨粒表面相对高度进行求解，根据公式(15)计算出高度图，然后进行傅里叶反变换即可计算出磨粒表面的相对高度，完成SFS算法恢复磨粒表面形貌，
公式(15)：FE(ω1,ω2)＝FZ(ω1,ω2)(-iω1)sx+FZ(ω1,ω2)(-iω2)sy
式中：FE(ω1,ω2)表示图像强度E(x,y)的傅里叶变换，FZ(ω1,ω2)表示高度函数Z(x,y)的傅里叶变换。
所述步骤5包括如下具体步骤：
步骤5.1：为建立空间点云模型中各个特征点之间的关系，采用Power Curst算法对点云模型进行网格化处理；
步骤5.2：进行中点加密，以步骤5.1中建立的网格每条边的中间点作为加密点，加密点的高度值根据SFS算法求解出的磨粒表面的高度与点云模型中特征点的二维位置计算而来，对网格化后的点云模型每条边均进行加密，加密完成后可再次进行，直到稠密重建效果满足要求即可。
本发明应用于机械设备磨损状态监测领域，具有以下有益效果：
(1)本发明将基于卡尔曼滤波法的自动跟踪方法与三维重建技术相结合，解决了磨粒分析技术中只能分析静态磨粒与只能分析单一磨粒表面的问题，实现了运动磨粒多表面三维形貌的获取，适用于磨粒分析技术中所有磨粒三维形貌的获取。
(2)本发明结合磨粒形态参数，利用卡尔曼滤波法可准确的预测目标运动磨粒所在区域，通过特征匹配自动确定目标运动磨粒，解决了以往运动速度不同造成多磨粒跟踪错误问题，为实现运动磨粒三维形貌的获取提供了运动磨粒不同视角下的图像序列。
(3)本发明采用SIFT算法提取与匹配运动磨粒图像的特征点对，进而通过求解投影矩阵计算特征点的空间坐标，建立了运动磨粒多表面的空间点云模型，然后将基于SFS算法获取的磨粒表面形貌与空间点云模型进行融合，获取了磨粒表面形貌的精细模型，提高了磨粒表面三维形貌的精度。
附图说明
图1为运动磨粒多表面三维形貌的获取方法的总体流程图。
图2为运动磨粒视频采集平台。
图3为基于卡尔曼滤波法自动跟踪流程图。
图4为目标运动磨粒多视角图像序列。
图5(a)为目标磨粒特征点匹配示意图。
图5(b)为目标磨粒空间点云模型。
图5(c)为基于SFS法恢复的磨粒表面三维形貌。
图5(d)为空间点云模型网格化图。
图6(a)为运动磨粒二维图。
图6(b)为运动磨粒多表面三维形貌重构结果俯视图。
图6(c)为运动磨粒多表面三维形貌重构结果斜视图。
具体实施方式
下面结合附图，对本发明作进一步说明。
参照图1，一种运动磨粒多表面三维形貌的获取方法，包括以下步骤：
步骤1：获取运动磨粒多表面三维形貌的首要步骤是采集运动磨粒多视角图像，通过多视角运动磨粒视频采集平台，使磨粒在特制的流道中运动翻滚，图像采集设备采集磨粒运动的视频，然后利用卡尔曼滤波方法实现磨粒的自动跟踪，并基于关键帧提取模型选择出目标磨粒的多视角图像序列。
所述步骤1包括如下具体步骤：
步骤1.1：建立运动磨粒多视角视频采集装置，如图2所示，该装置由二部分组成，一是含磨粒的油回路包括油样、数字微量泵、管道和特制的流道，二是磨粒运动视频采集部分，包括图像采集设备，光源和计算机。
步骤1.2：启动数字微量泵，磨粒在流道内作翻滚运动，图像采集设备拍摄磨粒运动的视频，传送给计算机。
步骤1.3：对步骤1.2中采集到的视频中的运动磨粒进行跟踪，跟踪流程图如图3所示，具体实现步骤如下：
步骤1.3.1：采用Matlab的视频工具箱中VideoRader函数按帧读取视频。
步骤1.3.2：根据磨粒的质心与形态参数能反映磨粒的实时运动的原理，选取磨粒的状态变量如公式(1)所示。
公式(1)：
式中x,y,l,w,A分别表示磨粒的质心横坐标、纵坐标、长轴、短轴、观测视角面积，vx,vy,vl,vw,vA分别表示上述值的变化速度。
步骤1.3.3：记当前时刻目标的状态向量为其卡尔曼预测为uk表示外部输入，将预测位置与外界输入相结合，得到当前状态如公式(2)所示。
公式(2)：
步骤1.3.4：记当前预测结果的不确定性为将之前的不确定性与本次的不确定性Q进行求和，如公式(3)所示。
公式(3)：
步骤1.3.5：记当前预测的卡尔曼增益(亦即权重)为Kk，由预测结果的不确定性和观测结果的不确定性R计算Kk，如公式(4)所示。
公式(4)：
步骤1.3.6：记当前的观测结果为zk，则当前的观测结果是预测结果和观测结果的加权平均所得到的如公式(5)所示。
公式(5)：
步骤1.3.7：记本次状态估计的不确定性为Pk，利用公式(6)更新Pk。
公式(6)：
公式(2)至公式(6)中A表示状态转移矩阵，B表示外部输入系数，H表示单位矩阵。
步骤1.3.8：利用Matlab对步骤1.3.2至步骤1.3.7进行编程，实现磨粒的自动跟踪。
步骤1.4：关键帧筛选方法，基于颜色、形状、特征点图像不变量构建颜色特征、形状特征和特征点对特征限制性因子，进而建立关键帧筛选模型，去除图像序列中冗余图像，提高重构效率。提取出被跟踪的磨粒的有效图像序列，如图4所示。
步骤2：由于本发明图像采集设备采用的是CCD摄像机，CCD摄像机拍摄的过程实际上是一个光学成像的过程，由于摄像机镜头内的透镜为凸透镜，只有穿过镜头中心部分的光线是直线传播的，因此拍摄的图像会出现不同程度的形变。为了提高获取的磨粒三维形貌的精度，需要对摄像机进行标定，对原始图像进行校正。
所述步骤2包括如下具体步骤：
步骤2.1：首先对棋盘格图像进行拍照，拍照时摄像机位置固定，通过旋转和平移棋盘格在不同角度下拍摄8幅棋盘格图像。
步骤2.2：启动Matlab相机标定工具箱Camera Calibration Toolbox，载入步骤2.1中拍摄的8幅棋盘格图像。
步骤2.3：对步骤2.1中所拍摄的8幅棋盘格图像进行角点检测。点击相机标定工具箱上Extract grid cornder按钮，输入棋盘格单元方格的尺寸，自动提取角点。然后依次循环操作，对每一幅棋盘格图像进行角点检测。
步骤2.4：点击相机标定工具箱上的Calibration按钮进行标定，最终得到摄像机内外参数与畸变参数，完成摄像机的标定过程。
步骤2.5：利用步骤2.4中得到的标定结果，对运动磨粒多视角图像进行校正。
步骤3：建立磨粒的空间点云模型是获取运动磨粒多表面的三维形貌的基础。磨粒的表面形貌信息在图像中以二维特征点的形式存在，只需在多幅图像中获得对应特征点的二维坐标，然后计算出特征点的三维坐标，即可建立空间点云模型。
所述步骤3包括如下具体步骤：
步骤3.1：利用SIFT算法对磨粒特征点进行提取，首先利用高斯函数进行尺度空间极值点的提取，然后利用2×2的Hessian矩阵剔除半圆效应较强的极值点，最后未剔除的点即为磨粒的特征点。
步骤3.2：通过公式(7)和公式(8)计算每个特征点在规定邻域内的梯度和方向，并以梯度直方图的形式表示。然后根据梯度直方图的信息计算出每个特征点方向，尺度和位置信息，建立128维的SIFT特征点描述子。
公式(7)：
公式(8)：
公式(7)和公式(8)中L(x,y)表示高斯尺度空间。
步骤3.3：以欧式距离作为相似性度量准则对SIFT特征点进行匹配，利用公式(9)计算步骤3.1中提取的磨粒特征点的欧式距离：
公式(9)：
式中siftp、siftq分别表示待匹配特征点p、q的128维描述子。
规定匹配度阈值T为0.75，如果最小欧式距离min(d)与次小欧式距离secmin(d)的比值满足公式(10)，则特征点匹配成功。以相邻两帧图像为例，匹配的特征点如图5(a)所示。
公式(10)：
步骤3.4：求解投影矩阵M，以步骤3.2中相邻两帧图像为例进行说明。获取磨粒图像时摄像机保持不动、磨粒运动，于是通过所假设磨粒不动，摄像机经过旋转与平移拍摄磨粒的另一视角的图像，进而求得匹配点的空间坐标。磨粒图像的投影矩阵定义如公式(11)所示。
公式(11)：
式中K为3×4摄像机内参矩阵(在步骤2中已求出)，R为3×3旋转矩阵，t为3×1平移矢量。
旋转矩阵R与平移矢量t均可由奇异值分解法求解，故投影矩阵已被求解。
步骤3.5：计算匹配点对所对应的空间坐标。设M1j和M2j(j＝1,2,3)分别是M1、M2的3个行向量，用齐次坐标表示的空间三维点w＝(Xw,Yw,Zw,1)T在两幅图像上的投影分别是p1＝(u1,v1,1)T和p2＝(u2,v2,1)T，则由摄像机的理论成像模型，对每一对匹配点，可得：
公式(12)：
通过最小二乘法求得每对匹配点所对应的空间坐标。图5(b)所示为得到的磨粒空间点云模型。
步骤4：利用阴影恢复形状法求解磨粒表面的相对高度，恢复磨粒表面的形貌。
所述步骤4包括如下具体步骤：
步骤4.1：对磨粒图像进行预处理，因为SFS算法需要的磨粒图像信息为亮度信息，所以需要将磨粒图像进行颜色空间转换，即由RGB颜色空间转换为HSI颜色空间，并提取HSI颜色空间的亮度分量I，如公式(13)所示。另外，在磨粒运动视频拍摄过程中总存在各种各样的噪声，因此还需要对磨粒图像进行平滑处理消除噪声。
公式(13)：
步骤4.2：光照模型的选择，选择磨粒的光照模型为朗伯模型。
步骤4.3：计算光源的方向，因为磨粒相对于光源十分的微小，所以直接估计光源的世界坐标(以摄像机镜头的光心为原点，光轴为z轴的坐标系)，并进行单位化处理，即光源方向的单位向量s：
公式(14)：s＝(sx,sy,sz)
步骤4.4：利用Pentland方法对磨粒表面相对高度进行求解，根据公式(15)可计算出高度图，然后进行傅里叶反变换即可计算出磨粒表面的相对高度，完成SFS算法恢复磨粒表面形貌，如图5(c)所示。
公式(15)：FE(ω1,ω2)＝FZ(ω1,ω2)(-iω1)sx+FZ(ω1,ω2)(-iω2)sy
式中：FE(ω1,ω2)表示图像强度E(x,y)的傅里叶变换，FZ(ω1,ω2)表示高度函数Z(x,y)的傅里叶变换。
步骤5：通过融合空间点云模型与阴影恢复形状算法恢复的磨粒表面，实现稠密重建，然后进行纹理映射，最终重构出运动磨粒多表面三维形貌。
所述步骤5包括如下具体步骤：
步骤5.1：为建立空间点云模型中各个特征点之间的关系，需对点云模型进行网格化处理，本发明采用的算法是Power Curst算法，效果如图5(d)所示。
步骤5.2：进行中点加密，以步骤5.1中建立的网格每条边的中间点作为加密点，加密点的高度值根据步骤4中SFS算法求解出的磨粒表面的高度与点云模型中特征点的二维位置计算而来。对网格化后的点云模型每条边均进行加密，加密完成后可再次进行，直到稠密重建效果达到最好为止。
(6)为使本发明获得的运动磨粒多表面三维形貌具有真实感，还需进行纹理映射，利用OpenGL纹理映射方法将原始图像的纹理映射到步骤(5)中稠密重建的结果上，实现磨粒的三维模型纹理的获取，完成运动磨粒的多表面三维形貌的获取。参照图6，将图6(a)所示的磨粒，进行以上处理后得到如图6(b)、(c)所示的三维形貌重构结果图。
Claims (6)
1.一种运动磨粒多表面三维形貌的获取方法，其特征在于，包括以下步骤：
步骤1：通过多视角运动磨粒视频采集平台，获取磨粒的运动视频，利用卡尔曼滤波法对其进行分析，实现运动磨粒的自动跟踪，并基于关键帧提取模型选择磨粒有效多视角图像序列；
步骤2：对摄像机进行标定，利用Matlab相机标定工具箱Camera Calibration Toolbox校准后得到的相机内外参数对步骤1中获得的图像序列进行校正；
步骤3：利用SIFT算法提取并匹配步骤2中校正后的运动磨粒图像中的特征点，然后通过求解图像的投影矩阵，结合最小二乘法计算出每对匹配点所对应的空间坐标，重建出空间点云；
步骤4：利用SFS算法计算磨粒表面的相对高度，恢复磨粒表面的三维形貌；
步骤5：利用Power Curst算法对点云进行网格化处理，然后通过融合算法将SFS算法重建的磨粒表面形貌与网格化后的空间点云进行融合，实现稠密重建；
步骤6：纹理映射，利用OpenGL纹理映射方法将原始图像的纹理映射到步骤5中稠密重建的结果上，获得运动磨粒的多表面的三维形貌。
2.根据权利要求1所述的一种运动磨粒多表面三维形貌的获取方法，其特征在于，步骤1中所述的利用卡尔曼滤波法对其进行分析，实现运动磨粒的自动跟踪，并基于关键帧提取模型选择磨粒有效多视角图像序列，具体如下：
步骤1.1：采用Matlab的视频工具箱中VideoRader函数按帧读取视频；
步骤1.2：根据磨粒的质心与形态参数能反映磨粒的实时运动的原理，选取磨粒的状态变量如公式(1)所示；
公式(1)：
式中x,y,l,w,A分别表示磨粒的质心横坐标、纵坐标、长轴、短轴、观测视角面积，vx,vy,vl,vw,vA分别表示上述值的变化速度；
步骤1.3：记当前时刻目标的状态向量为其卡尔曼预测为uk表示外部输入，将预测位置与外界输入相结合，得到当前状态如公式(2)所示；
公式(2)：
步骤1.4：记当前预测结果的不确定性为将之前的不确定性与本次的不确定性Q进行求和，如公式(3)所示；
公式(3)：
步骤1.5：记当前预测的卡尔曼增益亦即权重为Kk，由预测结果的不确定性和观测结果的不确定性R计算Kk，如公式(4)所示；
公式(4)：
步骤1.6：记当前的观测结果为zk，则当前的观测结果是预测结果和观测结果的加权平均所得到的如公式(5)所示；
公式(5)：
步骤1.7：记本次状态估计的不确定性为Pk，利用公式(6)更新Pk；
公式(6)：
公式(2)至公式(6)中A表示状态转移矩阵，B表示外部输入系数，H表示单位矩阵；
步骤1.8：利用Matlab对步骤1.3.2至步骤1.3.7进行编程，实现磨粒的自动跟踪；
步骤1.9：关键帧筛选方法，基于颜色、形状、特征点图像不变量构建颜色特征、形状特征和特征点对特征限制性因子，进而建立关键帧筛选模型，去除图像序列中冗余图像，提取出被跟踪的磨粒的有效图像序列。
3.根据权利要求1所述的一种运动磨粒多表面三维形貌的获取方法，其特征在于，所述步骤2包括如下具体步骤：
步骤2.1：首先对棋盘格图像进行拍照，拍照时摄像机位置固定，通过旋转和平移棋盘格在不同角度下不少于3幅棋盘格图像；
步骤2.2：启动Matlab相机标定工具箱Camera Calibration Toolbox，载入步骤2.1中拍摄棋盘格图像；
步骤2.3：对步骤2.1中所拍摄的棋盘格图像进行角点检测，点击相机标定工具箱上Extract grid cornder按钮，输入棋盘格单元方格的尺寸，自动提取角点；然后依次循环操作，对每一幅棋盘格图像进行角点检测；
步骤2.4：点击相机标定工具箱上的Calibration按钮进行标定，最终得到摄像机内外参数与畸变参数，完成摄像机的标定过程；
步骤2.5：利用步骤2.4中得到的标定结果，对运动磨粒多视角图像进行校正。
4.根据权利要求1所述的一种运动磨粒多表面三维形貌的获取方法，其特征在于，所述步骤3包括如下具体步骤：
步骤3.1：利用SIFT算法对磨粒特征点进行提取，首先利用高斯函数进行尺度空间极值点的提取，然后利用2×2的Hessian矩阵剔除半圆效应较强的极值点，最后未剔除的点即为磨粒的特征点；
步骤3.2：通过公式(7)和公式(8)计算每个特征点在规定邻域内的梯度和方向，并以梯度直方图的形式表示；然后根据梯度直方图的信息计算出每个特征点方向，尺度和位置信息，建立128维的SIFT特征点描述子；
公式(7)：
公式(8)：
公式(7)和公式(8)中L(x,y)表示高斯尺度空间；
步骤3.3：以欧式距离作为相似性度量准则对SIFT特征点进行匹配，利用公式(9)计算步骤3.1中提取的磨粒特征点的欧式距离：
公式(9)：
式中siftp、siftq分别表示待匹配特征点p、q的128维描述子；
规定匹配度阈值T为0.75，如果最小欧式距离min(d)与次小欧式距离secmin(d)的比值满足公式(10)，则特征点匹配成功；
公式(10)：
步骤3.4：求解投影矩阵M，获取磨粒图像时摄像机保持不动、磨粒运动，于是通过所假设磨粒不动，摄像机经过旋转与平移拍摄磨粒的另一视角的图像，进而求得匹配点的空间坐标；磨粒图像的投影矩阵定义如公式(11)所示；
公式(11)：
式中K为3×4摄像机内参矩阵，R为3×3旋转矩阵，t为3×1平移矢量；
旋转矩阵R与平移矢量t均可由奇异值分解法求解，故投影矩阵已被求解；
步骤3.5：计算匹配点对所对应的空间坐标；设M1j和M2j(j＝1,2,3)分别是M1、M2的3个行向量，用齐次坐标表示的空间三维点w＝(Xw,Yw,Zw,1)T在两幅图像上的投影分别是p1＝(u1,v1,1)T和p2＝(u2,v2,1)T，则由摄像机的理论成像模型，对每一对匹配点，可得：
公式(12)：
通过最小二乘法求得每对匹配点所对应的空间坐标；重建出空间点云。
5.根据权利要求1所述的一种运动磨粒多表面三维形貌的获取方法，其特征在于，所述步骤4包括如下具体步骤：
步骤4.1：对磨粒图像进行预处理，将磨粒图像进行颜色空间转换，即由RGB颜色空间转换为HSI颜色空间，并提取HSI颜色空间的亮度分量I，如公式(13)所示；另外，对磨粒图像进行平滑处理消除噪声；
公式(13)：
步骤4.2：选择磨粒的光照模型为朗伯模型；
步骤4.3：计算光源的方向，直接估计光源的世界坐标(以摄像机镜头的光心为原点，光轴为z轴的坐标系)，并进行单位化处理，即光源方向的单位向量s：
公式(14)：s＝(sx,sy,sz)
步骤4.4：利用Pentland方法对磨粒表面相对高度进行求解，根据公式(15)计算出高度图，然后进行傅里叶反变换即可计算出磨粒表面的相对高度，完成SFS算法恢复磨粒表面形貌，
公式(15)：FE(ω1,ω2)＝FZ(ω1,ω2)(-iω1)sx+FZ(ω1,ω2)(-iω2)sy
式中：FE(ω1,ω2)表示图像强度E(x,y)的傅里叶变换，FZ(ω1,ω2)表示高度函数Z(x,y)的傅里叶变换。
6.根据权利要求1所述的一种运动磨粒多表面三维形貌的获取方法，其特征在于，所述步骤5包括如下具体步骤：
步骤5.1：为建立空间点云模型中各个特征点之间的关系，采用Power Curst算法对点云模型进行网格化处理；
步骤5.2：进行中点加密，以步骤5.1中建立的网格每条边的中间点作为加密点，加密点的高度值根据SFS算法求解出的磨粒表面的高度与点云模型中特征点的二维位置计算而来，对网格化后的点云模型每条边均进行加密，加密完成后可再次进行，直到稠密重建效果满足要求即可。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910185038.3A CN109801291B (zh) | 2019-03-12 | 2019-03-12 | 一种运动磨粒多表面三维形貌的获取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910185038.3A CN109801291B (zh) | 2019-03-12 | 2019-03-12 | 一种运动磨粒多表面三维形貌的获取方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109801291A true CN109801291A (zh) | 2019-05-24 |
| CN109801291B CN109801291B (zh) | 2021-05-04 |
Family
ID=66562727
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910185038.3A Active CN109801291B (zh) | 2019-03-12 | 2019-03-12 | 一种运动磨粒多表面三维形貌的获取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109801291B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115326656A (zh) * | 2022-10-14 | 2022-11-11 | 西南交通大学 | 交通土建用颗粒材料松铺层颗粒粒径及级配无损测量方法 |
| CN115578523A (zh) * | 2022-11-21 | 2023-01-06 | 慧医谷中医药科技(天津)股份有限公司 | 一种多角度图像融合的舌部三维建模方法及系统 |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070084133A1 (en) * | 2005-10-18 | 2007-04-19 | 3M Innovative Properties Company | Agglomerate abrasive grains and methods of making the same |
| CN102890039A (zh) * | 2011-07-22 | 2013-01-23 | 浙江思博恩新材料科技有限公司 | 金刚线磨粒密度及分布的检测方法及装置 |
| CN103926276A (zh) * | 2014-03-25 | 2014-07-16 | 天津大学 | 一种在线油液磨粒监测装置及检测方法 |
| CN106941605A (zh) * | 2017-04-27 | 2017-07-11 | 华南理工大学 | 一种砂轮修锐修整的图像视觉监控装置及方法 |
| CN107657604A (zh) * | 2017-09-06 | 2018-02-02 | 西安交通大学 | 一种基于近场非标准光源的磨痕三维形貌原位获取方法 |
| CN107677216A (zh) * | 2017-09-06 | 2018-02-09 | 西安交通大学 | 一种基于光度立体视觉的多个磨粒三维形貌同步获取方法 |
| CN207039770U (zh) * | 2017-04-27 | 2018-02-23 | 华南理工大学 | 一种砂轮修锐修整的图像视觉监控装置 |
| KR20180136127A (ko) * | 2017-06-14 | 2018-12-24 | 엠.씨.케이 (주) | 연마제품 및 그 제조방법 |
-
2019
- 2019-03-12 CN CN201910185038.3A patent/CN109801291B/zh active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070084133A1 (en) * | 2005-10-18 | 2007-04-19 | 3M Innovative Properties Company | Agglomerate abrasive grains and methods of making the same |
| CN102890039A (zh) * | 2011-07-22 | 2013-01-23 | 浙江思博恩新材料科技有限公司 | 金刚线磨粒密度及分布的检测方法及装置 |
| CN103926276A (zh) * | 2014-03-25 | 2014-07-16 | 天津大学 | 一种在线油液磨粒监测装置及检测方法 |
| CN106941605A (zh) * | 2017-04-27 | 2017-07-11 | 华南理工大学 | 一种砂轮修锐修整的图像视觉监控装置及方法 |
| CN207039770U (zh) * | 2017-04-27 | 2018-02-23 | 华南理工大学 | 一种砂轮修锐修整的图像视觉监控装置 |
| KR20180136127A (ko) * | 2017-06-14 | 2018-12-24 | 엠.씨.케이 (주) | 연마제품 및 그 제조방법 |
| CN107657604A (zh) * | 2017-09-06 | 2018-02-02 | 西安交通大学 | 一种基于近场非标准光源的磨痕三维形貌原位获取方法 |
| CN107677216A (zh) * | 2017-09-06 | 2018-02-09 | 西安交通大学 | 一种基于光度立体视觉的多个磨粒三维形貌同步获取方法 |
Non-Patent Citations (2)
| Title |
|---|
| TONGHAI WU ET AL: "Morphological feature extraction based on multi-view images for wear debris analysis in on-line oil monitoring", 《TRIBOLOGY TRANSACTIONS》 * |
| 袁成清等: "磨粒的三维表面特征描述", 《摩擦学学报》 * |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115326656A (zh) * | 2022-10-14 | 2022-11-11 | 西南交通大学 | 交通土建用颗粒材料松铺层颗粒粒径及级配无损测量方法 |
| CN115326656B (zh) * | 2022-10-14 | 2022-12-23 | 西南交通大学 | 交通土建用颗粒材料松铺层颗粒粒径及级配无损测量方法 |
| CN115578523A (zh) * | 2022-11-21 | 2023-01-06 | 慧医谷中医药科技(天津)股份有限公司 | 一种多角度图像融合的舌部三维建模方法及系统 |
| CN115578523B (zh) * | 2022-11-21 | 2023-03-10 | 慧医谷中医药科技(天津)股份有限公司 | 一种多角度图像融合的舌部三维建模方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109801291B (zh) | 2021-05-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| WO2018028103A1 (zh) | 一种基于人眼视觉特性的电力线路无人机巡检方法 | |
| Kong et al. | Detecting abandoned objects with a moving camera | |
| CN111899172A (zh) | 一种面向遥感应用场景的车辆目标检测方法 | |
| Jeon et al. | Reconstruction-based pairwise depth dataset for depth image enhancement using CNN | |
| Wu et al. | Densely pyramidal residual network for UAV-based railway images dehazing | |
| CN112288627B (zh) | 一种面向识别的低分辨率人脸图像超分辨率方法 | |
| CN111539888B (zh) | 一种基于金字塔通道特征注意力的神经网络图像去雾方法 | |
| Gkitsas et al. | Panodr: Spherical panorama diminished reality for indoor scenes | |
| CN109801291A (zh) | 一种运动磨粒多表面三维形貌的获取方法 | |
| CN110222661A (zh) | 一种用于运动目标识别及跟踪的特征提取方法 | |
| CN115908354A (zh) | 一种基于双尺度策略和改进的yolov5网络的光伏面板缺陷检测方法 | |
| Maesumi et al. | Learning transferable 3D adversarial cloaks for deep trained detectors | |
| CN110111368B (zh) | 一种基于人体姿态识别的相似移动目标的检测跟踪方法 | |
| Atienza | A conditional generative adversarial network for rendering point clouds | |
| Yagi et al. | Diminished reality for privacy protection by hiding pedestrians in motion image sequences using structure from motion | |
| Li et al. | DAR-Net: dense attentional residual network for vehicle detection in aerial images | |
| CN113724273B (zh) | 一种基于神经网络区域目标分割的边缘光影融合方法 | |
| Jindal et al. | An ensemble mosaicing and ridgelet based fusion technique for underwater panoramic image reconstruction and its refinement | |
| CN112257505B (zh) | 一种应用于红外全景系统的虚假目标识别方法及存储介质 | |
| Jibrin et al. | Development of hybrid automatic segmentation technique of a single leaf from overlapping leaves image | |
| Wang et al. | PointGait: Boosting End-to-End 3D Gait Recognition with Point Clouds via Spatiotemporal Modeling | |
| Doshi et al. | Monochromatic Image Dehazing Using Enhanced Feature Extraction Techniques in Deep Learning | |
| CN105335399B (zh) | 一种信息处理方法及电子设备 | |
| Hu et al. | Dense 3D model reconstruction for digital city using computationally efficient multi-view stereo networks | |
| Zheng et al. | Piston Surface Defect Recognition Method Based on Image Processing |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20221208 Address after: 710061 Fengxi New Town, Xixian New District, Xi'an, Shaanxi Patentee after: Shaanxi Zhipu Weichuang Information Technology Co.,Ltd. Address before: Beilin District Xianning West Road 710049, Shaanxi city of Xi'an province No. 28 Patentee before: XI'AN JIAOTONG University |
|
| TR01 | Transfer of patent right |