专利名称：大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法
技术领域：
本发明涉及图像处理方法，特别是大鼠头部磁共振图像的处理方法。
背景技术：
通过分析入射到生物组织内部的光和由生物组织散射出的光，可以获得生物组织 内部的光学参数信息，其中迫切需要解决的问题就是研究光在复杂生物组织中的传播和分 布规律。蒙特卡洛方法(以下简称MC)是广泛用于研究光在生物组织中传播和分布的随机 统计方法。对于非均勻复杂组织中光的分布MC模拟，首先要根据不同的光学参数(主要是 吸收系数(absorption coefficient)、散射系数(scattering coefficient)、各向异性因 子(anisotropic factor))将生物组织分割成不同的区域，从而构成多组织MC仿真模型， 然后进行MC计算，获得相应的光场分布。在早期的研究中，仿真模型是由几个简单的几何 模型组成，与真实组织相比非常不精确，从而导致仿真结果出现较大偏差。随着医学影像技 术的成熟，通过对MRI图像进行分割，从而获得多组织MC仿真模型被证明是非常有效的。在 对于脑组织的研究中，根据光学参数的不同，比较公认的是把脑组织分割成为5个部分，包 括头皮、头骨、脑脊液、脑灰质和脑白质。从MRI图像如何准确的分割出大脑的5种不同组 织，成为利用MC仿真方法获得准确结果的关键。针对人脑组织提出了许多图像分割方法， 其中最具有代表性的方法包括阈值分割法、水平集分割法、模糊C均值聚类方法等。阈值分割法分为全局阈值法和局部阈值分割法，全局阈值分割方法在图像处理中 应用比较多，它在整幅图像内采用固定的阈值分割图像。如果图像中待分割组织的灰度信 息重叠较多的话，效果很差。水平集分割法的核心思想是把η维描述视为高一维(η+1)维的水平集，或者说是 把η维描述视为有η维变量的水平集函数f的水平集。这样一来就把求解η维描述的演化 过程转化为求解关于有η维变量的水平集函数f的演化所导致的水平集的演化过程。K均值聚类本身是迭代的，且不能确保它收敛于最优解。K均值算法的性能依赖于 聚类中心的初始位置。模糊C均值聚类方法，作为K均值聚类算法的一种改进，是用隶属度 确定每个数据点属于某个聚类的程度的一种聚类算法。模糊C均值聚类是目前分割MRI图 像最好的方法之一，它能很好的分割出颅脑图像中的灰质、白质和脑脊液等组织。由于大鼠模型广泛应用于科学研究中，因此建立满足MC仿真要求的大鼠大脑多 组织模型具有特别重要的意义。早期研究中，大鼠大脑多组织模型的建立方法，主要采用模 板匹配分割法。随着高场强小动物磁共振的出现，通过分割MRI获得大鼠大脑多组织MC仿 真模型，成为迫切的需要。大鼠大脑磁共振图像序列中，头皮组织的灰度分布范围广，与脑 组织(包括脂肪组织、脑脊液、脑白质、脑灰质)有很大重叠。因而直接通过灰度信息，利用 阈值分割或模糊聚类并不能够完成5中组织的分割；由于大鼠大脑体积小，颅骨很薄，利用 窄带水平集同样不能获得理想的结果。因此单独采用现有的分割方法都无法获得满足MC 仿真要求的多组织模型。
发明内容
本发明的目的在于克服现有技术的局限性，提供一种快速、准确的大鼠头部磁共 振图像蒙特卡洛仿真模型的获取方法。本发明的思路是针对大鼠头部磁共振图像中各组织的特点，选用最合适的图像分 割方法和相应的参数，从整体上通过多种图像分割方法的特定组合，将头皮、颅骨、脑脊液、 脑灰质、脑白质这5种组织图像准确分割出来，从而实现快速、准确的获取蒙特卡洛仿真模型。具体的说，本发明是通过如下各步骤实现的A、首先将图像分割为头皮、颅骨、脑脊液、脑灰质、脑白质这5种组织以及图像背 景共6部分图像；B、然后将上述6部分图像分别赋以不同的灰度值；C、最后将6部分图像合并，得到蒙特卡洛仿真模型。其中步骤A采用了多种图像分割方法的特定组合，具体包括下列各步骤Al、颅骨组织分割使用灰度阈值分割法对原始图像序列中的任意一张图像KO进 行分割，并利用膨胀腐蚀算法获得颅骨组织和不连续的噪声部分；然后利用二值图像最大 区域标记法去除噪声，即首先对图像不连续区域灰度值进行标记，然后计算各个标记区域 的面积，其中面积最大的连续区域即为需分割的组织，将其像素标记为1，其他部分为噪声 部分，其像素标记为0，得到颅骨图像Kl ；本步骤中的膨胀腐蚀算法是现有技术，具体内容参见文献(医学图像处理与分 析，罗述谦等编著，科学出版社，2003，第83-89页)A2、头皮组织分割首先对图像Kl利用边缘检测方法得到颅骨上边界的三维坐标 集；根据得到的颅骨上边界三维坐标对原始图像KO进行二值化处理，其中颅骨外边界上方 的区域像素赋为0，得到去除了部分头皮的图像K2 ；对图像K2利用窄带水平集方法提取包 含脑脊液、脑灰质、脑白质的脑组织部分的边界，得到图像K3;对图像K3进行二值化处理， 其中脑组织边界上以及边界内区域像素值赋为1，其他部分像素值赋为0，将处理后的结果 与原始图像KO相乘，得到图像K4 ；利用算式K5 = K0-K4-K1得到新的图像，并利用膨胀腐 蚀算法和二值图像最大区域标记法去除不连续的噪声部分，得到头皮图像K5 ；本步骤中所述边缘检测方法有多种现有技术可供选择，例如微分算子法等(医 学图像处理与分析，罗述谦等编著，科学出版社，2003，第68-75页)等；本步骤中使用的二值图像最大区域标记法与步骤Al中相同。A3、将原始图像KO去除颅骨部分和头皮部分，利用模糊聚类算法对图像进行分 害I]，对分割出的脑脊液、脑灰质、脑白质这三部分分别赋以不同的灰度值，得到图像K6;A4、对原始图像序列中的所有图像按照A1-A3步骤进行处理。本发明由于针对各组织的特点，分别选择了最适合的图像分割方法，且通过一定 的二值化处理方式，大大简化了图像处理的复杂度，因此，相比现有技术，具有准确、快速的 优点，具有良好的实际应用前景。
图1为本发明方法的流程图2为从原始大鼠头部磁共振图像序列中选取的第三层切片；图3为本发明方法得到的蒙特卡洛仿真模型。
具体实施例方式下面结合附图对本发明的技术方案进行详细说明如附图1所示，本发明按照如下步骤执行A、首先将图像分割为头皮、颅骨、脑脊液、脑灰质、脑白质这5种组织以及图像背景共6部分图像；Al、颅骨组织分割使用灰度阈值分割法对原始图像序列中的任意一张图像 KO (参见附图2)进行分割，并利用膨胀腐蚀算法获得颅骨组织和不连续的噪声部分；然后 利用二值图像最大区域标记法去除噪声，即首先对图像不连续区域灰度值进行标记，然后 计算各个标记区域的面积，其中面积最大的连续区域即为需分割的组织，将其像素标记为 1，其他部分为噪声部分，其像素标记为0，得到颅骨图像Kl ；A101、任意选取原始图像序列中的一张图像，手工标记出颅骨部分；A102、根据以下方程求取当0. 8彡S彡0. 85时的灰度值I 其中Pi是颅骨部分图像中灰度值为i的像素点个数，S为自定义的分割系数，Im 是颅骨部分图像中最大的灰度值；A103、将步骤A102中求得的I值作为灰度阈值，利用灰度阈值分割法对所述图像 KO进行分割；A104、对步骤A103得到的图像中所有像素点按照以下公式赋以不同的灰度值 M(i) 其中i是步骤A103得到的图像中任意点的灰度值；A105、选用大小为LXL的矩形模板对步骤A104得到的图像进行膨胀运算，其中L 取值为2-5 ；A106、对步骤A105得到的图像进行腐蚀运算，腐蚀模板为(L_l) X (L-I)的矩形；A107、将步骤A106得到的图像中各连续区域分别赋以不同灰度值，同一连续区域 中的所有像素点灰度值相同；包含像素点最多的灰度值所对应的连续区域即为颅骨部分；A108、将颅骨部分中所有像素点灰度值设为1，其他部分像素点灰度值设为0，得 到颅骨图像Kl。A2、头皮组织分割首先对图像Kl利用边缘检测方法得到颅骨上边界的三维坐标 集；根据得到的颅骨上边界三维坐标对原始图像KO进行二值化处理，其中颅骨外边界上方 的区域像素赋为0，得到去除了部分头皮的图像K2 ；对图像K2利用窄带水平集方法提取包含脑脊液、脑灰质、脑白质的脑组织部分的边界，得到图像K3;对图像K3进行二值化处理， 其中脑组织边界上以及边界内区域像素值赋为1，其他部分像素值赋为0，将处理后的结果 与原始图像KO相乘，得到图像Κ4 ；利用算式Κ5 = Κ0-Κ4-Κ1得到新的图像，并利用膨胀腐 蚀算法和二值图像最大区域标记法去除不连续的噪声部分，得到头皮图像Κ5 ； Α201、对步骤A得到的图像Κ1，利用边缘检测算法求得颅骨的上边界三维坐标集 (Χ，Υ，Ζ)，其中图像序列中第1幅图像的左上角为原点，X轴正方向指向右侧，Y轴正方向指 向下侧，坐标系为右手坐标系；Α202、根据下式将原始图像KO中坐标为(x，y，z)的像素点赋以相应的灰度值M(x， 7，2)，得到图像1(2: A203、对图像K2利用下式定义的水平集算法提取脑组织，得到包含脑脊液、脑灰 质、脑白质的脑组织边界的图像K3 上式 由文献(Level Set Evolution Without Re-initialization :A NewVariational Formulation)获得，为水平集进化函数。上式可通过文献(Fast DistancePreserving Level Set Evolution for Medical Image Segmentation)中窄带
方法快速计算。其中，Φ为水平集，t为时间，
为边界指示函数，G0是标准
差为σ的高斯核，λ为水平集曲线长度的权重，μ为内部能量权重，ν为收敛驱动力权重， μ，λ，ν的取值需要根据实际需要确定，本实施方案中，μ = 0.04, λ = 1. 21，ν = 12, τ = 5，迭代次数至少为50次；Α204、根据下式将图像Κ3中坐标为(x，y，ζ)的像素点赋以相应的灰度值M(x，y, ζ)，并将得到的结果与原始图像KO相乘，得到图像Κ4 ； Α205、根据算式Κ5 = Κ0-Κ4-Κ1得到新的图像；Α206、选用大小为LXL的矩形模板对步骤Α205得到的图像进行膨胀运算，其中L 取值为2-5 ；Α207、对步骤Α206得到的图像进行腐蚀运算，腐蚀模板为(L_l) X (L-I)的矩形；A208、将步骤A207得到的图像中各连续区域分别赋以不同灰度值，同一连续区域 中的所有像素点灰度值相同；提取包含像素点最多的灰度值所对应的连续区域，即得到头 皮图像K5。A3、将原始图像KO去除颅骨部分和头皮部分，利用模糊聚类算法对图像进行分 害I]，对分割出的脑脊液、脑灰质、脑白质这三部分分别赋以不同的灰度值，得到图像K6;A4、对原始图像序列中的所有图像按照A1-A3步骤进行处理。
B、将上述6部分图像分别赋以不同的灰度值；本实施方案中将图像ΚΙ、K5、K6中的头皮、颅骨、脑脊液、脑灰质和脑白质灰度值分别赋为1、2、3、4、5，背景为0 ；当然，也可取其他的灰度值，只要各部分组织及背景的灰度 不同即可。C、最后将6部分图像合并，得到蒙特卡洛仿真模型，如附图3所示。本发明可结合计算机系统并利用现有的图像处理软件完成。
权利要求
大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法，包括下述步骤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所述大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法，其特征在 于所述步骤Al包含以下步骤A101、任意选取原始图像序列中的一张图像，手工标记出颅骨部分；A102、根据以下方程求取当0. 8 < S < 0. 85时的灰度值I : 其中Pi是颅骨部分图像中灰度值为i的像素点个数，S为自定义的分割系数，Im是颅 骨部分图像中最大的灰度值；A103、将步骤A102中求得的I值作为灰度阈值，利用灰度阈值分割法对所述图像KO进 行分割；A104、对步骤A103得到的图像中所有像素点按照以下公式赋以不同的灰度值M(i) 其中i是步骤A103得到的图像中任意点的灰度值；A105、选用大小为LXL的矩形模板对步骤A104得到的图像进行膨胀运算，其中L取值 为 2-5 ；A106、对步骤A105得到的图像进行腐蚀运算，腐蚀模板为(L-I)X(L-I)的矩形；A107、将步骤A106得到的图像中各连续区域分别赋以不同灰度值，同一连续区域中的所有像素点灰度值相同；包含像素点最多的灰度值所对应的连续区域即为颅骨部分；A108、将颅骨部分中所有像素点灰度值设为1，其他部分像素点灰度值设为0，得到颅 骨图像K1。
3.如权利要求1所述大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法，其特征在 于所述步骤A2包括以下各步骤A201、对步骤A得到的颅骨图像K1，利用边缘检测算法求得颅骨的上边界三维坐标集 (X，Y，Z)，其中图像序列中第1幅图像的左上角为原点，X轴正方向指向右侧，Y轴正方向指 向下侧，坐标系为右手坐标系；A202、根据下式将原始图像KO中坐标为(X，y，ζ)的像素点赋以相应的灰度值M(x，y, z)，得到图像K2: 其他A203、对图像K2利用下式定义的水平集算法提取脑组织，得到包含脑脊液、脑灰质、脑 白质的脑组织边界的图像K3 其中，Φ为水平集，t为时间， 为边界指示函数，G。是标准差为σ的高斯核；λ为水平集曲线长度的权重，μ为内部能量权重，ν为收敛驱动力权重；Α204、根据下式将图像Κ3中坐标为(x，y，z)的像素点赋以相应的灰度值M(x，y，ζ)，并 将得到的结果与原始图像KO相乘，得到图像Κ4 ； Α205、根据算式Κ5 = Κ0-Κ4-Κ1得到新的图像；Α206、选用大小为LXL的矩形模板对步骤Α205得到的图像进行膨胀运算，其中L取值 为 2-5 ；Α207、对步骤Α206得到的图像进行腐蚀运算，腐蚀模板为(L-I)X(L-I)的矩形； Α208、将步骤Α207得到的图像中各连续区域分别赋以不同灰度值，同一连续区域中的 所有像素点灰度值相同；提取包含像素点最多的灰度值所对应的连续区域，即得到头皮图 像Κ5。
全文摘要
本发明公开了一种大鼠头部磁共振图像蒙特卡洛仿真模型的获取方法，属于组织光学技术领域。本发明基于多种图像分割方法的特定组合，首先通过分析确定灰度阈值，采用阈值分割、腐蚀和膨胀算法、二值图像最大区域标记法分割出颅骨组织；然后去除原始图像中颅骨上部所有的组织，采用特殊参数的窄带水平集分割出头皮组织；之后在去除头皮和颅骨的图像中，采用快速模糊聚类分割出脑脊液、脑灰质和脑白质；最后将分割出的5种组织以及背景部分，分别赋予不同的灰度值，获得多组织模型。本发明可以克服单一图像分割方法的局限性，快速、准确的获取大鼠头部磁共振图像蒙特卡洛仿真模型，具有良好的实际应用前景。
文档编号G06T7/00GK101847261SQ20101016748
公开日2010年9月29日 申请日期2010年5月10日 优先权日2010年5月10日
发明者李韪韬, 王新增, 王雪娜, 钱志余 申请人:南京航空航天大学