从MRI体数据提取大脑表面形态的新方法

从MRI体数据提取大脑表面形态的新方法北京生物医学工程2000年第4期第19卷论著作者：罗述谦闫华王萍单位：首都医科大学生物医学工程系100054关键词：MR图像；皮层与颅骨剔除；表面重建；形态；轮廓摘要人脑形态研究在许多领域都有日益广泛的应用。将人的大脑从MR图像体数据中提取出来是脑的可视化、形态学分析、脑图象配准等多种工作的重要基础。本文结合解剖学知识、用基于灰度梯度的极坐标搜索方法从人脑MR图像自动剔除大脑皮层和颅骨及其它非脑成分，仅保留大脑部分。此后，对于3D图像遍历，我们使用移动立方体法根据全部256种构型，通过查表法把位于大脑表面上的像素与大脑内部及图像背景分开，得到准确的大脑表面轮廓。然后计算该表面轮廓的各个小三角形面片法线，结合适当的光照模型实现大脑表面形态的真实再现。大脑的提取工作是在轴向的方向完成的。但与轴向图正交的冠状图和矢状图显示结果表明，本算法对大脑部分的提取无论从哪一个方向来看都是很成功的。AMethodofExtractionofBrainSurfaceMorphologyfromMRVolumeDataLuoShuqian，YanHua，WangPing(DepartmentofBiomedicalEngineering,CapitalUniversityofMedicalSciences,Beijing1000054)AbstractThestudyofmorphologyofhumnbrainhasincreasingapplicationsinmanyareas.ExtractionofcerebrumfromrawMRIvolumedataisanimportantworkforvisualizationofhumanbrain,analysisofmorphologyandimageregistration.Usingpriorsofanatomy,agreylevelgradientbasedsearchingmethodinapolarcoordinatesystemtoremovecortex,skullandotherno-cerebrumstructuresautomaticallyfromMRimageofhumanbrainisdescribedinthispaper.Then,For3Dimageroaming,amothodofcubicdisplacementwasusedbasedonall256possibleconfigurations,alook-uptableisusedtogettheaccuratecontourofbrainsurfacebydistinguishingpixelsofbrainsurfacefrompixelsofinternalbrainandthebackground.Then,normalsofallthetrianglefacetsofthesurfacecontourarecalculated,tore-displaythesurfacemorphology,combiningwithappropriateilluminationmodel.Althoughtheextractionofcerebrumisimplementedintransverseslices,thereconstructedcoronalandsagittalslicesshowtheworkisverysuccesful.Keywords:MRimage,Removalofcortexandskull,Surfacerendering,Morphology,Contour0引言脑表面形态的提取是脑的可视化、形态学分析、脑图像配准等多种工作的基础。神经解剖学表明，只要我们能够知道大脑表面各沟回(或皱折)的形状、大小和空间位置，就可以对脑内各个次级脑组织准确定位。显然，这对用医学图像做神经外科手术的计划和引导有着十分重要的意义。多模医学图像信息的融合、fMRI功能映射研究都要求直接显示大脑皮层的形态。最近，在认知科学研究中关于人类智能与大脑形态的关系提出很有趣的假设。有人将爱因斯坦的大脑与常人对照组比较(男35人，女56人)，发现爱因斯坦大脑下顶叶比常人宽15%，而且回间沟没有延续到顶叶。至于这是不是天才大脑特有的形态，一些学者建议应将爱因斯坦的大脑与活在世上的著名科学家、学者的大脑进行比较。如果我们能将活人做MR全脑扫描，经3D重建得到准确的大脑表面形态。就可以进行这项研究了。可见，人脑形态研究在许多领域都有重要应用。已有许多方法被用来解决这个问题。在这里，我们用一种基于一定的解剖知识［1，2］和灰度梯度的方法来提取脑组织，简单易行，取得了较好的效果。1方法1.1MR图像中大脑的提取以184×217×161的T1脑图像数据作为处理对象，其平行于XY平面的层面图像依次标记为T1—001到T1—161(从颅顶到颅底)。基本处理步骤为：(1)灰度直方图选取平行于XY平面的几个不同层面的图像作灰度直方图。应选取不太靠近颅底或颅顶的尽量含较多大脑组织的层面，以使直方图中灰、白质灰度的参与统计的数量多一些。在灰度为［30，200］之间自动检出两个最高峰，在峰值附近一段区域，分别为灰质和白质的灰度范围，以此初步确定下面处理中用到的灰度阈值。(2)去除颅骨、头皮及脂肪对每个平行于XY平面的层片，以图像坐标的中心为极坐标系的原点，从一足够大的半径外(头部之外低亮度处)逐渐向原点靠近，找到大脑表面点即停止，再沿逆时针转一小角度重复上述操作，直到遍历一周回到出发时的角度(图1)。在由外向内搜索的过程中，首先遇到的高亮度区判定为头皮及皮下脂肪层，继而遇到的低亮度区为颅骨，再次遇到的高亮度点即判定为脑表面点(一般为灰质)。为避免将板障误判为脑组织并抑制噪声，计算梯度变化时需要考虑邻近点的灰度，例如3×3邻域或径向相邻6点的邻域。将搜索到的脑表面点外的象素灰度置为0，即完成了剔除颅骨、头皮及皮下脂肪的工作。图1左图：原始的T1图像，右图：剔除头皮和颅骨算法示意图1初始轮郭，2头皮，3颅骨，4大脑(3)去除眼球对于有眼球组织的层片，从图中读出左右眼球中心坐标，与上述方法类似，只不过从内向外进行搜索，遇到低亮度象素点即停止，再转一小角度重复操作。将搜索到的边界点内部的象素灰度置为0。此过程中同样考虑邻近象素的灰度。(4)去除小脑小脑部分在T1—044到T1—077之间，形状、大小变化很大，比较难以准确去除。T1—077到T1—054中，小脑基本上包含在一近似椭圆的区域内，可以应用同一种方法处理。T1—077经颅骨去除后，在小脑区作一椭圆(大小需手动调整，尽量使其包含该层片所有小脑部分)，椭圆内的象素灰度置为零。其相邻层片T1—076中小脑部分会略大一些，将与T1—077中相同大小与位置的椭圆放置于其中，让椭圆向外膨胀，直至遇到低亮度区为止。类似地，依次处理较小序号的层片，直至T1—054。T1—044到T1—053中小脑与大脑的分界基本上为左右两条直线，可以应用另一种方法。以小脑中心为极坐标原点，在一定角度内，由内向外搜索到低亮度点。左右两边的点分别连成直线，将直线间的象素灰度置为0。(5)脑干的处理脑干也属于脑组织，应保留。其灰度变化不大，形状比较简单。T1—001到T1—053经去颅骨的处理后，仍由内向外，遇到低亮度区为止，将检测到的点连成封闭曲线，以曲线为边界作一掩模，使得在作小脑和其余非脑部分的去除过程中，掩模区域内不受影响。(6)其余非脑部分的去除对含脑组织少并且左右半脑分开的层片，分别以左右半脑的中心为极坐标的原点，仍按梯度法搜索脑表面。将得到的点分别连成两条封闭曲线，曲线外的区域(除了脑干掩模内部分)象素灰度全置为0。1.2提取轮廓的算法在剔除大脑皮层、颅骨和其它非脑成分之后，仅剩下大脑部分。由于我们感兴趣的是脑表面的形态而不考虑其内部的细节，因此，要把位于大脑表面上的像素与大脑内部分开，这个过程称做轮廓提取(contouring)。在3D体数据集中，所有的采样点都位于一个立体栅格系统中。其最小的单元(cell)是以8个相临顶点构成的立方体(cube)。3D图像的像素称做体素(voxel)。一个体素可以由一个或多个这样的单元组成。物体的表面实际上是一个闭合的灰度的等值面，其灰度值称做阈值。在该等值面的内部，所有的象素灰度值都大于这个阈值，在等值面的外部，所有的象素灰度值都小于这个阈值(或相反)，从而将物体与背景分开。显然，等值面上的体素内部灰度是不均匀的，即体素的一部分灰度大于这个阈值，另一部分灰度小于这个阈值。下面让我们以一个最简单的例子说明如何寻找物体和背景的边界的。假设图像的体素仅由一个单元构成。先从寻找2D图像轮廓线说起。一个单元有4个顶点。每两个顶点联接成一条边。每个顶点的灰度就是该数据点的数值。假设我们选取了一个灰度阈值G，根据各顶点的灰度与阈值G的关系这些顶点被分作两类，分别用黑、白两色圆点表示(记作1和0)。对正方形的四个边逐个判断，如果某一个边的两个顶点颜色相同，该边上不存在边缘点；否则，在此边上必有一个边缘点。用直线将不同边上的边缘点连接起来。这些连接线将正方形分割为两或三部分。对于2D图像的像素，共有16种构型(case)，如图2所示。图22D图像边缘点的16种构型这里有几点需要说明：(1)上述各构型只说明边缘线与哪几条边相交，并没有指明交点的具体位置。交点的位置应通过对该边的两个端点线性内插来实现。(2)如果我们遍历图像中所有的小正方形，并对公共边合并，就可以得到图像中物体的轮廓线了。(3)在某些构型中(例如case5和case10)，对边缘点的连接有两种不同的方法。这种联接的不确定性称做构型的二义性。上述方法可以直接推广至三维图像。这时，每个像素有8个顶点。根据这8个顶点与灰度阈值的关系一共有28=256种构型。2D图像的轮廓是由直线段连接而成，3D图像的轮廓则复杂的多。因为，3D图像的轮廓是由许许多多的小三角形面片镶嵌而成的。考虑到各构型的对称和互补性，图3给出简化后的15种基本构型。对于3D图像遍历，根据各体素的构型情况产生三角形面片镶嵌的表面轮廓的方法称作移动立方体法［4］。实际应用中要用到全部256种构型，因为仅靠15种基本构型的组合往往会在表面轮廓上产生空洞。为了方便起见，实用的遍历法是对每个体素用查表法。将体素的8个顶点与灰度阈值比较所产生的逻辑值依序构成一个8位的二进制编码索引值，全部256种构型的信息组成一个“构型—三角剖分”查找表。它包含256个索引项，每个索引项包含索引号以及指向该种三角剖分中的一个指针。通过查表可以直接得到轮廓段的拓扑信息、哪一个边与体素相交、应当使用那些顶点内插产生交点等。对于每个体素，根据它的索引号在“构型—三角剖分”查找表中确定其三角剖分形式。还要对相临正方形一致边合并。最终产生由小三角形面片镶嵌成的表面轮廓。图33D图像的15种基本构型1.3等值面的浓淡显示要想真实地显示物体表面的情况，我们采用等值面的浓淡显示。三角片的生成仅仅完成了等值面的构造，要真正显示出物体在一定光照条件的形态，还必须解决物体在特定的光照模型下的表面法向量的计算［3］。(1)光照模型这里采用的光照模型为I=Ia+(Is-Ia)*cos(θ)I：三角片的光强Ia：环境的光强Is：光源的光强θ：三角片指向物体外部的法向量与光线的夹角显然，三角片的光强与光源的方向和强度均有关。图4光照模型三角片的表面法向量的计算是真实、准确显示物体表面的关键问题。(2)表面法向量的计算基于灰度梯度的法向量估计方法是一种很有效的方法，首先，用灰度差分计算体素顶点(i,j,k)上的灰度梯度g=(gx,gy,gz)，其中gx=［s(i+1,j,k)-s(i-1,j,k)］/2gy=［s(i,j+1,k)-s(i,j-1,k)］/2gz=［s(i,j,k+1)-s(i,j,k-1)］/2对g进行归一化，得到(gx/｜g｜，gy/｜g｜，