CN102068281B - 一种占位性病变超声图像的处理方法 - Google Patents

Abstract

本发明公开了一种占位性病变超声图像的处理方法，对采集到的超声图像进行预处理，包括去除图像周边的文字信息、滤波、边界增强、确定有效信息区域等；然后自动定位病灶区域，确定占位病变的粗略轮廓，以所述粗略轮廓作为活动轮廓模型算法的初始轮廓来提取病灶区域的精确轮廓。其实现了占位性病变超声图像的自动分割、自动提取感兴趣区域的功能，以便可以自动的对占位性病变进行诊断，从而提高了临床诊断的客观性、准确性，对于占位性病变的诊断起到了很好的辅助作用。

Description

一种占位性病变超声图像的处理方法

技术领域

本发明涉及超声图像领域，尤其涉及一种占位性病变超声图像的处理技术。

背景技术

癌症是一大类恶性肿瘤的统称。癌症的早期往往症状很少，待发展到一定阶段后才渐渐表现出一系列症状和体征，这时再采取治疗措施往往为时已晚。因此，癌症的早期诊断和早期治疗尤为关键。目前癌症早期诊断的主要方法包括甲胎蛋白测定、铁蛋白检测、超声检查、CT检查、磁共振成像、选择性造影检查等。其中，超声检查具有无创、无辐射、价格低廉、实时、可重复操作且操作简易等特点，已成为癌症早期诊断的主要手段之一。

然而，在临床上，超声检查所采集的图像主要靠临床医生凭经验利用肉眼进行观测判断。这样一来，诊断结果就有很大的主观性，容易造成误诊、漏诊；同时，病灶区域的轮廓也主要靠医生进行手工勾勒，费时费力，并且结果重现性差。目前，国内外对占位性病变超声图像的处理方法的研究仍处于初级阶段。

因此，迫切地需要一种能够客观、准确、快速的实现占位性病变超声图像自动定位、提取的技术。

发明内容

本发明的目的是提供一种占位性病变超声图像的处理方法，目的在于克服目前临床医生在利用超声对占位性病变进行诊断时过分依赖经验，使得诊断结果具有极大主观性的弊端。

本发明的技术方案如下：一种占位性病变超声图像的处理方法，其中，包括以下步骤：S1，采集医用超声图像；S2，对所述医用超声图像进行预处理，依次包括去除无效文字信息、去除噪声、边界增强和确定有效信号区域；S3，对预处理后的超声图像进一步定位病灶区域；S4，确定病灶区域的粗略轮廓；S5，以所述粗略轮廓作为活动轮廓模型算法的初始轮廓来提取病灶区域的精确轮廓。

所述的占位性病变超声图像的处理方法，其中，所述步骤S2中综合采用中值滤波方法、各向同性扩散滤波、各向异性扩散滤波三种方法来实现去除噪声和边界增强。

所述的占位性病变超声图像的处理方法，其中，所述步骤S2中确定有效信号区域的方法进一步包括：S21，对熵处理后的超声图像进行区域生长，此时图像中存在两个区域：灰度值为0的背景区域和灰度值为255的目标区域；S22，自上而下扫描图像，停止条件是直至寻找到灰度值为255的点，记录下这两个点的坐标：左边的P点和右边的M点；S23，继续扫描图像，将符合如下条件的点的坐标分别记录下来：(1)位于点P左侧灰度值是0同时位于其上方或下方灰度值是0的点，(2)位于点M右侧灰度值是0同时位于其上方或下方灰度值也是0的点；S24，比较所有点组成的矩形包含有效回波信号区域面积的大小，取其中最合适的一个作为有效信号区域。

所述的占位性病变超声图像的处理方法，其中，所述步骤S3进一步包括：S31，设置一个同心圆模板，小圆半径r初始值为10，大圆半径R初始值为20，所述模板具有伸缩性，即可随病灶大小改变自身半径大小；S32，利用所述模板对确定的有效信号区域进行扫描，同时计算模板小圆和外环内的灰度平均值比、方差比；S33，扫描完成后，所述模板的圆心停留在灰度平均值比、方差比最小或最大的第一像素点处；S34，以第一像素点为中心点进行8方向的搜索病灶区域边界点，终止条件为所述方向上所有像素点的灰度值大于阈值T＝52，将得到的所述方向上的最后一个像素点保存下来，这样便得到8个方向上的8个点；S35，计算8个点与中心点的距离以及平均距离，若某点到中心点的距离大于平均距离的2倍，则去除该点；S36，采集剩下像素点的最大最小横、纵坐标，并以这四个值做一矩形，则所述矩形即为病灶区域的大致区域，实现了病灶区域的定位。

所述的占位性病变超声图像的处理方法，其中，所述步骤S4进一步包括：S41，在以第一像素点为中心的5×5第一邻域内，寻找与所述第一邻域内灰度平均值最接近的第二像素点作为种子点；S42，对第二像素点按8邻域进行区域生长，将符合第一生长条件的点合并到第一邻域中，并重新计算灰度平均值，将与新的灰度平均值最接近的点作为新的生长点；S43，标记被生长过的区域，当区域生长算法结束后，将最外层的标记点连接起来即是病灶区域的粗略轮廓。

所述的占位性病变超声图像的处理方法，其中，所述步骤S5进一步包括：S51，对病灶区域的粗略轮廓进行离散采样，得到第一离散点；S52，根据T-Snake模型，通过所述第一离散点内能外能的计算得到第二离散点；S53，将病灶区域的图像进行网格划分；S54，计算第二离散点与网格的焦点，得到新的蛇点，并重新从S52开始执行，直到计算得到的总能量不再变化；S55，将最终确定的蛇点连接形成新的轮廓，即为病灶区域的精确轮廓。

其中，所述步骤S53中采用ACID图像分解技术对病灶区域的图像进行网格划分。

本发明的有益效果：本发明的占位性病变超声图像的处理方法与传统的处理方法相比，实现了占位性病变超声图像的自动分割、自动提取感兴趣区域的功能，以便可以自动的对占位性病变进行诊断，从而提高了临床诊断的客观性、准确性，对于占位性病变的诊断起到了很好的辅助作用。

附图说明

图1是本发明占位性病变超声图像的处理方法的流程示意图。

图2是本发明实施例的占位性病变超声图像的处理方法中确定有效信号区域方法的示意图。

图3是本发明实施例的占位性病变超声图像的处理方法中定位病灶区域的流程示意图。

图4是本发明实施例的占位性病变超声图像的处理方法中确定病灶区域的粗略轮廓的方法示意图。

图5是本发明实施例的占位性病变超声图像的处理方法中确定病灶区域的精确轮廓的方法示意图。

图6是本发明实施例的ACID图像分解示意图。

具体实施方式

本发明提供的一种占位性病变超声图像的处理方法，目的在于克服目前临床医生在利用超声对占位性病变进行诊断时过分依赖经验，使得诊断结果具有极大主观性的弊端。

请参阅图1，图1为本发明占位性病变超声图像的处理方法的流程示意图。如图所示，本发明的技术方案为对采集到的超声图像进行预处理，包括去除图像周边的文字信息、滤波、边界增强、确定有效信息区域等；然后自动定位病灶区域，确定占位病变的粗略轮廓，以所述粗略轮廓作为活动轮廓模型算法的初始轮廓来提取病灶区域的精确轮廓。其流程包括以下步骤：S1，采集医用超声图像；S2，对所述医用超声图像进行预处理，依次包括去除无效文字信息、去除噪声、边界增强和确定有效信号区域；S3，对预处理后的超声图像进一步定位病灶区域；S4，确定病灶区域的粗略轮廓；S5，以所述粗略轮廓作为活动轮廓模型算法的初始轮廓来提取病灶区域的精确轮廓。

所述步骤S2中，采用计算图像局部区域的熵，再结合区域生长算法来去除无效文字信息。因为超声图像的四周分布着很多关于仪器型号、参数设置、扫描部位和方式以及病人信息，而真正有效的回声信息则分布在图像中间部位，如果不去掉这些信息，势必会影响对于超声图像的处理。因此，去除这些文字信息、提取有效的超声回波信息区域是预处理阶段必须做的工作。通过观察超声图像可以直观看出：文字信息比较亮，即灰度级较高，同时灰度级的波动范围较小。针对这一特征，本发明实施例中，采用计算图像局部区域的熵，再结合区域生长算法来区分并去除图像中的文字信息。

熵是图像数字特征的一种统计形式，反映了图像中所含信息量的多少，也即图像的混乱度。熵值计算表达式为： $E=-\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{p}_i{\log }_2{p}_i$ 对于超声图像(一般为8位灰度图)来说，N＝255，pi是每个灰度级像素数目与图像总的像素数之间的比值。E值大，说明图像包含的信息量大，图像灰度分布的离散度更大；E值小，说明图像包含的信息量小，图像灰度分布的离散度更小。超声图像中，文字区域的E值较小。根据熵值的这一特性，首先将图像***成20×20的小区域，按照公式分别计算每个子区域的熵值，并将该熵值取代子区域的像素值。熵值取代图像中的像素值后，对所有熵值进行归一化处理到[0，255]，根据经验设置一个阈值：M，M的范围可以为80到120，这里取M＝100，灰度值低于100的子区域赋值为0，即图像中的文字信息的区域，而大于100的区域赋值为255，即图像中的超声回波信息区域。

对图像进行了熵处理后，仍然会有部分文字信息没有被除掉，但是跟中间部分的有效信息区域已经有了一定的间隔，因此，接下来利用区域生长算法对超声图像进行进一步的处理。由于图像的中间部分是有效信息区域，且经过熵处理后，中间部分像素的灰度值是255，因此首先在图像中选取一个种子点，这里选择图像的中心点为种子点；然后以种子点为生长点进行8方向联通，生长的条件是邻域点的灰度值为255；重复上述过程，直到不再有满足条件的新结点加入集合为止。

计算图像局部区域的熵，再结合区域生长算法来去除无效文字信息，可以去除超声图像中的无效文字信息，得到有效的超声回波信息区域，能有效的节省后续的模板遍历时间。然而超声图像比较复杂且分辨率较低，噪声比较严重，因此还需要对图像进行去除噪声、边界增强处理。在本发明实施例中，采用中值滤波方法、各向同性扩散滤波、各向异性扩散滤波三种方法来实现去除噪声和边界增强。如果超声图像的梯度值过大，易产生随机噪声，则采用中值滤波；在图像梯度很小的均匀区域，采用各向同性扩散滤波，对图像进行平滑；而在感兴趣区域边界处，图像具有一定的梯度，采用各向异性扩散滤波，既平滑图像、去除噪声，又达到了保护图像边界的目的，这样以来，滤波速度快，去除噪声的效果比较好。

具体来说，滤波算法如下：

其中，Δμ是各向同性的拉普拉斯算子；

是

(x，y，t)的简写，代表梯度算子；div是散度算子；

是起扩散滤波作用的退化扩散项，μ是当前迭代的图像，I是原图像，(1-g)(μ-I)是保持初始图像不过度扩散的平衡控制项；g是用来检测边界和控制扩散速度的尺度函数：g＝1/(1+k×|G_σ|²)Gσ是高斯函数，作为对图像进行平滑的卷积核： $G_{\mathrm{\σ}}(x,y)=\frac{1}{{2\mathrm{\π\σ}}^2}\exp (-\frac{x^2+y^2}{{2\mathrm{\σ}}^2})$ 本***中依然沿用了算法中各个参数的设置：(1).g函数中参数的取值为：

其中高斯卷积核的参数取值为：σ＝2，迭代步长Δt＝0.48。

(2).λ的作用为调整平衡控制项，取值为1；(3).Y表示图像中像素点处的梯度值，取其经验值分别为Ymax＝120，Ymin＝1。

(4).迭代终止条件为： ${\left(\mathrm{ISC}\right)}_{k+1}=100\×\frac{|{\mathrm{RSNR}}_{k+1}-{\mathrm{RSNR}}_k|}{{\mathrm{RSNR}}_k}\≤\mathrm{\Ω}$ 当ISC满足条件时，图像自适应地停止迭代，其中Ω＝1。

以上步骤达到了对超声图像的噪声去除、边界增强的效果。然而，为了超声图像缩短处理时间、实现实时处理，还需要在超声图像内确定一个矩形，其包含尽可能多的回波信号以及尽可能少的无效信号，最为重要的一点是，病灶区域应该被包括在该矩形内，这里将这个矩形称为有效信号区域。

请参阅图2，在本发明实施例中，有效信号区域确定流程如下：S21，对熵处理后的超声图像进行区域生长，此时图像中存在两个区域：灰度值为0的背景区域和灰度值为255的目标区域；S22，自上而下扫描图像，停止条件是直至寻找到灰度值为255的点，记录下这两个点的坐标：左边的P点和右边的M点；S23，继续扫描图像，将符合如下条件的点的坐标分别记录下来：(1)位于点P左侧灰度值是0同时位于其上方或下方灰度值是0的点，(2)位于点M右侧灰度值是0同时位于其上方或下方灰度值也是0的点；S24，比较所有点组成的矩形包含有效回波信号区域面积的大小，取其中最合适的一个作为有效信号区域。其包含尽可能多的回波信号以及尽可能少的无效信号。

对于预处理后的超声图像，需要进一步进行定位病灶区域。在本发明的超声图像的处理方法中，采用一个同心圆模板遍历图像，以其内圆、外环的平均灰度比和方差比来定位病灶区域。

请继续参阅图3，其具体步骤如下：S31，设置一个同心圆模板，小圆半径r初始值可以为8到12，这里取r＝10，大圆半径R初始值可以为18到22，这里取R＝20，所述模板具有伸缩性，即可随病灶大小改变自身半径大小；S32，利用所述模板对确定的有效信号区域进行扫描，同时计算模板小圆和外环内的灰度平均值比、方差比；S33，扫描完成后，所述模板的圆心停留在灰度平均值比、方差比最小或最大的第一像素点处；S34，以第一像素点为中心点进行8方向的搜索病灶区域边界点，终止条件为所述方向上所有像素点的灰度值大于阈值T＝52，将得到的所述方向上的最后一个像素点保存下来，这样便得到8个方向上的8个点；S35，计算8个点与中心点的距离以及平均距离，若某点到中心点的距离大于平均距离的2倍，则去除该点；S36，采集剩下像素点的最大最小横、纵坐标，并以这四个值做一矩形，则所述矩形即为病灶区域的大致区域，实现了病灶区域的定位。

其中，同心圆模板在扫描过程中，会计算内圆、外环内的基于像素灰度值的平均值比、方差比：如果是弱回声类型，则会停留在比值最小的像素点处；如果是强回声类型，则会停留在比值最大的像素点处。

需要注意地是，对于等回声或混合回声类型，由于灰度特征不明显，因此模板设置成矩形，并对有效信息区域进行扫描，并计算矩形内的纹理特征值。超声图像的纹理特征主要有灰度共生矩阵、统计特征矩阵、傅里叶功率谱、灰度差分、Laws纹理能量等。灰度共生矩阵对超声图像具有很好的描述能力，本发明中超声图像的处理方法，利用熵、角二阶矩(能量)定位等回声病灶区域。

定义了病灶区域后，下面进一步来确定病灶区域的粗略轮廓。在本发明实施例中，采用区域生长算法来确定病灶区域的粗略轮廓，其算法容易实现，而且可以比较客观地描述出病灶区域的粗略轮廓。

请参阅图4，其方法流程如下：S41，在以第一像素点为中心的5×5第一邻域内，寻找与所述第一邻域内灰度平均值最接近的第二像素点作为种子点；S42，对第二像素点按8邻域进行区域生长，将符合第一生长条件的点合并到第一邻域中，并重新计算灰度平均值，将与新的灰度平均值最接近的点作为新的生长点；S43，标记被生长过的区域，当区域生长算法结束后，将最外层的标记点连接起来即是病灶区域的粗略轮廓。

得到了病灶区域的粗略轮廓后，以所述粗略轮廓作为活动轮廓模型算法的初始轮廓来提取病灶区域的精确轮廓。

在本发明实施例中，采用一种T-Snake模型。T.McInerney等提出了一种结合ACID(Affine Cell Image Decomposition)技术的T-Snake模型，该模型保留了传统模型的诸多优点，同时又实现了目标曲线的拓扑可变性。跟传统Snake模型类似，该模型仍定义为一系列节点组成的闭合轮廓，内力保持曲线的光滑性、连续性，并且可以通过弹性力及其他约束条件进行交互，膨胀力推动模型向图像边界方向进化，直到与外部图像力相平衡为止，模型的进化过程仍由拉格朗日运动方程描述；作为改进的地方，T-Snake的节点在模型的进化过程中并不是一成不变的，该模型将图像划分成网格状，在随着内外力作用的过程中，重新计算模型与网格的交点，找到一个新的点集作为模型再进化的新节点，从而实现了曲线的拓扑性。

T-Snake的进化过程可以认为是一种向前传输的过程，这一过程分为两个阶段：第一个阶段跟传统的活动轮廓模型类似；第二个阶段模型将计算与网格的焦点以形成新的点集，实现了拓扑结构的改变。因此T-Snake模型与传统轮廓模型的最大区别在于其能根据图像的边界信息对目标曲线进行拓扑结构上的改变，从而达到了检测比较复杂结构图像边界、以及实现了对多目标图像的分割。

T-Snake模型定义为一系列节点(设定其为N个)连接而成的闭合轮廓，其中第i(i＝0，1，2……N-1)个节点的运动方程为： ${\mathrm{\γ}}_i{\stackrel{\→}{x}}_i^{\′}\left(t\right)+a{\stackrel{\→}{\mathrm{\α}}}_i\left(t\right)+b{\stackrel{\→}{\mathrm{\β}}}_i\left(t\right)={\stackrel{\→}{\mathrm{\ρ}}}_i\left(t\right)+{\stackrel{\→}{f}}_i\left(t\right)$ 其中

是第i个节点的位置函数，

是其运动速度，γi为阻尼系数。方程左边后两项是模型的内部能量，右边两项是外部能量。假定节点的运动是没有惯性的，那么当模型达到内、外两种能量的平衡状态时，节点的运动速度变为零，此时的节点即是目标的边界点。由上式得到经过Δt时间后第i个节点的位置为： ${\stackrel{\→}{x}}_i(t+\mathrm{\Δt})={\stackrel{\→}{x}}_i\left(t\right)-\frac{\mathrm{\Δt}}{{\mathrm{\γ}}_i}(a{\stackrel{\→}{\mathrm{\α}}}_i\left(t\right)+b{\stackrel{\→}{\mathrm{\β}}}_i\left(t\right)-{\stackrel{\→}{\mathrm{\ρ}}}_i\left(t\right)-{\stackrel{\→}{f}}_i\left(t\right))$ T-Snake模型的内力结构中，

为第i个节点处的拉伸力，表示为： ${\stackrel{\→}{\mathrm{\α}}}_i\left(t\right)=2{\stackrel{\→}{x}}_i\left(t\right)-{\stackrel{\→}{x}}_{i-1}\left(t\right)-{\stackrel{\→}{x}}_{i+1}\left(t\right)$ 其是第i个节点位置函数二阶导数的离散近似，参数a控制着拉伸力的强度。是第i个节点处的弯曲力，表达式为： ${\stackrel{\→}{\mathrm{\β}}}_i\left(t\right)=2{\stackrel{\→}{\mathrm{\α}}}_i\left(t\right)-{\stackrel{\→}{\mathrm{\α}}}_{i-1}\left(t\right)-{\stackrel{\→}{\mathrm{\α}}}_{i+1}\left(t\right)$ 其是第i个节点位置函数的离散近似，参数b控制着弯曲力的强度。

T-Snake的外力结构中，

为第i个节点的膨胀力，表达式为： ${\stackrel{\→}{\mathrm{\ρ}}}_i\left(t\right)q\·\left(I\left({\stackrel{\→}{x}}_i\left(t\right)\right)\right)\·{\stackrel{\→}{n}}_i\left(t\right)$ 其中，

是模型在第i个节点处的单位法向量，q是膨胀力的幅值，控制着膨胀力的强度。将

记作F(I(x，y))，其为如下二值函数：

其中，I(x，y)是图像在(x，y)处的像素值，T为一个与图像像素值有关的阈值，函数F使T-Snake收缩(I(x，y)＜T)，从而阻止轮廓运动到背景上。考虑到图像像素值的统计特性，将F扩展为：

μ是图像像素值的平均值，σ是图像像素值的标准差，k为自定义参数，控制轮廓的膨胀范围。由以上膨胀力

的结构可以看出，

控制着轮廓点的移动方向及强度，能够使T-Snake快速收敛到图像边缘附近，但不能确保收敛到图像有效边缘上，因此引入另一种基于图像梯度特征的外力： ${\stackrel{\→}{f}}_i\left(t\right)=p\·\▿P\left({\stackrel{\→}{x}}_i\left(t\right)\right)$ 其中，p为该外力的幅值，满足： $P(x,y)=-c\left|\right|\▿[G_{\mathrm{\σ}}*I(x,y)]\left|\right|$ 图像I经过标准差为σ的高斯函数平滑滤波，去除噪声，而c为梯度势能的幅值。相当于梯度的变化。当

为正时，节点膨胀，反之则收缩。在T-snake中，往往设置p值略大于q，从而使得当T-snake到达边缘附近时，由于

大于

模型继续向着梯度变大的方向演变，直到节点处的梯度变为最大为止，此时，T-snake收敛于图像的真实边缘。即，

使得模型快速收敛到真实边缘附近，而

则保证模型收敛到真实边缘上。

请参阅图5，其是本发明实施例的占位性病变超声图像的处理方法中确定病灶区域的精确轮廓的方法示意图。所述确定精确轮廓的方法流程包括以下步骤：S51，对病灶区域的粗略轮廓进行离散采样，得到第一离散点；S52，根据T-Snake模型，通过所述第一离散点内能外能的计算得到第二离散点；S53，将病灶区域的图像进行网格划分；S54，计算第二离散点与网格的焦点，得到新的蛇点，并重新从S52开始执行，直到计算得到的总能量不在变化；S55，将最终确定的蛇点连接形成新的轮廓，即为病灶区域的精确轮廓。

具体来说，所述S51中，对得到的粗略轮廓进行离散采样，取得第一离散点，即为蛇点。本发明实施例中以同心圆模板的圆心为中心点，发射15条射线，位于射线方向上的点即为第一离散点，另外，第一离散点按逆时针方向进行存储。

S52中，所述第一离散点的内力外力达到平衡，能量函数取得最小值，其主要包括三方面：内能设定、外能设定、最小化能量设定。

(1)内能设定T-Snake模型的内能和传统活动轮廓模型的内能相同，起着保持轮廓光滑性、连续性的作用。其内力结构式为： $F_{\mathrm{int}}=a{\stackrel{\→}{\mathrm{\α}}}_i\left(t\right)+b{\stackrel{\→}{\mathrm{\β}}}_i\left(t\right)$ 转化为能量时，内能可以表示为连续能和曲率的线性组合。连续能一般为曲线的一阶导数积分，曲率则是轮廓曲线的二阶导数积分： $E_{\mathrm{int}}\left(r\left(s\right)\right)={\∫}_0^1\mathrm{\α}{\left|r^{\′}\left(s\right)\right|}^2\mathrm{ds}+{\∫}_0^1\mathrm{\β}{\left|r^{\′\′}\left(s\right)\right|}^2\mathrm{ds}$ 在计算过程当中，一般利用蛇点的差分格式代替微分，因此，当蛇点靠的比较近时，连续能的值较小。如果目标曲线只有连续能的作用，那么轮廓会收缩为一个点，无法得到目标的真实边界。为解决这个问题，在活动轮廓模型连续能项中加入了轮廓平均长度的控制，计算式为： $\mathrm{length}={\left(\underset{0}{\overset{1}{\∫}}{(x^{\′}{\left(s\right)}^2+y^{\′}{\left(s\right)}^2)}^{\frac{1}{2}}\mathrm{ds}\right)}^2$ 若节点个数为M，则目标曲线的平均长度为：c＝length/M。

因此，内部能量的表达式为： $E_{\mathrm{int}}=\mathrm{\α}\underset{0}{\overset{1}{\∫}}{|{\left|r^{\′}\left(s\right)\right|}^2-c|}^2\mathrm{ds}+\mathrm{\β}\underset{0}{\overset{1}{\∫}}{\left|r^{\′\′}\left(s\right)\right|}^2\mathrm{ds}$ 这样地话，加入了轮廓平均长度控制项后，当连续能在能量最小化过程当中趋于最小时，蛇点间的距离会趋向于平均长度c，从而避免了目标曲线收缩为一个点的情况，能更加准确地提取病灶区域的精确轮廓。

(2)外能设定内部能量仅仅保持曲线的连续性、光滑性，外部能量则将曲线推向目标的真实边界。

本发明中仍采用传统参数活动量轮廓模型中的图像力： $E_{\mathrm{image}}=-c{|\▿(G_{\mathrm{\σ}}(x,y)*I(x,y))|}^2$ 其中，Gσ是高斯滤波器，是梯度算子，*为卷积运算，c是梯度幅值的强度。单纯基于图像梯度的图像力只经过了高斯平滑滤波，对噪声仍然敏感，所以，为使目标曲线向着感兴趣区域的真实边界进化，加入了基于图像灰度的膨胀力以及基于图像区域的区域能量。

为使目标曲线进化到感兴趣区域的真实边界，一种基于区域的能量被引入到活动轮廓模型的能量函数中。基于图像区域的区域能量，其原理为假设图像中包含两个主要区域：感兴趣物体和背景区域。两部分区域象素有着不同的概率分布。当然两区域的假设也很容易推广到多个不同分布区域的情况，从而也适用于多目标的分割。在实现过程中，常假设物体和背景为均值和方差都不相同的高斯分布，根据似然函数设计区域能量为： $E_{\mathrm{region}}=-\underset{s\∈R}{\mathrm{\Σ}}\log P_R\left(I\left(s\right)\right)-\underset{s\∈R^{\′}}{\mathrm{\Σ}}\log P_{R^{\′}}\left(I\left(s\right)\right)$ 其中，R为目标区域，R’为背景区。为了计算方便，将上式改为： $E_{\mathrm{region}}=-\underset{s\∈R}{\mathrm{\Σ}}\log P_R\left(I\left(s\right)\right)-(C-\underset{s\∈R}{\mathrm{\Σ}}\log P_{R^{\′}}\left(I\left(s\right)\right))$ 其中C是对P_R′(I(s))在整个图像上求和的结果：C=∑logPR(I(S))。

由于C与轮廓位置无关，因此可以从能量函数中去掉，则区域能量可以进一步的简化为： $E_{\mathrm{region}}=-\underset{s\∈R}{\mathrm{\Σ}}\log \left(\frac{P_R\left(I\left(s\right)\right)}{P_{R^{\′}}\left(I\left(s\right)\right)}\right)$ 目标和背景区域的概率密度函数分别为： $P_R\left(I\left(s\right)\right)=\frac{1}{{\mathrm{\σ}}_R\sqrt{2\mathrm{\π}}}{e}^{\frac{-{(I\left(s\right)-{\mathrm{\μ}}_R)}^2}{2{{\mathrm{\σ}}_R}^2}}$ $P_{R^{\′}}\left(I\left(s\right)\right)=\frac{1}{{\mathrm{\σ}}_{R^{\′}}\sqrt{2\mathrm{\π}}}{e}^{\frac{-{(I\left(s\right)-{\mathrm{\μ}}_{R^{\′}})}^2}{{{2\mathrm{\σ}}_{R^{\′}}}^2}}$ 最后区域能量的表达式为： $E_{\mathrm{region}}=\underset{s\∈R}{\mathrm{\Σ}}((\frac{1}{{\mathrm{\σ}}_{R^{\′}}^2}-\frac{1}{{\mathrm{\σ}}_R^2})I^2\left(s\right)-2(\frac{u_{R^{\′}}}{{\mathrm{\σ}}_{R^{\′}}^2}-\frac{u_R}{{\mathrm{\σ}}_R^2})I\left(s\right)+\frac{{\mathrm{\μ}}_{R^{\′}}^2}{{\mathrm{\σ}}_{R^{\′}}^2}-\frac{{\mathrm{\μ}}_{R^{\′}}^2}{{\mathrm{\σ}}_R^2}+\log \left(\frac{{\mathrm{\σ}}_{R^{\′}}}{{\mathrm{\σ}}_R}\right))$ (3)最小化能量设定活动轮廓模型算法中曲线的进化问题，变成了一种能量函数最小化问题，将变分问题转为求最优解问题。我们最终设计的能量函数模型为： $E_{\mathrm{snake}}=\min \left(E_{\mathrm{total}}\right)=\min \underset{i=1}{\overset{N}{\mathrm{\Σ}}}[E_{\mathrm{int},i}+E_{\mathrm{image},i}+E_{\mathrm{inf\; lation},i}+E_{\mathrm{region},i})$ 综合考虑了计算复杂度、最优解等因素后，我们采用贪婪算法对蛇点进行优化，求解最小化能量函数。该算法是一种局部最优算法，只考虑当前蛇点在其邻域内的能量变化，而不考虑其他蛇点带来的影响，每次迭代中根据能量最小原则顺序调整每个蛇点的位置，直到所有蛇点的位置不再变化为止。

进一步地，所述步骤S53中采用ACID图像分解技术对病灶区域的图像进行网格划分。

ACID图像分解技术将二维空间图像分解成一系列的三角形网格，如图6所示。三角形网格避免了T-Snake模型在进化过程中出现二义性。

本发明提供的占位性病变超声图像的处理方法，实现了占位性病变超声图像的自动分割、自动提取感兴趣区域的功能，以便可以自动的对占位性病变进行诊断，从而提高了临床诊断的客观性、准确性，对于占位性病变的诊断起到了很好的辅助作用。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换、简单组合等多种变形，这些均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (1)

1.一种占位性病变超声图像的处理方法，其特征在于，包括以下步骤：

S1，采集医用超声图像；

S2，对所述医用超声图像进行预处理，预处理包括去除无效文字信息、去除噪声、边界增强和确定有效信号区域；

S3，对预处理后的超声图像进一步定位病灶区域；

S4，确定病灶区域的粗略轮廓；

S5，以所述粗略轮廓作为活动轮廓模型算法的初始轮廓来提取病灶区域的精确轮廓；

所述步骤S2中综合采用中值滤波方法、各向同性扩散滤波、各向异性扩散滤波三种方法来实现去除噪声和边界增强；

所述步骤S2中确定有效信号区域的方法包括：

S21，对熵处理后的超声图像进行区域生长，此时图像中存在两个区域：灰度值为0的背景区域和灰度值为255的目标区域；

S22，自上而下扫描图像，停止条件是直至寻找到灰度值为255的点，记录下这两个点的坐标：左边的P点和右边的M点；

S23，继续扫描图像，将符合如下条件的点的坐标分别记录下来：(1)位于点P左侧灰度值是0同时位于其上方或下方灰度值是0的点，(2)位于点M右侧灰度值是0同时位于其上方或下方灰度值也是0的点；

S24，比较所有点组成的矩形包含有效回波信号区域面积的大小，取其中最合适的一个作为有效信号区域；

所述步骤S3包括：

S31，设置一个同心圆模板，小圆半径r初始值为10，大圆半径R初始值为20，所述模板具有伸缩性，即可随病灶大小改变自身半径大小；

S32，利用所述模板对确定的有效信号区域进行扫描，同时计算模板小圆和外环内的灰度平均值比、方差比；

S33，扫描完成后，所述模板的圆心停留在灰度平均值比、方差比最小或最大的第一像素点处；

S34，以第一像素点为中心点进行8方向的搜索病灶区域边界点，终止条件为所述方向上所有像素点的灰度值大于阈值T＝52，将得到的所述方向上的最后一个像素点保存下来，这样便得到8个方向上的8个点；

S35，计算8个点与中心点的距离以及平均距离，若某点到中心点的距离大于平均距离的2倍，则去除该点；

S36，采集剩下像素点的最大最小横、纵坐标，并以这四个值做一矩形，则所述矩形即为病灶区域的大致区域，实现了病灶区域的定位；

所述步骤S4包括：

S41，在以第一像素点为中心的5×5第一邻域内，寻找与所述第一邻域内灰度平均值最接近的第二像素点作为种子点；

S42，对第二像素点按8邻域进行区域生长，将符合第一生长条件的点合并到第一邻域中，并重新计算灰度平均值，将与新的灰度平均值最接近的点作为新的生长点；

S43，标记被生长过的区域，当区域生长算法结束后，将最外层的标记点连接起来即是病灶区域的粗略轮廓；

所述步骤S5包括：

S51，对病灶区域的粗略轮廓进行离散采样，得到第一离散点；

S52，根据T-Snake模型，通过所述第一离散点内能外能的计算得到第二离散点；

S53，将病灶区域的图像进行网格划分；

S54，计算第二离散点与网格的焦点，得到新的蛇点，并重新从S52开始执行，直到计算得到的总能量不在变化；

S55，将最终确定的蛇点连接形成新的轮廓，即为病灶区域的精确轮廓；

其中，所述步骤S53中采用ACID图像分解技术对病灶区域的图像进行网格划分。

Images