一种三维医学图像配准方法、系统及计算机可读存储介质与流程

1.本发明涉及医学图像处理领域，尤其涉及一种三维医学图像配准方法、系统及计算机可读存储介质。


背景技术：

2.图像配准是图像处理的一个重要领域，配准指的是将两个或多个图像进行几何对齐，使浮动图像上的每一个点在目标图像上都有存在且唯一的点与其对应，其目标在于寻找不同图像之间对应的空间变换关系，从而尽可能地恢复待配准图像之间的几何不一致。图像配准是图像分析和处理的关键步骤,特别是在医学图像领域，对疾病诊断、手术引导和疾病治疗跟踪等方面具有重要应用价值，是图像融合、分析和目标识别的必要前提。医学图像配准在临床上有很多应用,图像引导手术、图像融合、器官图谱生成、肿瘤和骨骼生长监测等临床任务应用的关键技术，也是一个极具挑战性的问题。
3.随着医学成像技术的发展，可以获得不同模态的医学图像，例如：us，ct，mri，pet等，不同模态的图像具有不同的特点，us成像具有无创性、无放射性、实时性、便捷性和扫描成本低等优势，目前被广泛用于人体腹部器官和浅表器官疾病的筛查、诊断和与治疗，但是空间分辨率低、对比度低、穿透力有限、显像区域狭小、容易受医生操作技巧的影响等问题；ct成像逼真、清晰、空间分辨率高、显像区域完整，对于血管和骨头的造影效果好，可以用来突出解剖结构，将解剖结构与周围其他组织区分开来。但是ct成像过程中需要暴露在x射线下，有致癌的风险。mri核磁共振图像对人体的软组织器官有极佳的成像效果，而且其图像质量很好、分辨率高。但是体内带有金属异物的患者不能接受核磁共振检查，而且核磁共振有检查扫描时间长、器官的周期性运动易造成伪影、价格昂贵的缺点。pet可以检测人体的代谢情况，但其分辨率低、采集时间长。这些医学影像是无法相互取代的，若能够结合他们的优点将会对诊断治疗工作提供很大的帮助，因此配准技术的出现至关重要。配准技术可以提高检测治疗效果的效率，同时，该技术可以最大化地将不同模态或时间的医学图像融合，提高信息利用率和诊断的准确性。
4.目前图像配准的方式主要依靠手动校准对齐。例如：通过增强四期ct诊断有hcc风险的病人是否是hcc时，需要不断的调整四期图像，确保它们空间解剖位置对应起来，从而判断是否显示出对应的图像特征，最终根据这些特征进行诊断。但是在四期图像获取时，病人需要屏住呼吸，不同的呼吸程度会导致四期图像之间存在一定程度的形变，如果只是手动校准，无法到达解剖结构之间的对齐，从而影响诊断的精度。而且，人工校准是一个耗时的过程，因此通过计算机辅助将图像配准是必要的。
5.随着配准技术在图像分析领域的应用，越来越多的方法被用于医学图像领域，由于这个领域的特殊性，它对配准有更严格的要求：首先，配准过程产生的变换必须是一对一的，这样能保证浮动图像上的解剖点在目标图像能找到唯一对应的解剖点，这也保证配准过程是有意义的。其次，配准之后的图像与目标图像要尽可能的相似，即配准结束之后，两幅待配准的图像要尽可能的对齐。然后对变换有一些特定的要求，例如：四期配准时，变换
应该是体积保持的；对于不可压缩的组织配准时，变换应该满足不可压缩性；但是除去这些特定的要求，还必须得保证变换是没有折叠的，一旦变换出现网格折叠，就意味着这个变换是不合理、不能接受的。因此形变必须是没有网格折叠的。目前有几个方法可以达到这个要求，第一个方法是，通过加入时间维度，计算微分同胚的变化，例如：lddmm，但是这个方法对离散计算的格式要求很高，而且是一个非常耗时的过程。另一个方法是，在变换可微的情况下约束变换的雅可比行列式大于零。但是现有的约束雅可比行列式大于零的配准方法，例如：hyperelastic模型，它直接约束变换的雅可比行列式等于1，这只使用某些特定的应用，对于其他的应用会导致要求过于高而导致错误的匹配。
6.因此，需要一种新型的图像配准方法，其基于变分模型建立新的配准模型，以达到保拓扑配准的效果。


技术实现要素：

7.为了克服上述技术缺陷，本发明的目的在于提供一种三维医学图像配准方法、系统及计算机可读存储介质，实现三维医学图像的精确配准，且融合不同模态或时间的医学图像，提高信息利用率的准确性。
8.本发明公开了一种三维医学图像配准方法，包括以下步骤：
9.获取待配准图像，并对待配准图像作标准化预处理；
10.类比于复数空间的beltrami系数，于三维空间内定义一新测度，并基于新测度构建配准模型；
11.对配准模型作离散化处理，并使用高斯-牛顿迭代法对离散化处理后的配准模型优化，以得到最优变换模型；
12.基于最优变换模型构建待配准图像中的浮动图像，并将形变后的浮动图像与待配准图像中的目标图像融合显示。
13.优选地，获取待配准图像，并对待配准图像作标准化预处理的步骤包括：
14.搜集医学图像形成具有多个待配准图像的待配准图像集，并对每一待配准图像中的目标组织分割，形成分割图像；
15.对分割图像作灰度归一化处理、基于窗宽窗位作归一化处理，并对归一化处理后的分割图像作数字图像化重采样；
16.对分割图像裁剪，使得每一对待配准图像的尺寸一致。
17.优选地，类比于复数空间的beltrami系数，于三维空间内定义一新测度，并基于新测度构建配准模型的步骤包括：
18.类比于对复数空间的beltrami系数的定义，于三维空间内定义如下新测度：
19.其中表示雅可比行列式，sign()表示符号函数，且新测度的模小于1等价于变换的雅可比行列式大于0；
20.以新测度作为正则项，构建如下配准模型：
[0021][0022]
其中f:表示待配准图像中的目标图像，m:表示待配准图像中的浮动图像，y:r3→
r3表示空间变换，ω表示待配准图像的定义域，min
y l(y)表示最小化能量泛函，能量泛函中的第一项和第二项合起来表示数据保真项，化能量泛函，能量泛函中的第一项和第二项合起来表示数据保真项，表示归一化梯度算子，表示正常梯度算子，α1和α2表示光滑约束参数，β表示雅可比约束项的参数。
[0023]
优选地，以新测度作为正则项，构建如下配准模型的步骤包括：
[0024]
当待配准图像的模态为单模态时，数据保真项选取如下平方差之和ssd：ssd(y)＝，∫
ω
(m
°
y-f)2dx；
[0025]
当待配准图像的模态为多模态时，数据保真项选取如下梯度信息：
[0026]
优选地，对配准模型作离散化处理，并使用高斯-牛顿迭代法对离散化处理后的配准模型优化，以得到最优变换模型的步骤包括：
[0027]
基于中值公式对配准模型中的积分进行近似，并基于差分格式对配准模型中的导数进行离散，得到离散能量泛函；
[0028][0029]
其中φ1(rc)＝(φ1((rc)1)，
…
，φ1((rc)n))，φ(r)＝(φ((r)1)，
…
，φ((r)n))，
[0030]
对离散能量泛函求解出如下对应梯度：
[0031]
其中a和b分别对应一阶和二阶求导矩阵，un为离散后的位移场，h为离散单位，r(un)表示n(y)关于un离散的向量形式，rc(un)通过化简gf_tm得到；表示r(un)关于un的梯度向量，表示φ关于r(un)的梯度向量，表示φ1关于rc(un)的梯度向量；
[0032]
求解配准模型对应的如下海森矩阵：
[0033]
其中表示由φ1关于rc(un)的二阶导构成的对角矩阵，表示由φ关于r(un)的二阶导构成的对角矩阵；
[0034]
对海森矩阵求解可行方向，以求解对应的线性方程：
[0035][0036]
然后利用下述线性迭代求解满足需求的位移：
[0037]uk+1
＝uk+θkδuk[0038]yk+1
＝xn+u
k+1
，其中δu表示迭代方向，θ表示迭代步长，k表示第k次迭代，xn表示形变
图像的原始坐标；
[0039]
迭代后得到最细尺度对应的最优变换模型。
[0040]
本发明还公开了一种三维医学图像配准系统，包括：
[0041]
预处理模块，获取待配准图像，并对待配准图像作标准化预处理；
[0042]
建模模块，类比于复数空间的beltrami系数，于三维空间内定义一新测度，并基于新测度构建配准模型；
[0043]
优化模块，对配准模型作离散化处理，并使用高斯-牛顿迭代法对离散化处理后的配准模型优化，以得到最优变换模型；
[0044]
融合模块，基于最优变换模型构建待配准图像中的浮动图像，并将形变后的浮动图像与待配准图像中的目标图像融合显示。
[0045]
本发明又公开了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现如上所述的步骤。
[0046]
采用了上述技术方案后，与现有技术相比，具有以下有益效果：
[0047]
1、通过借助于模仿复数空间beltrami系数的定义来实现保拓扑结构，可以自动实现约束所有点的变换雅可比行列式大于0；
[0048]
2、借助于改进的高斯-牛顿方法，可以快速、稳定地计算最优解；
[0049]
3、所设计的归一化梯度向量场相似性测度适合多模态图像之间的配准，提升配准的精确度。
附图说明
[0050]
图1为符合本发明一优选实施例中三维医学图像配准方法的流程示意图；
[0051]
图2为符合本发明一优选实施例中四期ct图像的示意图。
具体实施方式
[0052]
以下结合附图与具体实施例进一步阐述本发明的优点。
[0053]
这里将详细地对示例性实施例进行说明，其示例表示在附图中。下面的描述涉及附图时，除非另有表示，不同附图中的相同数字表示相同或相似的要素。以下示例性实施例中所描述的实施方式并不代表与本公开相一致的所有实施方式。相反，它们仅是与如所附权利要求书中所详述的、本公开的一些方面相一致的装置和方法的例子。
[0054]
在本公开使用的术语是仅仅出于描述特定实施例的目的，而非旨在限制本公开。在本公开和所附权利要求书中所使用的单数形式的“一种”、“所述”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义。还应当理解，本文中使用的术语“和/或”是指并包含一个或多个相关联的列出项目的任何或所有可能组合。
[0055]
应当理解，尽管在本公开可能采用术语第一、第二、第三等来描述各种信息，但这些信息不应限于这些术语。这些术语仅用来将同一类型的信息彼此区分开。例如，在不脱离本公开范围的情况下，第一信息也可以被称为第二信息，类似地，第二信息也可以被称为第一信息。取决于语境，如在此所使用的词语“如果”可以被解释成为“在
……
时”或“当
……
时”或“响应于确定”。
[0056]
在本发明的描述中，需要理解的是，术语“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
[0057]
在本发明的描述中，除非另有规定和限定，需要说明的是，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是机械连接或电连接，也可以是两个元件内部的连通，可以是直接相连，也可以通过中间媒介间接相连，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。
[0058]
在后续的描述中，使用用于表示元件的诸如“模块”、“部件”或“单元”的后缀仅为了有利于本发明的说明，其本身并没有特定的意义。因此，“模块”与“部件”可以混合地使用。
[0059]
参阅图1，示出了符合本发明一优选实施例中三维医学图像的配准方法流程示意图，在该实施例中，三维图像配准方法包括以下步骤：
[0060]
s100：获取待配准图像，并对待配准图像作标准化预处理
[0061]
所获取的待配准图像可如图2所示，呈不同时期下病灶的ct图像。对这些待配准图像预处理后，待配准图像将具有相同的灰度和尺寸。
[0062]
s200：类比于复数空间的beltrami系数，于三维空间内定义一新测度，并基于所述新测度构建配准模型
[0063]
在复变函数空间中，给定一个函数μ，则而说明映射f是方向保持的，从而是拓扑保持的。其中μ:ω
→
d，f:ω
→
c，满足方程：
[0064][0065]
其中ω和d都是复数空间的子集，μ称为f的beltrami系数，对应的方程称为beltrami方程。然而beltrami系数在三维空间空没有相应的定义。因此，模仿它在二维空间中的定义，重新定义一个类beltrami系数的测度作为新测度，并以该新测度为正则项构建配准模型。
[0066]
s300：对所述配准模型作离散化处理，并使用高斯-牛顿迭代法对离散化处理后的配准模型优化，以得到最优变换模型
[0067]
为确定该配准模型的最优匹配效果，步骤s300中对配准模型作离散化处理，从而确定最小位移下(或称为最细尺度下)浮动图像与目标图像间的形变位移场。
[0068]
s400：基于所述最优变换模型构建待配准图像中的浮动图像，并将形变后的浮动图像与待配准图像中的目标图像融合显示。
[0069]
具有上述配置下，通过借助于类比复数空间beltrami系数的定义来实现保拓扑结构，可以自动实现约束所有点的变换雅可比行列式大于0。
[0070]
一优选实施例中，步骤s100包括：
[0071]
s110：搜集医学图像形成具有多个待配准图像的待配准图像集，并对每一待配准图像中的目标组织分割，形成分割图像；
[0072]
例如可从医院放射科收集待配准的图像，并从待配准图像中分割而出仅包含目标组织的部分，减少待配准区域。
[0073]
s120：对分割图像作灰度归一化处理到[0，1]，例如：mri和ct、基于窗宽窗位作归一化处理，并对归一化处理后的分割图像作数字图像化重采样，使得分割图像中的分辨率是各向同性的，各向同性指物体的物理、化学等方面的性质不会因方向的不同而有所变化的特性，即某一物体在不同的方向所测得的性能数值完全相同，亦称均质性；
[0074]
s130：对分割图像裁剪，使得每一对待配准图像的尺寸一致
[0075]
例如，将所有待配准图像的尺寸统一为256x256x128或128x128x128或224
×
224
×
128。再者可构建图像金字塔，从1到l，使得待配准图像的分辨率从粗精度到细精度。
[0076]
进一步优选或可选地，步骤s200包括：
[0077]
s210：类比于对复数空间的beltrami系数的定义，于三维空间内定义如下新测度：
[0078]
其中表示雅可比行列式，sign()表示符号函数，且所述新测度的模小于1等价于变换的雅可比行列式大于0；
[0079]
s220：以所述新测度作为正则项，构建如下配准模型：
[0080][0081]
其中f:表示待配准图像中的目标图像，m:表示待配准图像中的浮动图像，y:r3→
r3表示空间变换，ω表示待配准图像的定义域，min
y l(y)表示最小化能量泛函，最小化能量泛函中第一项和第二项合在一起表示数据保真项，化能量泛函，最小化能量泛函中第一项和第二项合在一起表示数据保真项，表示归一化梯度算子，表示正常梯度算子，α1和α2表示光滑约束参数，β表示雅可比约束项的参数。
[0082]
可以理解的是，当假设有f:表示固定图像(目标图像)，m:表示浮动图像(源图像)，y:r3→
r3表示空间变换，其中ω表示图像的定义域时，则基于变分模型配准的过程是优化如下能量泛函的过程：对于能量泛函的优化，便可得到满足条件的左右变换，其中s()表示数据保真项，即相似性测度，用于估计浮动图像和固定图像之间的距离(该值越小越好)，而r()表示正则项，即形变约束项，约束最终的变换满足一定的性质。由于医学图像的特殊性，对形变有更严格的要求。例如，变换应该是一对一的，必须保证浮动图像中的解剖点在固定图像中有且仅有一个对应的解剖位置。其次，形变之后的浮动图像与固定图像要尽可能的相似，也就是保真项的值越小越好。这两个要求在一般图像配准中也必须满足。但是医学图像显示的是一个解剖组织，在形变过程中，应该保证涉及到的解剖组织的拓扑结构是一致的，否则就会出现错误的匹配，造成信息的错误融合。因此，在医学图像配准中，形变网格不能出现折叠，即保拓扑结构，从而在步骤s220中新测度具有很好的性质，例如基于这个性质，可以通过约束n(y)＜1来达到约束变换是拓扑保持的。
[0083]
进一步优选实施例中，步骤s230包括：
[0084]
s231：当待配准图像的模态为单模态时，数据保真项选取如下平方差之和ssd：ssd(y)＝∫
ω
(m
°
y-f)2dx；
[0085]
s232：当待配准图像的模态为多模态时，数据保真项选取如下梯度信息(呈归一化表达)：
[0086][0087]
可以理解的是，选定了数据保真项，就完成了配准模型的构建，即能量泛函的构建。
[0088]
一优选实施例中，步骤s300包括：
[0089]
s310：基于中值公式对配准模型中的积分进行近似，并基于差分格式对配准模型中的导数进行离散，得到离散能量泛函(以多模态为例)；
[0090]
s320：对所述离散能量泛函求解出如下对应梯度：
[0091][0092]
其中φ1(rc)＝(φ1((rc)1)，
…
，φ1((rc)n))，φ(r)＝(φ((r)1)，
…
，φ((r)n))，其中a和b分别对应一阶和二阶求导矩阵，un为离散后的位移场，h为离散单位，r(un)表示n(y)关于un离散的向量形式，rc(un)通过化简gf_tm得到；表示r(un)关于un的梯度向量。表示φ关于r(un)的梯度向量，表示φ1关于rc(un)的梯度向量；
[0093]
s330：求解配准模型对应的如下海森矩阵：
[0094]
其中表示由φ1关于rc(un)的二阶导构成的对角矩阵，表示由φ关于r(un)的二阶导构成的对角矩阵。海森矩阵(hessian matrix或hessian)是一个多变量实值函数的二阶偏导数组成的方块矩阵；
[0095]
s340：对所述海森矩阵求解可行方向，即求解对应的线性方程，并且进行线性迭代求解满足需求的位移：
[0096][0097]uk+1
＝uk+θkδuk[0098]yk+1
＝xn+u
k+1
，其中δu表示迭代方向，θ表示迭代步长，k表示第k次迭代，xn表示形变图像的原始坐标；
[0099]
迭代后得到最细尺度对应的最优变换模型。
[0100]
本发明还公开了一种三维医学图像配准系统，包括：预处理模块，获取待配准图像，并对待配准图像作标准化预处理；建模模块，类比于复数空间的beltrami系数，于三维空间内定义一新测度，并基于新测度构建配准模型；优化模块，对配准模型作离散化处理，并使用高斯-牛顿迭代法对离散化处理后的配准模型优化，以得到最优变换模型；融合模块，基于最优变换模型构建待配准图像中的浮动图像，并将形变后的浮动图像与待配准图像中的目标图像融合显示。
[0101]
本发明又公开了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现如上所述的步骤。
[0102]
应当注意的是，本发明的实施例有较佳的实施性，且并非对本发明作任何形式的限制，任何熟悉该领域的技术人员可能利用上述揭示的技术内容变更或修饰为等同的有效实施例，但凡未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何修改或等同变化及修饰，均仍属于本发明技术方案的范围内。