CN105389811B - 一种基于多级阈值分割的多模态医学图像处理方法 - Google Patents

Abstract

本发明公开了一种基于多级阈值分割的多模态医学图像处理方法，包括如下步骤：步骤一：获取相应的医学图像灰度图像数据；步骤二：构建高斯尺度空间模型来同时完成去噪处理；步骤三：采用各向异性非线性扩散和0范式梯度最小化模型来完成医学图像的平滑处理；步骤四：缩小OTSU图像分割算法中最佳阈值的搜算范围，以提高搜算效率；步骤五：通过有效利用空间位置信息来提取完整的分割区域，从而得出更好的分割效果。本发明在图像预处理去噪阶段建立多尺度高斯模型，具有相比单一尺度模型下更好地去噪效果，在采用多阈值OTSU分割方法这一阶段，缩小了算法阈值搜素范围，提高了分割算法运行效率。因此本发明具有实时性好，分割效果准确这两大优点。

Description

一种基于多级阈值分割的多模态医学图像处理方法

技术领域

本发明涉及医学图像处理领域，特别涉及一种基于多级阈值分割的多模态医学图像处理方法。

背景技术

作为各种医学图像处理技术如配准、融合等的基础，同时鉴于在当前的临床医学应用如临床辅助诊断、图像引导的外科手术以及放射治疗中，医学图像分割技术显示出越来越重要的临床价值。医学图像分割技术是医学图像处理的研究基础，如多模态医学图像配准、多模态医学图像融合等，其分割质量的好坏将直接决定后续配准、融合等操作能否有效的进行。在当前的临床医学研究及应用如临床辅助诊断、图像引导的外科手术和放射治疗中，医学图像分割技术发挥着极其重要的临床研究及应用价值。更为重要的是，随着医学成像技术的快速发展，医学影像已成为医护工作者们的第三只眼睛，为临床医生们提供了有效的诊断线索。当前，各种成像模式如核磁共振成像(MRI)、计算机断层成像(CT)、超声成像(US)、正电子发射计算机断层成像(PET)等已在临床医学中得到了广泛应用，因此对多模态医学图像分割技术的研究已成必然趋势。而与现有的分割技术如基于偏微分方程的分割技术、基于图论的分割技术等相比，阈值分割技术具有计算简单、高效等特点，满足临床医学中的实时性要求。

发明内容

本发明的目的是提供一种基于多级阈值分割的多模态医学图像处理方法，对多级阈值进行分割，使阈值搜素范围减少从而提高医学图像分割效率，以获得更准确的医学图像分割结果。

本发明提供的技术方案为：

一种基于多级阈值分割的多模态医学图像处理方法，包括如下步骤：

步骤一：获取相应的医学图像灰度图像数据；

步骤二：构建高斯尺度空间模型来完成去噪处理；

步骤三：采用各向异性非线性扩散和0范式梯度最小化模型对步骤二处理后的医学图像的平滑处理；

步骤四：缩小OTSU图像分割算法中最佳阈值的搜算范围；

步骤五：利用空间位置信息来提取完整的分割区域，从而得出更好的分割效果。

优选的是，步骤二包括如下分步骤：

a、采用高斯函数与图像灰度函数的卷积得到高斯尺度空间模型：

L(x,y,σ)＝G(x,y,σ)*I(x,y)

其中，*表示卷积操作，σ为尺度因子，G(x,y,σ)为高斯函数，I(x,y)为所处理的灰度图像，x,y分别为相应像素点医学图像的横坐标和纵坐标位置；

b、计算第i层高斯尺度空间信息

L_i(x,y,σ_i)＝G(x,y,σ_i)*I(x,y)，

σ_i+1＝kσ_i，k＞1,i＝1,2,3...

其中，σ₀为初始尺度因子，σ_i为第i层尺度因子；k为平滑系控制参数；

c、将高斯尺度空间模型作为背景信息，采用背景差方法，获取第i层目标图像信息：

D_i(x,y)＝|I(x,y)-L_i(x,y,σ_i)|

d、对各层目标图像信息采用加权平均方式来获取最终目标图像信息：

其中ω_i为第i层目标图像信息所占权重，n为高斯尺度空间中图像总层数。

优选的是，步骤d中，所述第i层目标图像信息所占权重ω_i满足：

优选的是，步骤三中，采用各向异性非线性扩散模型进行图像的平滑处理采用的模型满足：

其中，I_t为平滑后图像灰度函数，t为平滑次数，为梯度算子，Δ表示拉普拉斯算子，div(·)表示散度算子，c(x,y,t)为扩散系数。

优选的是，所述扩散系数c(x,y,t)满足：

其中，

或

K为控制参数。

优选的是，步骤三中，使用0范式梯度最小化模型时引入一个约束项C(S)来避免弱边缘泄露现象，所述约束项C(S)满足：

S为平滑后图像，为x方向梯度，为y方向梯度，S_p为概率为p时的图像，p为对应像素点出现的概率。

优选的是，步骤四中，多级最佳阀值满足：

其中，为区域C_k概率总和，为区域C_k内的一阶矩总和，p_i为灰度图像中灰度值为i的像素点出现的概率，l为最佳阈值个数，L为灰度级。

优选的是，步骤五中，采用漫水法进行连通区域查找，同时将零碎小区域合并到周围与其最相似的大区域中，采用区域信息和边界相结合的计算策略，从而最终形成相似信息的分割区域。

本发明的有益效果是：本发明提供了一种基于多级阈值分割的多模态医学图像处理方法，在图像预处理去噪阶段建立多尺度高斯模型，由于在多尺度下图像各项信息相互结合，为图像去噪和非均匀光照矫正均提供了保证，具有相比单一尺度模型下更好地去噪效果。在采用多阈值OTSU分割方法这一阶段，对OTSU图像分割算法做出改进，缩小了算法阈值搜素范围，从而提高分割算法运行效率，使分割结果成像的实时性变得更好。本发明通过对多机阈值分割后的结果进行空间位置信息分析与处理，从而确定有效的分割联通域，进而获取更好的分割效果。因此本发明具有实时性好，分割效果准确这两大优点。

附图说明

图1为本发明所述的基于多级阈值分割的多模态医学图像处理方法流程图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

如图1所示，本发明提供了一种基于多级阈值分割的多模态医学图像处理方法，通过对医学图像建立多尺度高斯模型，采用各向异性非线性扩散和0范式梯度最小化模型进行去噪和平滑等预处理，并对多级阈值分割OTSU算法进行改进，使阈值搜素范围减少从而提高医学图像分割效率，并对分割后的图像进行空间位置信息分析与处理来获得更准确的医学图像分割结果。具体步骤如下：

步骤一：获取相应的医学图像灰度图像数据。

通过核磁共振成像(MRI)、计算机断层成像(CT)、超声成像(US)、正电子发射计算机断层成像(PET)等方式得到医学图像，读取原始图像各像素点灰度值，归一化为0-255，并以I(x,y)表示相应灰度图像，x,y分别为相应像素点医学图像的横坐标和纵坐标位置。即I(x,y)表示在数据矩阵中第x行第y列的像素灰度值。

步骤二S120：构建高斯尺度空间并进行去噪处理。

鉴于现有的非均匀光照矫正和去噪通常仅限于单一尺度下进行，并且分开进行，带来的问题不仅是去噪效果不够理想，而且非均匀光照信息也没有得到有效矫正。为了弥补该缺陷，本发明提出了一种高斯多尺度背景模型。

高斯函数在计算机视觉以及模式识别领域中有着十分广泛的应用，其定义见下式所示：

其中，σ为尺度因子。

对于灰度函数为I(x,y)的医学图像，采用高斯函数与该图像灰度函数的卷积来表示该医学图像对应的高斯尺度空间L(x,y,σ)：

L(x,y,σ)＝G(x,y,σ)*I(x,y)

其中，*表示卷积操作。

记初始尺度因子为σ₀，第i层尺度因子为σ_i，通过上述尺度空间的定义，则第i层尺度信息L_i(x,y,σ_i)计算如下：

L_i(x,y,σ_i)＝G(x,y,σ_i)*I(x,y)

σ_i+1＝kσ_i，k＞1,i＝1,2,3...

其中，k为平滑系控制参数，用来控制平滑速率。

通过将高斯尺度空间视作背景信息，则每层尺度信息L_i(x,y,σ)所对应的目标信息D_i(x,y)可采用背景差方法获取，如下式所示：

D_i(x,y)＝|I(x,y)-L_i(x,y,σ_i)|

为了能够得到一个更加准确的目标图像，并且保留更多的目标细节信息，对各层得到的目标图像采用加权平均方式来获取高斯滤波图像计算过程如下：

其中ω_i为各层目标图像所占权重，n为高斯尺度空间中图像总层数。由于随着尺度因子σ_i不断增大提取出的目标图像也更加准确，则它所占权重也将越大，因此采用下式计算各层目标图像的权重值：

在尺度空间建立过程中，本发明将根据差值图像的变化来控制尺度空间是否停止创建，另外，拟通过γ增强操作来突出目标图像中的目标信息。由于在各尺度上均进行了高斯模糊，从而为图像去噪提供了基础。由于多尺度信息的相互结合，为图像去噪和非均匀光照矫正均提供了保证。

步骤三S130：采用各向异性非线性扩散和0范式梯度最小化模型进行医学图像的平滑处理。

经过高斯滤波后的图像D(x,y)视为原始数据，再经t次平滑后的图像数据记为I_t，

D(x,y)＝u(x,y,0)

I_t＝u(x,y,t)

其中，t为平滑次数。

其中，为梯度算子，Δ表示拉普拉斯算子，div(·)表示散度算子。

记

则

其中，(i,j)为像素点坐标位置，λ为控制参数。

由于对于每个像素点都具有相同的处理方式，通常该模型被称为各向同性线性扩散模型。但该方法通常会造成边缘被过度平滑，理想的滤波器应该是使得区域内部尽量平滑而保持边缘信息尽可能不变，为此本发明拟采用各向异性非线性扩散模型来提高局部区域内部的平滑度，该模型定义如下：

其中，c(x,y,t)扩散系数，式中g函数通常有两种选取方式：

或

其中K为控制参数。

在弱边缘区域，各向异性非线性扩散模型容易造成边缘泄露，为此需要引入一个约束项C(S)来避免弱边缘泄露现象。其中S为平滑后图像。本约束项是通过限制平滑后图像中的梯度数量来达到平滑并保持边缘效果，定义如下：

s为平滑后图像，为x方向梯度，为y方向梯度，S_p为概率为p时的图像，p为对应像素点出现的概率。

最终能量公式定义为：

其中，p为对应像素点出现的概率，λ表示约束系数，λ＞1，用来控制约束强度。I_p为在概率p时的原始图像。

该公式意为每次把平滑过的I_t分别作为平滑后图像S带入上述能量公式，使该能量函数的值最小的一组I_t则为边缘保存最好的一个平滑结果，从而选取出来作为最终平滑结果。

综上所述，本发明通过结合各向异性非线性扩散模型和0范式梯度最小化模型来达到图像平滑并保持图像边缘这一目的，从而降低同质区域内细节信息对后续分割过程的影响。

步骤四S140：通过缩小OTSU图像分割算法中最佳阈值的搜算范围，以提高算法效率。

在一幅大小为M×N、灰度级为L(通常L＝256)的灰度图像中,灰度值为i的像素个数用n_i表示,总的像素个数用n表示,其中n＝n₀+n₁+...+n_L-1用pi表示灰度图像中灰度值为i的像素点出现的概率则有：

其中,将图像中的像素按灰度值用阈值t分成背景和目标两类C₀和C₁，其中C₀＝{0,1,...,t}，C₁＝{t+1,t+2,...,L-1}则C₀和C₁各自分布的概率分别为:

其中ω₁+ω₀＝1,令ω(t)＝ω₀,则ω₁＝1-ω(t)。

C₀和C₁各自分布的均值为μ₀μ₁，其中：

C₀和C₁各自分布的方差σ₀ ²和σ₁ ²分别为：

则C₀和C₁两类的类间方差σ_B ²和类内方差σ_W ²分别为：

在Otsu方法中，最大类间方差等价于最小类内方差，其目标函数如下：

O(t)＝(1-p_t)·σ_B ²＝(1-p_t)·[ω₀(μ₀)²+ω₁(μ₁)²]

最后最佳阈值t*通过下式获得：

由单阈值推广到多级阈值，假设最佳阈值个数为l，则最佳阈值利用下式可得：

其中，为区域C_k概率总和，为区域C_k内的一介矩总和。

传统Otsu法计算阈值时,对类间方差的穷举计算为最耗时部分,因此,求解Otsu阈值时若能减少对阈值的计算次数,则可以提高计算效率。有相关文献记载有关Otsu单阈值下的阈值性质，证明了Otsu方法找出的最佳阈值是用该阈值分割出的两类的均值的平均值，利用该阈值性质，不再单纯的计算所有阈值的类间方差，然后取使类间方差最大的那个阈值。而是判断分割出的两类的均值是否相等来确定最优阈值。假设一幅图像的Otsu阈值是100，传统情况下需要从整个灰度空间穷举计算，也就是从0到255都需计算一遍才能确定最优阈值100。但是利用Otsu阈值性质就可以这样：从0开始计算到100，这时可以确定100是Otsu阈值，然后结束计算。

在上述技术方案中，假设阀值T(T₁,T₂,…T_n-1)使得类内方差最小，即对任意的T_other≠T，有采用T将图像分成n类后，第0,1,2,…n-1类的均值分别为μ₀,μ₁,…μ_n-1，则一定有

证明过程如下：

设T′＝(T₁′,T₂′,…T_n-1′),并令

对于T和T′不妨设T₁＜T₁′＜T₂,T₂＜T₂′＜T₃，…，T_n-1＜T_n-1′＜L,对于其他情况同理可证。以T′为分割阀值得到的n类均值分别是μ₀′,μ₁′,…μ_n-1′，各类所占的概率为P₀′,P₁′,...,P_n-1′，对于i₁∈[T₁+1,…,T₁′],…,i_n-1∈[T_n-1+1,…,T_n-1′]，有2i_k≤2T_k′≤μ_k-1+μ_k,k＝1,2,...,n-1。因此i_k-μ_k-1≤μ_k-i_k，由于i_k≥T_k≥μ_k-1，得i_k-μ_k-1＞0。

故(i_k-μ_k-1)²≤(μ_k-i_k)²,k＝1,2,…,n-1。

由

得到

移项得

对于i₁∈[T₁+1,…,T₁′],…,i_n-1∈[T_n-1+1,…,T_n-1′]若是其中所包含的任意i都有P_i＝0，则有若存在P_i≠0，则有这两种情形都和是最小类间方差矛盾。

步骤五S150：通过有效利用空间位置信息来提取完整的分割区域，从而得出更好的分割效果。通过以上方法改进Otsu分割算法，缩小阈值搜索范围，从而提高算法效率。

由于阈值分割技术建立在统计方法之上，而完全忽略了空间位置信息，以导致分割区域间的不连续性，并且还存在着大量的零碎小区域。为此，本发明采用漫水法(Flooding漫水填充算法)进行连通区域查找，同时将那些零碎小区域合并到周围与其最相似的大区域中。其相似度的计算本发明采用区域信息和边界相结合的计算策略，从而最终形成信息的分割区域，这些分割区域具有空间连续性，从而有助于目标区域的提取。

由于像素区域就是原图像局部区域内像素点的一种聚集，这些像素点在亮度、颜色、纹理、梯度等信息上十分相似，而在特征上差异较大的相邻像素点则被划分到不同的超像素内。根据分析很容易发现，一幅图像超像素的形成可以看做图像局部区域像素点的聚类结果。基于此，本发明拟提出一种亮度、颜色、纹理、梯度等多信息融合的局部聚类算法。本算法的基本思路如下：首先指定像素区域个数l，并在每个像素区域中选定l个种子点，假设像素点总数为N，则两相邻种子点的间距可设置为为了防止种子点被设置在目标边缘或噪声点上，需要对其进行适当调整。调整方案为：对于当前每一个种子点，在其k×k邻域范围内，选择具有最小梯度值的像素点作为新的种子点。其中k可在范围{3,5,7,...}中选取。然后将以上设置好的种子点作为聚类中心，对图像所有像素点进行局部聚类操作。对于每一像素点，根据距离种子点的位置信息，以及在亮度、颜色、纹理、梯度等特征信息上的差异分别进行聚类。将每个像素点划分到坐标位置上距其最近、特征信息上差异最小的一类中。处理完所有像素点后，重新更新所有聚类中心，计算特征均值并将其作为新的聚类中心。重复以上过程直到聚类中心不在发生改变或改变很小时迭代停止。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (1)

1.一种基于多级阈值分割的多模态医学图像处理方法，其特征在于，包括如下步骤：

步骤一：获取医学图像灰度图像数据；

步骤二：构建高斯尺度空间模型来完成去噪处理；具体包括：

a、采用高斯函数与图像灰度函数的卷积得到高斯尺度空间模型：

L(x,y,σ)＝G(x,y,σ)*I(x,y)

其中，*表示卷积操作，σ为尺度因子，G(x,y,σ)为高斯函数，I(x,y)为所处理的灰度图像，x,y分别为相应像素点医学图像的横坐标和纵坐标位置；

b、计算第i层高斯尺度空间信息

L_i(x,y,σ_i)＝G(x,y,σ_i)*I(x,y)，

σ_i+1＝kσ_i，k＞1,i＝1,2,3...

其中，σ₀为初始尺度因子，σ_i为第i层尺度因子；k为平滑系控制参数；

c、将高斯尺度空间模型作为背景信息，采用背景差方法，获取第i层目标图像信息：

D_i(x,y)＝|I(x,y)-L_i(x,y,σ_i)|

d、对各层目标图像信息采用加权平均方式来获取最终目标图像信息：

其中ω_i为第i层目标图像信息所占权重，n为高斯尺度空间中图像总层数；

步骤三：采用各向异性非线性扩散模型进行图像的平滑处理，采用的模型满足：

其中，I_t为平滑后图像灰度函数，t为平滑次数，为梯度算子，Δ表示拉普拉斯算子，div(·)表示散度算子，c(x,y,t)为扩散系数；

所述扩散系数c(x,y,t)满足：

其中，

或

K为控制参数；

使用0范式梯度最小化模型时引入一个约束项C(S)来避免弱边缘泄露现象，所述约束项C(S)满足：

S为平滑后图像，为x方向梯度，为y方向梯度，S_p为概率为p时的图像，p为对应像素点出现的概率；

步骤四：缩小OTSU图像分割算法中最佳阈值的搜算范围，其中

多级最佳阈值满足：

其中，为区域C_k概率总和，为区域C_k内的一阶矩总和，p_i为灰度图像中灰度值为i的像素点出现的概率，l为最佳阈值个数，L为灰度级；

步骤五：采用漫水法进行连通区域查找，同时将零碎小区域合并到周围与其最相似的大区域中，采用区域信息和边界相结合的计算策略，从而最终形成相似信息的分割区域。