CN109408909B - 一种三维粗糙表面微凸体拟合方法 - Google Patents

一种三维粗糙表面微凸体拟合方法 Download PDF

Info

Publication number
CN109408909B
CN109408909B CN201811155174.XA CN201811155174A CN109408909B CN 109408909 B CN109408909 B CN 109408909B CN 201811155174 A CN201811155174 A CN 201811155174A CN 109408909 B CN109408909 B CN 109408909B
Authority
CN
China
Prior art keywords
fitting
discrete points
height
micro
convex body
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.)
Active
Application number
CN201811155174.XA
Other languages
English (en)
Other versions
CN109408909A (zh
Inventor
唐进元
温昱钦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN201811155174.XA priority Critical patent/CN109408909B/zh
Publication of CN109408909A publication Critical patent/CN109408909A/zh
Application granted granted Critical
Publication of CN109408909B publication Critical patent/CN109408909B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/003D [Three Dimensional] image rendering
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30108Industrial image inspection
    • G06T2207/30164Workpiece; Machine component

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Graphics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Length Measuring Devices By Optical Means (AREA)

Abstract

本发明公开了一种三维粗糙表面微凸体拟合方法,包括以下步骤:a.测量待拟合粗糙表面的表面形貌离散点高度矩阵z(i,j);b.计算得到表面平均高度面h0,取其为拟合参考面;c.对表面形貌离散点进行分割,划分为若干封闭区域;d.根据拟合参考面以及表面形貌离散点高度之间的相互关系,对每一个封闭区域内的离散点进行判断,将符合微凸体特征的离散点序列判定为微凸体;e.以步骤d中判定得到的微凸体,根据最小均方误差,用椭球面对每个微凸体的离散点进行拟合,得到拟合后椭球方程,据此计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等表面微凸体参数。本发明的一种三维粗糙表面微凸体拟合方法,具有测量精确的优点。

Description

一种三维粗糙表面微凸体拟合方法
技术领域
本发明涉及一种三维表面拟合方法,尤其涉及一种三维粗糙表面微凸体拟合方法。
背景技术
机械加工零件表面在微观尺度上都是凹凸不平的粗糙表面,当两表面接触的时候,实际产生接触的是一系列离散分布的接触区域,接触面积仅是光滑名义接触面积下的一部分,远小于名义接触面积,这会导致真实接触压力远大于名义接触压力。因此,接触表面的真实接触面积、接触压力等直接影响着传动部件的承载能力以及摩擦、磨损、疲劳等性能。
基于赫兹接触理论与数理统计原理,考虑微凸体间接触的粗糙表面接触模型目前广泛运用于粗糙表面接触分析,粗糙表面的表面微凸体参数:微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等对接触分析计算结果影响极大,然而考虑到光学仪器测量的表面形貌受到采样间距的影响,其微凸体参数计算同样将受到采样间距的巨大影响,通过经典三点法、五点法等计算微凸体参数得到的微凸体平均曲率半径通常过小,不能合适的用于粗糙表面微凸体接触模型计算,这极大影响了接触分析结果的准确性。为了提高参数计算的准确性,减小采样间距的影响,有研究者采用抛物线对表面微凸体进行拟合简化,再计算表面微凸体参数,在一定程度上提高了计算准确性。
然而现有拟合计算方法局限于二维粗糙轮廓的拟合简化,而粗糙表面接触模型是基于三维表面的接触发展而来,二维拟合简化得到的微凸体与实际三维微凸体存在显著的差别,简单的二维轮廓无法完全准确的表征粗糙表面的三维形貌,因此,现有方法的准确性还有待提高。
发明内容
本发明所要解决的技术问题是提供一种拟合误差小并且精确度高的三维粗糙表面微凸体拟合方法。
本发明是通过以下的技术方案实现的:
一种三维粗糙表面微凸体拟合方法,包括以下步骤,
a.测量待拟合粗糙表面的表面形貌离散点高度矩阵z(i,j);
b.以步骤a中测量得到的表面形貌离散点高度序列z(i,j)为基础,计算得到表面平均高度面H0,取其为拟合参考面;
c.基于分水岭算法,在Matlab软件中对表面形貌离散点进行分割,划分为若干封闭区域;
d.根据拟合参考面以及表面形貌离散点高度之间的相互关系,对每一个封闭区域内的离散点进行判断,将符合微凸体特征的离散点序列判定为微凸体;
e.以步骤d中判定得到的微凸体,根据最小均方误差,用椭球面对每个微凸体的离散点进行拟合,得到拟合后椭球方程,据此计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等表面微凸体参数。
优选的,a步骤中采用白光干涉仪进行测量。
优选的,b步骤中计算表面平均高度面H0包括以下步骤:
(b1)根据表面形貌高度离散点矩阵z(i,j)计算得到平均高度线h0,计算公式为
Figure GDA0003911420240000031
式中,N为表面形貌离散点垂直方向的个数,M为表面形貌离散点水平方向的个数;
(b2)令拟合参考面为高度为h=h0的平面。
优选的,c步骤中,包括以下步骤:
(c1)计算表面十点高度Sz,Sz定义为表面形貌离散点中,5个最高顶点的高度和5个最深谷点的深度的平均值,计算公式为
Figure GDA0003911420240000032
式中,zmaxk和zmink(k=1,2,...,5)为5个最高顶点的高度和5个最深谷点的深度,若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最大值,则该点为顶点,同理若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最小值,则该点为谷点;
(c2)根据拟合参考面对表面形貌离散点矩阵z(i,j)进行处理,仅对高于拟合参考面的表面进行分割处理,若z(i,j)>h,则令z(i,j)=z(i,j)-h,若z(i,j)<=h,则令z(i,j)=h;
(c3)使用Matlab对表面形貌离散点矩阵z(i,j)进行垂直方向和水平方向的Sobel算子滤波,然后求取模值,得到表面形貌离散点的梯度幅值矩阵gradz(i,j),计算公式为:
Figure GDA0003911420240000033
式中,Ix(i,j)为使用Sobel算子对表面形貌离散点矩阵进行水平方向边缘滤波的卷积结果,Iy(i,j)为使用Sobel算子对表面形貌离散点矩阵进行垂直方向边缘滤波的卷积结果;
(c4)根据表面十点高度Sz对梯度幅值矩阵gradz(i,j)进行处理,若gradz(i,j)<0.05Sz,则令gradz(i,j)=0;
(c5)使用Matlab中的watershed分水岭算法函数对梯度幅值矩阵进行分水岭分割,将表面分割为若干封闭区域。
优选的,d步骤中,包括以下步骤:
(d1)使用Matlab中的regiongrops图像处理函数对分水岭算法分割后的结果进行统计处理,得到每个封闭区域内所有离散点的坐标和区域边界上离散点的坐标;
(d2)针对单个封闭区域,若该区域内所有离散点高度的平均值大于该区域边界上离散点高度的平均值,且该区域内所有离散点高度的平均值大于拟合参考面高度h,则将该封闭区域内的离散点判定为组成一个微凸体。
优选的,e步骤中,包括以下步骤:
(e1)针对判定后单个微凸体的封闭区域,根据Matlab中的regiongrops图像处理函数可确定与该封闭区域具有相同标准二阶中心距的椭圆,令a为与该区域具有相同标准二阶中心矩的椭圆的长轴长度,b为与该区域具有相同标准二阶中心矩的椭圆的短轴长度,O点为该区域的质心;
(e2)针对(e1)中确定的椭圆区域,以O点为坐标系原点,椭圆长轴方向为x轴,椭圆短轴方向为y轴,高度方向为z轴建立坐标系,设区域内共有n个离散点,离散点高度分别为z(x1,y1),…,z(xi,yi),…,z(xn,yn);
(e3)用椭球面对该椭圆区域内离散点进行拟合,令椭球面方程为:
Figure GDA0003911420240000041
式中,c为椭球面系数;
(e4)运用最小均方根误差来判定最适合的拟合椭球面,拟合椭球面最小均方误差Er由以下公式确定:
Figure GDA0003911420240000051
Figure GDA0003911420240000052
即可计算得到每个椭球面的系数c:
Figure GDA0003911420240000053
将得到的椭球面系数c代入椭球面方程,即可得到微凸体离散点的拟合椭球面方程,确定拟合形貌。
优选的,e步骤中,包括以下步骤:
(e5)每一个椭球顶端曲率半径即为微凸体顶峰曲率半径,椭球在x轴方向等效曲率半径为Rx,在y轴方向等效曲率半径为Ry,计算公式如下:
Figure GDA0003911420240000054
Figure GDA0003911420240000055
设共有m个微凸体,则微凸体平均曲率半径R为:
Figure GDA0003911420240000056
(e6)微凸体峰点高度标准偏差σ、微凸体密度η为:
Figure GDA0003911420240000057
η=m/A (11)
式中,A为整个粗糙表面的名义面积。
有益效果是:
与现有技术相比,本发明的一种三维粗糙表面微凸体拟合方法通过实测工件的表面形貌,提取表面三维形貌离散点,基于分水岭算法对表面区域进行分割,通过离散点高度与初始参考面对分割区域进行判定,进一步确定拟合微凸体的区域,根据最小均方误差对每个微凸体进行椭球面拟合,方便计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等微凸体参数,拟合得到的椭球微凸体更加贴近实际接触的微凸体,使得拟合误差小,能更加准确的表征三维接触计算中的表面形貌,在一定程度上避免了采样间距对微凸体参数计算的巨大影响,使微凸体参数计算更加合理,以此计算得到微凸体参数能更加准确的表征接触微凸体的几何特征,提高了粗糙表面接触分析的准确性。
附图说明
以下结合附图对本发明的具体实施方式作进一步的详细说明,其中:
图1为工件表面形貌示意图;
图2为分割结果局部图;
图3为单个微凸体拟合前原始形貌示意图;
图4为单个微凸体拟合后形貌示意图;
图5为工件表面形貌拟合结果图。
具体实施方式
一种三维粗糙表面微凸体拟合方法,包括以下步骤:
a.测量待拟合粗糙表面的表面形貌离散点高度矩阵z(i,j);
b.以步骤a中测量得到的表面形貌离散点高度序列z(i,j)为基础,计算得到表面平均高度面H0,取其为拟合参考面;
c.基于分水岭算法,在Matlab软件中对表面形貌离散点进行分割,划分为若干封闭区域;
d.根据拟合参考面以及表面形貌离散点高度之间的相互关系,对每一个封闭区域内的离散点进行判断,将符合微凸体特征的离散点序列判定为微凸体;
e.以步骤d中判定得到的微凸体,根据最小均方误差,用椭球面对每个微凸体的离散点进行拟合,得到拟合后椭球方程,据此计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差等表面微凸体参数。
较佳的实施方式,a步骤中采用白光干涉仪进行测量。
较佳的实施方式,b步骤中计算表面平均高度面H0包括以下步骤:
(b1)根据表面形貌高度离散点矩阵z(i,j)计算得到平均高度线h0,计算公式为
Figure GDA0003911420240000071
式中,N为表面形貌离散点垂直方向的个数,垂直方向也即是x方向,M为表面形貌离散点水平方向的个数,水平方向也即是y方向,表面形貌图如图1所示;
(b2)令拟合参考面为高度为h=h0的平面。
较佳的实施方式,c步骤中,包括以下步骤:
(c1)计算表面十点高度Sz,Sz定义为表面形貌离散点中,5个最高顶点的高度和5个最深谷点的深度的平均值,计算公式为
Figure GDA0003911420240000072
Figure GDA0003911420240000081
式中,zmaxk和zmink(k=1,2,...,5)为5个最高顶点的高度和5个最深谷点的深度,若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最大值,则该点为顶点,同理若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最小值,则该点为谷点;
(c2)根据拟合参考面对表面形貌离散点矩阵z(i,j)进行处理,仅对高于拟合参考面的表面进行分割处理,若z(i,j)>h,则令z(i,j)=z(i,j)-h,若z(i,j)<=h,则令z(i,j)=h;
(c3)使用Matlab对表面形貌离散点矩阵z(i,j)进行垂直方向和水平方向的Sobel算子滤波,然后求取模值,得到表面形貌离散点的梯度幅值矩阵gradz(i,j),计算公式为:
Figure GDA0003911420240000082
式中,Ix(i,j)为使用Sobel算子对表面形貌离散点矩阵进行水平方向边缘滤波的卷积结果,Iy(i,j)为使用Sobel算子对表面形貌离散点矩阵进行垂直方向边缘滤波的卷积结果;
(c4)根据表面十点高度Sz对梯度幅值矩阵gradz(i,j)进行处理,若gradz(i,j)<0.05Sz,则令gradz(i,j)=0;
(c5)使用Matlab中的watershed分水岭算法函数对梯度幅值矩阵进行分水岭分割,将表面分割为若干封闭区域,表面分割结果局部如图2所示。
较佳的实施方式,d步骤中,包括以下步骤:
(d1)使用Matlab中的regiongrops图像处理函数对分水岭算法分割后的结果进行统计处理,得到每个封闭区域内所有离散点的坐标和区域边界上离散点的坐标;
(d2)针对单个封闭区域,若该区域内所有离散点高度的平均值大于该区域边界上离散点高度的平均值,且该区域内所有离散点高度的平均值大于拟合参考面高度h,则将该封闭区域内的离散点判定为组成一个微凸体。
较佳的实施方式,e步骤中,包括以下步骤:
(e1)针对判定后单个微凸体的封闭区域,根据Matlab中的regiongrops图像处理函数可确定与该封闭区域具有相同标准二阶中心距的椭圆,令a为与该区域具有相同标准二阶中心矩的椭圆的长轴长度,b为与该区域具有相同标准二阶中心矩的椭圆的短轴长度,O点为该区域的质心;
(e2)针对(e1)中确定的椭圆区域,以O点为坐标系原点,椭圆长轴方向为x轴,椭圆短轴方向为y轴,高度方向为z轴建立坐标系,设区域内共有n个离散点,离散点高度分别为z(x1,y1),…,z(xi,yi),…,z(xn,yn);
(e3)用椭球面对该椭圆区域内离散点进行拟合,令椭球面方程为:
Figure GDA0003911420240000091
式中,c为椭球面系数;
(e4)运用最小均方根误差来判定最适合的拟合椭球面,拟合椭球面最小均方误差Er由以下公式确定:
Figure GDA0003911420240000092
Figure GDA0003911420240000093
即可计算得到每个椭球面的系数c:
Figure GDA0003911420240000094
将得到的椭球面系数c代入椭球面方程,即可得到微凸体离散点的拟合椭球面方程,确定拟合形貌,表面拟合结果如图5所示。
较佳的实施方式,e步骤中,包括以下步骤:
(e5)每一个椭球顶端曲率半径即为微凸体顶峰曲率半径,椭球在x轴方向等效曲率半径为Rx,在y轴方向等效曲率半径为Ry,计算公式如下:
Figure GDA0003911420240000101
Figure GDA0003911420240000102
设共有m个微凸体,则微凸体平均曲率半径R为:
Figure GDA0003911420240000103
(e6)微凸体峰点高度标准偏差σ、微凸体密度η为:
Figure GDA0003911420240000104
η=m/A (11)式中,A为整个粗糙表面的名义面积。
以下,根据上述的方法提供具体算例:
用白光干涉仪Wyko NT9100对平板工件表面形貌进行测量,选择5×镜头倍数,采样间隔1μm,得到表面形貌离散点高度矩阵z(i,j),共472×632个离散点;
一、确定拟合参考面
根据表面形貌高度离散点矩阵z(i,j)计算表面平均高度h0,数据代入式(1):
Figure GDA0003911420240000105
二、分割表面区域
若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最大值,则该点为顶点,同理若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最小值,则该点为谷点,以此得到所有顶点与谷点深度,计算表面十点高度Sz,Sz为5个最高顶点的高度和5个最深谷点的深度的平均值,数据代入式(2):
Figure GDA0003911420240000111
再根据拟合参考面对表面形貌离散点矩阵z(i,j)进行处理,若z(i,j)>0.0075,则令z(i,j)=z(i,j)-0.0075,若z(i,j)<=0.0075,则令z(i,j)=0.0075;
然后,使用Matlab对表面形貌离散点矩阵z(i,j)进行垂直方向和水平方向的Sobel算子滤波,然后求取模值,得到表面形貌离散点的梯度幅值矩阵gradz(i,j):
Figure GDA0003911420240000112
式中,Ix(i,j)为使用Sobel算子对表面形貌离散点矩阵进行水平方向边缘滤波的卷积结果,Iy(i,j)为使用Sobel算子对表面形貌离散点矩阵进行垂直方向边缘滤波的卷积结果;
再根据表面十点高度Sz对梯度幅值矩阵进行处理,若gradz(i,j)<0.05Sz,即如果gradz(i,j)<0.275,则令gradz(i,j)=0;
最后使用Matlab中的watershed分水岭算法函数对梯度幅值矩阵进行分水岭分割,将表面分割为若干封闭区域,分割结果局部如图2所示。
三、微凸体判定
使用Matlab中的regiongrops图像处理函数对分水岭算法分割后的结果进行统计处理,得到每个封闭区域内所有离散点的坐标和区域边界上离散点的坐标。
然后针对单个封闭区域,若该区域内所有离散点高度的平均值大于该区域边界上离散点高度的平均值,且该区域内所有离散点高度的平均值大于拟合参考面高度h,则将该封闭区域内的离散点判定为组成一个微凸体。
四、微凸体拟合
首先,针对判定后单个微凸体的封闭区域,根据Matlab中的regiongrops图像处理函数可确定与该封闭区域具有相同标准二阶中心距的椭圆,令a为与该区域具有相同标准二阶中心矩的椭圆的长轴长度,b为与该区域具有相同标准二阶中心矩的椭圆的短轴长度,O点为该区域的质心。任取一微凸体为例,如一区域内表面微凸体形貌如图3所示,该区域质心为点O(279,37),与该区域具有相同标准二阶中心矩的椭圆的长轴长度a=8.7μm,与该区域具有相同标准二阶中心矩的椭圆的短轴长度b=6.4μm;
然后,针对上段中确定的椭圆区域,以点O(279,37)为坐标系原点,椭圆长轴方向为x轴,椭圆短轴方向为y轴,高度方向为z轴建立坐标系,则区域内共有46个离散点,离散点高度分别为z(x1,y1),…,z(xi,yi),…,z(x46,y46);
再用椭球面对该椭圆区域内离散点进行拟合,数据代入式(4),令椭球面方程为:
Figure GDA0003911420240000121
式中,c为椭球面系数;
然后,运用最小均方根误差来判定最适合的拟合椭球面,数据代入式(5),拟合椭球面最小均方误差Er由以下公式确定:
Figure GDA0003911420240000122
Figure GDA0003911420240000123
即可计算得到每个椭球面的系数c,数据代入式(6),可得:
Figure GDA0003911420240000124
Figure GDA0003911420240000131
将得到的椭球面系数c代入椭球面方程式(15),即可得到微凸体离散点的拟合椭球面方程,确定拟合形貌,单个微凸体拟合后形貌如图4所示,整体表面拟合结果如图5所示。
五、微凸体参数的计算
每一个椭球顶端曲率半径即为微凸体顶峰曲率半径,椭球在x轴方向等效曲率半径为Rx,在y轴方向等效曲率半径为Ry,数据分别代入式(7)、(8),如上述微凸体拟合后顶峰曲率半径如下:
Figure GDA0003911420240000132
Figure GDA0003911420240000133
分割并判定后共有859个微凸体,数据代入式(9),则微凸体平均曲率半径R为:
Figure GDA0003911420240000134
将数据分别代入式(10)、(11),即可获得微凸体峰点高度标准偏差σ、微凸体密度η为:
Figure GDA0003911420240000135
η=859/(472×632)=0.0029/μm2 (22)。
以上实施例仅用以说明本发明的技术方案而并非对其进行限制,凡未脱离本发明精神和范围的任何修改或者等同替换,其均应涵盖在本发明技术方案的范围内。

Claims (7)

1.一种三维粗糙表面微凸体拟合方法,其特征在于:包括以下步骤,
a.测量待拟合粗糙表面的表面形貌离散点高度矩阵z(i,j);
b.以步骤a中测量得到的表面形貌离散点高度序列z(i,j)为基础,计算得到表面平均高度面H0,取其为拟合参考面;
c.基于分水岭算法,在Matlab软件中对表面形貌离散点进行分割,划分为若干封闭区域;
d.根据拟合参考面以及表面形貌离散点高度之间的相互关系,对每一个封闭区域内的离散点进行判断,将符合微凸体特征的离散点序列判定为微凸体;
e.以步骤d中判定得到的微凸体,根据最小均方误差,用椭球面对每个微凸体的离散点进行拟合,得到拟合后椭球方程,据此计算微凸体平均曲率半径、微凸体密度、微凸体峰点高度标准偏差。
2.根据权利要求1所述的三维粗糙表面微凸体拟合方法,其特征在于:a步骤中采用白光干涉仪进行测量。
3.根据权利要求1所述的三维粗糙表面微凸体拟合方法,其特征在于,b步骤中计算表面平均高度面H0包括以下步骤:
(b1)根据表面形貌高度离散点矩阵z(i,j)计算得到平均高度线h0,计算公式为:
Figure FDA0003911420230000011
式中,N为表面形貌离散点垂直方向的个数,M为表面形貌离散点水平方向的个数;
(b2)令拟合参考面为高度为h=h0的平面。
4.根据权利要求1所述的三维粗糙表面微凸体拟合方法,其特征在于,c步骤中,包括以下步骤:
(c1)计算表面十点高度Sz,Sz定义为表面形貌离散点中,5个最高顶点的高度和5个最深谷点的深度的平均值,计算公式为
Figure FDA0003911420230000021
式中,zmaxk和zmink(k=1,2,...,5)为5个最高顶点的高度和5个最深谷点的深度,若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最大值,则该点为顶点,同理若z(i,j)在矩阵z(i-1:i+1,j-1:j+1)中为最小值,则该点为谷点;
(c2)根据拟合参考面对表面形貌离散点矩阵z(i,j)进行处理,仅对高于拟合参考面的表面进行分割处理,若z(i,j)>h,则令z(i,j)=z(i,j)-h,若z(i,j)<=h,则令z(i,j)=h;
(c3)使用Matlab对表面形貌离散点矩阵z(i,j)进行垂直方向和水平方向的Sobel算子滤波,然后求取模值,得到表面形貌离散点的梯度幅值矩阵gradz(i,j),计算公式为:
Figure FDA0003911420230000022
式中,Ix(i,j)为使用Sobel算子对表面形貌离散点矩阵进行水平方向边缘滤波的卷积结果,Iy(i,j)为使用Sobel算子对表面形貌离散点矩阵进行垂直方向边缘滤波的卷积结果;
(c4)根据表面十点高度Sz对梯度幅值矩阵gradz(i,j)进行处理,若gradz(i,j)<0.05Sz,则令gradz(i,j)=0;
(c5)使用Matlab中的watershed分水岭算法函数对梯度幅值矩阵进行分水岭分割,将表面分割为若干封闭区域。
5.根据权利要求1所述的三维粗糙表面微凸体拟合方法,其特征在于,d步骤中,包括以下步骤:
(d1)使用Matlab中的regiongrops图像处理函数对分水岭算法分割后的结果进行统计处理,得到每个封闭区域内所有离散点的坐标和区域边界上离散点的坐标;
(d2)针对单个封闭区域,若该区域内所有离散点高度的平均值大于该区域边界上离散点高度的平均值,且该区域内所有离散点高度的平均值大于拟合参考面高度h,则将该封闭区域内的离散点判定为组成一个微凸体。
6.根据权利要求1所述的三维粗糙表面微凸体拟合方法,其特征在于,e步骤中,包括以下步骤:
(e1)针对判定后单个微凸体的封闭区域,根据Matlab中的regiongrops图像处理函数可确定与该封闭区域具有相同标准二阶中心距的椭圆,令a为与该区域具有相同标准二阶中心矩的椭圆的长轴长度,b为与该区域具有相同标准二阶中心矩的椭圆的短轴长度,O点为该区域的质心;
(e2)针对(e1)中确定的椭圆区域,以O点为坐标系原点,椭圆长轴方向为x轴,椭圆短轴方向为y轴,高度方向为z轴建立坐标系,设区域内共有n个离散点,离散点高度分别为z(x1,y1),…,z(xi,yi),…,z(xn,yn);
(e3)用椭球面对该椭圆区域内离散点进行拟合,令椭球面方程为:
Figure FDA0003911420230000031
式中,c为椭球面系数;
(e4)运用最小均方根误差来判定最适合的拟合椭球面,拟合椭球面最小均方误差Er由以下公式确定:
Figure FDA0003911420230000041
Figure FDA0003911420230000042
即可计算得到每个椭球面的系数c:
Figure FDA0003911420230000043
将得到的椭球面系数c代入椭球面方程,即可得到微凸体离散点的拟合椭球面方程,确定拟合形貌。
7.根据权利要求6所述的三维粗糙表面微凸体拟合方法,其特征在于:e步骤中,包括以下步骤:
(e5)每一个椭球顶端曲率半径即为微凸体顶峰曲率半径,椭球在x轴方向等效曲率半径为Rx,在y轴方向等效曲率半径为Ry,计算公式如下:
Figure FDA0003911420230000044
Figure FDA0003911420230000045
设共有m个微凸体,则微凸体平均曲率半径R为:
Figure FDA0003911420230000046
(e6)微凸体峰点高度标准偏差σ、微凸体密度η为:
Figure FDA0003911420230000047
η=m/A (11)
式中,A为整个粗糙表面的名义面积。
CN201811155174.XA 2018-09-30 2018-09-30 一种三维粗糙表面微凸体拟合方法 Active CN109408909B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811155174.XA CN109408909B (zh) 2018-09-30 2018-09-30 一种三维粗糙表面微凸体拟合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811155174.XA CN109408909B (zh) 2018-09-30 2018-09-30 一种三维粗糙表面微凸体拟合方法

Publications (2)

Publication Number Publication Date
CN109408909A CN109408909A (zh) 2019-03-01
CN109408909B true CN109408909B (zh) 2022-12-13

Family

ID=65465833

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811155174.XA Active CN109408909B (zh) 2018-09-30 2018-09-30 一种三维粗糙表面微凸体拟合方法

Country Status (1)

Country Link
CN (1) CN109408909B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111862300B (zh) * 2020-06-16 2022-07-01 湖南大学 一种三维粗糙表面抛物面形微凸体拟合方法
CN114354541B (zh) * 2020-10-14 2024-07-16 中粮集团有限公司 大米加工精度检测方法及系统
CN112967256B (zh) * 2021-03-09 2023-11-24 扬州大学 一种基于空间分布的隧道椭圆化检测方法
CN113221892A (zh) * 2021-05-12 2021-08-06 佛山育脉科技有限公司 手掌图像确定方法、装置及计算机可读存储介质
CN117906529B (zh) * 2024-03-18 2024-05-28 板石智能科技(深圳)有限公司 倾斜空间平面自动平衡方法、装置、电子设备及存储介质

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2897455A1 (fr) * 2006-02-13 2007-08-17 Univ Hokkaido Nat Univ Corp Dispositif, procede et programme de segmentation de donnees de modele en treillis
CN106709207B (zh) * 2017-01-16 2020-06-16 东北大学 一种考虑粗糙表面微凸体相互作用影响的确定受载结合部法向接触刚度的方法
CN107423497B (zh) * 2017-07-12 2018-07-24 中南大学 一种粗糙表面微凸体拟合方法和系统

Also Published As

Publication number Publication date
CN109408909A (zh) 2019-03-01

Similar Documents

Publication Publication Date Title
CN109408909B (zh) 一种三维粗糙表面微凸体拟合方法
CN110807781B (zh) 一种保留细节与边界特征的点云精简方法
CN108986048B (zh) 基于线激光扫描三维点云快速复合滤波处理方法
CN102944174B (zh) 一种三维激光点云数据的去噪与精简方法及系统
CN108846888B (zh) 一种古木建筑构件精细尺寸信息自动化提取方法
CN109506580A (zh) 基于线激光三维扫描的锪孔质量检测方法
CN109147040B (zh) 基于模板的人体点云孔洞修补方法
KR20150058522A (ko) 타이어 트레드 파라미터를 분석하는 시스템 및 방법
KR20140020837A (ko) 표면 검사에 이용하기 위한 타이어 표면의 삼차원 이미지 전처리 방법
CN110288640B (zh) 基于凸密度极值的点云配准方法
KR20140009209A (ko) 연속된 b-스플라인 변형을 이용한 타이어 표면의 삼차원 이미지 전처리 방법
CN112785596B (zh) 基于dbscan聚类的点云图螺栓分割和高度测量方法
CN113052797B (zh) 基于深度图像处理的bga锡球三维检测方法
CN115984359A (zh) 基于球坐标积分的地基激光点云单木树冠体积提取方法
CN117057206A (zh) 一种三维模具的智能建模方法及系统
CN103218809A (zh) 一种珍珠长度参数的图像测量方法
CN114419140A (zh) 一种轨道激光测量装置光斑中心的定位算法
CN110310322B (zh) 一种10微米级高精度器件装配表面检测方法
CN107908841A (zh) 三维壁面可抓取位置判别算法
CN117928385A (zh) 一种基于远程无人机和传感器的工程施工智能测量方法
CN117152344B (zh) 基于照片重建点云的隧道围岩结构面解析方法及系统
CN117292181A (zh) 基于3d点云处理的钣金件孔组分类及全尺寸测量方法
CN112541897A (zh) 一种基于人工智能的面瘫程度评估系统
CN110207618B (zh) 三维扫描测量数据的表面线数据提取方法
CN107609255A (zh) 车身特征线消失位置的检测方法

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