CN102135606A - 一种基于knn分类算法的mr图像灰度不均匀性校正分割方法 - Google Patents

Abstract

一种基于KNN分类算法的MR图像灰度不均匀性校正分割方法，属于图像处理领域。本发明首先利用曲面拟合的的知识，用一组正交归一化的基函数构造灰度不均匀场模型，建立能量函数；然后依据能量函数最小化原理求解模型参数，实现灰度不均匀性校正和图像分割。其中解模型参数过程中采用迭代算法，并采用KNN算法求解隶属函数。这样可以在去除灰度不均匀场的同时，大大降低部分容积效应，同时降低了噪声对MR图像灰度不均匀性校正分割的影响。在KNN求解隶属函数时，先用核估计算法得到准确光滑的归一化直方图，然后用最大类间方差法分别求解脑脊液与灰质间的阈值TCG以及灰质与白质间的阈值TGW，利用这两个阈值对算法进行粗分类，最后采用传统的KNN算法对待定点进行精确分类。

Description

一种基于KNN分类算法的MR图像灰度不均匀性校正分割方法

技术领域

本发明属于图像处理技术领域，涉及MR图像的灰度不均匀性校正算法。

背景技术

为了提高磁共振(MR)图像的质量和图像分割、配准等处理的精确度，需要在处理前对磁共振图像的灰度不均匀性进行校正，传统的校正方法主要包括：

(1)基于滤波的方法：因为灰度不均匀场在整个图像区域中是缓慢变化的，所以可以把灰度不均匀场频谱归为低频谱。因此可以采用低通滤波器来将灰度不均匀场和真实图像分开。但是在磁共振图像中，真实图像和灰度不均匀场的频谱有重叠的部分，所以这种方法效果有限。

(2)基于统计信息的方法：该方法假设局部区域的中值等于这一类型组织的中值，而灰度不均匀场的存在使之偏移了组织中值，因此，可以利用局部图像的中值和某组织的中值之差来校正图像。但是这种假设并不是对所有MR图像都合适，有一定的局限性。

(3)基于分割的方法：这种方法在进行图像分割的同时对图像进行校正，获得不均匀场。这种方法能够有效利用像素点的信息，进行图像校正，能够自动获得校正图像和分割图像。例如基于EM(expectation-maximization，期望最大化)的灰度不均匀性校正和分割方法，该方法对初始值要求较高，需要手动选取各个组织的特征点作为初始值，分割结果受初始值影响较大；基于FCM(Fuzzy c-means clustering，模糊C均值聚类)的灰度不均匀性校正和分割方法，该方法引入一个平滑项来确保灰度不均匀场的平滑性，分割结果较好，但是算法复杂，计算量较大。

发明内容

本发明目的是提供一种基于能量最小和KNN(K-Nearest Neighbor，K最近邻)分类算法的MR图像灰度不均匀性校正分割方法。该方法属于一种基于分割的方法，可以确保分割结果不受初始值的影响，不用人工干预，而且在消除灰度不均匀性的同时，大大降低部分容积效应(partial volume effect)和噪声对分割算法的影响。

该方法首先基于曲面拟合原理，采用一组正交归一化的基函数构造灰度不均匀场模型，然后建立灰度不均匀场模型的能量函数(也就是原始MR图像和灰度不均匀场模型之间的误 差函数)；接下来通过对能量函数最小化的迭代计算，当原始MR图像和灰度不均匀场模型之间的误差达到最小时，就可以求出灰度不均匀场模型的模型参数；再将灰度不均匀场模型的模型参数代入灰度不均匀场模型中求出最终经灰度校正的MR图像分割结果。

为了方便理解和描述，先介绍本发明所涉及到的一些原理和数学模型：

(1)灰度不均匀场模型：现在广泛使用的模型是将灰度不均匀场简化为一个光滑的乘性场。以二维图像为例，即：I＝BJ。其中I为原始图像，B灰度不均匀场，J为未受灰度不均匀场污染的图像。MR图像的灰度不均匀场在空间上是缓慢变化的，表现在二维空间上就是一个光滑的曲面，可以用一组正交归一化基函数来拟合：

$B(x,y)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{w}_k{g}_k(x,y)$

对于J，同一组织的灰度值应该是恒定的，脑可以大致分为白质、灰质、脑脊液和背景4类，因此J为分段线性函数，由4个灰度值组成。

$J(x,y)=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{c}_i{u}_i(x,y)$

其中c_i为各个组织的灰度值，u_i(x，y)为隶属函数，其定义如下：

且 $\underset{i}{\mathrm{\Σ}}{u}_i(x,y)=1$

综上所述，灰度不均匀场模型如下：

$I(x,y)=B(x,y)J(x,y)=\left[\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{w}_k{g}_k(x,y)\right]\left[\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{c}_i{u}_i(x,y)\right]$

(2)最大类间方差法：最大类间方差法(Otsu方法)是阈值选取方法中较为成功的一种，通常可以用来作为初分割。其原理如下，假设一副图像的灰度值范围为[0，L)，灰度值为i的像素点个数为n₁，图像中像素点总数为n，对直方图归一化，可得：

$p_i=\frac{n_i}{n},$ p_i＞0且 $\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{p}_i=1$

假设阈值T将图像分成C₀与C₁两个类，C₀是一组灰度级为[0，1，...，T-1)的像素，C₁是一组灰 度级为[T，T+1，...，L)的像素，则C₀、C₁类出现的概率与均值分别为：

某像素点属于C₀类的概率： 

C₀类的灰度均值： 

某像素点属于C₁类的概率：w₁＝1-w₀；C₁类的灰度均值： 

其中，u为整体图像的灰度均值 

那么类C₀与C₁的类间方差为：

σ²(T)＝w₀(u₀-u)²+w₁(u₁-u)²＝w₀w₁(u₁-u₀)²

最佳的阈值应该使类间方差最大，即 

这个T^*就是我们需要找的最佳阈值。

本发明技术方案如下：

一种基于KNN分类算法的MR图像灰度不均匀性校正分割方法，如图1所示，包括以下步骤：

步骤1：建立灰度不均匀场模型。

所述灰度不均匀场模型为 

其中：B(x，y)为灰度不均匀场，J(x，y)为未受灰度不均匀场污染的MR图像；且 

其中g_k(x，y)为一组正交归一化的基函数，w_k为加权系数，k＝1，2，L，N，N为基函数的组数；c_i为MR图像中各个组织的灰度值，u_i(x，y)为隶属函数，i＝1，2，3或4；隶属函数u_i(x，y)相对于一个像素点(x，y)而言，定义为：

且 $\underset{i}{\overset{4}{\mathrm{\Σ}}}{u}_i(x,y)=1.$

步骤2：建立灰度不均匀场模型B(x，y)J(x，y)与原始MR图像I(x，y)之间的能量函数F(可以理解为误差函数)：

$F=\∫{|I(x,y)-[\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{w}_k{g}_k(x,y)\left]\right[\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{c}_i{u}_i(x,y)\left]\right|}^2\mathrm{dxdy}$

步骤3：将步骤2所得能量函数F写成向量形式，即：

F(U，C，W)＝∫|I(x，y)-(W^TG(x，y))(C^TU(x，y))|²dxdy            (1)

其中：

B(x，y)＝W^TG(x，y)，                    (2)

J(x，y)＝C^TU(x，y)，                    (3)

C＝(c₁，c₂，c₃，c₄)^T，W＝(w₁，w₂，...w_N)^T，

G(x，y)＝(g₁(x，y)，g₂(x，y)，...g_N(x，y))^T，

U(x，y)＝(u₁(x，y)，u₂(x，y)，u₃(x，y)，u₄(x，y))^T。

步骤4：任意确定灰度不均匀场模型B(x，y)J(x，y)的模型参数W和C的初始值W⁽¹⁾和C⁽¹⁾，由B(x，y)＝W^TG(x，y)计算出初始灰度不均匀场B⁽¹⁾(x，y)＝(W⁽¹⁾)^TG(x，y)，再由I(x，y)＝B(x，y)J(x，y)计算出初始未受灰度不均匀场污染的MR图像J⁽¹⁾(x，y)＝I(x，y)/B⁽¹⁾(x，y)，最后对初始未受灰度不均匀场污染的MR图像J⁽¹⁾(x，y)进行分类，求出初始隶属函数U⁽¹⁾(x，y)。

其中对初始未受灰度不均匀场污染的MR图像J⁽¹⁾(x，y)进行分类，求出初始隶属函数U⁽¹⁾(x，y)的具体过程包括以下步骤：

步骤4-1：定义灰度值概率密度函数 

其中p(q)表示灰度值为q的概率密度，n为MR图像的像素点总数，K(·)为0均值方差为1的高斯核函数，h为窗口宽度。

步骤4-2：在[0，50]像素宽度范围内确定一个合适的窗口宽度h，使得灰度值概率密度函 数p(q)曲线的波峰数量为4。灰度值概率密度函数p(q)曲线的4个波峰对应着MR图像中脑白质、脑灰质、脑脊液和背景4个区域的最大概率处的灰度值。

步骤4-3：采用最大类间方差法确定脑脊液与灰质间的阈值TCG以及灰质与白质间的阈值TGW。具体方法是：a)从步骤4-1和步骤4-2确定的灰度值概率密度函数p(q)曲线中确定一个灰度区间[q_min，q_max]，其中q_min为[5，10]之间的一个灰度值；b)在灰度区间[q_min，q_max]内，找出一个灰度值q_med，使得在灰度区间[q_min，q_med]内，利用最大类间方差法求出脑脊液与灰质间的阈值TCG，在灰度区间[TCG，q_max]内，利用最大类间方差法求出白质与灰质间的阈值TGW，且TGW＝q_med。

步骤4-4：粗分类(如图2所示)。利用步骤4-3)所得的阈值TCG与TGW，加上一定的偏差δ，将灰度区间[q_min，q_max]分为4类区域：灰度值为[q_min，TCG-δ]的脑脊液区域Ω₁；灰度值为[TCG+δ，TGW-δ]的脑灰质区域Ω₂；灰度值为[TGW+δ，q_max]的脑白质区域Ω₃；灰度值为[TCG-δ，TCG+δ]或[TGW-δ，TGW+δ]的待定区域Ω₀。

步骤4-5：采用KNN分类算法对待定区域Ω₀中的像素点进行分类。

最终得到初始隶属函数： 

其中：

i＝1，2，3或4

步骤5：利用下述迭代公式(4)(5)对灰度不均匀场模型模型参数C和W进行求解。

W(k+1)＝(A^(k))^-1V^(k)，                (4)

$c_i^{(k+1)}=\frac{\∫I(x,y)B^{\left(k\right)}(x,y)u_i^{\left(k\right)}(x,y)\mathrm{dxdy}}{\∫{\left(B^{\left(k\right)}\right)}^2(x,y)u_i^{\left(k\right)}(x,y)\mathrm{dxdy}},$ i＝1，2，3或4            (5)

其中：

B^(k)(x，y)＝(W^(k))^TG(x，y)，

J^(k)(x，y)＝(C^(k))^TU^(k)(x，y)，

V^(k)＝∫I(x，y)G(x，y)J^(k)(x，y)dxdy，

A^(k)＝∫G(x，y)G^T(x，y)(J^(k))²(x，y)dxdy。

步骤6：对隶属函数进行迭代求解。

首先由公式J^(k+1)(x，y)＝I(x，y)/B^(k)(x，y)计算出第k+1次迭代时未受灰度不均匀场污染的MR图像J^(k+1)(x，y)，接着按步骤4所述方法求出第k+1次迭代时的隶属函数U^(k+1)(x，y)。

步骤7：判断是否停止迭代。当迭代次数达到设定值k′时停止迭代(即经过k′次迭代就可以得到最优的灰度不均匀性模型模型参数)，将此时的灰度不均匀性模型模型参数分别记为记为W^*、U^*、c^* _i，即：

W^*＝W^(k′)                    (6)

${c^{*}}_i=c_i^{\left(k^{\′}\right)}$ i＝1，2，3，4                (7)

$U^{*}=(U_1^{*}(x,y),u_2^{*}(x,y),u_3^{*}(x,y),u_4^{*}(x,y){)}^T=U^{\left(k\right)}---\left(8\right)$

其中，迭代次数设定值k′可设为[10，20]之间。

步骤8：确定最终校正后的图像以及分割后的图像。

将步骤7中得到的最优灰度不均匀场模型模型参数c^* _i、U^*、W^*代入公式(1)计算出最终的灰度不均匀场B^*，然后由灰度不均匀场模型公式计算出最终校正后的MR图像J^*(x，y)，最后由公式 

计算出最终分割后的图像：背景区域的MR图像I_RACK；脑脊液区域的MR图像I_CSF；脑灰质区域的MR图像I_GM；脑白质区域的MR图像I_WM。

补充说明：

1.公式1中，选取勒让德多项式作为的基函数，以2-D为例，需3阶16个多项式。

2.步骤4-5中采用的KNN分类算法如下所述：

对于待定区域的某一个待定点(x，y)而言，搜索其K个近邻，当K个近邻中Ω_i(i＝1，2，3)类的像素点个数最多时，则将该待定点划分为Ω_i(i＝1，2，4)类。如果K个近邻中Ω₀类的像素点个数最多(仍然无法判断待定点属于哪一类)，则扩大搜索范围重复上述操作，直到能正确划分该待定点为止。

3.迭代公式(4)是通过计算能量函数F(公式(1))的偏导数 

，并令其等于0得到的；迭代公式(5)是通过计算能量函数F(公式(1))的偏导数 

，并令其等于0得到的。

4.A^-1表示A的逆矩阵，W^(K)、U^(k)此类符号表示第k次迭代W、U的值。

5.在步骤4-3中q_min取[5，10]之间的一个灰度值的目的是将背景与其他组织结构分割开来。在步骤4-4中只需要寻找脑白质区域、脑灰质区域与脑脊液区域之间的阈值。

6.本算法中建立灰度不均匀场、建立能量函数、求解迭代公式的过程均属于分析过程，实际的程序流程从初始化灰度不均匀场模型模型参数开始。

本发明的优点

(1)在医学图像中，普遍存在着部分容积效应(由于分辨率的原因，在人体组织器官交界处，一个像素点很可能对应于两个或多个组织器官，此时无法直观的对该像素点所属类别进行划分)，因此，本算法采用了KNN分类算法求解隶属函数U，对这些像素点进行划分。KNN分类算法能充分利用了像素点的领域信息，能较为准确的对这些像素点进行划分，较好好的抑制部分容积效应，同时降低噪声对不均匀性校正的影响。

(2)选取勒让德多项式作为基函数，大大减少曲线拟合所需参数个数，减少运算量。同时勒让德多项式在[-1，1]区间正交归一化，更加符合灰度不均匀场的分布情况，大大减少迭代次数。

(3)在求解灰度不均匀场模型模型参数时，假定迭代次数大于等于k′(通常设定为[5，10])次时，所得的灰度不均匀场模型模型参数的值即为最优解。这是由于本算法构造的能量函数是收敛的，当迭代次数大于一定值时，能量值(灰度不均匀场模型B(x，y)J(x，y)与原始MR图像I(x，y)之间的误差值)趋近于一个常量(误差最小)，经过大量反复的实验，对于MR图像来说，迭代次数定为[5，10]之间足以保证能量值最小(误差最小)，不需要在每次迭代后都 计算能量函数来判断收敛。大大降低了程序的复杂度和运算时间。

附图说明

图1是本发明总流程图。

图2是灰度区间初始划分图。

具体实施方式

本发明技术方案在实现时，首先使用Matlab语言编写仿真程序；然后使用MRI医学图像序列数据进行参数设定和程序优化处理；最后使用C++语言重写程序代码及交互界面框架，以提高程序性能。

Claims (3)

1.一种基于KNN分类算法的MR图像灰度不均匀性校正分割方法，包括以下步骤：

步骤1：建立灰度不均匀场模型；

所述灰度不均匀场模型为

其中：B(x，y)为灰度不均匀场，J(x，y)为未受灰度不均匀场污染的MR图像；且

其中g_k(x，y)为一组正交归一化的基函数，w_k为加权系数，k＝1，2，L，N，N为基函数的组数；c_i为MR图像中各个组织的灰度值，u_i(x，y)为隶属函数，i＝1，2，3或4；隶属函数u_i(x，y)相对于一个像素点(x，y)而言，定义为：

且 $\underset{i}{\overset{4}{\mathrm{\Σ}}}{u}_i(x,y)=1;$

步骤2：建立灰度不均匀场模型B(x，y)J(x，y)与原始MR图像I(x，y)之间的能量函数F：

$F=\∫{|I(x,y)-[\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{w}_k{g}_k(x,y)\left]\right[\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{c}_i{u}_i(x,y)\left]\right|}^2\mathrm{dxdy}$

步骤3：将步骤2所得能量函数F写成向量形式，即：

F(U，C，W)＝∫|I(x，y)-(W^TG(x，y))(C^TU(x，y))|²dxdy            (1)

其中：

B(x，y)＝W^TG(x，y)                        (2)

J(x，y)＝C^TU(x，y)                        (3)

C＝(c₁，c₂，c₃，c₄)^T，W＝(w₁，w₂，...w_N)^T

G(x，y)＝(g₁(x，y)，g₂(x，y)，...g_N(x，y))^T

U(x，y)＝(u₁(x，y)，u₂(x，y)，u₃(x，y)，u₄(x，y))^T

步骤4：任意确定灰度不均匀场模型B(x，y)J(x，y)的模型参数W和C的初始值W⁽¹⁾和C⁽¹⁾，由B(x，y)＝W^TG(x，y)计算出初始灰度不均匀场B⁽¹⁾(x，y)＝(W⁽¹⁾)^TG(x，y)，再由I(x，y)＝B(x，y)J(x，y)计算出初始未受灰度不均匀场污染的MR图像J⁽¹⁾(x，y)＝I(x，y)/B⁽¹⁾(x，y)，最后对初始未受灰度不均匀场污染的MR图像J⁽¹⁾(x，y)进行分类，求出初始隶属函数U⁽¹⁾(x，y)；

步骤5：利用下述迭代公式(4)(5)对灰度不均匀场模型模型参数C和W进行求解；

W^(k+1)＝(A^(k))^-1V^(k)，                    (4)

$c_i^{(k+1)}=\frac{\∫I(x,y)B^{\left(k\right)}(x,y)u_i^{\left(k\right)}(x,y)\mathrm{dxdy}}{\∫{\left(B^{\left(k\right)}\right)}^2(x,y)u_i^{\left(k\right)}(x,y)\mathrm{dxdy}},$ i＝1，2，3或4                (5)

其中：

B^(k)(x，y)＝(W^(k))^TG(x，y)

J^(k)(x，y)＝(C^(k))^TU^(k)(x，y)

V^(k)＝∫I(x，y)G(x，y)J^(k)(x，y)dxdy

A^(k)＝∫G(x，y)G^T(x，y)(J^(k))²(x，y)dxdy

步骤6：对隶属函数进行迭代求解；

首先由公式J^(k+1)(x，y)＝I(x，y)/B^(k)(x，y)计算出第k+1次迭代时未受灰度不均匀场污染的MR图像J^(k+1)(x，y)，接着按步骤4所述方法求出第k+1次迭代时的隶属函数U^(k+1)(x，y)；

步骤7：判断是否停止迭代；

当迭代次数达到设定值k′时停止迭代，将此时的灰度不均匀性模型模型参数分别记为记为W^*、U^*、c^* _i，即：

W^*＝W^(k′)                            (6)

${c^{*}}_i=c_i^{\left(k^{\′}\right)}$ i＝1，2，3，4                        (7)

$U^{*}=(U_1^{*}(x,y),u_2^{*}(x,y),u_3^{*}(x,y),u_4^{*}(x,y){)}^T=U^{\left(k\right)}---\left(8\right)$

步骤8：确定最终校正后的图像以及分割后的图像；

将步骤7中得到的最优灰度不均匀场模型模型参数c^* _i、U^*、W^*代入公式(1)计算出最终的灰度不均匀场B^*，然后由灰度不均匀场模型公式计算出最终校正后的MR图像J^*(x，y)，最后由公式

计算出最终分割后的图像：背景区域的MR图像I_BACK；脑脊液区域的MR图像I_CSF；脑灰质区域的MR图像I_GM；脑白质区域的MR图像I_WM。

2.根据权利要求1所述的基于KNN分类算法的MR图像灰度不均匀性校正分割方法，其特征在于，步骤4中对初始未受灰度不均匀场污染的MR图像J⁽¹⁾(x，y)进行分类，求出初始隶属函数U⁽¹⁾(x，y)的具体过程包括以下步骤：

步骤4-1：定义灰度值概率密度函数

其中p(q)表示灰度值为q的概率密度，n为MR图像的像素点总数，K(·)为0均值方差为1的高斯核函数，h为窗口宽度；

步骤4-2：在[0，50]像素宽度范围内确定一个合适的窗口宽度h，使得灰度值概率密度函数p(q)曲线的波峰数量为4；灰度值概率密度函数p(q)曲线的4个波峰对应着MR图像中脑白质、脑灰质、脑脊液和背景4个区域的最大概率处的灰度值；

步骤4-3：采用最大类间方差法确定脑脊液与灰质间的阈值TCG以及灰质与白质间的阈值TGW；具体方法是：a)从步骤4-1和步骤4-2确定的灰度值概率密度函数p(q)曲线中确定一个灰度区间[q_min，q_max]，其中q_min为[5，10]之间的一个灰度值；b)在灰度区间[q_min，q_max]内，找出一个灰度值q_med，使得在灰度区间[q_min，q_med]内，利用最大类间方差法求出脑脊液与灰质间的阈值TCG，在灰度区间[TCG，q_max]内，利用最大类间方差法求出白质与灰质间的阈值TGW，且TGW＝q_med；

步骤4-4：粗分类；

利用步骤4-3)所得的阈值TCG与TGW，加上一定的偏差δ，将灰度区间[q_min，q_max]分为4类区域：灰度值为[q_min，TCG-δ]的脑脊液区域Ω₁；灰度值为[TCG+δ，TGW-δ]的脑灰质区域Ω₂；灰度值为[TGW+δ，q_max]的脑白质区域Ω₃；灰度值为[TCG-δ，TCG+δ]或[TGW-δ，TGW+δ]的待定区域Ω₀；

步骤4-5：采用KNN分类算法对待定区域Ω₀中的像素点进行分类；

最终得到初始隶属函数：其中：

i＝1，2，3或4。

3.根据权利要求1所述的基于KNN分类算法的MR图像灰度不均匀性校正分割方法，其特征在于，步骤7中迭代次数设定值k′设为[10，20]之间。

Images