本发明公开一种肝脏CTA序列图像的血管分割方法，首先对输入三维肝脏序列图像进行对比度增强和平滑噪声预处理，然后，采用OOF和OFA算法增强肝脏血管及其边界，细化血管中心；根据血管的几何结构，自动搜索血管中心线的种子点，并提取肝脏血管的中心线、构建肝脏血管树；最后，结合快速行进法初步分割肝脏血管并计算相应的血管和背景灰度直方图，采用图割算法实现肝脏血管的精确分割。本发明针对低对比度、强噪声和模糊边界的CTA序列图像，充分利用血管的几何形状和灰度信息，能有效地准确分割肝脏血管。本发明的肝脏CTA序列图像的血管分割方法能够推广到其他的三维血管分割中。
1.一种腹部CTA序列图像的肝脏血管分割方法，其特征在于包括以下步骤：步骤1：对3D肝脏图像进行预处理，提高血管的对比度，平滑CTA图像噪声首先，分别调节3D肝脏图像的窗宽/窗位到[400～800]/[200～500]范围，并将灰度值归一化到[0,255]，得到图像I0(x)，其中x为三维空间位置坐标，然后采用各向异性滤波算法去除噪声：其中，t为扩散时间，I(x,t)为t时刻的图象，div(·)为散度操作， 为图像梯度，c(x,t)为非负的单调递减传导系数：其中，k为扩散门限，k为[50～80]之间的自然数；
步骤2：采用多尺度OOF和OFA算法检测血管及其相应的边界，细化血管中心OOF滤波器通过最小化球面通量寻找最佳的投影方向，构建血管增强滤波器，其中*为卷积运算，Br表示以x为球心,r为半径的球表面，为球表面的外向单位法向量， 为球表面Br上的位置坐标，dA为球表面无穷小的面积分量， 是方差为σ0的高斯函数，I(x)为经各向异性滤波后得到的图象，为方向向量，Qr,x为3×3通量矩阵,其第i行、第j列元素为其中，i,j∈{1,2,3}，br(x)为以x为球心、r为半径的球体阶跃函数函数，即在该球体内br(x)取1，否则br(x)取0；
通过对通量矩阵Qr,x的特征值λ1(x,r)、λ2(x,r)、λ3(x,r)和特征向量e1(x,r)、e2(x,r)、e3(x,r)分析，令|λ1(x,r)|≥|λ2(x,r)|≥|λ3(x,r)|，设计OOF血管增强算法fOOF(x,r)＝-(λ1(x,r)+λ2(x,r)) (3)结合OFA算法
为提高计算效率，所述多尺度OOF和OFA均采用快速傅里叶变换及其反变换进行优化，并对血管中心进行细化，计算过程包括：(1)设置最小半径rmin、最大半径rmax和半径增量rstep，初始化r＝rmin；
(2)计算通量矩阵Qr,x的元素 i,j∈{1,2,3}，其中，FFT(·)和FFT-1(·)分别是快速傅里叶变换和反快速傅里叶变换，σ为高斯函数的方差，u为频域的位置矢量；
(3)计算通量矩阵Qr,x的特征值λ1(x,r)、λ2(x,r)、λ3(x,r)和特征向量e1(x,r)、e2(x,r)、e3(x,r)；
fOOF(x,r)＝-(λ1(x,r)+λ2(x,r))(5)计算半径为r的fOFA(x,r)
(6)结合OOF和OFA算法，细化血管中心M(x,r)：M(x,r)＝max(0,fOOF(x,r)-fOFA(x,r))(7)如果r＜rmax，r＝r+rstep，重复(2)和(6)，否则，按照如下公式计算多尺度fOOF(x)，fOFA(x)和M(x)响应：其中，R为半径r的集合；获取体素x为圆心的最优血管截面半径 时，对应的通量矩阵 的特性向量e3即为体素x的切线方向
步骤3：根据CTA图像血管特征，自动选取血管中心线种子点并提取相应的中心线，结合血管的几何结构，构建肝脏血管树血管中心线种子点位于血管截面中心，考虑到细化图像中，血管中心的灰度值最大，首先选取M(x)大于阈值tseed的体素作为候选种子点，然后，以每个候选种子点为中心，构建边长为rseed的立方体区域，在该区域内选择M(x)最大的体素为种子点xi；
血管中心线由一系列血管截面中心点连接而成，运用以下步骤提取，即：(1)沿种子点xi的正切线 方向搜索邻域体素；
(2)如果该邻域体素 的 值大于阈值tline,且 与切线方向 的夹角小于θline，将该邻域体素标记为中心线体素，并作为新的种子点；
(3)重复(2)，当无法检测到符合条件的邻域体素，或该邻域体素已被标记为中心线体素时循环结束；
(5)如果该邻域体素 的 值大于阈值tline，且与 与负切线方向 的夹角小于θline，将该邻域体素标记为中心线体素，并作为新的种子点；
(6)重复(5)，当无法检测到符合条件的邻域体素，或该邻域体素已被标记为中心线体素时循环结束；
采用(1)～(6)步骤，遍历所有种子点，实现血管中心线的提取；所述的rseed为4～10之间的自然数；tseed选择为按M(x)从大到小排列在前4～10％的灰度值，tline选择前10～30％的值，nvoxels为10～30之间的自然树，tconn选择为0.05～1之间的值；
基于肝脏血管的几何结构，构建静脉期肝脏血管树：首先，统计中心线的体素个数，去除体素个数小于nvoxels的孤立中心线，然后，以中心线端点xe为中心，在边长为15～20的立方体区域内搜索其他中心线体素，计算xe与其他中心线体素xj的连接概率其中，d为xe到xj的距离，为以xe为圆心的最优血管截面半径，θtree为xe的切线方向与的夹角；当同时满足以下条件时，将xe与最大conn(xe,xj)的体素xj连接起来：a)conn(xe,xj)大于tconn；
最后，根据血管的连通性，选择体素数目最多的2条中心线作为肝静脉和肝门静脉血管树；
步骤4：结合快速行进法初步分割肝脏血管并计算相应的血管和背景灰度直方图，构建图割能量函数，实现肝脏血管的分割首先，以肝脏血管树作为初始目标区域，采用快速行进法初步分割肝脏血管，分别统计目标和背景的灰度直方图；然后，利用血管的灰度分布和边界信息，构建图割能量函数，实现对肝脏血管的准确分割：其中，α和β分别为区域项和边界项的权重系数；Lp为体素P的标签，Lp的取值为1和0，Rp为区域项， hobj和hbkj分别为归一化的初步分割目标和背景灰度直方图；P表示体素p的集合，Bpq为边界项，dist(p,q)为体素p到邻域体素q的距离，N为体素p的邻域体素的集体；最后，采用优化算法最小化该能量函数，实现肝脏血管的分割。
一种肝脏CTA序列图像的血管分割方法
技术领域
[0001] 本发明属于医学图像处理领域，涉及CTA序列图像中的肝脏血管增强、中心线提取和肝脏血管分割。
背景技术
[0002] 肝脏血管分割和三维重建有助于准确获取腹部肝脏血管组织的整体信息，是计算机辅助肝脏疾病诊断和肝脏手术规划的前提。CTA(计算机断层扫描血管造影)是一种无创伤的成像技术，具有密度分辨率高、对人体损伤小等优点，广泛应用于肝脏疾病的评估和诊断。由于肝脏血管结构复杂、血管之间相互缠绕、不同个体差异较大，肝脏血管分割面临巨大挑战。在临床应用中，为构建肝脏血管模型，放射科专家需根据血管的先验知识，对每一张CTA切片进行手动分割。由于肝脏CTA序列图像的切片数量很大(通常每位病人的CTA数据包含200～500张切片)，手动分割肝脏血管劳动强度大、耗时长、易产生错误分割结果，而且，分割的精度高度依赖于放射科专家的经验，具有很大的主观性。
[0003] 现有的三维肝脏血管分割方法一般可以分为基于灰度和梯度两大类。受血流速度、数据采样时间和血管狭窄处的影响，造影剂在血液中分布不均匀，导致在CTA图像中血管区域的灰度不一致。譬如，血管中心的灰度高于血管边界的灰度，主血管的灰度高于小血管的灰度。单一的基于灰度或梯度的分割方法，如三维区域生长、模糊聚类等，都无法有效地提取低对比度的肝脏血管。近年来，主动轮廓模型及其混合模型被广泛应用到三维血管分割中，但是这类模型的演化曲面容易越过血管的弱边界而导致严重的过分割，且需交互提供血管的初始区域。此外，以上方法对小血管的分割难以奏效。
发明内容
[0004] 针对现有方法的不足，本发明提出了一种肝脏CTA序列图像的血管分割方法，并且能够推广到其他三维血管的分割中。
[0005] 本发明的肝脏CTA序列图像的血管分割方法，包括以下步骤：
[0006] 步骤1：对原始三维肝脏图像进行预处理，提高血管的对比度，平滑CTA图像噪声。
[0007] 所述步骤1包括：
[0008] 首先，分别调节3D肝脏图像的窗宽/窗位到[400～800]/[200～500]范围，并将灰度值归一化到[0,255]，得到图像I0(x)，其中x为三维空间位置坐标。然后采用各向异性滤波算法去除噪声：
[0009]
[0010] 其中，t为扩散时间，I(x,t)为t时刻的图像，div(·)为散度操作， 为图像梯度，c(x,t)为非负的单调递减传导系数：
[0011]
[0012] 其中，k为扩散门限，k为[50～80]之间的自然数。
[0013] 步骤2：采用多尺度OOF(Optimal Orient Flux)和OFA(Oriented Flux Anti-symmetry)算法增强血管及其相应的边界，细化血管中心。
[0014] 所述多尺度算法指针对肝脏主血管和小血管的半径(尺度)差异，选择不同半径的血管滤波器，提取滤波器的最大响应值得到最终结果。
[0015] 所述OOF滤波器指通过最小化球面通量
[0016]
[0017] 寻找最佳的投影方向，构建血管增强滤波器，其中，*为卷积运算，Br表示以x为球心，r为半径的球表面，为球表面的外向单位法向量， 为球表面 上的位置坐标，dA为球表面无穷小的面积分量， 是方差为σ0的高斯函数，I(x)为经各向异性滤波后得到的图像，为方向向量，Qr,x为3×3通量矩阵，其第i行、第j列元素为
[0018]
[0019] 其中，i,j∈{1,2,3}，br(x)为以x为球心、r为半径的球体阶跃函数，即在该球体内，br(x)取1，否则，br(x)取0；
[0020] 通过对通量矩阵Qr,x的特征值λ1(x,r)、λ2(x,r)、λ3(x,r)和特征向量e1(x,r)、e2(x,r)、e3(x,r)分析，令|λ1(x,r)|≥|λ2(x,r)|≥|λ3(x,r)|，设计OOF血管增强算法[0021] fOOF(x,r)＝-(λ1(x,r)+λ2(x,r))
[0022] 结合OFA算法
[0023]
[0024] 构建血管边界检测算法
[0025]
[0026] 为提高计算效率，所述多尺度OOF和OFA均采用快速傅里叶变换及其反变换进行优化，并对血管中心进行细化，计算过程包括：
[0027] (1)设置最小半径rmin、最大半径rmax和半径增量rstep，初始化r＝rmin；
[0028] (2)计算通量矩阵Qr,x的元素 i,j∈{1,2,3}，
[0029]
[0030] 其中，FFT(·)和FFT-1(·)分别是快速傅里叶变换和反快速傅里叶变换，σ为高斯函数的方差，u为频域的位置矢量；
[0031] (3)计算通量矩阵Qr，x的特征值λ1(x,r)、λ2(x,r)、λ3(x,r)和特征向量e1(x,r)、e2(x,r)、e3(x,r)；
[0032] (4)计算半径为r的fOOF(x,r)
[0033] fOOF(x,r)＝-(λ1(x,r)+λ2(x,r))
[0034] (5)计算半径为r的fOFA(x,r)
[0035]
[0036] 其中， (6)结合OOF和OFA算法，细化血管中心M(x,r)：
[0037] M(x,r)＝max(0,fOOF(x,r)-fOFA(x,r))
[0038] (7)如果r＜rmax，r＝r+rstep，重复(2)至(6)，否则，按照如下公式计算多尺度fOOF(x)，fOFA(x)和M(x)响应：
[0039]
[0040]
[0041] 其中，R为半径r的集合；获取以体素x为圆心的最优血管截面半径此时，对应的通量矩阵 的特征向量e3的方向即为体素x的切线方向
[0042] 步骤3：根据CTA图像血管特征，自动选取血管中心线种子点并提取相应的中心线，结合血管的几何结构，构建肝脏血管树。
[0043] 所述血管中心线种子点位于血管截面中心。考虑到细化图像中，血管中心的灰度值最大，首先选取M(x)大于tseed的体素作为候选种子点，然后，以每个候选种子点为中心，构建边长为rseed的立方体区域，在该区域内选择M(x)最大的体素作为种子点xi。
[0044] 所述血管中心线由一系列血管截面中心点连接而成，运用以下追踪算法提取，即：
[0045] (1)沿种子点xi的正切线 方向搜索邻域体素。
[0046] (2)如果该邻域体素 的 值大于阈值tline(tline＜tseed)，且 与切线方向的夹角小于θline，将该邻域体素标记为中心线体素，并作为新的种子点。
[0047] (3)重复(2)，当无法检测到符合条件的邻域体素，或该邻域体素已被标记为中心线体素时循环结束。
[0048] (4)沿种子点xi的负切线方向 搜索邻域体素。
[0049] (5)如果该邻域体素 的 值大于阈值tline(tline＜tseed)，且与 与负切线方向 的夹角小于θline，将该邻域体素标记为中心线体素，并作为新的种子点。
[0050] (6)重复(5)，当无法检测到符合条件的邻域体素，或该邻域体素已被标记为中心线体素时循环结束。
[0051] 采用上述方法，遍历所有种子点，实现血管中心线的提取。
[0052] 基于肝脏血管的几何结构，构建静脉期肝脏血管树。首先，统计中心线的体素个数，去除体素个数小于nvoxels的孤立中心线，然后，以中心线端点xe为中心，在边长为15～20的立方体区域内搜索其他中心线体素，计算xe与其他中心线体素xj的连接概率[0053]
[0054] 其中，d为xe到xj的距离，为以xe为圆心的最优血管截面半径，θtree为xe的切线方向与 的夹角。当同时满足以下条件时，将xe与最大conn(xe,xj)的体素xj连接起来：
[0055] (a)conn(xe,xj)大于tconn；
[0056] (b)θtree小于
[0057] (c)d小于20。
[0058] 最后，根据血管的连通性，选择体素数目最多的2条中心线作为肝静脉和肝门静脉血管树。
[0059] 所述rseed优选4～10之间的自然数；tseed选择为按M(x)从大到小排列在前4～10％的灰度值，tline选择按M(x)从大到小排列在前10～30％的值，θline选择为 之间的值，nvoxels优选10～30之间的自然数，tconn选择为0.05～1之间的值。
[0060] 步骤4：结合快速行进法初步分割肝脏血管并计算相应的血管和背景灰度直方图，构建图割能量函数，实现肝脏血管的准确分割。
[0061] 首先，以肝脏血管树作为初始目标区域，采用快速行进法初步分割肝脏血管，分别统计其目标和背景的灰度直方图。然后，利用血管的灰度分布和边界信息，构建图割能量函数，实现肝脏血管的准确分割：
[0062]
[0063] 其中，α和β分别为区域项和边界项的权重系数；Lp为体素p的标签，Lp的取值为1和0，Rp为区域项， hobj和hbkj分别为归一化的初步分割目标和背景灰度直方图，P表示体素p的集合；Bpq为边界项，
dist(p,q)为体素p到邻域体素q的距离；N为体素p的邻域体素的集合，最后，采用优化算法最小化该能量函数，实现肝脏血管的分割。
[0064] 与现有方法相比，本发明具有如下优点：
[0065] (1)提供了一种自动的肝脏血管分割方法，针对每位病人腹部CTA数据采用OOF和OFA分别增强血管及其边界，优化血管的中心响应，利用血管的几何结构和灰度分布，提取血管的树形拓扑结构，采用图割算法实现准确分割，可有效处理复杂的肝脏血管，鲁棒性强；
[0066] (2)不需要手工选取中心线种子点，能自动提取血管中心线并构建血管树；
[0067] (3)不需要手工选取血管的初始分割区域；
[0068] (4)本发明的肝脏CTA序列图像的血管分割方法能推广到其他三维血管的分割中，如脑部MRI图像的血管分割。
附图说明
[0069] 图1为本发明的流程图；
[0070] 图2为本发明对4位病人肝脏CTA序列图像的处理结果图；
[0071] 图3为不同方法误分割体素到最近的手动分割血管表面的距离误差统计结果图。
具体实施方式
[0072] 图1是本发明实施的肝脏CTA序列图像的血管分割方法流程图。首先从输入的肝脏血管图像调节窗宽/窗位，提高血管的对比度，并采用各向异性滤波平滑噪声。然后，采用OOF和OFA方法增强血管及其边界，优化血管的中心响应。接下来，根据血管的几何结构，提取血管中心线并构建血管树。最后，运用快速行进法初步分割肝脏血管，结合其初始血管和背景的灰度分布，构建图割能量函数，优化该能量函数实现肝脏血管的准确分割。
[0073] 下面结合图1，以一实施例详细说明本发明的肝脏CTA序列图像的血管分割方法。
[0074] 1、预处理。首先输入肝脏CTA序列图像，其切片平面像素间距为0.5～1.0mm，层间距为0.7～2.0mm，调节该序列图像的窗宽/窗位到300/500，并将灰度值归一化到[0～255]，得到图像I0(x)，其中x为三维空间位置坐标，然后采用各向异性滤波方法去除噪声：
[0075]
[0076] 其中，t为扩散时间，I(x,t)为t时刻的图像，div(·)为散度操作， 为图像梯度，c(x,t)为非负的单调递减传导系数：
[0077]
[0078] 其中，k为扩散门限，k选取为70。
[0079] 2.采用多尺度OOF(Optimal Orient Flux)和OFA(Oriented Flux Anti-symmetry)算法增强血管及其相应的边界，细化血管中心。
[0080] 为提高计算效率，所述多尺度OOF和OFA均采用快速傅里叶变换及其反变换进行优化，并对血管中心进行细化，计算过程包括
[0081] (1)设置最小半径rmin＝0.5，最大半径rmax＝8.5和半径增量rstep＝1，初始化r＝rmin。
[0082] (2)计算通量矩阵Qr,x的元素 i,j∈{1,2,3}，
[0083]
[0084] 其中，FFT(·)和FFT-1(·)分别是快速傅里叶变换和反快速傅里叶变换，σ为高斯函数的方差，σ设置为0.5，u为频域的位置矢量；
[0085] (3)计算通量矩阵Qr,x的特征值λ1(x,r)、λ2(x,r)、λ3(x,r)和特征向量e1(x,r)、e2(x,r)、e3(x,r)。
[0086] (4)计算半径为r的fOOF(x,r)
[0087] fOOF(x,r)＝-(λ1(x,r)+λ2(x,r))
[0088] (5)计算半径为r的fOFA(x,r)
[0089]
[0090] 其中，
[0091] (6)结合OOF和OFA算法，细化血管中心M(x,r)：
[0092] M(x,r)＝max(0,fOOF(x,r)-fOFA(x,r))
[0093] (7)如果r＜rmax，r＝r+rstep，重复(2)至(6)；否则，按照如下公式计算最终多尺度fOOF(x)，fOFA(x)和M(x)响应：
[0094]
[0095]
[0096]
[0097] 其中，R为半径r的集合；获取以体素x为圆心的最优血管截面半径此时，对应的通量矩阵 的特征向量e3的方向即为体素x的切线方向
[0098] 3、根据CTA序列图像血管特征，自动选取血管中心线种子点并提取相应的中心线，再结合血管的几何结构，构建肝脏血管树。
[0099] 首先选取M(x)大于40的体素作为候选种子点，然后，以每个候选种子点为中心，构建边长为10的立方体区域，在该区域内选择M(x)最大的体素作为种子点xi，最后，运用以下算法提取血管中心线，即：
[0100] (1)沿种子点xi的正切线 方向搜索邻域体素。
[0101] (2)如果该邻域体素 的 值大于阈值20，且 与切线方向 的夹角小于π/3，将该邻域体素标记为中心线体素，并作为新的种子点。
[0102] (3)重复(2)，当无法检测到符合条件的邻域体素，或该邻域体素已被标记为中心线体素时循环结束。
[0103] (4)沿种子点xi的负切线方向 搜索邻域体素。
[0104] (5)如果该邻域体素 的 值大于阈值20，且与 与负切线方向 的夹角小于π/3，将该邻域体素标记为中心线体素，并作为新的种子点。
[0105] (6)重复(5)，当无法检测到符合条件的邻域体素，或该邻域体素已被标记为中心线体素时循环结束。
[0106] 采用上述方法，遍历所有种子点，实现血管中心线的提取。
[0107] 基于血管几何结构，构建静脉期肝脏血管树。首先，统计中心线的体素个数，去除中心线体素个数小于10的孤立中心线。然后，以中心线端点xe为中心，在边长为20的立方体区域内搜索其他中心线体素，计算xe与其他中心线体素xj的连接概率：
[0108]
[0109] 其中，d为xe到xj的距离，为以xe为圆心的最优血管截面半径，θtree为xe的切线方向与 的夹角。当同时满足以下条件时，将xe与最大conn(xe,xj)的体素xj连接起来：
[0110] (a)conn(xe,xj)大于0.1；
[0111] (b)θtree小于
[0112] (c)d小于20。
[0113] 最后，基于血管的连通性，选择体素数目最多的2条中心线作为肝静脉和肝门静脉血管树。
[0114] 4、结合快速行进法初步分割肝脏血管并计算相应的血管和背景灰度直方图，构建图割能量函数，实现肝脏血管的准确分割。
[0115] 首先，以肝脏血管树作为初始目标区域，采用快速行进法初步分割肝脏血管，分别统计其目标和背景的灰度直方图。然后，基于血管的灰度分布和边界信息，构建图割能量函数，实现肝脏血管的准确分割：
[0116]
[0117] 其中，α和β分别为区域项和边界项的权重系数，α＝1.2，β＝1；Lp为体素p的标签，Lp的取值为1和0，Rp为区域项， hobj和hbkj分别为归一化的初步分割目标和背景灰度直方图，P表示体素p的集合；Bpq为边界项，dist(p,q)为体素p到邻域体素q的距离；N为体素p的邻域体素的集合，最后，采用优化算法最小化该能量函数，实现肝脏血管的准确分割。
[0118] 为验证本发明的有效性，对6位病人的肝脏CTA序列图像进行了实验，其切片平面像素间距为0.5～1mm，层间距为0.7～2mm，血管的金标准由放射科专家手动分割产生。图2给出了随机抽取的4位病人肝脏CTA序列图像的血管处理结果图，其中第一行为调节窗宽/窗位的三维血管图像，第二行为基于OOF增强的血管图像，第三行为基于OFA增强的血管边界图像，第四行为细化的血管中心图像，第五行为提取的肝静脉和门静脉血管树，最后一行为采用图割算法分割的三维肝脏血管图像。从图中可以看出，本发明能有效地分割CTA序列图像肝脏血管，分割结果中包含较丰富的小血管，血管表面光滑，且连通性较好。
[0119] 为定量评价本发明的有效性，采用精确度(Accuracy)、敏感度(Sensitivity)和特性度(Specificity)对肝脏血管分割结果进行评价，它们的计算公式分别为：
[0120]
[0121]
[0122]
[0123] 其中，TP表示正确分割的血管体素数目，TN表示正确分割的背景体素数目，FP表示本来属于背景但被误分割为血管的体素数目，而FN则表示本来属于血管但被误分割为背景的体素数目。表1给出了本发明对6位病人的肝脏CTA序列图像的血管分割结果的精确度、敏感度和特异度，可以看出，本发明的精确度、敏感度和特异度都比较高，分割性能优异。
[0124] 表1本发明对6位病人肝脏CTA序列图像的血管分割结果的性能指标
[0125]
[0126] 为进一步评价本发明方法，对所有错误分割体素到最近的手动分割血管表面的距离进行统计，并与主动轮廓算法(GAC)和CV模型进行比较，统计结果如图3所示，其中，横坐标为被错误分割的体素到最近的手动分割血管表面的误差距离，纵坐标为相应的误分割体素数目，可以看出，在相同误差距离下，本发明方法的误分割体素数目小于主动轮廓和CV模型，从而进一步证明了本发明方法的优越性。
