基于单层表面跟踪的超大规模医学图像表面重建方法

专利名称：基于单层表面跟踪的超大规模医学图像表面重建方法
技术领域：
本发明涉及模式识别，特别涉及基于单层表面跟踪的超大规模医学图像表面重建方法。
但是，通过医学影像设备所得到的影像都是二维的，需要经过训练的放射科医生来做出判断。随着计算机可视化技术的发展，通过计算机的辅助，可以将一系列的二维医学图像经过处理，生成逼真的三维图像，使医生看得更好，看得更准，看得更方便。
医学图像的三维可视化是该学科的核心组成，它利用一系列的二维切片图像进行三维模型重建和显示，是进行定量分析的前提。三维可视化实现中有两种绘制技术表面绘制和直接体绘制。面绘制最大的特点是需要先对二维数据场进感兴趣对象进行分割和三维重建，生成对象边缘等值的曲面表示，再采用光照模型绘制图像。而体绘制是将三维体数据中的体素看成一个个半透明物质，并分类赋予其一定的颜色和阻光度，由光线穿过整个数据场，进行颜色合成，得到最终的绘制结果。
由于表面显示速度快，并且能为医生提供更加逼真的人体组织与器官的交互式显示，从而提高医疗诊断的准确性，所以表面重建和表面显示技术在医学图像处理与分析中应用的越来越广泛。其中最为著名的就是87年由W.Lorensen提出来的MC(Marching Cubes)算法，由于MC算法实现简单，并且能得到高解析度的输出，还能充分利用当前主流的图形显示硬件，因此使得它成为最流行的表面抽取算法。但是，MC算法在表面重建时要生成大量的三角面片去逼近目标曲面，这使得要处理大规模模型比较困难，另外，MC算法的速度比较慢。尽管研究人员已经开发出了大量的方法来解决这些问题，但是这些方法或者需要额外的计算，或者需要大量的辅助存储空间，都不适应于超大规模数据集的实时处理。
目前已经有一些商用的软件，如Mayo Clinic的ANALYZE、SensorSystems的MEDx、Advanced Visual Systems的AVS以及枢法模·丹历的Stealth Station等，它们都提供了大量的工具用来进行2D、3D图像显示、图像处理、分割、配准、定量分析等，并且都提供了表面显示功能。但是由于三维重建的计算量很大，这些商业软件往往要求运行在高档的工作站上，并且价格昂贵，因此应用范围收到了很大的限制。
由于现在的医学影像设备所产生的图像分辨率越来越高(如1024×1024)，数据越来越多，给三维重建也带来了更大的挑战。特别是自从美国的虚拟人体项目出现以来，海量的数据使得实时处理和显示变得更为困难，从而使得开发新的能够高效处理超大规模医学图像数据集的方法尤为重要。
为实现上述目的，基于单层表面跟踪的超大规模医学图像表面重建方法包括(1)分割步骤，从医学图像的二维切片中将感兴趣的部分分割出来；(2)表面提取步骤，使用单层表面跟踪将感兴趣的器官表面提取出来；(3)三角带生成步骤，对提取出来的三维表面模型进行处理，使其适合于快速绘制；(4)交互式显示步骤，对器官的三维表面模型进行真实感显示和实时交互。
本发明的目标平台是在我国普及率很高的普通PC，操作系统是界面友好的Windows系列操作系统，不仅成本低廉，并且容易操作。通过利用单层的表面跟踪和几何数据压缩技术，大大减少了算法的内存消耗量，同时还保证了非常快的重建速度。在医学图像数据量越来越大的今天，该发明在医学领域具有重要的应用价值，并且具有高可信度、可应用性和可采纳性。
图3是使用基于单层表面跟踪算法能正确处理的连接性示意图；在图3中，用粗线表示的和用细线表示的是一个连续曲面的两个部分，但是它们在一层切片内大部分是不连通的。通过使用下面描述的域值方法，它们最终可以被连在一个表面里。
图4是查找表中表项是如何组织的示意图；在图4(a)中，举出了一个简单的体素配置的例子，(b)图是六个方向的示意图，(c)图给出了在查找表中它对应的表项的值。
图5是所用数据缓冲区的示意图；其中，图5(a)示意了三个缓冲区的位置(Top Edges，Bottom Edges和Z Edges)，(b)图左边示意了数据集的一个截面图，右边示意缓冲区中的内容如何向上移一层。
图6是一个三角面片网格及其邻接图的示意；图7是虚拟人体的骨骼重建；图8是皮肤重建效果图。
发明的实施方式目前个人电脑在我国不断普及，并且价格也不断在下降，而PC机的硬件功能却以指数级飞速地发展。新型的高速CPU不断出现，现在已经达到了GHz，处理速度越来越快，并且不断有新的指令集出现来加速多媒体应用程序的执行。支持三维硬件加速技术的显卡频繁更新换代，现在已经出现了比CPU处理速度更快的图形处理器(GPU)，如NVidia的GeForce2等，可以使用硬件来完成诸如三维图形变换、光照、裁剪等原先需要软件来完成的功能。所有这些条件都使得在当前的主流PC机上开发出实用的可以高效处理超大规模数据集的方法成为可能。
下面详细描述本发明的超大规模医学图像的表面重建与显示方法。本实现方案由四个主要步骤组成，结构图可以参见

图1。这四个步骤分别是分割步骤、表面提取步骤、Triangle Strips生成步骤和交互式显示步骤，下面逐一进行介绍。
分割这一步的目的是为表面重建算法做预处理，将目标物体从背景中分割出来，也称为二值化的过程。分割对于高质量的三维重建是至关重要的，因为它决定着最终显示出来的物体是否是我们感兴趣的器官。
分割方法有很多不同的种类，每一种类都适合于不同的来源影像。比如阈值分割对CT比较有效，但是对于MRI图像来讲，由于人体内部结构复杂、生物组织的蠕动和MRI成像的特点，造成医学图像中目标物体不可避免的受到其它物体甚至是噪声的干扰，使得物体局部边缘特征模糊，用阈值分割就难以得到较好的效果。所以最好的方法就是将分割方法和三维重建方法结合起来，提供尽可能多的分割方法，针对不同的来源影像选用不同的分割方法，得到高精确度的分割结果以后再应用表面重建方法。
在这里我们介绍两种比较实用的分割方法域值法和区域增长法。阈值法的关键是阈值的选择，可以由用户选择区分背景与非背景的灰度阈值，也可用自动阈值法确定阈值。常见的自动阈值法有P-参数法，状态法，微分直方图法，判别分析法和可变阈值法。针对医学图像噪声多的特点，可以采用判别分析法。即在图像灰度值的直方图中，求得阈值t把灰度值的集合分成两组，使得两组得到最佳分离。最佳分离的标准是两组的平均值的方差和各组方差的比为最大。该方法在直方图中有两个波峰时，可作为状态法起作用；即使不存在波峰时也可求出阈值。设给定图像具有L级灰度值，阈值为k，k将图像的像素分成两组1，2。组1的像素数设为ω1(k)，平均灰度值为M1(k)，方差为σ1(k))；组2的像素数设为ω2(k)，平均灰度值为M2(k)，方差为σ2(k)。设全体像素的平均灰度值定为Mτ。则组内的方差σw2=ω1σ12+ω2σ22]]>组间的方差σB2=ω1(M1-Mτ)2+ω2(M2-Mτ)2=ω1ω2(M1-M2)2]]>最佳标准 值为最大，即 取最大值。
对于区域增长法，需要用户选择皮肤轮廓上的一个点作为种子点。区域生长的基本思想是将具有相似性质的象素集合起来构成区域，该方法需要先选取一个种子点，然后依次将种子象素周围的相似象素合并到种子象素所在的区域中。区域生长算法的研究重点一是特征度量和区域增长规则的设计，二是算法的高效性和准确性。我们使用对称区域增长算法，可以有效地弥补区域增长算法的两大弱点对初始种子点的选择敏感，以及内存占用过多，而且对3D连接对象标记和删除空洞的算法效率高。
表面提取要从大规模数据集里快速地重建出感兴趣的三维表面，有很多因素必须认真地考虑1)内存消耗。因为我们的目标平台是普通的PC机，对于象虚拟人体这样大规模的数据集来说，内存的合理使用是很重要的。
2)遍历速度。因为要处理的体素个数太多，如何加速遍历体素的速度也是很关键的。
3)绘制速度。产生的大量的三角面片将给绘制速度带来很大的挑战。
但是很不幸，在这些因素中，内存的消耗和遍历的速度是一对矛盾体。如果使用一些空间分割技术如八叉树等来加快体素的遍历速度，则除了数据集以外，在内存里还要存储八叉树等辅助结构，这对于大规模数据集显然是不实际的；如果要节省内存，则必然要顺序地遍历体素，这将降低重建的速度。
这里所提出的算法综合考虑了以上因素，它可以在消耗较少内存的情况下，快速地重建出感兴趣的表面，并实时地显示出来。
因为在六个方向的表面跟踪使得数据缓冲变得困难，并且需要更多的辅助存储空间，所以这里提出了单层的表面跟踪技术，即只在切片所在平面上的四个方向进行表面跟踪，在与切片相垂直的方向上按照从下到上的顺序进行处理，我们的实验证明此算法的效率是相当高的。
在标准的Marching Cubes算法里，体素遍历的次序是从左到右，从下到上，不考虑体素的邻居信息，因此算法的大部分时间花在对空体素的检测上。如果使用传统的表面跟踪算法，则需要将整个数据集读入内存，并且需要对整个数据集进行随机访问，这限制了它在大规模数据集上的应用。
我们的思想是既利用表面跟踪的优点，又尽量少消耗内存，并且实现对数据集按照切片顺序访问，从而在速度和内存消耗上取得一个最好的平衡点。很自然地，我们将表面跟踪限制在单层内，从而达到了这一目的。
在这里提出来的算法中，切片数据按照从下到上的顺序被分片读入内存，每两片切片组成了一层。算法在一层中只在四个方向进行表面跟踪，处理完一层后再读进下一片切片数据。
单层表面跟踪与六个方向的表面跟踪相比，可能会导致一个问题抽取出来的表面不完全。因为缺少两个方向的自由度，在三维空间中是相互连接的曲面可能在一张切片中不连接，如图2所示。为了解决这个问题，我们在第一层的处理中采用顺序扫描，同时将在上方向上有连接的立方体加入到一个种子点集合中，在处理下一层切片时，就从这个种子点集合出发，在四个方向上进行跟踪，同样的，记录在上方向上有连接的立方体。通过这样的种子点传播，不仅提高了运算速度，也可以部分解决上述问题。在图2中，用粗线表示的和用细线表示的部分在大部分切片里都不连通，但是通过这种方法，最终连在了一个表面里。但尽管这样，仍然损失了一个向下的自由度，当第一层中不包括所有的种子点时，算法仍然可能只搜索出部分表面，如图3所示，只有用粗线表示的部分被抽取出来。尽管这时可以再向下搜索，但一部分切片数据就要被重复读取，同时整个算法的复杂度也大为增加。这里我们使用了一个简单的方法，设置一个三角面片个数域值，如果某一层里抽取出来的三角面片个数少于这个域值，则在此层中再进行顺序扫描，得到此层中完整的种子点，然后再向上传播。实践证明，如果域值选择合适，这样作是比较有效的。
为了快速在单层数据中进行表面跟踪，可以构造一个查找表NeighbourTable，它对256种可能的体素配置方式都用一个字节来记录与邻居体素的连接信息。字节的最高三位代表在此体素邻居体素中有效体素的个数，范围为0-6，低五位从高位到低位分别代表上、前、后、左、右五个方向上的邻居体素是否有效，其中1代表有效，0代表无效。例如，对于图4所示的体素，它所对应的NeighbourTable里的表项应该是01111001。
另外考虑到输出三角面片的个数很庞大，如果采用传统的顶点表和面表的网格表达方式的话，将会耗费大量的内存。这里对于顶点数据的组织，我们采取了一种比较紧凑的数据结构。如果数据集的解析度是Ix×Iy×Iz，并且假设最大的一维是Ix，首先分配Iy×Iz个指针PointsList*，其中PointsList结构组织如下x坐标16bit法向量16bit法向量在一个单位立方体上被量化，每个面被分成100×100个点，总共6×100×100个不同的法向量，被量化成16bit，这样的精度并不会造成视觉上图像质量的下降。
采用这种数据结构后，要保存网格里M个点的数据只需要Iy×Iz×4+M×4个字节，而如果使用传统的坐标法向量表示法的话就需要M×(3+3)×4个字节，从我们的观察看，在大部分情况下新的方法只占用传统方法的不到1/4的存储空间。
算法在对一层进行跟踪时，首先从上层得到的种子点集合中取出一个点，如果此点没有被访问过，则将其加入到一个队列的尾部。算法每次都去队列的头部取出一个体素，得到此体素内的等值面与各个边的交点(此处直接采用中点坐标)，并且得到三角片的连接方式。然后通过查找NeighbourTable表得到此体素在四个方向上的连接信息，对于相连接的邻居体素，如果其没有被访问过，则将它加入到队列的尾部。同时还要判断此体素在上方向是否有连接，如果有，记录进种子点集合中。循环执行上面的过程，直到队列为空时，单层中相连接的等值面就被抽取出来了。
在表面跟踪的过程中，可以利用数据缓冲机制来避免一些重复的计算。对于一个体素来说，它与每个邻居体素都要共享一个面，这些面上的等值点实际上只需计算一次，保存在一个数据缓冲区中，下次就可以直接从缓冲中将计算结果直接取出来。要实现一个比较好的数据缓冲算法的难点就在于如何尽量减少数据缓冲区所占用的内存，并且要在一个数据没有用的时候将其释放掉，供其它的数据使用。
由于本算法将表面跟踪限制在单层里面，所以可以很容易地实现一种有效的数据缓冲机制。在这里使用了三个缓冲区来分别记录顶层边(TopEdges)、下层边(Bottom Edges)和中间的边(Z Edges)，如图5所示。假设切片数据的高和宽分别为ImageHeight和ImageWidth，则顶层边和下层边缓冲区的长度都为2×ImageHeight×ImageWidth，而中间边只需要ImageHeight×ImageWidth的长度，这样不仅保证了对每条边都有一个哈希表项所对应，同时也没有任何内存上的浪费。
对于等值面的表达，这里采用了传统的顶点表和面表的表达方式。顶点表VertexTable记录了所有的顶点的坐标和法向量，它的表达形式如下所示{x，y，z；nx，ny，nz}其中x，y，z是顶点的三个坐标，nx，ny，nz是顶点在三个方向上的法向量。而面表FacetTablei与传统的表达方式有所不同，它只记录了一层中的等值面的连接方式，它的表达形式如下所示
{index1，index2，index3}其中index1，index2，index3分别代表一个三角面片的三个顶点在VertexTable中的位置索引。当一层的等值面被抽取出来后，FacetTablei就被加入到一个链表中去。
三个缓冲区(Top Edges、Bottom Edges和Z Edges)里面的所有表项在最初始状态下全部被置为Empty，当处理一个体素的某条边时，如果这条边上有等值点，那么首先去对应的缓冲区里查找对应的位置上的表项是否为Empty，如果为Empty，那么就要计算这个点的坐标、法向量，并且要插入到顶点表VertexTable中，这时还要同时更新对应的表项，记录这个顶点在VertexTable中的索引；如果表项不为Empty，表示这个顶点在前面已经计算过，可以直接取出里面记录的顶点索引，加入到FacetTablei中即可。
在计算完一层的等值面后，缓冲区里面的内容也应该向上推进一层，可以把缓冲区看作为一个窗口，这个窗口整个往上移了一格。首先，把缓冲区里面第二层的内容赋值给第一层；然后将第二层里的所有表项内容全部清为Empty，这样就重复利用了有限的数据缓冲的空间，如图5所示。
三角带(Triangle Strips)生成通过上面的两个步骤，已经得到了等值面的三角面片网格表达。因为三个顶点才能表达一个三角面片，所以要想表达具有n个三角面片的网格，就必须有3×n个顶点；而如果使用三角带(Triangle Strips)的话，因为三角面片之间共享了顶点，所以可以仅使用n+2个顶点就能表达n个三角面片。这样可以大量地节省内存和显示卡的内存之间传输的带宽，极大地提高显示速度。
为了更快速地绘制生成的网格，这里采用了一个快速的三角带生成算法，对提取出来的等值面生成更紧凑的三角带表达。由于上面采用了数据缓冲机制，从而保证在顶点表VertexTable中没有重复的顶点，并且在面表FacetTablei中记录的是顶点的索引值，这样算法可以仅仅在面表中进行操作。下面是算法的描述。
首先，算法由面表创建出一个邻接图。所谓邻接图，实际上就是一个大的哈希表，它记录了每个三角片在每条边上的邻居信息，图6显示了一个三角网格及它所对应的邻接图。要创建邻接图，在面表FacetTablei中查找所有共享一条边的三角面片，并将它们插入到邻接图中。
其次，得到邻接图以后，从图中找到具有最小度数(邻居数)的三角片。如果有好几个三角片都具有同样的度数，那么再查找它们的邻居的度数，选中邻居度数最小的那个三角片；如果它们的邻居的度数仍然一样的话，那可以随便选中一个。
然后，从选中的这个三角片开始，循环生成三角带。将被选中的三角片加入到三角带中，并且在它的邻居中查找具有最小度数的三角片，将其加到三角带中，一直这样循环，直到所有的邻居都被访问过为止。如果碰到一个三角片，它没有任何邻居，那么此时应该创建一个新的三角带；如果碰上了有几个相同度数的三角片时，处理方法与上面所介绍的相同。为了避免无限循环，在将一个三角片加入到三角带之后，需要在邻接图中将它和所有指向它的指针全部去掉。
最后，检查生成的三角带是否有效地组织了所有的顶点，是否是有效的三角带。这时要对一个三角带中的所有三角片检测，看它的最后两个顶点是否是下一个三角片的前两个顶点。
交互式显示在得到要提取器官的三维模型以后，如何将它们交互式地显示出来也是一个相当重要的问题。为了尽可能地挖掘普通PC机的潜能，实现高速度的真实感显示，我们考虑了如下因素1)使用成熟的三维图形API OpenGL，利用硬件加速OpenGL是目前比较成熟的三维图形API，它可以使用当前PC机上所大量采用的主流的显卡(如nVidia公司的GeForce系列)所提供的硬件加速技术，已经成为当前三维显示API方面的工业标准。
因为我们在上一个步骤已经将三维模型表达成了三角带(TriangleStrips)的形式，并且现在的主流显卡都专门针对三角带的显示作了优化，所以我们使用了nVidia的OpenGL扩展NV_vertex_array_range和NV_fence，从我们的实验来看，在一个PIII800的机器上，一秒钟大概能绘制700万左右个三角面片，能满足大部分大规模数据的实时显示要求。
2)用软件实现连续LoD控制，软硬件结合实现实时绘制对于虚拟人体这样的超大规模数据集，抽取出来的三维模型可能包含有上千万甚至上亿的三角面片，如果纯粹只用OpenGL的能力来进行绘制的话，还是得不到满意的结果。所以必须用软件来实现LoD控制和消隐，这将大大减轻显卡的GPU的工作负担，并且可以完全充分地利用CPU和GPU的运算能力，使得它们同时并行工作。
对于抽取出来的三角面片网格，将其组织进一个分层的边界球中。在绘制过程中，使用分层边界球来进行可视性判断、LoD控制。通过从上到下遍历分层边界球，如果一个边界球是可见的，并且投影到屏幕上的面积低于一个设定的域值，那么此边界球被绘制。在用户使用鼠标交互控制的时候，为了维护一个比较高的帧速率，可以将面积域值设得稍微高一点；在系统空闲的时候，可以将面积域值设为一个像素，得到最高质量的图像，因此可以很自然地实现LoD控制。
通过使用上面描述的软硬件混合绘制方法，要处理超大规模的数据集，可以在普通的PC机上就能取得很好的效果。
下面说明利用基于单层表面跟踪的方法来处理超大规模医学图像数据集的具体实施过程。实验数据为美国国立医学图书馆网站上提供的虚拟人体数据集，这里使用的是CT数据集，数据集规模为512×512×522。因为在数据获取时一些参数的差异(比如fov，像素大小，切片间隔等)，整个数据集被分成了九个部分，因为最后一个部分(脚部的数据)与上一个部分的数据之间的间隔太大，我们只用了前八个部分。具体操作步骤如下1)首先通过数据接口读取数据。
2)点击“高级构造”按钮，进入分割界面。
3)系统提供了多种分割方法可供选择，有种子生长、腐蚀膨胀、模糊连接、域值、交互式分割等，此时可以选取一种分割方法对数据进行分割。因为处理的是CT数据，域值方法是比较有效的，用鼠标指定好低域值和高域值以后，系统就会将在此域值之内的物质分割出来。
4)分割以后，点“3D显示”按钮，系统就会调用本申请中所描述的算法进行三维重建，并将重建后的真实感图形显示出来，允许用户使用鼠标进行交互式观察。对虚拟人体数据我们分别重建了皮肤和骨骼两个模型，效果图如图7和图8所示。
上述实验结果与发明人对利用基于单层的表面跟踪来处理超大规模的医学图像数据集的理论分析结论一致，具有高可信度、可应用性和可采纳性。
权利要求
1.一种基于单层表面跟踪的超大规模医学图像表面重建方法，包括(1)分割步骤，从医学图像的二维切片中将感兴趣的部分分割出来；(2)表面提取步骤，使用单层表面跟踪将感兴趣的器官表面提取出来；(3)三角带生成步骤，对提取出来的三维表面模型进行处理，使其适合于快速绘制；(4)交互式显示步骤，对器官的三维表面模型进行真实感显示和实时交互。
2.按权利要求1所述的方法，其特征在于所述的单层表面跟踪包括步骤在切片所在平面上的四个方向进行表面跟踪；在与切片相垂直的方向上从下到上的顺序进行处理。
3.按权利要求2所述的方法，其特征在于还包括步骤在第一层的处理中采用顺序扫描；同时将在上方向上有连接的立方体加入到一个种子点集合中；在处理下一个切片时，从该种子点集合出发，在四个方向上进行跟踪；记录在上方向上有连接的立方体。
4.按权利要求3所述的方法，其特征在于包括步骤设置三角面片个数阈值，如果某一层里抽取出来的三角面片个数少于阈值，则在此层中再进行顺序扫描，得到此层中完整的种子点，然后再向上传播。
5.按权利要求2所述的方法，其特征在于所述的表面跟踪包括设置缓冲区，用于分别记录顶层边、下层边和中间的边。
6.按权利要求2所述的方法，其特征在于所述的缓冲区有三个。
7.按权利要求2所述的方法，其特征在于所采用的网格压缩表达方法如果数据集的解析度是Ix×Iy×Iz，并假设最大的一维是Ix，则首先分配Iy×Iz个指针PointsList。
8.按权利要求7所述的方法，其特征在于所述的PointsList数据结构为x坐标16bit；法向量16bit。
全文摘要
一种基于单层表面跟踪的超大规模医学图像表面重建方法，包括(1)分割步骤，从医学图像的二维切片中将感兴趣的部分分割出来；(2)表面提取步骤，使用单层表面跟踪将感兴趣的器官表面提取出来；(3)三角带生成步骤，对提取出来的三维表面模型进行处理，使其适合于快速绘制；(4)交互式显示步骤，对器官的三维表面模型进行真实感显示和实时交互。本发明的平台是在我国普及率很高的普通PC，操作系统是界面友好的Windows系列操作系统，不仅成本低廉，并且容易操作。通过利用单层的表面跟踪和几何数据压缩技术，大大减少了算法的内存消耗量，同时保证了非常快的重建速度。在医学图像数据量越来越大的今天，本发明在医学领域具有重要的应用价值。
文档编号G06T17/05GK1430185SQ01138698
公开日2003年7月16日 申请日期2001年12月29日 优先权日2001年12月29日
发明者田捷, 蒋永实, 常红星, 张晓鹏 申请人:田捷, 蒋永实, 常红星, 张晓鹏