一种基于多层P样条和稀疏编码的医学图像配准方法转让专利
申请号 : CN201710259589.0
文献号 : CN107025650B
文献日 : 2019-10-29
发明人 : 王丽芳 , 成茜 , 董侠 , 王雁丽 , 史超宇
申请人 : 中北大学
摘要 :
本发明公开了一种基于多层P样条和稀疏编码的医学图像配准方法，涉及图像处理技术领域，在P样条的基础上让网格控制节点由少变多，直至在某个网格密度下配准的误差最小，并使用基于图像块稀疏编码的相似性测度，不仅考虑了医学图像中存在灰度不均匀性造成的灰度偏移场，也考虑了像素之间的空间依赖性，同时使用K‑SVD算法，相比确定的字典，适用范围更广。
权利要求 :
1.一种基于多层P样条和稀疏编码的医学图像配准方法，其特征在于，包括：使用步长为1的滑动窗把含有灰度偏移场的参考图像R和浮动图像F划分成为大小为的图像块，把两幅图像的图像块分别记为：IR＝[R1,R2,…,Rn],IF＝[F1,F2,…,Fn] (1)训练集为两幅图像中图像块的集合：
初始化多层P样条变换模型参数和优化模块的参数；
使用K-SVD算法训练图像块，记录训练得到的分析字典Ω，寻找每个图像块的稀疏系数Yi，进而计算两幅图像中图像块的稀疏表示αRn和αFn；
通过L1SM计算图像块的相似性程度函数，将相似性程度函数L1SM作为优化模块的目标函数C(R,F)的第一部分Csim(R,F)，另一部分加入相关的平滑变换Csmooth，对目标函数迭代优化，进而更新变形场；
使用多层P样条变换更新浮动图像，并判断是否达到迭代次数，如果没有达到迭代次数，则重新使用K-SVD算法训练图像块，直到达到迭代次数；
输出配准后的浮动图像和控制网格。
2.如权利要求1所述的方法，其特征在于，多层P样条变换模型如下：其中，控制点位置为 [.]
表示取整数值运算，βl(u)表示三次P样条第l个基函数，总的惩罚P为：其中，P0，P1，P2，P3分别代表边缘基B0，B1，B2，B3的惩罚，λ0，λ1，λ2，λ3为平滑系数，Ik0，Ik1，Ik2，Ik3为单位矩阵。
3.如权利要求2所述的方法，其特征在于，步骤使用K-SVD算法训练图像块，记录训练得到的分析字典Ω，寻找每个图像块的稀疏系数Yi，进而计算两幅图像中图像块的稀疏表示αRn和αFn包括稀疏编码和字典更新两个阶段：在稀疏编码阶段，计算图像块 在分析字典Ω上的稀疏系数Yi：其中，||Y||0表示稀疏系数向量Y中非零元素的个数，ε表示允许偏差的精度，上述公式的求解过程即为稀疏编码；
在字典更新阶段，对字典中的每个原子进行更新，假定稀疏系数向量Y和分析字典Ω都是固定的，待更新字典的第k列为Ωk，令稀疏系数矩阵中Y与Ωk相乘的第k行为 则式(4)可以写为：公式(6)为除第k个原子以外其他原子产生的表示误差，定义为含有原子Ωk成分的图像块αi的索引所组成的集合，定义大小为N×wk的矩阵Dk，矩阵元素(wk(i),i)为1，其他矩阵元素均为0，定义向量 长度为|wk|，定义矩阵大小为n×|wk|，是去掉不受原子Ωk影响的样本中带来的误差，此时公式(6)等价于：将误差矩阵 进行奇异值SVD分解，得到 其中Δ中的奇异值是由大到小排列的，记 为矩阵U中的第一列，使用 来更新字典中的Ωk，同时将矩阵V中的第一列与Δ(1,1)的乘积来更新稀疏系数矩阵中的 实现字典的更新；
通过交替执行稀疏编码和字典更新两个阶段，即可得到分析字典Ω和两幅图像中图像块的稀疏表示αRn和αFn。
4.如权利要求3所述的方法，其特征在于，步骤通过L1SM计算图像块的相似性程度函数，将相似性程度函数L1SM作为优化模块的目标函数C(R,F)的第一部分Csim(R,F)，另一部分加入相关的平滑变换Csmooth，对目标函数迭代优化，进而更新变形场具体为：通过L1SM计算图像块的相似性程度函数：
L1SM＝||Ω(R-F)||1＝||αR-αF||1 (8)将相似性程度函数L1SM作为优化模块的目标函数C(R,F)的第一部分Csim(R,F)，另一部分加入相关的平滑变换Csmooth，使待配准的两幅图像之间保持权衡，并且变换光滑；
其中，A是区域的面积，λ是权重系数。
5.如权利要求4所述的方法，其特征在于，使用梯度下降法迭代优化目标函数，进而更新变形场，L1SM的导数表达式如下：其中，▽F为浮动图像的梯度，θ表示变换参数，sign(.)为符号函数。
说明书 :
一种基于多层P样条和稀疏编码的医学图像配准方法
技术领域
[0001] 本发明涉及图像处理技术领域，特别是涉及一种基于多层P样条和稀疏编码的医学图像配准方法。
背景技术
[0002] 医学图像配准是指对浮动图像寻求某种空间变换，使它与参考图像上的同一解剖点或者具有诊断意义的关键点达到空间上的一致。当两幅图像配准完成后，可以对它们进行比较和分析。医学图像配准作为医学图像处理中的一个重要组成部分，对手术定位，图像融合等具有重要意义。
[0003] 医学图像配准方法主要有刚性配准和非刚性配准两大类。刚性配准仅适用于不存在形变的配准。当浮动图像与参考图像存在较大形变时，就需要通过非线性的变换模型描述复杂的变换过程。非刚性配准方法主要有基于物理模型的变换和基于函数描述的变换。基于物理模型的变换将图像间的差异看作是由某种物理变形引起的，通常计算量较大，难以准确建立模型。基于函数描述的变换来源于插值和近似理论，它们采用具有少量参数的基函数描述复杂，稠密的非线性几何变换域。目前常用于非刚性配准的基函数有：径向基函数、B样条函数和小波函数等。基于B样条的自由形式变换(Free-Form Deformation，FFD)是目前最为流行的非刚性图像变换方法。Brian D.Marx等人在文献中在B样条函数的基础上附加了正则项，提出了P样条方法。Pradhan等人在文献中将P样条方法运用到医学图像配准领域。汪军等人在文献中将P样条和局部互信息用于非刚性医学图像配准。与B样条变换模型相比，P样条附加了额外的正则项来避免形变场奇异点和折叠效应，但单层P样条的配准精度仍受控制网格稀疏程度的影响。
[0004] 在图像配准中，相似性测度也是其中一个重要组成部分，常用的基于灰度的相似性测度有平方差(SSD)、绝对差(SAD)、相关系数(CC)、互信息(MI)等。在医学图像中，受时间和成像条件的影响，待配准的图像的灰度场可能发生显著的变化。例如，脑磁共振图像中存在缓慢变化的灰度偏移场，然而，现存的相似性测度对于含有灰度不均匀性的图像鲁棒性较差。为了解决这个问题，Myronenko等人在文献中提出在残差复杂性(RC)的分析上解决灰度级偏移场，RC的表达式中不含灰度偏移场，从而自适应地约束了灰度偏移场，但没有考虑残差图像自身特性，因而造成误配。卢振泰等人在文献中提出了局部方差和残差复杂性的医学图像配准方法，它考虑残差图像自身的特性，比RC更适用于非均匀医学图像配准，但它的计算仍以两幅图像之间的残差图像为基础，缺乏普适性。Ghaffari等人在文献中提出了稀疏诱导的相似性测度(SISM)，其是稀疏的方法，参考了像素之间的空间依赖性，但它的字典是基于DCT和小波变换的，利用了确定的字典，虽易快速实现，但表示能力有局限性。此外，Ghaffari等人又在文献中提出了Rank Induced相似性测度(RISM)，它把图像配准视为一个非线性和低秩矩阵分解问题，该方法可以产生较准确的配准结果，但适用范围较小，只适用于单模态医学图像配准。
发明内容
[0005] 本发明实施例提供了一种基于多层P样条和稀疏编码的医学图像配准方法，可以解决现有技术中存在的问题。
[0006] 一种基于多层P样条和稀疏编码的医学图像配准方法，包括：
[0007] 使用步长为1的滑动窗把含有灰度偏移场的参考图像R和浮动图像F划分成为大小为 的图像块，把两幅图像的图像块分别记为：
[0008] IR＝[R1,R2,…,Rn],IF＝[F1,F2,…,Fn] (1)
[0009] 训练集为两幅图像中图像块的集合：
[0010]
[0011] 初始化多层P样条变换模型参数和优化模块的参数；
[0012] 使用K-SVD算法训练图像块，记录训练得到的分析字典Ω，寻找每个图像块的稀疏系数Yi，进而计算两幅图像中图像块的稀疏表示αRn和αFn；
[0013] 通过L1SM计算图像块的相似性程度函数，将相似性程度函数L1SM作为优化模块的目标函数C(R,F)的第一部分Csim(R,F)，另一部分加入相关的平滑变换Csmooth，对目标函数迭代优化，进而更新变形场；
[0014] 使用多层P样条变换更新浮动图像，并判断是否达到迭代次数，如果没有达到迭代次数，则重新使用K-SVD算法训练图像块，直到达到迭代次数；
[0015] 输出配准后的浮动图像和控制网格。
[0016] 优选地，多层P样条变换模型如下：
[0017]
[0018] 其中，控制点位置为[.]表示取整数值运算，βl(u)表示三次P样条第l个基函数，总的惩罚P为：
其中，P0，P1，P2，P分别代表边缘基B0，B1，B2，
B3的惩罚，λ0，λ1，λ2，λ3为平滑系数，Ik0，Ik1，Ik2，Ik3为单位矩阵。
其中，P0，P1，P2，P分别代表边缘基B0，B1，B2，
B3的惩罚，λ0，λ1，λ2，λ3为平滑系数，Ik0，Ik1，Ik2，Ik3为单位矩阵。
[0019] 优选地，步骤使用K-SVD算法训练图像块，记录训练得到的分析字典Ω，寻找每个图像块的稀疏系数Yi，进而计算两幅图像中图像块的稀疏表示αRn和αFn包括稀疏编码和字典更新两个阶段：
[0020] 在稀疏编码阶段，计算图像块I＝[R1,R2,…,Rn,F1,F2,…,Fn]在分析字典Ω上的稀疏系数Yi：
[0021]
[0022] 其中，||Y||0表示稀疏系数向量Y中非零元素的个数，ε表示允许偏差的精度，上述公式的求解过程即为稀疏编码；
[0023] 在字典更新阶段，对字典中的每个原子进行更新，假定稀疏系数向量Y和分析字典Ω都是固定的，待更新字典的第k列为Ωk，令稀疏系数矩阵中Y与Ωk相乘的第k行为 则式(4)可以写为：
[0024]
[0025]
[0026] 公式 (6) 为除 第k个原子以 外其他原子产生的 表示误差 ，定 义为含有原子Ωk成分的图像块αi的索引所组成的集合，定义大小为
N×wk的矩阵Dk，矩阵元素(wk(i),i)为1，其他矩阵元素均为0，定义向量 长度为
|wk|，定义矩阵 大小为n×|wk|，是去掉不受原子Ωk影响的样本中带来的误差，
此时公式(6)等价于：
N×wk的矩阵Dk，矩阵元素(wk(i),i)为1，其他矩阵元素均为0，定义向量 长度为
|wk|，定义矩阵 大小为n×|wk|，是去掉不受原子Ωk影响的样本中带来的误差，
此时公式(6)等价于：
[0027]
[0028] 将误差矩阵 进行奇异值SVD分解，得到 其中Δ中的奇异值是由大到小排列的，记 为矩阵U中的第一列，使用 来更新字典中的Ωk，同时将矩阵V中的第一列与Δ(1,1)的乘积来更新稀疏系数矩阵中的 实现字典的更新；
[0029] 通过交替执行稀疏编码和字典更新两个阶段，即可得到分析字典Ω和两幅图像中图像块的稀疏表示αRn和αFn。
[0030] 优选地，步骤通过L1SM计算图像块的相似性程度函数，将相似性程度函数L1SM作为优化模块的目标函数C(R,F)的第一部分Csim(R,F)，另一部分加入相关的平滑变换Csmooth，对目标函数迭代优化，进而更新变形场具体为：
[0031] 通过L1SM计算图像块的相似性程度函数：
[0032] L1SM＝||Ω(R-F)||1＝||αR-αF||1 (8)
[0033] 将相似性程度函数L1SM作为优化模块的目标函数C(R,F)的第一部分Csim(R,F)，另一部分加入相关的平滑变换Csmooth，使待配准的两幅图像之间保持权衡，并且变换光滑；
[0034]
[0035] 其中，A是区域的面积，λ是权重系数。
[0036] 优选地，使用梯度下降法迭代优化目标函数，进而更新变形场，L1SM的导数表达式如下：
[0037]
[0038] 其中， 为浮动图像的梯度，θ表示变换参数，sign(.)为符号函数。
[0039] 本发明实施例提供的一种基于多层P样条和稀疏编码的医学图像配准方法，在P样条的基础上让网格控制节点由少变多，直至在某个网格密度下配准的误差最小，并使用基于图像块稀疏编码的相似性测度，不仅考虑了医学图像中存在灰度不均匀性造成的灰度偏移场，也考虑了像素之间的空间依赖性，同时使用K-SVD算法，相比确定的字典，适用范围更广。
附图说明
[0040] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0041] 图1为本发明实施例提供的一种基于多层P样条和稀疏编码的医学图像配准方法的流程图；
[0042] 图2为一组头颅矢状面图像的配准结果，a为参考图像，b为浮动图像，c为多层P样条+LMI，d为多层P样条+LMI网格，e为多层P样条+RC，f为多层P样条+RC网格，g为多层P样条+L1SM，h为多层P样条+L1SM网格；
[0043] 图3为一组脑部MR图像在不同网格密度下的配准结果，a为参考图像，b为浮动图像，c为L1SM(4×4)，d为L1SM(4×4)网格，e为L1SM(8×8)，f为L1SM(8×8)网格。
具体实施方式
[0044] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0045] 参照图1，本发明实施例中提供的一种基于多层P样条和稀疏编码的医学图像配准方法，包括以下步骤：
[0046] 步骤100，使用步长为1的滑动窗把含有灰度偏移场的参考图像R和浮动图像F划分成为大小为 的图像块，对于图像R和F都有 个图像块，把两幅图像的图像块分别记为：
[0047] IR＝[R1,R2,…,Rn],IF＝[F1,F2,…,Fn] (1)
[0048] 训练集为两幅图像中图像块的集合：
[0049]
[0050] 步骤110，初始化多层P样条变换模型参数和优化模块的参数，在本实施例中网格初始大小设为4×4，参考图像和浮动图像之间的配准误差为10-4，梯度下降算法的最大迭代次数为30次；多层P样条变换模型如下：
[0051]
[0052] 其中，控制点位置为[.]表示取整数值运算，βl(u)表示三次P样条第l个基函数。
[0053]
[0054] 总的惩罚P为：
[0055] 其中，P0，P1，P2，P分别代表边缘基B0，B1，B2，B3的惩罚，λ0，λ1，λ2，λ3为平滑系数，Ik0，Ik1，Ik2，Ik3为单位矩阵。
[0056] 步骤120，使用K-SVD算法训练图像块，记录训练得到的分析字典Ω，寻找每个图像块的稀疏系数Yi，本步骤包括稀疏编码和字典更新两个阶段：
[0057] 在稀疏编码阶段，计算图像块I＝[R1,R2,…,Rn,F1,F2,…,Fn]在分析字典Ω上的稀疏系数Yi：
[0058]
[0059] 其中，||Y||0表示稀疏系数向量Y中非零元素的个数，ε表示允许偏差的精度，上述公式的求解过程即为稀疏编码；
[0060] 在字典更新阶段，对字典中的每个原子进行更新，假定稀疏系数向量Y和分析字典Ω都是固定的，待更新字典的第k列为Ωk，令稀疏系数矩阵中Y与Ωk相乘的第k行为 则式(4)可以写为：
[0061]
[0062]
[0063] 公式 (6) 为除 第k个原子以 外其他原子产生的 表示误差 ，定 义为含有原子Ωk成分的图像块αi的索引所组成的集合，定义大小为
N×wk的矩阵Dk，矩阵元素(wk(i),i)为1，其他矩阵元素均为0，定义向量 (长度为
|wk|)，定义矩阵 (大小为n×|wk|，是去掉不受原子Ωk影响的样本中带来的误
差)，此时公式(6)等价于：
N×wk的矩阵Dk，矩阵元素(wk(i),i)为1，其他矩阵元素均为0，定义向量 (长度为
|wk|)，定义矩阵 (大小为n×|wk|，是去掉不受原子Ωk影响的样本中带来的误
差)，此时公式(6)等价于：
[0064]
[0065] 将误差矩阵 进行奇异值(SVD)分解，得到 (其中Δ中的奇异值是由大到小排列的)，记 为矩阵U中的第一列，使用 来更新字典中的Ωk，同时将矩阵V中的
第一列与Δ(1,1)的乘积来更新稀疏系数矩阵中的 实现字典的更新。
第一列与Δ(1,1)的乘积来更新稀疏系数矩阵中的 实现字典的更新。
[0066] K-SVD算法是一种迭代算法，通过交替执行稀疏编码和字典更新两个阶段，即可得到分析字典Ω和两幅图像中图像块的稀疏表示αRn和αFn；
[0067] 步骤130，通过L1SM计算图像块的相似性程度函数：
[0068] L1SM＝||Ω(R-F)||1＝||αR-αF||1 (8)
[0069] 将相似性程度函数L1SM作为优化模块的目标函数C(R,F)的第一部分Csim(R,F)，另一部分加入相关的平滑变换Csmooth，使待配准的两幅图像之间保持权衡，并且变换光滑。
[0070]
[0071] 其中，A是区域的面积，λ是权重系数。
[0072] 本实施例中使用梯度下降法迭代优化目标函数，进而更新变形场，L1SM的导数表达式如下：
[0073]
[0074] 其中， 为浮动图像的梯度，θ表示变换参数，sign(.)为符号函数。
[0075] 步骤140，使用多层P样条变换更新浮动图像，并判断是否达到迭代次数，如果没有达到迭代次数，则返回步骤120，直到达到迭代次数；
[0076] 步骤150，输出配准后的浮动图像和控制网格。
[0077] 实例
[0078] 为验证本发明方法的有效性，分别从两方面予以验证：
[0079] (1)相似性测度在配准中的影响
[0080] 图2给出了一组头颅矢状面图像的配准结果。2(a)、2(b)分别作为参考图像和浮动图像，分辨率为354×353。分别使用LMI，RC和本发明的L1SM三种测度进行对比实验，配准结果如图2所示。
[0081] 由图2(c)和2(d)可看出，使用LMI测度后的浮动图像在脑部多处均未配准，且形变网格在一些地方出现折叠现象；由图2(e)和2(f)可看出，使用RC测度只在颅骨内板，枕叶处没有配准，形变网格较为光滑；而由图2(g)和2(h)可看出，使用L1SM在各处均配准且形变网格光滑，并未出现折叠现象。
[0082] 为客观评价配准的效果，本发明采用均方根误差(RMSE)、配准时间(T/s)来定量评估配准的性能。均方根误差(RMSE)的数学表达式如公式16所示：
[0083]
[0084] 其中R(i,j)和F(i,j)分别是参考图像R，浮动图像F在点(i,j)处图像块的稀疏表示，M×N是图像的分辨率。RMSE值越小，配准效果越好。
[0085] 实验1中配准后图像和参考图像的RMSE值与配准运行时间如表1所示:
[0086] 表1相似性测度在配准中的影响
[0087]
[0088] 表1给出在多层P样条变换下局部互信息(LMI)、残差复杂性(RC)配准测度与本发明的L1SM配准结果的定量指标。表中数据是5次实验数据的平均值。依据对图像的分析和对评价指标均方根误差(RMSE)的分析，得出基于稀疏编码的L1SM相似性测度相比局部互信息(LMI)和残差复杂性(RC)对含有灰度偏移场的图像有着更好的配准表现；RMSE分别下降了
35.51％和31.74％。从配准时间上来看，基于稀疏编码的L1SM的效率要略低于LMI和RC的相似性测度。
35.51％和31.74％。从配准时间上来看，基于稀疏编码的L1SM的效率要略低于LMI和RC的相似性测度。
[0089] 实验2：网格稀疏程度对配准的影响
[0090] 图3给出一组脑部MR图像在不同网格密度下的配准结果，图像分辨率为557×583。其中，3(a)作为参考图像，3(b)作为浮动图像，图3(c)和3(d)显示了在4×4网格密度下的配准结果，图3(e)和3(f)显示了在8×8网格密度下的配准结果。
[0091] 由图3(d)和3(f)可看出，当控制网格稀疏时，由于控制顶点数目较少，每一个控制点影响的区域较大，导致配准精度受到影响。采用多层次P样条变换，控制网格由疏到密，模拟从全局到局部的变形，最终能够得到精确的配准结果。
[0092] 实验2中配准后图像和参考图像的RMSE值和配准运行时间如表2所示:
[0093] 表2网格稀疏程度对配准的影响
[0094]
[0095]
[0096] 表2给出在4×4，8×8不同网格密度下配准结果的定量指标。表中数据是5次实验数据的平均值。结果显示在8×8网格密度下RMSE较小，表明控制点越稠密，配准后两幅图像差异越小，配准效果越好，因此多层次配准可以通过逐层增加控制网格密度来准确选择网格大小，最终提升配准效果。
[0097] 依据上述两个实验的实验结果，表明了基于多层P样条和稀疏编码的非刚性医学图像配准方法降低了配准误差，提升了配准效果。
[0098] 本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
[0099] 本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0100] 这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0101] 这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0102] 尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
[0103] 显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。