CN116718137A - 面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法 - Google Patents
面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法 Download PDFInfo
- Publication number
- CN116718137A CN116718137A CN202310608783.0A CN202310608783A CN116718137A CN 116718137 A CN116718137 A CN 116718137A CN 202310608783 A CN202310608783 A CN 202310608783A CN 116718137 A CN116718137 A CN 116718137A
- Authority
- CN
- China
- Prior art keywords
- point
- distance
- function
- measuring
- point cloud
- 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.)
- Pending
Links
- 230000002159 abnormal effect Effects 0.000 claims abstract description 65
- 238000006243 chemical reaction Methods 0.000 claims abstract description 38
- 239000011159 matrix material Substances 0.000 claims abstract description 36
- 238000005259 measurement Methods 0.000 claims abstract description 35
- 238000000034 method Methods 0.000 claims abstract description 20
- 238000011156 evaluation Methods 0.000 claims abstract description 10
- 230000009466 transformation Effects 0.000 claims abstract description 4
- 239000013598 vector Substances 0.000 claims description 53
- 238000004364 calculation method Methods 0.000 claims description 12
- 238000013519 translation Methods 0.000 claims description 9
- 238000012937 correction Methods 0.000 claims description 6
- 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 3
- 238000010276 construction Methods 0.000 claims description 2
- 230000008569 process Effects 0.000 abstract description 4
- 238000012545 processing Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 230000036961 partial effect Effects 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 238000003754 machining Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005498 polishing Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
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/30—Measuring arrangements characterised by the use of optical techniques for measuring roughness or irregularity of surfaces
-
- 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
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- 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
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0004—Industrial image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
- G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30108—Industrial image inspection
- G06T2207/30164—Workpiece; Machine component
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
- Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
- Y02P90/30—Computing systems specially adapted for manufacturing
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Mathematical Analysis (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Computing Systems (AREA)
- Algebra (AREA)
- General Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Multimedia (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Evolutionary Computation (AREA)
- Artificial Intelligence (AREA)
- Quality & Reliability (AREA)
- Operations Research (AREA)
- Length Measuring Devices With Unspecified Measuring Means (AREA)
Abstract
本发明公开了一种面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,首先构建并输入目标点云P以及源点云Q;对目标点云和源点云进行配对,并分为正负测点;计算点对的法向距离,输入函数变形参数a,设置权重尺度因子c;构建鲁棒函数,根据点对距离将测点分为正常测点与异常测点,并分别施加不同程度的权重值;建立目标函数,求解精配准的刚体转换矩阵,利用计算出来的刚体转换矩阵作用于源点云、以更新源点云位置;定义误差评价函数;重新计算点对然后按照上述过程迭代一次,并每隔K次迭代更新一次c的值,直到达到收敛条件,最终获得精配准的刚体转换矩阵T。本发明能够在存在大量异常点云情况下,进行精配准,效率高,精度好。
Description
技术领域
本发明属于复杂构件测量-加工一体化制造领域,具体地,涉及一种点云精配准算法,特别涉及一种面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法(RFWVM)。
背景技术
航空发动机和燃气轮机叶片通常具有前后缘薄且曲率变化较大、型面扭曲,材料难加工等特点。其型面加工质量对发动机性能有决定性的影响。使用机器人对其磨抛过程中,在叶片的高曲率区域,如前后缘等地方磨抛时,容易出现过磨(材料去除量过多)的现象,而在低曲率区域容易出现欠磨(材料去除量较少)的现象,叶片的加工需要多次磨抛,在加工过程中,需要校正误差,确定加工余量,而这就依赖测量点云与模型点云的精配准。在面对这种结构复杂、余量分布不均且存在大量异常点云的情况,现有的点云配准算法难以实现有效配准,获取工件加工过程中的真实余量情况困难。文献“A method forregistration of 3-D shapes”(IEEE Transactions on Pattern Analysis and MachineIntelligence,14(1992)239-256)提出一种迭代最近点(ICP)算法,通过迭代计算最小化点到点距离平方和完成点云精配准,该算法未考虑到存在大量余量点云的情况,在用于上述复杂工件时会导致匹配失真。申请号为CN201510226138.8的授权发明专利提出了一种基于距离方差最小的工件点云匹配算法,该算法可在点云部分缺失以及点云密度不均等情况下取得良好的匹配效果,但对于存在大量不均匀余量和结构偏差的工件,该算法也会匹配失真。申请号为CN202110573411.X的授权发明专利提出一种基于加权正负余量方差最小化算法的复杂工件光学测量方法,该算法能够抑制部分余量不均引起的匹配失真,但在结构偏差程度不大但偏差点数量较多的情况下,本该被抑制的结构偏差点云却被视为正常余量点云,参与配准计算,使算法陷入局部最优解,无法有效配准。
发明内容
本发明针对上述问题,提出了一种面向复杂构件视觉测量的鲁棒函数加权方差最小化算法,应用于复杂构件的定位与测量,可有效解决传统算法因大量异常点云存在导致的匹配失真问题。
本发明的目的在于提供一种面向复杂构件视觉测量的鲁棒函数加权方差最小化(Robust Function Weighted Variance Minimization,RFWVM)算法,建立了一种可有效抑制结构偏差点云与余量不均点云的鲁棒权重函数,在面对大量异常点云时不易出现的匹配失真,算法具有极强的鲁棒性,能够抑制大量余量不均点云对精配准的影响,尤其适用于复杂构件定位测量,并通过对复杂构件的测量定位,验证了该算法的有效性。
为实现上述目的,本发明的技术方案如下:
一种面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,包括如下步骤:
步骤1、构建并输入目标点云P以及源点云Q;
步骤2、利用KD-Tree算法搜索目标点云和源点云中最近邻点对集合,目标点云中的一个测点pi与源点云中与该测点pi距离最近的一个点qj组成一对最近邻点对,并基于点对的位置关系将目标点云中的点区分为正负测点;
步骤3、计算最近邻点对集合中每个点对中pi到qj的法向距离,输入函数变形参数α值,设置权重尺度因子c值;
步骤4、构建鲁棒函数,根据点对的距离将目标点云Q中的点区分为正常测点与异常点云,并通过鲁棒函数对点对距离施加不同程度的权重值;
步骤5、基于点对法向距离与鲁棒函数建立鲁棒权重加权方差最小化的目标函数;
步骤6、根据鲁棒权重加权方差最小化的目标函数求解源点云Q和目标点云P之间精配准的刚体转换矩阵,利用计算出来的刚体转换矩阵作用于源点云、以更新源点云位置;
步骤7、定义余量不均点云、存在结构偏差点云时的误差评价函数;
步骤8、重复步骤2-7,并每隔K次迭代更新一次c的值,直到达到收敛条件,最终获得精配准的刚体转换矩阵T。
进一步地,步骤1中,源点云Q构建方法如下:
使用CAD模型离散出源点云Q,并通过PCL库中的NormalEstimation函数去估算源点云的法向量,设置质心为视点,并计算使每个法向量的方向相对质心一致向外,完成法向量定向;
所述目标点云P使用结构光三维扫描仪扫描工件得到。
进一步地,步骤2中,区分正负测点的方法为:目标点云中的测点位置减去源点云中的最近邻点位置得到的向量与源点云相应点的法向量求点积,若其值大于0,则为正余量测点,若小于0则为负余量测点。
进一步地,步骤3中具体包括以下步骤:
步骤3.1、点对法向距离diT为点对中点到平面距离,包含正向测点对法向距离dkT和负向测点对法向距离dsT,i∈[1,m+n],i为顺序角标,diT表示第i个测点的点对法向距离,m为正向测点的总数,n为负向测点的总数;
定义δ∈R3、t∈R3分别为旋转向量跟平移向量,R3表示三维向量空间,建立点对法向距离函数,如公式(1):
其中,pi表示目标点云P中的第i个测点,pi+为测点pi经过单步转换后的测点,qj表示源点云Q中与pi组成最邻近点对的点,F为距离函数,nj表示第j个源点云qj处的法向量,niT=nj,是转换向量;
步骤3.2、根据点云结构特征输入函数变形参数α值,设置权重尺度因子c值。
进一步地,步骤4中,将异常点云分为可控异常测点和完全异常测点,并对完全异常测点的权重赋值0,根据正向测点和负向测点分别计算和判断如下:
对于正向测点,构建正向测点对距离鲁棒函数如公式(2):
其中wk为第k个正向测点的距离权重函数,当wk<0时修正为wk=0,dkT为第k个正向测点到平面距离,是经单步转换后的所有正向测点到平面距离均值,若为第一次迭代的转换向量ξ=06×1,m是正向测点的总数量,k∈[1,m],α∈R为控制权重函数鲁棒性的变形参数;c∈(0,∞]为控制权重函数零点值的权重尺度因子;d0+表示正向测点权重曲线的零点,;
对于负向测点,构建负向测点对距离鲁棒函数如公式(3):
其中ws为第s个负向测点的距离权重函数,当ws<0时修正为ws=0,dsT为第s个负向测点到平面距离,是经单步转换后的负向测点到平面距离均值,若为第一次迭代的转换向量ξ=06×1,n是负向测点的总数量,s∈[1,n],α∈R为控制权重函数鲁棒性的变形参数;d0-表示负向测点权重曲线的零点。
进一步地,步骤5中,建立鲁棒函数加权方差最小化算法的目标函数如公式(4):
其中,为第k个正向测点的正向鲁棒函数加权法向距离偏差,满足 为第s个负向测点的负向鲁棒函数加权法向距离偏差,满足为所有正向测点鲁棒函数加权法向距离均值,/>为所有负向测点鲁棒函数加权法向距离均值;dkT为正向测点对法向距离,dsT为负向测点对法向距离;
考虑异常点云,将公式(4)改写表示为公式(5):
其中,与/>为正负向在正常范围内的测点的权值,其值为1;m1、m2分别表示正常正向测点总数和异常正向测点总数;n1、n2分别表示正常负向测点总数和异常负向测点总数;/>为所有异常测点的加权距离偏差平方和,包括正向异常测点与负向异常测点两部分,Df为第f个异常测点的加权距离偏差;所有正向测点鲁棒函数加权法向距离均值/>的计算如公式(6)所示:
其中,m21和m20分别为可控异常测点和完全异常测点的总数量,为第f1+个可控异常测点的鲁棒函数权重,/>为第f0+个完全异常测点的鲁棒函数权重,/>为第f1+个可控异常测点对应的点对法向距离,/>为第f0+个完全异常测点对应的点对法向距离,由步骤5中所建立的鲁棒权重函数,得到/>将其与/>代入公式(6)中,得到如公式(7)所示:
进一步地,步骤6具体包括以下步骤:
步骤6.1、将正向测点的鲁棒函数加权法向距离偏差平方和化简,如公式(10):
上式中,E+是6×6矩阵,F+是6×1矩阵,D+为标量;
步骤6.2、将负向测点的鲁棒函数加权法向距离偏差平方和进行同样化简,如公式(11):
其中,E-是6×6矩阵,F-是6×1矩阵,D-为标量;
步骤6.3、根据公式(4)所示目标函数G(R,t),使用目标函数对转换向量ξ求导并求解,如公式(12):
上式中, Ak,As均为1×6的矩阵,/>所有正向测点的点对距离均值,/>是未经单步转换的第k个正向测点的点对距离;/>所有负向测点的点对距离均值,/>是未经单步转换的第s个负向测点的点对距离;/>鲁棒函数加权负向平均向量,/>鲁棒函数加权正向平均向量;
故旋转矩阵与平移矩阵求得如公式(13):
e为自然底数,t为平移向量,为旋转向量δ的反对称矩阵;
步骤6.4、将R与T作用与源点云,更新源点云位置。
进一步地,步骤7中,误差评价函数如公式(14):
公式(14)中,pi表示配准后的测点中点对距离小于平均点对距离的测点,qj表示从源点云Q中搜到的最近点。
进一步地,步骤8中,K取值范围为2-8。
进一步地,步骤8中,收敛条件为迭代次数M大于设定值Mmax或连续两次迭代误差评价函数值相差小于设定阈值G。
本发明的有益效果:
本发明充分考虑了大量复杂异常点云导致的匹配失真的问题,将算法目标函数中施加鲁棒权重函数,并在迭代过程中更改权重函数曲线,逐步弱化异常点云对配准的影响,能够提高配准精度,避免匹配倾斜现象的发生,相比较于传统配准算法,本发明在存在大量异常点云的情况下,能够快速准确的得到最优解,不易陷入局部最优,能够有效解决传统算法在此情况下的匹配失真的问题,适用于存在大量异常点云的复杂工件的精配准,具有强稳定性和强抗异常干扰性等特点。
附图说明
图1鲁棒函数加权方差最小化算法精配准流程图。
图2为本发明实施例中正负余量点的鲁棒函数权重随迭代次数曲线变化图,其中图2中(a)为正向测点权重变化曲线,图2中(b)为负向测点权重变化曲线。
图3叶片CAD点云及将25.8%的点偏置2mm之后的模拟测量叶片点云初始位置图。
图4RFWVM算法与WPMAVM算法全局配准误差对比图。
图5为不同算法叶片配准误差色谱图,其中图5(a)为ICP算法叶片全局配准误差色谱图,图5(b)为WPMAVM算法叶片全局配准误差色谱图,图5(c)为RFWVM算法叶片全局配准误差色谱图。
图6为本发明实施例中ICP、WPMAVM和RFWVM算法全局配准定位误差图。
图7车身扫描点云及车身CAD点云初始位置图。
图8为本发明实施例中不同算法对冲压车身工件全局配准误差色谱图,其中图8(a)为ICP算法车身工件全局配准误差色谱图,图8(b)为WPMAVM算法车身工件全局配准误差色谱图,图8(c)为RFWVM算法车身工件全局配准误差色谱图。
具体实施方式
下面结合附图和实施例对本发明的实施方式作进一步详细描述。以下实施例用于说明本发明,但不能用来限制本发明的范围。
如图1所示,本发明提供一种面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,包括如下步骤:
步骤1、构建并输入目标点云P以及源点云Q。
示例性的:源点云Q构建方法如下:
使用CAD模型离散出源点云Q,并通过PCL库中的NormalEstimation函数去估算源点云的法向量,设置质心为视点,并计算使每个法向量的方向相对质心一致向外,完成法向量定向。
示例性的目标点云P使用结构光三维扫描仪扫描工件得到测量点云,将其作为目标点云P。
步骤2、利用KD-Tree算法搜索目标点云和源点云中最近邻点对集合,目标点云中的一个测点pi与源点云中与该测点pi距离最近的一个点qj组成一对最近邻点对,并基于点对的位置关系将目标点云中的点区分为正负测点。
更进一步,步骤2中,区分正负测点的方法为:目标点云中的测点位置减去源点云中的最近邻点位置得到的向量与源点云相应点的法向量求点积,若其值大于0,则为正余量测点,若小于0则为负余量测点;计算公式如下:
如果测点pi满足,则测点pi为正向测点,/>为源点云qj处的法向量;
如果测点pi满足,则测点pi为负向测点。
步骤3、计算最近邻点对集合中每个点对中P到Q的法向距离,输入函数变形参数α值,设置权重尺度因子c值。
示例性的,步骤3中具体包括以下步骤:
步骤3.1、点对法向距离diT为点对中点到平面距离(第i个测点对应的点对中,测点pi到源点云qj所在切平面之间的距离,后面简称点到平面距离或者点对法向距离),包含正向测点对法向距离dkT和负向测点对法向距离dsT,i∈[1,m+n],i为顺序角标,diT表示第i个测点的点对法向距离,m为正向测点的总数,n为负向测点的总数;
定义δ∈R3、t∈R3分别为旋转向量跟平移向量,R3表示三维向量空间,建立点对法向距离函数(测点pi到相应源点云qj所在切平面的距离),如公式(1):
其中,pi表示目标点云P中的第i个测点,pi+为测点pi经过单步转换后的测点,qj表示源点云Q中与pi组成最邻近点对的点,F为距离函数,nj表示第j个源点云qj处的法向量,niT=nj,是转换向量,第一次迭代时,转换向量ξ=06×1;
步骤3.2、根据点云结构特征输入函数变形参数α值,设置权重尺度因子c值。
为控制权重函数鲁棒性的变形参数,当输入目标点云噪声较高的时候,α设置为负数会更好,而当输入目标点云为低噪声时,α设置为0时会更好,一般情况下设置α=-2。
权重尺度因子的初值一般设置为0,即c=10,之后每隔4次迭代对c值进行一次退火处理,即c=c/2,当c<0.5时取c=0.5,c值越小,参与计算的异常点云就越少,迭代后期应使c值处于一个较小的值。
步骤4、构建鲁棒函数,根据点对的距离将目标点云Q中的点区分为正常测点与异常点云,并通过鲁棒函数对点对距离施加不同程度的权重值。
实施例性的,步骤4具体方法如下:
步骤4中,将异常点云分为可控异常测点和完全异常测点,判断方法如下:
计算法向距离与平均法向距离比值;
当该比值小于或等于1时,判断为正常测点,权重赋值为1,即wk=1;
当该比值大于1时但是小于或等于权值曲线的零点时,判断为可控异常测点,根据异常偏离程度权重对赋值;
当该比值大于权值曲线的零点时,完全异常测点,对完全异常测点的权重赋值0;
具体计算公式如下:
对于正向测点,构建正向测点对距离鲁棒函数如公式(2):
其中wk为第k个正向测点的距离权重函数,当wk<0时修正为wk=0,dkT为第k个正向测点到平面距离,是经单步转换后的所有正向测点到平面距离均值,若为第一次迭代则转换向量ξ=06×1,m是正向测点的总数量,k∈[1,m],α∈R为控制权重函数鲁棒性的变形参数;c∈(0,∞]为控制权重函数零点值的权重尺度因子;d0+表示正向测点权重曲线的零点,如图2中(a)所示,每个正向测点权重曲线与水平轴交点均为零点;
对于负向测点,构建负向测点对距离鲁棒函数如公式(3):
其中ws为第s个负向测点的距离权重函数,当ws<0时修正为ws=0,dsT为第s个负向测点到平面距离,是经单步转换后的负向测点到平面距离均值,若为第一次迭代则转换向量ξ=06×1,n是负向测点的总数量,s∈[1,n],α∈R为控制权重函数鲁棒性的变形参数;d0-表示负向测点权重曲线的零点,如图2中(b)所示,每个负向测点权重曲线与水平轴交点均为零。
步骤5、基于点对法向距离与鲁棒函数建立鲁棒权重加权方差最小化的目标函数。
实施例性的,建立鲁棒函数加权方差最小化算法的目标函数如公式(4):
其中,为第k个正向测点的正向鲁棒函数加权法向距离偏差,满足 为第s个负向测点的负向鲁棒函数加权法向距离偏差,满足为所有正向测点鲁棒函数加权法向距离均值,/>为所有负向测点鲁棒函数加权法向距离均值;dkT为正向测点对法向距离,dsT为负向测点对法向距离;
考虑异常点云,将公式(4)改写表示为公式(5):
其中,与/>为正负向在正常范围内的测点的权值,其值为1;m1、m2分别表示正常正向测点总数和异常正向测点总数;n1、n2分别表示正常负向测点总数和异常负向测点总数;/>为所有异常测点的加权距离偏差平方和,包括正向异常测点与负向异常测点两部分,Df为第f个异常测点的加权距离偏差;所有正向测点鲁棒函数加权法向距离均值/>的计算如公式(6)所示:
其中,m21和m20分别为可控异常测点和完全异常测点的总数量,为第f1+个可控异常测点的鲁棒函数权重,/>为第f0+个完全异常测点的鲁棒函数权重,/>为第f1+个可控异常测点对应的点对法向距离,/>为第f0+个完全异常测点对应的点对法向距离,由步骤5中所建立的鲁棒权重函数,得到/>将其与/>代入公式(6)中,得到如公式(7)所示:
对上述目标函数具有鲁棒性验证如下:
公式(7)为实际的表达式,设理想的正向加权距离均值/>并计算/>与/>的差值如公式(8)所示:
由公式(8)可知,假设全部的与/>值非常相近时,/>的值接近于0,即实际的加权平均距离与理想值非常接近;若/>为(0,1)之间的一个居中值,那么由上式可知,若存在过多的可控异常测点时,即算法将余量不大的大量偏差点云区域视为可控余量参与计算,m21>>m20甚至m21与m1是一个数量级的情况下,/>的值会非常大,此时理想正向收敛位置与实际收敛位置存在较大偏差,从而导致匹配失真。而在步骤5中设定的权重函数wk会随着迭代次数缩小可控余量点云的距离范围来减小m21的值,同时还能缩小/>与/>之间的差值,此时正向实际收敛位置/>会逐渐向正向理想收敛位置/>靠近,因此RFWVM算法(面向复杂构件视觉测量的鲁棒函数加权方差最小化算法)可以抑制大量异常点云引起的匹配失真。将Df中可控异常测点剔除得到/>假设RFWVM算法使/>同理,可证/>将其与/>一起代入/>中,如公式(9):
由公式(9)可以证明,在权重的作用下,不可控异常测点不会参与ξRFWVM的求解,且此时等价于理想正负余量均值距离,则最小化/>的目的与配准的目的是一致的,因此目标函数具有鲁棒性。
步骤6、根据鲁棒权重加权方差最小化的目标函数求解源点云Q和目标点云P之间精配准的刚体转换矩阵,利用计算出来的刚体转换矩阵作用于源点云、以更新源点云位置。
示例性的,具体步骤如下:
步骤6.1、将正向测点的鲁棒函数加权法向距离偏差平方和化简,如公式(10):
上式中,E+是6×6矩阵,F+是6×1矩阵,D+为标量;
步骤6.2、将负向测点的鲁棒函数加权法向距离偏差平方和进行同样化简,如公式(11):
其中,E-是6×6矩阵,F-是6×1矩阵,D-为标量;
步骤6.3、根据公式(4)所示目标函数G(R,t),使用目标函数对转换向量ξ求导并求解,如公式(12):
上式中, Ak,As均为1×6的矩阵,/>所有正向测点的点对距离均值,/>是未经单步转换的第k个正向测点的点对距离;/>所有负向测点的点对距离均值,/>是未经单步转换的第s个负向测点的点对距离;/>鲁棒函数加权负向平均向量,/>鲁棒函数加权正向平均向量;
故旋转矩阵与平移矩阵求得如公式(13):
e为自然底数,t为平移向量,为旋转向量δ的反对称矩阵;
当δ=[δx,δy,δz]时,
步骤6.4、将R与T作用与源点云,更新源点云位置。
步骤7、定义余量不均点云、存在结构偏差点云时的误差评价函数。
示例性,误差评价函数如公式(14):
公式(14)中,pi表示配准后的测点中点对距离小于平均点对距离的测点,qj表示从源点云Q中搜到的最近点。
步骤8、重复步骤2-7,并每隔4次迭代更新一次c的值,直到达到收敛条件,最终获得精配准的刚体转换矩阵T。
需要说明的是,迭代间隔次数K不限于上述取值4,可以才2-8之间根据需要任意取值,取值为4仅仅是本实施例中计算精度和算法的平衡。
示例性的,上述收敛条件为迭代次数M大于设定值Mmax或连续两次迭代误差评价函数值相差小于设定阈值G。
需要说明的是,Mmax设定值越大越好,但是考虑执行效率和所需算力问题,一般为100-500之前,可以满足本发明要求。阈值G取值越小精度越高,但是考虑执行效率和所需算力问题,一般设置在1×10-6左右;可以为1×10-5~1×10-7。
手动将叶片点云模型25.8%的点进行2mm的偏置,并对源点云进行一定的位姿转换,作为配准开始前两片点云的初始位置如图3所示;利用WPMAVM算法对叶片点云迭代配准,在第14次左右就陷入了局部最优解,最终其将正负余量权重换成RFWVM算法中的鲁棒权重函数后,算法可正常收敛,并且其/>十分接近于0,配准基本达到全局最优,其误差变化曲线如图4所示。后续加入ICP算法对同样的叶片点云进行配准计算,将各算法结果中的最近邻点距离映射为色谱图,如图5所示,可以看出,RFWVM算法除了手动偏置部分,其余正常点云均配准到正确位置,相对其他两种算法稳定性更强,其最终/>值对比如图6所示。
后续使用冲压件车身的模型点云与测量点云进行异源点云配准实验,其初始状态见附图7,同样使用三种算法进行配准,最终效果如图8所示,观察图8可知,RFWVM算法更能有效的区分出哪些点云为异常点云,并加以抑制,其配准效果相对其他两种算法更准确。
需要说明的是,本发明使用单步转换计算技术可以参考CN202110573411.X。
以上实施方式仅用于说明本发明,而非对本发明的限制。尽管参照实施例对本发明进行了详细说明,本领域的普通技术人员应当理解,对本发明的技术方案进行各种组合、修改或者等同替换,都不脱离本发明技术方案的精神和范围,均应涵盖在本发明的权利要求范围当中。
Claims (10)
1.一种面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,包括如下步骤:
步骤1、构建并输入目标点云P以及源点云Q;
步骤2、利用KD-Tree算法搜索目标点云和源点云中最近邻点对集合,目标点云中的一个测点pi与源点云中与该测点pi距离最近的一个点qj组成一对最近邻点对,并基于点对的位置关系将目标点云中的点区分为正负测点;
步骤3、计算最近邻点对集合中每个点对中pi到qj的法向距离,输入函数变形参数α值,设置权重尺度因子c值;
步骤4、构建鲁棒函数,根据点对的距离将目标点云Q中的点区分为正常测点与异常点云,并通过鲁棒函数对点对距离施加不同程度的权重值;
步骤5、基于点对法向距离与鲁棒函数建立鲁棒权重加权方差最小化的目标函数;
步骤6、根据鲁棒权重加权方差最小化的目标函数求解源点云Q和目标点云P之间精配准的刚体转换矩阵,利用计算出来的刚体转换矩阵作用于源点云、以更新源点云位置;
步骤7、定义余量不均点云、存在结构偏差点云时的误差评价函数;
步骤8、重复步骤2-7,并每隔K次迭代更新一次c的值,直到达到收敛条件,最终获得精配准的刚体转换矩阵T。
2.根据权利要求1所述面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,步骤1中,源点云Q构建方法如下:
使用CAD模型离散出源点云Q,并通过PCL库中的NormalEstimation函数去估算源点云的法向量,设置质心为视点,并计算使每个法向量的方向相对质心一致向外,完成法向量定向;
所述目标点云P使用结构光三维扫描仪扫描工件得到。
3.根据权利要求1所述面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,步骤2中,区分正负测点的方法为:目标点云中的测点位置减去源点云中的最近邻点位置得到的向量与源点云相应点的法向量求点积,若其值大于0,则为正余量测点,若小于0则为负余量测点。
4.根据权利要求1所述面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,步骤3中具体包括以下步骤:
步骤3.1、点对法向距离diT为点对中点到平面距离,包含正向测点对法向距离dkT和负向测点对法向距离dsT,i∈[1,m+n],i为顺序角标,diT表示第i个测点的点对法向距离,m为正向测点的总数,n为负向测点的总数;
定义δ∈R3、t∈R3分别为旋转向量跟平移向量,R3表示三维向量空间,建立点对法向距离函数,如公式(1):
其中,pi表示目标点云P中的第i个测点,pi+为测点pi经过单步转换后的测点,qj表示源点云Q中与pi组成最邻近点对的点,F为距离函数,nj表示第j个源点云qj处的法向量,niT=nj,是转换向量;
步骤3.2、根据点云结构特征输入函数变形参数α值,设置权重尺度因子c值。
5.根据权利要求4所述面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,步骤4中,将异常点云分为可控异常测点和完全异常测点,并对完全异常测点的权重赋值0,根据正向测点和负向测点分别计算和判断如下:
对于正向测点,构建正向测点对距离鲁棒函数如公式(2):
其中wk为第k个正向测点的距离权重函数,当wk<0时修正为wk=0,dkT为第k个正向测点到平面距离,是经单步转换后的所有正向测点到平面距离均值,第一次迭代的转换向量ξ=06×1,m是正向测点的总数量,k∈[1,m],α∈R为控制权重函数鲁棒性的变形参数;c∈(0,∞]为控制权重函数零点值的权重尺度因子;d0+表示正向测点权重曲线的零点;
对于负向测点,构建负向测点对距离鲁棒函数如公式(3):
其中ws为第s个负向测点的距离权重函数,当ws<0时修正为ws=0,dsT为第s个负向测点到平面距离,是经单步转换后的负向测点到平面距离均值,第一次迭代的转换向量ξ=06×1,n是负向测点的总数量,s∈[1,n],α∈R为控制权重函数鲁棒性的变形参数;d0-表示负向测点权重曲线的零点。
6.根据权利要求5所述面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,步骤5中,建立鲁棒函数加权方差最小化算法的目标函数如公式(4):
其中,为第k个正向测点的正向鲁棒函数加权法向距离偏差,满足/>dswsRM_为第s个负向测点的负向鲁棒函数加权法向距离偏差,满足为所有正向测点鲁棒函数加权法向距离均值,/>为所有负向测点鲁棒函数加权法向距离均值;dkT为正向测点对法向距离,dsT为负向测点对法向距离;
考虑异常点云,将公式(4)改写表示为公式(5):
其中,与/>为正负向在正常范围内的测点的权值,其值为1;m1、m2分别表示正常正向测点总数和异常正向测点总数;n1、n2分别表示正常负向测点总数和异常负向测点总数;为所有异常测点的加权距离偏差平方和,包括正向异常测点与负向异常测点两部分,Df为第f个异常测点的加权距离偏差;所有正向测点鲁棒函数加权法向距离均值/>的计算如公式(6)所示:
其中,m21和m20分别为可控异常测点和完全异常测点的总数量,为第f1+个可控异常测点的鲁棒函数权重,/>为第f0+个完全异常测点的鲁棒函数权重,/>为第f1+个可控异常测点对应的点对法向距离,/>为第f0+个完全异常测点对应的点对法向距离,由步骤5中所建立的鲁棒权重函数,得到/>将其与/>代入公式(6)中,得到如公式(7)所示:
7.根据权利要求6所述面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,步骤6具体包括以下步骤:
步骤6.1、将正向测点的鲁棒函数加权法向距离偏差平方和化简,如公式(10):
上式中,E+是6×6矩阵,F+是6×1矩阵,D+为标量;
步骤6.2、将负向测点的鲁棒函数加权法向距离偏差平方和进行同样化简,如公式(11):
其中,E_是6×6矩阵,F_是6×1矩阵,D_为标量;
步骤6.3、根据公式(4)所示目标函数G(R,t),使用目标函数对转换向量ξ求导并求解,如公式(12):
上式中, Ak,As均为1×6的矩阵,/>所有正向测点的点对距离均值,/>是未经单步转换的第k个正向测点的点对距离;/>所有负向测点的点对距离均值,dsT0是未经单步转换的第s个负向测点的点对距离;/>鲁棒函数加权负向平均向量,/>鲁棒函数加权正向平均向量;
故旋转矩阵与平移矩阵求得如公式(13):
e为自然底数,t为平移向量,为旋转向量δ的反对称矩阵;
步骤6.4、将R与T作用与源点云,更新源点云位置。
8.根据权利要求7所述面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,步骤7中,误差评价函数如公式(14):
公式(14)中,pi表示配准后的测点中点对距离小于平均点对距离的测点,qj表示从源点云Q中搜到的最近点。
9.根据权利要求1-8任意一项所述面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,步骤8中,K取值范围为2-8。
10.根据权利要求1-8任意一项所述面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法,其特征在于,步骤8中,收敛条件为迭代次数M大于设定值Mmax或连续两次迭代误差评价函数值相差小于设定阈值G。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310608783.0A CN116718137A (zh) | 2023-05-26 | 2023-05-26 | 面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310608783.0A CN116718137A (zh) | 2023-05-26 | 2023-05-26 | 面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN116718137A true CN116718137A (zh) | 2023-09-08 |
Family
ID=87867030
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202310608783.0A Pending CN116718137A (zh) | 2023-05-26 | 2023-05-26 | 面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN116718137A (zh) |
-
2023
- 2023-05-26 CN CN202310608783.0A patent/CN116718137A/zh active Pending
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104867136B (zh) | 一种基于距离方差最小的工件点云匹配算法 | |
CN104697462B (zh) | 一种基于中轴线的航空叶片型面特征参数提取方法 | |
US7619749B2 (en) | Methods and apparatus of aligning surfaces | |
CN109147040B (zh) | 基于模板的人体点云孔洞修补方法 | |
CN106446343B (zh) | 一种径流叶轮直纹叶片参数化型线的自动提取方法 | |
Makem et al. | A virtual inspection framework for precision manufacturing of aerofoil components | |
CN107091627B (zh) | 一种多表面系统的综合测量与评估方法 | |
CN113516695A (zh) | 激光轮廓仪平面度测量中的点云配准策略 | |
CN113910001B (zh) | 一种数控机床空间误差辨识方法 | |
CN109773593B (zh) | 一种基于余量约束条件下的磨削方法 | |
Lv et al. | WPMAVM: Weighted plus-and-minus allowance variance minimization algorithm for solving matching distortion | |
CN113029049B (zh) | 基于加权正负余量方差最小化算法的工件光学测量方法 | |
CN113094964B (zh) | 一种生成叶片加工坐标的方法和装置 | |
Chen et al. | Fast registration of 3D point clouds with offset surfaces in precision grinding of free-form surfaces | |
Pechenin et al. | Development of a method of ICP algorithm accuracy improvement during shaped profiles and surfaces control | |
CN116718137A (zh) | 面向复杂构件视觉测量的鲁棒函数加权最小化匹配算法 | |
RU2607867C2 (ru) | Процесс адаптивной обработки литых лопаток | |
CN109035238B (zh) | 一种面向自由曲面零件的加工余量离线分析方法 | |
Wu et al. | Rigid shape matching for 3-D robotic grinding measurement with applications to blades | |
Tan et al. | A fast and differentiated localization method for complex surfaces inspection | |
Wan et al. | The machining surface localization of free-form blade considering form tolerance | |
CN113570147B (zh) | 一种薄壁件多工况加工误差快速预测方法及设备 | |
Hsu et al. | An iterative coordinate setup algorithm for airfoil blades inspection | |
CN114898069A (zh) | 一种曲面蒙皮三维数据拼合方法 | |
CN115164809A (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 |