0 引言

核磁共振成像(Magnetic Resonance Imaging，MRI)技术是一种生物磁自旋成像技术，它可以直接做出各个方向的体层图像。同用X射线进行检查的电子计算机断层扫描(Computed Tomography，CT)图像相比，MRI以其非介入性、非损伤性、多角度多参数成像、高度的软组织分辨能力等优点，已成为脑疾病临床诊断的重要辅助手段^[1]。然而，由于噪声、场偏移效应和局部体效应等影响，获取的MRI图像常具有含噪高、对比度低、灰度不均匀、不同软组织之间或软组织与病灶之间边界模糊等特点，因此，MRI图像的精确分割是一大难题^[2]。

针对MRI图像分割的难题，国内外许多学者进行了深入的研究，提出了许多典型的分割算法，主要有基于阈值、基于模式识别、基于活动轮廓、基于马尔可夫随机场(Markov Random Field，MRF)和基于图切割这五个方面^[3]。其中，阈值分割以其简单快速的特点成为常用的一类MRI图像分割方法。目前，在已有的各类阈值选取方法^[4-6]中，以熵为目标函数的阈值分割方法备受关注。最初，一维最大Shannon熵法^[7]被用于图像阈值分割。此后，二维最大Shannon熵法^[8]及其粒子群优化算法^[9]也相继被提出，在提高算法的抗噪性的同时避免了计算量的膨胀。最大Shannon熵法选取最佳阈值时仅依照灰度级概率信息，忽视了类内灰度均匀性，导致一部分图像分割不准确；针对这一问题，文献[10]提出了Shannon灰度熵阈值选取方法，确保了分割后的图像中各类类内灰度的均匀性。上述最大熵法和最大灰度熵法建立在Shannon对数信息熵的基础上，对数运算存在零点无定义的问题, 而且仅依据一维灰度直方图，忽略了像素邻域信息，没有考虑各类类间的差异性。为此，文献[11]提出了二维倒数灰度熵阈值选取方法。由于MRI图像的成像质量较差，以上熵阈值方法用于MRI图像分割时存在以下问题：1) MRI图像具有含噪高、对比度低、灰度不均匀、不同软组织之间或软组织与病灶之间边界模糊等特点；而一维空间模型和简单的灰度-邻域平均灰度空间模型不能很好地反映图像的细节信息。2) MRI图像个体具有多样性和复杂性，目标对象可能为多个，且各目标或者背景所占灰度区间相异；单阈值分割无法体现组织间的层次性，分割得到的结果存在明显偏差。3)为减少计算量，基本粒子群优化算法因其实现简单、搜索速度较快而被广泛运用；然而，传统粒子群优化的熵阈值分割都以单个熵作为目标函数，不能同时考虑类间差异性与类内均匀性，且粒子群优化算法用于二维空间时计算量较大。

针对上述问题，本文提出了一种基于Curvelet变换和多目标粒子群(Multi-Objective Particle Swarm Optimization, MOPSO)的二维混合熵MRI图像多阈值分割方法。首先，为有效表示MRI图像中的纹理、边缘等目标细节信息，本文利用具有多分辨率、多方向性的Curvelet变换提取MRI图像的低频和高频信息，建立概貌-细节灰度级矩阵模型；其次，同时考虑目标与背景的类间差异性与类内均匀性，将提出的二维多阈值倒数熵和倒数灰度熵作为混沌多目标粒子群算法的目标函数，协同搜索得到最优的分割多阈值; 最后，为提高算法的求解速度，提出了二维倒数熵和倒数灰度熵多阈值选取的梯度算法。本文算法充分考虑了目标和背景的所有信息，多阈值选取能更好地适应MRI图像中目标的多样性和复杂性，递推算法大幅减少迭代过程中的重复计算，降低计算复杂度。最终实验结果表明本文算法用于分割各类MRI图像时，分割准确性最高。

1 概貌-细节灰度级矩阵模型

灰度-梯度直方图和灰度-邻域平均灰度空间模型是两种最常见的二维空间模型，这两种模型虽然利用了类内灰度均匀性，却忽视了高频细节信息和噪声的影响，鲁棒性较差。Curvelet变换是在小波变换的基础上发展起来的一种典型的多尺度几何分析方法，具有较好的线奇异性特点，可以保留更多的图像边缘信息。为此，本文提出基于Curvelet变换的概貌-细节灰度级矩阵模型，利用Curvelet变换具有较好的线奇异性、多分辨率和多方向性特点实现图像高频细节的提取和噪声的抑制。对于Curvelet变换而言，随着尺度的增加，纹理细节等有用信息和噪声衰减程度不同，故其对边缘、纹理及噪声有更强的区分能力，能够在有效抑制图像噪声的同时更好地保留目标的边缘、纹理信息。

1.1 Curvelet变换

为了有效地表示图像的边缘等特征, Candes和Donoho提出了第一代Curvelet变换，其构造思想是通过足够小的分块将曲线近似为每一个分块中的直线来看待，然后利用脊波(Ridgelet)变换理论来进行分析。该过程包括子带分解、平滑分块、正规化和Ridgelet分析等步骤，实现复杂，数据冗余量也比较大。Candes和Donoho在2004年提出了更为简洁的第二代Curvelet变换，其基本思想是将频域划分为楔形区域，再对这些楔形区域采用局部傅里叶基变换来实现。由于与脊波变换没有关系，第二代Curvelet变换更为简单，速度更快。实现快速离散的Curvelet变换有两种方法^[12]:一种是基于频域特殊采样的卷绕规则Wrapping方法; 另一种是基于非均匀采样的快速傅里叶变换(Unequally-Spaced Fast Fourier Transform, USFFT)方法。

1.2 基于Curvelet变换的概貌-细节灰度级矩阵模型

本文利用Curvelet变换中的Wrapping方法构建更为有效的MRI图像概貌-细节灰度级矩阵模型。具体步骤如下：

假设I是一幅待分割的MRI图像，大小为M×N，灰度级数为L-1。首先对其进行Curvelet变换，从低频到高频得到一系列不同尺度不同方向的子带；然后保留Curvelet低频子带系数，将高频子带系数置零，进行逆变换得到概貌图像f；再保留高频子带系数，将低频子带系数置零，进行逆变换得到细节图像g。将概貌图像和细节图像取绝对值后映射到；用r(i, j)表示概貌图像f的灰度级是i且细节图像g的灰度级是j出现的频数。定义联合概率p(i, j)=r(i, j)/(M×N)，则{p(i, j)}为该MRI图像的概貌一细节灰度级矩阵，如图 1(a)所示。

本文使用多阈值分割，假设含有n个目标，则概貌域选取n个阈值为t₁, t₂, …, t_n(0≤t₁ < t₂ < … < t_n-1 < t_n < L-1)，细节域因集中分布在较小灰度级，故取单阈值s₁。概貌-细节灰度级矩阵分割成图 1(b)所示，其中阴影部分表示目标和背景，其余部分表示背景和目标的边缘和噪声。

2 二维倒数熵和倒数灰度熵的多阈值选取

灰度熵选取的是使目标或背景内部的灰度分布尽可能均匀的最优阈值；最大倒数熵准则强调的是类间的差异性，应用在阈值化分割中就是搜索使目标和背景间的差异最大的最优阈值。

2.1 二维倒数熵的多阈值选取 2.1.1 倒数熵

若假设p_i表示事件i(i=0, 1, …, l-1)发生的概率$\sum\limits_{i = 0}^{l-1} {{p_i} = 1, 0 \le {p_i} \le 1} $，l为信源符号个数，自信息量为ΔI(p_i)=1/(1+p_i)，所有事件自信息量的期望即为一维倒数熵：

$ H = E(\Delta I) = \sum\limits_{i = 0}^{l - 1} {\frac{{{p_i}}}{{1 + {p_i}}}} $

(1)

一维倒数熵式(1)推广至二维倒数熵, 传统的二维空间为：灰度级-邻域平均灰度级，令p_ij表示该二维空间中二元对(i, j)发生的联合概率，${P_{ts}} = \sum\limits_{j = 0}^s {\sum\limits_{i = 0}^t {{p_{ij}}} } $, 则目标类和背景类的二维倒数熵H_o(t, s)和H_b(t, s)可表示为：

$ {H_o}(t,s) = \sum\limits_{j = 0}^s {\sum\limits_{i = 0}^t {\frac{{{p_{ij}}/{P_{ts}}}}{{1 + {p_{ij}}/{P_{ts}}}}} } = \sum\limits_{j = 0}^s {\sum\limits_{i = 0}^t {\frac{{{p_{ij}}}}{{{p_{ij}} + {P_{ts}}}}} } $

(2)

$ \begin{array}{c} {H_b}(t,s) = \sum\limits_{j = 0}^s {\sum\limits_{i = t + 1}^{L - 1} {\frac{{{p_{ij}}/(1 - {P_{ts}})}}{{1 + {p_{ij}}/(1 - {P_{ts}})}}} } = \\ \sum\limits_{j = 0}^s {\sum\limits_{i = t + 1}^{L - 1} {\frac{{{p_{ij}}}}{{{p_{ij}} + 1 - {P_{ts}}}}} } \end{array} $

(3)

总倒数熵为：

$ H(t,s) = {H_o}(t,s) + {H_b}(t,s) $

(4)

二维最大倒数熵选取方法的准则函数为：

$ ({t^*},{s^*}) = \arg \mathop {\max \max }\limits_{0 \le t \le L - 1,0 \le s \le L - 1} \left\{ {H(t,s)} \right\} $

(5)

当总倒数熵取最大值时，所对应的阈值向量(t^*, s^*)为二维最大倒数熵分割方法的最优阈值。

2.1.2 基于概貌-细节矩阵的二维倒数熵的多阈值选取

如图 1(b)所示，假设含有n个目标，则概貌域相应选取n个阈值为t₁, t₂, …, t_n，其中0≤t₁ < t₂ < … < t_n-1 < t_n < L-1，细节域因集中分布在较小灰度级，故取单阈值s₁。

本文基于概貌-细节矩阵模型将式(4)推广得到二维多阈值倒数熵H_R(t₁, t₂, …, t_n, s₁)：

$ \begin{array}{l} {H_R}({t_1},{t_2}, \cdots ,{t_n},{s_1}) = \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {\frac{{{p_{_{ij}}}}}{{{p_{_{ij}}} + \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {{p_{_{ij}}}} } }}} } + \\ \;\;\;\;\;\;\;\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_1} + 1}^{{t_2}} {\frac{{{p_{_{ij}}}}}{{{p_{_{ij}}} + \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_1} + 1}^{{t_2}} {{p_{_{ij}}}} } }}} } + \cdots + \\ \;\;\;\;\;\;\;\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {\frac{{{p_{_{ij}}}}}{{{p_{_{ij}}} + \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {{p_{_{ij}}}} } }}} } \end{array} $

(6)

最大倒数熵多阈值选取公式：

$ \begin{array}{l} (t_{R1}^*,t_{R2}^*, \cdots ,t_{Rn}^*,s_{R1}^*) = \\ \;\;\;\arg \mathop {\max }\limits_{0 \le {t_1} < {t_2} < \cdots < {t_{n - 1}} < {t_n} < L - 1} \mathop {\max }\limits_{0 \le {s_1} < L - 1} \{ {H_R}({t_1},{t_2}, \cdots ,{t_n},{s_1})\} \end{array} $

(7)

当总倒数熵H_R(t₁, t₂, …, t_n, s₁)取最大值时，所对应的n个灰度级阈值(t_R1^*, t_R2^*, …, t_Rn^*, s_R1^*)即为最大倒数熵多阈值分割方法的最佳阈值。

2.2 二维倒数灰度熵的多阈值选取 2.2.1 倒数灰度熵

若令${\mu _{oi}}(t, s) = \sum\limits_{i = 0}^t {\sum\limits_{j = 0}^s {h(i, j)} } i$, ${\mu _{oj}}(t, s) = \sum\limits_{i = 0}^t {\sum\limits_{j = 0}^s {h(i, j)} } j$, ${\mu _{bi}}(t, s) = \sum\limits_{i = t + 1}^{L-1} {\sum\limits_{j = s + 1}^{L-1} {h(i, j)} } $, ${\mu _{bj}}(t, s) = \sum\limits_{i = t + 1}^{L-1} {\sum\limits_{j = s + 1}^{L-1} {h(i, j)} } $，其中h(i, j)表示图像中灰度级-邻域平均灰度级二元对(i, j)出现的频数，则目标类的二维倒数灰度熵H_o(t, s)定义为：

$ \begin{array}{c} {\boldsymbol{H}_o}(t,s) = {({H_{oi}}(t,s),{H_{oj}}(t,s))^{\rm{T}}} = \\ {(\sum\limits_{j = 0}^s {\sum\limits_{i = 0}^t {\frac{{h(i,j)i}}{{{\mu _{oi}}(t,s) + i}}} } ,\sum\limits_{j = 0}^s {\sum\limits_{i = 0}^t {\frac{{h(i,j)j}}{{{\mu _{oj}}(t,s) + j}}} } )^{\rm{T}}} \end{array} $

(8)

背景类的二维倒数灰度熵H_b(t, s)定义为：

$ \begin{array}{c} {\boldsymbol{H}_b}(t,s) = {({H_{bi}}(t,s),{H_{bj}}(t,s))^{\rm{T}}} = \\ {(\sum\limits_{j = s + 1}^{L - 1} {\sum\limits_{i = t + 1}^{L - 1} {\frac{{h(i,j)i}}{{{\mu _{bi}}(t,s) + i}}} } ,\sum\limits_{j = s + 1}^{L - 1} {\sum\limits_{i = t + 1}^{L - 1} {\frac{{h(i,j)j}}{{{\mu _{bj}}(t,s) + j}}} } )^{\rm{T}}} \end{array} $

(9)

目标类和背景类的总二维倒数灰度熵H(t, s)为：

$ \boldsymbol{H}(t,s) = {\boldsymbol{H}_o}(t,s) + {\boldsymbol{H}_b}(t,s) $

(10)

定义：

$ \begin{array}{c} \xi (t,s) = \sum\limits_{j = 0}^s {\sum\limits_{i = 0}^t {\frac{{h(i,j)i}}{{{\mu _{oi}}(t,s) + i}}} } + \sum\limits_{j = s + 1}^{L - 1} {\sum\limits_{i = t + 1}^{L - 1} {\frac{{h(i,j)i}}{{{\mu _{bi}}(t,s) + i}} + } } \\ \sum\limits_{j = 0}^s {\sum\limits_{i = 0}^t {\frac{{h(i,j)j}}{{{\mu _{oj}}(t,s) + j}}} } + \sum\limits_{j = s + 1}^{L - 1} {\sum\limits_{i = t + 1}^{L - 1} {\frac{{h(i,j)j}}{{{\mu _{bj}}(t,s) + j}}} } \end{array} $

(11)

式中:ξ(t, s)为倒数灰度熵的阈值选取准则函数，当函数ξ(t, s)达到最大值时，获得最佳阈值向量(t^*, s^*)为:

$ ({t^*},{s^*}) = \arg \mathop {\max \max }\limits_{0 \le t \le L - 1,0 \le s \le L - 1} \left\{ {\xi (t,s)} \right\} $

(12)

2.2.2 基于概貌-细节矩阵的二维倒数灰度熵的多阈值选取

假设含有n个目标，则概貌域相应选取n个阈值为t₁, t₂, …, t_n(0≤t₁ < t₂ < …t_n-1 < t_n < L-1)，细节域因集中分布在较小灰度级，故取单阈值s₁。

本文在概貌-细节矩阵模型中将式(11)推广得二维多阈值灰度熵H_G(t₁, t₂, …, t_n, s₁):

$ \begin{array}{l} {H_G}({t_1},{t_2}, \cdots ,{t_n},{s_1}) = \\ \;\;\;\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {\frac{{h(i,j)i}}{{\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {h(i,j)i} } + i}}} } + \\ \;\;\;\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {\frac{{h(i,j)j}}{{\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {h(i,j)j} } + j}}} } + \cdots + \\ \;\;\;\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{n - 1}} + 1}^{{t_n}} {\frac{{h(i,j)i}}{{\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{n - 1}} + 1}^{{t_n}} {h(i,j)i} } + i}}} } + \\ \;\;\;\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{n - 1}} + 1}^{{t_n}} {\frac{{h(i,j)j}}{{\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{n - 1}} + 1}^{{t_n}} {h(i,j)j} } + j}}} } + \\ \;\;\;\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {\frac{{h(i,j)i}}{{\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {h(i,j)i} } + i}}} } + \\ \;\;\;\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {\frac{{h(i,j)j}}{{\sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {h(i,j)j} } + j}}} } \end{array} $

(13)

获得最佳二维倒数灰度多阈值向量(t_G1^*, t_G2^*, …, t_Gn^*, s₁^*)为:

$ \begin{array}{l} (t_{G1}^*,t_{G2}^*, \cdots ,t_{Gn}^*,s_1^*) = \\ \;\;\;\arg \mathop {\max }\limits_{0 \le {t_1} < {t_2} < \cdots < {t_{n - 1}} < {t_n} < L - 1} \mathop {\max }\limits_{0 \le {s_1} < L - 1} \{ {H_G}({t_1},{t_2}, \cdots ,{t_n},{s_1})\} \end{array} $

(14)

3 多目标粒子群混合熵多阈值医学分割 3.1 方法流程

本文算法的流程如图 2所示。

3.2 多目标粒子群混合熵的多阈值分割的具体步骤

本文依据多目标规划理论^[13]将倒数熵和倒数灰度熵同时作为多阈值搜索的优化准则，这里把这两种熵的结合定义为混合熵，即将倒数熵和倒数灰度熵作为两个目标函数协同优化粒子群。

算法步骤：

1) 混沌粒子群初始化：以改进Tent序列映射到概貌-细节灰度级模型内，从中选择N个二维向量作为粒子的初始位置X[i]，速度V[i]在[-Vm, Vm]上随机产生。

2) 将H_R(t₁, t₂, …, t_n, s₁)和H_G(t₁, t₂, …, t_n, s₁)作为两个目标函数，分别计算每个粒子的适应度值：

$ \begin{array}{l} Fit1[i] = {H_R}\left( {\boldsymbol{X}[i]} \right);i = 1, 2, \cdots, N\\ Fit2[i] = {H_G}\left( {\boldsymbol{X}[i]} \right);i = 1, 2, \cdots, N \end{array} $

3) 在目标函数H_R(t₁, t₂, …, t_n, s₁)和H_G(t₁, t₂, …, t_n, s₁)下分别对每个粒子求个体极值：

$ \begin{array}{l} \boldsymbol{pBest}[1, i] \leftarrow {H_R}\left( {\boldsymbol{X}[i]} \right);i = 1, 2, \cdots, N\\ \boldsymbol{pBest}[2, i] \leftarrow {H_G}\left( {\boldsymbol{X}[i]} \right);i = 1, 2, \cdots, N \end{array} $

4) 对目标函数H_R(t₁, t₂, …, t_n, s₁)和H_G(t₁, t₂, …, t_n, s₁)分别求全部极值：

$ \begin{array}{l} \boldsymbol{pBest}[1] \leftarrow {H_R}\left( {\boldsymbol{X}[i]} \right)\\ \boldsymbol{pBest}[2] \leftarrow {H_G}\left( {\boldsymbol{X}[i]} \right) \end{array} $

5) 计算全局矢量的均值gBest和距离dgBest：

$ \begin{array}{l} \boldsymbol{gBest} = Average(\boldsymbol{gBest}[1], \boldsymbol{gBest}[2]\\ dgBest = Distance(\boldsymbol{gBest}[1], \boldsymbol{gBest}[2]) \end{array} $

6) 计算每个粒子pBest[1, i]和pBest[2, i]之间的距离dpBest[i]：

$ \begin{array}{l} dpBest[i] = Dis\tan ce(\boldsymbol{pBest}[1, i], \boldsymbol{pBest}[2, i]);\\ \;\;\;\;\;\;\;\;i = 1, 2, \cdots, N \end{array} $

7) 对每个粒子计算更新速度V[i]和位置X[i]时，个体极值pBest[i]为：

For i=1 to N

 If (dpBest[i] < dgBest)

  pBest[i]=RandSelect(pBest[1, i], pBest[2, i]);

 else

  pBest[i]=Average(pBest[1, i], pBest[2, i]);

Next i

8) 用步骤5)和步骤7)所得的gBest和pBest[i]更新每个粒子速度V[i]和位置X[i]。

9) 若满足至最大迭代次数，停止迭代，否则返回2)(若优化停滞步数大于5时，进行Tent混沌变异，更新粒子的速度和位置)。

10) 迭代结束，得到最优阈值，使用最优阈值对图像进行多阈值分割。

3.3 二维灰度熵和倒数灰度熵多阈值选取的梯度算法

由式(6)可知，计算复杂度为O(L⁴)，计算量大，为提高运算效率、降低计算冗余度，可采用递推方式计算式(6)中的中间变量，计算复杂度为O(L²)。

令w(0, 0)=p₀₀；w(t, 0)=w(t-1, 0)+p_t0；w(0, s)=w(0, s-1)+p_0s；w(t, s)=w(t-1, s)+w(t, s-1)-w(t-1, s-1)+p_ts；则

$ \begin{array}{l} v({t_1}, 0, {s_1}) = w({t_1}, {s_1}) = \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {{p_{ij}}} } \\ v({t_2}, {t_1}, {s_1}) = w({t_2}, {s_1})-w({t_1}, {s_1}) = \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_1} + 1}^{{t_2}} {{p_{ij}}} } \\ v({t_n}, {t_{n-1}}, {s_1}) = w({t_n}, {s_1})-w({t_{n - 1}}, {s_1}) = \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{n - 1}} + 1}^{{t_n}} {{p_{ij}}} } \\ v(L - 1, {t_n}, {s_1}) = w(L - 1, {s_1}) - w({t_n}, {s_1}) = \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {{p_{ij}}} } \end{array} $

则式(6)二维多阈值倒数熵的递推形式为：

$ \begin{array}{l} {H_R}({t_1},{t_2}, \cdots ,{t_n},{s_1}) = \\ \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {\frac{{{p_{ij}}}}{{{p_{ij}} + v({t_1},{s_1})}} + } } \\ \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_1} + 1}^{{t_2}} {\frac{{{p_{ij}}}}{{{p_{ij}} + v({t_2},{t_1},{s_1})}} + } } \cdots + \\ \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{n - 1}} + 1}^{{t_n}} {\frac{{{p_{ij}}}}{{{p_{ij}} + v({t_n},{t_{n - 1}},{s_1})}}} } + \\ \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {\frac{{{p_{ij}}}}{{{p_{ij}} + v(L - 1,{t_n},{s_1})}}} } \end{array} $

(15)

同理，采用递推方式计算式(13)二维多阈值倒数灰度熵中的中间变量，则式(13)的递推形式为：

$ \begin{array}{l} {H_G}({t_1},{t_2}, \cdots ,{t_n},{s_1}) = \\ \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {\frac{{h(i,j)i}}{{{m_i}({t_1},{s_1}) + i}}} } + \\ \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = 0}^{{t_1}} {\frac{{h(i,j)j}}{{{m_j}({t_1},{s_1}) + j}}} } + \cdots + \\ \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{n - 1}} + 1}^{{t_n}} {\frac{{h(i,j)i}}{{{m_i}({t_n},{s_1}) + i}}} } \\ + \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{n - 1}} + 1}^{{t_n}} {\frac{{h(i,j)j}}{{{m_j}({t_n},{s_1}) + j}}} } + \\ \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {\frac{{h(i,j)i}}{{{m_i}(L - 1,{s_1}) + i}}} } + \\ \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_n} + 1}^{L - 1} {\frac{{h(i,j)j}}{{{m_j}(L - 1,{s_1}) + j}}} } \end{array} $

(16)

其中${m_i}({t_k}, {s_1}) = \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{k-1}} + 1}^{{t_k}} {h(i, j)} i} $, ${m_j}({t_k}, {s_1}) = \sum\limits_{j = 0}^{{s_1}} {\sum\limits_{i = {t_{k-1}} + 1}^{{t_k}} {h(i, j)} j} $。

4 实验结果与分析

为了验证本文算法的可行性和有效性，对大量MRI进行实验，其中选取MRI的脑部图像序列、腰椎病变图像、头颅切面图像以及含噪模拟图像的分割结果进行分析。与二维tsallis熵分割法^[14]、自动细菌觅食分割法(Adaptive Bacterial Foraging，ABF)^[15]和改进的Otsu分割法^[16]作对比实验，先对一组MRI脑部图像进行双阈值和三阈值分割，结果表明阈值数越多分割结果越好。为此，再对其他类型的MRI图像，如腰椎病变图像、头颅切面图像和含噪模拟图像进行三阈值分割。最后，分别采用差异性和均匀性的综合测度^[17]，分割准确性^[18]对分割结果进行定量分析，综合测度值越大表明分割效果越好。

实验环境为Intel Core i5处理器，4 GB内存，CPU为2.3 GHz，Windows 7系统，采用Matble2012编译语言实现。

为了公平起见，同时保证能使各自的算法优化效果最佳，且尽量达到稳定收敛，对各算法的参数设定如下：二维tsallis熵分割法：量子克隆进化法的种群规模为50，克隆控制系数为κ=16，混沌扰动系数λ₀=1.8，最大迭代次数(Maxgen)设置时，为了满足不同阈值选择的需要，Maxgen将随着阈值个数的增加而增加，设D为阈值个数，则Maxgen表示为Maxgen=2*D²+(D-1)*20+20，以达到自动调整Maxgen的目的；改进的Oust算法：因脑部图像相邻灰度差别较小，则需适当加大梯度影响，梯度权值a₁=0.3，距离权值a₂=0.1，灰度权值a₃=0.6，CHA=60；ABF算法：细菌个数20，最大迭代次数为Maxgen/2，趋向性操作次数为10，最大前进步数为4，复制和驱散操作次数都为2，细菌被驱散的概率为0.02；本文算法：MOPSO算法的种群规模和最大迭代次数都与量子克隆进化法一致，学习因子c₁和c₂均为2.1，惯性因子为0.6299。

4.1 分割结果的视觉效果评价

如图 4~9所示，本文先用一组MRI脑部图像与二维tsallis熵分割法、ABF分割法和改进的Otsu分割法作双阈值和三阈值分割的对比实验。由图 4~9可知，本文算法较其他算法分割效果更精准，脑沟壑更连通，局部散点更少而且边界清晰。这是由于本文算法充分发挥了Curvelet变换局部细节提取能力强的优势，弥补了PSO算法局部搜索能力不足的缺陷，又充分发挥了概貌-细节框架的全局把握能力，解决了ABF算法全局搜索能力较差的问题。对于复杂的多阈值分割，二维tsallis熵分割法丢失部分边缘信息，脑灰质无法分割出来。而改进Otsu算法不能自动获得梯度权值、距离权值和灰度权值的比值，需人为取经验值，无法做到自适应。无论是在全局收敛的可靠性上还是在自适应方面，本文算法都明显优于ABF算法和改进的Otsu算法。实验结果说明了本文算法的有效性。

如图 11~13所示，因由脑部图像分割可知三阈值优于双阈值，则本文再用其他类型的MRI与二维tsallis熵分割法、ABF分割法和改进的Otsu分割法作三阈值分割对比实验。图 11为MRI腰椎病变图像，三种对比算法均可在一定程度分割出右侧病变部分，但噪声明显，且无法分割出腰椎骨节部分。图 12为MRI头颅切面图像，由分割结果可知，三种对比算法都无法得到完整的脑灰质部分：而二维tsallis熵分割法分割的脑白质部分过细，部分边缘信息丢失；本文方法得到的脑灰质和脑白质较为完整、准确，细节信息丰富。图 13为医学加噪模拟图像，二维tsallis熵分割算法因图中含有噪声而分割失败，ABF分割法和改进的Otsu分割法无法克服噪声的影响，噪声干扰严重。本文方法对含噪的图像分割具有较好的适应性，分割结果准确。

二维tsallis熵分割算法分割效果不佳是因其存在无定义值和零值的缺陷，造成了分割后图像信息的丢失；ABF的全局搜索能力较差，导致类间差异性无法顾及；改进的Otsu算法先考虑类间的差异性，再利用梯度权值、距离权值和灰度权值的比值来考虑类内的均匀性，但该比值无法自适应且会影响差异性。本文方法将倒数熵和倒数灰度熵共同作为多目标粒子群的判断准则，不仅考虑反映了图像中目标类与背景类间的灰度级差异，而且同时考虑了均匀性。再者，使用高低频分层的概貌-细节空间，可以捕获目标的细节信息，滤除噪声，从而得到较为精确的边界分割结果。因此，本文方法能够适用于结构和层次复杂的含噪MRI图像，对肿瘤等异常部位的检测定位较为准确，实验结果证实了这一点。

4.2 分割结果的客观评价

为进一步准确评价算法的性能，本文分别采用差异性和均匀性的综合测度、分割准确性对分割结果进行定量分析。

1) 差异性和均匀性的综合测度。为了对比不同方法对类间差异性和类内均匀性的考虑情况，分别用以下几个参数来衡量。

①类间差异性测度(GC)用区域间对比度来衡量，定义为：

$ GC = \frac{{\left| {{f_1} - {f_3}} \right|}}{{{f_1} + {f_2} + {f_3}}} $

(17)

式中:f₁、f₂和f₃分别为原图中对应于分割图像的3个区域的平均灰度值。

②各区域类内均匀性测度(UM)定义为：

$ UM = 1 - \frac{1}{C}\sum\limits_i {\left\{ {{{\sum\limits_{(x,y) \in {R_i}} {\left[ {f(x,y) - \frac{1}{{{A_i}}}\sum\limits_{(x,y) \in {R_i}} {f(x,y)} } \right]} }^2}} \right\}} $

(18)

式中：R_i为分割图的第i个区域；A_i表示第i个区域像素个数；C为归一化参数；f(x, y)为(x, y)点的像素值。

③综合测度M为:

$ M = \alpha * GC + \beta * UM $

(19)

式中:α和β分别为GC和UM的权值系数，α+β=1。综合测度的值越大，分割效果越好。本文依据区域内部含有细节信息的多少，假定设综合测度α=0.4和β=0.6。

表 1为几种方法的最佳分割阈值。

从表 1可知，因脑部图像局部灰度差较小，所以4种算法所得阈值相差不大，但还是可以看出tsallis熵阈值与其他方法相比更高，ABF次之，这是因为这两种算法都为单一算法，考虑不全面。改进的Otsu算法与本文算法分割结果较近，但还是有较小波动，这与图 3~4的分割的主观视觉评价相符。同时也可看出这些对比算法没有同时考虑类间差异性和类内均匀性，而本文算法满足该条件，获得了比较合理的MRI图像分割阈值。

表 2为各算法的差异性、均匀性和综合测度，加粗部分为几种算法分割同一图片时各测度的最优值。由表 2可知，各算法中本文算法的差异性测度最大，改进的Otsu算法次之，这是由于改进的Otsu算法是虽以类间方差最大作为寻优准则，与差异性的测量原理相一致，但还是会受均匀性优化算法的影响。各算法中，本文算法的均匀性最大，ABF算法次之，ABF算法的局部搜索能力强，因而均匀性测度较高。而本文以灰度倒数熵和灰度熵作为寻优准则，类间差异性和均匀性可同时最优；由各方法的综合测度值可知，本文算法综合测度值最高，说明分割效果最好。

2) 分割准确性。为了进一步定量对比不同方法对MRI图像的分割准确性，通过测量面积将不同方法分割含噪模拟图(图 13)结果与去噪标准图进行对比。具体采用准确率R和检全率P加以衡量，其分别定义为：

$ R = ({N_{o1}} \cap {N_{s1}} + {N_{o2}} \cap {N_{s2}})/({N_{s1}} + {N_{s2}}) $

(20)

$ P = ({N_{o1}} \cap {N_{s1}} + {N_{o2}} \cap {N_{s2}})/({N_{o1}} + {N_{o2}}) $

(21)

其中：N_o1、N_o2分别表示标准分割的目标1、2的像素数; N_s1、N_s2为不同方法分割出的目标1、2的像素数。

为了便于分析，将R和P合并为一个参量F进行度量：

$ F = PR/[\alpha P + (1 - \alpha )R] $

(22)

其中：α为R和P的权衡系数，这里α取值为0.5。可以看出，F随着R和P的增大而增大，当R和P均较大，即F较大时，表示分割结果越接近标准结果，分割准确性越高。对比其他三种方法(图 14)，可以看出，本文方法具有较高准确率R和检全率P，参数F均大于80%，分割准确性最高。

5 结语

提出了基于Curvelet和MOPSO的混合熵MRI图像分割算法。首先，基于Curvelet变换构造高低频分层的概貌-细节空间，使得本文方法能有效提取高频细节信息，精确分割图像细节部分的效果尤其显著。再者，本文方法将改进的二维递推多阈值倒数熵和倒数灰度熵共同作为多目标粒子群的判断准则，弥补了其他熵无定义值和零值的缺陷，充分考虑所有图像信息，用递推算法降低熵的计算复杂度。采用多目标粒子群算法不仅能大幅减少计算量，而且同时满足类间差异性与类内均匀性。最后，与二维tsallis熵分割法、ABF分割法和改进的Otsu分割法实验比较可知，本文算法对灰度不均、低对比度和含噪的MRI图像分割效果最好，目标清晰准确，误分割最少。