基于特征点的图像配准方法、系统和计算机可读存储介质

1.本发明属于辅助医疗技术领域，尤其涉及一种基于特征点的图像配准方法、系统和计算机可读存储介质。


背景技术：

2.近年来随着成像技术的迅速发展，医学图像也被广泛用于辅助医学手术、临床诊断、患者治疗等方面。传统意义上的二维平面影像技术如计算机x线摄影(dr)等，在医学上帮助的局限性日益显现。正电子发射断层成像(positron emission tomography，pet)、计算机断层成像(computed tomography,ct)等三维成像方法解决了二维图像提供病征信息少的局限性，这种成像方法可以立体而直观的展示患者体内的器官的三维信息，给医生提供充足而精确的病灶信息数据。但是这类成像技术辐射较高、操作复杂、成像时间长，对工作环境要求严苛，也使之在临床医治、手术上的使用受限。通常在脊柱手术中，采用的是x光图像这类二维图像，这种成像方法辐射较低也能较好的适应手术中的环境，最主要是它成像速度较快，在术中能够即拍即成，不会耽误手术进程。因此若能将术中的二维xr透视数据与术前的三维ct数据进行配准，就能结合两者长处，对手术导航提供精准且快速的帮助。
3.达到维度统一的方法不止一种。配准结果通过最大化图像梯度、灰度、轮廓特征等的相似性来变换图像对应点上的空间位置。目前的配准算法大致分为两类，一类是基于灰度，另一类是基于特征。
4.基于灰度的配准方法是目前应用最广泛也是研究最多的方法。它需要利用数字重建放射影像(digital reconstruction radiography,drr)技术，将术前的三维ct图像投影成二维drr图像，与术中的透射图像进行基于灰度的相似性度量。其中drr生成算法，参数搜索算法以及目标函数的选择会对配准的精度以及速度产生重大影响，因此大家的研究方向大致也是围绕这三个点。比如通过基于gpu的drr生成算法，跳出传统基于cpu的方法，大大提高了drr生成速度。通过slab算法，从算法角度提高了drr生成速度，但是这种方法的匹配误差较大，相似性测度性能比较差。在目标函数的选择上，有人提出了融合归一化互信息和梯度差分的相似性测度指标作为目标函数来适应不同的xr图像，提高鲁棒性。在参数搜索算法上，除了针对最优化算法的改进，还有例如为了进一步降低2d/3d配准过程中的大量计算，有人提出空间参数解耦的方法来配准。使用梯度方向加权直方图提取图像的信息，计算直方图之间的相似性从而求得旋转参数和平移参数，完成配准。
5.基于特征的配准方法分为基于标记物和内部固有特征两类，标记物指的是利用打入患者体内的钛钉等作为配准的定位点，但是这类方法会对患者身体造成额外的伤害。基于内部固有特征的方法指的是提取目标组织已存在的特征点，例如轮廓，角点等，据此进行匹配，与上述提到的在目标函数中加入梯度差分项有异曲同工之妙，用sobel和laplas算子提取骨骼的边缘特征进行配准，这种方法利用了图像除灰度外的更多信息，对精度的提升有显著效果，但是他们的实验数据只有模型数据，所以图像的梯度信息与骨骼结构相吻合，但是真实人体数据中有很多其他器官的干扰，会带来很多干扰项，严重的会将骨骼的梯度
信息淹没在噪声里。还有很多研究使用的数据是头骨数据，没有其他器官的干扰，噪声更小。
6.值得一提的是近几年深度学习的迅速发展，使得很多研究者将神经网络应用到配准任务上来，特别在特征识别方面，卷积神经网络有着很高的精度及极快的速度。但是一方面，目前暂时没有公开的利用深度学习实现2d/3d刚性配准的网络，不同的网络结构会极大的影响配准效果，开源网络大多是端到端的非刚性配准例如voxelmorph，主要针对血管等组织的变化进行非刚性配准，但对于椎骨来说，骨性结构因为病人体位的变化导致的位移是刚性变化，现有的非刚性配准方法显然不可行。有人采用提取多张2d图的轮廓特征组合成3d轮廓，再提取ct的轮廓，配准两者的轮廓点云的方法，然而这种方法需要基于精准的轮廓提取方法，且正常情况难以即时获得多张2d图像。
7.另一方面，因为医疗数据的特殊性，训练数据极度匮乏，导致深度学习不是精度不如传统方法，就是泛化性很差，但是对于临床医学来说，精度的不足使得深度学习引以为傲的速度变的毫无意义。
8.有鉴于此，有必要设计一种新的配准方法，以解决上述问题。


技术实现要素：

9.本发明的目的是解决现有技术中对于特征点较多的器官，如何保证在提升精度的同时加快配准速度的问题。
10.为实现以上目的，本发明提供了一种基于特征点的图像配准方法，包括以下步骤：
11.步骤s1、将ct数据经数字重建放射影像技术生成drr图像；
12.步骤s2、将x射线透视图像与drr图像经过图像增强后输入特征匹配神经网络，得到粗配准的位置参数；
13.步骤s3、将粗配准得到的位置参数作为精配准的初始解，利用粒子群算法结合powell算法得到六个空间变换参数。
14.本发明的进一步改进在于，在步骤s3中，先使用随机优化算法粒子群算法将参数收敛至全局最优点附近，再使用局部最优powell算法将参数收敛至最优点。
15.本发明的进一步改进在于，步骤s3中粒子群算法包括以下步骤：
16.步骤s311、设置最大迭代次数、目标函数的自变量个数、粒子的最大速度，位置信息为整个搜索空间，在速度区间和搜索空间上随机初始化速度和位置，设置粒子群规模为m，每个粒子随机初始化一个速度；
17.步骤s312、定义适应度函数，个体极值为每个粒子找到的最优解，从这些最优解中找到一个全局值，该全局值为本次全局最优解；
18.步骤s313、与历史全局最优解进行比较，并更新，更新速度和位置的公式分别为：
[0019][0020]
其中，w为根据迭代次数线性下降的自适应权重，wmax为最大权重，wmin为最小权重，iter为当前迭代次数，ger为最大迭代次数；
[0021]
v＝w
×
v+c1
×
rand
×
(gbest-x)+c2
×
rand
×
(zbest-x)
[0022]
其中，v为更新的速度，rand为0～1之间的随机数，gbest为局部最优位置，zbest为
全局最优解，c1、c2为速度因子，分别控制向局部以及全局最优解的靠近速度；
[0023]
x＝x+v
[0024]
其中，x为更新的位置也就是输入参数，等于更新的速度与之前位置之和。
[0025]
本发明的进一步改进在于，步骤s3中powell算法的初始点x0由粒子群算法提供，该powell算法包括以下步骤：
[0026]
步骤s321、选取n个线性无关的初始搜索方向d0,d1,...,dn-1；给定允许误差err》0,令k＝0；
[0027]
步骤s322、进行基本搜索：令y0＝x(k),依次沿d0,d1,...,dn-1进行一维搜索；对一切j＝1,2,...,n，记f(y(j-1)+λ(j-1)d(j-1))＝min(f(y(j-1)+λd(j-1)))，y(j)＝f(y(j-1)+λ(j-1)d(j-1))；
[0028]
步骤s323、检查是否满足终止条件：取加速度方向d(n)＝y(n)-y(0)；若||d(n)||《err,则迭代终止，得yn为问题的最优解；否则转步骤s324；
[0029]
步骤s324、确定搜索方向，按照前述公式确定m，若验证式子成立，转步骤s325；否则，转步骤s326；
[0030]
步骤s325、调整搜索方向，从点yn出发沿方向dn进行一维搜索，求出λn,使得f(yn+λn*dn)＝min f(yn+λ*dn)；令x(k+1)＝yn+λn*dn；再令d(j)＝d(j+1)，j＝m,m+1,...,n-1.k＝k+1，返回步骤s322；
[0031]
步骤s326、不调整搜索方向，令x(k+1)＝yn；k＝k+1,返回步骤s322。
[0032]
本发明的进一步改进在于，在步骤s2中，特征提取网络的输入是xr和drr的二维图像，输出是两张图像中的特征点来计算对应位置关系；特征点检测网络为全卷积网络(fcnn)，且对训练数据进行旋转、缩放等数据增强；特征提取网络的编码部分采用vgg结构降低图像的维度后分为两部分：特征点检测和描述子检测；其中，特征点检测的特征图有65个通道，对应8*8个区域通道和一个无关键点的通道，通过softmax去除最后一个通道，reshape得到最后的h*w的特征图；描述子检测采用了ucn的model和双三次插值和l2归一化，并且一个关键点对应d维特征。
[0033]
本发明的进一步改进在于，在步骤s2中，粗配准选择在z轴上间隔1
°
生成drr，计算与xr图像的特征点相似度后得到z轴旋转的初始解rz，其他参数设为默认值，即r
x
＝-90
°
,ry＝0
°
,t
x
＝ty＝tz＝0。
[0034]
本发明的进一步改进在于，为了综合分析量化数据，定义配准误差为:
[0035][0036]
其中，rj,ri为配准后的参数，ri为需要计算误差的参数值，rj为其他参数值。为各切片相对于第一个切片的旋转参数真值，第一个切片的旋转参数为0
°
。
[0037]
本发明的进一步改进在于，将特征点的相似性度量sim_points与ncc结合，得到新的相似性度量np：
[0038]
np＝w1
×
0.5
×
(1-ncc)+w2
×k×
sim_points
[0039]
其中，w1和w2均为0.5，k为统一两者数量级的系数，取10000；
[0040]
其中，
[0041]
sim_points＝(cosine_distance([x
11
,x
12
…
x
1i
],[x
21
,x
22
…
x
2j
])
[0042]
+cosine_distance([y
21
,y
22
…y2i
],[y
21
,y
22
…y2j
]))/2
[0043]
其中，
[0044][0045]
其中，
[0046][0047]
x,y为特征点的坐标值。
[0048]
为了实现以上发明目的，本发明还提供了一种系统，执行前述任一项所述的方法。
[0049]
为了实现以上发明目的，本发明还提供了一种计算机可读存储介质，介质中存储有计算机程序，所述计算机程序被处理器进行加载，以执行前述任一项所述的方法中的步骤。
[0050]
本发明的有益效果如下：本发明提出了基于特征点的配准方法、系统以及计算机存储介质，利用一个特征提取网络识别图像的点特征而非物体的形状特征，对高噪声的图像具有更高的鲁棒性；采用自监督和半监督的训练模式，利用大量自然图像训练通用模型解决数据量问题，并添加了图像增强模块，以适应不同模态的图像质量；利用改进的粒子群算法结合powell算法进行精配准确定六个空间变换参数，在保证精度的同时，利用深度学习提高了配准速度；最后，同时利用灰度与脊柱特征信息，提出添加了特征点相似度项的目标函数，以此控制下降方向，提升精度并加快收敛速度。
附图说明
[0051]
图1是本发明提供的基于特征点的图像配准方法的流程图。
[0052]
图2是本发明提供的特征匹配网络配准流程示意图。
[0053]
图3是本发明提供的特征提取网络的网络结构示意图。
[0054]
图4是匹配图像对特征点匹配结果的效果图，该不匹配图的特征连线稀疏且交错。
[0055]
图5是匹配图像对特征点匹配结果的效果图，该匹配图特征点连线丰富且基本平行。
[0056]
图6是表1中第一组的四类收敛曲线图：(a)为有初始解，np；(b)为无初始解，np；(c)为有初始解，ncc；(d)为无初始解，ncc。
[0057]
图7是表1中第四组的配准效果演示图：左图为xr图像，右图为drr图像。
具体实施方式
[0058]
为了使本发明的目的、技术方案和优点更加清楚，下面结合附图和具体实施例对本发明进行详细描述。
[0059]
需要强调的是，在描述本发明过程中，各种公式和约束条件分别使用前后一致的标号进行区分，但也不排除使用不同的标号标志相同的公式和/或约束条件，这样设置的目的是为了更清楚的说明本发明特征所在。
[0060]
请参考图1、图2所示，本发明的配准方法主要包括图像预处理步骤、粗配准步骤以及精配准步骤三部分。
[0061]
步骤s1、将ct数据经数字重建放射影像技术生成drr图像。
[0062]
步骤s2、粗配准：将x射线透视图像与drr图像经过图像增强后输入特征匹配神经网络，得到粗配准的结果。
[0063]
步骤s3、精配准：将粗配准得到的位置参数作为精配准的初始解，利用优化的粒子群算法得到最终的六个空间变换参数。
[0064]
以下将对此进行详细说明。
[0065]
1、最优化方法
[0066]
基于投影策略的2d/3d医学图像配准需要求解六个空间变换的旋转平移参数，而对于没有标记物的图像对，无法构建对应方程，也就没有解析解，因此只能通过迭代搜索的方式求最优解。
[0067]
不同的最优化算法的寻优能力是不同的，目前被广泛应用的局部最优化算法包括：下降单纯形法(downhill simplex)、鲍威尔(powell)算法、levenberg-marquardt搜索法、拟牛顿法等，另外还有针对寻求全局最优提出来的随机优化算法，例如，例如遗传算法(genetic algorithm，ga)、模拟退火算法(simulate anneal，sa)、粒子群优化算法(particle swarm optimization，pso)等。
[0068]
随机优化算法容易跳出局部最优找到全局最优解，但是只能收敛至最优解附近，而局部最优化算法很容易陷入局部最优却能很快收敛至局部最优点。脊柱配准有六个维度的变量需要求解，很容易陷入局部最优，因此我们选择结合两种算法，先使用随机优化算法粒子群算法将参数收敛至全局最优点附近，再使用局部最优算法powell算法将参数收敛至最优点。
[0069]
1.1优化方法的初始解求解
[0070]
最优化方法是在解空间中按损失函数的下降方向迭代搜索全局最优解，因此初始点的选取会极大的影响求解的速度与精度，精确的初始点能避免算法陷入局部最优，更能加速搜索至最优点的速度，一般算法的初始点为全0，然后依靠多分辨率策略确定初始点，过于费时。因为z轴的旋转数据变化范围为360
°
，其他参数最大只有几十的变化，所以我们优先考虑z轴的初始解，我们将ct数据在z轴上以1
°
的步长生成360张drr图像，依靠特征提取网络，计算特征点间的余弦距离，选取最小距离的角度作为初始解。精度可达到2
°
。在此基础上进行最优化搜索，网络的具体结构图3。
[0071]
下面介绍这两种算法以及针对本发明的优化。
[0072]
1.2粒子群算法
[0073]
粒子群算法的基本概念源于对鸟群觅食行为的研究。用一种粒子来模拟上述的鸟类个体，每个粒子可视为n维搜索空间中的一个搜索个体，粒子的当前位置即为对应优化问题的一个候选解，粒子的飞行过程即为该个体的搜索过程.粒子的飞行速度可根据粒子历史最优位置和种群历史最优位置进行动态调整.粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子单独搜寻的最优解叫做个体极值，粒子群中最优的个体极值作为当前全局最优解。不断迭代，更新速度和位置。最终得到满足终止条件的最优解。
[0074]
算法流程如下：
[0075]
step1(步骤s311)：初始化
[0076]
首先，我们设置最大迭代次数，目标函数的自变量个数，粒子的最大速度，位置信息为整个搜索空间，我们在速度区间和搜索空间上随机初始化速度和位置，设置粒子群规模为m，每个粒子随机初始化一个速度。
[0077]
step2(步骤s312)：个体极值与全局最优解
[0078]
定义适应度函数，个体极值为每个粒子找到的最优解，从这些最优解找到一个全局值，叫做本次全局最优解。与历史全局最优比较，进行更新。
[0079]
step3(步骤s313)：更新速度和位置的公式
[0080][0081]
w为根据迭代次数线性下降的自适应权重，wmax为最大权重，wmin为最小权重，iter为当前迭代次数，ger为最大迭代次数。
[0082]
v＝w
×
v+c1
×
rand
×
(gbest-x)+c2
×
rand
×
(zbest-x)
[0083]
v为更新的速度，rand为0～1之间的随机数，gbest为局部最优位置，zbest为全局最优解，c1，c2为速度因子，分别控制向局部以及全局最优解的靠近速度，默认均为经验值2，这样更容易保持收敛速度和搜索效果的均衡。
[0084]
x＝x+v
[0085]
x为更新的位置也就是输入参数，等于更新的速度与之前位置之和。
[0086]
1.3 powell算法
[0087]
鲍威尔法又称方向加速法，它由powell于1964年提出，是利用共轭方向可以加快收敛速度的性质形成的一种搜索方法。该方法不需要对目标函数进行求导，当目标函数的导数不连续的时候也能应用，因此，鲍威尔算法是一种十分有效的直接搜索法。
[0088]
powell法可用于求解一般无约束优化问题，对于维数n《20的目标函数求优化问题，此法可获得较满意的结果。不同于其他的直接法，powell法有一套完整的理论体系，故其计算效率高于其他直接法。该方法使用一维搜索，而不是跳跃的探测步。然而，在原始的powell算法中，必须保持每次迭代中前n搜索方向线性无关，否则将永远得不到问题的最优解。当某一循环方向组中的矢量系出现线性相关的情况(退化，病态)时，搜索过程在降维的空间进行，致使计算不能收敛而失败。为了避免此种情况产生，我们使用修正的powell算法
[0089]
算法流程如下：
[0090]
初始点x0由粒子群算法提供。
[0091]
step1(步骤s321)：选取n个线性无关的初始搜索方向d0,d1,...,dn-1；给定允许误差err》0,令k＝0；
[0092]
step2(步骤s322)：进行基本搜索，令y0＝x(k),依次沿d0,d1,...,dn-1进行一维搜索。对一切j＝1,2,...,n，记f(y(j-1)+λ(j-1)d(j-1))＝min(f(y(j-1)+λd(j-1)))，y(j)＝f(y(j-1)+λ(j-1)d(j-1))；
[0093]
step3(步骤s323)：检查是否满足终止条件，取加速度方向d(n)＝y(n)-y(0)；若||d(n)||《err,迭代终止，得yn为问题的最优解；否则转step4。
[0094]
step4(步骤s324)：确定搜索方向，按照上式公式确定m，若验证式子成立，转
step5；否则，转step6。
[0095]
step5(步骤s325)：调整搜索方向，从点yn出发沿方向dn进行一维搜索，求出λn,使得f(yn+λn*dn)＝min f(yn+λ*dn)；令x(k+1)＝yn+λn*dn。再令d(j)＝d(j+1)，j＝m,m+1,...,n-1.k＝k+1,返回step2。
[0096]
step6(步骤s326)：不调整搜索方向，令x(k+1)＝yn；k＝k+1,返回step2。
[0097]
2、基于深度学习的特征提取网络
[0098]
深度学习之于传统算法的优点是耗时短，速度快，实际应用时只需要权重与输入数据间进行简单的运算，但是在精度上却不如传统算法那样精准可控，且需要大量训练数据，这在医学相关算法的研究中很难实现，所以本发明将深度学习应用于粗配准，提升整体的配准速度，在精配准中依旧使用以优化算法为基础的传统算法保证配准精度。为了解决数据量不足的问题，网络设计采用自监督和半监督形式训练通用模型，识别物体的点特征而不是形状特征，使得其在医学图像上也能有较好的识别精度。
[0099]
脊柱骨节和棘突部分的特征点较为密集，本发明计划构建的特征提取网络的输入是xr和drr的二维图像，输出是两张图像中的特征点来计算对应位置关系。特征点检测网络为全卷积网络(fcnn)，且对训练数据进行旋转，缩放等数据增强，因此特征点的检测具有旋转和平移不变性，切合脊柱在空间上的平移旋转变换的实际情形，适合应用于配准任务。
[0100]
如图3所示，特征提取网络包括编码部分和解码部分。
[0101]
编码部分采用vgg结构去降低图像的维度，然后分为两部分，一部分为特征点检测，另一部分为描述子检测。
[0102]
第一部分的特征图有65个通道，对应8*8个区域通道和一个无关键点的通道。通过softmax去除最后一个通道，reshape得到最后的h*w的特征图。
[0103]
第二部分采用了类似ucn的model和双三次插值和l2归一化。一个关键点对应d维特征。
[0104]
训练的数据采用ms-coco数据集，首先标注一部分训练数据作为预训练数据，然后以预训练模型预测其他数据的真值，自监督地训练模型，因为特征点多是拐点或者角点，所以对于一副图片来说并没有一个唯一准确的真值，因此自监督的真值质量与手工标注的质量接近一致。
[0105]
同时，不同模态不同角度的数据投影出的图片质量参差不齐，本发明在检测特征点之前增加图像增强模块，自动调整图像对的对比度，饱和度等，利于提高后续步骤的精度。
[0106]
3、相似性度量
[0107]
图像相似性的度量方法有ncc(normalized cross correlation)，归一化互相关算法，ssim((structural similarity))结构相似性指标等，ncc的相似性测度基于整幅图像的灰度，对于图像对i1，i2，ncc的计算公式如下：
[0108][0109]
ncc只针对图像的灰度做相似性度量，为了更大化的利用图像信息，脊柱的特征可
以同时作为相似性度量的衡量标准，考虑对于真实人体数据，其他器官的影响不利于脊柱梯度的提取，而特征匹配网络的特征提取结果对其他噪声的影响更加鲁棒，图像对的特征点相似性度量采用所有特征点坐标组成的向量的余弦距离。
[0110]
sim_points＝(cosine_distance([x
11
,x
12
…
x
1i
],[x
21
,x
22
…
x
2j
])
[0111]
+cosine_distance([y
21
,y
22
…y2i
],[y
21
,y
22
…y2j
]))/2
[0112]
其中，x,y为特征点的坐标值。
[0113]
本发明将特征点的相似性度量sim_points与ncc结合，得到新的相似性度量np：
[0114]
np＝w1
×
0.5
×
(1-ncc)+w2
×k×
sim_points
[0115]
本发明中w1和w2均为0.5，k为统一两者数量级的系数，取10000。
[0116]
本发明实验所用数据包括真实人体数据以及脊柱模型数据，其中ct数据大小为512*512*512，xr数据为n*1024*1024(n为张数)。
[0117]
初始解对于迭代搜索算法的精度与速度的提升有重要的作用，精确的初始解可以避免搜索陷入局部最优，同时缩短到达最优点的距离。在求初始解时，一般采用多分辨率策略，将低分辨率的图像配准结果作为初始解，此种方法依旧需要采用迭代搜索，考虑到数据集中z轴的旋转角度最大，需要考虑360
°
，x和y轴的旋转角度都相对较小，且特征识别网络对平移具有鲁棒性，因此粗配准选择在z轴上间隔1
°
生成drr，计算与xr图像的特征点相似度后得到z轴旋转的初始解rz，其他参数设为默认值，即r
x
＝-90
°
,ry＝0
°
,t
x
＝ty＝tz＝0。利用特征匹配网络计算xr和drr图像的特征点的匹配度，缩小drr配准集的范围，在粗配准阶段的耗时从低分辨率迭代搜索需要的90+s缩短至10+s，节省了88.9％的时间。
[0118]
在真实人体数据上本文实验了初始解对配准精度和速度的影响，从c臂上采集的ct数据只有各切片旋转参数相对值的真值，因此为了综合分析量化数据，我们定义配准误差为:
[0119][0120]
其中，rj,ri为配准后的参数，ri为需要计算误差的参数值，rj为其他参数值；为各切片相对于第一个切片的旋转参数真值(第一个切片的旋转参数为0
°
)。
[0121]
表1为真实人体数据各条件下配准结果(20次迭代)。图6为表1中第一组的四类收敛曲线：(a)为有初始解，np，(b)为无初始解，np，(c)为有初始解，ncc，(d)为无初始解，ncc。
[0122]
序号有初始解，np/e(
°
)无初始解，np/e(
°
)有初始解，ncc/e(
°
)无初始解，ncc/e(
°
)10.1785511192.7498407893.54122565435.3813584220.5312382712.8232172532.23229624710.5454285230.1945781927.4559841168.79786756544.1008851240.1785511193.4861541272.23229624719.3124691650.212387722.4719189812.40129284819.3124691660.2347114472.4719189812.3720808220.55129097平均0.2550029783.5765057083.59617656324.86731689
[0123]
表1真实人体数据各条件下配准结果(20次迭代)
[0124]
表1均匀取样了6个xr图像进行实验，从第二三列以及第四五列可以看出在20次迭
代的情况下，有初始解的配准的精度远高于无初始解的配准，当然，理论上在迭代次数足够的情况下，无初始解也能搜索至最优点，但是会耗费数倍的时间。从第一三列以及二四列可以看到，在20个迭代的情况下，np组的精度远高于ncc组，没有特征点辅助的情况下，ncc的第三组结果出现了较大偏差，进入了局部最优点未能跳出。所以特征匹配网络在粗配准和目标函数的引入能有效提高配准的精度。并且，从图6可以看出相似性度量np对收敛速度的影响比初始解更大，基本在第10个迭代就大致完成了收敛。表1中第四组配准效果演示如图7所示，左图为xr图像，右图为drr图像。
[0125]
其他研究者提取图像的梯度信息作为特征信息参与配准，例如用sobel算子，log(laplas+高斯滤波)算子等，但是他们都是在模型数据上开展的实验，或是加上了高斯噪声简单模拟实际场景等，并没有试验真实数据的效果。
[0126]
表2和表3对本文算法在真实数据上的效果做了展示。表2可看到针对两个梯度提取算法在存在大量噪声的真实人体数据上，无法良好地提取出需要的特征的问题，特征匹配网络在两种数据中都能稳定地提取到特征点并进行匹配，具有良好的鲁棒性。
[0127][0128]
表2真实人体数据和模型数据上的各算法效果对比
[0129]
自表3可看到在不同数据下，另外两种梯度提取算法在模型数据上的效果好于人体数据，sobel的配准精度更高于log。我们的方法在两种数据上的配准精度都是最高的，甚至在人体数据上的精度更高，非骨骼的特征点也能被提取用于辅助配准。说明我们的方法在不同的噪声水平下都能有稳定的高精度。
[0130][0131]
表3不同数据上不同方法的配准误差
[0132]
为了更好地实施本技术实施例中基于特征点的图像配准方法，在该方法的基础上，本技术实施例还提供了一种可实施以上基于特征点的图像配准方法的系统。该系统既可以为一个终端，也可以展现为一个配准装置。该终端或该配准装置包括一个或多个处理器、存储器以及一个或多个应用程序。一个或多个应用程序被存储于存储器中，配置为由处理器执行上述基于特征点的图像配准方法的步骤。
[0133]
综上所述，本发明提出了一套2d/3d配准方法、系统。利用一个特征提取网络识别图像的点特征而非物体的形状特征，所以对高噪声的图像具有更高的鲁棒性，采用自监督和半监督的训练模式，利用大量自然图像训练通用模型解决数据量问题，并添加了图像增强模块，以适应不同模态的图像质量。如此网络无需针对性的训练就能继续用于其他器官的2d-3d配准工作，特别是特征点较多的器官，例如其他部位的骨骼等。再利用改进的粒子群算法结合powell算法进行精配准确定六个空间变换参数，在保证精度的同时，利用深度学习提高了配准速度。同时利用灰度与脊柱特征信息，提出添加了特征点相似度项的目标函数，以此控制下降方向，提升精度并加快收敛速度。
[0134]
同时，本发明实施例还提供了一种计算机可读存储介质，该计算机可读存储介质可以包括：只读存储器(rom,read only memory)、随机存取记忆体(ram,random access memory)、磁盘或光盘等。该计算机可读存储介质中存储有多条指令，该指令能够被处理器进行加载，以执行本发明所提供的的、基于特征点的图像配准方法的步骤。例如，该指令可以执行如下步骤：
[0135]
步骤s1、将ct数据经数字重建放射影像技术生成drr图像；
[0136]
步骤s2、将x射线透视图像与drr图像经过图像增强后输入特征匹配神经网络，得到粗配准的位置参数；
[0137]
步骤s3、将粗配准得到的位置参数作为精配准的初始解，利用粒子群算法结合powell算法得到六个空间变换参数。
[0138]
或者，在步骤s3中，执行以下粒子群算法步骤：
[0139]
步骤s311、设置最大迭代次数、目标函数的自变量个数、粒子的最大速度，位置信息为整个搜索空间，在速度区间和搜索空间上随机初始化速度和位置，设置粒子群规模为
m，每个粒子随机初始化一个速度；
[0140]
步骤s312、定义适应度函数，个体极值为每个粒子找到的最优解，从这些最优解中找到一个全局值，该全局值为本次全局最优解；
[0141]
步骤s313、与历史全局最优解进行比较，并更新，更新速度和位置的公式分别为：
[0142][0143]
其中，w为根据迭代次数线性下降的自适应权重，wmax为最大权重，wmin为最小权重，iter为当前迭代次数，ger为最大迭代次数；v＝w
×
v+c1
×
rand
×
(gbest-x)+c2
×
rand
×
(zbest-x)
[0144]
其中，v为更新的速度，rand为0～1之间的随机数，gbest为局部最优位置，zbest为全局最优解，c1、c2为速度因子，分别控制向局部以及全局最优解的靠近速度；
[0145]
x＝x+v
[0146]
其中，x为更新的位置也就是输入参数，等于更新的速度与之前位置之和。
[0147]
以及/或者，在步骤s3中，执行以下步骤：
[0148]
步骤s321、选取n个线性无关的初始搜索方向d0,d1,...,dn-1；给定允许误差err》0,令k＝0；
[0149]
步骤s322、进行基本搜索：令y0＝x(k),依次沿d0,d1,...,dn-1进行一维搜索；对一切j＝1,2,...,n，记f(y(j-1)+λ(j-1)d(j-1))＝min(f(y(j-1)+λd(j-1)))，y(j)＝f(y(j-1)+λ(j-1)d(j-1))；
[0150]
步骤s323、检查是否满足终止条件：取加速度方向d(n)＝y(n)-y(0)；若||d(n)||《err,则迭代终止，得yn为问题的最优解；否则转步骤s324；
[0151]
步骤s324、确定搜索方向，按照前述公式确定m，若验证式子成立，转步骤s325；否则，转步骤s326；
[0152]
步骤s325、调整搜索方向，从点yn出发沿方向dn进行一维搜索，求出λn,使得f(yn+λn*dn)＝min f(yn+λ*dn)；令x(k+1)＝yn+λn*dn；再令d(j)＝d(j+1)，j＝m,m+1,...,n-1.k＝k+1，返回步骤s322；
[0153]
步骤s326、不调整搜索方向，令x(k+1)＝yn；k＝k+1,返回步骤s322。
[0154]
在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述的部分，可以参见其他实施例的相关描述。
[0155]
以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。

技术特征：
1.一种基于特征点的图像配准方法，其特征在于，包括以下步骤：步骤s1、将ct数据经数字重建放射影像技术生成drr图像；步骤s2、将x射线透视图像与drr图像经过图像增强后输入特征匹配神经网络，得到粗配准的位置参数；步骤s3、将粗配准得到的位置参数作为精配准的初始解，利用粒子群算法结合powell算法得到六个空间变换参数。2.根据权利要求1所述的方法，其特征在于：在步骤s3中，先使用随机优化算法粒子群算法将参数收敛至全局最优点附近，再使用局部最优powell算法将参数收敛至最优点。3.根据权利要求2所述的方法，其特征在于：步骤s3中粒子群算法包括以下步骤：步骤s311、设置最大迭代次数、目标函数的自变量个数、粒子的最大速度，位置信息为整个搜索空间，在速度区间和搜索空间上随机初始化速度和位置，设置粒子群规模为m，每个粒子随机初始化一个速度；步骤s312、定义适应度函数，个体极值为每个粒子找到的最优解，从这些最优解中找到一个全局值，该全局值为本次全局最优解；步骤s313、与历史全局最优解进行比较，并更新，更新速度和位置的公式分别为：其中，w为根据迭代次数线性下降的自适应权重，wmax为最大权重，wmin为最小权重，iter为当前迭代次数，ger为最大迭代次数；v＝w
×
v+c1
×
rand
×
(gbest-x)+c2
×
rand
×
(zbest-x)其中，v为更新的速度，rand为0～1之间的随机数，gbest为局部最优位置，zbest为全局最优解，c1、c2为速度因子，分别控制向局部以及全局最优解的靠近速度；x＝x+v其中，x为更新的位置也就是输入参数，等于更新的速度与之前位置之和。4.根据权利要求3所述的方法，其特征在于：步骤s3中powell算法的初始点x0由粒子群算法提供，该powell算法包括以下步骤：步骤s321、选取n个线性无关的初始搜索方向d0,d1,...,dn-1；给定允许误差err＞0,令k＝0；步骤s322、进行基本搜索：令y0＝x(k),依次沿d0,d1,...,dn-1进行一维搜索；对一切j＝1,2,...,n，记f(y(j-1)+λ(j-1)d(j-1))＝min(f(y(j-1)+λd(j-1)))，y(j)＝f(y(j-1)+λ(j-1)d(j-1))；步骤s323、检查是否满足终止条件：取加速度方向d(n)＝y(n)-y(0)；若||d(n)||＜err,则迭代终止，得yn为问题的最优解；否则转步骤s324；步骤s324、确定搜索方向，按照前述公式确定m，若验证式子成立，转步骤s325；否则，转步骤s326；步骤s325、调整搜索方向，从点yn出发沿方向dn进行一维搜索，求出λn,使得f(yn+λn*dn)＝min f(yn+λ*dn)；令x(k+1)＝yn+λn*dn；再令d(j)＝d(j+1)，j＝m,m+1,...,n-1.k＝k+1，返回步骤s322；
步骤s326、不调整搜索方向，令x(k+1)＝yn；k＝k+1,返回步骤s322。5.根据权利要求1所述的方法，其特征在于：在步骤s2中，特征提取网络的输入是xr和drr的二维图像，输出是两张图像中的特征点来计算对应位置关系；特征点检测网络为全卷积网络(fcnn)，且对训练数据进行旋转、缩放等数据增强；特征提取网络的编码部分采用vgg结构降低图像的维度后分为两部分：特征点检测和描述子检测；其中，特征点检测的特征图有65个通道，对应8*8个区域通道和一个无关键点的通道，通过softmax去除最后一个通道，reshape得到最后的h*w的特征图；描述子检测采用了ucn的model和双三次插值和l2归一化，并且一个关键点对应d维特征。6.根据权利要求5所述的方法，其特征在于：在步骤s2中，粗配准选择在z轴上间隔1
°
生成drr，计算与xr图像的特征点相似度后得到z轴旋转的初始解r
z
，其他参数设为默认值，即r
x
＝-90
°
,r
y
＝0
°
,t
x
＝t
y
＝t
z
＝0。7.根据权利要求6所述的方法，其特征在于：为了综合分析量化数据，定义配准误差为:其中，r
j
,r
i
为配准后的参数，r
i
为需要计算误差的参数值，r
j
为其他参数值。为各切片相对于第一个切片的旋转参数真值，第一个切片的旋转参数为0
°
。8.根据权利要求7所述的方法，其特征在于：将特征点的相似性度量sim_points与ncc结合，得到新的相似性度量np：np＝w1
×
0.5
×
(1-ncc)+w2
×
k
×
sim_points其中，w1和w2均为0.5，k为统一两者数量级的系数，取10000；其中，sim_points＝(cosine_distance([x
11
,x
12
…
x
1i
],[x
21
,x
22
…
x
2j
])+cosine_distance([y
21
,y
22
…
y
2i
],[y
21
,y
22
…
y
2j
]))/2其中，其中，x,y为特征点的坐标值。9.一种系统，其特征在于：执行权利要求1至8任一项所述的方法。10.一种计算机可读存储介质，其特征在于：存储有计算机程序，所述计算机程序被处理器进行加载，以执行权利要求1～8任一项所述的方法中的步骤。

技术总结
本发明提供了一种基于特征点的图像配准方法、实施该方法的系统以及计算机可读存储介质。所述方法包括以下步骤：将CT数据经数字重建放射影像技术生成DRR图像；将X射线透视图像与DRR图像经过图像增强后输入特征匹配神经网络，得到粗配准的位置参数；以及将粗配准得到的位置参数作为精配准的初始解，利用粒子群算法结合Powell算法得到六个空间变换参数。本发明利用一个特征提取网络识别图像的点特征而非物体的形状特征，提高了对高噪声图像的鲁棒性；采用自监督和半监督的训练模式，并添加了图像增强模块，以适应不同模态的图像质量；再利用改进的粒子群算法结合Powell算法进行精配准确定六个空间变换参数，在保证精度的同时，利用深度学习提高了配准速度。利用深度学习提高了配准速度。利用深度学习提高了配准速度。


技术研发人员：谢世朋 林泽
受保护的技术使用者：南京邮电大学
技术研发日：2022.06.30
技术公布日：2022/11/1