多相期增强CT扫描是肝病诊断的常见检查方法。通过注射造影剂，扫描可以获得多个相期下的CT图像序列，从而得到肝脏动脉、静脉血管的解剖信息。在基于多相期增强CT影像的肝脏计算机辅助诊断、治疗和手术规划系统中，从复杂的组织结构中单独提取血管有重要的意义，它是肝血管结构3维可视化及病理分析的基础。特别是在局部肝切除以及肝脏移植等手术中，能否保证剩余组织的足够供血是决定手术成败的关键因素之一，需要准确的血管系统结构及分布信息作为支撑^[1]。

在多相期增强CT影像中，组织血管的成像质量非常容易受到造影剂分布的影响，此外，图像本身的分辨率、噪声、伪影等都会提升血管分割的难度，特别是肝脏的血管系统包含了复杂的动脉、静脉和肝门静脉系统，使得肝内血管分割十分具有挑战。

利用Hessian矩阵特征进行管状结构增强是目前最常用的肝内血管分割方法，这一类的方法应用于噪声较多的肝脏CT图像具有比较好的效果^[2]。Frangi等人^[3]将血管增强当成是管状结构的滤波处理。通过分析管状、球状和片状等不同结构下Hessian矩阵3个特征值之间的关系，定义了一个多参数的度量函数用于管状血管结构的增强，同时提出采用多层滤波的方式来处理不同尺度的血管。这种增强方法对管状血管几何结构具有比较好的增强效果，但是增强的血管树可能存在一些断落而使得拓扑结构受到影响，并且增强结果还需要经过进一步的分割标记才能提取出最终的血管树^[4]。Mirhassani等人^[5]提出采用高阶阈值的方法去除一部分非血管噪声，然后使用霍夫变换和形态学运算的方法对断裂血管进行修复；Truc等人^[6]则提出加入方向信息来提高噪声环境下细小血管的增强效果，同时还可以解决血管断裂的问题。

本文提出了一种全自动的肝内血管分割方法。首先，对造影增强的肝脏CT图像进行各向异性的滤波处理，消除噪声影响；然后，改进了Frangi等人^[3]提出的血管度量函数，结合血管的灰度分布信息设计了一个新的度量项，可以在抑制噪声的同时增强血管树的连续性，避免节点断裂；随后，提出了一种自适应的区域增长算法，自动选取种子点并且在区域增长过程能自适应的调整阈值，这种方法对增强结果具有比较好的分割效果；最后，提出通过“由粗到细”的策略对多层增强融合结果进行迭代式分割，将上层的分割结果作为下一层的先验信息，从而更好地调整各层自适应区域增长算法的阈值，最终实现肝内血管的全自动分割。

因为处理主要针对肝脏及其内部的血管，为了消除临近组织对于分割效果的影响，本文采取了半自动的分割算法来提取感兴趣范围，首先将包含了血管的肝脏部分标记出来。

为了减弱肝脏CT图像上的噪声对血管分割的影响，需要在分割之前对图像进行平滑处理。选择各向异性滤波作为图像的降噪预处理^[7]，即

$ \left\{ \begin{align} & {{\partial }_{t}}u=\text{div}\left[ g\left( {{\left| \nabla u \right|}^{2}} \right)\cdot \nabla u \right] \\ & u\left| _{t=0}={{u}_{0}} \right. \\ \end{align} \right. $

(1)

式中，u表示图像的像素值，$ \nabla u$表示图像的梯度，div表示图像的散度，t表示热扩散的时间，g则表示扩散系数。这种方法能够在去除图像噪声的同时避免对边界区域的模糊。

虽然基于汉森矩阵特征值的管状结构增强方法很多。但是Frangi等人^[3]所提出的算法是其中最被广泛使用的。该方法首先分析了Hessian矩阵特征值与图像结构之间的关系。Hessian矩阵描述了3维图像上某一点的二阶导数。如果用λ₁、λ₂、λ₃分别表示Hessian矩阵的3个特征值，并且假设|λ₁|≤|λ₂|≤|λ₃|, Frangi等人所提出的血管度量函数为

$ \begin{array}{*{20}{l}} {\;\;\;\;\;\;\;\;\;\;\;\;\;\;f\left( \sigma \right) = }\\ {\left\{ {\begin{array}{*{20}{l}} {0\;\;\;\;\;{\lambda _2} > 0{\rm{或}}{\lambda _3} > 0}\\ {\left( {1 - \exp \left( {\frac{{ - R_A^2}}{{2{\alpha ^2}}}} \right)} \right)\exp \left( {\frac{{ - R_B^2}}{{2{\beta ^2}}}} \right)\left( {1 - \exp \left( {\frac{{ - {S^2}}}{{2{c^2}}}} \right)} \right)}\\ {\;\;\;\;\;\;{\rm{其他}}} \end{array}} \right.} \end{array} $

(2)

式中，σ为进行汉森矩阵特征值计算时所使用的高斯滤波器的标准差，$ {R_A} = \left| {{\lambda _2}} \right|/\left| {{\lambda _3}} \right|,{R_B} = \left| {{\lambda _1}} \right|/\sqrt {{\lambda _2}{\lambda _3}} ,S = \sqrt {\lambda _1^2 + \lambda _2^2 + \lambda _3^2} $。

在这个函数中，R_B用来将管状结构与球状结构区分开来；R_A是用来区分管状结构和薄片结构。但是Frangi方法指出，这两个度量值虽然有比较好的几何特性，可以使得结果不受到图像灰度变化的影响。然而CT血管造影图像和MR图像的随机噪声容易使增强结果受到影响，因此需要利用“血管比周围组织具有更高的灰度值且血管只是少数高亮体素”^[3]的先验信息来弥补。因此S作为补充项被提出，它利用了Frobenius矩阵范数来表征特征值的特性，认为在背景区域或者随机噪声处，由于对比不强烈，矩阵的特征值都比较小，所以S也比较小，而在血管结构上，由于其高亮的特性，特征值相对较大。

但是在实际的实验中，S对应的阈值c需要依据图像的灰度分布做调整，这给设计自动分割算法带来了难度。并且式(2)的度量函数没有充分考虑到血管主支或者血管节点处往往膨胀粗大，实际上不是规则的管状结构，所以它在这些部位的增强效果不佳，如图 1所示。但血管粗大的主支和交叉却又是血管拓扑分析最为关心的。

为了解决上述问题，提出了一种改进的血管度量函数。充分利用血管比周围组织具有更高的灰度值且血管只是少数高亮体素的特点，使用灰度信息来代替矩阵范数设计噪声抑制项。该血管度量函数为

$ \begin{array}{l} \;\;\;\;\;\;\;\;f\left( \sigma \right) = \\ \left\{ \begin{array}{l} 0\;\;\;\;\;{\lambda _2} > 0或{\lambda _3} > 0\\ \left( {1 - \exp \left( {\frac{{ - R_A^2}}{{2{\alpha ^2}}}} \right)} \right)\exp \left( {\frac{{ - R_B^2}}{{2{\beta ^2}}}} \right)\left( {1 - \exp \left( {\frac{{ - {I^2}}}{{2{c^2}}}} \right)} \right)\\ \;\;\;\;\;\;其他 \end{array} \right. \end{array} $

(3)

式中，I=(y(x)-y_back)/(y_max-y_back)。y(x)是当前像素点的灰度值，y_max是感兴趣区域像素点最大的灰度值，y_back为肝脏与血管的分割灰度值。如图 2所示，统计感兴趣区域的所有像素点得到灰度直方图，取直方图峰值对应的灰度值作为y_back，只有当像素点的灰度值大于y_back，时，才认为该点属于血管。

区域增长是图像分割领域常用的分割算法。由于血管树的联通特性，所以区域增长算法适用于血管的分割。区域增长中决定是否将种子点邻近区域像素合并的标准多种多样, 因此区域增长的方法也有很多。Orkisz等人^[8]提出了结合变分模型来描述区域增长，以便于将先验信息加入到区域增长过程中。他们采取以下的方式来描述血管分割中的先验信息：对于当前需要被判定的像素点x而言，如果它属于血管区域，那么它的像素值v应该更接近已经划分出来的血管区域(Ω_in)的平均像素值μ_in，而远离背景区域(Ω_out)的平均像素值μ_out。假设以v(x)表示当前待判断像素点x的响应值，如果有x∈Ω_in，需要满足$ \left| {v\left( \boldsymbol{x} \right) - {\mu _{{\rm{in}}}}} \right| < \left| {v\left( \boldsymbol{x} \right) - {\mu _{{\rm{out}}}}} \right|$。

将血管先验信息作为判定准则，结合肝脏血管增强图像的直方图特性，提出了一种自适应的区域增长算法实现对增强图像的全自动分割。自适应区域增长包括以下两个关键步骤：

1)初始种子点选取。初始种子点应该选取属于血管概率最高的点。依据血管增强的度量函数，单一尺度的增强结果中响应值越大，其属于血管的概率越高。如图 3所示，统计增强响应值得到直方图表现出了非常强的单调性，靠近直方图右边的点属于血管的概率更大。

为了防止随机噪声干扰种子点选取，如图 4所示，以当前所选择的增强尺度作为半径设置区域，统计区域响应值的总和，采用遍历法选取总和值最大的区域作为种子点的选择范围，然后再在区域中选择响应值最大的像素点作为初始种子点。

2)阈值自适应的区域增长。依据判定准则：以v(x)表示当前待判断像素点x的响应值，如果有x∈Ω_in，则$ \left| {v\left( x \right) - {\mu _{{\rm{in}}}}} \right| < \left| {v\left( x \right) - {\mu _{{\rm{out}}}}} \right|$。在这里μ_in和μ_out不是固定不变的，而是随着区域增长过程而收敛的。μ_in的初始值由初始种子点的响应值决定，μ_out的初始值则为0。对于当前待判断像素点y而言，有

$ \left\{ \begin{array}{l} {\mu _{{\rm{in}}}} = \frac{{\sum {v\left( y \right)} }}{N}\;\;\;\;\;\boldsymbol{y} \in {\boldsymbol{\Omega} _{{\rm{in}}}}\\ {\mu _{{\rm{out}}}} = \frac{{\sum {v\left( y \right)} }}{M}\;\;\;\;\boldsymbol{y} \in {\boldsymbol{\Omega} _{{\rm{out}}}} \end{array} \right. $

(4)

式中，N表示当前已经过判断属于Ω_in区域的像素点的数目，M表示当前已经过判断属于Ω_out区域的像素点的数目。

由于肝脏血管树中不同级别的分支半径分布不均，单一尺度增强效果无法满足整个血管树的分割。常用的方法是选取不同的参数对血管进行多尺度增强并将不同尺度的增强结果融合，最后从融合结果中提取血管树。这其中需要克服的关键问题是，不同尺度的增强图像存在响应强度不一致性，融合参数的选取将会直接影响到最终的分割结果。

为此，提出一种全自动的血管树分割算法，在多尺度血管融合的基础上，采用迭代的区域增长来对增强结果进行分割，将上一次迭代的分割结果补充为下一次迭代的初始种子点，以减弱图像融合不一致对阈值μ_in和μ_out的影响。本文分割算法流程图如图 5所示，主要包括以下几个关键步骤：

1)多尺度的图像增强。在标准差σ下，图像的二阶导数可以表示为

$ {I_{xx}}\left( \boldsymbol{r} \right) = \frac{{{\partial ^2}}}{{\partial {x^2}}}G\left( \boldsymbol{r} \right) * I\left( \boldsymbol{r} \right) $

(5)

$ G\left( \boldsymbol{r} \right) = \frac{1}{{\sqrt {2\pi } \sigma }}\exp \left( { - \frac{{\left\| r \right\|}}{{2{\sigma ^2}}}} \right) $

(6)

式中，I(r)表示像素点的灰度。因此标准差σ的选取将影Hessian矩阵的特征值，从而决定增强结果的尺度。

Li等人^[9]提出在多尺度融合时，尺度σ的取值范围为[σ_min，σ_max]，对应的血管值范围为[d_min，d_max]，则两者满足关系$ \sigma = \frac{d}{4}$，因此假设增强总次数为N，第n次增强尺度σ_n可表示为

$ {\sigma _n} = {\left( {\frac{{{d_{\max }}}}{{{d_{\min }}}}} \right)^{\frac{{n - 1}}{N}}}\frac{{{d_{\min }}}}{4},n = 1,2,3, \cdots ,N $

(7)

2)标准化的图像融合。Frangi等人^[3]提出引入标准化参数γ来解决多尺度融合中响应不一致的问题。对不同尺度，在式(6)中加入标准化项，即

$ G\left( \boldsymbol{r} \right) = {\sigma ^\gamma }\frac{1}{{\sqrt {2\pi } \sigma }}\exp \left( { - \frac{{\left\| \boldsymbol{r} \right\|}}{{2{\sigma ^2}}}} \right) $

(8)

以此为基础，在本文所提出的框架中，不同尺度图像的融合采用取层间最大值的方法，且每一次增强以后均对之前所有尺度的增强结果进行融合。从对应最粗血管的σ_N开始，对于以σ_n为尺度的增强而言，像素点x的增强结果v(x)表示为

$ v\left( \boldsymbol{x} \right) = \mathop {\max }\limits_{{\sigma _N} \le \sigma \le {\sigma _n}} v\left( {\boldsymbol{x},\boldsymbol{\sigma} } \right) $

(9)

3)迭代式的自适应区域增长。经过标准化的图像融合虽然可以一定程度解决不同尺度增强相应结果不一致的问题，但是对于不同的数据，融合参数的选取会对最终的分割效果造成差异，这给设计自动算法设计带来了难度。因此，在固定融合参数的基础上，采用迭代式区域增长的方式来对增强图像进行分割。

迭代从最大尺度σ_N开始，第1次增强以后，直接用遍历法选取初始种子点进行自适应的区域增长。之后，在以σ_n为尺度进行增强之后，均先进行标准化图像融合并选取初始种子点，然后将上一层区域增长的结果，全部标记成种子点补充进入本层，以此为基础进行区域增长。迭代式的自适应区域增长的结果如图 6所示。从图 6中可以看出，这种方法可以有效地保证血管主支的完整性。

选择5例临床数据(Case1、Case2、Case3、Case4、Case5)来评估本文方法。实验数据来源于华中科技大学同济医院提供的16层螺旋肝脏CT造影图像。图像的层间距为1.25 mm, 层内间距为0.703 125 mm。采用的双线性的方法对处理前的数据进行插值，从而得到层间和层内间距均为0.703 125 mm的图像。实验平台采用的计算机配置为：Inter^® Core^TM i5-4660 CPU@3.20 GHz，8 GB内存，NVIDIA GeForce GTX 980显卡，64位操作系统。

按照图 5所示的分割流程，首先对输入图像进行各项异性滤波降噪和肝脏模板提取预处理。预处理的结果如图 7所示。

为了对比本文提出的新的度量函数与Frangi方法^[3]在血管主支及血管交叉处的增强效果，在实验中使用了多个不同的尺度对原始图像进行了血管增强，并将不同尺度的增强结果融合。最终得到的对比结果如图 8所示，图中红色区域为增强效果。从图 8中可以看到，新的度量函数在血管分叉处拥有更好的连续性，并且对粗大的主支有更好的增强效果。

在血管增强结果的基础上，比较了本文使用的迭代式自适应区域增长算法与传统的连接门限法区域增长算法的效果。设增强结果像素的最大值为I_max，最小值为I_min，I(x)代表候选像素点的像素值，选取I(x)∈[I_min+ω×(I_max-I_min), I_max]作为连接门限。在实验中，分别对ω取0.2、0.15和0.1得到3组门限，2维的分割结果如图 9所示。图中红色标记表示血管分割结果与原始图像的叠加。从图 9(a)(b)中可以看出，ω取0.2和0.15时，血管处于欠分割状态，血管主支出现了明显的断裂；从图 9(c)中可以看出，当ω取0.1时，血管处于过分割状态，噪声被当成血管分割进来，但此时血管主支的分割仍然不完全。图 9(d)展示了本文提出的迭代式自适应区域增长的分割结果，其中血管主支部分分割完全，细支没有出现明显的过分割。

最后，对5组临床数据进行了完整分割，前4组数据选取了门脉期CT图像作为输入，获得了门静脉的分割结果，最后一组数据选择了静脉期的图像作为输入，同时获得了门静脉和肝静脉的分割结果。图 10展示了最终的分割结果。从上至下每一排展示每个临床病例的分割结果，从左至右分别是该病例的血管体绘制可视化结果、血管分割在2维切片上的标记以及用血管体绘制结果穿过2维切片的展示。从血管体绘制结果可以看出，本文方法所分割的血管具有比较完整的树形结构；从2维叠加结果及3维切片展示中可以看出，本文方法可以标记出绝大部分的血管，在主支和细支的分割上都具有比较好的效果。

在分割基础上，进一步对分割的结果进行了拓扑分析。其中Case1至Case4为门脉期图像的门脉分析，Case5为静脉期图像中单独对静脉的分析，各病例可检测到的分支级别，以及不同级别分支数目的统计结果如表 1所示。从表 1中可以看出，所有的数据都可以检测到比较完整的5级拓扑结构，其中Case1至Case3可以检测到第6级的细小分支。

提出了一种自动的肝内血管分割方法，适用于多相期CT影像中的肝内门静脉及静脉提取。改进了传统的基于汉森矩阵的图像增强算法，提出将图像灰度信息和汉森矩阵特征值相结合来设计增强滤波器，解决了传统方法在血管节点处增强不连续的问题，然后提出了一种迭代式的自适应区域增长算法来对多尺度的增强结果进行进一步分割，这种方法采用“由粗到细”的策略，将上层的分割结果作为下层先验信息，以完成血管树的完整分割。使用5例临床CT病例评估了本文方法，结果显示本文方法对肝内血管提取具有良好的效果。未来，将使用更多的临床数据对本文方法进行评估，检测算法的鲁棒性。