一种医学图像的非刚性配准方法

1.本发明属于医学图像处理技术领域，涉及一种医学图像的配准方法，特别是基于深度学习的非刚性、无监督医学图像配准方法。


背景技术：

2.医学图像配准，是通过寻找某种空间变换,使两幅图像上所有的解剖点或者是所有具有医学诊断意义的参考点都达到对齐或者匹配对应。具体可以根据变换关系分为刚性变换和柔性变换两大类配准方法。刚性变换本质是全局的线性变换，典型方法有旋转、缩放、仿射等，但其缺点是无法模拟局部的几何差异，同时由于模型参数简单而产生的形变不够鲁棒。非刚性变换则克服了上述缺点，将一个稠密的位移场(又称形变场、变换场)作为模型参数，此时每个像素或体素都有自己的位移向量，因此该种变换具有局部性的变形能力，在复杂且多变的医学图像配准中广为使用。
3.针对医学图像配准，学者们设计了许多传统的配准算法，使得这一技术能通过计算机来自动化完成。但是这些传统配准方法仍存在一些局限性，比如计算速度慢、整个配准流程操作复杂以及无法利用海量数据等缺点。相比于传统配准方法，基于深度学习的方法可以把较长的迭代时间转移到模型训练时间中去，其实际推理时间可以远少于传统方法的迭代时间；另一方面，端到端的深度回归网络设计，替代了传统配准方法的复杂流程，可以大大减轻配准的操作难度。还有就是深度学习网络可以通过对海量数据的拟合，来获得更泛化、更鲁棒的配准结果。
4.由于金标准标注稀缺，人工标注成本高，因此有监督医学图像配准的发展受到制约，而无监督的医学图像配准对标注的要求低或无需标注，有更好的发展前景。无监督深度学习的图像配准起源自jaderberg将空间变换器(stn)运用于配准网络，可以基于网络预测的形变场，对源图像进行插值和采样后得到形变后的图像，从而可以在训练过程中计算图像的相似性损失，构成损失函数而推动深度网络的反向传播，达成训练学习的目的。balakrishnan提出的配准框架voxelmorph，其利用3d u-net输出全局的3d形变场，且过程中无需生成位移参数的中间量。该方法一经提出，其优异的性能使得无监督深度学习的配准方法备受关注。然而，该框架对有些变形程度大的图像学习效果并不好，这是因为网络捕捉长距离形变的能力不够。为此，有许多基于“由粗到细”方式的方法被提出。zhao采用级联的方式，逐次对图像进行递进式配准，级联的每一阶段都是一个voxelmorph，这种方法虽然有效，但是模型参数较多，配准性能和速度都受到了限制。hu提出了一个名为dual-prnet的双流金字塔配准网络，其利用图像金字塔，从解码器产生的特征图直接产生多个不同尺度的配准场。


技术实现要素：

5.本发明要解决的技术问题是：当今方法对变形较大的图像配准不好，并且当今常用的“由粗到细”方式存在误差累积的问题。因此提出了一种基于自适应相关与不确定度估
计的医学图像的非刚性配准方法。
6.本发明的技术方案如下：
7.一种医学图像的非刚性配准方法，步骤如下：
8.步骤(1)金字塔特征编码
9.该步骤的输入为浮动图像和固定图像，然后分别送入一个共享参数的金字塔特征编码网络，从而得到不同采样倍率的浮动特征图或固定特征图。金字塔特征编码网络的作用是生成不同采样倍率的特征图，即一个特征金字塔。特征编码网络的主要组件是3d卷积层，并以resnet的残差块结构为基本单位进行堆叠，且通过设置卷积层的stride而得到不同的采样倍率。这种编码网络可以很好地将三维医学图像进行自编码，提取适合进行配准的特征。接下来的步骤则以每一种倍率的特征图为一个阶段，并于最高的采样倍率，即最小尺寸的特征图开始，逐个阶段地生成更大尺寸、更精密的场。
10.步骤(2)阶段预处理
11.从步骤(2)开始到步骤(4)生成位移场结束，定义为一个阶段。
12.该步骤输入首先将上一阶段的位移场φ、不确定度σ进行上采样，从而得到与特征图相同的尺寸；如果是第一个阶段，则初始化位移场与不确定度为零向量。接着将步骤(1)获得的当前阶段的浮动特征进行形变，只需通过上采样后的位移场φ即可完成，得到当前阶段的形变后的浮动特征，简称为形变特征：
13.feats
warp
(c,x,y,z)＝φ(x,y,z)
⊙
feats
mov
(c,x,y,z)
14.其中x、y、z表示体素位置，对应3d图像的三个坐标轴；feats
mov
和feats
warp
分别为浮动特征和形变特征；位移场φ包含了3个3维矩阵，每个矩阵尺寸与浮动特征一致，表示从浮动特征向固定特征进行变换的关系矩阵，
⊙
则为形变操作，且不同通道所做的形变操作是相同的。
15.步骤(3)生成位移空间
16.该步骤输入由两部分组成，一部分是当前阶段的浮动特征图，将浮动特征图送入多层卷积器，生成位移微调值，具体形式为一个四维张量：
17.disp
shift
(t,x,y,z)＝convs[feats
mov
(c,x,y,z)]
[0018]
其中disp
shift
表示位移微调值，而t表示位移空间中的点数；其中的c表示特征图的通道数；convs表示一系列的自网络结构，由多层的3d卷积、实例归一化和非线性激活函数组成。
[0019]
该步骤输入的另一部分是上一阶段的不确定度，将其通过一个线性变换而生成位移基数值。
[0020]
disp
space
(t，x，y，z)＝[α*σ(x，y，z)+β]*space(t)
[0021]
其中disp
space
表示位移基数值；σ表示不确定度，α和β表示线性变换的参数，space表示三维空间中一个立方体状的采样空间，其点数为t个；做第二个乘法前，需要将张量扩展到同一维度。
[0022]
将位移基数、位移微调值进行加和，则得到一个坐标向量，即最终的离散位移空间：
[0023]
disp(t，x，y，z)＝disp
space
(t，x，y，z)+disp
shift
(t，x，y，z)
[0024]
其中disp表示最终的位移空间，具体为一个坐标值向量。
[0025]
该位移空间是根据浮动图像的结构先验信息而进行动态变换的，因此使得配准拥有了绕开噪声、避免位移场不连续的能力；同时可以精准扩大采样范围，有效应对大变形的图像配准，从而提高初步匹配的精确度。
[0026]
步骤(4)多头场估计
[0027]
该步骤的输入由两部分组成，第一部分为通过步骤(2)获得的当前阶段的形变特征以及通过步骤(1)获得的固定特征图，另一部分为(3)步骤生成的自适应坐标(位移空间disp)，通过多头相关、代价聚合、归一化和位移融合后得到位移场，并同时估计该场的不确定度。具体而言，先将形变特征feats
warp
根据自适应坐标进行再采样，注意是进行t轮采样，因此生成的特征会增加一个t维度：
[0028]
feats
disp
(c，t，x，y，z)＝disp(t，x，y，z)
⊙
feats
warp
(c，x，y，z)
[0029]
其中feats
disp
表示采样得到的多位移特征。
[0030]
将多位移特征、固定特征feats
fix
分别划分成k个头，本质就是对通道进行分组，即k_divide操作：
[0031][0032][0033]
然后在维度上进行相关运算，得到k头相似代价体sim：
[0034][0035]
其中表示内积运算。
[0036]
然后通过k头卷积网络进行代价聚合，最终得到代价体cost：
[0037]
cost(k,t,x,y,z)＝k_convs[sim(k,t,x,y,z)]
[0038]
其中k_convs表示一系列的自网络结构，由多层的3d卷积、实例归一化、非线性激活函数组成，且卷积层的分组数为k。
[0039]
再对所有不同头生成的位移进行归一化运算，得到位移空间所对应的概率向量pro：
[0040]
pro(k,t,x,y,z)＝softmax[cost(k,t,x,y,z)]
[0041]
其中softmax表示归一化函数，此处在k和t两个维度上进行。
[0042]
最后将位移空间disp与概率向量pro进行内积，得到融合后的位移场φ：
[0043][0044]
其中表示内积运算，且在k和t两个维度上进行，且进行前将张量进行扩展。由于位移场φ本质是一阶矩结果，而二阶中心矩的根号结果就得到了不确定度σ。具体计算过程如下：
[0045][0046]
当生成了与浮动特征图同种尺度的位移场和不确定度后，则认为该阶段完成，即将开始下一阶段，跳转至(2)步骤；如果所有阶段均完成，则进入(5)步骤。
[0047]
步骤(5)构建无监督损失函数
[0048]
得到与原图相同尺度的位移场φ后，对浮动图像i
mov
进行形变，得到的形变图像即为最终的配准结果i
φ(mov)
：
[0049]iφ(mov)
＝φ(x,t,z)
⊙imov
[0050]
由于训练神经网络模型需要构建损失函数，本发明则使用无监督损失函数，其构建由两部分组成：一是相似性损失函数二是位移场的平滑性损失函数将二者加和得到最终的损失函数，λ则为二者的平衡参数，又称正则化系数：
[0051][0052]
本发明的有益效果：该方法通过自适应相关，即生成自适应的位移空间来精准控制特征的采样范围，可有效提高大变形的图像配准效果，还能应对噪声干扰、弱纹理和位移场不连续等情形，提高配准后的相似度；另外还使用不确定度估计来缓解误差累积的问题，以及使用多头机制深化配准。
附图说明
[0053]
图1是本发明方法的基本原理图。
[0054]
图2为本发明方法提供的整体流程图。
[0055]
图3为本发明方法提供的步骤(3)的详细示意图。
[0056]
图4为本发明方法提供的步骤(4)的详细示意图。
[0057]
图5为本发明方法提供的配准生成的位移场示例图。
[0058]
图6(a)和图6(b)为本发明方法提供的配准结果图像对比示例图，其中，图6(a)是固定图像(目标)，图6(b)是形变图像(结果)。
具体实施方式
[0059]
以下结合附图和技术方法，进一步说明本发明的具体实施方式。
[0060]
本发明的一种医学图像的非刚性配准方法的基本流程如图1和图2所示，基本步骤包括：(1)输入浮动、固定图像，送入金字塔特征编码器而生成不同采样倍率的特征图；(2)从最高的采样倍率(即最小尺寸)的特征开始，逐阶段地生成自适应的位移空间、位移场，并估计不确定度；(3)根据上一阶段的不确定度与当前阶段的浮动特征，分别形成位移基数与位移微调值，再加和得到自适应的位移空间；(4)根据当前阶段的形变特征、目标特征以及位移空间，通过多头相关、代价聚合、归一化和位移融合步骤后得到位移场，并同时估计该场的不确定度；(5)重复步骤(2)至(4)，直至获得与原图尺寸相同的位移场，形变后得到最终的配准结果。
[0061]
本实施例以肺部ct图像的呼吸双相配准，具体配准过程如下：
[0062]
步骤(1)金字塔特征编码
[0063]
该步骤输入肺部ct的吸气相和呼气相，分别作为浮动、固定图像(如图2中左侧所示)送入共享参数的特征编码器(如图2中顶端与底端的所示箭头)，得到输入图像尺寸的1/2、1/4、1/8、1/16的四种特征图，具体通过的是3d卷积层、in归一化层和pooling层，每次pooling得到下一种尺寸的特征图，而每次pooling前有2、3、4、6层3d卷积，且每层3d卷积都
有跳接结构，并跟着一个in归一化层，最终生成的四种特征图分别有12、24、24、24的通道维度。
[0064]
步骤(2)阶段预处理
[0065]
该步骤的输入首先将上一阶段的位移场φ、不确定度σ进行上采样，具体方式为三线性插值，从而得到与特征图相同的尺寸；如果是第一个阶段，则初始化位移场与不确定度为零向量。接着将步骤(1)获得的当前阶段的浮动特征进行形变(如图2中所示模块)，将上采样后的位移场φ作为坐标进行采样即可，得到当前阶段的形变后的浮动特征，简称为形变特征。
[0066]
步骤(3)生成位移空间
[0067]
如图2中ce模块，具体步骤如图3所示。具体选取位移空间中的点数为27个，对应着三维空间中每个轴上的取值分别为{-1，0，1}，总共构成33＝27个点的采样空间。位移微调值则由浮动特征通过一系列的自卷积结构，该结构由三个沙漏式的卷积模块堆叠而成，其中每个沙漏卷积模块有四层卷积，前三层卷积的stride分别为1、2、1，而第四层为stride为2的反卷积，并由跳跃结构连接，通过这四层之后特征尺度不变；位移基数则由不确定度通过线性变化而生成。最终将位移微调值与位移基数进行加和，得到位移空间。
[0068]
步骤(4)多头场估计
[0069]
如图2中de模块，具体步骤如图4所示。首先将形变特征根据位移空间进行再采样，例如空间点数为27，通道数为24时，则将24*h*w*d的向量变成了27*24*h*w*d。再对两种特征进行分组，例如将24的通道维度分成4组，每组有6的维度。得到相似代价体后，通过k头卷积网络，进行代价聚合。k头卷积网络定义为三层沙漏卷积模块(类似图3中的上侧网络结构)，且最后一层沙漏卷积包含一个拼接型连接，另外k头需定义卷积层的分组数为k。然后对所有不同头生成的特征进行归一化运算，得到位移空间对应的概率向量；最后将位移空间与概率向量进行内积，得到融合后的位移场。最后对位移场进行二阶中心矩分析，得到不确定度。
[0070]
当生成了与浮动特征图同种尺度的位移场和不确定度后，则认为该阶段完成，即将开始下一阶段，跳转至(2)步骤；如果所有阶段均完成，则进入(5)步骤。
[0071]
步骤(5)构建监督损失函数
[0072]
如图2所示，得到与原图尺度一致的位移场后，对浮动图像进行形变并与固定图像构建相似性损失函数，同时构建一个平滑性损失函数，二者加和得到最终的损失函数。
[0073]
步骤(6)训练、验证方法
[0074]
网络优化器采用adam，使用指数衰减的学习率曲线，初始学习率设为0.0005，并在每1万次迭代后将学习率减小为原来的0.95，共进行10万次迭代，并采用移动平均的方法更新模型参数。如果数据集样本不充裕时，则使用k-fold交叉验证的方法。
[0075]
步骤(7)评价指标
[0076]
由于ct图像中肺实质区域往往会占据大量的体素，因此使用地标点间的tre来评估肺内部的配准质量，也就是配准后的相似程度。
[0077]
[0078]
其中xf,xm分别表示浮动图像和固定图像间所对应的地标点集合，通常为手工标记或机器生成的具有解剖学价值的参考点。除了此之外还有平滑性指标，平滑性主要考察形变图像的折叠体素。在形变场φ而言，折叠体素出现的位置可以用雅可比行列式来求。定义坐标为(i,j,k)的体素，雅可比行列式值：
[0079][0080]
平滑性指标一般定义为雅克比行列式为负值的体素个数
[0081][0082]
本实施例获得的位移场如图5所示，配准结果图像对比如图6(a)和图6(b)所示。采用评价指标进行评价，在保证平滑性的情况下，即雅可比行列式负值个数为0，发现形变图像与固定图像相似程度高，配准结果较好，证明了本发明方法的有效性。

技术特征：
1.一种医学图像的非刚性配准方法，其特征在于，步骤如下：步骤(1)金字塔特征编码该步骤的输入为浮动图像和固定图像，然后分别送入一个共享参数的金字塔特征编码网络，从而得到不同采样倍率的浮动特征图或固定特征图；金字塔特征编码网络的作用是生成不同采样倍率的特征图，即一个特征金字塔；特征编码网络的主要组件是3d卷积层，并以resnet的残差块结构为基本单位进行堆叠，且通过设置卷积层的stride而得到不同的采样倍率；接下来的步骤以每一种倍率的特征图为一个阶段，并于最高的采样倍率，即最小尺寸的特征图开始，逐个阶段地生成更大尺寸、更精密的场；步骤(2)阶段预处理该步骤输入首先将上一阶段的位移场φ、不确定度σ进行上采样，从而得到与特征图相同的尺寸；如果是第一个阶段，则初始化位移场与不确定度为零向量；接着将步骤(1)获得的当前阶段的浮动特征进行形变，只需通过上采样后的位移场φ即可完成，得到当前阶段的形变后的浮动特征，简称为形变特征：feats
warp
(c,x,y,z)＝φ(x,y,z)
⊙
feats
mov
(c,x,y,z)其中x、y、z表示体素位置，对应3d图像的三个坐标轴；feats
mov
和feats
warp
分别为浮动特征和形变特征；位移场φ包含了3个3维矩阵，每个矩阵尺寸与浮动特征一致，表示从浮动特征向固定特征进行变换的关系矩阵，
⊙
则为形变操作；不同通道所做的形变操作是相同的；步骤(3)生成位移空间该步骤输入由两部分组成，一部分是当前阶段的浮动特征图，将浮动特征图送入多层卷积器，生成位移微调值，具体形式为一个四维张量：disp
shift
(t,x,y,z)＝convs[feats
mov
(c,x,y,z)]其中disp
shift
表示位移微调值，而t表示位移空间中的点数；其中的c表示特征图的通道数；convs表示一系列的自网络结构，由多层的3d卷积、实例归一化和非线性激活函数组成；该步骤输入的另一部分是上一阶段的不确定度，将其通过一个线性变换而生成位移基数值；disp
space
(t,x,y,z)＝[α*σ(x,y,z)+β]*space(t)其中disp
space
表示位移基数值；σ表示不确定度，α和β表示线性变换的参数，space表示三维空间中一个立方体状的采样空间，其点数为t个；做第二个乘法前，需要将张量扩展到同一维度；将位移基数、位移微调值进行加和，则得到一个坐标向量，即最终的离散位移空间：disp(t,x,y,z)＝disp
space
(t,x,y,z)+disp
shift
(t,x,y,z)其中disp表示最终的位移空间，具体为一个坐标向量；步骤(4)多头场估计该步骤的输入由两部分组成，第一部分为通过步骤(2)获得的当前阶段的形变特征以及通过步骤(1)获得的固定特征图，另一部分为(3)步骤生成的自适应坐标即位移空间disp，通过多头相关、代价聚合、归一化和位移融合后得到位移场，并同时估计该场的不确定度；具体而言，先将形变特征feats
warp
根据自适应坐标进行采样，是进行t轮采样，因此生成的特征会增加一个t维度：feats
disp
(c,t,x,y,z)＝disp(t,x,y,z)
⊙
feats
warp
(c,x,y,z)
其中feats
disp
表示采样得到的多位移特征；将多位移特征、固定特征feats
fix
分别划分成k个头，本质就是对通道进行分组，即k_divide操作：divide操作：然后在维度上进行相关运算，得到k头相似代价体sim：其中表示内积运算；然后通过k头卷积网络进行代价聚合，最终得到代价体cost：cost(k,t,x,y,z)＝k_convs[sim(k,t,x,y,z)]其中k_convs表示一系列的自网络结构，由多层的3d卷积、实例归一化、非线性激活函数组成，且卷积层的分组数为k；再对所有不同头生成的位移进行归一化运算，得到位移空间所对应的概率向量pro：pro(k,t,x,y,z)＝softmax[cost(k,t,x,y,z)]其中softmax表示归一化函数，此处在k和t两个维度上进行；最后将位移空间disp与概率向量pro进行内积，得到融合后的位移场φ：其中表示内积运算，且在k和t两个维度上进行，且进行前将张量进行扩展；由于位移场φ本质是一阶矩结果，而二阶中心矩的根号结果就得到了不确定度σ；具体计算过程如下：当生成了与浮动特征图同种尺度的位移场和不确定度后，则认为该阶段完成，即将开始下一阶段，跳转至(2)步骤；如果所有阶段均完成，则进入(5)步骤；步骤(5)构建无监督损失函数得到与原图相同尺度的位移场φ后，对浮动图像i
mov
进行形变，得到的形变图像即为最终的配准结果i
φ(mov)
：i
φ(mov)
＝φ(x,y,z)
⊙
i
mov
训练神经网络模型需要构建损失函数，使用无监督损失函数，其构建由两部分组成：一是相似性损失函数二是位移场的平滑性损失函数将二者加和得到最终的损失函数，λ则为二者的平衡参数，又称正则化系数：

技术总结
本发明提供了一种医学图像的非刚性配准方法，具体步骤包括：(1)输入浮动、固定图像，送入金字塔特征编码器而生成不同采样倍率的特征图；(2)从最高的采样倍率的特征开始，逐阶段地生成自适应的位移空间、位移场，并估计不确定度；(3)根据上一阶段的不确定度与当前阶段的浮动特征，分别形成位移基数与位移微调值，再加和得到自适应的位移空间；(4)根据当前阶段的形变特征、目标特征以及位移空间，通过多头相关、代价聚合、归一化和位移融合步骤后得到位移场，并同时估计该场的不确定度；(5)重复步骤(2)至(4)，直至获得与原图尺寸相同的位移场，形变图像后得到最终的配准结果。形变图像后得到最终的配准结果。形变图像后得到最终的配准结果。


技术研发人员：张立和 张恒宇 卢湖川
受保护的技术使用者：大连理工大学
技术研发日：2023.07.17
技术公布日：2023/10/15