CN114119626A - 一种基于统计学模型和多尺度滤波的脑血管图像分割方法 - Google Patents
一种基于统计学模型和多尺度滤波的脑血管图像分割方法 Download PDFInfo
- Publication number
- CN114119626A CN114119626A CN202111210464.1A CN202111210464A CN114119626A CN 114119626 A CN114119626 A CN 114119626A CN 202111210464 A CN202111210464 A CN 202111210464A CN 114119626 A CN114119626 A CN 114119626A
- Authority
- CN
- China
- Prior art keywords
- blood vessel
- point
- denoised
- image
- dimensional 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 72
- 238000001914 filtration Methods 0.000 title claims abstract description 32
- 238000003709 image segmentation Methods 0.000 title claims abstract description 32
- 238000013179 statistical model Methods 0.000 title claims abstract description 30
- 210000004204 blood vessel Anatomy 0.000 claims abstract description 209
- 230000011218 segmentation Effects 0.000 claims abstract description 36
- 210000004556 brain Anatomy 0.000 claims abstract description 30
- 238000012545 processing Methods 0.000 claims abstract description 15
- 238000007781 pre-processing Methods 0.000 claims abstract description 5
- 239000000203 mixture Substances 0.000 claims description 27
- 230000009191 jumping Effects 0.000 claims description 18
- 238000004364 calculation method Methods 0.000 claims description 16
- 210000003625 skull Anatomy 0.000 claims description 13
- 238000005381 potential energy Methods 0.000 claims description 12
- 238000003384 imaging method Methods 0.000 claims description 8
- 230000036961 partial effect Effects 0.000 claims description 8
- 239000011159 matrix material Substances 0.000 claims description 7
- 241001632427 Radiola Species 0.000 claims description 6
- 210000001519 tissue Anatomy 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 5
- 208000032843 Hemorrhage Diseases 0.000 claims description 4
- 238000007476 Maximum Likelihood Methods 0.000 claims description 3
- 210000005252 bulbus oculi Anatomy 0.000 claims description 3
- 230000004087 circulation Effects 0.000 claims description 3
- 230000009850 completed effect Effects 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 3
- 238000005315 distribution function Methods 0.000 claims description 3
- 230000002708 enhancing effect Effects 0.000 claims description 3
- 230000002401 inhibitory effect Effects 0.000 claims description 3
- 230000003313 weakening effect Effects 0.000 claims description 3
- 238000012512 characterization method Methods 0.000 claims description 2
- 239000000919 ceramic Substances 0.000 claims 1
- 238000000605 extraction Methods 0.000 abstract description 10
- 230000001502 supplementing effect Effects 0.000 abstract description 3
- 238000012216 screening Methods 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 44
- 208000026106 cerebrovascular disease Diseases 0.000 description 10
- 201000010099 disease Diseases 0.000 description 10
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 10
- 230000000694 effects Effects 0.000 description 9
- 238000010586 diagram Methods 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 238000002372 labelling Methods 0.000 description 4
- 208000006011 Stroke Diseases 0.000 description 3
- 210000005013 brain tissue Anatomy 0.000 description 3
- 230000036541 health Effects 0.000 description 3
- 230000002265 prevention Effects 0.000 description 3
- 230000002829 reductive effect Effects 0.000 description 3
- 238000007619 statistical method Methods 0.000 description 3
- 238000011282 treatment Methods 0.000 description 3
- 208000032382 Ischaemic stroke Diseases 0.000 description 2
- 230000032683 aging Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000017531 blood circulation Effects 0.000 description 2
- 210000000988 bone and bone Anatomy 0.000 description 2
- 230000034994 death Effects 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 208000028867 ischemia Diseases 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000000926 separation method Methods 0.000 description 2
- 230000006438 vascular health Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 206010002329 Aneurysm Diseases 0.000 description 1
- 208000024172 Cardiovascular disease Diseases 0.000 description 1
- 208000007530 Essential hypertension Diseases 0.000 description 1
- 102000009123 Fibrin Human genes 0.000 description 1
- 108010073385 Fibrin Proteins 0.000 description 1
- BWGVNKXGVNDBDI-UHFFFAOYSA-N Fibrin monomer Chemical compound CNC(=O)CNC(=O)CN BWGVNKXGVNDBDI-UHFFFAOYSA-N 0.000 description 1
- 208000016988 Hemorrhagic Stroke Diseases 0.000 description 1
- 206010020772 Hypertension Diseases 0.000 description 1
- 208000012902 Nervous system disease Diseases 0.000 description 1
- 208000025966 Neurological disease Diseases 0.000 description 1
- 206010033799 Paralysis Diseases 0.000 description 1
- 208000007536 Thrombosis Diseases 0.000 description 1
- 206010068268 Tumour compression Diseases 0.000 description 1
- 208000009443 Vascular Malformations Diseases 0.000 description 1
- 206010057469 Vascular stenosis Diseases 0.000 description 1
- 206010047115 Vasculitis Diseases 0.000 description 1
- 230000001154 acute effect Effects 0.000 description 1
- 210000001367 artery Anatomy 0.000 description 1
- 239000008280 blood Substances 0.000 description 1
- 210000004369 blood Anatomy 0.000 description 1
- 230000036995 brain health Effects 0.000 description 1
- 230000002490 cerebral effect Effects 0.000 description 1
- 210000001175 cerebrospinal fluid Anatomy 0.000 description 1
- 208000010877 cognitive disease Diseases 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000002526 effect on cardiovascular system Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 210000001508 eye Anatomy 0.000 description 1
- 229950003499 fibrin Drugs 0.000 description 1
- 210000004884 grey matter Anatomy 0.000 description 1
- 230000008076 immune mechanism Effects 0.000 description 1
- 208000020658 intracerebral hemorrhage Diseases 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000004660 morphological change Effects 0.000 description 1
- 230000017074 necrotic cell death Effects 0.000 description 1
- 230000000926 neurological effect Effects 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 230000003449 preventive effect Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000011321 prophylaxis Methods 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 239000013589 supplement Substances 0.000 description 1
- 208000024891 symptom Diseases 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 238000011277 treatment modality Methods 0.000 description 1
- 230000002792 vascular Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
- 210000004885 white matter Anatomy 0.000 description 1
Images
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
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
- G06N3/084—Backpropagation, e.g. using gradient descent
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration using local operators
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration using local operators
- G06T5/30—Erosion or dilatation, e.g. thinning
-
- 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
- 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/30016—Brain
-
- 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)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computational Linguistics (AREA)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Health & Medical Sciences (AREA)
- Evolutionary Computation (AREA)
- General Health & Medical Sciences (AREA)
- Molecular Biology (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Image Analysis (AREA)
Abstract
本发明提出了一种基于统计学模型和多尺度滤波的脑血管图像分割方法，属于医学计算及图像处理技术领域。所述分割方法将统计模型和多尺度滤波的相结合，以在高层与低层模型相结合提取脑血管图像并分割，包括：S1，数据预处理，得到去噪后脑部区域图像；S2，拟合S1图像背景及血管类的高斯模型，搭建FMM；S3，设置FMM初始值，估计参数及计算低层能量函数；S4，对图像进行分割；S5，结合脊线性质对种子点进行筛选，在3D空间中提取血管中心线，沿远端血管中心线进行膨胀，对低层模型进行补充。所述方法提高了远端血管图像的分割精度，能实现在较为粗糙的血管图像分割结果下得到精准的血管中心线提取。
Description
技术领域
本发明涉及一种基于统计学模型和多尺度滤波的脑血管图像分割方法，属于医学计算及图像处理技术领域。
背景技术
血管狭窄、动脉瘤和血管畸形等脑血管疾病是威胁人类健康的重大疾病。随着我国老龄化速度加快，老年脑血管病人数量也在逐年上升，及早发现脑血管形态变化，对脑血管疾病进行有效预防和诊断极为重要。基于脑血管影像的血管数字化分割提取是其中的关键技术。
脑血管病指的是在大脑中部分血管由于病理原因无法正常运作，使得脑组织缺血而受损的疾病，比较典型的有免疫机制相关的血管炎、肿瘤压迫等导致的血栓、高血压造成的纤维蛋白坏死等，这些疾病大多能够引发病患脑部组织出现缺血或出血，最终导致患者死亡或残疾。据调查，这些疾病的临床患者以中老年人居多，且大多是急性发病。“中风”就是一种脑血管病，临床称之为脑卒中。临床将该疾病划分为两类：出血性卒中与缺血性卒中。调查发现，该类疾病是导致我国成年人残疾甚至死亡的最常见原因。该疾病具有发病率高、复发率高、致残率以及死亡率较高的特点，会对病人家庭造成较重的经济负担，其症状常表现为瘫痪、言语认知出现障碍以及抑郁症等。患病率以及发病率相比于2015年也在持续上升，对于人类的健康造成了巨大的威胁。上文提到的两种脑卒中，治疗方式并不相同。考虑到当前有效治疗方式的缺乏，预防是最为合适的应对措施。由此，对于缺血性卒中等脑血管疾病的相应治疗措施而言，是否能够将完整血管进行三维可视化显得至关重要。
目前基于深度学习方法进行的血管图像分割已经取得了较好的效果，但其依赖大量的标注数据进行训练，同时需要标准的数据采集规范，限制了模型在临床上的推广应用。而统计分析模型依据像素灰度建立有限混合模型，提取血管过程中无需标注数据，对于小样本的医学数据具有优势，但受其模型限制对于微小血管图像分割效果不理想。本申请致力于解决上述问题。
发明内容
本发明的目的是为解决由于依据像素灰度建立有限混合模型的统计分析模型对于高强度区域能有效分割，但对于低强度的小血管图像分割的效果较差这一技术缺陷，提出了一种基于统计学模型和多尺度滤波的脑血管图像分割方法，该方法基于脊线的3D中心线提取，将统计模型和多尺度滤波相结合，通过高层与低层模型相结合实现对脑血管图像提取，具体为：利用统计学模型在高强度像素点上分割提取，结合脊线性质对种子点进行筛选；应用脊线性质，在3D空间中提取血管中心线，沿远端血管中心线进行膨胀，对低层模型结果进行补充，实现了多尺度滤波的脑血管图像分割。
为了达到上述目的，本发明采取如下技术方案。
所述脑血管图像分割方法，包括如下步骤：
步骤1：数据预处理，得到去噪后脑部区域图像，具体为：采集大脑部分血管原始图像、并对得到的血管三维图像进行去噪，得到去噪后血管三维图像，并根据去噪后血管三维图像构建血管三维图像的灰度的分布直方图；具体为：
步骤1.1：通过TOF-MRA及SNAP-MRA采集出大脑部分血管的原始图像；
其中，TOF-MRA，即飞行时间磁共振血管成像；SNAP-MRA即非对比增强血管与斑块内出血同时成像；血管三维图像即颅骨剥离后的脑部区域图像；
步骤1.2：使用FSL软件的bet命令对步骤1.1采集的原始图像进行颅骨剥离，将原始数据去除眼球脂肪等不相干的图像，得到血管三维图像；
其中，步骤1.2中对原始图像进行处理的原因是：TOF-MRA及SNAP-MRA采集时由于设备限制，无法单独采集脑部区域；
步骤1.3：选用大小为n*n*n，方差为se的高斯核对血管三维图像颅骨剥离后的脑部区域图像进行高斯滤波,滤除图像中的无关噪声，得到去噪后脑部区域图像，n的取值范围是2到10，se的取值范围是0.2-1.0，具体包括如下子步骤：
步骤1.3.1：构建血管三维图像的黑塞矩阵，其中，血管三维图像中每一个像素点的一阶导数表示灰度的梯度，二阶导数是灰度梯度的梯度，即灰度的曲率；
步骤1.3.2：求出血管三维图像中每一个像素点的特征向量后，将该像素点的曲率进行分解；判断出当前像素点的曲率主方向及该方向对应的照其绝对值大小顺序排列的特征向量λ1，λ2和λ3；
步骤1.3.3：建立针对某一像素点(x,y,z)的高斯卷积核G(x,y,z)：
步骤1.3.4：将高斯卷积核G(x,y,z)与对应像素点(x,y,z)的体素强度 I(x,y,z)进行卷积得到，在尺度σ下的体素强度I′(x，y，z，σ)；
步骤1.3.5重复步骤1.3.3到1.3.4，得到在尺度σ下血管三维图像中所有像素点的体素强度；
步骤1.3.6：使用Hessian滤波器检测出血管三维图像中的管状结构，再对管状结构进行增强，对其余结构信号进行抑制减弱，以增强血管组织的对比度，得到增强后的血管三维图像，具体为：如果特征值λ2和λ3中都小于0，那么我们引入超参数α用于区分板状结构与管状结构；β用于区分斑点状结构与管状结构；c是用于区分板状结构与斑点状结构，计算对于某一像素点s在尺度σ下增强之后的像素表征V0(s，σ)；
步骤1.3.7：重复步骤1.3.6，在每一个尺度下σ，进行计算该尺度时各个像素点计算得到的像素表征V0(s，σ)；
步骤1.3.8：将步骤1.3.7中各个尺度下得到的所有的V0(s，σ)中的最大值作为去噪后血管三维图像该像素点的像素值V(s)；
其中，去噪后血管三维图像即去噪后脑部区域图像；
步骤1.3.9：根据步骤1.3.8得到的去噪后的血管三维图像的像素值V(s)的灰度分布情况构建灰度分布直方图Iinit；
步骤2：通过分别拟合去噪后脑部区域图像背景及血管类的高斯模型，搭建有限混合模型FMM，具体为：
步骤2.1：使用两个高斯模型G1(x)、G2(x)拟合步骤1.3.8得到的去噪后血管三维图像中的背景类；
步骤2.2：使用一个高斯模型G3(x)拟合步骤1.3.8得到的去噪后血管三维图像中的血管类；
步骤2.3：根据步骤2.1与步骤2.2拟合出来的三个高斯模型G1(x)、G2(x) 以及G3(x)搭建有限混合模型FMM，其FMM公式描述的灰度分布his()是步骤2.1与步骤2.2拟合出来的三个高斯模型的线性组合，his()利用超参数α1，α2，α3调节三个高斯模型G1(x)、G2(x)以及G3(x)间各自的权重；
步骤3：使用最大期望方法对有限混合模型FMM设置初始值，对初始值进行参数估计，并计算低层能量函数，具体为：为有限混合模型设置初始值，并对有限混合模型执行最大期望方法，估计该有限混合模型FMM设置初始值，具体为：
步骤3.1：为有限混合模型设置初始值，得到三个高斯模型的参数集θ，具体包括如下子步骤：
步骤3.1.1：初始化步骤2.2拟合得到的血管类的高斯分布G3(x)，将其在步骤2.3中得到的his(x)表达式中的占比α3设定为区间0.01～0.05中的常数值；
步骤3.1.2：使用经过步骤1.3.8中去噪的血管三维图像中最亮的前M％的像素点计算得出经过步骤2.2拟合出来的高斯分布函数G3(x)的初始均值μ3及方差σ32；
其中，M为取值1-5之间的常数；
步骤3.1.3：将经过步骤2.1拟合出来的G2(x)的初始均值μ2设定为步骤 1.3.9中得到的灰度直方图Iinit的最高点的灰度；
步骤3.1.4：将经过步骤2.1拟合出来的G1(x)的初始均值μ1设立为G2(x)的初始均值μ2的t倍；
其中，t为常数，取值范围为0.1-0.5；
步骤3.1.5：计算G1(x)贡献在步骤1.3.9中的得到的灰度直方图中Iinit占比作为α1的估计值；
步骤3.1.6：根据步骤3.1.5计算得到的G1(x)的体素占比α1，计算中去掉 G1(x)贡献后的剩余分布直方图Ileft；
步骤3.1.7：取步骤3.1.6得到的剩余分布直方图Ileft在区间[(μ2-μ1)/ 2，(μ2+μ1)/2]上的灰度点进行最大似然估计计算G2(x)函数的方差σ22；
步骤3.1.8：根据步骤3.1.1的求解的α3和步骤3.1.5的求解的α1，G2(x) 的占比α2使用1-α1-α3得到；
均值参数的集合θμ、方差参数的的集合θσ及体素占比参数的集合θα，统称为参数集θ；
步骤3.2：对有限混合模型执行最大期望方法，具体步骤如下：
步骤3.2.1：设定循环参数t＝1；
步骤3.2.2：根据三个高斯模型的参数集θ，计算隐变量γ其联合概率的条件概率P；
步骤3.2.3：根据步骤3.2.2计算出来的隐变量γ其联合概率的条件概率P，计算期望E；
步骤3.2.4：根据步骤3.2.3计算得到的期望E，计算使得其最大的参数取值，作为t+1次循环的参数集，将体素占比的和，即∑θα设为1，得到新拉格朗日函数；
步骤3.2.5：使用步骤3.2.4得到的新拉格朗日函数对λ、αz、μz、σz求偏导，令其偏导数为0，得到新的参数集；
步骤3.2.6：判断t是否大于等于2，若是，则判断第t次及第t-1两次循环的所有参数的变化是否均小于q％，若是，将步骤3.2.5得到的新的参数集作为 his(ys)函数的初始参数的估计，并跳至步骤4；否则判断t是否等于N，若是，将步骤3.2.5得到的新的参数集作为his(ys)函数的初始参数的估计，并跳至步骤4，否则将t加1，跳至步骤3.2.2；
其中，q的取值范围在0.5-1.5之间，N的取值范围为25到75；
至此，完成了对his(ys)函数的初始参数θ的估计；
步骤4：使用步骤2搭建的有限混合模型FMM与步骤1中去噪后血管三维图像的像素值，计算得到去噪的血管三维图像分割的初步结果，具体步骤如下：
步骤4.1：通过步骤3得到的his(ys)函数的初始参数θ的估计，提出低层能量函数U1(xs)满足U1(xs)∝his(ys)，计算对于去噪后的血管三维图像的单个像素点s的模型G1、G2以及G3的概率表达P(ys|G1)、P(ys|G2)以及P(ys|G3)；
步骤4.2：针对去噪后的血管三维图像的单个像素，计算血管类的条件概率函数P(Y|V)和背景类的条件概率函数P(Y|B)；
步骤4.3：判断P(Y|V)是否大于P(Y|B)，若是，则认为该像素点属于血管类，否则认为该像素点属于背景类，得到该像素标签；
步骤4.4：重复上述步骤4.1-4.3，遍历去噪后血管三维图像的每一个像素，为每一个像素生成其对应的标签；
步骤4.5：将步骤4.4中针对每一个像素得到的相应的标签进行汇总，进而得到去噪的血管三维图像分割的初步结果Mask0；
步骤5：使用步骤1中去噪后血管三维图像进行基于脊线性质的血管中心线提取，具体为：提取去噪后血管三维图像中的种子点，并且对提取的种子点生成初始方向，之后针对生成的初始方向对每一个种子点沿正向与负向进行延伸以实现中心线生长，具体为：
步骤5.1：提取去噪的血管三维图像中的种子点,具体步骤如下：
步骤5.1.1：遍历经过步骤1中去噪后血管三维图像中的所有像素点，对于每一个像素点s将该点强度值Is与其邻域内的点N(s)强度值进行比较，判断像素点s的灰度是否大于邻域内任意一点灰度，如果是则认为该像素点是一个极值点，否则不认为该像素点是一个极值点；
步骤5.1.2：遍历步骤5.1.1得到的所有极值点，针对每个强度为I0，梯度方向为
的极值点s0(x0，y0，z0)进行一阶微分方向修正，找到点 s1(x1，y1，z1)，其强度为I1；
步骤5.1.3：对于通过步骤5.1.2得到的一阶微分方向修正过的点s1(x1，y1，z1)，计算其黑塞矩阵较大的两个特征值对应的特征向量
步骤5.1.4：根据步骤5.1.3针对修正过的点s1(x1，y1，z1)及
构成平面Ω，并以s1为圆心r为半径的范围内在平面Ω寻找灰度极值点作为种子点；
步骤5.2：针对提取的种子点生成初始方向，具体步骤如下：
步骤5.2.1：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点 s0沿着
寻找其邻域内强度最大点s1，判断s1强度是否小于s0，若是则结束此次初始方向的生成，否则将其作为正向传播的下一点s1，称为后沿点；
其中，将从种子点s0指向后沿点s1的方向作为正向传播的方向
步骤5.2,2：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s0沿着
寻找其邻域内强度最大点s-1，判断s-1强度小于s0，若是则结束此次初始方向的生成，否则将其作为负向传播的下一点s-1，称为前驱点；
其中，将从种子点s0指向前驱点s-1的方向作为负向传播的方向
步骤5.2.3：判断步骤5.2.1和5.2.2处理得到的正向传播与负向传播夹角是否小于阈值θ，若是则将s-1从s0的邻域点集内移除，跳至步骤5.2.1；否则完成该种子点的初始方向的生成，跳至步骤5.3；
步骤5.3：针对步骤5.2生成的初始方向对每一个种子点沿正向与负向进行延伸实现中心线生长；
步骤5.4：对于步骤5.1提取的每一个种子点重复步骤5.2—5.3，得到该血管三维图像的血管中心线图像；
步骤5.5：对于步骤5.4中提取到的血管中心线做膨胀处理，得到处理之后的中心线图像；
步骤6：引入血管特征场及血管形状场，计算势能函数，并且根据势能函数计算高层能量函数，具体步骤如下：
步骤6.1：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一个像素点s，计算该点血管特征场对于高层能量函数的贡献
其中，像素点s的邻域点集为N(s)，针对邻域内一点任意一点s′∈N(s)使用
表示该像素点处梯度，
表示其邻域内梯度最大值；
步骤6.2：重复步骤7.1，对去噪的血管三维图像数据的每一个像素点求解血管特征场；
步骤6.3：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一像素点s，计算该点血管方向场对于高层能量函数的贡献
其中，像素点s的邻域点集为N(s)，使用表
示s处梯度，其邻域内一点s′∈ N(s)，使用
表示s′处最小特征值对应特征向量，则s′处血管方向场对于高层能量函数的贡献为
步骤6.4：重复步骤6.3，对去噪的血管三维图像数据的每一个像素点求解血管方向场；
步骤6.5：基于步骤6.2和6.4求解的针对所有像素点的血管特征场与血管方向场，对于各个像素点，通过步骤4.5得到的初始标签Mask0计算出该像素点与其邻域内一点s′所对应的势能函数Ita(xs,xs′)，其为血管特征场与血管方向场的线性组合，用ε1,ε2控制各自的权重；
步骤6.6：基于步骤6.5的运算得到的势能函数Ita(xs,xs′)，根据结合s′处在Mask0中的标签计算对xs的能量函数的贡献值E(xs,xs')；
步骤6.7：对s的邻域内所有点进行计算加和，计算该点最终的高层能量函数U2(xs)；
步骤7：计算当前分割结果，直到图像分割结果Mask变化小于设立阈值，具体步骤如下:
步骤7.1：结合步骤3得到的低层能量函数U1(xs)与步骤6.7得到的高层能量函数U2(xs)以及步骤5.5处理之后的中心线图像可以得到点s处的标签概率 P(xs)
步骤7.2：判断此次去噪后血管三维图像分割的结果与上一次去噪后血管三维图像的分割结果间的相似性系数是否低于阈值c，若是，则输出分割的标签 Mask0作为最终的结果，否则Mask0设置为此次去噪的血管三维图像分割的结果并且跳转到对步骤5；
其中，阈值c的取值范围为0.01到0.05。
有益效果
本发明一种基于统计学模型和多尺度滤波的TOF-MRA脑血管图像分割方法，与已有技术相比较，具有如下有益效果：
1.所述血管图像分割方法针对统计学模型中应用脊线性质，在3D空间中提取血管中心线，沿远端血管中心线进行膨胀，对低层模型结果进行补充，提高低层模型对远端血管图像分割精度，解决了低层模型对于远端血管图像分割效果较差的问题；
2.所述方法基于脊线的3D中心线提取，即利用统计学模型在高强度像素点上分割提取的优势，结合脊线性质对种子点进行筛选，可以实现在较为粗糙的血管图像分割结果下得到精准的血管中心线提取。
附图说明
图1为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法的中心线膨胀示意图；
图2为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法的流程图；
图3为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法中地面真实标注图像的示意图；
图4为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法一张切片下原图像(a)与增强图像(b)的对比示意图；
图5中应用本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法给出了颅骨剥离后的血管图像灰度分布直方图Iinit；
图6中应用本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法低层模型结果对比；
图7中应用本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法种子点延伸示意图对比；
图8为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法中心线分割的两例结果对比图；
图9为本发明一种基于统计学模型和多尺度滤波的脑血管图像分割方法针对Brave数据集的运行结果。
具体实施方式
下面结合附图和实例，对本发明提供的一种基于统计学模型和多尺度滤波的脑血管图像分割方法作详细地说明。
实施例1
脑血管疾病是导致人类死亡的主要疾病之一，对其早期预防及治疗极为重要。通过影像技术进行的无创伤脑血管病变的检测，是一种有效的诊断和监测脑血管疾病的方法，也是最能被人所接受而不会产生副作用的技术。此项研究为脑血管疾病的诊断与预防提供技术支持，相关技术可以推广到多模式图像配准和三维可视化等其它领域。由于我国每年参加体检和预防性保健人数规模庞大，此研究成果的应用可产生广泛的社会效益和经济效益。
本实施例中以Brave数据集作为我们的目标对象；Brave，即Brain and VascularHealth in the Elderly，此数据集研究的目的是使用多模式核磁功能成像来研究血管健康和大脑老化之间的相互关系。这项研究连续招募了200名老年男性，他们患有原发性高血压，没有临床上明显的神经、心脑血管疾病。所有与会者都提供了书面知情同意。
数据集在进行采集时，由于设备原因使用了分段扫描并拼接，使得整张图象在成像时呈现了明显的灰度阶梯，这部分噪声使得最终成像效果较差，进而影响到了分割结果，使得整体参数较低。由于医疗数据获取地面真实困难，对于一例 TOF-MRA的血管进行标注就需要一位有经验的外科医生3-4天的工作量，本课题所用数据集对于血管分割部分仅标注了关心的血管结果，对于其余血管并未进行标注。具体标注结果如附图3所示，在附图3中左侧图为TOF-MRA切片，右侧图为对应切片的地面真实与切片叠加展示，使用红色表示地面真实的标注，本实施例在Brave数据集中随机选取了20例数据进行试验；结合实际情况，对于血管分割结果本课题使用了真实阳性值作为核心指标。
地面真实，也即Ground Truth，即为真值、真实的有效值或者是标准的答案。
本实施例中的一种基于统计学模型和多尺度滤波的TOF-MRA脑血管图像分割方法，其操作流程如附图2所示，具体实现步骤为：
步骤1：数据预处理，得到去噪后脑部区域图像，具体为：采集大脑部分血管原始图像、并对得到的血管三维图像进行去噪，得到去噪后血管三维图像，并根据去噪后血管三维图像构建血管三维图像的灰度的分布直方图；具体为：
步骤1.1：通过TOF-MRA及SNAP-MRA采集出大脑部分血管的原始图像；
其中，TOF-MRA，即飞行时间磁共振血管成像；SNAP-MRA即非对比增强血管与斑块内出血同时成像；血管三维图像即颅骨剥离后的脑部区域图像；
步骤1.2：使用FSL软件的bet命令对步骤1.1采集的原始图像进行颅骨剥离，将原始数据去除眼球脂肪等不相干的图像，得到血管三维图像；
其中，步骤1.2中对原始图像进行处理的原因是：TOF-MRA及SNAP-MRA采集时由于设备限制，无法单独采集脑部区域；
步骤1.3：选用大小为3×3×3，方差为0.4的高斯核对血管三维图像颅骨剥离后的脑部区域图像进行高斯滤波,滤除图像中的无关噪声，得到去噪后脑部区域图像，具体包括如下子步骤：
步骤1.3.1：构建血管三维图像的黑塞矩阵，其中，血管三维图像中每一个像素点的一阶导数表示灰度的梯度，二阶导数是灰度梯度的梯度，即灰度的曲率，其黑塞矩阵的表达形式为公式(1)：
步骤1.3.2：求出血管三维图像中每一个像素点的特征向量后，将该像素点的曲率进行分解；判断出当前像素点的曲率主方向及该方向对应的照其绝对值大小顺序排列的特征向量λ1，λ2和λ3；
步骤1.3.3：建立针对某一像素点(x,y,z)的高斯卷积核G(x,y,z)，如公式 (2)所示：
步骤1.3.4将高斯卷积核G(x,y,z)与对应像素点(x,y,z)的体素强度 I(x,y,z)进行卷积得到，在尺度σ下的体素强度I′(x，y，z，σ)；
步骤1.3.5重复步骤1.3.3到1.3.4，得到在尺度σ下血管三维图像中所有像素点的体素强度，计算方法如公式(3)所示：
I′(x,y,z,σ)＝G(x,y,z,σ)*I(x,y,z) (3)
步骤1.3.6：使用Hessian滤波器检测出血管三维图像中的管状结构，再对管状结构进行增强，对其余结构信号进行抑制减弱，以增强血管组织的对比度，得到增强后的血管三维图像，具体为：如果特征值λ2和λ3中都小于0，那么我们引入超参数α用于区分板状结构与管状结构；β用于区分斑点状结构与管状结构；c是用于区分板状结构与斑点状结构，计算对于某一像素点s在尺度σ下增强之后的像素表征V0(s，σ)，如公式(4)所示：
步骤1.3.7：重复步骤1.3.6，在每一个尺度下σ，进行计算该尺度时各个像素点计算得到的结果；
步骤1.3.8：将步骤1.3.7中各个尺度下得到的所有的V0(s，σ)中的最大值作为去噪后血管三维图像该像素点的像素值V(s)，如公式(5)所示：
其中，去噪后血管三维图像即去噪后脑部区域图像；
附图4给出了图像中一张切片下原图像(a)与增强图像(b)的对比。图像经过增强后，可以鲜明的看到其中的血管部分被增强，但其中也有少量的非血管组织也被错误的增强；同时由于血管尺度不同，其增强的血管会有部分变得比原图像更加细，而由于这样的原因多尺度滤波器更多的是作为血管图像分割时图像的预处理手段而不是直接作为分割结果；
步骤1.3.9：根据步骤1.3.8得到的处理之后的血管三维图像的像素图V(s) 的灰度分布情况构建灰度分布直方图Iinit；
附图5中给出了颅骨剥离后的血管图像灰度分布直方图Iinit：纵坐标表示频率，横坐标表示图像像素的灰度；在灰度分布直方图中，最低强度级别主要对应于脑脊液、骨骼和背景空气；中等强度级别对应于脑组织，包括灰质和白质以及眼睛的部分；高强度级别由邻近皮肤的脂肪和动脉组成；基于TOF-MRA的成像原理，血液流速越快的部分强度越高，即大血管理论上是图像中最亮的部分，与周围组织对比度较高，也更加容易通过灰度的差异来进行分离；在血管末梢部位血管变细，血液流速降低，对应部位在图像中的强度、与周围组织的对比度也会下降，较难通过灰度差异分离；在低层图像处理时，使用的是颅骨剥离后的图像，这是由于骨骼部分对应的强度范围内几乎不会有血管体素的存在，然而低强度区域的灰度分布在拟合时会引入新的误差，进而影响中强度区域及高强度区域的拟合效果；同时去掉颅骨区域后，减少了灰度拟合时的模型个数，更容易拟合；
步骤2：通过分别拟合去噪后脑部区域图像背景及血管类的高斯模型，搭建有限混合模型FMM，具体为：
步骤2.1：使用两个高斯模型G1(x)、G2(x)拟合步骤1.3.8得到的去噪后血管三维图像中的背景类的背景类，其中α1、μ1、σ12代表G1(x)中的权重、均值与方差，α2、μ2、σ22代表G2(x)中的权重、均值与方差。对于一个高斯模型Gn(x)，其某像素点x(i)处的公式表达为公式(6)：
步骤2.2：单独使用一个高斯模型G3(x)拟合其血管类，其参数表示为α3、μ3、σ32，其高斯函数的表达形式如公式(6)所示。
步骤2.3：根据步骤2.1与步骤2.2拟合出来的三个高斯模型G1(x)、G2(x) 以及G3(x)搭建有限混合模型FMM以拟合灰度分布his()，公式如(7)所示：
his(x)＝α1×G1(x)+α2×G2(x)+α3×G3(x) (7)
其中超参数α1，α2，α3用于调节步骤2.1与步骤2.2拟合出来的三个高斯模型之间的权重；
步骤3：使用最大期望方法对有限混合模型FMM设置初始值，对初始值进行参数估计，并计算低层能量函数，具体为：为有限混合模型设置初始值，并对有限混合模型执行最大期望方法，估计该有限混合模型FMM设置初始值，具体为：
步骤3.1：为有限混合模型设置初始值，得到三个高斯模型的参数集θ，具体包括如下子步骤：
步骤3.1.1：首先初始化步骤2.2拟合得到的血管类的高斯分布G3(x)，将其在步骤2.3中得到的his(x)表达式中的占比α3设定为区间0.03；
其中，步骤3.1.1中将超参数α3设定为0.03的原因是：人的大脑中血管体素所占比例通常为0.01～0.05，将其设立为0.03可以较好的进行拟合；
步骤3.1.2：将经过步骤2.2拟合出来的高斯分布函数G3(x)的初始均值μ3及方差σ32使用经过步骤1.3.8的中经去噪的血管三维图像中的最亮的前3％点进行计算得出；
步骤3.1.3：将经过步骤2.1拟合出来的G2(x)的初始均值μ2设定为步骤 1.3.9中得到的灰度直方图Iinit的最高点的灰度；
步骤3.1.4：将经过步骤2.1拟合出来的G1(x)的初始均值μ1设立为G2(x)的初始均值μ2的0.3倍；
步骤3.1.5：计算G1(x)贡献在步骤1.3.9中的得到的灰度直方图中Iinit占比作为α1的估计值，引入超参数k制衡G1(x)的峰值高度，计算方法参照公式(8)：
步骤3.1.6：根据步骤3.1.5计算得到的结果，计算中去掉G1(x)贡献后的剩余分布直方图Ileft，计算方法参照公式(9)：
步骤3.1.7：取步骤3.1.6得到的剩余分布直方图Ileft在区间[(μ2-μ1)/ 2，(μ2+μ1)/2]上的灰度点进行最大似然估计计算方差σ22；
步骤3.1.8：根据步骤3.1.1的求解的α3和步骤3.1.5的求解的α1，G2(x) 的占比α2使用1-α1-α3得到；
均值参数的集合θμ、方差参数的的集合θσ及体素占比参数的集合θα，统称为参数集θ；
步骤3.2：对有限混合模型执行最大期望方法，具体步骤如下：
步骤3.2.1：设定循环参数t＝1；
步骤3.2.2：根据三个高斯模型的参数集θ，计算隐变量γ其联合概率的条件概率P，计算方法参照公式(10)；
步骤3.2.3：根据步骤3.2.2计算出来的隐变量γ其联合概率的条件概率P，计算期望E，如公式(11)所示：
步骤3.2.4：根据步骤3.2.3计算得到的期望E，计算使得其最大的参数取值，作为t+1次循环的参数集，将体素占比的和，即∑θα设为1，得到新拉格朗日函数L，如公式(12)所示：
步骤3.2.4：使用步骤3.2.3得到的L对λ、αz、μz、σz求偏导，令其偏导数为0，可得新的参数集θ′，如公式(13)所示；
步骤3.2.5：使用步骤3.2.4得到的新拉格朗日函数对λ、αz、μz、σz求偏导，令其偏导数为0，得到新的参数集；
步骤3.2.6：判断t是否大于等于2，若是，则判断第t次及第t-1两次循环的所有参数的变化是否均小于1％，若是，将步骤3.2.5得到的新的参数集作为 his(ys)函数的初始参数的估计，并跳至步骤4；否则判断t是否等于50，若是，将步骤3.2.5得到的新的参数集作为his(ys)函数的初始参数的估计，并跳至步骤4，否则将t加1，跳至步骤3.2.2；
至此，完成了对his(ys)函数的初始参数θ的估计；
步骤4：使用步骤2搭建的有限混合模型FMM与步骤1中去噪后血管三维图像的像素值，计算得到去噪的血管三维图像分割的初步结果，具体步骤如下：
步骤4.1：通过步骤3得到的his(ys)函数的初始参数θ的估计，提出低层能量函数U1(xs)满足U1(xs)∝his(ys)，计算对于去噪后的血管三维图像的单个像素点s的模型G1、G2以及G3的概率表达P(ys|G1)、P(ys|G2)以及P(ys|G3)。对于单个像素点s，计算其属于第i个模型Gi的概率表达P(ys|Gi)，如公式(14)所示：
步骤4.2：针对去噪后的血管三维图像的单个像素，计算血管类的条件概率函数P(Y|V)＝P(ys|G3)和背景类的条件概率函数P(Y|B)＝P(ys|G1)+P(ys|G2)；
步骤4.3：判断P(Y|V)是否大于P(Y|B)，若是，则认为该像素点属于血管类，否则认为该像素点属于背景类，得到该像素标签；
步骤4.4：重复上述步骤4.1-4.3，遍历去噪后血管三维图像的每一个像素，为每一个像素生成其对应的标签；
步骤4.5：将步骤4.4中针对每一个像素得到的相应的标签进行汇总，进而得到去噪的血管三维图像分割的初步结果Mask0；
附图6给出了两例病人的通过底层模型的血管分割的对比图，每一行代表一例图像，左侧图是三高斯模型的生成结果、中间图为分段三高斯模型结果，右侧图为标准结果；图像结果均使用最大密度投影图呈现；
步骤5：使用步骤1中去噪后血管三维图像进行基于脊线性质的血管中心线提取，具体为：提取去噪后血管三维图像中的种子点，并且对提取的种子点生成初始方向，之后针对生成的初始方向对每一个种子点沿正向与负向进行延伸以实现中心线生长，具体为：
步骤5.1：提取去噪的血管三维图像中的种子点,具体步骤如下：
步骤5.1.1：遍历经过步骤1中去噪后血管三维图像中的所有像素点，对于一个像素点s，将该点强度值Is与其邻域内的点N(s)强度值进行比较，判断像素点 s的灰度是否大于邻域内任意一点灰度，如果是则认为该像素点是一个极值点，否则不认为该像素点是一个极值点；
步骤5.1.2：遍历步骤5.1.1得到的所有极值点，针对每个强度为I0，梯度方向为
的极值点s0(x0，y0，z0)进行一阶微分方向修正，找到点 s1(x1，y1，z1)，其强度为I1，如公式(15)所示；
步骤5.1.3：对于通过步骤5.1.2得到的一阶微分方向修正过的点s1(x1，y1，z1)，计算其黑塞矩阵较大的两个特征值对应的特征向量
步骤5.1.4：根据步骤5.1.3针对修正过的点s1(x1，y1，z1)及
构成平面Ω，并以s1为圆心r为半径的范围内在平面Ω寻找灰度极值点作为种子点，计算方法参考公式(16)；
步骤5.2：针对提取的种子点生成初始方向，具体步骤如下：
步骤5.2.1：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点 s0沿着
寻找其邻域内强度最大点s1，若s1强度小于s0，则结束此次初始方向的生成，否则将其作为正向传播的下一点s1，称为后沿点；
其中，将从种子点s0指向后沿点s1的方向作为正向传播的方向
步骤5.2.2：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s0沿着
寻找其邻域内强度最大点s-1，若s-1强度小于s0，则结束此次初始方向的生成，否则将其作为负向传播的下一点s-1，称为前驱点；
其中，将从种子点s0指向前驱点s-1的方向作为负向传播的方向
步骤5.2.3：判断步骤5.2.1和5.2.2处理得到的正向传播与负向传播夹角是否小于阈值θ，若是则将s-1从s0的邻域点集内移除，跳至步骤5.2.1；否则完成该种子点的初始方向的生成，跳至步骤5.3；
种子点初始方向生成过程如附图7所示，其中P0点表示当前点，对当前点有正向传播及负向传播，直到延伸至强度达到生长阈值或是强度开始下降，则停止生长。
步骤5.3：针对步骤5.2生成的初始方向对每一个种子点沿正向与负向进行延伸实现中心线生长；
其中，步骤5.2-5.3应用了脊线性质，在3D空间中提取血管中心线，沿远端血管中心线进行膨胀，对低层模型结果进行补充，提高低层模型对远端血管图像分割精度，解决了低层模型对于远端血管图像分割效果较差的问题，以达到我们预期的有益效果；
步骤5.4：对于步骤5.1提取的每一个种子点重复步骤5.2—5.3，得到该血管三维图像的血管中心线图像；
附图8给出了两例结果对比图，其中第一列为本专利所述中心线提取方法，第三列为中心线的标准分割结果，每一行代表了一例不同的图像。
步骤5.5：对于步骤5.4中提取到的血管中心线做膨胀处理，得到处理之后的中心线图像，对于步骤5.4中提取到的血管中心线与血管图像分割的结果相交点，如附图1中点P，若其在实际图像中对应点P0处半径为R，则对中心线PQ 做膨胀处理，使得血管半径沿PQ逐渐变小，直到汇聚于PQ延长线。若PQ长度为D，对于PQ中与P距离为d的一点，其则其膨胀核尺寸r满足
因中心线点不会出现在血管末梢，因此对于r加上限制r＝max(r,2)；
步骤6：引入血管特征场及血管形状场，计算势能函数，并且根据势能函数计算高层能量函数，具体步骤如下：
步骤6.1：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一个像素点s，计算该点血管特征场对于高层能量函数的贡献
其中，像素点s的邻域点集为N(s)，针对邻域内一点任意一点s′∈N(s)使用
表示该像素点处梯度，
表示其邻域内梯度最大值；
步骤6.2：重复步骤7.1，对去噪的血管三维图像数据的每一个像素点求解血管特征场；
步骤6.3：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一像素点s，计算该点血管方向场对于高层能量函数的贡献
其中，像素点s的邻域点集为N(s)，使用表
示s处梯度，其邻域内一点s′∈ N(s)，使用
表示s′处最小特征值对应特征向量，则s′处血管方向场对于高层能量函数的贡献为
步骤6.4：重复步骤6.3，对去噪的血管三维图像数据的每一个像素点求解血管方向场；
步骤6.5：基于步骤6.2和6.4求解的针对所有像素点的血管特征场与血管方向场，对于各个像素点，通过步骤4.5得到的初始标签Mask0计算出该像素点与其邻域内一点s′所对应的势能函数Ita(xs,xs′)，其为血管特征场与血管方向场的线性组合，用ε1,ε2控制各自的权重，如公式(17)所示：
步骤6.6：基于步骤6.5的运算得到的势能函数Ita(xs,xs′)，根据结合s′处在Mask0中的标签计算对xs的能量函数的贡献值E(xs,xs')，如公式(18)所示：
步骤6.7：对s的邻域内所有点进行计算加和，计算该点最终的高层能量函数U2(xs)，通过设置常数βc其随着团结构的不同及图像质量不同而变化。对图像内所有点进行计算，分别记录其对应Lv以及LB的能量函数，并将其用于最终血管图像分割结果的计算，如公式(19)所示；
步骤7：计算当前分割结果，直到图像分割结果Mask变化小于设立阈值，具体步骤如下:
步骤7.1：结合步骤3得到的低层能量函数U1(xs)与步骤6.7得到的高层能量函数U2(xs)以及步骤5.5处理之后的中心线图像可以得到点s处的标签概率 P(xs)，如公式(20)所示:
步骤7.2：判断此次去噪后血管三维图像分割的结果与上一次去噪后血管三维图像的分割结果间的相似性系数是否低于阈值0.03，若是，则输出分割的标签Mask0作为最终的结果，否则Mask0设置为此次去噪的血管三维图像分割的结果并且跳转到对步骤5；
针对Brave数据集，上述方法的运行结果如附图9所示：
图9中的Brave数据集实验结果展示了对比第一列为本课题所提方法结果，第二列血管提取的标准结果。分割结果图像均使用最大密度投影图给出。针对 Brave数据集对比结果，可以发现针对图像质量较低、信噪比较差图像，本专利所提方法对血管图像分割准确性有实际提升，尤其对于统计学模型表现不佳的细小血管部分，由于中心线方法的引入，使得该部分结果相较于传统统计学模型有较大的改善，细小血管的分割结果变得更粗，证明了其有效性。
上述描述对本发明的特征和方法进行了具体的说明，但应了解，在所述权利要求中定义的本发明并不局限于所述的具体特征或方法。本领域人员可在权利要求的范围内做出修改，并不影响本发明的实质内容。
Claims (10)
1.一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：包括如下步骤：
步骤1：数据预处理，得到去噪后脑部区域图像，具体为：采集大脑部分血管原始图像、并对得到的血管三维图像进行去噪，得到去噪后血管三维图像，并根据去噪后血管三维图像构建血管三维图像的灰度的分布直方图；
步骤2：通过分别拟合去噪后脑部区域图像背景及血管类的高斯模型，搭建有限混合模型FMM，具体为：
步骤2.1：使用两个高斯模型G1(x)、G2(x)拟合步骤1的去噪后血管三维图像中的背景类；
步骤2.2：使用一个高斯模型G3(x)拟合步骤1.3.8得到的去噪后血管三维图像中的血管类；
步骤2.3：根据步骤2.1与步骤2.2拟合出来的三个高斯模型G1(x)、G2(x)以及G3(x)搭建有限混合模型FMM，其FMM公式描述的灰度分布his()是步骤2.1与步骤2.2拟合出来的三个高斯模型的线性组合，his()利用超参数α1，α2，α3调节三个高斯模型G1(x)、G2(x)以及G3(x)间各自的权重；
步骤3：使用最大期望方法对有限混合模型FMM设置初始值，对初始值进行参数估计，并计算低层能量函数，具体为：为有限混合模型设置初始值，并对有限混合模型执行最大期望方法，估计该有限混合模型FMM设置初始值；
步骤4：使用步骤2搭建的有限混合模型FMM与步骤1中去噪后血管三维图像的像素值，计算得到去噪的血管三维图像分割的初步结果，具体步骤如下：
步骤4.1：通过步骤3得到的his(ys)函数的初始参数θ的估计，提出低层能量函数U1(xs)满足U1(xs)∝his(ys)，计算对于去噪后的血管三维图像的单个像素点s的模型G1、G2以及G3的概率表达P(ys|G1)、P(ys|G2)以及P(ys|G3)；
步骤4.2：针对去噪后的血管三维图像的单个像素，计算血管类的条件概率函数P(Y|V)和背景类的条件概率函数P(Y|B)；
步骤4.3：判断P(Y|V)是否大于P(Y|B)，若是，则认为该像素点属于血管类，否则认为该像素点属于背景类，得到该像素标签；
步骤4.4：重复上述步骤4.1-4.3，遍历去噪后血管三维图像的每一个像素，为每一个像素生成其对应的标签；
步骤4.5：将步骤4.4中针对每一个像素得到的相应的标签进行汇总，进而得到去噪的血管三维图像分割的初步结果Mask0；
步骤5：使用步骤1中去噪后血管三维图像进行基于脊线性质的血管中心线提取，具体为：提取去噪后血管三维图像中的种子点，并且对提取的种子点生成初始方向，之后针对生成的初始方向对每一个种子点沿正向与负向进行延伸以实现中心线生长，具体为：
步骤5.1：提取去噪的血管三维图像中的种子点，具体步骤如下：
步骤5.1.1：遍历经过步骤1中去噪后血管三维图像中的所有像素点，对于每一个像素点s将该点强度值Is与其邻域内的点N(s)强度值进行比较，判断像素点s的灰度是否大于邻域内任意一点灰度，如果是则认为该像素点是一个极值点，否则不认为该像素点是一个极值点；
步骤5.1.2：遍历步骤5.1.1得到的所有极值点，针对每个强度为I0，梯度方向为
的极值点s0(x0，y0，z0)进行一阶微分方向修正，找到点s1(x1，y1，z1)，其强度为I1；
步骤5.1.3：对于通过步骤5.1.2得到的一阶微分方向修正过的点s1(x1，y1，z1)，计算其黑塞矩阵较大的两个特征值对应的特征向量
步骤5.1.4：根据步骤5.1.3针对修正过的点s1(x1，y1，z1)及
陶成平面Ω，并以s1为圆心r为半径的范围内在平面Ω寻找灰度极值点作为种子点；
步骤5.2：针对提取的种子点生成初始方向，具体步骤如下：
步骤5.2.1：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s0沿着
寻找其邻域内强度最大点s1，判断s1强度是否小于s0，若是则结束此次初始方向的生成，否则将其作为正向传播的下一点s1，称为后沿点；
其中，将从种子点s0指向后沿点s1的方向作为正向传播的方向
步骤5.2，2：对于通过步骤5.1.3提取的所有的种子点，对于每一个种子点s0沿着
寻找其邻域内强度最大点s-1，判断s-1强度小于s0，若是则结束此次初始方向的生成，否则将其作为负向传播的下一点s-1，称为前驱点；
其中，将从种子点s0指向前驱点s-1的方向作为负向传播的方向
步骤5.2.3：判断步骤5.2.1和5.2.2处理得到的正向传播与负向传播夹角是否小于阈值θ，若是则将s-1从s0的邻域点集内移除，跳至步骤5.2.1；否则完成该种子点的初始方向的生成，跳至步骤5.3；
步骤5.3：针对步骤5.2生成的初始方向对每一个种子点沿正向与负向进行延伸实现中心线生长；
步骤5.4：对于步骤5.1提取的每一个种子点重复步骤5.2-5.3，得到该血管三维图像的血管中心线图像；
步骤5.5：对于步骤5.4中提取到的血管中心线做膨胀处理，得到处理之后的中心线图像；
步骤6：引入血管特征场及血管形状场，计算势能函数，并且根据势能函数计算高层能量函数，具体步骤如下：
步骤6.1：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一个像素点s，计算该点血管特征场对于高层能量函数的贡献
其中，像素点s的邻域点集为N(s)，针对邻域内一点任意一点s′∈N(s)使用
表示该像素点处梯度，
表示其邻域内梯度最大值；
步骤6.2：重复步骤6.1，对去噪的血管三维图像数据的每一个像素点求解血管特征场；
步骤6.3：针对步骤1.3.8得到的经过预处理的去噪的血管三维图像数据的每一像素点s，计算该点血管方向场对于高层能量函数的贡献
其中，像素点s的邻域点集为N(s)，使用表
示s处梯度，其邻域内一点s′∈N(s)，使用
表示s′处最小特征值对应特征向量，则s′处血管方向场对于高层能量函数的贡献为
步骤6.4：重复步骤6.3，对去噪的血管三维图像数据的每一个像素点求解血管方向场；
步骤6.5：基于步骤6.2和6.4求解的针对所有像素点的血管特征场与血管方向场，对于各个像素点，通过步骤4.5得到的初始标签Mask0计算出该像素点与其邻域内一点s′所对应的势能函数Ita(xs，xs′)，其为血管特征场与血管方向场的线性组合，用ε1，ε2控制各自的权重；
步骤6.6：基于步骤6.5的运算得到的势能函数Ita(xs，xs′)，根据结合s′处在Mask0中的标签计算对xs的能量函数的贡献值E(xs，xs′)；
步骤6.7：对s的邻域内所有点进行计算加和，计算该点最终的高层能量函数U2(xs)；
步骤7：计算当前分割结果，直到图像分割结果Mask变化小于设立阈值，具体步骤如下：
步骤7.1：结合步骤3得到的低层能量函数U1(xs)与步骤6.7得到的高层能量函数U2(xs)以及步骤5.5处理之后的中心线图像可以得到点s处的标签概率P(xs)
步骤7.2：判断此次去噪后血管三维图像分割的结果与上一次去噪后血管三维图像的分割结果间的相似性系数是否低于阈值c，若是，则输出分割的标签Mask0作为最终的结果，否则Mask0设置为此次去噪的血管三维图像分割的结果并且跳转到对步骤5；
其中，阈值c的取值范围为0.01到0.05。
2.依据权利要求1所述的一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：步骤1，具体为：
步骤1.1：通过TOF-MRA及SNAP-MRA采集出大脑部分血管的原始图像；
步骤1.2：使用FSL软件的bet命令对步骤1.1采集的原始图像进行颅骨剥离，将原始数据去除眼球脂肪等不相干的图像，得到血管三维图像；
步骤1.3：选用大小为n*n*n，方差为se的高斯核对血管三维图像颅骨剥离后的脑部区域图像进行高斯滤波，滤除图像中的无关噪声，得到去噪后脑部区域图像，具体包括如下子步骤：
步骤1.3.1：构建血管三维图像的黑塞矩阵，其中，血管三维图像中每一个像素点的一阶导数表示灰度的梯度，二阶导数是灰度梯度的梯度，即灰度的曲率；
步骤1.3.2：求出血管三维图像中每一个像素点的特征向量后，将该像素点的曲率进行分解；判断出当前像素点的曲率主方向及该方向对应的照其绝对值大小顺序排列的特征向量λ1，λ2和λ3；
步骤1.3.3：建立针对某一像素点(x，y，z)的高斯卷积核G(x，y，z)：
步骤1.3.4：将高斯卷积核G(x，y，z)与对应像素点(x，y，z)的体素强度I(x，y，z)进行卷积得到，在尺度σ下的体素强度I′(x，y，z，σ)；
步骤1.3.5重复步骤1.3.3到1.3.4，得到在尺度σ下血管三维图像中所有像素点的体素强度；
步骤1.3.6：使用Hessian滤波器检测出血管三维图像中的管状结构，再对管状结构进行增强，对其余结构信号进行抑制减弱，以增强血管组织的对比度，得到增强后的血管三维图像，具体为：如果特征值λ2和λ3中都小于0，那么我们引入超参数α用于区分板状结构与管状结构；β用于区分斑点状结构与管状结构；c是用于区分板状结构与斑点状结构，计算对于某一像素点s在尺度σ下增强之后的像素表征V0(s，σ)；
步骤1.3.7：重复步骤1.3.6，在每一个尺度下σ，进行计算该尺度时各个像素点计算得到的像素表征V0(s，σ)；
步骤1.3.8：将步骤1.3.7中各个尺度下得到的所有的V0(s，σ)中的最大值作为去噪后血管三维图像该像素点的像素值V(s)；
其中，去噪后血管三维图像即去噪后脑部区域图像；
步骤1.3.9：根据步骤1.3.8得到的去噪后的血管三维图像的像素值V(s)的灰度分布情况构建灰度分布直方图Iinit。
3.依据权利要求2所述的一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：步骤1.1中，TOF-MRA，即飞行时间磁共振血管成像；SNAP-MRA即非对比增强血管与斑块内出血同时成像；血管三维图像即颅骨剥离后的脑部区域图像；步骤1.2中对原始图像进行处理的原因是：TOF-MRA及SNAP-MRA采集时由于设备限制，无法单独采集脑部区域。
4.依据权利要求2所述的一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：步骤1.3中，n的取值范围是2到10，se的取值范围是0.2-1.0。
5.依据权利要求2所述的一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：步骤2中，去噪后脑部区域图像即去噪后血管三维图像。
6.依据权利要求1所述的一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：步骤3，具体为：
步骤3.1：为有限混合模型设置初始值，得到三个高斯模型的参数集θ，具体包括如下子步骤：
步骤3.1.1：初始化步骤2.2拟合得到的血管类的高斯分布G3(x)，将其在步骤2.3中得到的his(x)表达式中的占比α3设定为区间0.01～0.05中的常数值；
步骤3.1.2：使用经过步骤1.3.8中去噪的血管三维图像中最亮的前M％的像素点计算得出经过步骤2.2拟合出来的高斯分布函数G3(x)的初始均值μ3及方差σ32；
步骤3.1.3：将经过步骤2.1拟合出来的G2(x)的初始均值μ2设定为步骤1.3.9中得到的灰度直方图Iinit的最高点的灰度；
步骤3.1.4：将经过步骤2.1拟合出来的G1(x)的初始均值μ1设立为G2(x)的初始均值μ2的t倍；
步骤3.1.5：计算G1(x)贡献在步骤1.3.9中的得到的灰度直方图中Iinit占比作为σ1的估计值；
步骤3.1.6：根据步骤3.1.5计算得到的G1(x)的体素占比α1，计算中去掉G1(x)贡献后的剩余分布直方图Ileft；
步骤3.1.7：取步骤3.1.6得到的剩余分布直方图Ileft在区间[(μ2-μ1)/2，(μ2+μ1)/2]上的灰度点进行最大似然估计计算G2(x)函数的方差σ22；
步骤3.1.8：根据步骤3.1.1的求解的α3和步骤3.1.5的求解的α1，G2(x)的占比α2使用1-α1-α3得到；
均值参数的集合θμ、方差参数的的集合θσ及体素占比参数的集合θα，统称为参数集θ；
步骤3.2：对有限混合模型执行最大期望方法，具体步骤如下：
步骤3.2.1：设定循环参数t＝1；
步骤3.2.2：根据三个高斯模型的参数集θ，计算隐变量γ其联合概率的条件概率P；
步骤3.2.3：根据步骤3.2.2计算出来的隐变量γ其联合概率的条件概率P，计算期望E；
步骤3.2.4：根据步骤3.2.3计算得到的期望E，计算使得其最大的参数取值，作为t+1次循环的参数集，将体素占比的和，即∑θα设为1，得到新拉格朗日函数；
步骤3.2.5：使用步骤3.2.4得到的新拉格朗日函数对λ、αz、μz、σz求偏导，令其偏导数为0，得到新的参数集；
步骤3.2.6：判断t是否大于等于2，若是，则判断第t次及第t-1两次循环的所有参数的变化是否均小于q％，若是，将步骤3.2.5得到的新的参数集作为his(ys)函数的初始参数的估计，并跳至步骤4；否则判断t是否等于N，若是，将步骤3.2.5得到的新的参数集作为his(ys)函数的初始参数的估计，并跳至步骤4，否则将t加1，跳至步骤3.2.2；
至此，完成了对his(ys)函数的初始参数θ的估计。
7.依据权利要求6所述的一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：步骤3.1.2中，M为取值1-5之间的常数。
8.依据权利要求6所述的一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：步骤3.1.4中，t为常数，取值范围为0.1-0.5。
9.依据权利要求6所述的一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：步骤3.2.6中，q的取值范围在0.5-1.5之间。
10.依据权利要求6所述的一种基于统计学模型和多尺度滤波的脑血管图像分割方法，其特征在于：步骤3.2.6中，N的取值范围为25到75。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111210464.1A CN114119626B (zh) | 2021-10-18 | 2021-10-18 | 一种基于统计学模型和多尺度滤波的脑血管图像分割方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111210464.1A CN114119626B (zh) | 2021-10-18 | 2021-10-18 | 一种基于统计学模型和多尺度滤波的脑血管图像分割方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114119626A true CN114119626A (zh) | 2022-03-01 |
| CN114119626B CN114119626B (zh) | 2024-07-19 |
Family
ID=80376250
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111210464.1A Active CN114119626B (zh) | 2021-10-18 | 2021-10-18 | 一种基于统计学模型和多尺度滤波的脑血管图像分割方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114119626B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115375719A (zh) * | 2022-08-18 | 2022-11-22 | 华科精准（北京）医疗科技有限公司 | 一种血管图像分割方法及装置 |
| CN115953413A (zh) * | 2023-03-13 | 2023-04-11 | 同心智医科技(北京)有限公司 | 一种mra图像分割方法、装置及存储介质 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2008091583A2 (en) * | 2007-01-23 | 2008-07-31 | Dtherapeutics, Llc | Image-based extraction for vascular trees |
| CN106296675A (zh) * | 2016-08-04 | 2017-01-04 | 山东科技大学 | 一种高噪声灰度不均匀图像的分割方法 |
| CN109102511A (zh) * | 2018-07-06 | 2018-12-28 | 深圳先进技术研究院 | 一种脑血管分割方法、系统及电子设备 |
| CN110853051A (zh) * | 2019-10-24 | 2020-02-28 | 北京航空航天大学 | 基于多注意密集连接生成对抗网络的脑血管影像分割方法 |
| WO2020177126A1 (zh) * | 2019-03-07 | 2020-09-10 | 深圳先进技术研究院 | 图像处理方法、系统、计算设备及存储介质 |
| CN111815663A (zh) * | 2020-06-29 | 2020-10-23 | 浙江工贸职业技术学院 | 一种基于Hessian矩阵和灰度法的肝血管分割系统 |
| CN112102327A (zh) * | 2019-06-18 | 2020-12-18 | 中国科学院深圳先进技术研究院 | 一种图像处理方法、装置及计算机可读存储介质 |
-
2021
- 2021-10-18 CN CN202111210464.1A patent/CN114119626B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2008091583A2 (en) * | 2007-01-23 | 2008-07-31 | Dtherapeutics, Llc | Image-based extraction for vascular trees |
| CN106296675A (zh) * | 2016-08-04 | 2017-01-04 | 山东科技大学 | 一种高噪声灰度不均匀图像的分割方法 |
| CN109102511A (zh) * | 2018-07-06 | 2018-12-28 | 深圳先进技术研究院 | 一种脑血管分割方法、系统及电子设备 |
| WO2020177126A1 (zh) * | 2019-03-07 | 2020-09-10 | 深圳先进技术研究院 | 图像处理方法、系统、计算设备及存储介质 |
| CN112102327A (zh) * | 2019-06-18 | 2020-12-18 | 中国科学院深圳先进技术研究院 | 一种图像处理方法、装置及计算机可读存储介质 |
| CN110853051A (zh) * | 2019-10-24 | 2020-02-28 | 北京航空航天大学 | 基于多注意密集连接生成对抗网络的脑血管影像分割方法 |
| CN111815663A (zh) * | 2020-06-29 | 2020-10-23 | 浙江工贸职业技术学院 | 一种基于Hessian矩阵和灰度法的肝血管分割系统 |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115375719A (zh) * | 2022-08-18 | 2022-11-22 | 华科精准（北京）医疗科技有限公司 | 一种血管图像分割方法及装置 |
| CN115953413A (zh) * | 2023-03-13 | 2023-04-11 | 同心智医科技(北京)有限公司 | 一种mra图像分割方法、装置及存储介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114119626B (zh) | 2024-07-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gul et al. | Deep learning techniques for liver and liver tumor segmentation: A review | |
| Zhang et al. | Mapping population-based structural connectomes | |
| Gebru et al. | Detection of cerebrovascular changes using magnetic resonance angiography | |
| Liu et al. | Automatic whole heart segmentation using a two-stage u-net framework and an adaptive threshold window | |
| CN107403201A (zh) | 肿瘤放射治疗靶区和危及器官智能化、自动化勾画方法 | |
| Rueda et al. | Feature-based fuzzy connectedness segmentation of ultrasound images with an object completion step | |
| CN110415234A (zh) | 基于多参数磁共振成像的脑部肿瘤分割方法 | |
| CN114119626B (zh) | 一种基于统计学模型和多尺度滤波的脑血管图像分割方法 | |
| Jo et al. | Segmentation of the main vessel of the left anterior descending artery using selective feature mapping in coronary angiography | |
| Kamiya | Deep learning technique for musculoskeletal analysis | |
| Bermejo-Peláez et al. | Emphysema classification using a multi-view convolutional network | |
| Hasan | A hybrid approach of using particle swarm optimization and volumetric active contour without edge for segmenting brain tumors in MRI scan | |
| Kabir | Early stage brain tumor detection on MRI image using a hybrid technique | |
| Goswami et al. | An analysis of image segmentation methods for brain tumour detection on MRI images | |
| Bui et al. | An automatic random walk based method for 3D segmentation of the heart in cardiac computed tomography images | |
| Brahim et al. | A 3D network based shape prior for automatic myocardial disease segmentation in delayed-enhancement MRI | |
| CN114926487A (zh) | 多模态影像脑胶质瘤靶区分割方法、系统及设备 | |
| Subudhi et al. | Extraction of brain from MRI images by skull stripping using histogram partitioning with maximum entropy divergence | |
| CN111429404B (zh) | 一种用于心脑血管检测的成像系统和方法 | |
| Jalab et al. | Fractional Renyi entropy image enhancement for deep segmentation of kidney MRI | |
| Hasan et al. | Lung Segmentation from Chest X-Ray Images Using Deeplabv3plus-Based CNN Model | |
| Delmoral et al. | Segmentation of pathological liver tissue with dilated fully convolutional networks: A preliminary study | |
| WO2019210124A1 (en) | Atlas for segmentation of retina layers from oct images | |
| Hussain et al. | Segmentation of brain MRI with statistical and 2D wavelet features by using neural networks | |
| Dachraoui et al. | Computerized image segmentation of multiple sclerosis lesions using fuzzy level set model |
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 |