CN101763644A - 肺结节三维分割与特征提取方法及系统 - Google Patents

Abstract

本发明属于医学图像的计算机分析技术的应用领域，为一种基于胸部CT的肺结节三维分割与特征提取方法及系统。本发明采用基于三维空间的水平集方法在三维感兴趣区域内分割出肺结节，并结合肺结节边界曲面和感兴趣区域的特性提取结节特征。本发明包括肺区域分割模块、三维感兴趣区域提取模块、肺结节三维分割模块、结节特征提取及量化模块和输出模块。本发明充分利用了HRCT图像三维空间信息，可以减少结节周围正常组织对其分割和特征提取的干扰，提高了特征提取的准确性。

Description

肺结节三维分割与特征提取方法及系统

技术领域

本发明属于医学图像的计算机分析技术的应用领域，具体涉及一种基于胸部CT的肺结节三维分割与特征提取方法及系统。

背景技术

随着空气污染造成环境的不断恶化及吸烟人口的大量增加，肺部疾病的发病率及死亡率也在逐年上升。在全世界大部分国家，肺癌是男性、女性癌症死亡排名第一的原因。而在我国，过去的5年中肺癌病患者增加了约12万；每4个癌症死亡患者中就有1人是肺癌患者。由于肺部发生的病变在早期难以被人察觉，等到有临床症状出现而就医时，大多数患者已到达了病程的中晚期，治愈的机率将大大降低。而肺结节是肺癌早期最常见的表现形式，在肺癌计算机辅助诊断技术中也主要是对肺结节的检测与量化方法进行研究。

在肺结节的检测中，CT被认为是最好的影像学手段，肺结节所具有的三维特征，比如形态近似球形、灰度近似空间正态分布的特点，在高分辨率CT(High Resolution Computed Tomography，HRCT)中更为明显。虽然，目前许多学者在肺结节的检测和特征量化方面做了大量研究，相关研究也已取得了一定的成果，但大部分研究都集中在对二维CT横断面图像(原始的CT图像)进行处理，没有充分利用CT三维空间信息，而且大部分特征提取方法或者只考虑感兴趣区域的特性而忽略了肺结节自身的特性，或者只考虑肺结节自身的特性而忽略了感兴趣区域的特性，没有充分结合两者的特性。如公开号CN101556650A的专利文献在没有获得肺结节边界的情况下，直接从感兴趣区域提取特征，如分叶征、空洞征等，但感兴趣区域内还包含了其他正常组织，这些组织特性对感兴趣区域的特征提取会造成干扰。

由于支气管树通过主支气管入口与心脏、胸廓等高密度外围组织相连，在对肺实质进行分割时，肺的表面会沿着与外围高密度组织粘连的主支气管通道发生凹陷，位于凹陷内的支气管段将会被漏分，因此使用几何凸点搭桥法填补整个肺表面的凹陷，从而得到完整的肺内区域(具体参见KawataY.，N.Niki，H.Ohmatsu，et al.“Computerized analysis of 3-D pulmonarynodule images in surrounding and internal structure feature spaces，”2001International Conference on Image Processing，889-892(2001).)。

发明内容

本发明的目的在于提供一种肺结节三维分割与特征提取方法，该方法能够从大量的二维胸部CT横断面图像中利用三维空间信息分割出肺结节，能够准确可靠地提取及量化肺结节特征；本发明还提供了实施该方法的系统，系统稳定可靠。

本发明提供的一种肺结节三维分割与特征提取方法，其特征在于，该方法包括下述步骤：

第1步在输入的所有二维胸部CT横断面图像上自动分割出肺实质后，使用几何凸点搭桥法修补肺表面的凹陷，将肺结节所在的肺内区域与肺外区域分隔开；

第2步构建三维空间中不同尺度的三种形状选择性增强滤波器，分别仅对球形、管状及面状物体具有响应。利用这三种形状选择性增强滤波器分别初步甄别出由所有二维胸部CT横断面图像构建的三维CT图像中的疑似肺结节、疑似粗血管及疑似厚支气管壁，提取出肺内区域中候选结节区域，再利用Fisher线性分类器对候选结节区域的特征值进行分类从而得到三维感兴趣区域；

第3步在第2步提取的三维感兴趣区域内使用基于三维空间的水平集方法分割出肺结节，确定肺结节边界曲面；

第4步根据第3步得到的肺结节边界曲面的信息，结合感兴趣区域的总体特性，提取相关特征并量化特征值；

第5步根据第3步得到的肺结节边界曲面信息，将肺结节的位置和范围以三维形式显示，同时以用户需要的形式存储及显示第4步得到的特征及其量化值。

本发明提供的一种肺结节三维分割与特征提取系统，其特征在于：该系统包括肺区域分割模块、三维感兴趣区域提取模块、肺结节三维分割模块、肺结节特征提取及量化模块和输出模块；

肺区域分割模块用于将输入的所有二维胸部CT横断面图像中肺实质分割出来后，使用几何凸点搭桥法修补肺表面的凹陷，并将分割出的肺内区域传送给三维感兴趣区域提取模块；

三维感兴趣区域提取模块用于从肺内区域中提取候选结节区域，利用Fisher线性分类器对候选结节区域的特征值进行分类从而获取三维感兴趣区域，并将感兴趣区域的边界信息传送给肺结节三维分割模块和结节特征提取及量化模块；

肺结节三维分割模块使用基于三维空间的水平集方法从三维感兴趣区域中将肺结节分割出来，并将得到的肺结节边界曲面信息传送给结节特征提取及量化模块和输出模块；

结节特征提取及量化模块用于接收肺结节三维分割模块获得的肺结节的边界曲面信息，结合三维感兴趣区域提取模块得到的三维感兴趣区域的总体特性，提取肺结节的相关特征，计算特征值，并将其特征及其量化值传送给输出模块；

输出模块用于接收肺结节三维分割模块得到的肺结节边界曲面信息，并将其位置和范围以三维形式显示，同时接收结节特征提取及量化模块得到的特征及其量化值，将其结果根据需要存储及显示。

本发明首先批量地分割出所有二维胸部CT横断面图像中肺内区域，利用三维空间信息提取感兴趣区域，然后在感兴趣区域内对肺结节进行三维分割，根据得到的肺结节边界曲面信息，结合感兴趣区域的总体特性提取并量化结节特征，最后以三维形式显示结节位置和范围，并根据用户需要存储及显示肺结节特征值。总之，本发明通过一个肺结节三维分割与特征提取方法及系统自动地对大量的二维胸部CT横断面图像中的肺结节进行三维分割和检测，显示肺结节位置及计算一系列相关特征值。

本发明在感兴趣区域提取时，根据肺结节、粗血管、厚支气管壁或者其交叉区域在三维空间中的形状特点，使用了三维空间的多尺度形状选择性增强滤波器，可以减少粗血管、厚支气管壁或者其交叉区域对候选结节区域的干扰，能降低假阳性结节出现的概率。同时，本发明结合了肺结节边界曲面和感兴趣区域的特性进行特征提取，既考虑到肺结节所在区域的总体特性，又减少了肺结节周围正常组织对其特征的干扰，提高了特征提取的准确性，此外还可以提取肺结节的体积、边界清晰度和钙化含量等特征。

附图说明

图1为本发明肺结节三维分割与特征提取方法的流程图；

图2为本发明肺结节三维分割与特征提取系统的功能结构图；

图3为本发明实施例的感兴趣区域提取的特征列表；

图4为本发明实施例中肺结节三维分割的过程图；

图5为本发明实施例中毛刺征提取时肺结节边界正方形感兴趣区域的选取示意图。

具体实施方式

如图1所示，本发明提供的肺结节三维分割与特征提取方法，包括以下步骤：

(1)在输入的所有二维胸部CT横断面图像上自动分割出肺实质后，使用几何凸点搭桥法修补肺表面的凹陷，将肺结节所在的肺内区域与肺外区域分隔开。

采用阈值分割方法对每层二维胸部CT横断面图像进行自动分割，即根据OTSU阈值法得到的肺实质分割阈值T₁，对图像进行二值化处理，其中肺实质内充满空气、密度小，它与背景区域的CT值都小于T₁，则被转换成二值图像中黑色区域；而脂肪、肌肉等组织密度大，其CT值都高于T₁，则被转换成二值图像中白色区域。

由于支气管树通过主支气管入口与心脏、胸廓等高密度外围组织相连，在通过上述阈值分割方法对肺实质进行分割时，肺的表面会沿着与外围高密度组织粘连的主支气管通道发生凹陷，位于凹陷内的支气管段将会被漏分。本发明使用几何凸点搭桥法修补每层二维胸部CT横断面图像的肺边界，来对整个肺表面的凹陷进行填补，从而得到完整的肺内区域。

(2)构建三维空间中不同尺度的三种形状选择性增强滤波器，分别仅对球形、管状及面状物体具有响应。利用这三种形状选择性增强滤波器分别初步甄别出由所有二维胸部CT横断面图像构建的三维CT图像中的疑似肺结节、疑似粗血管及疑似厚支气管壁，提取出肺内区域中候选结节区域，再利用Fisher线性分类器对候选结节区域的特征值进行分类从而得到三维感兴趣区域。

一些呈毛玻璃状的结节与周围区域对比度很低，并且大部分肺结节跟肺部血管和气管的灰度值一样以及肺结节可能与血管和气管粘连在一起等原因，使得不能通过简单的阈值分割得到感兴趣区域。考虑到三维空间中肺结节的形状与粗血管、厚支气管壁有着比较明显的差异，本发明方法基于三维局部空间的Hessian矩阵，构建三种形状选择性增强滤波器，分别用来初步甄别CT图像中的疑似肺结节、疑似粗血管及疑似厚支气管壁。设f(x，y，z)为三维CT图像中像素(x，y，z)处的CT值，提取三维感兴趣区域的具体实现过程如下：

(2.1)确定形状选择性增强滤波器的尺度范围；

虽然结节的准确尺度在分割之前很难估计，但根据结节的定义，其直径是在[d₀＝3mm，d₁＝30mm]之间。根据肺结节大小的范围[d₀，d₁]和候选尺度个数N₁，N₁由用户预先设定，其取值范围可以为[3，7]，同时根据公式(1)确定形状选择性增强滤波器的尺度：

${\mathrm{\σ}}_1=d_0/4,{\mathrm{\σ}}_2=r\×{\mathrm{\σ}}_1,{\mathrm{\σ}}_3=r^2\×{\mathrm{\σ}}_1,......,{\mathrm{\σ}}_{N_1}=r^{N_1-1}\×{\mathrm{\σ}}_1=d_1/4---\left(1\right)$

其中，σ_i表示第i(1≤i≤N₁)个形状选择性增强滤波器的尺度， $r={(d_1/d_0)}^{1/(N_1-1)}.$

(2.2)对每个像素(x，y，z)，使用步骤(2.1)得到的形状选择性增强滤波器的各尺度分别计算球形、管状和面状选择性增强滤波输出值z_nodule(x，y，z)、z_vessel(x，y，z)和z_airway(x，y，z)，形成三维多尺度形状选择性增强滤波图像；

分别应用尺度为σ_i(1≤i≤N₁)的三维高斯球形滤波器对三维CT图像进行卷积得到卷积图像，即得到像素(x，y，z)处卷积后的CT值

并计算其每个像素(x，y，z)的9个CT值差分

和

$f_{\mathrm{xx},{\mathrm{\σ}}_i}=f_{{\mathrm{\σ}}_i}(x+1,y,z)+f_{{\mathrm{\σ}}_i}(x-1,y,z)-2\×f_{{\mathrm{\σ}}_i}(x,y,z)$

$f_{\mathrm{yy},{\mathrm{\σ}}_i}=f_{{\mathrm{\σ}}_i}(x,y+1,z)+f_{{\mathrm{\σ}}_i}(x,y-1,z)-2\×f_{{\mathrm{\σ}}_i}(x,y,z)$

$f_{\mathrm{zz},{\mathrm{\σ}}_i}=f_{{\mathrm{\σ}}_i}(x,y,z+1)+f_{{\mathrm{\σ}}_i}(x,y,z-1)-2\×f_{{\mathrm{\σ}}_i}(x,y,z)$

$f_{\mathrm{xy},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x+1,y+1,z)+f_{{\mathrm{\σ}}_i}(x-1,y-1,z)-f_{{\mathrm{\σ}}_i}(x+1,y-1,z)-f_{{\mathrm{\σ}}_i}(x-1,y+1,z))/4$

$f_{\mathrm{xz},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x+1,y,z+1)+f_{{\mathrm{\σ}}_i}(x-1,y,z-1)-f_{{\mathrm{\σ}}_i}(x+1,y,z-1)-f_{{\mathrm{\σ}}_i}(x-1,y,z+1))/4---\left(2\right)$

$f_{\mathrm{yx},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x+1,y+1,z)+f_{{\mathrm{\σ}}_i}(x-1,y-1,z)-f_{{\mathrm{\σ}}_i}\left(x-1,y+1,z\right)-f_{{\mathrm{\σ}}_i}(x+1,y-1,z))/4$

$f_{\mathrm{yz},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x,y+1,z+1)+f_{{\mathrm{\σ}}_i}(x,y-1,z-1)-f_{{\mathrm{\σ}}_i}(x,y+1,z-1)-f_{{\mathrm{\σ}}_i}(x,y-1,z+1))/4$

$f_{\mathrm{zx},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x+1,y,z+1)+f_{{\mathrm{\σ}}_i}(x-1,y,z-1)-f_{{\mathrm{\σ}}_i}(x-1,y,z+1)-f_{{\mathrm{\σ}}_i}(x+1,y,z-1))/4$

$f_{\mathrm{zy},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x,y+1,z+1)+f_{{\mathrm{\σ}}_i}(x,y-1,z-1)-f_{{\mathrm{\σ}}_i}(x,y-1,z+1)-f_{{\mathrm{\σ}}_i}(x,y+1,z-1))/4$

将这9个CT值差分代入如公式(3)所示的Hessian矩阵，计算其三个特征值

和

$H_{{\mathrm{\σ}}_i}=\left[\begin{array}{ccc}{f}_{\mathrm{xx},{\mathrm{\σ}}_i}& f_{\mathrm{xy},{\mathrm{\σ}}_i}& f_{\mathrm{xz},{\mathrm{\σ}}_i}\\ f_{\mathrm{yx},{\mathrm{\σ}}_i}& f_{\mathrm{yy},{\mathrm{\σ}}_i}& f_{\mathrm{yz},{\mathrm{\σ}}_i}\\ f_{\mathrm{zx},{\mathrm{\σ}}_i}& f_{\mathrm{zy},{\mathrm{\σ}}_i}& f_{\mathrm{zz},{\mathrm{\σ}}_i}\end{array}\right]---\left(3\right)$

据此构造出尺度为σ_i的球形、管状和面状选择性增强滤波器如公式(4)～(6)所示：

其中，和

分别表示像素(x，y，z)在尺度σ_i下的球形、管状和面状选择性增强滤波器的输出函数值。

取所有不同尺度σ_i(1≤i≤N₁)下

和

中的最大值，分别用z_nodule(x，y，z)、z_vessel(x，y，z)和z_airway(x，y，z)表示，即

$z_{\mathrm{vessel}}(x,y,z)=\underset{i}{\max }[z_{\mathrm{vessel},{\mathrm{\σ}}_i}(x,y,z)\×{\mathrm{\σ}}_i^2],$ $z_{\mathrm{airway}}(x,y,z)=\underset{i}{\max }[z_{\mathrm{airway},{\mathrm{\σ}}_i}(x,y,z)\×{\mathrm{\σ}}_i^2],$ 1≤i≤N₁，从而形成三个三维多尺度形状选择性增强滤波图像，分别为三维多尺度球形选择性增强滤波图像、三维多尺度管状选择性增强滤波图像和三维多尺度面状选择性增强滤波图像，分别用来初步甄别CT图像中的疑似肺结节、疑似粗血管和疑似厚支气管壁。

(2.3)在(2.2)得到的三维多尺度球形选择性增强滤波图像中，去除假阳性结节，得到候选结节区域。

由于(2.2)得到的三维多尺度球形选择性增强滤波图像既包括真阳性结节，也包括假阳性结节，其中假阳性结节大多位于粗血管、厚支气管壁或者其交叉区域。因此本发明方法采用CT值作为分割的主要依据来分割候选结节区域，再从中去除假阳性结节，其具体实现方法如下：

①选取一个滤波阈值T₂，T₂越低，真结节会被检出的越多，即结节的灵敏度会提高，但同时假阳性个数也会随之上升；反之，T₂越高，真结节会被检出得越少，即结节的灵敏度会降低，但同时假阳性个数也会随之减少，理论上T₂可以为任何大于0的实数，但取[10，50]区间内的值比较合理。将三维多尺度球形选择性增强滤波图像上所有滤波输出值大于T₂的连通区域均作为种子区域，记种子区域的个数为N₂，第j个种子区域记为G_j，1≤j≤N₂，设生长次数t的初始值为1；

②计算三维CT图像中对应于种子区域G_j内部的CT值均值

和方差|G_j|为G_j包含的像素数；

③将满足下述生长条件的像素(x，y，z)加入到种子区域G_j内部，生长条件为：a)像素(x，y，z)与G_j是26邻接的；b)f(x，y，z)∈|m_j-3×v_j，m_j+3×v_j]；

④为了避免与血管粘连的结节生长进入血管区域内，选取一个合适的生长次数阈值T₃，(T₃原则上可以取任何正整数，但取值范围为[3，7]时效果更好)，判断生长次数t，若t＜T₃，则t＝t+1，转步骤②；否则，结束生长，将此时的种子区域G_j作为候选结节区域Z_j；

⑤判断j，若j＜N₂，则令j＝j+1，t＝1，然后转入步骤②；否则，利用Fisher线性分类器对所有候选结节区域Z_j特征值进行分类，1≤j≤N₂，去除分类结果中的假阳性结节。

Fisher线性分类器是通过降维的思想来寻找最优线性判别函数，将两类特征向量(特征空间样本)进行分类。本发明对候选结节区域Z_j(1≤j≤N₂)，从三维CT图像、三维多尺度球形选择性增强滤波图像、三维多尺度管状选择性增强滤波图像以及三维多尺度面状选择性增强滤波图像中提取20个特征，具体描述如图3所示。

每个候选结节区域对应的特征形成一个N₃维特征向量(N₃为提取的特征个数)，将N₃维空间的特征向量投影到某个方向的直线上，使用Fisher评价函数来描述两类特征向量投影的分离程度，使特征向量的投影能最好地分开。Fisher评价函数是指特征向量投影的平均值间距与各类内的方差的比值，当Fisher评价函数最大时，两类特征向量在该特征空间中能够被最好地分开。

以保留下来的候选结节区域的质心为中心的83×83×83(默认设置)立方体就是所需的三维感兴趣区域。

(3)在步骤(2)提取的三维感兴趣区域内使用基于三维空间的水平集方法分割出肺结节，确定肺结节边界曲面，本发明具体实现过程如图4所示：

(3.1)令迭代次数k的初始值为0，将步骤(2)提取的三维感兴趣区域的闭合边界作为初始边界曲面C^(k)，并根据距边界的有符号离散距离函数得到初始水平集函数g^(k)；

距边界的有符号离散距离函数是指对于给定的闭合边界曲面，空间上某点的函数值大小为该点到给定边界曲面的距离，其符号根据该点在闭合曲面内部还是外部决定(本发明中定义闭合曲面内部的点对应的函数值符号为正，外部为负)。

(3.2)计算三维CT图像上边界曲面C^(k)内部区域与外部区域的CT值均值c₁和c₂，并得到下一代水平集函数g^(k+1)，具体实现步骤如下：

①本实施例使用公式(7)计算边界曲面C^(k)内部区域Ω₁的CT值均值c₁：

$c_1\left(C^{\left(k\right)}\right)=\frac{\underset{(x,y,z)\∈{\mathrm{\Ω}}_1}{\mathrm{\Σ}}f(x,y,z)}{\mathrm{Volume}\left({\mathrm{\Ω}}_1\right)}---\left(7\right)$

其中，Volume(Ω₁)为边界曲面C^(k)内部区域Ω₁的体积。

②计算边界曲面C^(k)外部区域的CT值均值c₂；

首先判断边界曲面C^(k)所围区域是高对比度结节区域还是低对比度结节区域。本发明使用的判别方法为：选取一个CT值阈值T₄，取值范围为[-700HU，300HU]，对于给定的结节区域，取其中心点及其周围26邻域CT值均值，记为m_CT，若m_CT小于T₄，则判为低对比度结节区域，否则判为高对比度结节区域。然后使用公式(8)计算c₂：

其中，Ω₂＝{(x，y，z)|(x，y，z)∈肺内区域曲面C^(k)内部区域}，Ω₃＝{(x，y，z)|(x，y，z)∈Ω₂∩f(x，y，z)＜c₁}，即Ω₂是所有属于肺内区域中曲面C^(k)外部区域的点组成的区域，Ω₃是所有属于Ω₂区域且CT值小于c₁的点组成的区域。Volume(Ω₂)和Volume(Ω₃)分别表示区域Ω₂和Ω₃的体积，LowVal为用户预先设定的一个小于肺实质CT值的常量参数，取值范围可以为[-1200HU，-900HU]，V_ext表示低CT值像素的体积V_ext＝((1-p)/p)×Volume(Ω₂)，p定义为max(min(β₁×Volume(Ω₁)+β₂，1)，0.5)(该值的取值范围为[0.5，1])，β₁与β₂均为预先设定的常量参数，其取值范围可以为[0，1]。

③使用公式(9)得到下一代水平集函数g^(k+1)：

$g^{(k+1)}={\mathrm{\δ}}_{\mathrm{\ϵ}}\left(g^{\left(k\right)}\right)[\mathrm{\μ}\mathrm{div}\left(\frac{\▿g^{\left(k\right)}}{|\▿g^{\left(k\right)}|}\right)-v-{\mathrm{\λ}}_1{(f-c_1)}^2+{\mathrm{\λ}}_2{(f-c_2)}^2]---\left(9\right)$

其中，预先设定的常量参数μ值越大，图像的分割结果越平滑，但实际操作时，由于数据离散化的影响，μ值设置过大会导致图像演化的异常，因此μ的取值范围为[5，25]。预先设定的常量参数v作为分割区域面积项的系数，控制曲线向内收缩，由于曲线所围面积与其长度有一定的相关性，因此v被设置为0。预先设定的常量参数λ₁和λ₂的作用主要有表现为以下两方面：一，控制分割结果的圆形度，μ过大会导致分割异常，因此只能通过同时降低λ₁和λ₂的值来等价的提高曲率力在演化过程中的权重，从而达到平滑曲面的作用；二，使分界面偏向某一聚类中心，通过调节

大小来达到此目的。因此，λ₁和λ₂的取值范围均可以为(1.0/(2200×2200)，1)。f就是指f(x，y，z)，

表示向量

的散度，

为水平集函数g^(k)梯度，为水平集函数g^(k)梯度的范数，δ_ε(g^(k))用于限制水平集函数的演化范围，定义为如公式(10)所示的K_ε(g^(k))函数的导数：

$K_{\mathrm{\ϵ}}\left(g^{\left(k\right)}\right)=\frac{1}{2}(1+\frac{2}{\mathrm{\π}}\arctan \left(\frac{g^{\left(k\right)}}{\mathrm{\ϵ}}\right))---\left(10\right)$

其中ε为预先设定的常量参数，取值范围为[-1，1]。

(3.3)由g^(k+1)的零水平集得到新的边界曲面C^(k+1)；

(3.4)比较C^(k)与C^(k+1)，判断|C^(k+1)-C^(k)|是否收敛，如果|C^(k+1)-C^(k)|小于预先设定的阈值T₅，则收敛，C^(k+1)即为所求的肺结节边界曲面；否则，令k＝k+1，转步骤(3.2)。

阈值T₅的取值范围可以设为[0，10]。

(4)根据步骤(3)得到的肺结节边界曲面信息，结合感兴趣区域的总体特性，提取相关特征并量化特征值；

本实施例结合了肺结节边界曲面信息和感兴趣区域的总体特性，提取了包括最大直径、体积、球形度、边界清晰度、最大密度、最小密度、平均密度、密度方差、钙化含量、脂肪含量、毛刺征、分叶征、胸膜凹陷征在内的13个特征并进行了量化，其中最大直径、体积、球形度、边界清晰度、最大密度、最小密度、平均密度、密度方差、钙化含量、脂肪含量是在肺结节三维分割基础上提取，毛刺征、分叶征、胸膜凹陷征的提取结合了肺结节和感兴趣区域在二维和三维空间中的特性。

本发明方法对上述各特征的量化方法并没有特别的规定，下面对每一种特征列举一种方法进行特征量化，但也可以采用其它的特征量化方法对各特征进行量化。

设分割后第q个肺结节的边界曲面为E_q，其内部区域为S_q，具体特征量化方法如下：

(4.1)最大直径

最大直径是描述肺结节大小的一个特征，本发明方法使用公式(11)定义的肺结节边界曲面中两点间的最大距离，表示第q个肺结节的最大直径：

(4.2)毛刺征

毛刺为近似垂直于肺结节边界生长的“蟹足式”针状结构。毛刺与周围肺实质的对比度比肺结节边界与肺实质的对比度低得多，且形状狭窄，因此很难将其分割出来。本发明根据毛刺生长的特点，直接在三维感兴趣区域的中心层图像上，通过引入边界法线方向作为参照，分析边界附近区域梯度方向的规律来量化毛刺水平。具体方法如下：

①如图5所示，将三维感兴趣区域中心层所对应的结节区域的边界线标记为edge，从边界edge开始向外宽P_最大直径(q)/4的带状区域为第1层，标记为A₁；而从第1层的外边界开始向外宽P_最大直径(q)/4的带状区域则为第2层，标记为A₂；

②使用公式(12)计算edge上所有边界点(x，y)的梯度-法线正交指数GN(x，y)：

GN(x，y)＝ρ₁×GN₁(x，y)+ρ₂×GN₂(x，y)       (12)

其中，GN₁(x，y)和GN₂(x，y)分别表示A₁和A₂区域的梯度-法线正交指数，其定义分别如公式(13)和(14)所示，ρ₁和ρ₂分别为GN₁(x，y)和GN₂(x，y)的权值，由于毛刺大部分集中在第1层区域，而离边界较远的第2层区域中干扰信息也更多，因此ρ₁＞ρ₂≥0。

W(x，y)是以edge上的边界点(x，y)为顶点，并覆盖该点法线向量的边界正方形感兴趣区域，其边长为P_最大直径(q)/2；点(x′，y′)是指W(x，y)与A₁相交的区域中CT值在[-850HU，300HU]范围内的点，点(x″，y″)是指W(x，y)与A₂相交区域中CT值在[-850HU，300HU]范围内的点；

和

分别为点(x′，y′)和点(x″，y″)的梯度方向与点(x，y)的法线方向的夹角；当夹角越接近π/2，其

值越大。|{(x′，y′)|(x′，y′)∈W(x，y)∩A₁，f(x′，y′)∈[-850HU，300HU]}|表示W(x，y)与A₁相交的区域中CT值在[-850HU，300HU]范围内的点的个数，|{(x″，y″)|(x″，y″)∈W(x，y)∩A₂，f(x″，y″)∈[-850HU，300HU]}|表示W(x，y)与A₂相交的区域中CT值在[-850HU，300HU]范围内的点的个数。

③由于毛刺分散生长于边界线上，并不是所有边界线edge上的点都具有大的梯度-法线正交指数，因此本发明方法仅根据其中具有较大的梯度-法线正交指数的点来量化第q个肺结节的毛刺征P_毛刺征(q)如公式(15)所示：

其中，边界曲面E_q中所有边界点的梯度-法线正交指数GN(x，y)由大到小排序后，排在前5％～10％的那些边界点组合的集合表示为E_max，且|E_max|表示E_max中包含的像素数。

P_毛刺征(q)值越大，说明该结节的毛刺征越明显。本发明根据P_毛刺征(q)取值，将肺结节的毛刺征水平分为三类：如果P_毛刺征(q)∈[0，25)，说明该结节无毛刺；如果P_毛刺征(q)∈[25，40)，说明该结节有短毛刺；如果P_毛刺征(q)∈[40，80]，说明该结节有长毛刺。

(4.3)体积

体积是肺结节量化的基本特征，常用于估计肺结节的倍增时间，本发明方法使用公式(16)计算第q个肺结节的体积：

P_体积(q)＝|S_q|×(PS_a×PS_c×PS_s)            (16)

其中，|S_q|表示第q个肺结节内部区域S_q包含的像素数；PS_a、PS_c及PS_s分别为CT图像中横断面、冠状面及矢状面中的像素间距，在等方性处理后PS_a＝PS_c＝PS_s。

(4.4)球形度

球形度是描述肺结节形状近似球形的程度，在一定程度上描述了结节的规则程度，它是0到1的数，当球形度越接近1，则说明该肺结节越类似于球形。本发明方法使用公式(17)定义第q个肺结节的球形度：

其中，D_q为第q个肺结节的表面积。

(4.5)边界清晰度

边界清晰度是描述肺结节边界清晰/模糊程度的参数，本发明方法使用公式(18)定义第q个肺结节的边界清晰度：

其中，|grad(x，y，z)|为结节边界曲面E_q上像素(x，y，z)的梯度grad(x，y，z)的绝对值。P_{边界清晰度}(q)值越大，说明该肺结节的边界越清晰。

(4.6)分叶征

分叶征是由于肿瘤生长异向性产生的表面凹陷征，分叶征的凹陷可能发生在两个相邻的凸点之间。

以三维肺结节分割结果中的第layer层为例，从结节区域的中心点向每个边界点按逆时针引线段，如果某个边界点到结节区域的中心点的线段长度均大于与其直接相邻的两个边界点到结节区域的中心点的线段长度，则该边界点为凸点；如果某个边界点到结节区域的中心点的线段长度均小于与其直接相邻的两个边界点到结节区域的中心点的线段长度，则该边界点为凹点，因此可将边界点划分为凸点与凹点两个集合，分别记为fp_layer与sp_layer。本发明方法使用公式(19)来量化第q个肺结节的分叶征P_分叶征(q)：

其中，|sp_layer|表示该肺结节第layer层凹点的个数，|fp_layer|表示该肺结节第layer层凸点的个数，该肺结节第layer层的第n(1≤n≤|fp_layer|-1)对相邻凸点的分叶指数LB_q，layer，n定义如公式(20)所示：

${\mathrm{LB}}_{q,\mathrm{layer},n}=\min (\frac{\mathrm{CDis}}{\mathrm{CLen}},1)---\left(20\right)$

其中，CLen为这两个相邻凸点间线段的长度，CDis为该结吊区域的中心点到这两个凸点连线间的距离与到凸点间每个边界点的长度之差的最大值。

P_分叶征(q)值越大，说明该肺结节的分叶程度越明显。本发明根据P_分叶征(q)取值对肺结节的分叶程度进行分类：如果P_分叶征(q)∈[0，0.1]，说明该肺结节没有分叶；如果P_分叶征(q)∈(0.1，0.2]，说明该肺结节有浅分叶；如果P_分叶征(q)∈(0.2，0.4]，说明该肺结节有中分叶；如果P_分叶征(q)∈(0.4，1]，说明该肺结节有深分叶。

(4.7)胸膜凹陷征

典型的胸膜凹陷征是近脏层胸膜受到附近肺结节的牵引，其表面呈现小三角形影或小喇叭状阴影，三角形的底部在胸壁，尖指向结节，结节与三角形影之间可为线状影相连。

本发明方法从肺边界的某点开始按照逆时针顺序搜索第一个可能位于胸膜凹陷三角形底部的肺边界点，该边界点满足两个条件：a)它与某个结节区域中心的距离在[P_最大直径(q)/2，2×P_最大直径(q)]区间；b)其CT值高于生长阈值T₆，取值范围是[-800HU，-900HU]。以该边界点为种子点进行区域生长，所得的连通区域所包含的像素数目在范围[10，100]内，且其区域的形状为以肺边界线为底的类三角形，则该连通区域发生了胸膜凹陷；反之，该连通区域不是胸膜凹陷区域，则继续搜索下一个可能的边界点，直到搜索完整个肺边界。

(4.8)最大密度

本发明方法使用公式(21)定义第q个肺结节的最大密度：

(4.9)最小密度

本发明方法使用公式(22)定义第q个肺结节的最小密度：

(4.10)平均密度

本发明方法使用公式(23)定义第q个肺结节的平均密度：

(4.11)密度方差

根据密度方差大小，肺结节可分为均匀结节和非均匀结节，本发明方法使用公式(24)定义第q个肺结节的密度方差：

(4.12)钙化含量

钙化含量是指在肺结节内部钙化组织占肺结节总体积的百分比，本发明方法使用公式(25)定义第q个肺结节的钙化含量：

其中，|{(x，y，z)|(x，y，z)∈S_q∩f(x，y，z)∈[100HU，3000HU]}|表示第q个肺结节内部区域S_q中CT值在[100HU，3000HU]范围内的像素数。肺结节内部钙化组织越多，P_钙化含量(q)取值越大；反之，越小。

(4.13)脂肪含量

脂肪含量是指在肺结节内部脂肪组织占肺结节总体积的百分比，本发明方法使用公式(26)定义第q个肺结节的脂肪含量：

其中，|{(x，y，z)|(x，y，z)∈S_q∩f(x，y，z)∈[-90HU，-50HU]}|表示第q个肺结节内部区域S_q中CT值在[-90HU，-50HU]范围内的像素数。肺结节内部脂肪组织越多，P_脂肪含量(q)取值越大；反之，越小。

(5)根据步骤(3)得到的肺结节边界曲面信息，将肺结节的位置和范围以三维形式显示，同时以用户需要的形式存储及显示步骤(4)得到的特征及其量化值。

如图2所示，本发明提供的一种肺结节三维分割与特征提取系统，包括肺区域分割模块1、三维感兴趣区域提取模块2、肺结节三维分割模块3、肺结节特征提取及量化模块4和输出模块5。

肺区域分割模块1，按照上述步骤(1)所述的方法，将输入的所有二维胸部CT横断面图像中肺实质分割出来后，使用几何凸点搭桥法修补肺表面的凹陷，并将分割出的肺内区域传送给三维感兴趣区域提取模块2；

三维感兴趣区域提取模块2，按照上述步骤(2)所述的步骤从肺内区域中提取候选结节区域，利用Fisher线性分类器对候选结节区域的特征值进行分类从而获取三维感兴趣区域，并将感兴趣区域的边界信息传送给肺结节三维分割模块3和结节特征提取及量化模块5；

肺结节三维分割模块3按照上述步骤(3)所述的过程，使用基于三维空间的水平集方法从三维感兴趣区域中将肺结节分割出来，并将得到的肺结节边界曲面信息传送给结节特征提取及量化模块4和输出模块5；

结节特征提取及量化模块4用于接收肺结节三维分割模块3获得的肺结节的边界曲面信息，结合三维感兴趣区域提取模块2得到的三维感兴趣区域的总体特性，提取肺结节的相关特征，计算特征值，并将其特征及其量化值传送给输出模块5；

输出模块5用于接收肺结节三维分割模块3得到的肺结节边界曲面信息，并将其位置和范围以三维形式显示，同时接收结节特征提取及量化模块4得到的特征及其量化值，将其结果根据需要存储及显示。

实施例：

本发明提出的肺结节三维分割与特征提取系统中涉及到若干参数，这些参数要针对具体处理的数据特点进行综合调节设定以达到整体系统的良好性能，此处列出针对本发明处理数据集合而设定的参数：

步骤(2.1)在确定形状选择性增强滤波器的尺度范围和候选尺度个数时，待增强物体直径范围为[d₀＝3mm，d₁＝30mm]，候选尺度个数N₁取5。

步骤(2.3)选择三维多尺度球形选择性增强滤波图像上的种子区域时，选取的滤波阈值T₂为30。同时，生长次数阈值T₃取值为5。

步骤(3.2)判断边界曲面所围区域是高对比度结节区域还是低对比度结节区域时，选取的CT阈值T₄＝-430HU；常量参数LowVal取值为-1200HU，常量参数β₁与β₂取值分别为0.00019075和0.39849711；求解下一代水平集函数g^(k+1)时，常量参数ε，μ，v，λ₁和λ₂的取值分别为0.0001，10，0，1.0/2200/1.5和1.0/2200/1.5。

步骤(3.4)判断|C^(k+1)-C^(k)|是否收敛时，阈值T₅的取值为5。

步骤(4.2)A₁和A₂两个区域的梯度-法线正交指数的权值ρ₁和ρ₂取值分别为1.5和1。

步骤(4.7)判断胸膜凹陷征时，生长阈值T₆取值为-850HU。

本发明方法通过一种肺结节三维分割与特征提取系统自动地对胸部二维CT横断面图像中肺内三维感兴趣区域进行分析和处理，提供肺结节三维边界曲面和CT值，并根据需要提供肺结节相关的一系列特征参数，此方法利用了三维局部空间Hessian矩阵的多尺度形状选择性增强滤波器，可以减少粗血管、厚支气管壁或者其分叉区域对候选结节区域的干扰，能降低假阳性结节出现的概率。在特征提取时，结合了感兴趣区域和结节边界曲面信息，可以减少其周围正常组织对结节特征的干扰，提高了特征提取的准确性，还可以提取结节的体积、边界清晰度、钙化含量等特征。本发明的实现并不局限于上述实例所公开的范围，可以采用不同于上述实例的方式实现上述技术方案。

Claims (9)

1.一种肺结节三维分割与特征提取方法，其特征在于，该方法包括下述步骤：

第1步在输入的所有二维胸部CT横断面图像上自动分割出肺实质后，使用几何凸点搭桥法修补肺表面的凹陷，将肺结节所在的肺内区域与肺外区域分隔开；

第2步构建三维空间中不同尺度的三种形状选择性增强滤波器，分别仅对球形、管状及面状物体具有响应；利用这三种形状选择性增强滤波器分别初步甄别出由所有二维胸部CT横断面图像构建的三维CT图像中的疑似肺结节、疑似粗血管及疑似厚支气管壁，提取出肺内区域中候选结节区域，再利用Fisher线性分类器对候选结节区域的特征值进行分类从而得到三维感兴趣区域；

第3步在第2步提取的三维感兴趣区域内使用基于三维空间的水平集方法分割出肺结节，确定肺结节边界曲面；

第4步根据第3步得到的肺结节边界曲面的信息，结合感兴趣区域的总体特性，提取相关特征并量化特征值；

第5步根据第3步得到的肺结节边界曲面信息，将肺结节的位置和范围以三维形式显示，同时以用户需要的形式存储及显示第4步得到的特征及其量化值。

2.根据权利要求1所述的肺结节三维分割与特征提取方法，其特征在于：第2步中，设f(x，y，z)为三维CT图像中像素(x，y，z)处的CT值，提取三维感兴趣区域的具体实现过程如下：

(2.1)确定每一种形状选择性增强滤波器的尺度：

定义N₁为用户预先设定的候选尺度的个数，肺结节大小的范围为[d₀，d₁]，利用公式(I)计算形状选择性增强滤波器的第i个尺度σ_i，1≤i≤N₁，则，

${\mathrm{\σ}}_1=d_0/4,{\mathrm{\σ}}_2=r\×{\mathrm{\σ}}_1,{\mathrm{\σ}}_3=r^2\×{\mathrm{\σ}}_1,......,{\mathrm{\σ}}_{N_1}=r^{N_1-1}\×{\mathrm{\σ}}_1=d_1/4---\left(I\right)$

其中， $r={(d_1/d_0)}^{1/(N_1-1)};$

(2.2)对每个像素(x，y，z)，使用步骤(2.1)得到的形状选择性增强滤波器的各尺度分别计算球形、管状和面状选择性增强滤波输出值z_nodule(x，y，z)、z_vessel(x，y，z)和z_airway(x，y，z)，形成三维多尺度形状选择性增强滤波图像；

(2.3)在步骤(2.2)得到的三维多尺度球形选择性增强滤波图像中，去除假阳性结节，得到候选结节区域。

3.根据权利要求2所述的肺结节三维分割与特征提取方法，其特征在于：步骤(2.2)的具体过程如下：

分别应用尺度为σ_i的三维高斯球形滤波器对三维CT图像进行卷积得到卷积图像，即得到像素(x，y，z)处卷积后的CT值

并计算其每个像素(x，y，z)的9个CT值差分

和

$f_{\mathrm{xx},{\mathrm{\σ}}_i}=f_{{\mathrm{\σ}}_i}(x+1,y,z)+f_{{\mathrm{\σ}}_i}(x-1,y,z)-2\×f_{{\mathrm{\σ}}_i}(x,y,z)$

$f_{\mathrm{yy},{\mathrm{\σ}}_i}=f_{{\mathrm{\σ}}_i}(x,y+1,z)+f_{{\mathrm{\σ}}_i}(x,y-1,z)-2\×f_{{\mathrm{\σ}}_i}(x,y,z)$

$f_{\mathrm{zz},{\mathrm{\σ}}_i}=f_{{\mathrm{\σ}}_i}(x,y,z+1)+f_{{\mathrm{\σ}}_i}(x,y,z-1)-2\×f_{{\mathrm{\σ}}_i}(x,y,z)$

$f_{\mathrm{xy},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x+1,y+1,z)+f_{{\mathrm{\σ}}_i}(x-1,y-1,z)-f_{{\mathrm{\σ}}_i}(x+1,y-1,z)-f_{{\mathrm{\σ}}_i}(x-1,y+1,z))/4$

$f_{\mathrm{xz},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x+1,y,z+1)+f_{{\mathrm{\σ}}_i}(x-1,y,z-1)-f_{{\mathrm{\σ}}_i}(x+1,y,z-1)-f_{{\mathrm{\σ}}_i}(x-1,y,z+1))/4$

$f_{\mathrm{yx},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x+1,y+1,z)+f_{{\mathrm{\σ}}_i}(x-1,y-1,z)-f_{{\mathrm{\σ}}_i}(x-1,y+1,z)-f_{{\mathrm{\σ}}_i}(x+1,y-1,z))/4$

$f_{\mathrm{yz},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x,y+1,z+1)+f_{{\mathrm{\σ}}_i}(x,y-1,z-1)-f_{{\mathrm{\σ}}_i}(x,y+1,z-1)-f_{{\mathrm{\σ}}_i}(x,y-1,z+1))/4$

$f_{\mathrm{zx},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x+1,y,z+1)+f_{{\mathrm{\σ}}_i}(x-1,y,z-1)-f_{{\mathrm{\σ}}_i}(x-1,y,z+1)-f_{{\mathrm{\σ}}_i}(x+1,y,z-1))/4$

$f_{\mathrm{zy},{\mathrm{\σ}}_i}=(f_{{\mathrm{\σ}}_i}(x,y+1,z+1)+f_{{\mathrm{\σ}}_i}(x,y-1,z-1)-f_{{\mathrm{\σ}}_i}(x,y-1,z+1)-f_{{\mathrm{\σ}}_i}(x,y+1,z-1))/4$

将这9个CT值差分代入如公式(II)所示的Hessian矩阵，计算其三个特征值

和

$H_{{\mathrm{\σ}}_i}=\left[\begin{array}{ccc}{f}_{\mathrm{xx},{\mathrm{\σ}}_i}& f_{\mathrm{xy},{\mathrm{\σ}}_i}& f_{\mathrm{xz},{\mathrm{\σ}}_i}\\ f_{\mathrm{yx},{\mathrm{\σ}}_i}& f_{\mathrm{yy},{\mathrm{\σ}}_i}& f_{\mathrm{yz},{\mathrm{\σ}}_i}\\ f_{\mathrm{zx},{\mathrm{\σ}}_i}& f_{\mathrm{zy},{\mathrm{\σ}}_i}& f_{\mathrm{zz},{\mathrm{\σ}}_i}\end{array}\right]---\left(\mathrm{II}\right)$

据此构造出尺度为σ_i的球形、管状和面状选择性增强滤波器如公式(III)～(V)所示：

其中，

和

分别表示像素(x，y，z)在尺度σ_i下的球形、管状和面状选择性增强滤波器的输出函数值；

取所有不同尺度σ_i下

和

中的最大值，1≤i≤N₁，分别用z_nodule(x，y，z)、z_vessel(x，y，z)和z_airway(x，y，z)表示，形成三维多尺度球形选择性增强滤波图像、三维多尺度管状选择性增强滤波图像和三维多尺度面状选择性增强滤波图像。

4.根据权利要求2所述的肺结节三维分割与特征提取方法，其特征在于：步骤(2.3)中，去除假阳性结节的具体过程为：

①选取一个滤波阈值T₂，将三维多尺度球形选择性增强滤波图像上所有滤波输出值大于T₂的连通区域均作为种子区域，记种子区域的个数为N₂，第j个种子区域记为G_j，1≤j≤N₂，设生长次数t的初始值为1；

②计算三维CT图像中对应于种子区域G_j内部的CT值均值和方差

|G_j|为G_j包含的像素数；

③将满足下述生长条件的像素(x，y，z)加入到种子区域G_j内部，生长条件为：a)像素(x，y，z)与G_j是26邻接的；b)f(x，y，z)∈[m_j-3×v_j，m_j+3×v_j]；

④选取一个能够避免与血管粘连的结节生长进入血管区域内的生长次数阈值T₃，判断生长次数t，若t＜T₃，则t＝t+1，转步骤②；否则，结束生长，将此时的种子区域G_j作为候选结节区域Z_j；

⑤判断j，若j＜N₂，则令j＝j+1，t＝1，然后转入步骤②；否则，利用Fisher线性分类器对所有候选结节区域Z_j特征值进行分类，1≤j≤N₂，去除分类结果中的假阳性结节。

5.根据权利要求1所述的肺结节三维分割与特征提取方法，其特征在于：第3步包括下述具体过程：

(3.1)令迭代次数k的初始值为0，将第2步提取的三维感兴趣区域的闭合边界作为初始边界曲面C^(k)，并根据距边界的有符号离散距离函数得到初始水平集函数g^(k)；

距边界的有符号离散距离函数是指对于给定的闭合边界曲面，空间上某点的函数值大小为该点到给定边界曲面的距离，其符号根据该点在闭合曲面内部还是外部决定；

(3.2)计算三维CT图像上边界曲面C^(k)内部区域与外部区域的CT值均值c₁和c₂，并得到下一代水平集函数g^(k+1)；

(3.3)由g^(k+1)的零水平集得到新的边界曲面C^(k+1)；

(3.4)比较C^(k)与C^(k+1)，判断|C^(k+1)-C^(k)|是否收敛，若收敛，C^(k+1)即为所求的肺结节边界曲面；否则，令k＝k+1，再转步骤(3.2)。

6.根据权利要求5所述的肺结节三维分割与特征提取方法，其特征在于：步骤(3.2)中，CT值均值c₁、c₂和下一代水平集函数g^(k+1)的计算过程为：

①使用公式(VI)计算边界曲面C^(k)内部区域Ω₁的CT值均值c₁：

$c_1\left(C^{\left(k\right)}\right)=\frac{\underset{(x,y,z)\∈{\mathrm{\Ω}}_1}{\mathrm{\Σ}}f(x,y,z)}{\mathrm{Volume}\left({\mathrm{\Ω}}_1\right)}---\left(\mathrm{VI}\right)$

其中，Volume(Ω₁)为边界曲面C^(k)内部区域Ω₁的体积；

②计算边界曲面C^(k)外部区域的CT值均值c₂；

对于给定的结节区域，取其中心点及其周围26邻域CT值均值，记为m_CT，对于预先设定的CT值阈值T₄，使用公式(VII)计算c₂：

其中，

Ω₃＝{(x，y，z)|(x，y，z)∈Ω₂∩f(x，y，z)＜c₁}，Volume(Ω₂)和Volume(Ω₃)分别表示区域Ω₂和Ω₃的体积，LowVal为预先设定的一个小于肺实质CT值的常量参数，低CT值像素的体积V_ext＝((1-p)/p)×Volume(Ω₂)，p定义为max(min(β₁×Volume(Ω₁)+β₂，1)，0.5)，β₁与β₂均为预先设定的常量参数；

③使用公式(VIII)得到下一代水平集函数g^(k+1)：

$g^{(k+1)}={\mathrm{\δ}}_{\mathrm{\ϵ}}\left(g^{\left(k\right)}\right)[\mathrm{\μ}\mathrm{div}\left(\frac{\▿g^{\left(k\right)}}{|\▿g^{\left(k\right)}|}\right)-v-{\mathrm{\λ}}_1{(f-c_1)}^2+{\mathrm{\λ}}_2{(f-c_2)}^2]---\left(\mathrm{VIII}\right)$

其中，μ，v，λ₁和λ₂均为预先设定的常量参数，f就是指f(x，y，z)，

表示向量

的散度，

为水平集函数g(k)的梯度，为

的范数，δ_ε(g^(k))为K_ε(g^(k))函数的导数，K_ε(g^(k))定义如下：

$K_{\mathrm{\ϵ}}\left(g^{\left(k\right)}\right)=\frac{1}{2}(1+\frac{2}{\mathrm{\π}}\arctan \left(\frac{g^{\left(k\right)}}{\mathrm{\ϵ}}\right))$

其中ε为预先设定的常量参数，取值范围为[-1，1]。

7.根据权利要求1所述的肺结节三维分割与特征提取方法，其特征在于：第4步提取的特征包括最大直径、体积、球形度、边界清晰度、最大密度、最小密度、平均密度、密度方差、钙化含量、脂肪含量、毛刺征、分叶征和胸膜凹陷征。

8.根据权利要求7所述的肺结节三维分割与特征提取方法，其特征在于：按照下述过程量化毛刺征：

①将三维感兴趣区域中心层所对应的结节区域的边界线标记为edge，从边界edge开始向外宽P_最大直径(q)/4的带状区域为第1层，标记为A₁；而从第1层的外边界开始向外宽P_最大直径(q)/4的带状区域则为第2层，标记为A₂；其中，P_最大直径(q)表示第q个肺结节的最大直径；

②计算edge上所有边界点(x，y)的梯度-法线正交指数GN(x，y)：

GN(x，y)＝ρ₁×GN₁(x，y)+ρ₂×GN₂(x，y)

其中，GN₁(x，y)和GN₂(x，y)分别表示A₁和A₂区域的梯度-法线正交指数，其定义分别如公式(IX)和(X)所示，ρ₁和ρ₂分别为GN₁(x，y)和GN₂(x，y)的权值，ρ₁＞ρ₂≥0；

W(x，y)是以edge上的边界点(x，y)为顶点，并覆盖该点法线向量的边界正方形感兴趣区域，其边长为P_最大直径(q)/2，点(x′，y′)是指W(x，y)与A₁相交的区域中CT值在[-850HU，300HU]范围内的点，点(x″，y″)是指W(x，y)与A₂相交区域中CT值在[-850HU，300HU]范围内的点；

和

分别为点(x′，y′)和点(x″，y″)的梯度方向与点(x，y)的法线方向的夹角；

③第q个肺结节的毛刺征量化结果P_毛刺征(q)为：

其中，边界曲面E_q中所有边界点的梯度-法线正交指数GN(x，y)由大到小排序后，排在前5％～10％的那些边界点组合的集合表示为E_max，E_q为分割后第q个肺结节的边界曲面，且|E_max|表示E_max中包含的像素数。

9.一种实现权利要求1所述肺结节三维分割与特征提取方法的系统，其特征在于：该系统包括肺区域分割模块(1)、三维感兴趣区域提取模块(2)、肺结节三维分割模块(3)、肺结节特征提取及量化模块(4)和输出模块(5)；

肺区域分割模块(1)用于将输入的所有二维胸部CT横断面图像中肺实质分割出来后，使用几何凸点搭桥法修补肺表面的凹陷，并将分割出的肺内区域传送给三维感兴趣区域提取模块(2)；

三维感兴趣区域提取模块(2)用于从肺内区域中提取候选结节区域，利用Fisher线性分类器对候选结节区域的特征值进行分类从而获取三维感兴趣区域，并将感兴趣区域的边界信息传送给肺结节三维分割模块(3)和结节特征提取及量化模块(5)；

肺结节三维分割模块(3)使用基于三维空间的水平集方法从三维感兴趣区域中将肺结节分割出来，并将得到的肺结节边界曲面信息传送给结节特征提取及量化模块(4)和输出模块(5)；

结节特征提取及量化模块(4)用于接收肺结节三维分割模块(3)获得的肺结节的边界曲面信息，结合三维感兴趣区域提取模块(2)得到的三维感兴趣区域的总体特性，提取肺结节的相关特征，计算特征值，并将其特征及其量化值传送给输出模块(5)；

输出模块(5)用于接收肺结节三维分割模块(3)得到的肺结节边界曲面信息，并将其位置和范围以三维形式显示，同时接收结节特征提取及量化模块(4)得到的特征及其量化值，将其结果根据需要存储及显示。

Images