徐 林,陈明生,孙 建,金 贵,张海生,宁 旭,许 佳,王 凤,庄 伟,白泽霖,秦明新
(陆军军医大学生物医学工程与影像医学系电子工程学教研室,重庆 400038)
脑脊液是一种无色透明的液体,在蛛网膜下腔、脑室和脊髓中央管间循环流动,其具有保护大脑免受外部轻微撞击带来的损伤、为中枢神经系统提供营养以及稳定颅内化学调控等作用。大多数情况下,许多神经紊乱(如脑水肿和阿尔茨海默症)会影响脑脊液的体积和压力,因此不论在临床上还是科学研究中,监测、分析和研究脑脊液的变化都十分重要[1]。
而目前还缺乏在体量化分析脑脊液的方法。文献[2] 利用三维速度编码梯度回波相位对比成像来评估人脑内脑脊液的流动情况;文献[3] 提出了一种常规的MRI计算机数值化分析方法评估和诊断自发性脊髓脑脊液漏出症(spontaneous spinal cerebrospinal fluid leaks,SSCFL);因为腰骶脑脊液体积影响脊髓麻醉的临床反应,一些学者研究了基于MRI的腰骶脑脊液体积估算,文献[4] 利用了一种基于阈值的区域增长算法评估腰骶脑脊液体积;文献[5] 研究了量化分析腰骶脑脊液和颅内脑脊液时的阈值选择标准。
上述脑脊液评估和量化的方法大多数是针对人脑的。在研究动物模型脑脊液对脑出血检测[6-7]的影响时,出血和脑脊液的变化希望能够同步进行观测,这使得量化分析动物模型颅内的脑脊液变化非常重要。但是观测和量化动物模型颅内的脑脊液不同于人脑。首先,动物大脑的结构不同于人脑;其次,动物的脑组织在MRI图像上亮度和对比度都不同于人脑,并且动物脑组织图像的尺寸和分辨力不同于人脑;最后,大多数MRI处理算法都是针对人脑开发的且没有在动物上进行测试。当前也出现了一些针对动物模型的基于计算机辅助的图像分析技术。为分析家兔大脑的神经发育,文献[8] 建立了基于图谱的家兔脑组织MRI自动分割方法,但在该分割方法中脑脊液并未被单独分割出来;文献[9] 构建了一种针对小型动物的肺感染疾病的自动肺部图像分析框架;文献[10] 基于在相似性传播聚类框架内利用灰度相似性度量,提出了一种最优阈值选择方法,并应用于基于PET图像分割的小型动物肺结核的量化分析中。
本课题组在研究基于磁感应相移(magnetic induction phase shift,MIPS)[6,11-14]非接触检测家兔脑出血模型时,希望获得MIPS信号与家兔颅内主要组织(特别是脑脊液,因为脑脊液在颅内主要组织中的电导率最高)间的关系。针对这一目的,需要对家兔脑出血模型中的脑脊液体积进行量化,并描述脑脊液体积随着脑出血进程发展而变化的过程。然而当前还缺乏针对家兔的准确脑脊液量化方法,且上述关于动物模型的计算机辅助图像分析技术[8-10]也难以满足这一要求。
本文通过结合马尔科夫随机场(Markov random field,MRF)和模糊聚类(fuzzy clustering method,FCM)建立一种新的图像分割算法[15-16]。该方法中,FCM被用来提取和利用图像中的灰度信息,MRF被用来提取和利用图像中的空间信息和抑制噪声给图像处理带来的干扰。
首先,通过自动选取待分类类别数目以取得更准确的脑脊液量化结果;然后,通过该分割算法获得潜在的脑脊液区域,并与已经准备好的对应家兔脑脊液模板做掩码获得最终准确的家兔脑脊液区域;最后,通过家兔脑脊液分割区域结合MRI成像参数计算得到该家兔脑脊液体积的量化结果,量化分析流程如图1所示。
图1 MRI图像处理和脑脊液量化计算步骤
MRI图像中包含了2种重要信息:亮度信息和纹理信息。其中像素自身灰度值代表了该点亮度信息,反映了该点物理属性;像素与邻域内相邻像素构成了纹理信息。分割就是充分利用上述信息实现图像中相似性质区域的提取。对于亮度信息,不同组织在MRI图像中有不同的亮度,分析利用亮度值并基于亮度值进行分类的方法有很多,如直方图、混合高斯模型以及FCM等。对于纹理信息,在MRI图像中,当前点的属性不仅由自身亮度决定,还由周围相邻像素决定。如果周围都是类别A,那么当前像素将有极大概率属于类别A。MRF是对图像中邻近像素间相互影响建模的强有力工具,它可以充分利用图像自身包含的纹理信息进行图像分割。
本文建立了一个结合FCM和MRF的MRI图像分割方法[17-22]:在最大后验概率马尔科夫随机场(maximum a posteriori Markov random field,MAP-MRF)理论框架下,定义空间约束势能函数,进一步引入FCM隶属度信息提高分割性能。此外也对原始图像进行多尺度分解,从其中的低频图像提取FCM隶属度信息,进一步减轻噪声干扰和提高鲁棒性。
MRF理论提供了一个便利和一致的方式使得能够对上下文依赖关系建模,如图像的相邻像素以及相关特征[23]。MAP-MRF框架由Geman等[24]提出用来解决图像统计分析问题。对于图像分割问题,可以看作为标签问题。通过求取标签问题的最大后验概率估计获得最终的分割结果。
对于二维图像(高度M、宽度N)定义二维网格:
对于二维网格S,元素i的邻域像素定义为以i为中心为半径的圆以内的像素。定义像素i的邻域像素集合Ni:
其中,pixeli和pixeli'分别表示图像中的第i和第i'个像素。
S中的基团c定义为S的子集。单点基团表示为c={i},双点基团表示为c={i,i'},三点基团表示为c={i,i',i''}(i,i',i''∈S,且彼此相邻)。所有单点基团、双点基团和三点基团的集合表示为C1、C2和C3,所有的基团为 C=C1∪C3∪C3…
在MAP-MRF框架内,MRI图像分割问题可以转换为标签问题。如在脑组织MRI图像分割中将白质、灰质、脑脊液等分别赋予不同的标签。MRI图像中的每个像素标签的联合概率通过MRF建模,联合概率表示为
其中,fi(fi=fM×N)为第i个像素在相应标签下的概率。依据Hammersley-Clifford理论:
其中,Z和T是常数;d(d=d1,d2,…,dM×N)为观测值;U(f|d)为能量函数,表示为
能量函数是基团势能函数Vc(f|d)在c中所有可能的基团的和。由此,原始的图像分割问题或标签问题表示为
其中,P(f)为图像中的每个像素标签的联合概率,U(f)为能量函数,f为分割标签,argmax为取最大值,argmin为取最小值。
如果在所有基团上进行计算求解标签,计算量将十分巨大。实际上,在大多应用中仅在单点基团和双点基团上进行计算。或者说,通常采用二阶邻域系统,像素i的标签仅由像素i自身及其邻域像素集合Ni={pixeli'∈S|[dist(pixeli',pixeli)]2≤2,i'≠i}决定。这时,上述模型也称作auto-model[20]。
本文中,定义MAP-MRF下的势能函数:
其中,VC1(fi,di)是在单点基团上的势能函数,VC2(fi,fi')是在双点基团上的势能函数,所有相邻像素点都包含在双点基团中并被计算势能函数。
1.1.1 单点基团势能函数
单点基团势能函数从像素i自身为其分割提供支持信息,主要是利用图像的灰度信息为分割提供支撑。本文中,主要是通过FCM隶属度从像素i提取支持信息。本文中单点基团势能函数定义为
其中,λ1和λ2是系数,Ucoarse为隶属度矩阵。
(1)在原始图像上的单点基团势能函数。
单点基团势能函数定义为
其中,ui,j是像素i属于类别j的FCM隶属度。Lk是类别标签,k={1,2,3,4},L={BG,GM,WM,CSF},BG、GM、WM和CSF分别表示背景、灰质、白质和脑脊液。
(2)在低频图像上的单点基团势能函数。
对原始图像进行小波变换并提取低频图像,在该低频图像上应用FCM产生另外一个FCM隶属度矩阵Ucoarse,以进一步抑制噪声带来的影响。
其中,μcoarse是模糊隶属度矩阵Ucoarse中的元素。
1.1.2 双点基团势能函数
双点基团势能函数从像素i的相邻像素为像素i的分割提供决策支持信息,以充分利用图像的空间信息为分割提供支撑。本文中双点基团势能函数定义为
其中,λ3和 λ4是系数。
(1)基于相邻像素FCM隶属度的势能函数。
基于空间一致性假设,每个像素值都应该与其邻近像素的值接近,如当前像素属于某一类别,其隶属度也应该与它相邻像素的隶属度一致。定义相邻像素隶属度的势能函数为
其中,Ni是像素i的邻域像素集;abs(·)是求绝对值函数;di,i'是像素i和它的邻域像素i'间的距离,并且作为权重反映了不同距离的邻域像素带来的不同影响。
(2)针对邻域像素类别的势能函数。
当前像素i最有可能与其相邻像素属于同一类别,或者说像素i有最大概率与它的相邻像素属于同一分割区域,因此定义针对邻域像素类别的势能函数为
其中,Ii是像素i的灰度值,Li是像素i的标签(或类别)。
在家兔脑部MRI成像过程中,因为口腔、食道里的唾液以及其他颅内在MRI图像上成像灰度接近脑脊液的组织会使得分割结果中包含这些组织对应的区域,从而影响分割准确度。在通过1.1章节所述的图像分割算法获得潜在脑脊液区域后,需要再与预先建立的家兔脑脊液模板图像进行掩码处理以获取精确的脑脊液分割结果。
本研究所用的家兔脑脊液模板库是先由图像自动分割提取脑脊液,然后由专家人工修正并划定脑脊液大致存在区域得到的。模板库中,每一个切片对应一个原始图像IT和一个脑脊液区域划定好的图像IM。
图像配准和掩码处理过程:首先,将待处理图像IO与模板库中对应的切片图像IT配准,并将配准参数θ(R)保存。图像配准有刚性配准和非刚性配准,在配准精度上有的已经发展到亚像素级配准,但这里所用到的配准可以只涉及刚性配准,并且对配准精度要求不高。然后,将待处理图像利用1.1章节所述的图像分割算法进行分割,得到所有的候选脑脊液区域IP。最后,将IP与和IT对应的IM结合配准参数θ(R)进行掩码,得到仅包含脑脊液的图像IC。图像配准和掩码处理流程如图2所示。
图2 图像配准和掩码处理流程
目前,在针对人脑的MRI图像自动分割算法中,通常是将图像自动分割为背景、白质、灰质和脑脊液4个类别。但本课题组发现针对兔脑的MRI图像自动分割和脑脊液提取实验分为4个类别并不能得到最准确的脑脊液提取结果,因此设计了FCM分割类别自动选取方法,如图3所示。
图3 自动确定类别数量流程
实验步骤:首先,制作家兔脑出血模型,并在脑出血过程中进行MRI成像,获得脑出血不同阶段的家兔脑组织MRI图像;其次,对MRI图像分割得到家兔脑脊液区域,并结合MRI成像参数计算此时脑脊液的体积;最后,结合脑出血不同阶段家兔脑脊液体积的量化值,实现观察家兔脑脊液随着脑出血进程发展的变化。
家兔脑出血模型[25]采用自体血注入法[19]。首先用乌拉坦(25%,5 ml/kg)通过耳缘静脉注射麻醉。其次用立体定位仪确定大脑内囊出血位置,沿头骨中线做一个小的皮肤切口确定立体定位位置,以前卤作为参考点,冠状平面中轴线左右1 mm处即为内囊所在点。然后进行穿刺,穿刺点在矢状缝右侧6 mm,冠状缝偏下1 mm的位置,使用直径1 mm的钻头钻入颅骨,埋置直径0.7 mm的注射针插入约13 mm,而后用牙科水泥予以封闭确保颅腔保持密封状态;再用家兔的自体动脉血缓缓注入到实验组家兔的大脑中。最后将麻醉后的兔子固定于实验板上,各时间点分别做可变反转角快速自旋回波(sampling perfection with application optimized contrasts using different flip angle evolutions,SPACE)序列扫描,按照 3 ml/54 min的速度进行不间断注血,同时进行整个头颅扫描,扫描采用矢状位和俯卧位进行分层扫描。
本研究中家兔脑出血MIPS监测和MRI观测实验中所用家兔由陆军军医大学实验动物中心提供。所选家兔体质量(2.2±0.2)kg,颅径偏差范围±5%,实验室温度、湿度等环境条件严格控制。通过实验,获得了8只家兔的脑组织影像以及脑脊液量化分析结果。
家兔脑组织MRI成像在陆军军医大学第一附属医院实现,实验机器为西门子3.0T MRI扫描仪,线圈为18膝关节线圈,采用SPACE序列。SPACE序列是T2权重快速自旋回波序列,适合于获得脑脊液与其他相邻脑组织的清晰边界。MRI扫描仪的参数设置:重复时间(repetition time,TR)1 300 ms,回波链长度(echo train length,ETL)49,回波时间(echo time,TE)44 ms,矩阵 320×275,视野(field of view,FOV)160mm×160mm,层数192,层厚0.5mm,层间距0 mm。
图4为其中一只家兔MRI影像矢状方向第95帧图像的手动分割结果和本文所用算法分割结果,2种分割结果的脑脊液分割区域和量化分析值都非常接近。在实验中,依据手动分割结果对脑脊液体积进行量化分析,得到脑脊液体积约为1.18 ml;依据图像自动分割算法和量化计算对脑脊液体积进行量化分析,得到脑脊液体积约为1.10 ml(颅外椎管中的脑脊液未计算)。手动分割和图像自动分割算法的误差约为6.8%。
图4 家兔脑脊液量化计算实验
图5为第一只脑出血家兔针对脑脊液MRI成像后第95帧图像。表1为从8只家兔中得到的脑脊液量化分析结果。图6则显示了以8只家兔平均值反映的脑脊液体积随着出血量增加的变化趋势。其中兔子4和兔子5相比其他兔子变化趋势更缓慢,但随着注血量增加脑脊液体积减小趋势仍然存在。可能原因有两方面:(1)个体差异,兔子4和兔子5身体比较健壮,代偿能力强;(2)可能是由实验造成的,注血管插入太深,血液进入基底节并弥散开来,引起测量脑脊液体积减小。
图5 脑脊液观测实验结果
家兔的动物模型在脑出血和脑缺血等大脑医学研究中有非常广泛的应用。在这些研究中,往往对动物模型中脑脊液的变化十分关注[26]。当前针对脑脊液的医学影像分割方法主要集中在人脑,而对家兔脑脊液的分割提取还缺乏比较有效的方法。
本研究利用MRI成像和图像处理算法对家兔脑脊液进行了量化分析。
本研究结合FCM和MRF建立了脑脊液分割算法。该算法通过FCM提取利用图像中包含的灰度信息,通过MRF提取利用图像中包含的像素间空间信息。在充分利用灰度和空间信息的基础上,与预先建立的家兔脑脊液模板进行掩码处理,实现了对脑脊液的准确自动分割。
为实现定量观察和分析家兔脑出血过程中脑脊液的变化,以8只家兔脑出血模型为研究对象,以本文图像分割算法为手段,结合MRI成像参数计算家兔在不同出血阶段的脑脊液体积。在图6实验结果中,可以显著观察到脑脊液随着家兔脑出血量的增加不断减少,本文通过影像学手段观察了脑脊液的代偿过程,随着颅内出血的增加,脑脊液主动减少分泌增加吸收,达到减少颅内脑脊液体积以维持正常颅内压的目的;并且由图6还可以看到,当家兔颅内出血2 ml后,脑脊液体积变化的曲线变得平坦,即脑脊液减少速率变慢,这正说明2 ml前脑脊液处于正常的代偿阶段,2 ml后脑脊液代偿能力接近末期,大脑已经难以通过减少脑脊液体积维持颅内压。
表1 8只脑出血家兔的脑脊液量化计算结果
图6 脑脊液体积随着出血量增加变化趋势
本文通过利用MRI对家兔脑出血模型在脑出血过程中脑脊液的变化进行了定量分析。首先结合FCM和MRF建立了MRI图像分割方法,并结合掩码处理实现对家兔脑脊液的准确分割;基于脑脊液分割结果,结合MRI成像参数实现对脑脊液的定量分析,得到了随着脑出血的增加,家兔脑脊液的变化量化值,并观察到脑出血时脑脊液是不断减少的;通过定量实验说明了脑脊液代偿机制,为脑出血时颅内病生理变化研究提供了有效手段,为本课题组成功实现基于磁感应相移非接触检测家兔脑出血状况提供了有力的数据支撑。