CN107274399B - 一种基于Hession矩阵和三维形状指数的肺结节分割方法 - Google Patents
一种基于Hession矩阵和三维形状指数的肺结节分割方法 Download PDFInfo
- Publication number
- CN107274399B CN107274399B CN201710464681.0A CN201710464681A CN107274399B CN 107274399 B CN107274399 B CN 107274399B CN 201710464681 A CN201710464681 A CN 201710464681A CN 107274399 B CN107274399 B CN 107274399B
- Authority
- CN
- China
- Prior art keywords
- dimensional
- nodule
- lung
- image
- shape index
- 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.)
- Active
Links
- 239000011159 matrix material Substances 0.000 title claims abstract description 64
- 230000011218 segmentation Effects 0.000 title claims abstract description 42
- 238000000034 method Methods 0.000 title claims abstract description 38
- 230000002685 pulmonary effect Effects 0.000 title claims abstract description 18
- 238000001514 detection method Methods 0.000 claims abstract description 52
- 210000004072 lung Anatomy 0.000 claims abstract description 51
- 206010056342 Pulmonary mass Diseases 0.000 claims abstract description 40
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 6
- 210000004204 blood vessel Anatomy 0.000 claims description 20
- 230000000694 effects Effects 0.000 claims description 19
- 238000001914 filtration Methods 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 11
- 210000004224 pleura Anatomy 0.000 claims description 11
- 238000013178 mathematical model Methods 0.000 claims description 7
- 210000001519 tissue Anatomy 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 5
- 238000009499 grossing Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- 210000000621 bronchi Anatomy 0.000 claims description 3
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 238000009795 derivation Methods 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 3
- 239000013598 vector Substances 0.000 claims description 3
- 230000008859 change Effects 0.000 claims description 2
- 238000003384 imaging method Methods 0.000 claims description 2
- 238000002591 computed tomography Methods 0.000 description 17
- 238000003745 diagnosis Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 239000005337 ground glass Substances 0.000 description 3
- 206010058467 Lung neoplasm malignant Diseases 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 201000005202 lung cancer Diseases 0.000 description 2
- 208000020816 lung neoplasm Diseases 0.000 description 2
- 230000003211 malignant effect Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- 230000035945 sensitivity Effects 0.000 description 2
- 206010028980 Neoplasm Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000002360 explosive Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000003902 lesion Effects 0.000 description 1
- 230000036210 malignancy Effects 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 208000024891 symptom Diseases 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000002792 vascular Effects 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/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G06T5/70—
-
- 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
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/50—Depth or shape recovery
- G06T7/55—Depth or shape recovery from multiple images
-
- 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/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
-
- 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/20—Special algorithmic details
- G06T2207/20024—Filtering details
-
- 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
- G06T2207/30064—Lung nodule
Abstract
本发明公开了一种基于Hession矩阵和三维形状指数的肺结节分割方法。该方法充分利用医学CT影像，首先根据序列CT影像的灰度值利用最佳阈值分割序列肺实质，并构建三维肺实质的体数据；然后计算肺实质三维体数据中个每个体素点的Hession矩阵特征值；接着根据三维结节模型的形状特征，并结合Hession矩阵特征值和二维形状指数对三维形状指数进行构建；最终构建的三维类球形滤波器‑3D shape index结节检测函数，对肺实质三维体数据进行结节检测，并利用检测得到结节区域作为区域增长的多种子点，进行结合置信度的区域增长算法将结节进行三维分割。本发明操作简单，可以实现对不同类型的疑似肺结节进行自动的检测和分割，具较强的稳定性和准确性。
Description
技术领域
本发明属于医学图像处理领域，具体涉及了基于CT序列影像中Hession矩阵和三维形状指数的多尺度三维空间特征的肺结节自动分割方法。
背景技术
肺癌是发病率和死亡率最高的恶性肿瘤之一。其主要原因是细小的病变特征难以被发现，临床漏诊率高。临床症状出现时，大多数患者往往已经到了病程的中晚期，治愈率非常低。肺癌早期的检出对提高治愈率起着至关重要的作用。计算机X射线断层扫描(CT)具有较高的组织分辨率，被广泛地应用于肺结节的筛查。随着病灶区域成像精度的不断提高，CT扫描厚度不断减小，影像数据呈现爆炸式增长，由于肺部组织结构复杂和结节形状大小各异，医生对多层二维CT影像筛查结节并排除血管等其他组织的干扰的工作难度很大，容易出现诊断疲劳，甚至可能造成漏诊和误诊。因此，如何能够在不降低精度要求的前提下，快速高效的检出肺结节并分割，对后续肺结节良恶性自动诊断具有重大的意义。但是肺结节检测系统的检测准确率和效率并不能满足临床的需求，其主要原因是在肺结节的检测过程中，由于血管的密度、CT值等都与结节相似，而且在病变区域内，部分结节和血管都会交叉，这些都会导致结节的检测过程会产生敏感性不高和较高的假阳性，影响结节的检测准确率。
发明内容
本发明的目的在于克服现有肺结节检测技术中存在的缺陷，提供一种基于Hession矩阵和三维形状指数的肺结节分割方法。
本发明所采用的技术方案为：
一种基于Hession矩阵和三维形状指数的肺结节分割方法，包括以下步骤：
步骤A；肺实质三维体数据的构建：首先利用各向异性滤波对CT影像进行预处理，去除噪声图像的影像，然后利用最佳阈值法对图像中的肺实质区域进行提取，按照序列依次提取单张CT影像的肺实质区域，最终构建肺实质三维体数据；
步骤B；在三维图像中将结节、血管、胸膜分别近似看作球体、圆柱体和曲面，通过分析典型的球、圆柱和曲面的基本特征，构造三个理想的数学模型，并计算肺实质三维体数据中每个体素点的Hession矩阵特征值，并分析每个模型的三维特征；
步骤C；分析三维模型的形状特征，利用肺结节的三维特征将其与血管等疑似区域进行区分，通过将二维形状指数扩展到三维体数据上，并结合Hession矩阵特征值对三维形状指数进行构建，最终构建能够较为准确的区分结节和血管的三维形状指数；
步骤D；肺结节灰度值呈Gaussian分布，不同尺度空间上的高斯滤波函数对于复杂的医学影像具有高效的增强效果，对不同尺度不同类型的疑似肺结节，选取不同的尺度参数σ进行多尺度增强，得到多尺度的三维类球形滤波器-3D shape index结节检测函数；
步骤E；利用步骤D中得到的3D shape index结节检测函数对肺实质三维体数据进行结节检测，并利用检测得到结节区域作为区域增长的多种子点，结合置信连接的区域增长算法将结节进行三维分割，最终得到分割完整的肺结节区域。
所述的肺结节分割方法，所述步骤A的具体步骤为：
A1；利用各向异性滤波对CT影像进行预处理，通过各向异性滤波去噪不仅能够去除光滑区域的噪声，且能较好保留边缘，有助于之后的检测和分割；
A2；利用最佳阈值法对CT影像分割，得到序列肺实质图像；
A21，在包含序列肺部区域的影像中，求出每一张影像的最大灰度值和最小灰度值，分别记为Gmax和Gmin，令初始阈值T0＝(Gmax+Gmin)/2；
A22，根据阈值T(k)(k＝0，1，2…，k)将CT影像分割为前景和背景，分别求出两者的平均灰度值H1和H2；
A23，求出新阈值T(k+1)＝(H1+H2)/2；计算新的阈值
A24，若T(k)＝T(k+1)，则所得即为阈值；否则转A22，迭代计算；
A25，最终分割去每张影像的肺实质区域，并对包括支气管等在内的其他组织进行去除；
A3；根据A2得到的序列肺实质图像，构建三维肺实质体数据。
所述的肺结节分割方法，所述步骤B的具体步骤为：
B1；在三维图像中将结节、血管、胸膜分别近似看作球体、圆柱体和曲面；为此，构造三个理想模型，在三维空间中分别代表球体，柱状，面状，表达式为：
B2；在三维体数据中，Hessian矩阵描述了体素点附近的局部强度变化的二阶结构，三维图像I(x，y，z)中目标体素点为P，HP为P点的Hessian矩阵，由点P的二阶导数构建成3维图像的Hessian矩阵为一个3×3的实对称矩阵，具有3个实特征值；因此由Hessian矩阵的定义可得点P处的Hessian矩阵表达式为：
B3；通过计算可得，H是三阶对称矩阵，且其六个混合偏导fxy，fxz，fyx，fyz，fzx，fzy的值都为0，fxx，fyy，fzz计算公式为：
Hession的特征值和特征向量可以描述图像I(P)的二阶导数的大小与方向；最大特征值相对应的特征向量代表着P点处曲率的最大方向，最小特征值相对应的特征向量则代表着P点处曲率的最小方向；
B4；对每个Hessian矩阵都按照以下公式(4)进行矩阵分解，其中λ1，λ2，λ3是矩阵的特征值，且|λ1|>＝|λ2|>＝|λ3|，其对应的特征向量分别为e1，e2和e3；
Hession矩阵特征值对应的特征向量两两正交，而且特征向量的方向对应着三维椭球各轴的主方向，特征值的大小对应着各轴的长度，它们共同反映了物体的形状和大小；
B5；根据上述B3，B4计算可得，不同理想模型对应的Hession矩阵的特征值关系总结如下：
球状结构：λ1≈λ2≈λ3＜0
管状结构：λ1≈λ2＜0,λ3≈0
面状结构：λ1＜0,λ2≈λ3≈0。
所述的肺结节分割方法，所述步骤C的具体步骤为：
C1；曲面是描述形状的常用方式，而曲率是曲面最重要的几何特征之一，经典的表面曲率测量，例如高斯曲率和平均曲率，不能很好的指示局部形状，结合两个主曲率提供信息，构建形状指标，曲面的局部的纯几何结构进行量化的指标，不同的数值代表着不同的形状，形状指数定义为：
高斯曲率：K(P)，主曲率：H(p)，
k1(p)，k2(p)为Hessian计算出来的最大、最小特征值；
C2；形状指数由二维推导过程如下：
C21，计算二维形状指数SI的主曲率k1，k2，并将其单位化后变化到极坐标：
C22，在极坐标系中由角度关系可得：
C23，由C22中的角度关系可知：
C24，最终θ1的值为：
C25，最终可以推导出二维形状指数：
C3；形状指数由二维扩展到三维过程如下：
C31，根据C2形状指数二维的推导计算可得
C32，根据arctan和arccos关系可知：
C33，在三维图像中，球坐标是三维坐标系的一种，用以确定三维空间中点、线、面以及体的位置，它以坐标原点为参考点，由方位角、仰角和距离构成如球坐标；设P(x，y，z)为空间内一点，球坐标系(r，θ，
)与直角坐标系(x，y，z)的转换关系如下面所示：
C34，在三维体素中一点P(x，y，z)，有三个主曲率k1，k2，k3，将其单位化后变化到球坐标为：
C35，根据C24计算可得：
C36，最终可得三维形状指数：
所述的肺结节分割方法，所述步骤D的具体步骤为：
D1；由于肺结节灰度值呈Gaussian分布，对不同尺度不同类型的疑似肺结节，选取不同的尺度参数σ进行多尺度增强；根据高斯模版具有三维旋转对称性可得：σx＝σy＝σz＝σ，三维高斯模板函数为：
D2；在三维体数据上使用服从高斯分布的球状模型模拟结节，利用高斯函数来建立肺结节类球型b(x，y，z)数学模型定义为：
D3；将Hessian矩阵的差分运算与高斯卷积相结合，通过改变高斯函数的标准偏移量来获得不同尺度σ下的线条滤波图像；
D31；三维高斯函数的表达式为：
D32；根据高斯函数的卷积性质得图像I为：
D4；带检测的肺结节直径范围为[d0，d1]，为了检测在这个范围内的所有物体，将高斯滤波器的尺度范围设定在[d0/4，d1/4]之间；然后选取N个不同的尺度参数σ值，分别对图像进行卷积和增强运算；尺度参数σ计算公式为：
D5；将得到的高斯滤波图像结合Hessian矩阵的表达式为：
D6；根据所述Hession矩阵的特征值和所述三维形状指数，最终构建类球形滤波器结节检测函数V为：
D7；结节检测的具体步骤为
D71，使用各向异性滤波对序列图像进行去噪；
D72，根据最佳阈值法分割序列肺实质区域，并构建三维体数据I；
D73，确定经验范围内的高斯滤波尺度σ以及σ个数，并计算每个σ的值；
D74，对每一个尺度σn分别进行步骤D75-D710；
D75，使用高斯滤波器平滑中三维体数据Iin；
D76，对每一个体素点分别进行步骤D77-D79；
D77，计算Hessian矩阵以及其对应的特征值λ1，λ2，λ3；
D78，计算三维形状指数:DSI；
D79，计算类球形结节检测函数V，并输出三维体数据Iout；
D710，终止体数据循环；
D711，终止尺度循环；
D712，最终输出不同尺度检测到的结节区域位置。
所述的肺结节分割方法，所述步骤E的具体步骤为：
E1，将通过三维类球形滤波器-3D shape index结节检测肺实质三维体数据得到的结节区域的多体素点坐标作为区域增长的多种点；
E2，计算肺实质体数据中体素灰度值的平均值：
和标准差：
E3，根据指定的惩罚因子，乘以标准差，I(X)∈[m-fσ,m+fσ]计算以期望为中心的灰度值范围，X是图像I中的像素点，m和σ分别是当前区域灰度的平均值和标准差；如果邻域的灰度值位于这个范围内是话就被包含在种子区域，否自就被排除；
E4，更新各部分肺实质体数据体素权重，遍历所有体素点；
E5，最终完成结节的三维分割。
与现有的肺结节分割技术相比：
1、本发明的目的在于克服现有检测和分割技术中存在的缺陷，提供一种简单、自动化的肺结节检测分割方法。
2、利用本发明技术，对肺结节可以准确的检测以及分割，具有稳定性和可再现性特点。
3、为后续肺结节的良恶性分类诊断提供了良好的基础。
附图说明
图1是本发明的肺结节检测算法总体流程图；
图2是本发明的序列肺实质分割区域；
图3是本发明的结节、血管、胸膜的三维效果图及其对应的数学模型和三维结构；
图4是本发明的SI值的意义；
图5是本发明的极坐标和球坐标的对应和转换关系；
图6是本发明不同类型结节在不同尺度下的检测效果；
图7是本发明不同类型的结节检测实验结果；
具体实施方式
以下结合具体实施例，对本发明进行详细说明。
参照图1，主要流程包括：CT影像预处理：各向异性滤波去噪，和序列肺实质分割、类球形滤波器构建：Hession矩阵特征值计算，3D形状指数函数的构建，类球形滤波器的构建、三维肺结节分割：结合区域增长的置信连接，三维肺结节分割等步骤。本发明方法的具体实施方式如下：
A.三维肺实质区域体数据的构建
A1；利用各向异性滤波对CT影像进行预处理，通过各向异性去噪不仅能够去除光滑区域的噪声，且能较好保留边缘，有助于之后的分割和检测。
A2；利用最佳阈值法对CT影像分割，得到序列肺实质图像。
第1步，在包含序列肺部区域的影像中，求出每一张影像的最大灰度值和最小灰度值，分别记为Gmax和Gmin，令初始阈值T0＝(Gmax+Gmin)/2；
第2步，根据阈值T(k)(k＝0，1，2…，k)将CT影像分割为前景和背景，分别求出两者的平均灰度值H1和H2；
第3步，求出新阈值T(k+1)＝(H1+H2)/2；计算新的阈值
第4步，若T(k)＝T(k+1)，则所得即为阈值；否则转2，迭代计算。
第5步，最终分割去每张影像的肺实质区域，并对支气管等其他组织进行去除。
A3；根据A2得到的序列肺实质图像，构建三维肺实质体数据。
参考图2，图2为经过上述操作得到的序列肺实质区域。
B.分析结节的三维特征，计算Hession矩阵特征值
在检测疑似肺结节时，由于结节、血管、胸膜等都是三维的实体，而且在形态上与标准的球、圆柱和曲面有一定的相似性，其中结节、血管、胸膜等三维实体对应的理想球体模型、圆柱体模型和曲面模型。参考图3，图3为结节、血管、胸膜的三维效果图及其对应的数学模型和三维结构，通过对典型的球、圆柱和曲面的基本特征的分析，选择能够区分球、圆柱、和曲面的特征提取出来，并将这些特征引入到结节、血管和胸膜的三维模型中来，作为区分结节、血管等的重要依据。具体步骤如下：
B1；在三维图像中将结节、血管、胸膜分别近似看作球体、圆柱体和曲面。为此，构造三个理想模型，在三维空间中分别代表球体，柱状，面状，表达式为：
B2；在三维体数据中，Hessian矩阵描述了体素点附近的局部强度变化的二阶结构，三维图像I(x，y，z)中目标体素点为P，HP为P点的Hessian矩阵，由点P的二阶导数构建成3维图像的Hessian矩阵为一个3×3的实对称矩阵，具有3个实特征值。因此由Hessian矩阵的定义可得点P处的Hessian矩阵表达式为：
B3；通过计算可得，H是三阶对称矩阵，且其六个混合偏导fxy，fxz，fyx，fyz，fzx，fzy的值都为0，fxx，fyy，fzz计算公式为：
Hession的特征值和特征向量可以描述图像I(P)的二阶导数的大小与方向。最大特征值相对应的特征向量代表着P点处曲率的最大方向，最小特征值相对应的特征向量则代表着P点处曲率的最小方向；
B4；对每个Hessian矩阵都可以按照以下方式公式(4)进行矩阵分解，其中λ1，λ2，λ3是矩阵的特征值，且|λ1|>＝|λ2|>＝|λ3|，其对应的特征向量分别为e1，e2和e3。
Hession矩阵特征值对应的特征向量两两正交，而且特征向量的方向对应着三维椭球各轴的主方向，特征值的大小对应着各轴的长度，它们共同反映了物体的形状和大小。
B5；根据上述B3，B4计算可得，不同理想模型对应的海森矩阵特征值关系总结如下：
a；球状结构：λ1≈λ2≈λ3＜0
b；管状结构：λ1≈λ2＜0,λ3≈0
c；面状结构：λ1＜0,λ2≈λ3≈0
C.构建三维形状指数
单独使用Hession矩阵特征值来构建增强滤波器包含较多的假阳性，而且仅用最大和最小曲率构建的SI，在二维图像上检测效果较好，没有充分利用Hession矩阵的特征值，面对复杂的肺部CT影像并不能很好的检测肺结节，而且存在较高的假阳性。三维体数据中，肺结节检测根据三维特征可以很好的区分结节和血管以及其他组织；因此本文将形状指数SI扩展到三维体数据，通过体数据计算三维SI。使用三维曲率，来更加完善的表示三维特征，能够较为准确的区分结节和血管，构建三维形状指数具体步骤如下：
C1；曲面是描述形状的常用方式，而曲率是曲面最重要的几何特征之一，经典的表面曲率测量，例如高斯曲率和平均曲率，不能很好的指示局部形状，结合两个主曲率提供信息，构建形状指标，曲面的局部的纯几何结构进行量化的指标，不同的数值代表着不同的形状，形状指数定义为公式(5)，所代表的值如图4所示：
高斯曲率：K(P)主曲率：H(p)，
k1(p)，k2(p)为Hessian计算出来的最大最小特征值。
C2；形状指数由二维推导过程如下：
第1步，计算二维形状指数SI的主曲率k1，k2，并将其单位化后变化到极坐标，参照图5(a)：
第2步，参照图5(b)在极坐标系中由角度关系可得：
第3步，由第2步中的角度关系可知：
第4步，最终θ1的值为：
第5步，最终可以推到出二维形状指数：
C3；形状指数由二维扩展到三维过程如下：
第1步，根据C2形状指数二维的推导计算可得
第2步，根据arctan和arccos关系可知：
第3步，在三维图像中，球坐标是三维坐标系的一种，用以确定三维空间中点、线、面以及体的位置，它以坐标原点为参考点，由方位角、仰角和距离构成如球坐标。参照图5(c)，设P(x，y，z)为空间内一点，球坐标系(r，θ，
)与直角坐标系(x，y，z)的转换关系如下面所示：
第4步，在三维体素中一点P(x，y，z)，有三个主曲率k1，k2，k3，将其单位化后变化到球坐标为参照图5(d)：
第5步，根据二维推导第4步计算可得：
第6步，最终可得三维形状指数：
D.多尺度的三维类球形滤波器
在肺部CT图像中，结节的尺寸大小具有不确定性，而且存在许多图像噪声，而体素的二次偏导计算过程对图像噪声具有很强的敏感性，因此，如果将增强滤波器直接应用在图像上，将不会产生好的结果。为了有效检测出不同大小的结节，本发明采用的是一种基于高斯函数的多尺度滤波方法。利用高斯函数对图像进行卷积运算，去除图像中的噪声同时平滑图像，再对图像进行多尺度增强，可以增强不同尺度大小的结节图像，具体步骤如下：
D1；由于肺结节灰度值呈Gaussian分布，对不同尺度不同类型的疑似肺结节，选取不同的尺度参数σ进行多尺度增强。根据高斯模版具有三维旋转对称性可得：σx＝σy＝σz＝σ，三维高斯模板函数为：
D2；在三维体数据上使用服从高斯分布的球状模型模拟结节，利用高斯函数来建立肺结节类球型b(x，y，z)数学模型定义为：
D3；将Hessian矩阵的差分运算与高斯卷积相结合，通过改变高斯函数的标准偏移量来获得不同尺度σ下的线条滤波图像；
a；三维高斯函数的表达式为：
b；根据高斯函数的卷积性质得图像I为：
D4；带检测的肺结节直径范围为[d0，d1]，为了检测在这个范围内的所有物体，将高斯滤波器的尺度范围设定在[d0/4，d1/4]之间。然后选取N个不同的σ值，分别对图像进行卷积和增强运算。尺度计算公式为：
D5；将得到的高斯滤波图像结合Hessian矩阵的表达式为：
D6；根据权利要求3 B5中得到的Hession矩阵的特征值和权利要求4C3第六步中推导出来三维形状指数，最终构建类球形滤波器结节检测函数V为：
D7；结节检测的具体步骤为
第1步，使用各向异性滤波对序列图像进行去噪。
第2步，根据最佳阈值法分割序列肺实质区域，并构建三维体数据I。
第3步，确定经验范围内的高斯滤波尺度σ以及σ个数，并计算每个σ的值。
第4步，对每一个尺度σn分别进行步骤5-10。
第5步，使用高斯滤波器平滑中三维体数据Iin。
第6步，对每一个体素点分别进行步骤7-9。
第7步，计算Hessian矩阵以及其对应的特征值λ1，λ2，λ3。
第8步，计算三维形状指数:DSI。
第9步，计算类球形结节检测函数V，并输出三维体数据Iout。
第10步，终止体数据循环。
第11步，终止尺度循环。
第12步，最终输出不同尺度检测到的结节区域位置。
E.三维结节分割
三维分割是针对的三维体数据算法，算法在三维体数据上直接进行分割，执行1次即可完成整体的分割，效率高，并且充分利用切片之间的空间、纹理等信息，提高分割的准确性，具体步骤如下：
第1步，将通过三维类球形滤波器-3D shape index结节检测肺实质三维体数据得到的结节区域的多体素点坐标作为区域增长的多种点；
第2步，计算肺实质体数据中体素灰度值的平均值：
和标准差：
第3步，根据指定的惩罚因子，乘以标准差，I(X)∈[m-fσ,m+fσ]计算以期望为中心的灰度值范围，X是图像I中的像素点，m和σ分别是当前区域灰度的平均值和标准差。如果邻域的灰度值位于这个范围内是话就被包含在种子区域，否自就被排除；
第4步，更新各部分肺实质体数据体素权重，遍历所有体素点；
第5步，最终完成结节的三维分割。
参考图6是本发明对其中不同类型肺结节的多尺度检测结果。其中，第一列中包含孤立性肺结节，第二列中包含血管粘连结节，第三列中包含胸膜牵拉结节，第三列中包含磨玻璃结节，第(5)列中包含多个结节类型；(a)列为三维肺实质体数据，(b)列为高斯尺度σ为5时的效果，(c)列为高斯尺度σ为5.5时的效果，(d)列为高斯尺度σ为6时的效果，(e)列为高斯尺度σ为6.5时的效果，(f)列为高斯尺度σ为7时的效果。
参考图7为本文方法对四类肺结节的检测结果和分割结果。(1)列是医生标注的结节区域；(2)列是肺实质区域的三维重建部分；(3)列是本文构建的球型滤波器检测结果，并将其对应到原来的三维重建区域；(4)列是结合置信连接的区域生长的分割效果图；(5)列是结节最终的分割结果。(a)行是孤立性结节的检测和分割效果；(b)行是血管牵拉结节的检测和分割效果；(c)行是胸膜牵拉结节的检测和分割效果；(d)行是磨玻璃结节的检测和分割效果；(e)行是多结节检测和分割效果。
本发明可以相对完整的检测出不同类型的结节，不仅对较易检测的孤立性肺结节有很好的检测效果，而且对胸膜牵拉型结节、血管粘连型结节和磨玻璃结节都有很好的检测效果。尤其是对磨玻璃这种特征不明显难以检测的结节也有很好的检测效果。
应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。
Claims (3)
1.一种基于Hession矩阵和三维形状指数的肺结节分割方法，其特征在于，包括以下步骤：
步骤A；肺实质三维体数据的构建：首先利用各向异性滤波对CT影像进行预处理，去除噪声图像的影像，然后利用最佳阈值法对图像中的肺实质区域进行提取，按照序列依次提取单张CT影像的肺实质区域，最终构建肺实质三维体数据；
步骤B；在三维图像中将结节、血管、胸膜分别近似看作球体、圆柱体和曲面，通过分析典型的球、圆柱和曲面的基本特征，构造三个理想的数学模型，并计算肺实质三维体数据中每个体素点的Hession矩阵特征值，并分析每个模型的三维特征；
步骤C；分析三维模型的形状特征，利用肺结节的三维特征将其与包括血管在内的疑似区域进行区分，通过将二维形状指数扩展到三维体数据上，并结合Hession矩阵特征值对三维形状指数进行构建，最终构建能够较为准确的区分结节和血管的三维形状指数；
步骤D；肺结节灰度值呈Gaussian分布，不同尺度空间上的高斯滤波函数对于复杂的医学影像具有高效的增强效果，对不同尺度不同类型的疑似肺结节，选取不同的尺度参数σ进行多尺度增强，得到多尺度的三维类球形滤波器-3D shape index结节检测函数；
步骤E；利用步骤D中得到的3D shape index结节检测函数对肺实质三维体数据进行结节检测，并利用检测得到结节区域作为区域增长的多种子点，结合置信连接的区域增长算法将结节进行三维分割，最终得到分割完整的肺结节区域；
所述步骤A的具体步骤为：
A1；利用各向异性滤波对CT影像进行预处理，通过各向异性滤波去噪不仅能够去除光滑区域的噪声，且能较好保留边缘，有助于之后的检测和分割；
A2；利用最佳阈值法对CT影像分割，得到序列肺实质图像；
A21，在包含序列肺部区域的影像中，求出每一张影像的最大灰度值和最小灰度值，分别记为Gmax和Gmin，令初始阈值T0＝(Gmax+Gmin)/2；
A22，根据阈值T(k),k＝0，1，2…,将CT影像分割为前景和背景，分别求出两者的平均灰度值H1和H2；
A23，求出新阈值T(k+1)＝(H1+H2)/2；计算新的阈值；
A24，若T(k)＝T(k+1)，则所得即为阈值；否则转A22，迭代计算；
A25，最终分割出每张影像的肺实质区域，并对包括支气管在内的其他组织进行去除；
A3；根据A2得到的序列肺实质图像，构建三维肺实质体数据；
所述步骤B的具体步骤为：
B1；在三维图像中将结节、血管、胸膜分别近似看作球体、圆柱体和曲面；为此，构造三个理想模型，在三维空间中分别代表球体，柱状，面状，表达式为：
B2；在三维体数据中，Hessian矩阵描述了体素点附近的局部强度变化的二阶结构，三维图像I(x，y，z)中目标体素点为P，HP为P点的Hessian矩阵，由点P的二阶导数构建成3维图像的Hessian矩阵为一个3×3的实对称矩阵，具有3个实特征值；因此由Hessian矩阵的定义可得点P处的Hessian矩阵表达式为：
B3；通过计算可得，HP是三阶对称矩阵，且其六个混合偏导fxy，fxz，fyx，fyz，fzx，fzy的值都为0；fxx，fyy，fzz计算公式为：
Hession的特征值和特征向量可以描述图像I(P)的二阶导数的大小与方向；最大特征值相对应的特征向量代表着P点处曲率的最大方向，最小特征值相对应的特征向量则代表着P点处曲率的最小方向；
B4；对每个Hessian矩阵都按照以下公式(4)进行矩阵分解，其中λ1，λ2，λ3是矩阵的特征值，且|λ1|>＝|λ2|>＝|λ3|，其对应的特征向量分别为e1，e2和e3；
Hession矩阵特征值对应的特征向量两两正交，而且特征向量的方向对应着三维椭球各轴的主方向，特征值的大小对应着各轴的长度，它们共同反映了物体的形状和大小；
B5；根据上述B3，B4计算可得，不同理想模型对应的Hession矩阵的特征值关系总结如下：
球状结构：λ1≈λ2≈λ3＜0
管状结构：λ1≈λ2＜0,λ3≈0
面状结构：λ1＜0,λ2≈λ3≈0；
所述步骤C的具体步骤为：
C1；形状指数定义为：
高斯曲率：K(P)，主曲率：H(p)，
k1(p)，k2(p)为Hessian计算出来的最大、最小特征值；
C2；形状指数由二维推导过程如下：
C21，计算二维形状指数SI的主曲率k1，k2，并将其单位化后变化到极坐标：
C22，在极坐标系中由角度关系可得：
C23，由C22中的角度关系可知：
C24，最终θ1的值为：
C25，最终可以推导出二维形状指数：
C3；形状指数由二维扩展到三维过程如下：
C31，根据步骤C2中形状指数的二维推导过程，计算可得
C32，根据arctan和arccos关系可知：
C33，在三维图像中，球坐标是三维坐标系的一种，用以确定三维空间中点、线、面以及体的位置，它以坐标原点为参考点，由方位角、仰角和距离构成球坐标；设P(x，y，z)为空间内一点，球坐标系
与直角坐标系(x，y，z)的转换关系如下面所示：
C34，在三维体素中一点P(x，y，z)，有三个主曲率k1，k2，k3，将其单位化后变化到球坐标为：
C35，根据C24计算可得：
C36，最终可得三维形状指数：
2.根据权利要求1所述的肺结节分割方法，其特征在于，所述步骤D的具体步骤为：
D1；由于肺结节灰度值呈Gaussian分布，对不同尺度不同类型的疑似肺结节，选取不同的尺度参数σ进行多尺度增强；根据高斯模板具有三维旋转对称性可得：σx＝σy＝σz＝σ，三维高斯模板函数为：
D2；在三维体数据上使用服从高斯分布的球状模型模拟结节，利用高斯函数来建立肺结节类球型b(x，y，z)数学模型定义为：
D3；将Hessian矩阵的差分运算与高斯卷积相结合，通过改变高斯函数的标准偏移量来获得不同尺度σ下的线条滤波图像；
D31；三维高斯函数的表达式为：
D32；根据高斯函数的卷积性质得图像I为：
D4；带检测的肺结节直径范围为[d0，d1]，为了检测在这个范围内的所有物体，将高斯滤波器的尺度范围设定在[d0/4，d1/4]之间；然后选取N个不同的尺度参数σ值，分别对图像进行卷积和增强运算；尺度参数σ计算公式为：
D5；将得到的高斯滤波图像结合Hessian矩阵的表达式为：
D6；根据所述Hession矩阵的特征值和所述三维形状指数，最终构建类球形滤波器结节检测函数V为：
D7；结节检测的具体步骤为：
D71，使用各向异性滤波对序列图像进行去噪；
D72，根据最佳阈值法分割序列肺实质区域，并构建三维体数据I；
D73，确定经验范围内的高斯滤波尺度σ以及σ个数，并计算每个σ的值；
D74，对每一个尺度σn分别进行步骤D75-D710；
D75，使用高斯滤波器平滑中三维体数据Iin；
D76，对每一个体素点分别进行步骤D77-D79；
D77，计算Hessian矩阵以及其对应的特征值λ1，λ2，λ3；
D78，计算三维形状指数:DSI；
D79，计算类球形结节检测函数V，并输出三维体数据Iout；
D710，终止体数据循环；
D711，终止尺度循环；
D712，最终输出不同尺度检测到的结节区域位置。
3.根据权利要求2所述的肺结节分割方法，其特征在于，所述步骤E的具体步骤为：
E1，将通过三维类球形滤波器-3D shape index结节检测肺实质三维体数据得到的结节区域的多体素点坐标作为区域增长的多种点；
E2，计算肺实质体数据中体素灰度值的平均值：
和标准差：
E3，根据指定的惩罚因子，乘以标准差，I(X)∈[m-cσ,m+cσ]计算以期望为中心的灰度值范围，X是图像I中的像素点，m和σ分别是当前区域灰度的平均值和标准差，如果邻域的灰度值位于这个范围内就被包含在种子区域，否则就被排除；
E4，更新各部分肺实质体数据体素权重，遍历所有体素点；
E5，最终完成结节的三维分割。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710464681.0A CN107274399B (zh) | 2017-06-19 | 2017-06-19 | 一种基于Hession矩阵和三维形状指数的肺结节分割方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710464681.0A CN107274399B (zh) | 2017-06-19 | 2017-06-19 | 一种基于Hession矩阵和三维形状指数的肺结节分割方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107274399A CN107274399A (zh) | 2017-10-20 |
| CN107274399B true CN107274399B (zh) | 2020-03-13 |
Family
ID=60069031
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710464681.0A Active CN107274399B (zh) | 2017-06-19 | 2017-06-19 | 一种基于Hession矩阵和三维形状指数的肺结节分割方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107274399B (zh) |
Families Citing this family (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108109159B (zh) * | 2017-12-21 | 2020-07-28 | 东北大学 | 一种基于hessian矩阵和区域增长相结合的视网膜血管分割系统 |
| CN108447066B (zh) * | 2018-03-19 | 2021-03-09 | 青岛海信医疗设备股份有限公司 | 胆道图像分割方法、终端、存储介质 |
| CN108510478B (zh) * | 2018-03-19 | 2021-01-26 | 青岛海信医疗设备股份有限公司 | 肺气道图像分割方法、终端、存储介质 |
| CN109191436A (zh) * | 2018-08-15 | 2019-01-11 | 复旦大学 | 基于视觉显著性谱残差方法的低剂量ct肺结节检测算法 |
| CN109410170B (zh) * | 2018-09-14 | 2022-09-02 | 东软医疗系统股份有限公司 | 图像数据处理方法、装置及设备 |
| CN109583444B (zh) * | 2018-11-22 | 2021-08-10 | 博志生物科技有限公司 | 空洞区域定位方法、装置及计算机可读存储介质 |
| CN111353328B (zh) * | 2018-12-20 | 2023-10-24 | 核动力运行研究所 | 一种超声三维体数据在线显示及分析方法 |
| CN111435532B (zh) * | 2019-01-14 | 2021-06-22 | 湖南大学 | 一种数字图像中树状结构末梢点的检测方法 |
| CN109727260A (zh) * | 2019-01-24 | 2019-05-07 | 杭州英库医疗科技有限公司 | 一种基于ct影像的三维肺叶分割方法 |
| CN110599482B (zh) * | 2019-09-18 | 2022-06-03 | 青岛海信医疗设备股份有限公司 | 肺裂曲面确定方法、装置、电子设备和存储介质 |
| CN112184659B (zh) * | 2020-09-24 | 2023-08-25 | 上海健康医学院 | 一种肺部图像处理方法、装置及设备 |
| CN114419077B (zh) * | 2022-03-31 | 2022-07-05 | 深圳市一图智能科技有限公司 | 一种基于ct图像的肺支气管自动分割方法及系统 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1836240A (zh) * | 2003-08-13 | 2006-09-20 | 美国西门子医疗解决公司 | Cad(计算机辅助决策)支持系统和方法 |
| US9014456B2 (en) * | 2011-02-08 | 2015-04-21 | University Of Louisville Research Foundation, Inc. | Computer aided diagnostic system incorporating appearance analysis for diagnosing malignant lung nodules |
-
2017
- 2017-06-19 CN CN201710464681.0A patent/CN107274399B/zh active Active
Non-Patent Citations (5)
| Title |
|---|
| "3D shape analysis for early diagnosis of malignant lung nodules";Ayman El-Baz等;《Med Image Comput Comput Assist Interv. 2011》;20111231;第772-783页 * |
| "Model-based quantitation of 3-D magnetic resonance angiographic images";Alejandro F. Frangi等;《IEEE Transactions on Medical Imaging》;19991031;第18卷(第10期);第946-956页 * |
| "Surface shape and curvature scales";Jan J Koenderink等;《Image and Vision Computing》;19921031;第10卷(第8期);第557-564页 * |
| "基于PET/CT的孤立性肺结节的自动分割方法";强彦等;《清华大学学报（自然科学版）》;20130215;第53卷(第2期);第200-204页 * |
| "基于分层模版种子点的分水岭分割孤立性肺结节的方法";赵涓涓等;《清华大学学报（自然科学版）》;20140715;第54卷(第7期);第910-916页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107274399A (zh) | 2017-10-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107274399B (zh) | 一种基于Hession矩阵和三维形状指数的肺结节分割方法 | |
| US11710242B2 (en) | Methods and systems for image segmentation | |
| Nailon | Texture analysis methods for medical image characterisation | |
| Liu et al. | Deformable registration of cortical structures via hybrid volumetric and surface warping | |
| CN108052942B (zh) | 一种飞机飞行姿态的视觉图像识别方法 | |
| US9129390B2 (en) | Method of segmenting anatomic entities in 3D digital medical images | |
| US6766043B2 (en) | Pleural nodule detection from CT thoracic images | |
| CN105719295B (zh) | 一种基于三维超体素的颅内出血区域分割方法及系统 | |
| WO2001078005A2 (en) | System and method for three-dimensional image rendering and analysis | |
| CN101116089A (zh) | 在三维图像容积中检测息肉的方法 | |
| CN107633514B (zh) | 一种肺结节周边血管量化评估系统及方法 | |
| CN104021547A (zh) | 肺部 ct 的三维配准方法 | |
| CN110738701A (zh) | 一种肿瘤三维定位系统 | |
| CN109300113A (zh) | 一种基于改进凸包方法的肺结节辅助检测系统及方法 | |
| Li et al. | Segmentation of pulmonary nodules using adaptive local region energy with probability density function-based similarity distance and multi-features clustering | |
| EP2006802A1 (en) | Method of constructing a grey value model and/or a geometric model of an anatomic entity in a 3D digital medical image | |
| Zhou et al. | Automatic segmentation of liver from CT scans with CCP–TSPM algorithm | |
| CN109872353B (zh) | 基于改进迭代最近点算法的白光数据与ct数据配准方法 | |
| Cheng et al. | A new strategy for boundary surface detection in 3D biomedical images | |
| Bin et al. | Rapid multimodal medical image registration and fusion in 3D conformal radiotherapy treatment planning | |
| Zhang et al. | Automatic liver segmentation from CT scans based on a statistical shape model | |
| Zhang et al. | Feature Point Extraction and Motion Tracking of Cardiac Color Ultrasound under Improved Lucas–Kanade Algorithm | |
| TWI548401B (zh) | 血管三維結構重建方法 | |
| CN117876833A (zh) | 一种用于机器学习的肺部ct图像特征提取方法 | |
| Zhou et al. | 3D-3D tubular organs registration based on bifurcations for the CT images |
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 |