CN103093474B - 基于同质体和局部能量的三维乳腺超声图像分割方法 - Google Patents

Abstract

本发明公开了一种基于同质体和局部能量的三维乳腺超声图像分割方法，包括如下步骤：输入待处理的三维乳腺超声图像；对待处理的图像进行边缘检测；对检测出的边缘曲面进行拟合；定义同质体；定义总能量泛函；数值计算；本发明具有以下优点：（1）解决了传统方法中三维邻域窗口难以选择的难题，使邻域可以随体素的空间位置变化而变化，同质体在一定程度上保证了邻域不会跨越组织，可以有效区分局部外表相似却属于不同组织的像素，有利于提高肿瘤轮廓提取精度。（2）克服了传统基于全局信息的水平集模型难以克服灰度不均匀的缺点，又避免了传统基于局部区域信息的水平集模型陷入局部极值的问题，提高了对三维乳腺肿瘤轮廓提取的精度。

Description

基于同质体和局部能量的三维乳腺超声图像分割方法

技术领域

本发明涉及医学超声图像处理技术领域，特别是涉及一种基于同质体和局部能量的三维乳腺超声图像分割方法。

背景技术

乳腺癌的发病率和死亡率居女性疾病的首位，已成为导致女性死亡的杀手。在乳腺肿瘤的检测中，超声检查凭借其无创伤、无辐射、适合大规模检查等优势而备受青睐。相对于二维乳腺超声图像而言，三维乳腺超声成像能提供更全面，更丰富的图像信息，更有利于医生的观察与分析。然而，医生对超声图像的阅读存在个体差异，尤其对于数据量庞大的三维图像来说，医生的阅读尤为繁琐。为了减少医生诊断的主观性和提高诊断效率，临床应用中迫切需要发展计算机辅助诊断系统。

乳腺超声图像中的肿瘤轮廓提取是计算机辅助诊断系统的重要组成部分。然而，由于受射频场不均匀性和超声成像设备本身对超声等影响，超声图像具有灰度不均匀性、低对比度、斑点噪声等特点。其次，乳腺肿瘤的大小、形状个体差异显著，肿瘤的浸润性质造成了边缘的模糊甚至缺失。以上这些使得乳腺超声图像的肿瘤处理变得非常困难。

水平集处理方法因具有自适应拓扑变化和数值解的稳定性高等优势，被广泛地应用于超声图像的处理中。利用水平集处理方法研究乳腺超声图像的肿瘤处理问题，关键在于如何定义合适的能量泛函。通过能量泛函的构造中可以引入图像的边缘信息，全局和局部区域信息，对于灰度不均匀和受斑点噪声影响的乳腺超声图像可以取得较准确的肿瘤轮廓提取结果。

传统的水平集处理方法包括基于边缘的方法和基于区域的方法两种。经典的测地活动轮廓模型仅利用图像的灰度梯度信息对图像进行划分，由于超声图像中的边缘通常非常模糊甚至缺失，导致该模型容易从弱边界处发生泄露。Chan和Vese提出的CV模型(T.Chan,L.Vese.Active contours without edges.IEEETransactions on Image processing,2001,10(2):266-277)不利用边缘信息，而利用全局信息来处理，可以避免测地活动轮廓模型容易在弱边界处发生泄露的不足。但CV模型假设图像由前景和背景两个灰度均匀的区域构成，由于乳腺超声图像常存在灰度不均匀现象，因此很难满CV模型的假设。Lankton等人提出的基于局部区域的水平集处理方法(S.Lankton,A.Tannenbaum.Localizing region-basedactive contours.IEEE Transactions on Image Processing,2008,17(11):2029-2039)通过对区域灰度信息建模来处理图像，改善了CV模型难以处理灰度不均匀图像的不足。但该方法利用的是基于固定邻域的灰度信息，难以处理乳腺超声图像中局部结构复杂的情形，此外，该方法由于仅考虑了局部信息，容易陷入局部最优。

目前，国内、外在三维乳腺超声图像肿瘤轮廓提取方面的研究还不多见。国际上知名的有美国密歇根大学的Sahiner等人和台湾大学的陈达人等人，他们主要采用基于边缘的水平集处理方法。而国内还未见三维乳腺超声图像肿瘤轮廓提取方面的研究。

发明内容

本发明的目的在于将二维图像上所定义的同质片概念推广至三维空间上，定义一种可以保证局部邻域体素相似的，可随体素空间位置变化的新的三维邻域系统，即同质体。在此基础上，结合图像灰度，纹理和边缘信息，定义合适的能量泛函，发挥水平集方法在曲面演化数值解法方面的优势，得到可以准确处理出三维乳腺超声图像肿瘤的基于同质体和局部能量的水平集轮廓提取方法。

本发明的技术方案为：基于同质体和局部能量的三维乳腺超声图像分割方法，其特征在于，包括如下步骤：

(1)输入待处理的三维乳腺超声图像；

(2)将基于相位的边缘检测方法扩展至三维空间上，从三维乳腺超声图像中检测出边缘体素，生成原图像的初始边缘图；

(3)对每个体素三维邻域的边缘信息进行局部拟合处理，以得到每个体素点的边缘能量，生成原图像的最终边缘图；边缘拟合定义包括如下步骤：

(3.1)拟合窗口的选择；(3.2)局部坐标变化；(3.3)系数确定；(3.4)计算拟合后的边缘能量；

(4)定义同质体；同质体的定义包括如下步骤：

(4.1)三维搜索窗口的确定；(4.2)扫描线搜索；(4.3)能量赋值；(4.4)同质体的定义；

(5)定义总能量函数；总能量函数的定义包括如下步骤：

(5.1)总能量泛函的定义；(5.2)局部能量泛函的定义；(5.3)全局能量泛函的定义；(5.4)正则化项的定义；

(6)数值计算；数值计算包括如下步骤：

(6.1)对于步骤(5)定义的总能量泛函，利用梯度下降流法得到水平集演化方程；(6.2)离散化水平集演化方程；

其中，步骤(5.2)中所述的局部能量泛函的定义包括如下步骤：

假设ξ表示窄带，以窄带ξ内体素y为中心，根据水平集函数φ的正负，将图像划分为两个区域Ω₁和Ω₂；结合步骤(4.4)定义的同质体Λ(y)，Λ(y)与两个区域Ω₁和Ω₂的相交区域分别为Λ(y)∩Ω₁和Λ(y)∩Ω₂，用下式计算对应相交区域的灰度统计特征：

$u_i\left(y\right)=\frac{{\∫}_{x\∈\mathrm{\Λ}\left(y\right)\∩{\mathrm{\Ω}}_i}{(I\left(x\right)-I\left(y\right))}^2\mathrm{dx}}{{\∫}_{x\∈\mathrm{\Λ}\left(y\right)\∩{\mathrm{\Ω}}_i}\mathrm{dx}},i=\mathrm{1,2}---\left(4\right)$

利用下式计算对应相交区域的纹理统计特征：

$v_i\left(y\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\∫}_{x\∈\mathrm{\Λ}\left(y\right)\∩{\mathrm{\Ω}}_i}{\left|\right|I\left(x\right)-m_i\left(k\right)\left|\right|}^2,i=\mathrm{1,2}---\left(5\right)$

其中，K表示纹理基元类别数；m₁(k)和m₂(k)表示区域Λ(y)∩Ω₁和Λ(y)∩Ω₂上第k个纹理基元通道内平均灰度信息。

进一步地，步骤(2)中所述的边缘检测包括：

针对三维乳腺超声图像的特点，采用一种非对称方法进行边缘检测，其表达式为：

其中，v为三维图像中的体素，其坐标为(x,y,z)，为球坐标系统下三维LogGabor过滤器的方向，和θ分别为天顶角和方位角，s为过滤器尺度，T_s为与尺度相关的噪声阈值，odd_s(v)和even_s(v)分别为点v在尺度s上，奇滤波器和偶滤波器的输出， $A_s\left(v\right)=\sqrt{{\left({\mathrm{odd}}_s\left(v\right)\right)}^2+{\left(\mathrm{even}\left(v\right)\right)}^2}.$

进一步地，步骤(3)的边缘拟合包括如下步骤：

步骤(3.1)中所述的拟合窗口的选择包括如下步骤：以当前体素p为中心，选择大小为7×7×7的拟合窗口；

步骤(3.2)中所述的局部坐标变化包括如下步骤：

以体素p在拟合窗口的位置为坐标原点，建立新坐标系统(u,v,h)，假设φ为体素p的梯度方向，将拟合窗口的邻域体素投影至φ上，则u为邻域体素在φ上的位置，v为邻域体素与体素p在原坐标系统下z方向上的相对距离，h为原坐标系统下邻域体素的边缘能量；

步骤(3.3)中所述的二次曲面函数的系数确定包括如下步骤：

将拟合窗口的邻域体素按步骤(3.2)进行局部坐标变化后，得到一组(u_i,v_i,h_i)的集合；结合以二次曲面函数描述的拟合面，利用最小二乘法确定其系数；

步骤(3.4)中所述的拟合后的边缘能量的计算包括如下步骤：

利用二次函数的系数，求得二次函数的极值后，将其代入二次函数，得到当前体素p拟合后的边缘能量。

进一步地，步骤(4)的同质体定义包括如下步骤：

步骤(4.1)中所述的搜索窗口的确定包括如下步骤：以当前体素点p(x,y,z)为中心，选择立方体形窗口为搜索窗口W_p；

步骤(4.2)中所述的扫描线搜索包括如下步骤：以当前体素点p(x,y,z)为中心体素，对其和W_p 6个网格面中某网格面上的体素p'间的三维直线pp'，利用快速体素遍历算法获得沿直线pp'的网格体素集合，遵循类似的方法处理该网格面上的其余体素以及其余5个网格面上的所有体素；

步骤(4.3)中所述的能量赋值包括如下步骤：根据最终边缘图，对于三维直线pp'上的网格点p_i，将其赋值为pp'上所遍历网格体素的最大边缘能量；

步骤(4.4)中所述的同质体定义为：

Λ(p)＝{(p_i,m^p(p_i)):m^p(p_i)≥τ,p_i∈W_p}    (2)

其中，m^p(·)表示网格点p_i属于同质体Λ(p)的隶属度，将其定义为1减去步骤(4.3)计算后的p_i的能量，τ为阈值。

进一步地，步骤(5)的总能量泛函的定义包括如下步骤：

步骤(5.1)中所述的总能量泛函的定义为：

E^HVLELS＝λE^L+(1-λ)E^G+E^R    (3)

其中，E^HVLELS，E^L，E^G，E^R分别表示总能量，局部能量，全局能量，正则化项；λ为调节因子，起到调节局部能量和全局能量平衡的作用。

结合公式(4)计算的灰度统计特征和公式(5)计算的纹理统计特征，设计一种基于同质体的局部能量泛函，函数表达式如下：

E^L＝∫_y∈ξ(u₁(y)-u₂(y))·(v₁(y)-v₂(y))dy    (6)

步骤(5.3)中所述的全局能量泛函的定义为：

E^G＝∫_ξ(I(y)-c₁(y))²dy-∫_ξ(I(y)-c₂(y))²dy    (7)

其中，c₁(y)和c₂(y)分别表示窄带ξ内的点y在区域Ω₁和Ω₂中的平均灰度值；

步骤(5.4)中所述的正则化项的定义为：

E^R＝∫_ξδ(φ(y))|▽φ(y)|dy    (8)。

进一步地，步骤(6)的数值求解包括如下步骤：

步骤(6.1)中所述的水平集演化方程的构造包括如下步骤：

对公式(3)应用梯度下降法，得到水平集演化方程：

$\frac{\∂\mathrm{\φ}}{\∂t}=|\▿\mathrm{\φ}|(F+\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right))---\left(9\right)$

其中，F＝F^L+F^G，F^L和F^G分别为公式(6)和公式(7)去掉积分后的表达式，▽为梯度算子；

步骤(6.1)中所述的水平集演化方程的离散化包括如下步骤：

对水平集演化方程(9)进行离散化，得到如下表达式：

φⁿ⁺¹＝φⁿ-Δt[max(Fⁿ,0)▽⁺+min(Fⁿ,0)▽^-+κ|▽φ|]    (10)

其中，Δt为迭代步长，κ为曲面的曲率，φⁿ⁺¹和φⁿ分别为第n+1次和第n次迭代后的水平集函数。

具体地，本发明包括如下步骤：

(1)输入待处理的三维乳腺超声图像；

(2)对待处理的三维图像进行边缘检测，得到初始边缘图；

针对三维乳腺超声图像的特点，采用一种非对称方法进行边缘检测，其表达式为：

其中，v为三维图像中的体素，其坐标为(x,y,z)，为球坐标系统下三维LogGabor过滤器的方向，和θ分别为天顶角和方位角，s为过滤器尺度，T_s为与尺度相关的噪声阈值，odd_s(v)和even_s(v)分别为点v在尺度s上，奇滤波器和偶滤波器的输出， $A_s\left(v\right)=\sqrt{{\left({\mathrm{odd}}_s\left(v\right)\right)}^2+{\left(\mathrm{even}\left(v\right)\right)}^2}.$

Log Gabor滤波器公式如下：

其中，为滤波器的中心频率，ρ_i，和θ_i为第i个滤波器的中心径向频率，天顶角和方位角，σ_ρ和σ_α为标准差， ω为频率空间中某点。

(3)对初始边缘图，进行局部拟合处理，得到最终边缘图；

步骤(2)中检测出的边缘面可能存在空洞或噪声，为避免此现象，采用边缘拟合方法，即对各体素的邻域点进行局部逼近，用拟合后产生的曲面的极值作为各体素的边缘能量，以此得到最终边缘图。采用二元二次函数，即h＝au²+buv+cv²+du+ev+f描述拟合面。

边缘拟合分为以下三个步骤：

(3.1)拟合窗口的选择：以当前体素p为中心，选择大小为7×7×7的拟合窗口；

(3.2)局部坐标的变换：以体素p在拟合窗口的位置为坐标原点，建立新坐标系统(u,v,h)。假设φ为体素p的梯度方向，将拟合窗口的邻域体素投影于φ上，则u为邻域体素在φ上的位置，v为邻域体素与体素p在原坐标系统下z方向上的相对距离，h为原坐标系统下邻域体素的边缘能量。

(3.3)将拟合窗口的邻域体素按步骤(3.2)进行局部坐标变化后，得到一组(u_i,v_i,h_i)的集合。结合以二次曲面函数描述的拟合面，利用最小二乘法确定其系数。

(3.4)利用二次函数的系数，求得二次函数的极值后，将其代入二次函数，得到当前体素p拟合后的边缘能量。类似地，计算得到每个体素点拟合后的边缘能量，生成原图像的最终边缘图，从而完成边缘检测；

(4)定义同质片；利用步骤(3)中生成的最终边缘图为每个图像体素定义同质体，所述同质体为一种连通的可以保证相邻体素相似性的可变三维邻域；同质体的构造包括如下步骤：

(4.1)搜索窗口的确定：以当前体素点p(x,y,z)为中心，选择立方体形窗口为搜索窗口W_p；

(4.2)扫描线搜索：以当前体素点p(x,y,z)为中心体素，对其和W_p 6个网格面中某网格面上的体素p'间的三维直线pp'，利用快速体素遍历算法获得沿直线pp'的网格体素集合，遵循类似的方法处理该网格面上的其余体素以及其余5个网格面上的所有体素；

(4.3)能量赋值：根据最终边缘图，对于三维直线pp'上的网格点p_i，将其赋值为pp'上所有网格体素的最大边缘能量。

(4.4)同质体的定义：假设Λ(p)表示当前体素点p的同质体，则Λ_p中的每个体素p_i必须满足下式：

Λ(p)＝{(p_i,m^p(p_i)):m^p(p_i)≥τ,p_i∈W_p}    (3)

其中，m^p(·)表示网格点p_i属于同质体Λ(p)的隶属度，将其定义为1减去步骤(4.3)计算后的p_i的能量，τ为阈值。

(5)定义总能量函数；

传统的水平集轮廓提取方法或采用全局信息，或采用局部信息来定义能量函数，但对于灰度不均匀的三维乳腺超声肿瘤图像的肿瘤轮廓提取并不合适。在本发明中，采用同质体替代固定邻域，在此基础上利用同质体中的局部统计信息来定义局部能量；同时结合图像中的全局信息定义全局能量。此外，为增强曲面的光滑性，增加了一个正则化项。假设E^HVLELS，E^L，E^G，E^R分别表示总能量，局部能量，全局能量，正则化项，则总能量函数的表达式如下：

E^HVLELS＝λE^L+(1-λ)E^G+E^R    (4)

其中，λ为调节因子，起到调节局部能量和全局能量平衡的作用。局部能量，全局能量，正则化项的计算包括以下三部分：

(5.1)局部能量

假设ξ表示窄带，以窄带ξ内体素y为中心，根据水平集函数φ的正负，将图像划分为两个区域Ω₁和Ω₂。结合步骤(4.4)定义的同质体Λ(y)，Λ(y)与两个区域Ω₁和Ω₂的相交区域分别为Λ(y)∩Ω₁和Λ(y)∩Ω₂，用下式计算对应相交区域的灰度统计特征：

$u_i\left(y\right)=\frac{{\∫}_{x\∈\mathrm{\Λ}\left(y\right)\∩{\mathrm{\Ω}}_i}{(I\left(x\right)-I\left(y\right))}^2\mathrm{dx}}{{\∫}_{x\∈\mathrm{\Λ}\left(y\right)\∩{\mathrm{\Ω}}_i}\mathrm{dx}},i=\mathrm{1,2}---\left(5\right)$

利用下式计算对应相交区域的纹理统计特征：

$v_i\left(y\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\∫}_{x\∈\mathrm{\Λ}\left(y\right)\∩{\mathrm{\Ω}}_i}{\left|\right|I\left(x\right)-m_i\left(k\right)\left|\right|}^2,i=\mathrm{1,2}---\left(6\right)$

其中，K表示纹理基元类别数；m₁(k)和m₂(k)表示区域Λ(y)∩Ω₁和Λ(y)∩Ω₂上第k个纹理基元通道内平均灰度信息。

结合公式(5)计算的灰度统计特征和公式(6)计算的纹理统计特征，设计一种基于同质体的局部能量泛函，函数表达式如下：

E^L _＝∫_y∈ξ(u₁(y)-u₂(y))·(v₁(y)-v₂(y))dy    (7)

(5.2)全局能量

利用图像的全局信息有利于避免陷入局部最优，采用简化的CV模型，全局能量的计算仅在演化曲面而非图像域上，函数表达式如下：

E^G＝∫_ξ(I(y)-c₁(y))²dy-∫_ξ(I(y)-c₂(y))²dy    (8)

其中，c₁(y)和c₂(y)分别表示窄带ξ内的点y在区域Ω₁和Ω₂中的平均灰度值。

(5.3)正则化项

利用正则化有利于增强曲面的光滑性，函数表达式如下：

E^R＝∫_ξδ(φ(y))|▽φ(y)|dy    (9)

(6)数值求解；

基于同质体和局部能量的水平集三维乳腺超声图像肿瘤轮廓提取方法是基于几何形变模型的曲线演变方法。几何形变模型的基本思想是将图像数据和形变速度结合起来，使演化曲面停止于图像的边缘处。几何形变模型的数值求解通常用水平集方法。

(6.1)对公式(5)应用梯度下降法，得到水平集演化方程：

$\frac{\∂\mathrm{\φ}}{\∂t}=|\▿\mathrm{\φ}|(F+\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right))---\left(10\right)$

其中，F＝F^L+F^G，F^L和F^G分别为公式(7)和公式(8)去掉积分后的表达式，▽为梯度算子。

(6.2)对水平集演化方程(10)离散化，得到如下表达式：

φⁿ⁺¹＝φⁿ-Δt[max(Fⁿ,0)▽⁺+min(Fⁿ,0)▽^-+κ|▽φ|]    (11)

其中，Δt为迭代步长，κ为曲面的曲率，φⁿ⁺¹和φⁿ分别为第n+1次和第n次迭代后的水平集函数。假设D⁺和D^-分别是向前差分和向后差分算子，则

▽⁺＝[max(D^-x,0)²+min(D^+x,0)²+max(D^-y,0)²+min(D^+y,0)²+max(D^-z,0)²+min(D^+z,0)²]^1/2

▽^-＝[max(D^+x,0)²+min(D^-x,0)²+max(D^+y,0)²+min(D^-y,0)²+max(D^+z,0)²+min(D^-z,0)²]^1/2。

与现有技术相比，本发明具有以下有益效果：

(1)本发明不是采用固定的三维邻域(如立方体，球体)，而是利用拟合后的边缘图定义了一种可变三维邻域，即同质体。该同质体可以随体素的空间位置变化而变化，具有自适应性；更为重要的是，该同质体为能保证体素相似性的连通区域，具有同质性。一方面，同质体既解决了三维邻域窗口选择的难题。另一方面，同质体在一定程度上保证了邻域不会跨越组织，可以有效区分局部外表相似却属于不同组织的像素，从而有利于提高肿瘤轮廓提取的精度。

(2)传统水平集方法中能量函数的定义或采用全局能量，即整个图像域上的平均灰度值，或采用局部能量，即对固定邻域内的灰度统计信息建模。与传统水平集方法不同的是，本发明综合考虑了全局能量和局部能量。其中，对于全局能量，利用演化曲面而非整个图像域上的点进行计算。对于局部能量，利用同质体中的灰度信息和纹理信息进行计算，而非仅考虑固定邻域中的灰度信息。能量函数的修改使本发明基本克服了CV模型难以处理灰度不均匀的不足，也改善了Lankton方法容易陷入局部最优的现象，提高了对三维乳腺肿瘤轮廓提取的精度。

附图说明

图1为基于同质体和局部能量的三维乳腺超声图像分割方法结构框图；

图2为某体素的局部边缘曲面拟合示意图；

图3为边缘曲面拟合前在某个断层(12层)上的边缘图；

图4为边缘曲面拟合后在某个断层(12层)上的边缘图；

图5为同质体的体素搜索示意图；

图6为光滑区域某体素的同质体在二维断层上的显示；

图7为非光滑区域某体素的同质体在二维断层上的显示；

图8为本发明对三维乳腺超声囊肿图像的轮廓提取结果在二维断层上的显示与现有两种方法轮廓提取结果对比图；

图9为本发明对图8的三维轮廓提取结果与现有两种方法轮廓提取结果对比图；

图10为本发明对三维乳腺超声良性肿瘤图像的轮廓提取结果在二维断层上的显示与现有两种方法轮廓提取结果对比图；

图11为本发明对图10的三维轮廓提取结果与现有两种方法轮廓提取结果对比图；

图12为本发明对三维乳腺超声恶性肿瘤图像的轮廓提取结果在二维断层上的显示与现有两种方法轮廓提取结果对比图；

图13为本发明对图12的三维轮廓提取结果与现有两种方法轮廓提取结果对比图。

具体实施方式

下面将结合附图及具体实施方式对本发明作进一步的描述。

参见图1，本发明的基于同质体和局部能量的三维乳腺超声图像分割方法，包括以下步骤：

步骤1，输入待处理的三维乳腺超声图像，对感兴趣区域进行裁剪。裁剪的方法是确保感兴趣区域包含整个乳腺肿瘤，其目的是减少肿瘤轮廓提取的计算量。

步骤2，对裁剪后的三维图像，利用公式(1)进行边缘检测，得到初始边缘图。本实例中选取三维体数据的一个断层切片(12断层)，如图2所示，右侧为左侧图像局部放大示意图。图中像素点的灰度值越接近1，表示该点为边缘点可能性越大。

步骤3，对初始边缘图中某个体素进行边缘拟合，如图4所示，绿色圆点为邻域体素，可以看出，能较好地进行局部拟合。类似地，对所有体素进行局部拟合，达到恢复缺失边缘信息的目的，得到最终边缘图。图3为对图2进行拟合后所得到的边缘图，右侧为左侧图像局部放大示意图。

步骤4，采用步骤3中的最终边缘图为每个体素定义同质体。

同质片体构造的基本过程为：

(4a)搜索窗口的确定：以边缘图上的体素点p(x,y,z)为中心，确定大小为13×13×13的搜索窗口W_p。

(4b)扫描线搜索：对中心体素，图5所示为三维扫描线搜索示意图。图中p(x,y,z)代表中心体素，立方体代表搜索窗口W_p，p'为W_p 6个网格面上的任意体素，pp'为连接p和p'的三维空间直线(红色)，记录pp'所遍历的网格体素(紫色)，p_i为其中的某个遍历体素。按类似方法，得到该网格面上的其余体素以及其余5个网格面所有体素和中心体素三维直线上所遍历的空间体素。

(4c)能量赋值：根据最终边缘图，对于遍历过的体素p_i，将其能量赋值为pp'所遍历的网格体素的最大边缘能量。

(4d)同质体定义：对于p_i，利用公式(3)进行判断，确定其是否属于同质体。本实例中τ＝0.2。

如图6所示，白色区域代表光滑区域某体素的同质体在二维断层上的显示；如图7所示，白色区域代表非表光滑区域某体素的同质体在二维断层上的显示。

步骤5，利用公式(11)对三维乳腺超声图像进行轮廓提取，获得乳腺肿瘤曲面。

本发明的效果可以通过以下实验进一步说明：

1.实验条件

硬件平台为：Intel(R)Core(TM)i5CPU 7502.67GHz,Windows XPProfessional。

软件平台为：Matlab7.0和Visual C++6.0。

2.实验内容与结果

本发明对25组来自于临床的三维乳腺超声评估图像(恶性肿瘤4例，良性肿瘤21例)进行研究，用本发明所提出的基于同质体(Homogeneous Volume,HV)和局部能量(Local Energy,LE)的水平集(Level Set,LS)算法(HVLE-LS)与现有的基于全局信息的水平集轮廓提取方法(Chan-Vese，CV)和Lanton等人提出的基于局部信息的水平集方法(以下简称为Lanton)进行比较。

2.1定性分析

图8、图10和图12为三组在二维断层切片上的轮廓提取实例。图9、图11和图13为对应的三维的轮廓提取结果实例。

图8(a)显示了一个灰度较均匀，具有规则形状的乳腺良性囊肿，图8(b)为手工轮廓提取结果，图8(c)为CV模型的轮廓提取结果，图8(d)为Lanton方法的轮廓提取结果，图8(e)为本发明的方法的轮廓提取结果。由图8可知，三维体数据中间几个断层上乳腺囊肿与周围组织对比度较大，且囊肿内部灰度均匀，但对比度和灰度均匀性随着中间断层向上和向下逐渐降低。CV模型仅利用了全局灰度信息，没有考虑局部信息，容易从灰度不均匀造成的弱边界处泄露。Lanton方法虽然考虑了局部信息，但没考虑全局信息，容易陷入灰度不均匀导致的局部极值。本发明不仅考虑了全局灰度信息，而且通过同质体考虑了局部区域内的灰度信息和纹理信息，因此可以得到理想的轮廓提取结果，既避免了CV模型不能处理灰度不均匀的问题，也避免了Lanton方法局部错误的轮廓提取问题。

图10(a)显示了一个灰度不均匀，图像对比度低的乳腺良性纤维瘤，图10(b)为手工轮廓提取结果，图10(c)为CV模型的轮廓提取结果，图10(d)为Lanton方法的轮廓提取结果，图10(e)为本发明的方法的轮廓提取结果。由于受射频场不均匀性，超声成像设备本身及肿瘤浸润效应的影响，乳腺超声图像灰度均匀性差，且肿瘤与周边组织的对比度低。从图10可知，CV模型由于假设目标和背景为灰度均匀的区域，因此容易受灰度不均性和低对比度的影响，致使轮廓提取结果中包含了大量的正常组织，产生了一个错误的轮廓提取结果。Lanton方法仅利用图像的局部灰度信息，由于乳腺超声图像中存在由噪声或灰度不均匀导致的局部极值，因此容易陷入到局部最优。而且这两种方法都没有考虑局部邻域的纹理信息和边缘信息，而这两种信息对于灰度不均匀的超声图像是非常有利的。本发明通过构造同质体和局部能量综合利用了边缘，灰度和纹理信息，有效地克服了CV模型难以处理灰度不均匀的缺点，同时避免了Lanton方法陷入局部最优的缺点，得到的轮廓提取结果具有较好的主观视觉效果。此外，由图11也可以看出，CV模型受图像不均匀的影响，提取出的肿瘤体积较大，而本发明的处理结果与手工处理结果最为接近。

图12(a)显示了一个灰度不均匀，形状不规则的浸润性小叶癌，图12(b)为手工轮廓提取结果，图12(c)为CV模型的轮廓提取结果，图12(d)为Lanton方法的轮廓提取结果，图12(e)为本发明的方法的轮廓提取结果。由图12可知，CV模型根据全局灰度均值对图像进行划分，能提取出中间几个断层切片上灰度较均匀的乳腺肿瘤。然而，在第1行和第6行的断层上，肿瘤面积逐渐变小，灰度均匀性和图像对比度变差，CV模型受到灰度不均匀的影响，误认为非肿瘤区域中的像素为肿瘤，致使产生了一些假阳性结果。Lanton方法只依靠局部灰度信息进行划分，很容易陷入图像的局部极值，难以准确地提取出乳腺肿瘤。本发明的HVLE-LS算法充分考虑了全局和局部图像信息，因此对图12(a)取得了令人满意的轮廓提取效果。由图13(a)的三维轮廓提取结果看出，本发明的HVLE-LS算法提取的肿瘤曲面比其余两种水平集轮廓提取方法都细致。

这几组实例充分证明了本发明的方法的有效性。

2.2定量分析

除主观视觉效果以外，将三种不同水平集轮廓提取方法所提取的肿瘤区域与专家手绘制的肿瘤区域的相似性度量来评价轮廓提取效果。假设A_a是算法的轮廓提取结果标记集合，A_m是手工标准轮廓提取的标记集合，则Jaccard和Dice相似性定义如下：

$\begin{array}{c}\mathrm{JS}=\frac{|A_m\∩A_a|}{|A_m\∪A_a|}\\ \mathrm{DSC}=2\frac{|A_m\∩A_a|}{\left|A_m\right|+{|A}_a|}\end{array}$

JS和DSC值越大，轮廓提取总体性能越好。JS＝1或DSC＝1为手工轮廓提取的结果和医生轮廓提取的结果完全匹配。

CV模型对前述的三组乳腺超声图像Jaccard相似性度量分别为：74.66％，69.54％，70.22％；DSC度量分别为85.49％，82.03％，82.50％。Lanton方法对这三组图像的Jaccard相似性度量分别为：57.83％，73.85％，57.13％；DSC度量分别为73.28％，84.96％，72.72％。本发明的HVLE-LS算法对这三组图像的Jaccard相似性度量结果分别为：84.92％，87.41％，84.02％；DSC度量分别为91.84％，93.28％，91.32％。

CV模型，Lanton方法和本发明的HVLE-LS算法对25组临床三维乳腺超声图像进行提取的Jaccard相似性度量分别为65.24％，57.36％，83.23％，DSC度量分别为77.67％，72.29％，90.79％。从以上分析可知，本发明提出的HVLE-LS算法可以获得更好的轮廓提取性能。

由以上的临床实验可以说明，针对临床三维乳腺超声图像的肿瘤轮廓提取，本发明在三维乳腺超声图像的肿瘤轮廓提取上存在一定的优势，克服了现有CV模型容易受灰度不均匀影响和Lanton方法容易陷入局部最优的不足，不论是定性分析还、还是定量分析，本发明均优于现有的CV模型和Lanton方法。

综上所述，本发明针对临床三维乳腺超声图像的肿瘤轮廓提取效果明显优于现有的CV模型和Lanton方法。

Claims (6)

1.基于同质体和局部能量的三维乳腺超声图像分割方法，其特征在于，包括如下步骤：

(1)输入待处理的三维乳腺超声图像；

(2)将基于相位的边缘检测方法扩展至三维空间上，从三维乳腺超声图像中检测出边缘体素，生成原图像的初始边缘图；

(3)对每个体素三维邻域的边缘信息进行局部拟合处理，以得到每个体素点的边缘能量，生成原图像的最终边缘图；边缘拟合定义包括如下步骤：

(3.1)拟合窗口的选择；(3.2)局部坐标变化；(3.3)系数确定；(3.4)计算拟合后的边缘能量；

(4)定义同质体；同质体的定义包括如下步骤：

(4.1)三维搜索窗口的确定；(4.2)扫描线搜索；(4.3)能量赋值；(4.4)同质体的定义；

(5)定义总能量函数；总能量函数的定义包括如下步骤：

(5.1)总能量泛函的定义；(5.2)局部能量泛函的定义；(5.3)全局能量泛函的定义；(5.4)正则化项的定义；

(6)数值计算；数值计算包括如下步骤：

(6.1)利用梯度下降流法得到水平集演化方程；(6.2)离散化水平集演化方程；

其中，步骤(5.2)中所述的局部能量泛函的定义包括如下步骤：

假设ξ表示窄带，以窄带ξ内体素y为中心，根据水平集函数φ的正负，将图像划分为两个区域Ω₁和Ω₂；结合步骤(4.4)定义的同质体Λ(y)，Λ(y)与两个区域Ω₁和Ω₂的相交区域分别为Λ(y)∩Ω₁和Λ(y)∩Ω₂，用下式计算对应相交区域的灰度统计特征：

$u_i\left(y\right)=\frac{{\∫}_{x\∈\mathrm{\Λ}\left(y\right)\∩{\mathrm{\Ω}}_i}{(I\left(x\right)-I\left(y\right))}^2\mathrm{dx}}{{\∫}_{x\∈\mathrm{\Λ}\left(y\right)\∩{\mathrm{\Ω}}_i}\mathrm{dx}},i=\mathrm{1,2}---\left(4\right)$

利用下式计算对应相交区域的纹理统计特征：

$v_i\left(y\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\∫}_{x\∈\mathrm{\Λ}\left(y\right)\∩{\mathrm{\Ω}}_i}{\left|\right|I\left(x\right)-m_i\left(k\right)\left|\right|}^2,i=\mathrm{1,2}---\left(5\right)$

其中，K表示纹理基元类别数；m₁(k)和m₂(k)表示区域Λ(y)∩Ω₁和Λ(y)∩Ω₂上第k个纹理基元通道内平均灰度信息。

2.根据权利要求1所述的基于同质体和局部能量的三维乳腺超声图像分割方法，其特征在于：

步骤(2)中所述的边缘检测包括：

针对三维乳腺超声图像的特点，采用一种非对称方法进行边缘检测，其表达式为：

其中，v为三维图像中的体素，其坐标为(x,y,z)，为球坐标系统下三维LogGabor过滤器的方向，和θ分别为天顶角和方位角，s为过滤器尺度，T_s为与尺度相关的噪声阈值，odd_s(v)和even_s(v)分别为点v在尺度s上，奇滤波器和偶滤波器的输出， $A_s\left(v\right)=\sqrt{{\left({\mathrm{odd}}_s\left(v\right)\right)}^2+{\left(\mathrm{even}\left(v\right)\right)}^2}.$

3.根据权利要求1所述的基于同质体和局部能量的三维乳腺超声图像分割方法，其特征在于：步骤(3)的边缘拟合包括如下步骤：

步骤(3.1)中所述的拟合窗口的选择包括如下步骤：以当前体素p为中心，选择大小为7×7×7的拟合窗口；

步骤(3.2)中所述的局部坐标变化包括如下步骤：

以体素p在拟合窗口的位置为坐标原点，建立新坐标系统(u,v,h)，假设φ为体素p的梯度方向，将拟合窗口的邻域体素投影至φ上，则u为邻域体素在φ上的位置，v为邻域体素与体素p在原坐标系统下z方向上的相对距离，h为原坐标系统下邻域体素的边缘能量；

步骤(3.3)中所述的二次曲面函数的系数确定包括如下步骤：

将拟合窗口的邻域体素按步骤(3.2)进行局部坐标变化后，得到一组(u_i,v_i,h_i)的集合；结合以二次曲面函数描述的拟合面，利用最小二乘法确定其系数；

步骤(3.4)中所述的拟合后的边缘能量的计算包括如下步骤：

利用二次函数的系数，求得二次函数的极值后，将其代入二次函数，得到当前体素p拟合后的边缘能量。

4.根据权利要求1所述的基于同质体和局部能量的三维乳腺超声图像分割方法，其特征在于：步骤(4)的同质体定义包括如下步骤：

步骤(4.1)中所述的搜索窗口的确定包括如下步骤：以当前体素点p(x,y,z)为中心，选择立方体形窗口为搜索窗口W_p；

步骤(4.2)中所述的扫描线搜索包括如下步骤：以当前体素点p(x,y,z)为中心体素，对其和W_p 6个网格面中某网格面上的体素p'间的三维直线pp'，利用快速体素遍历算法获得沿直线pp'的网格体素集合，遵循类似的方法处理该网格面上的其余体素以及其余5个网格面上的所有体素；

步骤(4.3)中所述的能量赋值包括如下步骤：根据最终边缘图，对于三维直线pp'上的网格点p_i，将其赋值为pp'上所遍历网格体素的最大边缘能量；

步骤(4.4)中所述的同质体定义为：

Λ(p)＝{(p_i,m^p(p_i)):m^p(p_i)≥τ,p_i∈W_p}       (2)

其中，m^p(·)表示网格点p_i属于同质体Λ(p)的隶属度，将其定义为1减去步骤(4.3)计算后的p_i的能量，τ为阈值。

5.根据权利要求1所述的基于同质体和局部能量的三维乳腺超声图像分割方法，其特征在于：步骤(5)的总能量泛函的定义包括如下步骤：

步骤(5.1)中所述的总能量泛函的定义为：

E^HVLELS＝λE^L+(1-λ)E^G+E^R         (3)

其中，E^HVLELS，E^L，E^G，E^R分别表示总能量，局部能量，全局能量，正则化项；λ为调节因子，起到调节局部能量和全局能量平衡的作用；

结合公式(4)计算的灰度统计特征和公式(5)计算的纹理统计特征，设计一种基于同质体的局部能量泛函，函数表达式如下：

E^L＝∫_y∈ξ(u₁(y)-u₂(y))·(v₁(y)-v₂(y))dy       (6)

步骤(5.3)中所述的全局能量泛函的定义为：

E^G＝∫_ξ(I(y)-c₁(y))²dy-∫_ξ(I(y)-c₂(y))²dy       (7)

其中，c₁(y)和c₂(y)分别表示窄带ξ内的点y在区域Ω₁和Ω₂中的平均灰度值；

步骤(5.4)中所述的正则化项的定义为：

E^R＝∫_ξδ(φ(y))|▽φ(y)|dy         (8)。

6.根据权利要求1所述的基于同质体和局部能量的三维乳腺超声图像分割方法，其特征在于：步骤(6)的数值求解包括如下步骤：

步骤(6.1)中所述的水平集演化方程的构造包括如下步骤：

对公式(3)应用梯度下降法，得到水平集演化方程：

$\frac{\∂\mathrm{\φ}}{\∂t}=|\▿\mathrm{\φ}|(F+\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right))---\left(9\right)$

其中，F＝F^L+F^G，F^L和F^G分别为公式(6)和公式(7)去掉积分后的表达式，▽为梯度算子；

步骤(6.1)中所述的水平集演化方程的离散化包括如下步骤：

(6.2)对水平集演化方程(9)进行离散化，得到如下表达式：

φⁿ⁺¹＝φⁿ-Δt[max(Fⁿ,0)▽⁺+min(Fⁿ,0)▽^-+κ|▽φ|]   (10)

其中，Δt为迭代步长，κ为曲面的曲率，φⁿ⁺¹和φⁿ分别为第n+1次和第n次迭代后的水平集函数。