CN101847261B - 大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法 - Google Patents

Abstract

本发明公开了一种大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法，属于组织光学技术领域。本发明基于多种图像分割方法的特定组合，首先通过分析确定灰度阈值，采用阈值分割、腐蚀和膨胀算法、二值图像最大区域标记法分割出颅骨组织；然后去除原始图像中颅骨上部所有的组织，采用特殊参数的窄带水平集分割出头皮组织；之后在去除头皮和颅骨的图像中，采用快速模糊聚类分割出脑脊液、脑灰质和脑白质；最后将分割出的5种组织以及背景部分，分别赋予不同的灰度值，获得多组织模型。本发明可以克服单一图像分割方法的局限性，快速、准确的获取大鼠头部磁共振图像蒙特卡洛仿真模型，具有良好的实际应用前景。

Description

大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法

技术领域

本发明涉及图像处理方法，特别是大鼠头部磁共振图像的处理方法。

背景技术

通过分析入射到生物组织内部的光和由生物组织散射出的光，可以获得生物组织内部的光学参数信息，其中迫切需要解决的问题就是研究光在复杂生物组织中的传播和分布规律。蒙特卡洛方法(以下简称MC)是广泛用于研究光在生物组织中传播和分布的随机统计方法。对于非均匀复杂组织中光的分布MC模拟，首先要根据不同的光学参数(主要是：吸收系数(absorption coefficient)、散射系数(scattering coefficient)、各向异性因子(anisotropic factor))将生物组织分割成不同的区域，从而构成多组织MC仿真模型，然后进行MC计算，获得相应的光场分布。在早期的研究中，仿真模型是由几个简单的几何模型组成，与真实组织相比非常不精确，从而导致仿真结果出现较大偏差。随着医学影像技术的成熟，通过对MRI图像进行分割，从而获得多组织MC仿真模型被证明是非常有效的。在对于脑组织的研究中，根据光学参数的不同，比较公认的是把脑组织分割成为5个部分，包括头皮、头骨、脑脊液、脑灰质和脑白质。从MRI图像如何准确的分割出大脑的5种不同组织，成为利用MC仿真方法获得准确结果的关键。针对人脑组织提出了许多图像分割方法，其中最具有代表性的方法包括阈值分割法、水平集分割法、模糊C均值聚类方法等。

阈值分割法分为全局阈值法和局部阈值分割法，全局阈值分割方法在图像处理中应用比较多，它在整幅图像内采用固定的阈值分割图像。如果图像中待分割组织的灰度信息重叠较多的话，效果很差。

水平集分割法的核心思想是把n维描述视为高一维(n+1)维的水平集，或者说是把n维描述视为有n维变量的水平集函数f的水平集。这样一来就把求解n维描述的演化过程转化为求解关于有n维变量的水平集函数f的演化所导致的水平集的演化过程。

K均值聚类本身是迭代的，且不能确保它收敛于最优解。K均值算法的性能依赖于聚类中心的初始位置。模糊C均值聚类方法，作为K均值聚类算法的一种改进，是用隶属度确定每个数据点属于某个聚类的程度的一种聚类算法。模糊C均值聚类是目前分割MRI图像最好的方法之一，它能很好的分割出颅脑图像中的灰质、白质和脑脊液等组织。

由于大鼠模型广泛应用于科学研究中，因此建立满足MC仿真要求的大鼠大脑多组织模型具有特别重要的意义。早期研究中，大鼠大脑多组织模型的建立方法，主要采用模板匹配分割法。随着高场强小动物磁共振的出现，通过分割MRI获得大鼠大脑多组织MC仿真模型，成为迫切的需要。大鼠大脑磁共振图像序列中，头皮组织的灰度分布范围广，与脑组织(包括脂肪组织、脑脊液、脑白质、脑灰质)有很大重叠。因而直接通过灰度信息，利用阈值分割或模糊聚类并不能够完成5中组织的分割；由于大鼠大脑体积小，颅骨很薄，利用窄带水平集同样不能获得理想的结果。因此单独采用现有的分割方法都无法获得满足MC仿真要求的多组织模型。

发明内容

本发明的目的在于克服现有技术的局限性，提供一种快速、准确的大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法。

本发明的思路是针对大鼠头部磁共振图像中各组织的特点，选用最合适的图像分割方法和相应的参数，从整体上通过多种图像分割方法的特定组合，将头皮、颅骨、脑脊液、脑灰质、脑白质这5种组织图像准确分割出来，从而实现快速、准确的获取蒙特卡洛仿真模型。

具体的说，本发明是通过如下各步骤实现的：

A、首先将图像分割为头皮、颅骨、脑脊液、脑灰质、脑白质这5种组织以及图像背景共6部分图像；

B、然后将上述6部分图像分别赋以不同的灰度值；

C、最后将6部分图像合并，得到蒙特卡洛仿真模型。

其中步骤A采用了多种图像分割方法的特定组合，具体包括下列各步骤：

A1、颅骨组织分割：使用灰度阈值分割法对原始图像序列中的任意一张图像K0进行分割，并利用膨胀腐蚀算法获得颅骨组织和不连续的噪声部分；然后利用二值图像最大区域标记法去除噪声，即首先对图像不连续区域灰度值进行标记，然后计算各个标记区域的面积，其中面积最大的连续区域即为需分割的组织，将其像素标记为1，其他部分为噪声部分，其像素标记为0，得到颅骨图像K1；

本步骤中的膨胀腐蚀算法是现有技术，具体内容参见文献(医学图像处理与分析，罗述谦等编著，科学出版社，2003，第83-89页)

A2、头皮组织分割：首先对图像K1利用边缘检测方法得到颅骨上边界的三维坐标集；根据得到的颅骨上边界三维坐标对原始图像K0进行二值化处理，其中颅骨外边界上方的区域像素赋为0，得到去除了部分头皮的图像K2；对图像K2利用窄带水平集方法提取包含脑脊液、脑灰质、脑白质的脑组织部分的边界，得到图像K3；对图像K3进行二值化处理，其中脑组织边界上以及边界内区域像素值赋为1，其他部分像素值赋为0，将处理后的结果与原始图像K0相乘，得到图像K4；利用算式K5＝K0-K4-K1得到新的图像，并利用膨胀腐蚀算法和二值图像最大区域标记法去除不连续的噪声部分，得到头皮图像K5；

本步骤中所述边缘检测方法有多种现有技术可供选择，例如：微分算子法等(医学图像处理与分析，罗述谦等编著，科学出版社，2003，第68-75页)等；

本步骤中使用的二值图像最大区域标记法与步骤A1中相同。

A3、将原始图像K0去除颅骨部分和头皮部分，利用模糊聚类算法对图像进行分割，对分割出的脑脊液、脑灰质、脑白质这三部分分别赋以不同的灰度值，得到图像K6；

A4、对原始图像序列中的所有图像按照A1-A3步骤进行处理。

本发明由于针对各组织的特点，分别选择了最适合的图像分割方法，且通过一定的二值化处理方式，大大简化了图像处理的复杂度，因此，相比现有技术，具有准确、快速的优点，具有良好的实际应用前景。

附图说明

图1为本发明方法的流程图；

图2为从原始大鼠头部磁共振图像序列中选取的第三层切片；

图3为本发明方法得到的蒙特卡洛仿真模型。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

如附图1所示，本发明按照如下步骤执行：

A、首先将图像分割为头皮、颅骨、脑脊液、脑灰质、脑白质这5种组织以及图像背景共6部分图像；

A1、颅骨组织分割：使用灰度阈值分割法对原始图像序列中的任意一张图像K0(参见附图2)进行分割，并利用膨胀腐蚀算法获得颅骨组织和不连续的噪声部分；然后利用二值图像最大区域标记法去除噪声，即首先对图像不连续区域灰度值进行标记，然后计算各个标记区域的面积，其中面积最大的连续区域即为需分割的组织，将其像素标记为1，其他部分为噪声部分，其像素标记为0，得到颅骨图像K1；

A101、任意选取原始图像序列中的一张图像，手工标记出颅骨部分；

A102、根据以下方程求取当0.8≤S≤0.85时的灰度值I：

$S=\frac{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{p}_i}{\underset{i=1}{\overset{\mathrm{Im}}{\mathrm{\Σ}}}{p}_i}$

其中p_i是颅骨部分图像中灰度值为i的像素点个数，S为自定义的分割系数，Im是颅骨部分图像中最大的灰度值；

A103、将步骤A102中求得的I值作为灰度阈值，利用灰度阈值分割法对所述图像K0进行分割；

A104、对步骤A103得到的图像中所有像素点按照以下公式赋以不同的灰度值M(i)：

其中i是步骤A103得到的图像中任意点的灰度值；

A105、选用大小为L×L的矩形模板对步骤A104得到的图像进行膨胀运算，其中L取值为2-5；

A106、对步骤A105得到的图像进行腐蚀运算，腐蚀模板为(L-1)×(L-1)的矩形；

A107、将步骤A106得到的图像中各连续区域分别赋以不同灰度值，同一连续区域中的所有像素点灰度值相同；包含像素点最多的灰度值所对应的连续区域即为颅骨部分；

A108、将颅骨部分中所有像素点灰度值设为1，其他部分像素点灰度值设为0，得到颅骨图像K1。

A2、头皮组织分割：首先对图像K1利用边缘检测方法得到颅骨上边界的三维坐标集；根据得到的颅骨上边界三维坐标对原始图像K0进行二值化处理，其中颅骨外边界上方的区域像素赋为0，得到去除了部分头皮的图像K2；对图像K2利用窄带水平集方法提取包含脑脊液、脑灰质、脑白质的脑组织部分的边界，得到图像K3；对图像K3进行二值化处理，其中脑组织边界上以及边界内区域像素值赋为1，其他部分像素值赋为0，将处理后的结果与原始图像K0相乘，得到图像K4；利用算式K5＝K0-K4-K1得到新的图像，并利用膨胀腐蚀算法和二值图像最大区域标记法去除不连续的噪声部分，得到头皮图像K5；

A201、对步骤A得到的图像K1，利用边缘检测算法求得颅骨的上边界三维坐标集(X，Y，Z)，其中图像序列中第1幅图像的左上角为原点，X轴正方向指向右侧，Y轴正方向指向下侧，坐标系为右手坐标系；

A202、根据下式将原始图像K0中坐标为(x，y，z)的像素点赋以相应的灰度值M(x，y，z)，得到图像K2：

A203、对图像K2利用下式定义的水平集算法提取脑组织，得到包含脑脊液、脑灰质、脑白质的脑组织边界的图像K3：

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\μ}[\mathrm{\Δ\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)]+\mathrm{\λ\δ}\left(\mathrm{\φ}\right)\mathrm{div}\left(g\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+\mathrm{vg\δ}\left(\mathrm{\φ}\right)$

上式由文献(Level Set Evolution Without Re-initialization：A NewVariational Formulation)获得，为水平集进化函数。上式可通过文献(Fast DistancePreserving Level Set Evolution for Medical Image Segmentation)中窄带方法快速计算。其中，φ为水平集，t为时间，

为边界指示函数，G_σ是标准差为σ的高斯核，λ为水平集曲线长度的权重，μ为内部能量权重，v为收敛驱动力权重，μ，λ，v的取值需要根据实际需要确定，

本实施方案中，μ＝0.04，λ＝1.21，v＝12，τ＝5，迭代次数至少为50次；

A204、根据下式将图像K3中坐标为(x，y，z)的像素点赋以相应的灰度值M(x，y，z)，并将得到的结果与原始图像K0相乘，得到图像K4；

A205、根据算式K5＝K0-K4-K1得到新的图像；

A206、选用大小为L×L的矩形模板对步骤A205得到的图像进行膨胀运算，其中L取值为2-5；

A207、对步骤A206得到的图像进行腐蚀运算，腐蚀模板为(L-1)×(L-1)的矩形；

A208、将步骤A207得到的图像中各连续区域分别赋以不同灰度值，同一连续区域中的所有像素点灰度值相同；提取包含像素点最多的灰度值所对应的连续区域，即得到头皮图像K5。

A3、将原始图像K0去除颅骨部分和头皮部分，利用模糊聚类算法对图像进行分割，对分割出的脑脊液、脑灰质、脑白质这三部分分别赋以不同的灰度值，得到图像K6；

A4、对原始图像序列中的所有图像按照A1-A3步骤进行处理。

B、将上述6部分图像分别赋以不同的灰度值；

本实施方案中将图像K1、K5、K6中的头皮、颅骨、脑脊液、脑灰质和脑白质灰度值分别赋为1、2、3、4、5，背景为0；当然，也可取其他的灰度值，只要各部分组织及背景的灰度不同即可。

C、最后将6部分图像合并，得到蒙特卡洛仿真模型，如附图3所示。

本发明可结合计算机系统并利用现有的图像处理软件完成。

Claims (3)

1.大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法，包括下述步骤：A、首先将图像分割为头皮、颅骨、脑脊液、脑灰质、脑白质这5种组织以及图像背景共6部分图像；B、然后将上述6部分图像分别赋以不同的灰度值；C、最后将6部分图像合并，得到蒙特卡洛仿真模型；其特征在于：所述步骤A具体包括以下步骤：

A1、颅骨组织分割：使用灰度阈值分割法对原始图像序列中的任意一张图像K0进行分割，并利用膨胀腐蚀算法获得颅骨组织和不连续的噪声部分；然后利用二值图像最大区域标记法去除噪声，即首先对图像不连续区域灰度值进行标记，然后计算各个标记区域的面积，其中面积最大的连续区域即为需分割的组织，将其像素标记为1，其他部分为噪声部分，其像素标记为0，得到颅骨图像K1；

A2、头皮组织分割：首先对颅骨图像K1利用边缘检测方法得到颅骨上边界的三维坐标集；根据得到的颅骨上边界三维坐标对原始图像K0进行二值化处理，其中颅骨外边界上方的区域像素赋为0，得到去除了部分头皮的图像K2；对图像K2利用窄带水平集方法提取包含脑脊液、脑灰质、脑白质的脑组织部分的边界，得到图像K3；对图像K3进行二值化处理，其中脑组织边界上以及边界内区域像素值赋为1，其他部分像素值赋为0，将处理后的结果与原始图像K0相乘，得到图像K4；利用算式K5＝K0-K4-K1得到新的图像，并利用膨胀腐蚀算法和二值图像最大区域标记法去除不连续的噪声部分，得到头皮图像K5；

A3、将原始图像K0去除颅骨部分和头皮部分，利用模糊聚类算法对图像进行分割，对分割出的脑脊液、脑灰质、脑白质这三部分分别赋以不同的灰度值，得到图像K6；

A4、对原始图像序列中的所有图像按照A1-A3步骤进行处理。

2.如权利要求1所述大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法，其特征在于：所述步骤A1包含以下步骤：

A101、任意选取原始图像序列中的一张图像，手工标记出颅骨部分；

A102、根据以下方程求取当0.8≤S≤0.85时的灰度值I：

$S=\frac{\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{p}_i}{\underset{i=1}{\overset{\mathrm{Im}}{\mathrm{\Σ}}}{p}_i}$

其中p_i是颅骨部分图像中灰度值为i的像素点个数，S为自定义的分割系数，I_m是颅骨部分图像中最大的灰度值；

A103、将步骤A102中求得的I值作为灰度阈值，利用灰度阈值分割法对所述图像K0进行分割；

A104、对步骤A103得到的图像中所有像素点按照以下公式赋以不同的灰度值M(i)：

其中i是步骤A103得到的图像中任意点的灰度值；

A105、选用大小为L×L的矩形模板对步骤A104得到的图像进行膨胀运算，其中L取值为2-5；

A106、对步骤A105得到的图像进行腐蚀运算，腐蚀模板为(L-1)×(L-1)的矩形；

A107、将步骤A106得到的图像中各连续区域分别赋以不同灰度值，同一连续区域中的所有像素点灰度值相同；包含像素点最多的灰度值所对应的连续区域即为颅骨部分；

A108、将颅骨部分中所有像素点灰度值设为1，其他部分像素点灰度值设为0，得到颅骨图像K1。

3.如权利要求1所述大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法，其特征在于：所述步骤A2包括以下各步骤：

A201、对步骤A得到的颅骨图像K1，利用边缘检测算法求得颅骨的上边界三维坐标集(X，Y，Z)，其中图像序列中第1幅图像的左上角为原点，X轴正方向指向右侧，Y轴正方向指向下侧，坐标系为右手坐标系；

A202、根据下式将原始图像K0中坐标为(x，y，z)的像素点赋以相应的灰度值M(x，y，z)，得到图像K2：

A203、对图像K2利用下式定义的水平集算法提取脑组织，得到包含脑脊液、脑灰质、脑白质的脑组织边界的图像K3：

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\μ}[\mathrm{\Δ\φ}-\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)]+\mathrm{\λ\δ}\left(\mathrm{\φ}\right)\mathrm{div}\left(g\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+\mathrm{vg\δ}\left(\mathrm{\φ}\right)$

其中，φ为水平集，t为时间，

为边界指示函数，G_σ是标准差为σ的高斯核；λ为水平集曲线长度的权重，μ为内部能量权重，v为收敛驱动力权重；

A204、根据下式将图像K3中坐标为(x，y，z)的像素点赋以相应的灰度值M(x，y，z)，并将得到的结果与原始图像K0相乘，得到图像K4；

A205、根据算式K5＝K0-K4-K1得到新的图像；

A206、选用大小为L×L的矩形模板对步骤A205得到的图像进行膨胀运算，其中L取值为2-5；

A207、对步骤A206得到的图像进行腐蚀运算，腐蚀模板为(L-1)×(L-1)的矩形；

A208、将步骤A207得到的图像中各连续区域分别赋以不同灰度值，同一连续区域中的所有像素点灰度值相同；提取包含像素点最多的灰度值所对应的连续区域，即得到头皮图像K5。

Images