CN111462202B - 非刚性注册方法及系统 - Google Patents

非刚性注册方法及系统 Download PDF

Info

Publication number
CN111462202B
CN111462202B CN202010269838.6A CN202010269838A CN111462202B CN 111462202 B CN111462202 B CN 111462202B CN 202010269838 A CN202010269838 A CN 202010269838A CN 111462202 B CN111462202 B CN 111462202B
Authority
CN
China
Prior art keywords
model
function
initial model
point
points
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
CN202010269838.6A
Other languages
English (en)
Other versions
CN111462202A (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.)
University of Science and Technology of China USTC
Original Assignee
University of Science and Technology of China USTC
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 University of Science and Technology of China USTC filed Critical University of Science and Technology of China USTC
Priority to CN202010269838.6A priority Critical patent/CN111462202B/zh
Publication of CN111462202A publication Critical patent/CN111462202A/zh
Application granted granted Critical
Publication of CN111462202B publication Critical patent/CN111462202B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • G06T7/344—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods involving models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/14—Transformations for image registration, e.g. adjusting or mapping for alignment of images
    • G06T3/147—Transformations for image registration, e.g. adjusting or mapping for alignment of images using affine transformations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种非刚性注册方法及系统，相关方案通过建立嵌入的采样图模型，来降低优化时间，采用拟牛顿方法求解替代函数使得求解更快；通过参数值的调整策略逐渐增加稀疏性使得模型更加鲁棒，总而言之，本发明能够有效的解决基于迭代最近点的非刚性注册时，求解速度慢、精度不高等缺点。

Description

非刚性注册方法及系统
技术领域
本发明涉及三维技术领域，尤其涉及一种非刚性注册方法及系统。
背景技术
由于三维扫描技术的飞速发展，通过扫描得到深度数据来重建三维物体或场景以及三维物体的动态重建及跟踪等方面应用越来越广泛。由于深度扫描设备每次只能测量得到物体局部坐标系下的部分表面，并且可能出现平移错位或旋转错位，以及有大量的噪声，因此为了得到物体完整表面数据，需要对这些局部点云数据进行整合和配准。
对于不同视角下的静止物体，一般通过点云刚性注册来寻找一种三维空间中的刚性变换使得在不同视角下的三维坐标点共同部分能够正确的匹配和拼接。然而对于动态的物体，由于局部变形的存在，仅仅通过一个刚性变换难以将两个或多个表面数据良好匹配，因此发展了非刚性注册技术，即对一个表面数据采用多个变换矩阵更精准的表达。
目前的非刚性注册方法很多，一些方法根据两个表面的弧度、曲率等几何特性全局搜索点或局部表面的对应关系，以得到一个全局的变换将二者良好匹配，但由于几何特征的歧义性，寻找全局准确的匹配点仍然是一个很难的问题；一些方法基于最近点迭代的思想，将最近点作为对应点迭代找到最优解，这种方法为局部注册方法，依赖于初始状态，并易受噪声、异常值等影响。由于非刚性注册问题对一个表面赋予了很多个局部变换，由此带来的自由度高使得表面可能扭曲很大，因此为了使得表面依然具有连续性和光滑性，一般引入正则项来约束变换的求解，也因此带来目标函数复杂度高的问题，无法让求解更加高效。
因此，如何有效的解决非刚性注册中求解速度慢、精度不高等缺点，是一项亟待解决的问题。
发明内容
本发明的目的是提供一种非刚性注册方法及系统，能够提高求解速度，并提升精度。
本发明的目的是通过以下技术方案实现的：
一种非刚性注册方法，包括：
步骤S101、基于输入的初始模型，根据给定的采样半径，得到嵌入简化图模型；
步骤S102、基于所述初始模型、嵌入简化图模型与输入的目标模型，进行最近邻搜索，得到匹配的对应点，构造原目标函数的替代函数；
步骤S103、采用拟牛顿方法迭代求解替代函数，得到非刚性变换矩阵；
步骤S104、判断所述非刚性变换矩阵是否收敛，若未收敛，则返回步骤S102更新替代函数；若收敛，则进入步骤S105；
步骤S105、判断是否达到参数值下界，若未达到参数值下界，将参数值减小后，返回步骤S102更新替代函数；若达到参数值下界，则转入步骤S106；
步骤S106、终止迭代，得到非刚性变换矩阵并将之应用到初始模型点上得到变换后的初始模型。
由上述本发明提供的技术方案可以看出，通过建立嵌入的采样图模型，来降低优化时间，采用拟牛顿方法求解替代函数使得求解更快；通过参数值的调整策略逐渐增加稀疏性使得模型更加鲁棒，总而言之，本发明能够有效的解决基于迭代最近点的非刚性注册时，求解速度慢、精度不高等缺点。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。
图1为本发明实施例提供的一种非刚性注册方法的流程图；
图2为本发明实施例提供的一种非刚性注册系统的流程图。
具体实施方式
下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。
本发明实施例提出一种非刚性注册方法，如图1所示，其主要包括：
步骤S101、基于输入的初始模型，根据给定的采样半径，得到嵌入简化图模型。
本发明实施例中，初始模型可以是可以计算测地线距离的模型，例如网格模型。
由于对初始模型中每个点都赋予一个仿射变换常常导致变量数多，优化时间长，效率低。因此，本发明实施例中，通过给定的采样半径从所述初始模型中提取一部分的节点构成一个变形图模型，也即嵌入简化图模型，再通过嵌入简化图模型中节点的变换带动所述初始模型中所有点的变换，主要过程如下：
1)根据给定的半径Ra在所述初始模型中迭代的提取采样点：对于初始模型中所有的点，将执行PCA提取协方差矩阵的最大特征值所在方向作为一个轴，将初始模型中所有的点投影到相应轴上并对投影值进行从小到大排序，根据排序添加采样点：如果采样点集合为空集或者初始模型中当前点到采样点集合中所有已存在点的测地线距离最小值大于给定半径Ra，则将当前点加入采样点集合，构成嵌入简化图的节点集。
2)连接测地线距离小于给定半径的任意两个节点构成嵌入简化图的边集，由嵌入简化图的节点集与边集构成嵌入简化图模型。
3)将嵌入简化图的点集中第j个节点处的仿射变换表示为(Aj,tj)；其中，
Figure BDA0002442730500000031
表示变换矩阵，
Figure BDA0002442730500000032
表示平移部分。如果初始模型中的第i个点在第j个节点的影响范围内，则用Aj(vi-pj)+pj+tj表示在变换(Aj,tj)作用下的新位置，其中
Figure BDA0002442730500000033
为第j个节点的三维坐标，
Figure BDA0002442730500000034
为第i个初始模型点的三维坐标，则第i个初始模型点在变形后的位置表示为：
Figure BDA0002442730500000035
其中，I(vi)表示影响第i个初始模型点的嵌入简化图节点集合，I(vi)＝{pj|D(vi,pj)＜Ra}，其中D(vi,pj)表示vi与pj之间的测地线距离；wij为权重，它是依赖距离的归一化权重：
Figure BDA0002442730500000036
步骤S102、基于所述初始模型、嵌入简化图模型与输入的目标模型，进行最近邻搜索，得到匹配的对应点，构造原目标函数的替代函数。
本步骤中，基于所述初始模型、嵌入简化图模型与目标模型，通过对所述目标模型建立分割k维数据空间的数据结构进行最近邻搜索，得到匹配的对应点，并构造原目标函数的替代函数。
由于传统的基于迭代最近点的非刚性算法存在易受异常值和噪声影响的问题，而改进的方法虽然可以有效减少噪声影响，但是求解速度很慢，不能实时优化。为此，本发明实施例提出采用Welsch函数改进逐点对齐项以及正则项，并用MM算法(Majorization-minimization，最大化-最小化方法)求解，首先找到目标函数的上界函数(也即替代函数)，其次采用L-BFGS方法求解上界函数，并且动态调整Welsch函数的参数，以逐渐增强稀疏性，有效改善传统的非刚性算法。
目标函数表示为：
min Ealign(X)+αEreg(X)+βErot(X)
其中，α和β分别为正则项和旋转项的权重，均为正数，用来平衡三项能量。具体实施中，默认设置
Figure BDA0002442730500000041
其中n为初始模型的点数，r为嵌入简化图节点数，|Eg|为嵌入简化图边数，可根据输入的初始模型和目标模型的形状及位置差异来调节，若差异较大，则将β值调大，否则调小，若需要的光滑程度高，则将α值调大，否则调小。
Figure BDA0002442730500000042
表示对齐项，
Figure BDA0002442730500000043
的计算方式参照前文，即可以通过嵌入简化图模型计算，va为对齐项参数，使得对应点尽可能接近，同时对于距离太远的点有忽略作用以使得对应点的准确性提高；
Figure BDA0002442730500000044
为正则项，表示相邻第j个和第l个嵌入图节点的变换差Djl尽可能一致，Dji可以通过嵌入简化图模型得到，νr为正则项参数，
Figure BDA0002442730500000045
希望保持局部表面的刚性；则上式可以表示为：
Figure BDA0002442730500000046
其中，
Figure BDA0002442730500000047
v为指定参数，即在上述优化问题中的va和vr，后面给出具体的选取方式，x表示上式中的
Figure BDA0002442730500000048
||Djl||；
Figure BDA00024427305000000411
分别为初始模型和目标模型，
Figure BDA00024427305000000412
分别为初始模型
Figure BDA00024427305000000413
和目标模型
Figure BDA00024427305000000414
上点的三维坐标的集合，m为目标模型的点数；初始模型
Figure BDA00024427305000000415
上第i个点vi与目标模型
Figure BDA00024427305000000416
上的第ρ(i)个点uρ(i)对应，该对应关系是通过对所述目标模型建立分割k维数据空间的数据结构，并根据初始模型上点的三维坐标进行最近邻搜索，找到目标模型上欧式距离最近的点为对应点。
Figure BDA0002442730500000049
表示待求解的非刚性变换矩阵，
Figure BDA00024427305000000410
Figure BDA0002442730500000051
表示仿射变换，r为节点数，uρ(i)表示vi的对应点，即目标模型上与vi欧式距离最近的点；Djl＝Aj(pl-pj)+pj+tj-(pl+tl)表示第j个节点与相邻第l个节点的变换差，其中pl,
Figure BDA0002442730500000052
分别为第l个节点和第j个节点的三维坐标，tl表示第l个节点处的仿射变换的平移部分；N(pj)为第j个节点的邻域集合；
Figure BDA0002442730500000053
为旋转群，其中，
Figure BDA0002442730500000054
为一个单位阵，det(R)表示求解的行列式值；||·||F表示Frobenius范数，
Figure BDA0002442730500000055
表示投影到SO3的投影算子。
为更好的求解，采用最大化-最小化方法进行求解，具体在于，首先构造待求解函数的上界函数，其次优化上界函数，以达到优化原函数的目的，二者交替迭代得到最优解，下面具体阐述求解过程。
由于Welsch函数该函数的任意一点y都可以找到一个凸二次上界函数：
Figure BDA0002442730500000056
使得
当x≠y时，
Figure BDA0002442730500000057
当x＝y时，
Figure BDA0002442730500000058
其中，x为函数的自变量，y为任意选取的点，v与上述一致，为Welsch函数参数；忽略与自变量无关的量，即
Figure BDA0002442730500000059
记为C，则上式等价于
Figure BDA00024427305000000510
忽略C，将之应用到对齐项，设置该项的Welsch参数为va,并根据第k步迭代值设置
Figure BDA00024427305000000511
则可得到
Figure BDA00024427305000000512
将之应用到正则项，设置该项的Welsch参数为vr，并根据第k步迭代值设置
Figure BDA00024427305000000513
则可得到
Figure BDA00024427305000000514
则可构造上界函数作为替代函数：
Figure BDA00024427305000000515
步骤S103、采用拟牛顿方法迭代求解替代函数，得到非刚性变换矩阵，即新一步的迭代值为：
Figure BDA00024427305000000516
前文构造的替代函数中的Erot(X)为高度的非线性函数，因此具体采用拟牛顿方法中的L-BFGS进行求解：
1)计算替代函数的梯度G(X)，公式如下：
Figure BDA0002442730500000061
矩阵形式可以表示为
Figure BDA0002442730500000062
其中
Figure BDA0002442730500000063
将
Figure BDA0002442730500000064
写为矩阵形式为：
Figure BDA0002442730500000065
其中
Figure BDA0002442730500000066
Figure BDA0002442730500000067
为存储权重，
Figure BDA0002442730500000068
是一个块矩阵，且当pj∈I(vi)时，
Figure BDA0002442730500000069
否则Fij＝[0,0,0,0]；
Figure BDA00024427305000000610
Figure BDA00024427305000000611
同样的，
Figure BDA00024427305000000612
其中，
Figure BDA00024427305000000613
Wr为存储权重
Figure BDA00024427305000000614
Y每行存储
Figure BDA00024427305000000615
B对应于xl位置的块矩阵为
Figure BDA00024427305000000616
和对应于xj位置的块矩阵为[0，0，0，1]。
则可得到替代函数
Figure BDA00024427305000000617
的梯度为：
Figure BDA00024427305000000618
固定替代函数作为目标函数，将当前第k步迭代的迭代值Xk作为拟牛顿求解的初值
Figure BDA00024427305000000619
通过循环计算梯度以及下降方向并更新迭代值
Figure BDA00024427305000000620
的过程，直至满足
Figure BDA00024427305000000621
时终止，将得到的迭代值
Figure BDA00024427305000000622
即作为最优解Xk+1。具体计算如下：
根据初始迭代值
Figure BDA00024427305000000623
得到初始近似Hessian矩阵为：
Figure BDA00024427305000000624
H0用来计算下降方向，进而得到
Figure BDA00024427305000000625
不断的循环计算，在第t次循环计算过程中，得到第t次迭代值
Figure BDA00024427305000000626
同时，根据第t-1次迭代值
Figure BDA00024427305000000627
计算梯度值
Figure BDA00024427305000000628
并通过L-BFGS方法计算下降方向，即结合历史迭代值
Figure BDA00024427305000000629
和梯度值
Figure BDA00024427305000000630
得到第t次迭代的近似Hessian矩阵Ht，即：
Figure BDA00024427305000000631
Figure BDA0002442730500000071
其中，
Figure BDA0002442730500000072
Figure BDA0002442730500000073
I为单位阵，Tr为矩阵迹函数。
得到下降方向
Figure BDA0002442730500000074
并通过线搜索得到最优步长λ＞0并满足Armijo条件
Figure BDA0002442730500000075
得到t+1步新的非刚性变换矩阵
Figure BDA0002442730500000076
不断迭代直到达到终止条件时终止，得到固定当前替代函数下的最优非刚性变换矩阵Xk+1，默认设置阈值∈1＝10-3，预设值mmax＝5。
步骤S104、判断当前迭代计算的非刚性变换矩阵是否收敛，若未收敛，则返回步骤S102，即根据当前迭代计算的非刚性变换矩阵更新替代函数；若收敛，则进入步骤S105。
利用相邻两次迭代(第k、第k+1次)构造的替代函数所求解出的最优非刚性变换矩阵对初始模型所有顶点的进行变换，逐一计算变换后顶点位置欧式距离最大值，并判断：
Figure BDA0002442730500000077
是否减小到小于一个用户指定的值∈2，若是则认为收敛；否则，未收敛，一般设置阈值∈2＝10-3；未收敛时，需要范围步骤S102，利用上一次迭代得到的非刚性变换矩阵来构造新的替代函数。
步骤S105、判断是否达到参数值下界，若未达到参数值下界，将参数值减小后，返回步骤S102更新替代函数；若达到参数值下界，则转入步骤S106。
本发明实施例中，由于在初始的最近点作为对应点非常不准确，由于迭代次数的增加，最近点作为匹配点准确性越来越高，因此我们将对齐项Welsch函数参数va逐渐减小，对于正则项，同样随着迭代次数的增加，相邻嵌入简化图节点变换的误差越来越小，因此我们也将vr逐渐减小以增加函数的稀疏性，具体策略为：
首先令
Figure BDA0002442730500000078
固定va和vr达到收敛之后将令
Figure BDA0002442730500000079
返回步骤S102重新构建目标函数的替代函数，不断迭代直至
Figure BDA00024427305000000710
时迭代终止。其中，设置对齐项参数值va上界为
Figure BDA00024427305000000711
Figure BDA00024427305000000712
为开始时初始模型中所有点到目标模型中的对应点距离的中值，其中目标模型中的对应点为距离初始模型中点的欧氏距离最近的点，va下界通过初始模型的平均边长
Figure BDA00024427305000000713
决定，设为
Figure BDA00024427305000000714
正则项参数vr最大值设为
Figure BDA00024427305000000715
其中
Figure BDA00024427305000000716
为用户指定的参数，默认设置
Figure BDA00024427305000000717
步骤S106、终止迭代，得到非刚性变换矩阵并将之应用到初始模型点上得到变换后的初始模型。
需要说明的，本发明实施例所给出各项参数的具体数值仅为示例，并非构成限制；在实际应用中，用户可以根据实际情况或者经验调整各项参数的具体数值。
本发明实施例上述方案中，通过采样方法得到嵌入的嵌入简化图模型，通过对应点对齐项，相邻节点变换正则项，节点变换旋转项三项的线性组合来构建目标函数，通过最大化-最小化方法(MM方法)构造替代函数并通过L-BFGS方法求解，不断迭代得到最优非刚性变换矩阵。
由于传统目标函数采用变换后的初始模型和目标模型逐点之间的L2范数作为距离度量，这样使得注册结果很容易受到噪声和异常值的影响，由于Welsch函数具有更加稀疏的特性，能够使得距离较远的对应点对对目标函数整体的影响较小，因此发明公开的方法实施例采用Welsch函数作为逐点之间的距离度量；由于相邻节点之间可能存在大的局部变形，比如人体模型的关节处，因此也采用Welsch函数作为相邻节点变换差的度量以允许局部大变形的存在；
对于本发明实施例提供的目标函数，使用最大化-最小化(MM)算法来构造替代函数：1)找到当前目标函数的一个上界函数(在当前值处与原函数值相同，其他点处均大于原函数值)；2)极小化替代函数以达到极小化原目标函数的效果，由于待求函数的非线性性，采用L-BFGS方法来求解。
由于Welsch函数的目标函数需要设置参数，而不同规模的模型需要的参数值不同，寻找适合该模型的参数值比较麻烦，而且随着注册过程的进行，固定的参数不能达到很好的效果，因此本发明提出动态调整参数的方法来解决该问题。首先初始参数根据初始的对应点逐点距离的中值来选取，当迭代收敛时，将参数值减半，继续迭代，直到降到终止值及以下迭代终止，终止值根据初始模型的平均边长来选取。
由上述发明公开的方法实施例提供的技术方案可以看出：发明公开的方法实施例通过改变原有的L2范数的能量函数为更加稀疏的目标函数，使得模型不易受异常值及噪声影响以及易于保持局部的尖锐特征；采用拟牛顿方法求解替代函数使得求解更快；通过参数值的调整策略逐渐增加稀疏性使得模型更加鲁棒。
本发明实施例还提供一种非刚性注册系统，该系统用来实现前述的方法，如图2所示，该系统主要包括：
采样与模型构建模块，用于执行步骤S101，即基于输入的初始模型，根据给定的采样半径，得到嵌入简化图模型；
函数构造模块，用于执行步骤S102，即基于所述初始模型、嵌入简化图模型与输入的目标模型，进行最近邻搜索，得到匹配的对应点，构造原目标函数的替代函数；
求解模块，用于执行步骤S103，即采用拟牛顿方法迭代求解替代函数，得到非刚性变换矩阵；
判断模块，用于执行步骤S104，即判断所述非刚性变换矩阵是否收敛，若未收敛，则返回步骤S102更新替代函数；若收敛，则进入步骤S105；
参数值判断与调整模块，用于执行步骤S105，即判断是否达到参数值下界，若未达到参数值下界，将参数值减小后，返回步骤S102更新替代函数；若达到参数值下界，则转入步骤S106；还用于执行步骤S106，即终止迭代，得到非刚性变换矩阵并将之应用到初始模型点上得到变换后的初始模型。
本发明实施例中，通过给定的采样半径从所述初始模型中提取一部分的节点构成一个变形图模型，也即嵌入简化图模型，再通过嵌入简化图模型中节点的变换带动所述初始模型中所有点的变换；
根据给定的半径Ra在所述初始模型中迭代的提取采样点：对于初始模型中所有的点，将执行PCA提取协方差矩阵的最大特征值所在方向作为一个轴，将初始模型中所有的点投影到相应轴上并对投影值进行排序，根据排序添加采样点：如果采样点集合为空或者初始模型中当前点到采样点集合中所有已存在点的测地线距离最小值大于给定半径Ra，则将当前点加入采样点集合，构成嵌入简化图的节点集；
连接测地线距离小于给定半径的任意两个节点构成嵌入简化图的边集，由嵌入简化图的节点集与边集构成嵌入简化图模型；
将嵌入简化图的点集中第j个节点处的仿射变换表示为(Aj,tj)；其中，
Figure BDA0002442730500000091
表示变换矩阵，
Figure BDA0002442730500000092
表示平移部分；如果初始模型中的第i个点在第j个节点的影响范围内，则用Aj(vi-pj)+pj+tj表示在变换(Aj,tj)作用下的新位置，
Figure BDA0002442730500000093
为第j个节点的三维坐标，
Figure BDA0002442730500000094
为第i个初始模型点的三维坐标，则第i个初始模型点在变形后的位置表示为：
Figure BDA0002442730500000101
其中，I(vi)表示影响第i个初始模型点的嵌入简化图节点集合，wij为权重wij。
本发明实施例中，所述基于所述初始模型、嵌入简化图模型与目标模型，进行最近邻搜索，得到匹配的对应点，构造原目标函数的替代函数包括：
基于所述初始模型、嵌入简化图模型与目标模型，通过对所述目标模型建立分割k维数据空间的数据结构进行最近邻搜索，得到匹配的对应点，并构造如下目标函数：
minEalign(X)+αEreg(X)+βErot(X)
其中，
Figure BDA0002442730500000102
表示对齐项，n为初始模型中点的数量；
Figure BDA0002442730500000103
表示正则项，Djl＝Aj(pl-pj)+pj+tj-(pl+tl)表示第j个节点与相邻第l个节点的变换差；
Figure BDA0002442730500000104
表示旋转项；α和β分别为正则项和旋转项的权重；
Figure BDA0002442730500000105
v为指定参数，即va和vr，x表示
Figure BDA0002442730500000106
||Djl||；；
Figure BDA00024427305000001023
分别为初始模型和目标模型，
Figure BDA00024427305000001024
Figure BDA00024427305000001025
分别为初始模型
Figure BDA00024427305000001026
和目标模型
Figure BDA00024427305000001027
上点的三维坐标的集合，
Figure BDA0002442730500000107
表示待求解的非刚性变换矩阵，
Figure BDA0002442730500000108
表示第j个节点处的仿射变换，
Figure BDA0002442730500000109
表示变换矩阵，
Figure BDA00024427305000001010
表示平移部分；N(pj)为第j个节点的邻域；
Figure BDA00024427305000001011
为旋转群，其中，
Figure BDA00024427305000001012
为一个单位阵；det(R)表示求解的行列式值；||·||F表示Frobenius范数，
Figure BDA00024427305000001013
表示投影到SO3的投影算子。
采用Welsch函数改进对齐项Ealign(X)以及正则项Ereg(X)得到
Figure BDA00024427305000001014
与
Figure BDA00024427305000001015
从而构造上述目标函数的上界函数，并作为替代函数，替代函数表示为：
Figure BDA00024427305000001016
本发明实施例中，采用拟牛顿方法迭代求解替代函数，得到非刚性变换矩阵：
Figure BDA00024427305000001017
其中，Xk+1表示第k+1步迭代获得的非刚性变换矩阵；计算步骤包括：
固定替代函数
Figure BDA00024427305000001018
作为目标函数，将第k步迭代获得的非刚性变换矩Xk作为拟牛顿求解的初值
Figure BDA00024427305000001019
并计算初始的近似Hessian矩阵H0，进而得到非刚性变换矩
Figure BDA00024427305000001020
不断的循环计算，在第t次循环计算过程中，得到迭代值
Figure BDA00024427305000001021
同时，根据第t-1次迭代值
Figure BDA00024427305000001022
计算梯度值
Figure BDA0002442730500000111
结合历史迭代值和历史梯度值得到第t次循环计算的近似Hessian矩阵Ht，并得到下降方向
Figure BDA0002442730500000112
再通过线性搜索得到最优步长λ，从而计算：
Figure BDA0002442730500000113
如果满足
Figure BDA0002442730500000114
终止循环计算，将得到的迭代值
Figure BDA0002442730500000115
作为k+1步迭代获得的非刚性变换矩阵Xk+1，t为循环计算的次数，∈1为设定的阈值。
本发明实施例中，对用于对齐项va和正则项参数vr采用动态调整的策略。首先令
Figure BDA0002442730500000116
固定va和vr达到收敛之后将令
Figure BDA0002442730500000117
vr＝0.5vr，返回步骤S102重新构建目标函数的替代函数，不断迭代直至
Figure BDA0002442730500000118
时迭代终止。其中，设置对齐项参数值νa上界为
Figure BDA0002442730500000119
Figure BDA00024427305000001110
为开始时初始模型中每个点到目标模型中的对应点距离的中值，其中目标模型中的对应点为距离初始模型中点的欧氏距离最近的点，νa下界通过初始模型的平均边长
Figure BDA00024427305000001111
决定，设为
Figure BDA00024427305000001112
正则项参数vr最大值设为
Figure BDA00024427305000001113
其中
Figure BDA00024427305000001114
为用户指定的参数。
上述系统中相关的计算部分在之前的方法实施例中进行了详细的说明，故不再赘述。
所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，仅以上述各功能模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能模块完成，即将系统的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM，U盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (8)

1.一种非刚性注册方法，其特征在于，包括：
步骤S101、基于输入的初始模型，根据给定的采样半径，得到嵌入简化图模型；
步骤S102、基于所述初始模型、嵌入简化图模型与输入的目标模型，进行最近邻搜索，得到匹配的对应点，构造原目标函数的替代函数；
步骤S103、采用拟牛顿方法迭代求解替代函数，得到非刚性变换矩阵；
步骤S104、判断所述非刚性变换矩阵是否收敛，若未收敛，则返回步骤S102更新替代函数；若收敛，则进入步骤S105；
步骤S105、判断是否达到参数值下界，若未达到参数值下界，将参数值减小后，返回步骤S102更新替代函数；若达到参数值下界，则转入步骤S106；
步骤S106、终止迭代，得到非刚性变换矩阵并将之应用到初始模型点上得到变换后的初始模型；
其中，通过给定的采样半径从所述初始模型中提取一部分的节点构成一个变形图模型，也即嵌入简化图模型，再通过嵌入简化图模型中节点的变换带动所述初始模型中所有点的变换；
根据给定的半径Ra在所述初始模型中迭代的提取采样点：对于初始模型中所有的点，将执行PCA提取协方差矩阵的最大特征值所在方向作为一个轴，将初始模型中所有的点投影到相应轴上并按投影值从小到大排序，根据排序添加采样点：如果采样点集合为空集或者初始模型中当前点到采样点集合中所有已存在点的测地线距离最小值大于给定半径Ra，则将当前点加入采样点集合，构成嵌入简化图的节点集；
连接测地线距离小于给定半径的任意两个节点构成嵌入简化图的边集，由嵌入简化图的节点集与边集构成嵌入简化图模型；
将嵌入简化图的点集中第j个节点处的仿射变换表示为(Aj,tj)；其中，
Figure FDA0003699000410000011
表示变换矩阵，
Figure FDA0003699000410000012
表示平移部分；如果初始模型中的第i个点在第j个节点的影响范围内，则用Aj(vi-pj)+pj+tj表示在变换(Aj,tj)作用下的新位置，
Figure FDA0003699000410000013
为第j个节点的三维坐标，
Figure FDA0003699000410000014
为第i个初始模型点的三维坐标，则第i个初始模型点在变形后的位置表示为：
Figure FDA0003699000410000015
其中，I(vi)表示影响第i个初始模型点的嵌入简化图节点集合，wij为权重；
所述基于所述初始模型、嵌入简化图模型与目标模型，进行最近邻搜索，得到匹配的对应点，构造原目标函数的替代函数包括：
基于所述初始模型、嵌入简化图模型与目标模型，通过对所述目标模型建立分割k维数据空间的数据结构进行最近邻搜索，得到匹配的对应点，并构造如下目标函数：
minEalign(X)+αEreg(X)+βErot(X)
其中，
Figure FDA0003699000410000021
表示对齐项，va为对齐项参数，n为初始模型中点的数量；
Figure FDA0003699000410000022
表示正则项，vr为正则项参数，Djl＝Aj(pl-pj)+pj+tj-(pl+tl)第j个节点与相邻节点l的变换差；
Figure FDA0003699000410000023
Figure FDA0003699000410000024
表示旋转项；α和β分别为正则项和旋转项的权重；
Figure FDA0003699000410000025
Figure FDA0003699000410000026
ν为指定参数，即νa和νr，x表示
Figure FDA0003699000410000027
与||Djl||；
Figure FDA0003699000410000028
分别为初始模型和目标模型，
Figure FDA0003699000410000029
分别为初始模型和目标模型上的点集，r为节点数，m为目标模型的点数，uρ(i)表示vi的对应点；
Figure FDA00036990004100000210
Figure FDA00036990004100000211
表示待求解的非刚性变换矩阵，
Figure FDA00036990004100000212
表示第j个节点处的仿射变换，
Figure FDA00036990004100000213
表示变换矩阵，
Figure FDA00036990004100000214
表示平移部分；N(pj)为第j个节点的邻域；
Figure FDA00036990004100000215
为旋转群，其中，
Figure FDA00036990004100000216
为一个单位阵，det(R)表示求解的行列式值；||·||F表示Frobenius范数，
Figure FDA00036990004100000217
表示投影到SO3的投影算子；
采用Welsch函数改进对齐项Ealign(X)以及正则项Ereg(X)得到
Figure FDA00036990004100000218
与
Figure FDA00036990004100000219
从而构造上述目标函数的上界函数，并作为替代函数，替代函数表示为：
Figure FDA00036990004100000220
2.根据权利要求1所述的一种非刚性注册方法，其特征在于，固定替代函数采用拟牛顿方法迭代求解替代函数，得到非刚性变换矩阵：
Figure FDA00036990004100000221
其中，Xk+1表示第k+1步迭代获得的非刚性变换矩阵；计算步骤包括：
固定替代函数
Figure FDA00036990004100000222
作为目标函数，将第k步迭代获得的非刚性变换矩Xk作为拟牛顿求解的初值
Figure FDA00036990004100000223
并计算初始的近似Hessian矩阵H0，进而得到非刚性变换矩
Figure FDA00036990004100000224
不断的循环计算，在第t次循环计算过程中，得到迭代值
Figure FDA00036990004100000225
同时，根据第t-1次迭代值
Figure FDA00036990004100000226
计算梯度值
Figure FDA0003699000410000031
结合历史迭代值和历史梯度值得到第t次循环计算的近似Hessian矩阵Ht，并得到下降方向
Figure FDA0003699000410000032
再通过线性搜索得到最优步长λ，从而计算：
Figure FDA0003699000410000033
如果满足
Figure FDA0003699000410000034
终止循环计算，将得到的迭代值
Figure FDA0003699000410000035
作为k+1步迭代获得的非刚性变换矩阵Xk+1，t为循环计算的次数，∈1为设定的阈值。
3.根据权利要求1所述的一种非刚性注册方法，其特征在于，对对齐项参数νa和正则项参数νr采用动态调整的策略：首先令
Figure FDA0003699000410000036
固定νa和νr达到收敛之后将令
Figure FDA0003699000410000037
νr＝0.5vr，返回步骤S102重新构建目标函数的替代函数，不断迭代直至
Figure FDA0003699000410000038
时迭代终止；其中，设置对齐项参数值va上界为
Figure FDA0003699000410000039
Figure FDA00036990004100000310
Figure FDA00036990004100000311
为开始时初始模型中每个点到目标模型中的对应点距离的中值，其中目标模型中的对应点为距离初始模型中点的欧氏距离最近的点，va下界通过初始模型的平均边长
Figure FDA00036990004100000312
决定，设为
Figure FDA00036990004100000313
正则项参数vr最大值设为
Figure FDA00036990004100000314
其中
Figure FDA00036990004100000315
为用户指定的参数。
4.一种非刚性注册系统，其特征在于，用于实现权利要求1～3任一项所述的方法，该系统包括：
采样与模型构建模块，用于执行步骤S101，即基于输入的初始模型，根据给定的采样半径，得到嵌入简化图模型；
函数构造模块，用于执行步骤S102，即基于所述初始模型、嵌入简化图模型与输入的目标模型，进行最近邻搜索，得到匹配的对应点，构造原目标函数的替代函数；
求解模块，用于执行步骤S103，即采用拟牛顿方法迭代求解替代函数，得到非刚性变换矩阵；
判断模块，用于执行步骤S104，即判断所述非刚性变换矩阵是否收敛，若未收敛，则返回步骤S102更新替代函数；若收敛，则进入步骤S105；
参数值判断与调整模块，用于执行步骤S105，即判断是否达到参数值下界，若未达到参数值下界，将参数值减小后，返回步骤S102更新替代函数；若达到参数值下界，则转入步骤S106；还用于执行步骤S106，即终止迭代，得到非刚性变换矩阵并将之应用到初始模型点上得到变换后的初始模型。
5.根据权利要求4所述的一种非刚性注册系统，其特征在于，通过给定的采样半径从所述初始模型中提取一部分的节点构成一个变形图模型，也即嵌入简化图模型，再通过嵌入简化图模型中节点的变换带动所述初始模型中所有点的变换；
根据给定的半径Ra在所述初始模型中迭代的提取采样点：对于初始模型中所有的点，将执行PCA提取协方差矩阵的最大特征值所在方向作为一个轴，将初始模型中所有的点投影到相应轴上并按投影值从小到大排序，根据排序添加采样点：如果采样点集合为空集或者初始模型中当前点到采样点集合中所有已存在点的测地线距离最小值大于给定半径Ra，则将当前点加入采样点集合，构成嵌入简化图的节点集；
连接测地线距离小于给定半径的任意两个节点构成嵌入简化图的边集，由嵌入简化图的节点集与边集构成嵌入简化图模型；
将嵌入简化图的点集中第j个节点处的仿射变换表示为(Aj,tj)；其中，
Figure FDA0003699000410000041
表示变换矩阵，
Figure FDA0003699000410000042
表示平移部分；如果初始模型中的第i个点在第j个节点的影响范围内，则用Aj(vi-pj)+pj+tj表示在变换(Aj,tj)作用下的新位置，
Figure FDA0003699000410000043
为第j个节点的三维坐标，
Figure FDA0003699000410000044
为第i个初始模型点的三维坐标，则第i个初始模型点在变形后的位置表示为：
Figure FDA0003699000410000045
其中，I(vi)表示影响第i个初始模型点的嵌入简化图节点集合，wij为权重。
6.根据权利要求4所述的一种非刚性注册系统，其特征在于，所述基于所述初始模型、嵌入简化图模型与目标模型，进行最近邻搜索，得到匹配的对应点，构造原目标函数的替代函数包括：
基于所述初始模型、嵌入简化图模型与目标模型，通过对所述目标模型建立分割k维数据空间的数据结构进行最近邻搜索，得到匹配的对应点，并构造如下目标函数：
minEalign(X)+αEreg(X)+βErot(X)
其中，
Figure FDA0003699000410000046
表示对齐项，va为对齐项参数，n为初始模型中点的数量；
Figure FDA0003699000410000047
表示正则项，vr为正则项参数，Djl＝Aj(pl-pj)+pj+tj-(pl+tl)第j个节点与相邻节点l的变换差；
Figure FDA0003699000410000048
Figure FDA0003699000410000049
表示旋转项；α和β分别为正则项和旋转项的权重；
Figure FDA00036990004100000410
Figure FDA00036990004100000411
v为指定参数，即va和vr，x表示
Figure FDA00036990004100000412
与||Djl||；
Figure FDA00036990004100000413
分别为初始模型和目标模型，
Figure FDA00036990004100000414
分别为初始模型和目标模型上的点集，r为节点数，m为目标模型的点数，uρ(i)表示vi的对应点；
Figure FDA00036990004100000415
Figure FDA0003699000410000051
表示待求解的节点处的非刚性变换矩阵，
Figure FDA0003699000410000052
表示第j个节点处的仿射变换，
Figure FDA0003699000410000053
表示变换矩阵，
Figure FDA0003699000410000054
表示平移部分；N(pj)为第j个节点的邻域；
Figure FDA0003699000410000055
为旋转群，其中，
Figure FDA0003699000410000056
为一个单位阵，det(R)表示求解的行列式值；||·||F表示Frobenius范数，
Figure FDA0003699000410000057
表示投影到SO3的投影算子；
采用Welsch函数改进对齐项Ealign(X)以及正则项Ereg(X)得到
Figure FDA0003699000410000058
与
Figure FDA0003699000410000059
从而构造上述目标函数的上界函数，并作为替代函数，替代函数表示为：
Figure FDA00036990004100000510
7.根据权利要求6所述的一种非刚性注册系统，其特征在于，采用拟牛顿方法迭代求解替代函数，得到非刚性变换矩阵：
Figure FDA00036990004100000511
其中，Xk+1表示第k+1步迭代获得的非刚性变换矩阵；计算步骤包括：
固定替代函数
Figure FDA00036990004100000512
作为目标函数，将第k步迭代获得的非刚性变换矩Xk作为拟牛顿求解的初值
Figure FDA00036990004100000513
并计算初始的近似Hessian矩阵H0，进而得到非刚性变换矩
Figure FDA00036990004100000514
不断的循环计算，在第t次循环计算过程中，得到迭代值
Figure FDA00036990004100000515
同时，根据第t-1次迭代值
Figure FDA00036990004100000516
计算梯度值
Figure FDA00036990004100000517
结合历史迭代值和历史梯度值得到第t次循环计算的近似Hessian矩阵Ht，并得到下降方向
Figure FDA00036990004100000518
再通过线性搜索得到最优步长λ，从而计算：
Figure FDA00036990004100000519
如果满足
Figure FDA00036990004100000520
终止循环计算，将得到的迭代值
Figure FDA00036990004100000521
作为k+1步迭代获得的非刚性变换矩阵Xk+1，t为循环计算的次数，∈1为设定的阈值。
8.根据权利要求6所述的一种非刚性注册系统，其特征在于，对对齐项va和正则项参数vr采用动态调整的策略：首先令
Figure FDA00036990004100000522
固定va和vr达到收敛之后将令
Figure FDA00036990004100000523
vr＝0.5vr，返回步骤S102重新构建目标函数的替代函数，不断迭代直至
Figure FDA00036990004100000524
时迭代终止；其中，设置对齐项参数值va上界为
Figure FDA00036990004100000525
Figure FDA00036990004100000526
为开始时初始模型中每个点到目标模型中的对应点距离的中值，其中目标模型中的对应点为距离初始模型中点的欧氏距离最近的点，va下界通过初始模型的平均边长
Figure FDA00036990004100000527
决定，设为
Figure FDA00036990004100000528
正则项参数vr最大值设为
Figure FDA00036990004100000529
其中
Figure FDA00036990004100000530
为用户指定的参数。
CN202010269838.6A 2020-04-08 2020-04-08 非刚性注册方法及系统 Active CN111462202B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010269838.6A CN111462202B (zh) 2020-04-08 2020-04-08 非刚性注册方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010269838.6A CN111462202B (zh) 2020-04-08 2020-04-08 非刚性注册方法及系统

Publications (2)

Publication Number Publication Date
CN111462202A CN111462202A (zh) 2020-07-28
CN111462202B true CN111462202B (zh) 2022-09-02

Family

ID=71680460

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010269838.6A Active CN111462202B (zh) 2020-04-08 2020-04-08 非刚性注册方法及系统

Country Status (1)

Country Link
CN (1) CN111462202B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103268631A (zh) * 2013-05-23 2013-08-28 中国科学院深圳先进技术研究院 点云骨架提取方法及装置
CN105631877A (zh) * 2015-12-28 2016-06-01 天津大学 基于加权双稀疏约束的非刚性表面配准方法
CN106934824A (zh) * 2017-03-16 2017-07-07 天津大学 可变形物体的全局非刚性配准与重建方法
CN108376408A (zh) * 2018-01-30 2018-08-07 清华大学深圳研究生院 一种基于曲率特征的三维点云数据快速加权配准方法
CN109202912A (zh) * 2018-11-15 2019-01-15 太原理工大学 一种基于单目深度传感器和机械臂配准目标轮廓点云的方法
CN110197503A (zh) * 2019-05-14 2019-09-03 北方夜视技术股份有限公司 基于增强型仿射变换的非刚性点集配准方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20150201910A1 (en) * 2014-01-17 2015-07-23 Centre For Imaging Technology Commercialization (Cimtec) 2d-3d rigid registration method to compensate for organ motion during an interventional procedure
US10360469B2 (en) * 2015-01-15 2019-07-23 Samsung Electronics Co., Ltd. Registration method and apparatus for 3D image data

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103268631A (zh) * 2013-05-23 2013-08-28 中国科学院深圳先进技术研究院 点云骨架提取方法及装置
CN105631877A (zh) * 2015-12-28 2016-06-01 天津大学 基于加权双稀疏约束的非刚性表面配准方法
CN106934824A (zh) * 2017-03-16 2017-07-07 天津大学 可变形物体的全局非刚性配准与重建方法
CN108376408A (zh) * 2018-01-30 2018-08-07 清华大学深圳研究生院 一种基于曲率特征的三维点云数据快速加权配准方法
CN109202912A (zh) * 2018-11-15 2019-01-15 太原理工大学 一种基于单目深度传感器和机械臂配准目标轮廓点云的方法
CN110197503A (zh) * 2019-05-14 2019-09-03 北方夜视技术股份有限公司 基于增强型仿射变换的非刚性点集配准方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Non-rigid point set registration via coherent spatial mapping;Jun Chen等;《Signal Processing》;20150131;第106卷;第62-72页 *
Robust single-view geometry and motion reconstruction;Hao Li等;《ACM Transactions on Graphics》;20091231;第28卷(第5期);第1-10页 *
一种面向大尺度变形的非刚性注册算法;李俊等;《计算机学报》;20110331;第34卷(第3期);第539-547页 *
三维几何模型的重建与结构优化;刘源;《信息科技辑》;20150915(第09期);第15-93页 *

Also Published As

Publication number Publication date
CN111462202A (zh) 2020-07-28

Similar Documents

Publication Publication Date Title
CN111986115A (zh) 激光点云噪声和冗余数据的精准剔除方法
CN106709947B (zh) 一种基于rgbd相机的三维人体快速建模系统
Li et al. Feature-preserving 3D mesh simplification for urban buildings
WO2015139574A1 (zh) 一种静态物体重建方法和系统
CN111795704A (zh) 一种视觉点云地图的构建方法、装置
CN111815686B (zh) 基于几何特征由粗到细点云配准方法
Li et al. Symmetry and template guided completion of damaged skulls
Zhu et al. AdaFit: Rethinking learning-based normal estimation on point clouds
CN113327275B (zh) 一种基于多约束点到局部曲面投影的点云双视角精配准方法
Das et al. Scan registration with multi-scale k-means normal distributions transform
CN105701455A (zh) 基于asm算法的人脸特征点采集及三维人脸建模方法
CN113034694B (zh) 自适应于流场特征的度量正交网格自动生成方法及装置
CN109034131B (zh) 一种半自动人脸关键点标注方法及存储介质
CN110009732A (zh) 基于gms特征匹配的面向复杂大尺度场景三维重建方法
CN110163902B (zh) 一种基于因子图的逆深度估计方法
Li et al. Three-dimensional point cloud registration based on normal vector angle
CN112184869A (zh) 基于绝对高斯曲率估计的保持几何特征的点云简化方法
CN110942077B (zh) 基于权重局部变化度和l1中值优化的特征线提取方法
CN112200915A (zh) 一种基于靶标三维模型纹理影像的前后形变量检测方法
CN110415339B (zh) 计算输入三维形体间的匹配关系的方法和装置
Yao et al. Fast and robust non-rigid registration using accelerated majorization-minimization
CN113313200B (zh) 一种基于法向约束的点云精匹配方法
Zang et al. Density-adaptive and geometry-aware registration of TLS point clouds based on coherent point drift
CN111462202B (zh) 非刚性注册方法及系统
CN116883590A (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