论文部分内容阅读
【摘 要】本文提出一种改进的使用主动轮廓模型分割脑部肿瘤图像的方法。采用改进的Greedy法，对轮廓线蛇点的邻域提出了新的搜索方法，并给出判定的准则。同时使用遗传算法处理图像的二维直方图，得到阈值后对图像进行二值化处理，并使用静电场模型对二值图像的梯度场计算外力。对比试验表明了该算法的有效性。
【关键词】主动轮廓模型 贪婪法 二維直方图 静电场 分割
一、引言
由于医学图像的数据量大，目标区域的解剖结构复杂多变，并且图像中伪影和各种噪声的影响使得目标的边缘往往模糊而不连续。这使得准确分割医学图像非常困难，由于问题的复杂性和重要性，如何构造有效的分割算法是一个很有挑战的工作。
被称为Snake的主动轮廓模型（Active Contour Model）由Kass等人提出，并被广泛的应用在一系列的图像视觉问题中。经过初始轮廓的估计后，Snake能够在轮廓线内力和图像外力的作用下自动的收敛能量的极小状态，最终运动到目标的边缘附近，从而完成图像的分割。近年来学者提出种种新方法，对轮廓线蛇点邻域的搜索和Snake的内力和外力的计算进行了各种改进，试图在各种约束条件下找到一个平衡点[1-2][4][5]，本文使用遗传算法首先对脑部MRI图像的二维直方图进行处理以提取的阈值，对图像进行二值处理后，求得二值图像的梯度边缘，并使用静电力模型对初始轮廓线内部的区域计算图像的外力，从而使轮廓线收敛，达到分割脑部MRI图像的目的。
二、使用主动轮廓模型进行图像分割
（一）主动轮廓模型的描述
参数化的Snake是一条动态的轮廓线，其参数方程为，整个轮廓线的能量定义为：
其中，，分别Snake的内力强度，外力强度，和人为交互力的强度。如果不需要人为干预Snake的行为的话，通常可以取。轮廓线的内力强度为：
式中的第一项为弹性势能，当轮廓线长度缩短时，该项积分取较小值，第二项为刚性能量，在轮廓线趋于平滑时，该项的积分取较小值，参数分别是弹性能量和刚性能量的控制系数，当两参数分别为零时，允许蛇点在此处的弹性能量和刚性能量取一个极大值，即允许轮廓线出现不连续，或曲率较大。所以，把这两个参数作为时变参数的话，可以近一步的影响蛇的行为。一般处理中为了简单，通常可以令这两个参数为常量。
（二）图像的外力
图像的外力一般可以取下面的形式：
即图像灰度和高斯函数卷积的梯度模值平方的负数，其中为控制系数。简化的外力还可以取公式（1）中高斯函数取极值的情况
但基于梯度的外力往往是一种短程力，难以捕获远处的轮廓线，并且轮廓线在这种外力的作用下难以收敛到目标区域的凹陷处。
三、基于静电场模型的Snake算法
（一）图像的二值化处理
本文使用Snake模型计算脑部MRI图像外力前，首先对图像进行二值化处理。进行二值化处理可以利用图像的直方图。图1（b）是图1左图脑部水肿图像的直方图，由于直方图的双峰间的峰谷往往过于平坦，使用一维的Otsu方法往往难以得到理想的灰度阈值。
为了得到较优的灰度阈值，可以使用二维直方图，本文使用基于原始图和其八邻域平滑图得到的二维直方图：设图像的大小为，在象素处的灰度为，象素点的八邻域均值为：
设原始图和平滑图为，联合象素满足的概率为，其中，则以为自变量，函数为联合图像的二维灰度直方图。图2为图1（a）的二维直方图的灰度图。
其中高亮度的区域对应二维直方图中取值大的部分，由图可以见，灰度二维直方图取值较大的部分集中在xy平面的对角线附近。 取x轴阈值s，和y轴阈值t，可以把二维直方图分为如下四个部分。
其中和区对应图像中的噪声和图边缘部分，概率较小，从而阈值划分问题转化为求阈值矢量（s，t），使得熵函数
对于公式（3）中最优参数，可以采用各种优化算法进行计算 [7][8]，例如可以使用遗传算法计算得到该类图像灰度阈值，由脑部MRI图像的先验知识可知，这类图像的阈值分布在100左右，可以尽一步缩小优化算法的搜索空间，加快算法的收敛速度。
本文使用的遗传算法采用直接对二维阈值矢量进行编码，设阈值分布以100为中心，幅度为±32的区间内，即[69，132]，对二维阈值矢量可以采用两基因（gene）的染色体编码，每个基因采用采用6比特编码，每个染色体（chromosome）为12比特，前面的6个比特为第一个基因，为二维直方图第一个坐标阈值的编码。后面第二个基因代表第二个坐标阈值。此时搜索空间中的点为61X61。染色体采用公式（2）作为适应度函数。
基本步骤：（a）初始染色体种群，种群数量为N=40。（b）对染色体解码，并计算适应度。（c）淘汰适应度最差的20个染色体，并对最优的20个染色体进行交叉，产生新的染色体，使种群恢复到原来的规模，并采用0.1概率对现有染色体进行变异。（d）如果近10代染色体种群的最优适应值保持不变，迭代结束，否则转（b）。
由于使用了图像的先验知识，大大减少了计算的工作量。否则要假设阈值分布在[0，255]，寻优空间的大小为255X255。实际上满足其他条件（例如图像的尺寸，和脑部图像与背景比例等）的脑部MRI图像的阈值将分布在更窄取值范围内，阈值矢量的搜索空间可以近一步的缩小。搜索空间足够小时，还可以使用全局搜索的方法，从而得到最佳的阈值。图3为使用遗传算法计算得到二维最大熵阈值，并使用该阈值对图1中的左图进行二值化的结果。
（二）基于静电场模型的图像外力的计算
由于主动轮廓模型传统外力的限制，学者提出了种种改进的外力模型[4][5]，本文使用文献[5]的思想，使用二值图像的梯度图在初始轮廓线内部区域的静电场力作为主动轮廓模型的图像外力。 设图像中每个象素的位置放置负电荷
，M是图像象素的个数，是象素处的灰度值，每个蛇点带单位正电荷，由于图像边缘的电场强度较大，蛇点将在电场的作用下向物体的边缘进行移动。同时，由于蛇点之间斥力的作用，蛇点将趋于均匀分布。设蛇点的数量为N，则轮廓线在静电场中的势能为：
其中分别为标号为的蛇点的矢量坐标，为标号为的象素的矢量坐标。带单位正电荷的蛇点在静电场中所受的引力为：
所受其他蛇点的斥力为：
为了突出二值图像边缘的影响，对图3进行求梯度的运算，并以梯度模值表征图像的边缘。图4（a）为是梯度运算的结果。图4（b）为使用静电场模型对图4（a）的目标区域（矩形框中）进行计算后得到的图像外力的结果。
由图4（b）可以看出，在水肿区域的边缘处形成较强的外力，即使目标区域的凹陷区也有外力的存在。
四、蛇点邻域的快速搜索算法
本文采用了改进的Greedy法在邻域中移动蛇点。为了使蛇点间的距离更加均匀，使目标蛇点在相邻蛇点的中垂线上移动是合理的。当前蛇点为，相邻蛇点为，两个相邻蛇点的中垂线为，如果和的邻域相交，则相交的邻域点为贪婪法的搜索点（图5中阴影部分），否则，把当前蛇点移动到离中垂线距离最近的邻域点。判断邻域点和是否相交，可以采用下面的准则：
设为定义好的一个比较小的值，设的邻域点为，如果
则判定在中垂线附近，并把归为搜索的邻域点之一。如果，则上述算法和传统的贪婪法相同。如果所有的邻域点均不满足公式（4），则待搜索的点为，并满足：
五、結果和讨论
使用上述算法对脑部MRI图像进行处理，得到的水肿区域的分割结果如图6所示。其中图6（a）是待分割图像的初始轮廓线。图6（b）是采用5X5邻域的Greedy法进行搜索的轮廓线的收敛结果，该结果使用了图4作为原始图像的外力，如前述，梯度力是一种短程力，难以捕捉较远范围的轮廓线，从而易使轮廓线收敛到局部极小值，并且难以收敛到目标区域的凹陷处，而且由于目标区域的边界较弱，导致轮廓线最终的收敛结果和目标区域仍然有较大的差别。图6（c）是使用本文算法的收敛的结果。使用静电场模型的图像外力是一种远路力，目标区域边缘的每一个象素点对图像上的任何地方均有力的作用，这种力反比与距离的平方。同时目标区域边缘的凹陷处也有力的作用，轮廓线在这种力的作用下容易收敛到真实的边界，从而得到较好的结果。
使用静电场模型计算当前位置图像外力的时候需要对目标区域各个象素点的贡献进行叠加。计算量较大，但是由于本文采用的初始轮廓线为目标区域的外轮廓线，所以不必计算轮廓线外部的电场力分布，同时本文采用改进的Greedy法对蛇点邻域进行搜索，故也不必计算蛇点邻域中所有点的静电力的分布，从而减小了计算量，由于蛇点邻域采用判据（4），（5）得到的搜索点的个数远少于邻域中所有点的个数，同时判据（4），（5）计算量远少于（3）的计算量，所以使用本文的算法可以近一步加快收敛的速度。
本研究使用二维直方图求得到图像二值化阈值，对二值图像的梯度图使用静电场模型计算图像外力，并提出蛇点邻域搜索的新准则，初始轮廓线在内力和外力的驱动下，最终收敛到的脑部图像的水肿瘤区域的边缘。将这种方法应用于实际的脑部MRI图像，得到了满意的结果。
参考文献：
[1] Amini A A， Weymouth TE， Jain R C. Using dynamic programming for solving variational problem in vision [J]， IEEE Transactions on Pattern Analysis and Machine Intelligence， 1990， 12 （9） ： 855～867
[2] Williams D J， Shab M. A fast algorithms， for active contours and curvature estimation [J]， CVGIP： Image Understanding， 1992， 55（1）： 14～26
[3] Kenji Suzuki， Isao Horiba， Boboru Sugie， Michio Nanki， Contour Extraction of Left Ventricular Cavity [J]， Digital systems and Computer in Japan Vol.34， No.2，2003
[4] C. Xu and J.L Price. Snakes， Shapes， and Gradient Vector Flow [J]， IEEE Trans. Pattern Analysis and Machine Intelligence， Vol.7， No.3， pp.359-369， Mar.1998
[5] Andrei C，Jalba， Michael H F.W， Jos B.T.M R. CPM： A Deformable Model for Shape Recovery and Segmentation Based on Charged Particles [J] ， IEEE Trans. Pattern Analysis and Machine Intelligence， VOL26， No.10， Oct， 2004
[6] Otsu N. A threshold selection method from gray-level histogram [J]. Automatica， 1979， 9； 62-69
[7] 陈果，左洪福. 图像分割的二维最大熵遗传算法[J]， 计算机辅助设计与图形学学报， 2002，14 （6） ： 530 - 534.
[8] 周露芳，古乐野. 基于量子遗传算法的二维最大熵图像分割[J] ， 计算应用， 2005， 25（8）
【关键词】主动轮廓模型 贪婪法 二維直方图 静电场 分割
一、引言
由于医学图像的数据量大，目标区域的解剖结构复杂多变，并且图像中伪影和各种噪声的影响使得目标的边缘往往模糊而不连续。这使得准确分割医学图像非常困难，由于问题的复杂性和重要性，如何构造有效的分割算法是一个很有挑战的工作。
被称为Snake的主动轮廓模型（Active Contour Model）由Kass等人提出，并被广泛的应用在一系列的图像视觉问题中。经过初始轮廓的估计后，Snake能够在轮廓线内力和图像外力的作用下自动的收敛能量的极小状态，最终运动到目标的边缘附近，从而完成图像的分割。近年来学者提出种种新方法，对轮廓线蛇点邻域的搜索和Snake的内力和外力的计算进行了各种改进，试图在各种约束条件下找到一个平衡点[1-2][4][5]，本文使用遗传算法首先对脑部MRI图像的二维直方图进行处理以提取的阈值，对图像进行二值处理后，求得二值图像的梯度边缘，并使用静电力模型对初始轮廓线内部的区域计算图像的外力，从而使轮廓线收敛，达到分割脑部MRI图像的目的。
二、使用主动轮廓模型进行图像分割
（一）主动轮廓模型的描述
参数化的Snake是一条动态的轮廓线，其参数方程为，整个轮廓线的能量定义为：
其中，，分别Snake的内力强度，外力强度，和人为交互力的强度。如果不需要人为干预Snake的行为的话，通常可以取。轮廓线的内力强度为：
式中的第一项为弹性势能，当轮廓线长度缩短时，该项积分取较小值，第二项为刚性能量，在轮廓线趋于平滑时，该项的积分取较小值，参数分别是弹性能量和刚性能量的控制系数，当两参数分别为零时，允许蛇点在此处的弹性能量和刚性能量取一个极大值，即允许轮廓线出现不连续，或曲率较大。所以，把这两个参数作为时变参数的话，可以近一步的影响蛇的行为。一般处理中为了简单，通常可以令这两个参数为常量。
（二）图像的外力
图像的外力一般可以取下面的形式：
即图像灰度和高斯函数卷积的梯度模值平方的负数，其中为控制系数。简化的外力还可以取公式（1）中高斯函数取极值的情况
但基于梯度的外力往往是一种短程力，难以捕获远处的轮廓线，并且轮廓线在这种外力的作用下难以收敛到目标区域的凹陷处。
三、基于静电场模型的Snake算法
（一）图像的二值化处理
本文使用Snake模型计算脑部MRI图像外力前，首先对图像进行二值化处理。进行二值化处理可以利用图像的直方图。图1（b）是图1左图脑部水肿图像的直方图，由于直方图的双峰间的峰谷往往过于平坦，使用一维的Otsu方法往往难以得到理想的灰度阈值。
为了得到较优的灰度阈值，可以使用二维直方图，本文使用基于原始图和其八邻域平滑图得到的二维直方图：设图像的大小为，在象素处的灰度为，象素点的八邻域均值为：
设原始图和平滑图为，联合象素满足的概率为，其中，则以为自变量，函数为联合图像的二维灰度直方图。图2为图1（a）的二维直方图的灰度图。
其中高亮度的区域对应二维直方图中取值大的部分，由图可以见，灰度二维直方图取值较大的部分集中在xy平面的对角线附近。 取x轴阈值s，和y轴阈值t，可以把二维直方图分为如下四个部分。
其中和区对应图像中的噪声和图边缘部分，概率较小，从而阈值划分问题转化为求阈值矢量（s，t），使得熵函数
对于公式（3）中最优参数，可以采用各种优化算法进行计算 [7][8]，例如可以使用遗传算法计算得到该类图像灰度阈值，由脑部MRI图像的先验知识可知，这类图像的阈值分布在100左右，可以尽一步缩小优化算法的搜索空间，加快算法的收敛速度。
本文使用的遗传算法采用直接对二维阈值矢量进行编码，设阈值分布以100为中心，幅度为±32的区间内，即[69，132]，对二维阈值矢量可以采用两基因（gene）的染色体编码，每个基因采用采用6比特编码，每个染色体（chromosome）为12比特，前面的6个比特为第一个基因，为二维直方图第一个坐标阈值的编码。后面第二个基因代表第二个坐标阈值。此时搜索空间中的点为61X61。染色体采用公式（2）作为适应度函数。
基本步骤：（a）初始染色体种群，种群数量为N=40。（b）对染色体解码，并计算适应度。（c）淘汰适应度最差的20个染色体，并对最优的20个染色体进行交叉，产生新的染色体，使种群恢复到原来的规模，并采用0.1概率对现有染色体进行变异。（d）如果近10代染色体种群的最优适应值保持不变，迭代结束，否则转（b）。
由于使用了图像的先验知识，大大减少了计算的工作量。否则要假设阈值分布在[0，255]，寻优空间的大小为255X255。实际上满足其他条件（例如图像的尺寸，和脑部图像与背景比例等）的脑部MRI图像的阈值将分布在更窄取值范围内，阈值矢量的搜索空间可以近一步的缩小。搜索空间足够小时，还可以使用全局搜索的方法，从而得到最佳的阈值。图3为使用遗传算法计算得到二维最大熵阈值，并使用该阈值对图1中的左图进行二值化的结果。
（二）基于静电场模型的图像外力的计算
由于主动轮廓模型传统外力的限制，学者提出了种种改进的外力模型[4][5]，本文使用文献[5]的思想，使用二值图像的梯度图在初始轮廓线内部区域的静电场力作为主动轮廓模型的图像外力。 设图像中每个象素的位置放置负电荷
，M是图像象素的个数，是象素处的灰度值，每个蛇点带单位正电荷，由于图像边缘的电场强度较大，蛇点将在电场的作用下向物体的边缘进行移动。同时，由于蛇点之间斥力的作用，蛇点将趋于均匀分布。设蛇点的数量为N，则轮廓线在静电场中的势能为：
其中分别为标号为的蛇点的矢量坐标，为标号为的象素的矢量坐标。带单位正电荷的蛇点在静电场中所受的引力为：
所受其他蛇点的斥力为：
为了突出二值图像边缘的影响，对图3进行求梯度的运算，并以梯度模值表征图像的边缘。图4（a）为是梯度运算的结果。图4（b）为使用静电场模型对图4（a）的目标区域（矩形框中）进行计算后得到的图像外力的结果。
由图4（b）可以看出，在水肿区域的边缘处形成较强的外力，即使目标区域的凹陷区也有外力的存在。
四、蛇点邻域的快速搜索算法
本文采用了改进的Greedy法在邻域中移动蛇点。为了使蛇点间的距离更加均匀，使目标蛇点在相邻蛇点的中垂线上移动是合理的。当前蛇点为，相邻蛇点为，两个相邻蛇点的中垂线为，如果和的邻域相交，则相交的邻域点为贪婪法的搜索点（图5中阴影部分），否则，把当前蛇点移动到离中垂线距离最近的邻域点。判断邻域点和是否相交，可以采用下面的准则：
设为定义好的一个比较小的值，设的邻域点为，如果
则判定在中垂线附近，并把归为搜索的邻域点之一。如果，则上述算法和传统的贪婪法相同。如果所有的邻域点均不满足公式（4），则待搜索的点为，并满足：
五、結果和讨论
使用上述算法对脑部MRI图像进行处理，得到的水肿区域的分割结果如图6所示。其中图6（a）是待分割图像的初始轮廓线。图6（b）是采用5X5邻域的Greedy法进行搜索的轮廓线的收敛结果，该结果使用了图4作为原始图像的外力，如前述，梯度力是一种短程力，难以捕捉较远范围的轮廓线，从而易使轮廓线收敛到局部极小值，并且难以收敛到目标区域的凹陷处，而且由于目标区域的边界较弱，导致轮廓线最终的收敛结果和目标区域仍然有较大的差别。图6（c）是使用本文算法的收敛的结果。使用静电场模型的图像外力是一种远路力，目标区域边缘的每一个象素点对图像上的任何地方均有力的作用，这种力反比与距离的平方。同时目标区域边缘的凹陷处也有力的作用，轮廓线在这种力的作用下容易收敛到真实的边界，从而得到较好的结果。
使用静电场模型计算当前位置图像外力的时候需要对目标区域各个象素点的贡献进行叠加。计算量较大，但是由于本文采用的初始轮廓线为目标区域的外轮廓线，所以不必计算轮廓线外部的电场力分布，同时本文采用改进的Greedy法对蛇点邻域进行搜索，故也不必计算蛇点邻域中所有点的静电力的分布，从而减小了计算量，由于蛇点邻域采用判据（4），（5）得到的搜索点的个数远少于邻域中所有点的个数，同时判据（4），（5）计算量远少于（3）的计算量，所以使用本文的算法可以近一步加快收敛的速度。
本研究使用二维直方图求得到图像二值化阈值，对二值图像的梯度图使用静电场模型计算图像外力，并提出蛇点邻域搜索的新准则，初始轮廓线在内力和外力的驱动下，最终收敛到的脑部图像的水肿瘤区域的边缘。将这种方法应用于实际的脑部MRI图像，得到了满意的结果。
参考文献：
[1] Amini A A， Weymouth TE， Jain R C. Using dynamic programming for solving variational problem in vision [J]， IEEE Transactions on Pattern Analysis and Machine Intelligence， 1990， 12 （9） ： 855～867
[2] Williams D J， Shab M. A fast algorithms， for active contours and curvature estimation [J]， CVGIP： Image Understanding， 1992， 55（1）： 14～26
[3] Kenji Suzuki， Isao Horiba， Boboru Sugie， Michio Nanki， Contour Extraction of Left Ventricular Cavity [J]， Digital systems and Computer in Japan Vol.34， No.2，2003
[4] C. Xu and J.L Price. Snakes， Shapes， and Gradient Vector Flow [J]， IEEE Trans. Pattern Analysis and Machine Intelligence， Vol.7， No.3， pp.359-369， Mar.1998
[5] Andrei C，Jalba， Michael H F.W， Jos B.T.M R. CPM： A Deformable Model for Shape Recovery and Segmentation Based on Charged Particles [J] ， IEEE Trans. Pattern Analysis and Machine Intelligence， VOL26， No.10， Oct， 2004
[6] Otsu N. A threshold selection method from gray-level histogram [J]. Automatica， 1979， 9； 62-69
[7] 陈果，左洪福. 图像分割的二维最大熵遗传算法[J]， 计算机辅助设计与图形学学报， 2002，14 （6） ： 530 - 534.
[8] 周露芳，古乐野. 基于量子遗传算法的二维最大熵图像分割[J] ， 计算应用， 2005， 25（8）