CN108492300A - 管状结构增强与能量函数结合的肺部血管树分割方法 - Google Patents
管状结构增强与能量函数结合的肺部血管树分割方法 Download PDFInfo
- Publication number
- CN108492300A CN108492300A CN201810222634.XA CN201810222634A CN108492300A CN 108492300 A CN108492300 A CN 108492300A CN 201810222634 A CN201810222634 A CN 201810222634A CN 108492300 A CN108492300 A CN 108492300A
- Authority
- CN
- China
- Prior art keywords
- tubulose
- response
- tree
- formula
- image
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 64
- 230000002792 vascular Effects 0.000 title claims abstract description 50
- 230000002685 pulmonary effect Effects 0.000 title claims abstract description 45
- 230000004044 response Effects 0.000 claims abstract description 58
- 230000011218 segmentation Effects 0.000 claims abstract description 46
- 230000006870 function Effects 0.000 claims abstract description 34
- 210000004072 lung Anatomy 0.000 claims abstract description 19
- 238000009792 diffusion process Methods 0.000 claims abstract description 17
- 230000002708 enhancing effect Effects 0.000 claims abstract description 17
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 14
- 238000004364 calculation method Methods 0.000 claims abstract description 8
- 230000014509 gene expression Effects 0.000 claims description 19
- 239000011159 matrix material Substances 0.000 claims description 19
- 239000013598 vector Substances 0.000 claims description 12
- 238000001914 filtration Methods 0.000 claims description 10
- 238000009499 grossing Methods 0.000 claims description 6
- 230000035945 sensitivity Effects 0.000 claims description 5
- 235000004338 Syringa vulgaris Nutrition 0.000 claims description 4
- 244000297179 Syringa vulgaris Species 0.000 claims description 4
- 238000009795 derivation Methods 0.000 claims description 4
- 238000002474 experimental method Methods 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 3
- 230000010339 dilation Effects 0.000 claims description 3
- 230000009977 dual effect Effects 0.000 claims description 3
- 230000000877 morphologic effect Effects 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 239000008280 blood Substances 0.000 claims description 2
- 210000004369 blood Anatomy 0.000 claims description 2
- 238000003325 tomography Methods 0.000 claims description 2
- 210000004204 blood vessel Anatomy 0.000 abstract description 13
- 238000004458 analytical method Methods 0.000 description 3
- 238000004445 quantitative analysis Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 208000029523 Interstitial Lung disease Diseases 0.000 description 1
- 206010028980 Neoplasm Diseases 0.000 description 1
- 208000010378 Pulmonary Embolism Diseases 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 208000020816 lung neoplasm Diseases 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 230000010412 perfusion Effects 0.000 description 1
- 210000001147 pulmonary artery Anatomy 0.000 description 1
- 208000002815 pulmonary hypertension Diseases 0.000 description 1
- 210000003492 pulmonary vein Anatomy 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004043 responsiveness Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 210000003462 vein Anatomy 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/11—Region-based segmentation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/155—Segmentation; Edge detection involving morphological operators
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30061—Lung
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30101—Blood vessel; Artery; Vein; Vascular
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明涉及一种管状结构增强与能量函数结合的肺部血管树分割方法，利用Pock函数计算管状结构响应度，从而检测出潜在的血管区域。然后采用基于扩散张量的管状结构增强算法对原始图像进行增强，降低噪声对原始图像的影响并增强血管区域。最后将Pock函数计算结果和图像增强结果相结合构建区域描述算子，并利用最小化能量分割方法VRG法对肺部血管进行精细分割。其分割结果显示该方法在分割出肺部主支血管的同时，提取出了大量的细小血管，且分割结果受噪声的影响较小。本方法特异性高，敏感性较强，同时能够区分血管与气管壁区域，进一步提高了分割结果的准确性。
Description
技术领域
本发明涉及一种医学图像处理技术，特别涉及一种管状结构增强与能量函数结合的肺部血管树分割方法。
背景技术
肺部血管由肺动脉与肺静脉组成，是人体各组织器官中最为复杂的血管结构之一。从肺主动脉和肺主静脉开始，肺部血管逐级分支形成类树状的血管树结构。在临床诊断中，准确获取肺部血管树的解剖结构信息，是评估肺动脉高压风险的重要参考依据，也是实现肺栓塞自动检测的基础，同时也有利于减少肺结节检测的假阳性率。在临床研究中，有效地分离出肺部血管树，对于肺灌注研究、间质性肺疾病研究以及肿瘤体积量化分析具有重要的临床意义。在图像处理领域，提取出的肺部血管树，还可用于引导肺气管和肺叶组织的分割。因此有效准确地分割出肺部血管区域，具有重要的临床意义与研究价值。而由于肺部血管树分布范围广、血管半径跨度大且分支数目庞大，完整有效地分割出肺部血管树难度依旧较大。针对这一难题，国内外提出的肺部血管分割方法相对较少。其中能够完整分割出肺部血管区域的方法又占少数，并且用于评判分割结果的量化指标不够全面。主要的分割方法包括基于区域生长的方法、基于水平集的方法以及基于管状滤波增强的方法。基于区域生长的方法能够有效地分割出管半径较大的肺部血管，但对细小血管的分割效果较差，且容易错分出肺部气管壁区域。基于水平集的分割方法，具有较高的分割精度，但由于引入了水平集函数，使得完整分割肺部血管的计算量较大。基于管状滤波函数的分割方法种类较多，该类方法主要通过对Hessian矩阵进行分析，从而提取出潜在的肺部血管区域。但该类方法容易在血管分叉处产生断裂，从而影响最终的分割结果。
发明内容
本发明是针对肺部血管树分割中存在的细小血管分割不完全、错分出气管壁区域以及计算量大的问题，提出了一种管状结构增强与能量函数结合的肺部血管树分割方法，以实现对肺部血管树的有效分割。
本发明的技术方案为：一种管状结构增强与能量函数结合的肺部血管树分割方法，具体包括如下步骤：
1)输入待分割DICOM格式的胸部CT序列断层图像；
2)利用文献《Automatic Lung Segmentation for Accurate Quantitation ofVolumetric X-Ray CT Images》中的阈值法与形态学修补法，对步骤1)的胸部CT序列断层图像肺部区域进行分割并获取掩膜Mask1；
3)利用文献《Two-pass region growing combined morphology algorithm forsegmenting airway tree from CT chest scans》中的双路径区域生长结合形态学重建的方法对步骤1)的胸部CT序列断层图像中气管树进行分割，获取不含气管壁的气管树；
4)对步骤3中分割得到的肺部气管树，利用多尺度球形结构元素进行形态学开操作，从而剥离出含有不同等级分枝的气管树；然后对剥离出的各个气管树，使用对应尺度的球形结构元素进行膨胀操作，从而获取含有气管壁区域的气管树Airway；
5)从步骤2)所得肺部区域掩膜Mask1中去除气管树Airway对应的区域，从而获取肺部血管分割掩膜Mask；
6)在分割掩膜Mask区域内，利用Pock管状响应函数，对步骤1)中的图像进行多尺度管状响应计算，并获取管状响应结果；
6-1)设定尺度集合σi，σmin≤σi≤σmax，其中最小尺度σmin等于0.5，最大尺度σmax等于8，尺度间隔为0.5，同时设定Pock响应函数对称敏感系数σw为0.2，噪声抑制项δ设置为0.2；
6-2)对步骤1)中的图像进行高斯核标准差大小为σi的高斯滤波，i＝0.5，1，1.5，......，8，然后求取滤波结果体素点的Hessian矩阵，用求取的Hessian矩阵乘以然后计算Hessian矩阵的特征值与特征向量；
6-3)利用Pock管状响应函数公式(1)至(8)，计算在尺度σ下的管状响应；
其中R+(x，σ，θ)为对称性约束的管状结构响应；x代表图像中的任意体素点；σ代表当前尺度；θ代表待检测管状结构半径与尺度σ间的比例系数；R-(x，σ)代表体素点x在尺度σ下的梯度幅值；δ代表噪声抑制项；
N为离散圆周点个数，w(bi)为对称系数，bi为原始边界响应，
I表示步骤1)中胸部CT序列断层图像；bi表示在尺度σ下，第i个圆周点处的梯度幅值，即边界响应，它表示了该点属于边界点的可能性大小；Vαi表示第i个圆周点对应的旋转相量，用于估算圆周点的位置；
V1和V2分别表示Hessian矩阵主曲率方向对应的两个特征向量，对应的特征值关系为|λ1|＞|λ2|；αi代表第i个圆周点的离散相量角；
其中为圆周点的平均边界响应；σω为检测对称性的敏感系数，取值范围为(0，1]；
6-4)利用管状响应规范化公式，对尺度σi下的管状响应进行响应规范化，规范化公式如下所示：
Rnorm(x，σi，θ)＝σi γR(x，σi，θ)
其中Rnorm(x，σi，θ)代表规范化后的管状响应，R(x，σi，θ)表示未经规范化的管状响应，γ代表规范化系数，取值为1；
6-5)重复子步骤6-2)至6-4)以求取各尺度下的规范化管状响应，并利用多尺度管状响应公式求取最终的管状响应，多尺度管状响应公式如下所示：
7)使用基于三维扩散张量的图像增强方法，对Mask掩膜内的区域进行增强；
8)利用步骤6)与步骤7)的计算结果，依据公式(21)构建区域描述算子k(x)；
P(x)表示体素点x的管状响应；apin表示被判定为血管区域体素点的平均管状响应；对应的apout表示非血管区域体素点的平均管状响应；MP表示掩膜Mask内的最大管状响应；以此类推，T(x)表示体素点x增强后的响应，与分别表示被划分为血管区域与非血管区域体素点的平均增强响应值；MT表示掩膜Mask内的最大增强响应值；
9)利用步骤8)中构建的区域描述算子以及初始化的分割区域，在掩膜Mask内利用VRG算法对肺部血管进行精细分割：
9-1)求取Pock函数管状响应结果的局部极大值区域，利用局部极大值区域作为初始分割种子集Seeds，并设置种子点集对应体素点的初始状态φ(x)为1；
9-2)从种子集Seeds开始，对内外边界点依据公式(22)和(23)进行能量变化计算，从而纳入或排除血管区域内的体素点，分割进程将持续至能量状态稳定即∑φn+1(x)＝∑φn(x)；
VRG算法纳入或排除体素点的能量方程为：
ΔJ(φn+1(x))＝(1-2φn(x))k(x) (22)
ΔJ代表能量变化项，φ(x)代表体素点x的状态值，它的表达式如公式(23)所示，n+1与n表示迭代次数：
Ωin表示血管区域，Ωout表示非血管区域，依据能量最小化原则，仅当体素点的能量变化值为负值时，该体素点将被纳入或排除出血管区域。
所述步骤4)实现的具体步骤如下：
4-1)设置开操作结构元素半径其中半径为毫米单位，共取七个尺度；
4-2)利用具有不同半径的结构元素对步骤3)中获取的气管树进行形态学开操作，从而剥离出含有不同等级分枝的气管树集合
4-3)设置膨胀操作结构元素的半径为开操作结构元素的二分之一，即
4-4)对步骤4-2)中获取的气管树集合利用对应的膨胀操作结构元素进行形态学膨胀操作，从而获取气管树集合将所得气管树集合进行并集，从而获取包含气管壁区域的气管树Airway。
所述步骤7)实现的具体步骤如下：
7-1)设置实验参数，设置ρ＝1以及σ＝1分别用于归一化图像的高斯平滑滤波与平滑三维结构张量的高斯平滑滤波；同时设置扩散模型相关参数，C＝3.31488，λc＝0.02，λe＝0.02，λh＝0.5。迭代的次数为5，步长为0.5；
7-2)利用高斯核标准差为1的高斯滤波器对归一化图像进行滤波，并对滤波结果进行求导：
对归一化的图像进行高斯滤波，然后对滤波后的图像分别使用模板与在三个方向上进行滤波，并获取滤波结果ux，uy和uz。其中三个模板都是3*3*3的三维矩阵，对应每一层的二维矩阵表达式分别如公式(9)，(10)以及(11)所示；
公式(9)中，表示x方向模板第一层矩阵表达式，为第二层矩阵表达式，为第三层表达式，以此类推；
7-3)依据公式(12)构建三维结构张量，对三维结构张量进行高斯核标准差为1的高斯滤波，然后计算其对应的特征值与特征向量；
7-4)依据公式(13)至(18)，构建基于连续转换的混合扩散HDCS模型的三维扩散张量：
表示一致增强扩散CED模型；表示边缘增强扩散EED模型；ε为权重值，用于控制EED模型与CED模型的连续转化；
α为经验值取值为0.001；λc为管状结构对比参数取值为0.02；κ＝(μ2/(α+μ3))4；代表图像梯度；则代表σ尺度下的图像梯度；μ2与μ3代表三维结构张量的两个特征值；μ1，μ2与μ3分别代表三维结构张量的三个特征值，且三个特征值的绝对值大小关系为：μ3≥μ2≥μ1；
C与m分别为阈值参数与经验取值。C取值为3.31488，m取值为4。λe为平面结构对比参数，并且取经验值为0.02，
λn为利用HDCS模型构建的特征值，υn1、υn2、υn3为三维结构张量所对应的三个特征向量；
7-5)依据公式(19)至(20)，进行迭代计算并获取最终的增强结果：
代表图像梯度；u表示归一化后的原始图像；div表示散度算子，*符号代表卷积，·符号代表点乘；同时M表示空间内核，内核值p及排列方式由Kroon等人在文献《Optimized anisotropic rotational invariant diffusion scheme on cone-beam CT》中提出的数值优化方法计算得出；迭代计算方程由公式(20)所示，其中k表示迭代次数，τ表示迭代时间步长；
本发明的有益效果在于：本发明管状结构增强与能量函数结合的肺部血管树分割方法，其分割结果显示该方法在分割出肺部主支血管的同时，提取出了大量的细小血管，且分割结果受噪声的影响较小。本方法特异性高，敏感性较强，同时能够区分血管与气管壁区域，进一步提高了分割结果的准确性。
附图说明
图1为本发明管状结构增强与能量函数结合的肺部血管树分割方法流程图；
图2为针对VESSEL12竞赛案例采用本发明方法分割所得的结果图。
具体实施方式
管状结构增强与能量函数结合的肺部血管树分割方法，利用Pock函数计算管状结构响应度，从而检测出潜在的血管区域。然后采用基于扩散张量的管状结构增强算法对原始图像进行增强，降低噪声对原始图像的影响并增强血管区域。最后将Pock函数计算结果和图像增强结果相结合构建区域描述算子，并利用最小化能量分割方法VRG法对肺部血管进行精细分割。
如图1所示管状结构增强与能量函数结合的肺部血管树分割方法流程图，包括如下步骤：
步骤1，输入待分割DICOM格式的胸部CT序列断层图像(原始图像)。
步骤2，利用文献《Automatic Lung Segmentation for Accurate Quantitationof Volumetric X-Ray CT Images》中的阈值法与形态学修补法，对步骤1的胸部CT序列断层图像肺部区域进行分割并获取掩膜Mask1。
步骤3，利用文献《Two-pass region growing combined morphology algorithmfor segmenting airway tree from CT chest scans》中的双路径区域生长结合形态学重建的方法对步骤1的胸部CT序列断层图像中气管树进行分割，获取不含气管壁的气管树。
步骤4，对步骤3中分割得到的肺部气管树，利用多尺度球形结构元素进行形态学开操作，从而剥离出含有不同等级分枝的气管树；然后对剥离出的各个气管树，使用对应尺度的球形结构元素进行膨胀操作，从而获取含有气管壁区域的气管树Airway。
4-1：设置开操作结构元素半径其中半径为毫米单位，共取七个尺度。
4-2：利用具有不同半径的结构元素对步骤3中获取的气管树进行形态学开操作，从而剥离出含有不同等级分枝的气管树集合
4-3：设置膨胀操作结构元素的半径为开操作结构元素的二分之一，即
4-4：对子步骤4-2中获取的气管树集合利用对应的膨胀操作结构元素进行形态学膨胀操作，从而获取气管树集合将所得气管树集合进行并集，从而获取包含气管壁区域的气管树Airway。
步骤5，从步骤2所得肺部区域掩膜Mask1中去除气管树Airway对应的区域，从而获取肺部血管分割掩膜Mask。
步骤6，在分割掩膜Mask区域内，利用Pock管状响应函数，对步骤1中的图像进行多尺度管状响应计算，并获取管状响应结果。
在掩膜Mask与步骤1中图像的乘积中计算管状响应。掩膜Mask是一个由数值0和1组成的矩阵，1代表肺实质区域。0代表背景区域，以及步骤3，4中提取的气管与气管壁区域。因此用掩膜Mask与原始图像相乘，就能把后续计算限定在Mask中值为1的区域(Mask其他区域为0，与原始图像相乘后还是0)。从而减少计算量，同时排除气管壁对血管分割精度的影响。因为在CT图像上，血管和气管壁紧密相连，亮度也几乎一样。
6-1：设定尺度集合σi，σmin≤σi≤σmax，其中最小尺度σmin等于0.5，最大尺度σmax等于8，尺度间隔为0.5。同时设定Pock响应函数对称敏感系数σw为0.2，噪声抑制项δ设置为0.2。
6-2：对步骤1中的图像进行高斯核标准差大小为σi(i＝0.5，1，1.5，......，8)的高斯滤波，然后求取滤波结果体素点的Hessian矩阵。用求取的Hessian矩阵乘以然后计算Hessian矩阵的特征值与特征向量。
6-3：利用Pock管状响应函数公式(1)至(8)，计算在尺度σ下的管状响应；
其中R+(x，σ，θ)为对称性约束的管状结构响应。x代表图像中的任意体素点，σ代表当前尺度，θ代表待检测管状结构半径与尺度σ间的比例系数。R-(x，σ)代表体素点x在尺度σ下的梯度幅值；δ代表噪声抑制项。
N为离散圆周点个数，w(bi)为对称系数，bi为原始边界响应，
I表示原始图像。bi表示在尺度σ下，第i个圆周点处的梯度幅值，即边界响应。它表示了该点属于边界点的可能性大小。表示第i个圆周点对应的旋转相量，用于估算圆周点的位置。
V1和V2分别表示Hessian矩阵主曲率方向对应的两个特征向量，对应的特征值关系为|λ1|＞|λ2|。αi代表第i个圆周点的离散相量角。
其中为圆周点的平均边界响应，表达式如公式(8)所示。σω为检测对称性的敏感系数，取值范围为(0，1]。
对称系数w(bi)权衡了各圆周点的边界响应对最终管状响应的贡献程度，从而使得具有高对称性的管状结构响应强，低对称性的管状结构响应弱。
6-4：利用管状响应规范化公式，对尺度σi下的管状响应进行响应规范化，规范化公式如下所示：
Rnorm(x，σi，θ)＝σi γR(x，σi，θ)
其中Rnorm(x，σi，θ)代表规范化后的管状响应，R(x，σi，θ)表示未经规范化的管状响应，γ代表规范化系数，取值为1。
6-5：重复子步骤6-2至6-4以求取各尺度下的规范化管状响应，并利用多尺度管状响应公式求取最终的管状响应。多尺度管状响应公式如下所示：
步骤7，使用基于三维扩散张量的图像增强方法，对Mask掩膜内的区域进行增强：
7-1：设置实验参数，设置ρ＝1以及σ＝1分别用于归一化图像的高斯平滑滤波与平滑三维结构张量的高斯平滑滤波。同时设置扩散模型相关参数，C＝3.31488，λc＝0.02，λe＝0.02，λh＝0.5。迭代的次数为5，步长为0.5。
7-2：利用高斯核标准差为1的高斯滤波器对归一化图像进行滤波，并对滤波结果进行求导，求导内核如公式(9)，(10)与(11)所示。
对归一化的图像进行高斯滤波，然后对滤波后的图像分别使用模板与在三个方向上进行滤波，并获取滤波结果ux，uy和uz。其中三个模板都是3*3*3的三维矩阵，对应每一层的二维矩阵表达式分别如公式(9)，(10)以及(11)所示。
公式(9)中，表示x方向模板第一层矩阵表达式，为第二层矩阵表达式，为第三层表达式，以此类推。
7-3：依据公式(12)构建三维结构张量，对三维结构张量进行高斯核标准差为1的高斯滤波，然后计算其对应的特征值与特征向量。
7-4，依据公式(13)至(18)，构建基于连续转换的混合扩散模型(HDCS)的三维扩散张量。
表示一致增强扩散模型(CED)，表达式如公式(14)所示。表示边缘增强扩散模型(EED)，表达式如公式(15)所示，ε为权重值，用于控制EED模型与CED模型的连续转化。
α为经验值取值为0.001。λc为管状结构对比参数取值为0.02，κ＝(μ2/(α+μ3))4。代表图像梯度(则代表σ尺度下的图像梯度)，μ2与μ3代表三维结构张量的两个特征值。μ1，μ2与μ3分别代表三维结构张量的三个特征值，且三个特征值的绝对值大小关系为：μ3≥μ2≥μ1。
C与m分别为阈值参数与经验取值。C取值为3.31488，m取值为4。λe为平面结构对比参数，并且取经验值为0.02，
λn为利用HDCS模型构建的特征值，υn1、υn2、υn3为三维结构张量所对应的三个特征向量。
7-5，依据公式(19)至(20)，进行迭代计算并获取最终的增强结果。
代表图像梯度；u表示归一化后的原始图像；div表示散度算子，*符号代表卷积，·符号代表点乘。同时M表示空间内核，内核值p及排列方式由Kroon等人在文献《Optimized anisotropic rotational invariant diffusion scheme on cone-beam CT》中提出的数值优化方法计算得出。迭代计算方程由公式(20)所示，其中k表示迭代次数，τ表示迭代时间步长。
步骤8，利用步骤6与步骤7的计算结果(步骤6中，利用Pock管状响应函数得到的滤波结果；步骤7中利用HDCS扩散滤波得到的管状结构增强结果)，依据公式(21)构建区域描述算子k(x)，并用于最小化能量分割VRG算法。
P(x)表示体素点x的管状响应，apin表示被判定为血管区域体素点的平均管状响应。对应的apout表示非血管区域体素点的平均管状响应，MP表示掩膜Mask内的最大管状响应。以此类推，T(x)表示体素点x增强后的响应，与分别表示被划分为血管区域与非血管区域体素点的平均增强响应值。MT表示掩膜Mask内的最大增强响应值。
步骤9，利用步骤8中构建的区域描述算子以及初始化的分割区域，在掩膜Mask内利用VRG算法(文献《VARIATIONAL REGION GROWING》中的变化区域生长算法)对肺部血管进行精细分割：
9-1：求取Pock函数管状响应结果的局部极大值区域，利用局部极大值区域作为初始分割种子集Seeds，并设置种子点集对应体素点的初始状态φ(x)为1。
9-2：从种子集Seeds开始，对内外边界点依据公式(22)和(23)进行能量变化计算，从而纳入或排除血管区域内的体素点。分割进程将持续至能量状态稳定即∑φn+1(x)＝∑φn(x)。
以Pock函数计算结果的局部最大值为种子点集合S。
VRG算法纳入或排除体素点的能量方程为：
ΔJ(φn+1(x))＝(1-2φn(x))k(x) (22)
ΔJ代表能量变化项，φ(x)代表体素点x的状态值，它的表达式如公式(23)所示，n+1与n表示迭代次数：
Ωin表示血管区域，Ωout表示非血管区域。因此，依据能量最小化原则，仅当体素点的能量变化值为负值时，该体素点将被纳入(x初始状态值为0)或排除出血管区域(x初始状态值为1)。
本实验采用了国际血管分割竞赛VESsel Segmentation in the Lung 2012(VESSEL12)中的图像数据，每幅断层图像都是512×512的16位DICOM格式图像。用于本发明实施例的竞赛数据可在官网上下载获取，并且该分割实验是在Matlab R2015a与VisualStudio 2013平台下进行的。实验环境为：Windows 8.1，4核Intel(R)i5-3470CPU 3.20GHz，且每个分割案例的分割时间不超过40分钟。
图2是本发明针对VESSEL12竞赛案例采用管状结构增强与能量函数结合的肺部血管树分割方法分割所得的结果图。如图2所示，对上述肺部CT序列断层图像，采用本发明所述的方法对肺部血管树进行了分割。从图2中可以看出，本发明的方法在分割出主支血管的同时，分割出了丰富的细小血管区域。
为了定量地评价本发明分割方法的效果，本发明所述方法分割获取的二值结果被上传至官方网站。通过定制的统一评价标准，竞赛主办方对本方法的分割结果进行了定量分析，并将分析结果在官方网站上发布。
本发明方法针对VESSEL12案例分割结果的定量分析结果，图中的定量评价标准由竞赛官方制定，并且评价标准与分析结果都可在竞赛官方网站查询得出。
利用本发明方法对20组肺部血管树进行了分割，并通过上传分割结果的方式，对本发明方法分割结果进行了全面的评估。从定量分析结果的总评分可得出，ROC曲线下面积指标AUC＝0.897，此评分在竞赛现有的二值结果中，排名第三。同时可看出本发明方法的敏感性为0.860，特异性为0.965。本发明方法在20组数据中表现平稳，AUC最低评分为0.863，最高评分为0.944，且能够排除气管壁对血管分割结果的影响。因此，使用本发明方法可以有效、准确地分割出肺部血管树。
所发明涉及的管状结构增强与能量函数结合的肺部血管树分割方法，因为将Pock管状响应函数与HDCS增强方法相互结合，并利用VRG算法对肺部血管进行精细分割。此方法能够检测出更多的潜在肺部血管区域，并受噪声的影响较小，分割结果较为准确。而且本实施例的管状结构增强与能量函数结合的肺部血管树分割方法能够分割出更多的细小血管区域。且本方法特异性高，敏感性较强，同时能够区分血管与气管壁区域，进一步提高了分割结果的准确性。本实施例所涉及的方法在VESSEL12竞赛20组案例中都有较好的表现。因此证明本发明涉及的管状结构增强与能量函数结合的肺部血管树分割方法能够有效、准确的分割出肺部血管树。
Claims (3)
1.一种管状结构增强与能量函数结合的肺部血管树分割方法，其特征在于，具体包括如下步骤：
1)输入待分割DICOM格式的胸部CT序列断层图像；
2)利用文献《Automatic Lung Segmentaion for Accurate Quantitaion ofVolumetric X-Ray CT Images》中的阈值法与形态学修补法，对步骤1)的胸部CT序列断层图像肺部区域进行分割并获取掩膜Mask1；
3)利用文献《Two-pass region growing combined morphology algorithm forsegmenting airway tree from CT chest scans》中的双路径区域生长结合形态学重建的方法对步骤1)的胸部CT序列断层图像中气管树进行分割，获取不含气管壁的气管树；
4)对步骤3中分割得到的肺部气管树，利用多尺度球形结构元素进行形态学开操作，从而剥离出含有不同等级分枝的气管树；然后对剥离出的各个气管树，使用对应尺度的球形结构元素进行膨胀操作，从而获取含有气管壁区域的气管树Airway；
5)从步骤2)所得肺部区域掩膜Mask1中去除气管树Airway对应的区域，从而获取肺部血管分割掩膜Mask；
6)在分割掩膜Mask区域内，利用Pock管状响应函数，对步骤1)中的图像进行多尺度管状响应计算，并获取管状响应结果；
6-1)设定尺度集合σi，σmin≤σi≤σmax，其中最小尺度σmin等于0.5，最大尺度σmax等于8，尺度间隔为0.5，同时设定Pock响应函数对称敏感系数σw为0.2，噪声抑制项δ设置为0.2；
6-2)对步骤1)中的图像进行高斯核标准差大小为σi的高斯滤波，i＝0.5，1，1.5，......，8，然后求取滤波结果体素点的Hessian矩阵，用求取的Hessian矩阵乘以然后计算Hessian矩阵的特征值与特征向量；
6-3)利用Pock管状响应函数公式(1)至(8)，计算在尺度σ下的管状响应；
其中R+(x，σ，θ)为对称性约束的管状结构响应；x代表图像中的任意体素点；σ代表当前尺度；θ代表待检测管状结构半径与尺度σ间的比例系数；R-(x，σ)代表体素点x在尺度σ下的梯度幅值；δ代表噪声抑制项：
N为离散圆周点个数，w(bi)为对称系数，bi为原始边界响应，
I表示步骤1)中胸部CT序列断层图像；bi表示在尺度σ下，第i个圆周点处的梯度幅值，即边界响应，它表示了该点属于边界点的可能性大小：表示第i个圆周点对应的旋转相量，用于估算圆周点的位置；
V1和V2分别表示Hessian矩阵主曲率方向对应的两个特征向量，对应的特征值关系为|λ1|＞|λ2|；αi代表第i个圆周点的离散相量角；
其中为圆周点的平均边界响应；σω为检测对称性的敏感系数，取值范围为(0，1]；
6-4)利用管状响应规范化公式，对尺度σi下的管状响应进行响应规范化，规范化公式如下所示：
Rnorm(x，σi，θ)＝σi γR(x，σi，θ)
其中Rnorm(x，σi，θ)代表规范化后的管状响应，R(x，σi，θ)表示未经规范化的管状响应，γ代表规范化系数，取值为1：
6-5)重复子步骤6-2)至6-4)以求取各尺度下的规范化管状响应，并利用多尺度管状响应公式求取最终的管状响应，多尺度管状响应公式如下所示：
7)使用基于三维扩散张量的图像增强方法，对Mask掩膜内的区域进行增强；
8)利用步骤6)与步骤7)的计算结果，依据公式(21)构建区域描述算子k(x)；
P(x)表示体素点x的管状响应；表示被判定为血管区域体素点的平均管状响应；对应的表示非血管区域体素点的平均管状响应；MP表示掩膜Mask内的最大管状响应；以此类推，T(x)表示体素点x增强后的响应，与分别表示被划分为血管区域与非血管区域体素点的平均增强响应值；MT表示掩膜Mask内的最大增强响应值；
9)利用步骤8)中构建的区域描述算子以及初始化的分割区域，在掩膜Mask内利用VRG算法对肺部血管进行精细分割：
9-1)求取Pock函数管状响应结果的局部极大值区域，利用局部极大值区域作为初始分割种子集Seeds，并设置种子点集对应体素点的初始状态φ(x)为1；
9-2)从种子集Seeds开始，对内外边界点依据公式(22)和(23)进行能量变化计算，从而纳入或排除血管区域内的体素点，分割进程将持续至能量状态稳定即∑φn+1(x)＝∑φn(x)；
VRG算法纳入或排除体素点的能量方程为：
ΔJ(φn+1(x))＝(1-2φn(x))k(x) (22)
ΔJ代表能量变化项，φ(x)代表体素点x的状态值，它的表达式如公式(23)所示，n+1与n表示迭代次数：
Ωin表示血管区域，Ωout表示非血管区域，依据能量最小化原则，仅当体素点的能量变化值为负值时，该体素点将被纳入或排除出血管区域。
2.根据权利要求1所述管状结构增强与能量函数结合的肺部血管树分割方法，其特征在于，所述步骤4)实现的具体步骤如下：
4-1)设置开操作结构元素半径ri o＝2(i+1)，其中半径为毫米单位，共取七个尺度；
4-2)利用具有不同半径ri o的结构元素对步骤3)中获取的气管树进行形态学开操作，从而剥离出含有不同等级分枝的气管树集合
4-3)设置膨胀操作结构元素的半径为开操作结构元素的二分之一，即ri d＝i+1，
4-4)对步骤4-2)中获取的气管树集合利用对应的膨胀操作结构元素进行形态学膨胀操作，从而获取气管树集合将所得气管树集合进行并集，从而获取包含气管壁区域的气管树Airway。
3.根据权利要求1或2所述管状结构增强与能量函数结合的肺部血管树分割方法，其特征在于，所述步骤7)实现的具体步骤如下：
7-1)设置实验参数，设置ρ＝1以及σ＝1分别用于归一化图像的高斯平滑滤波与平滑三维结构张量的高斯平滑滤波；同时设置扩散模型相关参数，C＝3.31488，λc＝0.02，λe＝O.02，λh＝0.5。迭代的次数为5，步长为0.5；
7-2)利用高斯核标准差为1的高斯滤波器对归一化图像进行滤波，并对滤波结果进行求导：
对归一化的图像进行高斯滤波，然后对滤波后的图像分别使用模板与在三个方向上进行滤波，并获取滤波结果ux，uy和uz。其中三个模板都是3*3*3的三维矩阵，对应每一层的二维矩阵表达式分别如公式(9)，(10)以及(11)所示；
公式(9)中，表示x方向模板第一层矩阵表达式，为第二层矩阵表达式，为第三层表达式，以此类推；
7-3)依据公式(12)构建三维结构张量，对三维结构张量进行高斯核标准差为1的高斯滤波，然后计算其对应的特征值与特征向量；
7-4)依据公式(13)至(18)，构建基于连续转换的混合扩散HDCS模型的三维扩散张量：
表示一致增强扩散CED模型；表示边缘增强扩散EED模型；ε为权重值，用于控制EED模型与CED模型的连续转化；
α为经验值取值为0.001；λc为管状结构对比参数取值为0.02；κ＝(μ2/(α+μ3))4；代表图像梯度；则代表σ尺度下的图像梯度；μ2与μ3代表三维结构张量的两个特征值；μ1，μ2与μ3分别代表三维结构张量的三个特征值，且三个特征值的绝对值大小关系为：μ3≥μ2≥μ1；
C与m分别为阈值参数与经验取值。C取值为3.31488，m取值为4。λe为平面结构对比参数，并且取经验值为0.02，
λn为利用HDCS模型构建的特征值，vn1、vn2、vn3为三维结构张量所对应的三个特征向量；
7-5)依据公式(19)至(20)，进行迭代计算并获取最终的增强结果：
du1＝u*Mx·(D11*Mx+D12*My+D13*Mz)
du2＝u*My·(D12*Mx+D22*My+D23*Mz)
du2＝u*Mz·(D13*Mx+D23*My+D33*Mz)
u11＝u*Mxx
u22＝u*Myy
u33＝u*Mzz
u12＝u*Mxy
u13＝u*Mxz
代表图像梯度；u表示归一化后的原始图像；div表示散度算子，*符号代表卷积，·符号代表点乘；同时M表示空间内核，内核值p及排列方式由Kroon等人在文献《Optimizedanisotropic rotational invariant diffusion scheme on cone-beam CT》中提出的数值优化方法计算得出；迭代计算方程由公式(20)所示，其中k表示迭代次数，τ表示迭代时间步长；
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810222634.XA CN108492300B (zh) | 2018-03-16 | 2018-03-16 | 管状结构增强与能量函数结合的肺部血管树分割方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810222634.XA CN108492300B (zh) | 2018-03-16 | 2018-03-16 | 管状结构增强与能量函数结合的肺部血管树分割方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108492300A true CN108492300A (zh) | 2018-09-04 |
| CN108492300B CN108492300B (zh) | 2021-07-13 |
Family
ID=63339721
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810222634.XA Expired - Fee Related CN108492300B (zh) | 2018-03-16 | 2018-03-16 | 管状结构增强与能量函数结合的肺部血管树分割方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108492300B (zh) |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109431531A (zh) * | 2018-12-25 | 2019-03-08 | 上海联影医疗科技有限公司 | 基于灌注成像的血管分割方法及装置和计算机装置 |
| CN109472807A (zh) * | 2018-11-30 | 2019-03-15 | 北京师范大学 | 基于深度神经网络的血管模型提取方法 |
| CN109584223A (zh) * | 2018-11-20 | 2019-04-05 | 北京中科研究院 | Ct图像中肺部血管分割方法 |
| CN109584167A (zh) * | 2018-10-24 | 2019-04-05 | 深圳市旭东数字医学影像技术有限公司 | 基于二阶特征的ct图像肝内血管增强与分割方法及系统 |
| CN111080556A (zh) * | 2019-12-23 | 2020-04-28 | 山东师范大学 | 一种强化ct图像气管壁增强方法、系统、设备及介质 |
| CN111626974A (zh) * | 2019-02-28 | 2020-09-04 | 苏州润迈德医疗科技有限公司 | 冠状动脉造影图像序列的质量评分方法和装置 |
| CN112184659A (zh) * | 2020-09-24 | 2021-01-05 | 上海健康医学院 | 一种肺部图像处理方法、装置及设备 |
| CN112651933A (zh) * | 2020-12-21 | 2021-04-13 | 山东省人工智能研究院 | 基于测地线距离图和程函方程的血管分割方法 |
| CN114240970A (zh) * | 2021-12-21 | 2022-03-25 | 北京适创科技有限公司 | 一种用于ct数据的自动摆正与自动强化方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102324109A (zh) * | 2011-09-26 | 2012-01-18 | 上海理工大学 | 基于模糊隶属度模型的非实质性肺结节三维分割方法 |
| CN104616307A (zh) * | 2015-02-12 | 2015-05-13 | 河北大学 | 一种肺部ct图像粘连血管型结节检测方法 |
| CN106097305A (zh) * | 2016-05-31 | 2016-11-09 | 上海理工大学 | 双行程区域生长结合形态学重建的肺部气管树分割方法 |
-
2018
- 2018-03-16 CN CN201810222634.XA patent/CN108492300B/zh not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102324109A (zh) * | 2011-09-26 | 2012-01-18 | 上海理工大学 | 基于模糊隶属度模型的非实质性肺结节三维分割方法 |
| CN104616307A (zh) * | 2015-02-12 | 2015-05-13 | 河北大学 | 一种肺部ct图像粘连血管型结节检测方法 |
| CN106097305A (zh) * | 2016-05-31 | 2016-11-09 | 上海理工大学 | 双行程区域生长结合形态学重建的肺部气管树分割方法 |
Non-Patent Citations (3)
| Title |
|---|
| CHRISTIANBAUER等: "Segmentation of interwoven 3d tubular tree structures utilizing shape priors and graph cuts", 《MEDICAL IMAGE ANALYSIS》 * |
| 操时力: "基于CT图像的肝脏门静脉分割方法研究", 《中国优秀硕士学位论文全文数据库》 * |
| 步蕊蕊等: "基于三维结构张量的CT肺血管树增强", 《中国医学物理学杂志》 * |
Cited By (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109584167A (zh) * | 2018-10-24 | 2019-04-05 | 深圳市旭东数字医学影像技术有限公司 | 基于二阶特征的ct图像肝内血管增强与分割方法及系统 |
| CN109584223A (zh) * | 2018-11-20 | 2019-04-05 | 北京中科研究院 | Ct图像中肺部血管分割方法 |
| CN109472807A (zh) * | 2018-11-30 | 2019-03-15 | 北京师范大学 | 基于深度神经网络的血管模型提取方法 |
| CN109472807B (zh) * | 2018-11-30 | 2021-11-26 | 北京师范大学 | 基于深度神经网络的血管模型提取方法 |
| CN109431531A (zh) * | 2018-12-25 | 2019-03-08 | 上海联影医疗科技有限公司 | 基于灌注成像的血管分割方法及装置和计算机装置 |
| CN111626974A (zh) * | 2019-02-28 | 2020-09-04 | 苏州润迈德医疗科技有限公司 | 冠状动脉造影图像序列的质量评分方法和装置 |
| CN111626974B (zh) * | 2019-02-28 | 2024-03-22 | 苏州润迈德医疗科技有限公司 | 冠状动脉造影图像序列的质量评分方法和装置 |
| CN111080556B (zh) * | 2019-12-23 | 2023-06-13 | 山东师范大学 | 一种强化ct图像气管壁增强方法、系统、设备及介质 |
| CN111080556A (zh) * | 2019-12-23 | 2020-04-28 | 山东师范大学 | 一种强化ct图像气管壁增强方法、系统、设备及介质 |
| CN112184659A (zh) * | 2020-09-24 | 2021-01-05 | 上海健康医学院 | 一种肺部图像处理方法、装置及设备 |
| WO2022063198A1 (zh) * | 2020-09-24 | 2022-03-31 | 上海健康医学院 | 一种肺部图像处理方法、装置及设备 |
| CN112184659B (zh) * | 2020-09-24 | 2023-08-25 | 上海健康医学院 | 一种肺部图像处理方法、装置及设备 |
| CN112651933A (zh) * | 2020-12-21 | 2021-04-13 | 山东省人工智能研究院 | 基于测地线距离图和程函方程的血管分割方法 |
| CN114240970A (zh) * | 2021-12-21 | 2022-03-25 | 北京适创科技有限公司 | 一种用于ct数据的自动摆正与自动强化方法 |
| CN114240970B (zh) * | 2021-12-21 | 2024-05-24 | 北京适创科技有限公司 | 一种用于ct数据的自动摆正与自动强化方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108492300B (zh) | 2021-07-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108492300A (zh) | 管状结构增强与能量函数结合的肺部血管树分割方法 | |
| Sharma et al. | Dermatologist-level classification of skin cancer using cascaded ensembling of convolutional neural network and handcrafted features based deep neural network | |
| Kooi et al. | Large scale deep learning for computer aided detection of mammographic lesions | |
| CN108364006B (zh) | 基于多模式深度学习的医学图像分类装置及其构建方法 | |
| Liu et al. | Automatic organ segmentation for CT scans based on super-pixel and convolutional neural networks | |
| Sori et al. | Multi-path convolutional neural network for lung cancer detection | |
| CN105741251B (zh) | 一种肝脏cta序列图像的血管分割方法 | |
| CN103714536B (zh) | 基于稀疏表示的多模态磁共振图像的分割方法及装置 | |
| CN101706843B (zh) | 一种乳腺cr图像交互式读片方法 | |
| Tian et al. | Multi-path convolutional neural network in fundus segmentation of blood vessels | |
| CN109389584A (zh) | 基于cnn的多尺度鼻咽肿瘤分割方法 | |
| Panda et al. | New binary Hausdorff symmetry measure based seeded region growing for retinal vessel segmentation | |
| CN105139377B (zh) | 一种腹部ct序列图像肝脏的快速鲁棒自动分割方法 | |
| CN109325942A (zh) | 基于全卷积神经网络的眼底图像结构分割方法 | |
| Liu et al. | Automatic whole heart segmentation using a two-stage u-net framework and an adaptive threshold window | |
| CN110265095A (zh) | 用于hcc复发及rfs的预测模型和诺模图的构建方法及应用 | |
| CN112270666A (zh) | 基于深度卷积神经网络的非小细胞肺癌病理切片识别方法 | |
| CN104299242B (zh) | 基于ngc‑acm的荧光造影眼底图像提取方法 | |
| CN108038513A (zh) | 一种肝脏超声图像的特征分类方法 | |
| CN101714153A (zh) | 基于视觉感知的交互式乳腺钼靶图像检索方法 | |
| Jony et al. | Detection of lung cancer from CT scan images using GLCM and SVM | |
| Mahapatra et al. | A novel framework for retinal vessel segmentation using optimal improved frangi filter and adaptive weighted spatial FCM | |
| Mahapatra et al. | Visual saliency based active learning for prostate mri segmentation | |
| CN112308846B (zh) | 血管分割方法、装置及电子设备 | |
| de Sousa Costa et al. | Classification of malignant and benign lung nodules using taxonomic diversity index and phylogenetic distance |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210713 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |