CN104835112A - 一种肝脏多相期ct图像融合方法 - Google Patents

Abstract

本发明公开了一种肝脏多相期CT图像融合方法，首先利用基于联合直方图的多分辨率CT图像配准方法对源图像序列进行粗配准，接着结合置信连接的区域生长算法实现了肝脏的图像自动分割和基于梯度向量流snake模型的肝脏图像分割，有效的提取肝脏的边缘信息；再对肝脏图像进行基于定向区域生长算法的血管提取，接着对肝实质图像进行基于B样条自由形变变换和基于空间加权互信息的肝脏非刚性配准，精确找到空间同一位置的图像对；最后基于小波变换进行图像融合。本发明针对肝脏CT图像的特点，将图像分割和图像配准过程结合到图像的融合过程中，从而大大提高了融合的精度。

Description

一种肝脏多相期CT图像融合方法

技术领域

本发明涉及医学图像融合技术领域，特别涉及一种肝脏多相期CT图像融合方法。

背景技术

图像融合是指将多源信道或多次采集到的关于同一目标的图像数据经过图像处理，最大限度的提取各自的有利信息，最后综合成高质量的图像，从而利于后续分析。在医学图像处理领域中，为获得人体某一部位更全面、更清晰、特征更显著的图像往往采用图像融合的方法来提高图像质量。这一技术将几组图像放在一起分析，将细节信息不同的图像综合成一副细节信息完整的图像，并保持与原始图像近似的细节部分，从而为后续的医学诊断提供帮助。进行图像融合时，待融合图像需要在空间上已经对齐，即已配准好，配准效果的好坏直接影响融合的效果。

肝脏CT增强扫描指的是往患者体内注射造影剂，造影剂随血流进入肝脏，在不同的时间段，对患者进行CT扫描。一般将肝脏增强扫描分为动脉期、门脉期及延迟期。CT扫描时，检查者呈仰卧位并屏气不动，操作者按照一定的流速，将一定剂量的造影剂注射到患者体内，给药后一段时间先进行全肝动脉期扫描；再间隔一段时间后进行全肝门脉期扫描；最后再间隔一定的时间对病灶区行延迟扫描。在不同时间进行扫描获得的多组的图像称为多相期图像，不同相期的图像是造影剂在流经肝脏内部不同的血管时CT扫描所获取的图像，具有不同的特点。动脉期图像是造影剂流经肝脏动脉时获取的图像，因此动脉期图像中肝动脉呈现高亮度显示，容易从图像中区分出来，但肝脏内部的其余管道，如肝静脉、门静脉等管道中不存在造影剂，因此这些管道与肝实质在CT数值上较为接近，不易从图像中区分出来。同理，门脉期图像中，容易将肝静脉、门静脉管道从图像中区分出来。

在肝病的诊断中，往往需要同时观察到同一位置肝内不同的组织结构以及它们之间的位置关系，给诊断提供依据。此时就需要把这些组织结构在同一幅图像中显示出来，并进行对比和分析。但如前所述，不同组织结构的图像是在不同相期扫描获得的，它们是在不同的图像中显示各自的特征，要把这些特征在同一张图像中显示出来，这就需要进行图像融合，即将图像细节信息不同的图像综合出一副细节信息完整的图像，并保持与原始图像近似的细节部分。

然而，由于CT扫描时其起始位置未必完全一致，以及患者在扫描过程中可能存在呼吸运动或其他轻微的身体位置的移动，肝脏不同相期的扫描图像序列，它们在空间上的位置关系并不是一一对应的，不同相期间的图像并不是完全对齐的。例如动脉期的第50张图像未必就对应着门脉期的第50张图像，它们可能表示肝脏的不同位置。目前，医生只能根据主观经验来选择、匹配同一位置的不同图像，这就给诊断分析带来了困难和不便。因此在对肝脏CT增强扫描多相期图像融合之前，首先需要将这些图像进行对齐，即需要对肝脏CT增强扫描多相期图像进行图像配准。配准结果的好坏直接影响图像融合的质量，进而甚至可能对医生的诊断带来影响。

由于肝脏CT增强扫描多相期图像中还存在胃、脊柱、肾脏等其他组织器官，图像中存在这些组织器官可能对肝脏图像配准带来误差甚至引起错误。因此，一种较为精确的配准方式，在其处理时，应将图像中的肝脏器官进行图像分割后，提取出相同的部分后再进行图像配准。

此外，由于人的腹腔没有骨骼包围(上腹部有部分肋骨包围)，是柔软且极易发生形变的。腹腔除了本身柔软的特点之外，还受到呼吸运动的影响。在增强扫描的过程中，病人需要躺在一个设备上大约5分钟，因此可以知道在扫描的过程中病人的正常呼吸也在进行，这将会导致胸腔体积的连续变化并使胸腔发生形变。刚体变化、仿射变换的一些传统线性变换往往作用于图像全局或者头部等有硬骨骼包被的主支气管的配准。对于腹部图像，这样传统的全局变换往往不能进行精确的局部建模，或者使用这些全局变换无法建模一些局部的形变，在实际的临床应用中，还需要更精确的建模方法来实现图像的配准。同时，由于非线性形变的复杂性，在非刚性配准领域，基于互信息的方法很难满足要求，对基于互信息的方法进行改进依然是一个值得深入研究的课题。本发明首先进行全局刚性变换配准方法对CT图像序列进行粗配准，经过对图像中的肝脏器官图像分割后，提取出肝实质部分，再采用基于B样条的自由形变变换进行精确的局部建模和基于空间加权互信息的配准方法实现图像的非刚性配准。

现有学者对肝脏CT增强扫描多相期图像融合方法研究主要有：厦门大学黄晓阳《在体肝脏图像配准方法及应用研究》以及厦门大学王振环《肝脏CT图像配准、融合技术的研究与实现》，但这些方法中都未进行分割图像，将肝脏器官提取后再进行配准，这有可能给图像配准带来误差。同时在非刚性配准的研究中，引入空间信息是现阶段基于像素强度医学图像配准领域的研究方向之一，Studholme等将空间位置信息作为与参考图像和浮动图像地位等同的第三维信息引入到图像配准中，提出基于区域互信息RMI的配准方法，Loeckx等对Studholme的方案进行了改进，提出了条件互信息cMI的概念。本发明采用空间加权互信息方案，结合基于B样条的自由变换模型对图像进行空间分块处理，计算形变模型B样条曲线对不同位置像素点的贡献加权，该方案配准精度明显优于传统的互信息方法。此外，配准框架中常用的鲍威尔算法是一种局部的确定性优化算法，收敛速度快。由于互信息函数通常不光滑，而鲍威尔算法的寻优结果与初始点的选择有关，所以鲍威尔算法容易陷入局部极值造成误配准。周其慧等人提出了一种改进的鲍威尔优化算法与多分辨率策略相结合的配准方法，降低配准算法陷入局部极值的可能。赵海峰等人提出了一种将改进的人工鱼群算法与鲍威尔算法相结合的配准方法，从而避免鲍威尔算法寻优过程中陷入局部极值的困境。基于上述考虑，本发明在刚性配准中结合了多种优化算法的特点，不仅将高斯多分辨率策略有效整合到配准算法中，并在不同分辨率层使用了不同优化算法，有效克服配准过程中局部极值的影响造成误配的问题，具有一定的实用价值。而且在刚性配准中采用基于联合直方图的相似性测度，改进了常用的互信息测度，大大提高了配准速度。但这些配准方法都未能针对肝脏CT图像的特点展开研究。

发明内容

为了解决上述问题，本发明公开了一种肝脏多相期CT图像融合方法，首先利用基于联合直方图的多分辨率CT图像配准方法对源图像序列进行粗配准，接着结合置信连接的区域生长算法实现了肝脏的图像自动分割和基于梯度向量流snake模型的肝脏图像分割，有效的提取肝脏的边缘信息；再对肝脏图像进行基于定向区域生长算法的血管提取，接着对肝实质图像进行基于B样条自由形变变换和基于空间加权互信息的肝脏非刚性配准，精确找到空间同一位置的图像对；最后基于小波变换进行图像融合。本发明针对肝脏CT图像的特点将图像分割和图像配准过程结合到图像的融合过程中，从而大大提高了融合的精度。

为达到上述目的，本发明的技术方案为：

一种肝脏多相期CT图像融合方法，包括如下步骤：

步骤一、输入动脉期图像序列A和静脉期图像序列B；

步骤二、利用基于联合直方图的多分辨率CT图像配准方法进行粗配准，对于图像序列A中每一张CT图像，在图像序列B中找到空间位置相同的CT图像，即得到与源图像序列A对应的图像序列C；

步骤三、对图像序列A进行基于置信连接的自动肝脏图像分割，得到分割后肝脏图像序列D；

步骤四、利用步骤三分割后所得肝脏的轮廓线作为图像序列C的初始轮廓；对图像序列C进行基于梯度向量流snake模型的肝脏图像分割，得到分割后图像序列E；

步骤五、分别对肝脏图像序列D，图像序列E进行基于定向区域生长算法的血管提取；

5a)对肝脏图像灰度分析选取血管根部的种子点；

5b)使用定向区域生长算法对血管进行分割；

5c)输出经血管提取后的肝实质图像序列F和肝实质图像序列G；

步骤六、利用基于B样条的自由形变变换和基于空间加权互信息的非刚性配准方法进行细配准，对上面步骤分割提取后的肝实质进行处理，即对于肝实质图像序列F中每一张CT图像，在肝实质图像序列G中找到空间位置相同的CT图像；

步骤七、基于小波变换进行图像融合。

进一步的，所述步骤二的详细步骤为：

2a)输入待配准的两幅源图像，分别记为参考图像和浮动图像；

2b)通过金字塔算法对待参考图和浮动图进行多分辨率处理，将图像分为三层，构成图像金字塔，顶层图像分辨率最小，最底层为原始图像；

2c)对参考图金字塔顶层图像的指定区域进行几何坐标刚体变换得到新的区域；

2d)通过三次B样条插值方法得到浮动图金字塔顶层图像在2c)新区域的坐标；

2e)利用基于联合直方图的相似性测度计算参考图金字塔顶层图像和插值图之间的相似度，得到一个关于几何变换参数的函数；

2f)相似度函数输入到粒子群算法优化算法中进行最优化计算得到最优变换参数，这个过程通过迭代来实现，即重复2c)～2e)步直到取得最大值；

2g)将通过粒子群算法搜索的变换参数作为下一阶段优化算法搜索的起始点，对图像金字塔上的最后两层图像进行处理，原理同2c)～2f),优化算法使用了收敛速度更快的改进鲍威尔算法。

2h)输出浮动图在最终变换下的配准图像。经粗配准，即得到与源图像序列A对应的图像序列C。

进一步的，所述步骤2e)中，对两幅图像的联合直方图的定义为：两幅图像中对应像素灰度值的统计概率分布，对于离散图像，其概率密度为：

P(m,n)＝N(m,n)/M

其中，M为一幅图像所包含的像素总数；N(m,n)为两幅图像中对应像素灰度值分别为m和n的像素总数，以一幅图像的灰度值为横坐标，另一幅图像的灰度值为纵坐标，画出的概率分布图称为联合直方图。

基于联合直方图的配准相似性测度如下：

$N(\mathrm{tx},\mathrm{ty},\mathrm{tz},\mathrm{\&theta;})=\frac{\underset{|i-j|<h}{\mathrm{\&Sigma;}}{p}^{\mathrm{XY}}(i,j)}{\underset{i,j}{\mathrm{\&Sigma;}}{p}^{\mathrm{XY}}(i,j)}$

其中阈值h取图像中像素取值范围的10％能获得较好的效果。

在实现时，计算阈值区域内的点的个数，在求解时，就能够计算相似函数中的值。

进一步的，所述步骤2f)中，设定一个搜索问题，规模为R的粒子群在D维空间优化，粒子i(1≤i≤R)在第n(1≤n≤Nmax)代的速度表示为Vi(n)＝(vi1,vi2,…,viD)，位置表示为Xi(n)＝(xi1,xi2,…,xiD)；为了保持每个粒子的速度惯性，引入了惯性常数ω这个量；ω值较小时，能使得粒子最终收敛到最佳位置；ω值较大时，能使得粒子在全局范围内搜索能力提高；改进的粒子速度更新公式如下所示：

v_ij(n)＝ωv_ij(n-1)+α₁β₁[p_ij-x_ij(n-1)]+α₂β₂[p_uj-x_ij(n-1)]

其中n为迭代次数，α1、α2为加速常数，β1、β2的取值为[0，1]范围内的随机数。p_ij为粒子i的历史最佳位置在j维上的值，而p_uj为某一局部区域内所有粒子的历史最佳位置在j维上的值。若该区域所包含的粒子为整个群体，则p_u代表的是整个粒子群最佳历史位置,否则p_u代表局部的最佳位置。当算法迭代次数达到最大值时或者当粒子搜索的最佳历史位置低于设定的阈值以及在参数空间中粒子之间的距离接近设定的阈值时，算法将终止迭代。

进一步的，所述步骤三的详细步骤为：

3a)利用改进的曲线各向异性扩散滤波对CT图像进行平滑预处理；

3b)对预处理图像的灰度分析自动选取序列种子点；

3c)使用置信连接的区域生长算法对肝脏进行分割；

3d)利用空洞填充法改善分割结果，得到源图像序列A分割后肝脏图像序列D。

进一步的，所述步骤四的详细步骤为：

4a)利用步骤三分割后所得肝脏的轮廓线作为图像序列C的初始轮廓；

4b)得到需分割图像的边界图f(x,y)；

4c)得到边界图f(x,y)的梯度向量流；

4d)得到轮廓线的内部能量，即弹性能量和弯曲能量；

4e)综合所有能量进行曲线的变化，并得到轮廓线的最终收敛结果；即得到图像序列C分割后肝脏图像序列E。

进一步的，所述步骤六的详细步骤为：

6a)输入待配准的两幅肝实质图像，分别记为参考图像和浮动图像；

6b)通过自由形变变换对参考图的肝实质进行几何建模，并采用B样条插值方法得到浮动图网格新坐标；

6c)利用空间加权方式改进的互信息作为相似性测度计算参考图和插值图之间的相似度，得到一个关于几何变换参数的函数；

6d)将相似度函数输入到有限存储LBFGS搜索算法优化器中进行最优化计算得到最终变换参数；

6e)输出浮动图在最优变换下的配准图像。

进一步的，所述步骤七的详细步骤为：

7a)分别对步骤六中的参考图和配准图像进行小波分解，得到参考图和配准图像的不同频带系数；

7b)选择低频系数融合规则进行融合处理；

7c)对融合后所得系数进行小波逆变换，所得的重构图像即为融合图像。

相对于现有技术，本发明的有益效果为：

本发明首先利用基于联合直方图的多分辨率CT图像配准方法对源图像序列进行粗配准，接着结合置信连接的区域生长算法实现了肝脏的全自动分割和基于梯度向量流snake模型的肝脏分割，有效的提取肝脏的边缘信息；再对肝脏图像进行基于定向区域生长算法的血管提取，接着对肝实质图像进行基于B样条自由形变变换和基于空间加权互信息的肝脏非刚性配准，精确找到空间同一位置的图像对；最后基于小波变换进行图像融合。对肝脏多相期图像融合前，首先在刚性配准中结合了多种优化算法的特点，在高斯多分辨率策略不同分辨率层使用了不同优化算法，有效克服配准过程中局部极值的影响造成误配的问题。而且采用基于联合直方图的相似性测度，改进了常用的互信息测度，大大提高了配准速度。然后结合置信连接的区域生长算法和基于梯度向量流snake模型进行肝脏分割，有效的提取肝脏的信息，并在此基础上进行血管提取，最终对肝实质进行配准，采用该方法最终是对多相期图像中的共同部分(肝实质部分)进行配准，从而大大提高了配准的精度，提高了图像融合的精度。

本发明的图像分割、配准以及融合几个步骤在功能上彼此相互支持，粗配准的目的是为了提高分割的效果，分割的目的是为了提高后续精细配准的精度，精细配准是为了提高融合的质量。即步骤二的粗配准，其目的在于步骤4a)，步骤三、四、五的目的在于提高步骤六的精度，经过步骤六的处理，目的在于提高了步骤七的融合效果。这些步骤的组合取得了新的技术效果，组合后的技术效果比每个技术特征效果的总和更优越。

附图说明

图1为本发明的原理流程图；

图2为本发明实施例的肝动脉期图像；

图3为本发明实施例的肝静脉期图像；

图4为本发明肝脏图像配准的流程图；

图5为本发明肝脏图像配准中两幅图像的联合直方图；

其中a)两幅相同图像的联合直方图；

b)平移3象素时的联合直方图；

c)平移9象素时的联合直方图；

图6为本发明基于置信连接的全自动肝脏分割的流程图；

图7为本发明基于梯度向量流snake模型的肝脏分割的流程图；

图8为本发明肝脏图像融合的流程图。

具体实施方式

下面结合附图和具体实施方式对本发明技术方案做进一步详细描述：

如图1所示：一种肝脏多相期CT图像融合方法，包括如下步骤：

步骤1，输入动脉期图像序列A和静脉期图像序列B；

1a)将待融合的动脉期图像序列A和静脉期图像序列B分别存入计算机硬盘；

1b)分别读取步骤1a)中存储在计算机硬盘空间中待融合的动脉期图像序列A和静脉期图像序列B。

步骤2，利用基于联合直方图的多分辨率CT图像配准方法进行粗配准，对于动脉期图像序列A中每一张CT图像，在静脉期图像序列B中找到空间位置相同的CT图像。具体步骤如下：

2a)输入待配准的两幅图像，分别记为参考图像和浮动图像。

2b)通过金字塔算法对待参考图和浮动图进行多分辨率处理，将图像分为三层，构成图像金字塔，顶层图像分辨率最小，最底层为原始图像。多分辨率金字塔算法通常对图像进行平滑以及下采样处理，经过金字塔算法处理后，形成了不同空间分辨率的序列图像。这些图像自上而下分辨率逐渐变大，最底层是原始图像。通常使用的金字塔算法有：高斯金字塔、拉普拉斯金字塔、小波变换金字塔等。在本专利中，使用高斯金字塔算法对图像进行多分辨率处理,G_K层图像是通过G_K-1层图像与高斯核卷积所得，对图像处理的计算公式如下所示：

$G_K(i,j)=\underset{m}{\mathrm{\&Sigma;}}\underset{n}{\mathrm{\&Sigma;}}W(m,n)G_{K-1}(2i+m,2j+n);1\&le;K\&le;N$

其中,G_k(i,j)为第k层高斯金字塔图像,k为分解层数,w(m,n)＝h(m)h(n)为5＊5的低通特性窗口函数,h为高斯密度分布函数,i,j为图像的行列数。这样通过多分辨率金字塔算法对原始图像进行迭代计算可以得到其他层图像。

2c)对参考图金字塔顶层图像的指定区域进行几何坐标刚体变换得到新的区域；一个空间中的点(x,y,z)经过变换后的新坐标为(x*,y*,z*),用矩阵表示如下：

[x* y* z* 1]＝[x y z 1]·Rx·Ry·Rz·T

其中R_X，R_Y，R_Z分别表示绕X，Y，Z轴的旋转变换矩阵，T为平移变换矩阵，它们分别为：

$\mathrm{Rx}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \mathrm{rx}& \sin \mathrm{rx}& 0\\ 0& -\sin \mathrm{rx}& \cos \mathrm{rx}& 0\\ 0& 0& 0& 1\end{array}\right]$ $\mathrm{Ry}=\left[\begin{array}{cccc}\cos \mathrm{ry}& 0& -\sin \mathrm{ry}& 0\\ 0& 1& 0& 0\\ \sin \mathrm{ry}& 0& \cos \mathrm{ry}& 0\\ 0& 0& 0& 1\end{array}\right]$

$\mathrm{Rz}=\left[\begin{array}{cccc}\cos \mathrm{rz}& \sin \mathrm{rz}& 0& 0\\ -\sin \mathrm{rz}& \cos \mathrm{rz}& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$ $T=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\\ \mathrm{tx}& \mathrm{ty}& \mathrm{tz}& 1\end{array}\right]$

2d)通过三次B样条插值方法得到浮动图金字塔顶层图像在2c)新区域的坐标。B样条函数是分段n次多项式，它是一组节点序列{xi}确定的B样条基函数及顶点系列{vi}构成的特征多边形二者的线性组合。这里取n＝3，即三次B样条函数做为插值函数。在坐标区间[a,b]上，取分割a＝x₀≤x₁≤...≤x_n＝b为节点，构造B样条基函数为

其中B样条基函数N_i,k(x)的第一个下标i表示序号，第二个下标k表示基函数的次数。该递推公式表明，欲确定第i个k次B样条基，需要用到x_i,x_i+1...x_i+k+1个节点。我们称区间[x_i,x_i+k+1]为N_i,k(x)的支撑区间，也就是说仅在这个区间内的值不为零。设有一组节点序列{x_i}(i＝0,1,2...,n)，由其确定的B样条基函数N_i,k(x)，有一顶点系列{v_i}(i＝0,1,2...,n)构成的特征多边形，将N_i,k(x)与v_i组合，得到k次B样条曲线，其方程为

$r\left(x\right)=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{N}_{i,k}\left(x\right)V_i,a\&le;x\&le;b$

设V_i(i＝0，1，2，…，n)为B样条曲线的控制顶点，则三次均匀B样条曲线的矩阵表达式为

$P_i\left(u\right)=\frac{1}{6}\left[\begin{array}{cccc}1& u& u^2& u^3\end{array}\right]\left[\begin{array}{cccc}1& 4& 1& 0\\ -3& 0& 3& 0\\ 3& -6& 3& 0\\ -1& 3& -3& 1\end{array}\right]\left[\begin{array}{c}{V}_i\\ V_{i+1}\\ V_{i+2}\\ V_{i+3}\end{array}\right]0\&le;u\&le;1,i=\mathrm{0,1}...n-3$

通过插值方法获得浮动图在上一步新区域的坐标，得到插值图。

2e)利用基于联合直方图的相似性测度计算参考图金字塔顶层图像和插值图之间的相似度，是一个关于几何变换参数的函数。相似性测度定量化地衡量了两幅图像匹配的效果。

本发明对两幅图像的联合直方图的定义为：两幅图像中对应像素灰度值的统计概率分布。对于离散图像，其概率密度为：

P(m,n)＝N(m,n)/M

其中，M为一幅图像所包含的像素总数；N(m,n)为两幅图像中对应像素灰度值分别为m和n的像素总数。以一幅图像的灰度值为横坐标，另一幅图像的灰度值为纵坐标，画出的概率分布图称为联合直方图。

本发明使用的基于联合直方图的配准相似性测度如下：

$N(\mathrm{tx},\mathrm{ty},\mathrm{tz},\mathrm{\&theta;})=\frac{\underset{|i-j|<h}{\mathrm{\&Sigma;}}{p}^{\mathrm{XY}}(i,j)}{\underset{i,j}{\mathrm{\&Sigma;}}{p}^{\mathrm{XY}}(i,j)}$

其中阈值h取图像中像素取值范围的10％能获得较好的效果。

在实现时，计算阈值区域内的点的个数，在求解时，就可以计算相似函数中的值。使用这种方法来计算代价函数，其时间复杂度约为O(n²)，比使用最大互信息法来配准图像的时间复杂度降低了2×n倍。

2f)相似度函数输入到粒子群算法优化算法中进行最优化计算得到最优变换参数，这个过程通过迭代来实现，即重复2c)～2e)步直到取得最大值。粒子群算法优化算法的执行是通过维持一群粒子来遍历参数空间，寻找到最佳变换参数。在迭代搜索时，参数空间中的每个粒子都要不断的更新三个属性值：适应值、位置和速度。而一个粒子的当前位置和速度以及该粒子当前最佳位置和粒子群体中当前最佳位置决定着该粒子的下一位置。算法迭代结束后，群体中的最佳位置就是该迭代算法的最终的优化结果。设定一个搜索问题，规模为R的粒子群在D维空间优化，粒子i(1≤i≤R)在第n(1≤n≤Nmax)代的速度表示为Vi(n)＝(vi1,vi2,…,viD)，位置表示为Xi(n)＝(xi1,xi2,…,xiD)。为了保持每个粒子的速度惯性，引入了惯性常数ω这个量。ω值较小时，能使得粒子最终收敛到最佳位置；ω值较大时，能使得粒子在全局范围内搜索能力提高。改进的粒子速度更新公式如下所示：

v_ij(n)＝ωv_ij(n-1)+α₁β₁[p_ij-x_ij(n-1)]+α₂β₂[p_uj-x_ij(n-1)]

其中n为迭代次数，α1、α2为加速常数，β1、β2的取值为[0，1]范围内的随机数。p_ij为粒子i的历史最佳位置在j维上的值，而p_uj为某一局部区域内所有粒子的历史最佳位置在j维上的值。若该区域所包含的粒子为整个群体，则p_u代表的是整个粒子群最佳历史位置,否则p_u代表局部的最佳位置。当算法迭代次数达到最大值时或者当粒子搜索的最佳历史位置低于设定的阈值以及在参数空间中粒子之间非常接近时，算法将终止迭代。

2g)将通过粒子群算法搜索的变换参数作为下一阶段优化算法搜索的起始点，对图像金字塔上的最后两层图像进行处理，原理同2c)～2f),优化算法使用了收敛速度更快的改进鲍威尔算法。优化起始点为[X Y R S]＝[0 0 0 1]，4个参数的含义分别为X表示X轴方向的平移，Y表示Y轴方向的平移，R表示旋转角度，S表示缩放系数。4个参数的搜索步长对应为[1 1 1 0.05]，搜索的范围为X[-20 20]、Y[-20 20]、R[-20 20]、S[-1 1]，迭代精度为e≤0.0001。

Powell优化算法步骤如下：

(1)选取起始点x⁽⁰⁾和一组线性无关的向量e⁽ⁱ⁾(i＝1,2,...,N)，N为变量个数，此例子中N＝4。e⁽ⁱ⁾取N个坐标轴的方向，即e⁽¹⁾(1,0,...0)、e⁽²⁾(0,1,...0)……；

(2)从x⁽⁰⁾沿方向e⁽ⁱ⁾(i＝1,2,...,N)依次进行N次一维搜索，得到：

x⁽ⁱ⁾＝x^(i-1)+λ_i ^*e⁽ⁱ⁾，i＝1，2，…N

$f\left(x^{\left(i\right)}\right)=\underset{\mathrm{\&lambda;}}{\min }f(x^{(i-1)}+\mathrm{\&lambda;}{e}^{\left(i\right)}),i=\mathrm{1,2},...N$

在完成这N次一维搜索后，得到x^(N)。

(3)计算最快上升方向上函数的变化：

$\mathrm{DEL}=\underset{1\&le;i\&le;N}{\max }|f\left(x^{\left(i\right)}\right)-f\left(x^{\left(0\right)}\right)|=|f\left(x^{\left(\mathrm{IBIG}\right)}\right)-f\left(x^{\left(0\right)}\right)|;$

(4)引进方向

e≡x^(N)-x⁽⁰⁾，PET≡2x^(N)-x⁽⁰⁾；

(5)计算

f_E≡f(PET)≡f(2x^(N)-x⁽⁰⁾)

若

f_E≥f₀

或f_E≤f₀且

2(f₀-2f_N+f_E)[(f₀-f_N)-DEL]²≥(f₀-f_E)²DEL

则将x^(N)作为新的起点，沿上面一组旧方向e⁽ⁱ⁾(i＝1,2,...,N)重复上面步骤，即转到步骤(2)。

若(5)中的条件不满足，沿方向e≡x^(N)-x⁽⁰⁾以x^(N)作为起始点进行搜索得到目标函数在此方向上的极小值点P。将原来的方向e^(IBIG)去掉而保留其余原有的N－1个方向，加上方向e仍得到N个方向：e⁽¹⁾，e⁽²⁾，…e^(N)，以此时的P作为新起始点，重复上面步骤，即转步骤(2)。

结束的标志有两个：

若则停止计算，其中ε为相对误差精度。

若上面的过程进行到等于某指定的迭代最大次数时，停止迭代，表示迭代失败。

2h)输出浮动图在最优变换下的配准图像，利用上一步计算出的最优变换参数输出浮动图在最优变换下的配准图像。经粗配准，即得到与动脉期图像序列A对应的静脉期图像序列C。

步骤3，对动脉期图像序列A进行基于置信连接的全自动肝脏分割。具体步骤如下：

3a)利用改进的曲线各向异性扩散滤波对动脉期图像序列A中CT图像进行平滑去噪处理。改进的曲率各向异性扩散方程为：

$\frac{\&PartialD;f}{\&PartialD;t}=|\&dtri;f|\&dtri;\&CenterDot;c\left(\right|\&dtri;f\left|\right)\frac{\&dtri;f}{|\&dtri;f|};$

其中，f＝f(x,y,t)为t时刻扩散滤波后的图像，f＝f(x,y,0)＝I(x,y)是输入图像，为梯度算子，是单调递减电导函数，电导系数k控制边缘对比处理的灵敏性，改进的曲率电导系数

3b)对预处理图像的灰度分析自动选取序列种子点。首先对平滑去噪后的图像序列A选取第3*n/8(取整)幅CT图像I进行直方图统计,因为在此幅图像可以保证肝脏所属像素点灰度占最大比例。在此直方图中有属于肌肉组织和属于肝脏2个峰值，其中，n为图像序列中CT图像的幅数，较大的峰值h属于肝脏。然后在每幅CT图像上按x,y方向依次统计各行各列中灰度值属于[h-10,h+10]的像素点数，得到拥有最多该灰度范围体素点的行列号x_max和y_max。最后，判断点(x_max,y_max)的灰度值是否在肝脏所属灰度范围[h-10,h+10]内。若是，则将此点作为种子点；否则舍弃此点。根据相邻CT图像上的种子点进行分割。

3c)使用置信连接的区域生长算法对肝脏进行分割。对上步自动选取种子点后的每幅CT图像，首先，利用种子点周围的一个小邻域计算包含在此区域中的所有像素灰度的平均值和标准差。通过一个给定控制亮度范围大小的乘法因子l乘以标准差来定义相似灰度的范围,即I(X)∈[m-l_δ,m+l_δ]；其中X是图像I中的像素点，m和δ分别是当前区域灰度的平均值和标准差；相邻像素点的灰度值在这个范围内的将被包含到当前区域中。然后，对包含在当前区域内的所有像素点再次计算其灰度平均值和标准差，从而定义一个新的灰度范围，用来查看当前区域的邻域中灰度值落在这个范围内的像素点，并将其包含进当前区域。重复这个过程，直到没有新的像素加进来或者达到指定的最大迭代次数。

3d)利用空洞填充法改善上步的分割结果，得到平滑的肝脏轮廓，最终得到动脉期图像序列A分割后的动脉期肝脏图像序列D。

步骤4，对静脉期图像序列C进行基于梯度向量流snake模型的肝脏分割。具体步骤如下：

4a)利用步骤3分割后所得肝脏的轮廓线作为图像序列C的初始轮廓。因为经过粗配准，两图像序列在空间位置上已经一一对齐，所以直接复制步骤3所得的肝脏轮廓线作为对应静脉期图像序列C的初始轮廓。

4b)得到需分割图像的边界图f(x,y)。设灰度图像I(x,y)的一个边界图 $f(x,y)=-\&dtri;E_{\mathrm{ext}}(x,y).$

4c)得到边界图f(x,y)的梯度向量流。梯度向量流(Gradient VectorFlow，GVF)将图像的梯度矢量向外扩散，比图像力更加有序、更能体现物体边界的宏观走势。该算法设GVF定义为向量场V(x,y)＝(μ(x,y),v(x,y))并满足下式的最小值：

$\mathrm{\&epsiv;}=\&Integral;\&Integral;\mathrm{\&mu;}(u_x^2+u_y^2+v_x^2+v_y^2)+{|\&dtri;f|}^2{|V-\&dtri;f|}^2\mathrm{dxdy}$ 根据上式，通过变分法，解下列欧拉方程就可以得到梯度向量流场向量场V(x,y)：

$\mathrm{\&mu;}{\&dtri;}^{2}u-(u-f_x)(f_x^2+f_y^2)=0$

$\mathrm{\&mu;}{\&dtri;}^{2}u-(v-f_y)(f_x^2+f_y^2)=0$

计算u场和v场的二阶微分后，根据下面的迭代公式：

$u_{i,j}^{n+1}=(1-b_{i,j}\&times;\mathrm{\&Delta;t})\&times;u_{i,j}^n+r\&times;{\&dtri;}^{2}u(i,j)\&times;4+c_{i,j}^1\&times;\mathrm{\&Delta;t}$

$v_{i,j}^{n+1}=(1-b_{i,j}\&times;\mathrm{\&Delta;t})\&times;v_{i,j}^n+r\&times;{\&dtri;}^{2}v(i,j)\&times;4+c_{i,j}^2\&times;\mathrm{\&Delta;t}$

计算出u场和v场的值。

4d)得到轮廓线的内部能量，即弹性能量和弯曲能量。内部能量由弹性能量和弯曲能量组成，这两个能量都是由轮廓自身决定的，和图像的特征无关。弹性能量能抑制曲线的伸长，弯曲能量能抑制曲线弯曲，合称为内能。在没有外能的情况下，内能的作用是使曲线平滑地收缩。本发明使用的是Williams和Shah提出的内部能量的计算方法，该方法对于控制点邻域上的每个点计算内部能量的公式为：

$E_{\mathrm{Elastic}}\left[x\right]\left[y\right]=|\mathrm{avg}\_\mathrm{dis}\tan \mathrm{ce}-\sqrt{{(x-x_{k-1})}^2+{(y-y_{k-1})}^2}|$

$E_{\mathrm{Bending}}\left[x\right]\left[y\right]=\sqrt{{(x_{k+1}-2x+x_{k-1})}^2+{(y_{k+1}-2y+y_{k-1})}^2}$

avg_distance为轮廓线上所有控制点的平均距离，随着控制点位置的变化而变化。通过这样计算距离偏差的方法可以使得每个控制点获得不同的弹性能量，进而在迭代的过程中趋近于均匀分布，保证曲线的连续性；同时，这种在离散的状态下计算曲率的方法较为方便，能提高算法的执行效率。

4e)综合所有能量进行曲线的变化，并得到轮廓线的最终收敛结果。曲线的演化就是在内外力共同作用下，朝着目标边界收敛的过程。在每一次曲线变化时，都会先根据控制点之间的间距，动态的添加和删除控制点，以保证曲线的正常收敛。然后根据Williams和Shah提出的贪婪算法，运用每次变化状态下每种能量的值，在每个控制点的领域上寻找合适的更新点，将控制点移动过去，该算法是将原来仅对控制点的计算转换到控制点领域上的每一点，这样通过计算每个轮廓点局部的能量最小值，进而得到整个轮廓线的最小值，完成轮廓线的变化。曲线演化的方程为：

X_t(s,t)＝αX”(s,t)-βX”(s,t)-γV。

对于曲线演化的迭代，本发明采用的是3x3邻域，在控制点邻域上查找能量最小点并更新控制点。最终即得到静脉期图像序列C分割后静脉期肝脏图像序列E。

步骤5，分别对肝脏图像序列D，E进行基于定向区域生长算法的血管提取。具体步骤如下：

5a)对肝脏图像灰度分析选取血管根部的种子点。

5b)使用定向区域生长算法对血管进行分割。

算法的输入是灰度图像I(共有M个体素),和一个种子点S0；算法的输出是边的集合E(初始值为空集)。对于体素x∈I，N(x)表示x的邻点的集合,称N(x)中的每个体素和x相邻。对于体素集N(S)表示不属于S并和S中某一体素相邻的体素的集合:

$N\left(S\right)=\left\{x\right|x\&Element;I,x\&NotElement;S,$

使得x∈N(y)}

算法是一个迭代过程，每次迭代过程中，向已生长区域加入一个体素，向边集加入一条边。设第m次迭代后已生长区域为Sm，边集为Em。第m+1次迭代中，将N(Sm)中最亮的体素g加入Sm；记Sm中和g相邻的最亮体素为s,将边(s,g)加入Em:

g＝max(N(S_m))

s＝max(N(g)∩S_m)

S_m+1＝S_m∪g

E_m+1＝E_m∪(s，g)

经M-1次迭代后,所有体素全部加入已生长区域:SM-1＝I。与常规的区域生长法不同，本发明算法在生长过程中，每次只生长到Sm的邻域中最亮的一个点，即为生长过程规定了某种方向性，我们称之为“定向”区域生长。这种生长方式的特点是在越过图像中的低灰度区域后,可以迅速回到高灰度区域(即血管)继续生长,再通过后续处理就有可能将被低灰度区域分断开的血管末梢包含进血管分割结果。生长结束时整幅图像形成了一棵以种子点S0为根结点的树,边集E＝EM-1存储了树中所有的边。算法在实际使用时没有必要生长至全图,可设定一个阈值,若g的灰度低于该阈值则停止生长,以节省运算时间。

接下来利用边集E分割感兴趣区域(ROI)。此时需要设定3个阈值：两个灰度阈值T1、T2和一个长度阈值T3，其中T1<T2。对于图像中任意体素x，存在一条从S0至x的路径:(S0,x1),(x1,x2),…,(xm,x)。则x属于ROI的条件是:

x的灰度不低于T2；

x1,x2,…，xm的灰度均不低于T1；

路径中灰度低于T2的连续体素个数不超过T3，即若xi,xi+1,…,xi+k-1的灰度均低于T2,则k≤T3。

此外,若x属于ROI,则从S0至x的路径上所有体素均属于ROI相对于常规的区域生长法,本发明算法在ROI的判定条件上有所放松:在从S0到x的路径上如果有连续若干个体素灰度较低,只要这些体素个数不是太多,仍然认为x属于ROI这样,算法就有可能跨越微细管道结构中的低灰度区域,将管道末梢包含进分割结果。

5c)输出经血管提取后的动脉期肝实质图像序列F和静脉期肝实质图像序列G。

步骤6，利用基于B样条自由形变变换和基于空间加权互信息的非刚性配准方法进行细配准，对上面步骤分割提取后的肝实质进行处理，即对于肝实质图像序列F中每一张CT图像，在肝实质图像序列G中找到空间位置相同的CT图像。具体步骤如下：

6a)输入待配准的两幅肝实质图像，分别记为参考图像和浮动图像。

6b)通过自由形变变换对参考图的肝实质进行几何建模，并采用B样条插值方法得到浮动图在控制点网格的新坐标。自由形变变换是一种对固体进行几何建模的方法，其主要目标是对固体的形变进行模拟和建模。自由形变变换的其核心思想是:变形操作并不直接作用于物体，而是作用于物体所嵌入的变形空间，通过改变变形空间来使物体发生形变。也就是说，如果变形空间被改变了，则嵌入其中的物体自然随之改变。把自由形变变换应用在二维的图像变换时候，通过引入一个覆盖在二维图像上的控制点网格来完成对图像的形变，当调整控制点的位置时，图像的形状也会发生变化。自由形变变换中位于控制点网格上控制点的位置是由外界指定的，而控制点网格内部的非控制点位置的像素则根据一定邻域内控制点的位置，通过样条插值来决定。本发明的算法在设计的时候，自由形变变换插值时采用基于B样条插值方法。

通过插值方法获得浮动图在控制点网格的新坐标，得到插值图。

6c)利用空间加权方式改进的互信息作为相似性测度计算参考图和插值图之间的相似度，是一个关于几何变换参数的函数。相似性测度定量化地衡量了两幅图像匹配的效果。

互信息可以作为图像配准相似测度标准的依据是，当两幅图像完全对准时，其互信息值达到最大。在图像配准过程中，由于几何变换的存在，待配准图像的重叠面积可能会发生改变，而互信息的大小与图像重叠面积相关联，为了规避互信息值对图像重叠面积的依赖性，Studholme提出了使用归一化的互信息作为图像配准的相似性度量标准，可用公式表述为:

$I(R,F)=\frac{H\left(R\right)+H\left(F\right)}{H(R,F)}$

在前面介绍基于B样条的自由变换模型时我们注意到，每个控制点φ_i,jj对图像形变的影响都是局部性的，其范围仅限于以其为中心的4×4个网格以内，且影响程度由近及远逐渐减小。因此，本发明方案考虑将浮动图像进行空间分块处理，针对每一个控制点φ_i,j，将受其影响的最大范围，也就是一个以φ_i,j为中心面积为4Δ×4Δ的图像块单独处理，计算其与参考图像中对应区域的归一化互信息。对大小为n_x×n_y的控制点网格，共要将待配准图像分成n_x×n_y个相互重叠的图像块，分别计算对应的局部区域归一化互信息，总的互信息根据每个图像块面积占图像总面积的比重pi，j依概率加权得到。

$I(R,F)=\underset{i,j}{\mathrm{\&Sigma;}}{p}_{i,j}I(R_{i,j},F_{i,j})$

同时，考虑到控制点坐标变化对其周围图像形变的影响程度随距离衰减，且衰减程度受B样条基函数控制，因此在计算浮动图像上每一个图像块与对应参考图像块的互信息时，与当前控制点距离不同的点对待配准图像的直方图和联合直方图的贡献不等同考虑。本发明方案采用与形变模型相对应的二维3阶B样条曲线对不同位置像素点的贡献加权。

考虑参考点φ_i,j对应的图像块，设φ_i,j的坐标为(x_i，j，y_i，j)，对浮动图像上的强度为f的点(x_f，y_f)，经过函数T(x，y)变形后，坐标为(T(x_f)，T(y_f))，对应的参考图像上的点的强度值为r，则当前待配准浮动图像块与对应参考图像块的联合直方图可表示为:

$h(r,f;{\mathrm{\&phi;}}_{i,j})=\underset{x_f\&Element;F_{i,j}}{\mathrm{\&Sigma;}}w(T\left(x_f\right)-x_{i,j},T\left(y_f\right)-y_{i,j})$

其中w(x，y)为二维空间权函数，即本方案所选用的二维3阶B样条曲线，其一维形式在[－2，2]上定义为:

$w\left(u\right)=\left\{\begin{array}{cc}{(u+2)}^3/6& u\&Element;[-2,-1]\\ (4-6u^2-3u^3)/6& u\&Element;[-\mathrm{1,0}]\\ (4-6u^2+3u^3)/6& u\&Element;\left[\mathrm{0,1}\right)\\ {(2-u)}^3& u\&Element;\left[\mathrm{1,2}\right]\end{array}\right.;$

假设待配准图像总的像素点数为N，则当前局部互信息对整幅图像互信息的贡献率：

$p_{i,j}=\frac{\underset{r\&Element;R_{i,j}}{\mathrm{\&Sigma;}}\underset{f\&Element;F_{i,j}}{\mathrm{\&Sigma;}}h(r,f;{\mathrm{\&phi;}}_{i,j})}{N}.$

6d)相似度函数输入到L-BFGS(L-BFGS即有限内存BFGS，“limited memory BFGS”)搜索算法优化器中进行最优化计算得到最终变换参数，这个过程通过迭代来实现，即重复6b)～6c)步直到取得最大值。L-BFGS是一种拟牛顿的最优化算法，L-BFGS使用了BFGS算法来估计海森矩阵。

BFGS是从牛顿的最优化算法衍生出来，属于登山算法的一种，寻找梯度为零的驻点。牛顿法假定目标函数在局部可以通过二阶泰勒公式近似，然后通过目标函数的一阶和二阶导数来确定驻点。拟牛顿算法中目标函数二阶导数的汉森矩阵不用被计算，而是通过分析梯度向量被更新。拟牛顿算法其实是寻找多维目标函数一阶导数的根的割线法的推广。

在BFGS中，第k次迭代的搜索方向d_k通过解拟牛顿方程得到:

$B_k{d}_k=-\&dtri;f\left({\mathrm{\&mu;}}_k\right)$

其中B_k是对海森矩阵的估计，且在每次迭代中被更新，是目标函数在x_k的梯度。得到了方向d_k后，再通过沿着d_k方向的线性搜索得到下一次的参数值μ_k+1。之后，BFGS不需要直接计算在μ_k+1的完整的海森矩阵来得到B_k+1，而是通过把第k次迭代得到的B_k上与两个矩阵相加得到：

B_k+1＝B_k+U_k+V_k

矩阵U_k和V_k的秩都为1但是有不同的基，所以U_k和V_k构成一个秩为2的更新矩阵，使得BFGS算法相对于梯度下降算法在缩放问题上更为鲁棒。BFGS的算法过程如下所示。

假设初始参数值为初始海森矩阵为B₀，反复执行如下步骤指导收敛到问题的解:

l)通过解 $B_k{d}_k=-\&dtri;f\left({\mathrm{\&mu;}}_k\right)$ 得到方向向量d_k；

2)使用线性搜索沿着方向向量方向得到一个最优步长a_k，更新参数μ_k+1＝μ_k+a_kd_k；

3)计算 $y_k=\&dtri;f\left({\mathrm{\&mu;}}_{k+1}\right)-\&dtri;f\left({\mathrm{\&mu;}}_k\right);$

4)计算 $B_{k+1}=B_k+\frac{y_k{y}_k^T}{y_k^T{d}_k}-\frac{B_k{d}_k{\left(B_k{d}_k\right)}^T}{d_k^T{B}_k{d}_k}.$

其中f(μ_k)是需要被最优化的目标函数，收敛过程可以通过检查梯度的模进行。在具体使用中，B0可以被初始化为单位矩阵即B0＝I，因此第一步搜索相当于一次梯度下降，之后的迭代则因为B_k逐步逼近海森矩阵而越来越精确。

上述算法的第一步对d_k的求解需要计算B_k的逆矩阵，而且这个计算过程一般是通过应用Sherman-Morrison公式来替代上述第四步来直接进行的:

$B_{k+1}^{-1}=B_k^{-1}+\frac{\left(d_k{d}_k^T\right)(d_k^T{y}_k+y_k^T{B}_k^{-1}{y}_k)}{{\left(d_k^T{y}_k\right)}^2}-\frac{B_k^{-1}{y}_k{d}_k^T+d_k{y}_k^T{B}_k^{-1}}{d_k^T{y}_k};$

L-BFGS的实现类似于任何一个直接的拟牛顿算法，目标是求解方向向量d_k＝-H_kg_k，其中H_k从为海森矩阵的逆。具体的，对于给定的第k次迭代参数μ_k和为目标函数f(μ_k)的梯度，保存增量s_k＝μ_k-μ_k-1和y_k＝g_k-g_k-1定义 是第k次迭代开始时对海森矩阵的初始估计。然后，方向向量可以按如下方法计算:

q＝g_k

For i＝m-1...0

$a_i=r_i{s}_i^{T}q$

q＝q-a_iy_i

$z=H_k^{0}q$

For i＝0...m-1

$b=r_i{y}_i^{T}z$

z＝z+s_i(a_i-b_i)

H_kg_k＝z

余下的计算步骤可见BFGS算法。

6e)输出浮动图在最优变换下的配准图像。

步骤7，基于小波变换进行图像融合。具体步骤如下：

7a)分别对上步的参考图和配准图像进行小波分解，得到参考图和配准图像的不同频带系数。对参考图和配准图像进行N层小波分解，最终有(3N+1)个不同频带，其中包含3N个高频带和一个低频带。

7b)按照融合规则进行融合处理。具体的融合规则为：

(1)融合图像的低频部分，取参考图和配准图像分解后的加权平均，即

C_N,F＝(C_N,A+C_N,B)/2

其中，CN,A、CN,B分别表示参加融合的参考图和配准图像在小波分解尺度N上的低频分量，CN,F表示融合图像F在小波分解尺度N上的低频分量。

(2)在最高分解层上，比较参考图和配准图像的3个方向高频分量的小波系数，取绝对值大的小波系数作为融合图像F的小波系数，即

DiN,F＝DiN,A if|DiN,A|≥|DiN,B|

DiN,F＝DiN,B else

其中，DiN,A、DiN,B分别表示参加融合的参考图和配准图像在小波分解尺度N上i方向上的小波系数，DiN,F表示融合图像F在小波分解尺度N上i方向上的小波系数。

(3)在中间分解层上，理想像素为中心的局部区域(这里取3×3)的均值方差最大的参考图或配准图像的小波系数作为融合图像F对应的小波系数，即

Dij,F＝Dij,A if MSEA≥MSEB

Dij,F＝Dij,B else

其中分解尺度j取1到N-1；MSEA、MSEB分别表示参考图和配准图像在分解尺度上方向上对应局部区域上的方差。方差MSE定义为：

$\mathrm{MSE}=\frac{1}{M*N}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{(X_{i,j}-\stackrel{\&OverBar;}{X})}^2$

其中，M、N分别为局部区域的行数和列数(这里为3)；xi,j为当前局部区域内的一个像素的灰度值，x为当前局部区域像素灰度值的平均值。

7c)对融合后所得系数进行小波逆变换，所得的重构图像即为融合图像。

本发明的原理流程图如图1所示，包括有输入原始数据、CT图像粗配准，动脉期肝脏自动分割、静脉期snake分割、肝脏图像血管提取、肝实质图像细配准、肝图像小波变换、输出融合图像等几个步骤。

本发明肝脏图像动脉期图像如图2所示，肝脏图像静脉期图像如图3所示。本发明肝脏图像配准的流程图如图4所示，肝脏图像配准包括输入待配准的两幅图像(参考图和浮动图)，参考图几何空间刚体变换，浮动图三次B样条插值，相似性测度计算参考图和插值图之间的相似度，Powell优化算法最优化变换参数，输出最优变换下的插值图像。

本发明肝脏图像粗配准中两幅图像的联合直方图如图5所示。a)为两幅相同图像的联合直方图；b)为平移3象素时的联合直方图；c)为平移9象素时的联合直方图。

本发明基于置信连接的全自动肝脏分割的流程图如图6所示，肝脏自动分割依次进行滤波平滑预处理，自动获取种子点，置信连接的区域生长法，后处理。

本发明基于梯度向量流snake模型的肝脏分割的流程图如图7所示，肝脏分割依次进行图像灰度化处理、获取边界图、获取梯度向量流、获取轮廓线内部能量、综合所有能量进行曲线变化。

本发明肝脏图像融合的流程图如图8所示，肝脏图像融合包括图像小波分解，低频系数融合规则进行融合处理，小波逆变换重构图像。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何不经过创造性劳动想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书所限定的保护范围为准。

Claims (8)

1.一种肝脏多相期CT图像融合方法，其特征在于，包括如下步骤：

步骤一、输入肝脏动脉期图像序列A和肝脏静脉期图像序列B；

步骤二、利用基于联合直方图的多分辨率CT图像配准方法进行粗配准，对于图像序列A中每一张CT图像，在图像序列B中找到空间位置相同的CT图像；即得到与源图像序列A对应的图像序列C；

步骤三、对图像序列A进行基于置信连接的自动肝脏图像分割，得到分割后肝脏图像序列D；

步骤四、利用步骤三分割后所得肝脏的轮廓线作为图像序列C的初始轮廓；对图像序列C进行基于梯度向量流snake模型的肝脏分割，得到分割后图像序列E；

步骤五、分别对肝脏图像序列D，图像序列E进行基于定向区域生长算法的血管提取；

5a)对肝脏图像灰度分析选取血管根部的种子点；

5b)使用定向区域生长算法对血管进行分割；

5c)肝脏图像序列D经血管提取后得到肝实质图像序列F，肝脏图像序列E经血管提取后得到肝实质图像序列G，输出图像序列F和图像序列G；

步骤六、利用基于B样条的自由形变变换和基于空间加权互信息的非刚性配准方法进行细配准，对上面步骤分割提取后的肝实质进行处理，即对于肝实质图像序列F中每一张CT图像，在肝实质图像序列G中进行配准；

步骤七、基于小波变换进行图像融合。

2.根据权利要求1所述的方法，其特征在于，所述步骤二的详细步骤为：

2a)输入待配准的两幅源图像，分别记为参考图像和浮动图像；

2b)通过金字塔算法对待参考图和浮动图进行多分辨率处理，将图像分为三层，构成图像金字塔，顶层图像分辨率最小，最底层为原始图像；

2c)对参考图金字塔顶层图像的指定区域进行几何坐标刚体变换得到新的区域；

2d)通过三次B样条插值方法得到浮动图金字塔顶层图像在2c)新区域的坐标；

2e)利用基于联合直方图的相似性测度计算参考图金字塔顶层图像和插值图之间的相似度，得到一个关于几何变换参数的函数；

2f)将2e)计算的相似度函数输入到粒子群算法优化算法中进行最优化计算得到最优变换参数，这个过程通过迭代来实现，即重复2c)～2e)步直到取得最大值；

2g)将通过粒子群算法搜索的变换参数作为下一阶段优化算法搜索的起始点，对图像金字塔上的最后两层图像进行处理，原理同2c)～2f),优化算法使用了鲍威尔算法；

2h)输出浮动图在最终变换下的配准图像；经粗配准，即得到与源图像序列A对应的图像序列C。

3.根据权利要求2所述的方法，其特征在于，所述步骤2e)中，对两幅图像的联合直方图的定义为：两幅图像中对应像素灰度值的统计概率分布，对于离散图像，其概率密度为：

P(m,n)＝N(m,n)/M

其中，M为一幅图像所包含的像素总数；N(m,n)为两幅图像中对应像素灰度值分别为m和n的像素总数，以一幅图像的灰度值为横坐标，另一幅图像的灰度值为纵坐标，画出的概率分布图称为联合直方图。

基于联合直方图的配准相似性测度如下：

$N(\mathrm{tx},\mathrm{ty},\mathrm{tx},\mathrm{\&theta;})=\frac{\underset{|i-j|<h}{\mathrm{\&Sigma;}}{p}_{\mathrm{XY}}(i,j)}{\underset{i,j}{\mathrm{\&Sigma;}}{p}_{\mathrm{XY}}(i,j)}$

其中阈值h取图像中像素取值范围的10％能获得较好的效果，

在实现时，计算阈值区域内的点的个数，在求解时，就能够计算相似函数中的值。

4.根据权利要求2所述的方法，其特征在于，所述步骤2f)中，设定一个搜索问题，规模为R的粒子群在D维空间优化，粒子i(1≤i≤R)在第n(1≤n≤Nmax)代的速度表示为Vi(n)＝(vi1,vi2,…,viD)，位置表示为Xi(n)＝(xi1,xi2,…,xiD)；为了保持每个粒子的速度惯性，引入了惯性常数ω这个量；ω值较小时，能使得粒子最终收敛到最佳位置；ω值较大时，能使得粒子在全局范围内搜索能力提高；改进的粒子速度更新公式如下所示：

v_ij(n)＝ωv_ij(n-1)+α₁β₁[p_ij-x_ij(n-1)]+α_2β2[p_uj-x_ij(n-1)]

其中n为迭代次数，α1、α2为加速常数，β1、β2的取值为

[0，1]范围内的随机数，p_ij为粒子i的历史最佳位置在j维上的值，而p_uj为某一局部区域内所有粒子的历史最佳位置在j维上的值，若该区域所包含的粒子为整个群体，则p_u代表的是整个粒子群最佳历史位置,否则p_u代表局部的最佳位置。

当算法迭代次数达到最大值时或者当粒子搜索的最佳历史位置低于设定的阈值以及在参数空间中粒子之间的距离接近设定的阈值时，算法将终止迭代。

5.根据权利要求1所述的方法，其特征在于，所述步骤三的详细步骤为：

3a)利用改进的曲线各向异性扩散滤波对CT图像进行平滑预处理；

3b)对预处理图像的灰度分析自动选取序列种子点；

3c)使用置信连接的区域生长算法对肝脏进行分割；

3d)利用空洞填充法改善分割结果，得到源图像序列A分割后肝脏图像序列D。

6.根据权利要求1所述的方法，其特征在于，所述步骤四的详细步骤为：

4a)利用步骤三分割后所得肝脏的轮廓线作为图像序列C的初始轮廓；

4b)得到需分割图像的边界图f(x,y)；

4c)得到边界图f(x,y)的梯度向量流；

4d)得到轮廓线的内部能量，即弹性能量和弯曲能量；

4e)综合所有能量进行曲线的变化，并得到轮廓线的最终收敛结果；即得到图像序列C分割后肝脏图像序列E。

7.根据权利要求1所述的方法，其特征在于，所述步骤六的详细步骤为：

6a)输入待配准的两幅肝实质图像，分别记为参考图像和浮动图像；

6b)通过自由形变变换对参考图的肝实质进行几何建模，并采用B样条插值方法得到浮动图网格新坐标；

6c)利用空间加权方式改进的互信息作为相似性测度计算参考图和插值图之间的相似度，得到一个关于几何变换参数的函数；

6d)将6c计算得到的相似度函数输入到有限存储LBFGS搜索算法优化器中进行最优化计算得到最终变换参数；

6e)输出浮动图在最优变换下的配准图像。

8.根据权利要求5所述的方法，其特征在于，所述步骤七的详细步骤为：

7a)分别对步骤六中的参考图和配准图像进行小波分解，得到参考图和配准图像的不同频带系数；

7b)按照融合规则进行融合处理；

7c)对融合后所得系数进行小波逆变换，所得的重构图像即为融合图像。