用于局部肺结构分类以进行计算机辅助结节检测的系统和方法

专利名称：用于局部肺结构分类以进行计算机辅助结节检测的系统和方法
技术领域：
本发明针对在数字化医学图像中局部结构类型的分类。
背景技术：
去年仅在美国就有超过160,000人死于肺癌。虽然预防肺癌最好的方法是不吸烟，但是早期检测是提高患者预后的关键。当在早期检测出癌症并进行手术，非小细胞肺癌I期患者的五年存活率是60％-80％。然而，不动手术的患者面临的五年存活率仅为10％.1。
成像技术(例如计算机断层摄影(CT)扫描)提供了非侵入的、灵敏的早期检测方法。在胸部CT扫描中对肺结节的计算机辅助检测和诊断(CAD)对于更有效和标准化诊断过程降低了人为误差的可能性。在CT扫描中，肺结节呈现出各种形状和尺寸的致密块。它们可能与其它组织例如血管或胸膜分离或粘连。
近来，已经提出了许多技术以在薄层(thin-slice)CT中对结节进行自动检测和分类，包括区域生长和自动阈值确定，与高斯(Gaussian)结节模型的模板匹配，使用3D结节选择和噪声抑制过滤器，结节匹配，以及可变几何和强度模板。然而，CAD系统的这些现状的缺点在于区分结节和其它致密结构(例如血管)。由于在大多数系统中所使用的圆形假设，弯曲的脉管及其交叉点经常被错误地检测为结节，从而导致假阳性(falsepositive)(FP)病例。
为了减少这类FP的数量，以前已经提出了两种解决方案基于相关性的过滤器用以增强模糊形状分析的感兴趣区域以进行脉管树重建，和利用由基于海森(Hessian)的分析所提供的脉管中轴跟踪。前一方法的缺点包括它的不灵活性。在研究中所使用的简单结构模板将不能处理许多复杂的脉管形状和拓扑。另一方面，后一方法尽管能够处理较多不规则的结构但是计算量非常大。

发明内容
这里所描述的本发明典型实施例一般包括用于在胸部CT扫描中对局部结构类型(例如结节、脉管和交叉点)进行分类的方法和系统。该分类对于肺结节的计算机辅助检测(CAD)是有用的，并且可以被用作任何肺部CAD系统的后处理组成部分以便减少由脉管和交叉点所产生的假阳性(FP)。该分类因此假定由这种CAD系统或从放射科医师的报告中提供阳性候选者，从而集中于减少FP的问题。在该方案中，分类结果提供了一种排除由脉管和交叉点引起的假阳性的有效手段从而改善了整体性能。
依据本发明实施例的方法将各种3D拓扑结构的分类转换成更加简单的2D数据聚类分类，对其在文献中可用更普遍和灵活的解决方案，并且其更适用于可视化。除了计算上的好处以外，这种方法具有以下优势，即更普遍和灵活的分析技术清单(inventory)以及更直观的可视化可能性，这在与放射科医师的可能交互作用环境中是有用的。
给定结节候选者，首先稳健地使各向异性高斯拟合于数据。所得到的高斯中心和扩展参数被用来仿射标准化数据域，以便将拟合的各向异性椭圆体弯曲成尺寸固定的各向同性球体。一种自动方法可以提取含有靶结构的适当边界表面的3D球形流形。通过数据驱动熵最小化法来进行尺度选择。采用例如利用自动模数估计的EM聚类、方向统计以及利用修正的Bhattacharyya距离的分级聚类等技术针对对应于突出结构的高强度群(cluster)来分析流形。高强度群的估计数量显性地确定肺结构的类型结节(0)，粘连结节(1)，脉管(2)，交叉点(＞3)。利用方向统计、尤其是多变量包绕高斯建模(multivariate wrappedGaussian modeling)，根据本发明一个实施例的方法扩展了高斯拟合法，包括自动模数选择。
在肺部CAD的范围之外，根据本发明一个实施例的分类方法可以被用于在各种领域(例如血管造影术)中提供血管结构的有意义的信息。这些局部程序比现有技术更灵活和有效，并且将有助于改善通常的肺部CAD系统的精确性。此外，在建模部分中所选择的感兴趣体积(VOI)表示具有有益的可视化能力，例如帮助用户(放射科医师)交互作用的展开的2D边界流形(boundingmanifold)。
在该领域中得出证明有利的分类的对所选择的胸部CT图像例子的定性研究。根据本发明一个实施例的算法可以稳健地对结节、粘连结节、脉管以及脉管交叉点的例子进行分类。
根据本发明的一个方面，提供了一种用于在数字化图像中对肺部结构进行分类的方法，包括在数字化三维(3D)图像中提供一个或多个靶结构的近似的靶结构位置，在每个所述近似的靶位置周围拟合各向异性高斯模型来产生更精确的3D靶模型和所述一个或多个靶结构的中心位置，使每个所述3D靶模型弯曲成3D球体，在每个所述弯曲成的3D球体周围构建边界流形，并且在所述边界流形上识别群，其中对所述一个或多个靶结构进行分类。
根据本发明的另一方面，数字化图像包括与三维栅格上的点域相对应的多个强度。
根据本发明的另一方面，在近似的靶位置周围拟合各向异性高斯模型包括采用高斯尺度-空间均值漂移分析和基于散度的詹森-香农(Jensen-Shannon)自动带宽选择，产生所述靶结构的3D椭圆体模型，其中所述3D椭圆体的中心和尺寸与所述高斯模型的中心和协方差相对应。
根据本发明的另一方面，弯曲所述3D靶模型包括仿射标准化所述3D椭圆体，其中缩放方向和系数从所述各向异性高斯模型的结构协方差获得。
根据本发明的另一方面，构建边界流形进一步包括将弯曲成的球体的3D表面展开成2D表示，并且确定适当的边界流形的半径。
根据本发明的另一方面，将弯曲成的球体的3D表面展开成2D表示包括将所述弯曲成的球体的表面转换成球坐标(θ，)，其中∈[-π，π]和θ∈[-π，π]。
根据本发明的另一方面，确定适当的边界流形的半径包括在所述弯曲成的球体周围构建多个不同半径的球形流形；将每个球形流形展开成2D表示；将在每个所述展开的球形流形上的强度值分布标准化；为每个所述展开的球形流形计算强度熵，其中强度值被看作概率值，其中熵分布被定义；以及找到使所述熵分布最小化的半径。
根据本发明的另一方面，识别群包括使用期望最大化使向量变量Θ＝(θ1，…，θF)T的多变量包绕高斯分布NWP(Θ)的混合Nw(Θ)=Σp=1pcpNwp(Θ)]]>拟合于通过所述边界流形突出的对象，所述边界流形遵循最小描述长度准则；其中混合分量概率cp在期望最大化中被估计，其中在每一维中，θi满足θ＝χw＝χmod 2π∈(-π，π]，NWP(Θ)满足Nwp(Θ)=Σk1=-∞∞···ΣkF=-∞∞Np(Θ+2πk1e1+···+2πkFeF)]]>，其中ek＝(0，...，0，1，0，...，0)T是为第k个欧几里得基础向量，在第k项记为1和其它地方记为0，其中混合分量p的估计 和 是在期望最大化内从样本集X＝{θ(1)，…，9(M)}中基于方向平均值(μ^θ)f=arg(1MΣm=1Mexp(iθf(m)))]]>和协方差Σ^θ=1M-1Σm=1MΘ(m)′Θ(m)′T]]>获得的，其中θ(m)′=(θf(m)-(μ^θ)f)mod2π]]>，并且其中观察(observation)X直接从2D展开图像I(θ，)中提取，其中每个取样的(θm，m)∈(-π，π]×(-π，π]的出现次数与相应的图像矩阵值I(θm，m)成比例地被设定。
根据本发明的另一方面，所述方法包括使用凝聚分级聚类来使相互的预定距离内的群合并，使用用于多变量包绕高斯分布对的等于18((μ2-μ1)mos2π)T(Σ1+Σ22)-1]]>((μ2-μ1)mod2π)+12ln|Σ1+Σ2||Σ1||Σ2|]]>的距离度量，其中μ1和μ2是高斯分布对的平均值，以及∑1和∑2是其各自的方差。
根据本发明的另一方面，肺结构的类由与靶结构相关的包绕高斯分量群的数量确定，其中孤立的结节有0个群，粘连结节有2个群，脉管有4个群，而血管交叉点有6个或更多群。
根据本发明的另一方面，提供了一种计算机可读的程序存储设备，有形地体现由计算机可执行的指令程序以实施在数字化图像中对肺结构进行分类的方法步骤。


图1(a)-(c)图示了根据本发明一个实施例的肺结构分类的方法。
图2(a)-(g)图示了根据本发明一个实施例的不同半径r的展开椭圆体和相应的图像强度直方图熵的效果。
图3图示了根据本发明一个实施例的具有方向数据的聚类(clustering)。
图4(a)-(d)和图5(a)-(d)图示了用于胸部CT扫描的本发明一个实施例的肺结构分类方法的图示例子。
图6(a)-(b)图示了根据本发明一个实施例的方向数据的例子。
图7是根据本发明一个实施例的分类方法的流程图。
图8是根据本发明一个实施例用于执行分类方法的典型的计算机系统的框图。
具体实施例方式
这里所描述的本发明典型的实施例一般包括用于在胸部扫描中对局部结构类型进行分类的系统和方法。因此，尽管本发明容许有各种修改和可代替的形式，它的特殊实施例通过在附图中举例来表示并且在这里将详细描述。然而，应该明白这并不意味着本发明局限于所公开的特定形式，而相反，本发明覆盖了落在本发明精神和范围内的所有修改、等效和替代方案。
如这里所使用的，术语“图像”指的是由离散的图像元素(例如，2-D图像的像素和3-D图像的体素)组成的多维数据。例如，图像可以是通过计算机断层摄影、磁共振成像、超声或本领域技术人员所公知的任何其他医学成像系统所采集的对主体(subject)的医学图像。图像也可以例如像这些遥感系统、电子镜检等由非医学环境提供。尽管图像可以看作是从R3至R的函数，然而本发明的方法并不限制于这些图像，并且可以应用于任意维数的图像，例如2-D图片和3-D体积。对于2或3维图像而言，图像的域典型地是2或3维图像矩形阵列，其中每个像素或体素都可以参照一组2或3相互正交的轴被寻址。适当地，本文所用的术语“数字的”或“数字化的”指的是在通过数字获取系统或从模拟图像的转化而得到的数字或数字化格式的图像或体积。
根据本发明一个实施例的分类系统包括(1)用于各向异性高斯拟合的模块，(2)在肺结构边界处的2D流形的构建，以及(3)对所述流形的稳健的聚类分析(cluster analysis)。第(2)部分使用基于熵最小化的数据驱动尺度选择。第(3)部分使用统计分析方法，例如用自动模数选择(automatic mode numberselection)的基于期望最大化(EM)的聚类、方向数据建模、以及基于Bhattacharyya距离变量的分级聚类。在该分析中高强度群的数量将直接决定肺结构类。与如脉管树重建之类的其它球形方法不同，这种方法允许对肺结构进行灵活的局部化检查。
在结节检测应用背景中，不正确检测的和分割的脉管和脉管分支结构表示假阳性(FP)病例。根据本发明一个实施例的分类方法拒绝所有这样的非结节结构，并且作为副产物，推断出研究中的肺结构类型的种类，即结节、粘连结节、脉管或脉管交叉点。而且，根据本发明一个实施例的分类解决方案可以在肺部CAD系统内起后处理过滤器的作用，以减少脉管和交叉点所引起的FP。
根据本发明一个实施例的肺部分类方法的流程图表示在图7中。根据本发明一个实施例的方法假定肺结构的近似位置例如从CAD系统、放射科医师的手册读物或报告等被提供。一按(one-click)结节分割算法可以被用来定位和分割靶结构，所述靶结构包括结节、粘连结节、脉管或脉管交叉点。现在参见该图，既定在步骤71提供结节候选位置用作该半自动分割解决方案的初始化。在步骤72，使各向异性高斯模型拟合于靶结构强度，形成靶的更精确的3D椭圆体模型。在步骤73，将这些椭圆体弯曲成3D球体。在步骤74，从弯曲成的3D球体构建边界流形。根据本发明一个实施例，这种构建包括将球体的3D表面展开成2D球坐标表示，接着确定适当的边界流形的半径。在步骤75，执行边界流形的聚类分析，接着在步骤76进行后处理。下面将描述这些步骤的细节。
参见步骤72，根据本发明一个实施例的算法基于使用高斯尺度-空间均值漂移分析和基于散度的詹森-香农自动带宽选择来稳健地使基于高斯的各向异性强度模型拟合于数据。这些技术从不精确的CAD或手动初始化提供对靶中心的精确估计。从靶结构边界中提取出3D椭圆体流形。然而，椭圆体拟合通常并不是不重要的，通过对局部结构分割的选择而减轻该任务，所述局部结构分割根据高斯参数平均值和协方差给出结节的中心和椭圆体形状的精确估计。该估计的稳健性还允许对非结节区域(例如感兴趣的脉管和脉管交叉点/分支)进行分割。
为了简化数学表示，在步骤73，对原始的感兴趣体积(VOI)进行仿射标准化。这包括使VOI弯曲以将分割的各向异性椭圆体转换成尺寸固定的各向同性球体，其被设置在VOI的中心。仿射标准化的参数、也就是说缩放方向和系数可以从由分割模块所估计的结构协方差的特征值分析中直接获得。
图1(a)-(c)图示了根据本发明一个实施例的示例性肺结构分类。图1(a)表示原始的感兴趣体积(VOI)和分割的结节候选者、用椭圆体拟合的结节结构，在这里是脉管。椭圆体拟合是从各向异性高斯拟合模块获得的。图1(b)表示原始的VOI的仿射标准化，其中椭圆体被弯曲成各向同性球体。图1(c)表示在距离rbound处分割结构的边界流形，其被展开成2D图像以及被球极坐标θ和参数化。通过三线性内插法获得图像灰度值。
再次参见图7，在步骤74，通过对适当的流形进行聚类分析来确定研究中的肺结构类型的种类，所述适当的流形从靶结构的边界区域被计算。根据本发明一个实施例的球形流形是从仿射标准化的3D图像中构建的。在几何学上，目的是表示稍微超出靶结构边界表面的球形层，使得含有关于通过表面的突出的对象(object)的信息。其形状在初始VOI中假定是椭圆体，尤其是与从基于高斯的各向异性分割中所获得的椭圆体成比例。因此，在仿射标准化的表示中，也与由中心点(abound，bbound)和半径rbound限定的各向同性球形状相对应。尽管中心点与分割的椭圆体的中心点相同，但球半径rbound将会以在下文所解释的数据驱动方式来确定。
假设一固定的rbound，边界流形表示可以从笛卡尔(x，y，z)被转换成球坐标(θ，)。这里，θ指的是方位角，指的是极角。结果是仿射标准化的椭圆体展开表示是2D图像矩阵I(θ，)。图(1c)图示了肺结构例子的结果。应注意的是，与一般惯例相反，极角包括间隔Interval＝2π(代替π)，即∈[-π，π]，导致笛卡尔体素的双倍出现。引入这种冗余的原因是以下将要介绍的聚类在两个参数方面在其相应的间隔Intervalθ和Interval上要求I(θ，)的周期性行为，也就是，I(θ+Intervalθ，+Interval)＝I(θ，)。在球坐标情况下，如果Interval＝π，则这显然不能实现。
根据本发明的一个实施例，在强度分布的熵的基础上，使用数据驱动方法来确定适当的半径rbound。为了推动这种方法，考虑图2(a)-(f)，每幅图都图示了在具有不同半径r的(θ，)域中展开的椭圆体表示，如图中所示。图2(g)表示对于半径r∈{1，…，32}而言关于图象强度所计算的相应的图像直方图熵Er。典型的图像熵可以根据 从图像强度I(θ，)计算出来。在适当地对图像强度值分布进行标准化以后，展开的流形图像被看作是2D似然函数。然后，直接用被解释为概率值的标准化强度值来计算强度熵。半径的选择包括自动地选定半径，使得高强度群由于突出的结构而在相应的流形中表现最特殊。这种由几个如图2(d)中所示的群所组成的流形图像与具有较小和较大半径的图像相比应该有较低的熵，其原因在于以下直观的论点。如图2(a)-(b)所示，较小半径使得相应的边界椭圆体穿过靶结构内部，导致具有更平的似然性的高熵值。在另一方面，如图2(e)-(f)所示，由于出现在附近的其它“非靶”结构，较大半径也引起高熵值。因此，适当的半径rbound形成熵分布Er的局部最小化。在这一方面，选定rbound位于正差商 第一次出现处，也就是，rbound=minr{r|Er+1>Er}.]]>将3D肺结构的部分转换成2D图像以后，可以应用经充分研究的、有效的和易于可视化的2D图像分析技术。如同可以从图1(c)所看到的，边界流形含有对肺结构分类有价值的信息。实际上，高强度群的数量揭露肺结构的类型，其等同于通过所限定的边界的突出的对象的数量。本发明一个实施例的分类建立在这一观察的基础上，想要进行以下域假设·边界流形中的0个群指示没有相连接的邻近结构，因此，分割的结构与孤立的结节相对应；·边界流形中的1个群指示至粘连结构的单一连接，这在很多情况下源自于附着至更大结构(例如肺壁等)的结节；·2个群指示两个连接，这对于血管最常见；以及·＞3个群指示脉管交叉点。
根据本发明的一个实施例，高强度群的数量通过聚类分析来识别，在图7的步骤75来执行。本发明的一个实施例的聚类策略的基础是基于期望最大化(EM)的高斯拟合。除了标准EM高斯聚类特性以外，本发明的一个实施例的聚类算法需要反映在(θ，)域中在2D边界流形图像的边缘处出现的连续性。特别是，通过球形角变量(θ，)被参数化的边界流形表示对应于方向数据。为了图示方向数据，参见图3的简化图示。在方向(θ，)域中适当的聚类算法应该重新获得单一群。然而，用线性代替方向建模，三种可观察到的结构中的每一种都将形成独立的群。而且，本发明一个实施例的聚类算法应该能够自动地确定模数。
方向数据可以在圆的圆周上以两维被可视化为超球面的表面上的点。图6(a)-(b)图示了方向数据的例子。产生方向数据的情况如下文所述。对于圆变量θ，加法“a+b”变成“(a+b)mod2π”，这里角度以间隔(-π，π]表示。应注意的是，在这种假定的情况下，mod(求模)算符也映射到(-π，π]。假设变量μθc和Vθc表示平均值和方差的圆对应部分。在用θ~=(θ-v)mod2π]]>表示的零方向漂移的情况下，对μθc和Vθc的合理定义应该维持不变。圆变量的不变性应该为μθ~c=(μθc-v)mod2π,]]>
Vθ~c=Vθc.]]>然而，从图6的例子中可以容易证明，所期望的不变性被违反。图6表示圆观察的简单集，其对应于ν＝-0.5π，在图6(a)中Θ＝{0.1π，0.25，0.6s}，而图6(b)中Θ~={0.6π,0.7π,-0.9π}]]>对于这些观察，对平均值和方差的无偏性最大似然(ML)估计可以被计算为μ^θ=0.3π,]]>σ^θ≈0.26π,]]>μ^θ~≈0.13π]]>和σ^θ~≈0.90π,]]>这明显违反了漂移不变性。在图中，平均值和方差的值分别用黑圆点的位置和附随弧(accompanying arc)的长度来图示。因此，对圆数据而言，平均值和方差的线性定义大大依赖于零方向，这对于适当的处理而言是不相称的行为和要求。
为了处理该情况，假定圆随机变量θ具有PDFp(θ)。在符合标准统计特性的情况下，PDF应该满足p(θ)≥0和∫-ππp(θ)dθ=1]]>变量θ表示复数eiθ，并且采用通过ρθexp(iμθc)=E[exp(iθ)]]]>限定的圆平均方向μθc和圆方差Vθc的符号，其中Vθc=1-ρθ.]]>参量pθ称为合成长度。形象地说，μθc是eiθ期望相位而ρθ是eiθ期望长度。Vθc∈
]]>测量离差(dispersion)的量。可以表明，平均值和方差的这些定义满足所期望的漂移不变性，并且能够被用作线性平均值和方差的适当的对应部分。
已经提出许多模型用于方向数据的统计建模。根据本发明的一个实施例，使用多变量包绕高斯分布，其是包绕高斯分布的扩展。线上的变量x的高斯分布N(x)可以“包绕(wrapped)”着单位半径的圆的圆周。也就是，包绕变量θ＝χw＝mod2π∈(-π，π]的包绕高斯分布Nw(θ)为Nw(θ)=Σk=-∞∞N(θ+2kπ).]]>向量变量Θ＝(θ1，…，θF)T的多变量包绕高斯分布可以类似地被定义为(1),Nw(Θ)=Σk1=-∞∞···ΣkF=-∞∞N(Θ+2πk1e1+···+2πkFeF),]]>其中ek＝(0，...，0，1，0，...，0)T是第k个欧几里得基础向量，第k项元素记为1和其它地方记为0。图3图示了两维多变量包绕高斯的例子。
已经表明，给定方向变量中的适当的小方差，可以利用(2),(μ^θ)f=arg(1MΣm=1Mexp(iθf(m)))]]>
和(3),Σ^θ=1M-1Σm=1MΘ(m)′Θ(m)′T]]>从样本集X＝{θ(1)，…，θ(M)}中获得方程式(1)的精确的平均值和协方差估计 和 ，其中θ(m)′=(θf(m)-(μ^θ)f)mod2π]]>，i2＝-1，以及“arg”是复数的相位。为了简化，对于Θ中的所有维数f都隐性地假定2π的周期和θf∈(-π，π]的范围。
期望最大化(EM)算法是用于在概率统计模型中找到参数的最大似然估计的一类统计方法，其中模型取决于未观察到的潜在的变量。EM在执行期望(E)步骤和最大化(M)步骤之间交替，其中期望步骤通过包括潜在的变量好像它们已经被观察到来计算似然性的期望，最大化步骤通过将在E步骤所找到的期望似然最大化来计算参数的最大似然估计。在步骤M中找到的参数接着被用来开始另一个E步骤，并且重复上述过程。EM算法将以迭代的方式改善初始估计θ0和构建新的估计θ1，…，θn。
如果y表示由可观察到的变量的值组成的不完全数据，而x表示缺少的数据，那么x和y合在一起就形成一个完全数据集。假设p是具有由向量θ所给定的参数的完全数据的联合概率分布函数p(y，x|θ)。这个函数提供了完全数据似然性。然后，利用Bayes规则，给定未观察到的变量的条件分布的期望是p(x|y,θ)=p(y,x|θ)p(y|θ)=p(y|x,θ)p(x|θ)∫p(y|x^,θ)p(x^|θ)dx^.]]>该公式表示仅要求关于给定未观察到的数据的观察似然性p(y|x，θ)、和未观察到的数据的概率p(x|θ)的知识。从θn得到θn+1的单独最大化步骤是θn=1=argmaxθEx[logp(y,x|θ)|y],]]>其中Ex[]表示在x的条件分布中θ被固定在θn时log p(y，x|θ)的条件期望。对数似然性log p(y，x|θ)经常被用来替代真正的似然性p(y，x|θ)，其原因在于它产生了更加容易的公式，但是仍然在与似然性相同的点处获得其最大值。换句话，在先前的参数值情况下给定所观察到的变量，θn+1是使完全数据对数似然性的条件期望(E)最大化(M)的值。通常，通过形成对数似然性的拉格朗日(Lagrangian)函数，然后相对于平均值和协方差估算导数来得到最大值。
在EM聚类算法中，规则的、线性高斯模型可以用上述的多变量包绕高斯模型代替。尤其是，一方面EQ.(1)而另一方面EQ.(2)和(3)分别在E和M步骤中替代初始线性当量。
根据本发明一个实施例，可以从采用高斯估计的有限混合和遵循最小描述长度(MDL)准则的模型选择的修正EM聚类算法来得到模数选择的模型，用以使混合中的分量的数量最小化。通常，EM聚类算法的输入是观察的样本集X＝{(θ1，1)，…，(θM，M)}，而当前数据是2D(图像)矩阵I(θ，)。为了克服这种不相容性，观察X直接从I(θ，)中被提取，其中每个取样的(θm，m)∈(-π，π]×(-π，π]的出现次数与相应的图像矩阵值I(θm，m)成比例地被设定。
当边界流形中的实际突出结构的形状之一与椭圆高斯形状不对应时，与高斯EM聚类相关的一个方面发生。在这种情形下，期望EM算法拟合具有一组高斯分量的该结构。这种效应将明显对分类产生不利影响，其中分量的数量起到了积分作用。
再次参见图7，根据本发明的一个实施例，在步骤76，应用后处理来合并适当的分量。尤其是，本发明的一个实施例的这种后处理可以被看作第二聚类分析，所述第二聚类分析对所有EM拟合的高斯分量集体行分析并将子集合并为单个群，直到某一尺度。用于这种情形的本领域公知的一种技术是凝聚分级聚类(agglomerative hierarchical clustering)。在分级聚类中，聚类空间用其元素的距离表示。在该情况下，元素是多变量包绕高斯函数，并且统计描述符被用于几何形状。高斯分布的合适的(和可计算分析的)统计距离量度是Bhattacharyya距离DBhatt(μ1,Σ1,μ2,Σ2)=18(μ2-μ1)T(Σ1+Σ22)-1(μ2-μ1)+121n|Σ1+Σ2||Σ1||Σ2|.]]>然而，DBhatt没有考虑包绕高斯的方向特征。因而，根据本发明的一个实施例提出了DBhatt的修正变型，“包绕Bhattacharyya距离”DBhattw(μ1,Σ1,μ2,Σ2)=18((μ2-μ1)mod2π)T(Σ1+Σ22)-1((μ2-μ1)mod2π)+121n|Σ1+Σ2||Σ1||Σ2|.]]>最后，包绕高斯分量群的数量决定了肺结构的类0代表孤立的结节，2×1＝2代表粘连结节，2×2＝4代表脉管，以及＞2×3＝6代表肪管交叉点。系数取2是由于在极坐标中的两倍间隔，如上面所讨论的。
本发明实施例的该方法的局限性在于以下事实，即在(θ，)域内尺度是依赖于位置的。根据本发明的一个实施例的一个替代方案将会是用von Mises-Fisher分布对方向数据建模可以避开该问题。对于d维单位随机向量的vonMises-Fisher分布采取形式f(x|μ，k)＝cd(k)exp(kμTx)，其中‖μ‖＝1，k≥0，和d≥2。标准化常数cd(k)由cd(k)=kd/2-1(2π)d/2Id/2-1(k)′]]>给出，其中Ir()表示r阶第一类修正贝塞耳(Bessel)函数。分布用平均方向(mean direction)μ和集中参数(concentration parameter)k来参数化，它表征根据f(x|μ，k)所提取的单位向量是如何强有力地集中在平均方向μ周围。然而，von Mises-Fisher分布的参数估计通常包括对贝塞耳(Bessel)函数比率的隐含方程式求解，对其不存在解析解。
根据本发明的一个实施例，对于所提出的肺结构分类进行定性试验。图4和5展示了胸部CT图像的分类图示，对每个类“结节”、“粘连结节”、“脉管”、“脉管交叉点”都有两个例子表示。
图4(a)-(d)和5(a)-(d)描述了用于胸部CT扫描的本发明一个实施例的肺结构分类方法的图示例子。每一行与一个例子的分割和确认相对应。图4的前两行关于结节对象、后两行关于粘连于肺壁的结节，而图5的行1和2显示了血管结构，行3和4是脉管交叉点。列(a)图示了在三个正交横截面中的CT VOI。分割结果用椭圆体图示。列(b)表示了初始VOI的仿射标准化，使得3D椭圆体被弯曲成球体。列(c)显示了在(θ，)域中展开的所构建的边界流形。然而要注意的是，已经引入了附加的强度阈值。采用这个步骤作为消除低强度结构的快速、简易手段，所述低强度结构可能会干扰高斯EM聚类。列(d)中的图形显示了由基于EM的算法拟合高斯混合模型的结果。虚线椭圆形与基于EM的聚类高斯分量相对应，而实线椭圆形描述了经过后处理的群。
如列(a)中所示，3D分割可以如图4中所示分割所有孤立和粘连的结节，也可以如图5中所示分割假阳性血管和脉管交叉点。在列(d)中，边界流形图像被转换成取样数据集X，如在2.2.1部分所描述的。此外，列(d)显示了基于EM的包绕高斯聚类的结果，也就是，k个分量的平均值和协方差用虚线椭圆形表示。特别要注意的是，在图4的行3、4和在图5的行3、4中(θ，)域的边界处的连续性。为了可视化的目的，已经包括了分级聚类后处理的图示。从该后处理中得到的群用k2实线椭圆形表示，其中心点和扩展对应于在一个经后处理的群中从所有包绕高斯的平均值计算得出的平均值和协方差。注意的是，如果群的基数较低，则这个图示可以形成退化的椭圆形，例如在图5的列2中。从分量数量k2推断结构类，可以验证当前的分类对所有8个例子给出了正确答案。其它情况获得了类似结果。
值得指出的是该分类的局限性，这在一些情形下可以导致误分类。在流形3D球体(对应于θ＝0和＝π)的极点处的结构在展开以后在2D图像的θ-维中不相称地变大。这种情形可以与从绘图学的现象相比较，其中北极和南极区域在2D墨卡托(Mercator)投影世界地图上比在3D球形世界地球仪上占据较大的区域。在上面图示的例子中，这种行为可以在图5的行4中观察到，其中在≈π的高强度结构遍布θ的整个范围(-π，π]。因而，当从展开的流形的尺度关系中尤其是对这些极点区域得出结论时建议谨慎。实际上这是包绕高斯建模、尤其是展开的缺点。在这一点上，应该注意的是，von Mises-Fisher建模避开了这种现象，因为假设没有展开。
应该理解的是，本发明可以以硬件、软件、固件、特殊目的的过程及其组合的各种形式实现。在一个实施例中，本发明可以以软件形式被实现为有形体现在计算机可读程序存储设备上的应用程序。该应用程序可以被上载到并且由包括任何适当的构造的机器执行。
图8是根据本发明一个实施例的用于执行分类方法的典型计算机系统的框图。现在参见图8，用于实现本发明的计算机系统81尤其可以包括中央处理单元(CPU)82、存储器83和输入/输出(I/O)接口84。计算机系统81一般通过I/O接口84与显示器85和各种输入设备86(例如鼠标和键盘)耦合。支持电路可以包括如高速缓冲存储器、电源、时钟电路和通讯总线之类的电路。存储器83可以包括随机存储器(RAM)、只读存储器(ROM)、磁盘驱动器、磁带驱动器等及其组合。本发明可以作为存储在存储器83中的例行程序87来实现，并通过CPU82执行以处理来自信号源88的信号。同样地，计算机系统81是通用计算机系统，当执行本发明的例行程序87时变成专用计算机系统。
计算机系统81还包括操作系统和微指令代码。这里所描述的各种过程和功能可以是通过操作系统执行的微指令代码的一部分或应用程序的一部分(或及其组合)。此外，其它各种外围设备可以被连接到计算机平台，例如辅助数据存储设备和打印设备。
应该进一步理解，因为在附图中所描述的一些组成系统部件和方法步骤可以以软件形式实现，所以根据本发明被编程的方式，系统部件(或过程步骤)之间的实际连接可以不同。给定这里所提供的本发明的教导，相关领域的普通技术人员将能够设想本发明的这些和类似的实施方式或配置。
尽管对本发明已经参照优选实施例进行了详细描述，本领域技术人员因该意识到，在不偏离如所附的权利要求书所阐述的本发明的精神和范围内，可以对其进行各种修改和替换。
权利要求
1.一种用于在数字化图像中对肺结构进行分类的方法，包括以下步骤在数字化三维(3D)图像中提供一个或多个靶结构的近似靶结构位置；在所述近似的靶位置周围拟合各向异性高斯模型，以产生更精确的3D靶模型和所述一个或多个靶结构的中心位置；使每个所述3D靶模型弯曲成3D球体；在每个所述弯曲成的3D球体周围构建边界流形；以及在所述边界流形上识别群，其中已经对所述一个或多个靶结构进行了分类。
2.如权利要求1所述的方法，其中所述数字化图像包括与三维栅格上的点域相对应的多个强度。
3.如权利要求1所述的方法，其中在近似的靶位置周围拟合各向异性高斯模型包括使用高斯尺度空间均值漂移分析和基于散度的詹森-香农自动带宽选择来产生所述靶结构的3D椭圆体模型，其中所述3D椭圆体的中心和尺寸与所述高斯模型的中心和协方差相对应。
4.如权利要求3所述的方法，其中弯曲所述3D靶模型包括仿射标准化所述3D椭圆体，其中缩放方向和系数是从所述各向异性高斯模型的结构协方差获得的。
5.如权利要求1所述的方法，其中构建边界流形进一步包括将弯曲成的球体的3D表面展开成2D表示，并且确定适当的边界流形的半径。
6.如权利要求5所述的方法，其中将弯曲成的球体的3D表面展开成2D表示包括将所述弯曲成的球体的表面转换成球坐标(θ，)，其中∈[-π，π]和θ∈[ -π，π]。
7.如权利要求5所述的方法，其中确定适当的边界流形的半径包括在所述弯曲成的球体的周围构建多个不同半径的球形流形；将每个球形流形展开成2D表示；对每个所述展开的球形流形上的强度值分布进行标准化；对每个所述展开的球形流形计算强度熵，其中强度值被看作概率值，其中熵分布被定义；以及找到使所述熵分布最小化的半径。
8.如权利要求1所述的方法，其中识别群包括使用期望最大化使向量变量Θ＝(θ1，…，θF)T的多变量包绕高斯分布Nwp(Θ)的混合Nw(Θ)=Σp=1pcpNwp(Θ)]]>拟合于通过所述边界流形突出的对象，所述边界流形遵循最小描述长度准则；其中混合分量概率cp在期望最大化中被估计，其中在每一维中，θi满足θ＝xw＝xmod 2π∈(-π，π]，Nwp(Θ)满足Nwp(Θ)=Σk1=-∞∞···ΣkF=-∞∞Np(Θ+2πk1e1+...+2πkFeF)]]>，其中ek＝(0，...，0，1，0，...，0)T是第k个欧几里得基础向量，在第k项记为1和其它地方记为0，其中混合分量p的估计 和 是在期望最大化内从样本集X＝{θ(1)，…，θ(M)}中基于方向平均值(μ^θ)f=arg(1MΣm=1Mexp(iθf(m)))]]>和协方差Σ^θ=1M-1Σm=1MΘ(m)′Θ(m)′T]]>获得的，其中θ(m)′=(θf(m)-(μ^θ)f)mod2π,]]>并且其中观察X直接从2D展开图像I(θ，)提取，其中每个取样的(θm，m)∈(-π，π]×(-π，π]的出现次数与相应的图像矩阵值I(θm，m)成比例地被设定。
9.如权利要求1所述的方法，其中进一步包括使用凝聚分级聚类来使相互的预定距离内的群合并，采用用于多变量包绕高斯分布对的等于18((μ2-μ1)mod2π)T(Σ1+Σ22)-1]]>((μ2-μ1)mod2π)+12ln|Σ1+Σ2||Σ1||Σ2|]]>的距离度量，其中μ1和μ2是高斯分布对的平均值，以及∑1和∑2是其各自的方差。
10.如权利要求1所述的方法，其中肺结构的类由与靶结构相关的包绕高斯分量群的数量确定，其中孤立的结节有0个群，粘连结节有2个群，脉管有4个群，而脉管交叉点有6个或更多群。
11.一种用于在数字化图像中对肺结构进行分类的方法，包括以下步骤在数字化三维(3D)图像中提供一个或多个3D球体的靶位置，所述图像包括与三维栅格上的点域相对应的多个强度，在所述图像中每个3D球体都表示靶结构；在所述3D球体周围构建多个不同半径的球形流形；对每个所述球形流形计算强度熵，其中强度值被看作概率值，其中熵分布被定义；找到使所述熵分布最小化的半径，其中所述最小化半径限定边界流形；将边界流形的表面展开成2D球坐标(θ，)表示，其中∈[-π，π]和θ∈[-π，π]；使用期望最大化来拟合向量变量Θ＝(θ1，…，θF)T的多变量包绕高斯分布Nw(Θ)的混合Nw(Θ)=Σp=1pcpNwp(Θ)]]>，其中混合分量概率cp在期望最大化中被估计，其中θi＝(θi，i)为通过所述边界流形突出的靶结构的群，以及其中肺结构通过多个突出的群被分类。
12.如权利要求11所述的方法，进一步包括使在所述多个球形流形的每一个上的强度分布标准化。
13.如权利要求11所述的方法，其中提供一个或多个3D球体的靶位置包括在所述数字化三维(3D)图像中提供所述一个或多个靶结构的近似靶结构位置；在所述近似的靶位置周围拟合各向异性高斯模型来产生更精确的3D椭圆体靶模型和所述椭圆体模型的中心位置；以及将所述椭圆体模型仿射标准化成3D球体。
14.如权利要求11所述的方法，其中在每一维中，θi满足θ＝xw＝x mod 2π∈(-π，π]，Nw(Θ)满足Nw(Θ)=Σk1=-∞∞···ΣkF=-∞∞N(Θ+2πk1e1+···+2πkFeF)]]>，其中ek＝(0，...，0，1，0，...，0)T是第k个欧几里得基础向量，在第k项记为1和其它地方记为0，其中混合分量p的估计 和 是在期望最大化中从样本集X＝{θ(1)，…，θ(M)}中基于方向平均值(μ^θ)f=arg(1MΣm=1Mexp(iθf(m)))]]>和协方差Σ^θ=1M-1Σm=1MΘ(m)′Θ(m)′T]]>获得的，其中θ(m)′=(θf(m)-(μ^θ)f)mod2π]]>，并且其中观察X直接从2D展开图像I(θ，)提取，其中每个取样的(θm，m)∈(-π，π]×(-π，π]的出现次数与相应的图像矩阵值I(θm，m)成比例地被设定。
15.计算机可读的程序存储设备，有形地体现由计算机可执行的指令程序，来实施用于在数字化图像中对肺结构进行分类的方法步骤，所述方法包括以下步骤在数字化三维(3D)图像中提供一个或多个靶结构的近似靶结构位置；在所述近似的靶位置周围拟合各向异性高斯模型来产生更精确的3D靶模型和所述一个或多个靶结构的中心位置；使每个所述3D靶模型弯曲成3D球体；在每个所述弯曲成的3D球体周围构建边界流形；以及在所述边界流形上识别群，其中已经对所述一个或多个靶结构进行了分类。
16.如权利要求15所述的计算机可读的程序存储设备，其中所述数字化图像包括与三维栅格上的点域相对应的多个强度。
17.如权利要求15所述的计算机可读的程序存储设备，其中在近似的靶位置周围拟合各向异性高斯模型包括使用高斯尺度空间均值漂移分析和基于散度的詹森-香农自动带宽选择来产生所述靶结构的3D椭圆体模型，其中所述3D椭圆体的中心和尺寸与所述高斯模型的中心和协方差相对应。
18.如权利要求17所述的计算机可读的程序存储设备，其中弯曲所述3D靶模型包括仿射标准化所述3D椭圆体，其中缩放方向和系数是从所述各向异性高斯模型的结构协方差获得的。
19.如权利要求15所述的计算机可读的程序存储设备，其中构建边界流形进一步包括将弯曲成的球体的3D表面展开成2D表示，并且确定适当的边界流形的半径。
20.如权利要求19所述的计算机可读的程序存储设备，其中将弯曲成的球体的3D表面展开成2D表示包括将所述弯曲成的球体的表面转换成球坐标(θ，)，其中∈[-π，π]和θ∈[-π，π]。
21.如权利要求19所述的计算机可读的程序存储设备，其中确定适当的边界流形的半径包括在所述弯曲成的球体周围构建多个不同半径的球形流形；将每个球形流形展开成2D表示；将在每个所述展开球形流形上的强度值分布标准化；对每个所述展开的球形流形计算强度熵，其中强度值被看作概率值，其中熵分布被定义；以及找到使所述熵分布最小化的半径。
22.如权利要求15所述的计算机可读的程序存储设备，其中识别群包括使向量变量Θ＝(θ1，…，θF)T的多变量包绕高斯分布Nwp(Θ)的混合Nw(Θ)=Σp=1pcpNwp(Θ)]]>拟合于通过所述边界流形突出的对象，所述边界流形遵循最小描述长度准则；其中混合分量概率cp在期望最大化内被估计，其中在每一维中，θi满足θ＝xw＝x mod 2π∈(-π，π]，Nwp(Θ)满足Nwp(Θ)=Σk1=-∞∞···ΣkF=-∞∞Np(Θ+2πk1e1+...+2πkFeF)]]>，其中ek＝(0，...，0，1，0，...，0)T是第k个欧几里得基础向量，在第k项记为1和其它地方记为0，其中混合分量p的估计 和 是在期望最大化内从样本集X＝{θ(1)，…，θ(M)}中基于方向平均值(μ^θ)f=arg(1MΣm=1Mexp(iθf(m)))]]>和协方差Σ^θ=1M-1Σm=1MΘ(m)′Θ(m)′T]]>获得的，其中θ(m)′=(θf(m)-(μ^θ)f)mod2π,]]>并且其中观察X直接从2D展开图像I(θ，)提取，其中每个取样的(θm，m)∈(-π，π]×(-π，π]的出现次数与相应的图像矩阵值I(θm，m)成比例地被设定。
23.如权利要求15所述的计算机可读的程序存储设备，所述方法进一步包括使用凝聚分级聚类来使相互的预定距离内的群合并，使用用于多变量包绕高斯分布对的等于18((μ2-μ1)mod2π)T(Σ1+Σ22)-1]]>((μ2-μ1)mod2π)+12ln|Σ1+Σ2||Σ1||Σ2|]]>的距离度量，其中μ1和μ2是高斯分布对的平均值，以及∑1和∑2是其各自的方差。
24.如权利要求15所述的计算机可读的程序存储设备，其中肺结构的类由与靶结构相关的包绕高斯分量群的数量确定，其中孤立的结节有0个群，粘连结节有2个群，脉管有4个群，而脉管交叉点有6个或更多群。
全文摘要
一种用于在数字化图像中对肺结构进行分类的方法，包括在数字化三维(3D)图像中提供(71)一个或多个靶结构的近似的靶结构位置；在所述近似的靶位置周围拟合(72)各向异性高斯模型来产生更精确的3D靶模型和所述一个或多个靶结构的中心位置；使每个所述3D靶模型弯曲(73)成3D球体；在每个所述弯曲成的3D球体周围构建(74)边界流形；以及在已经对所述一个或多个靶结构进行了分类的所述边界流形上识别(75)群。
文档编号G06T7/00GK101071506SQ20071008604
公开日2007年11月14日 申请日期2007年1月25日 优先权日2006年1月25日
发明者C·巴尔曼, X·李, 冈田和典 申请人:美国西门子医疗解决公司