CN1296874C - 基于点重建的超大规模医学影像三维可视化方法 - Google Patents

Abstract

本发明涉及计算机图形学和模式识别技术领域，特别涉及一种利用点的重建来进行超大规模医学影像三维可视化的方法。其包括：a.分割步骤，从医学影像的二维切片中将感兴趣的部分分割出来；b.基于点的重建步骤，使用基于点的重建方法提取并重建出感兴趣的器官的三维表面，并以点的形式保存起来；c.三维绘制步骤，对提取出来的基于点的三维表面模型使用基于点的绘制方法来进行交互式地显示。本发明的方法，在配有一块中档显卡的普通的PC(个人电脑)上，实现与虚拟人体数据集相当规模的超大规模数据集的快速三维重建与实时显示。本发明具有高可信度、可应用性和可采纳性，并且在医学领域具有重要的应用价值。

Description

基于点重建的超大规模医学影像三维可视化方法

技术领域

本发明涉及计算机图形学和模式识别技术领域，特别涉及一种利用点的重建来进行超大规模医学影像三维可视化的方法。

背景技术

自从X射线发明以来，CT(计算机断层成像)、MRI(核磁共振成像)、CR(计算机X线成像)、B超、电子内窥镜等现代医学影像设备先后出现，使传统的医学诊断方式发生了革命性的变化。随着现代计算机科学技术的发展，医学影像处理作为一门新兴的交叉学科也随之出现，给医学诊断带来了新的曙光。但是传统的看片方式都是二维的，需要经过训练的放射科医生来看片并做出判断，随着计算机可视化技术的发展，医学影像的三维可视化也变得可能。医生通过计算机的辅助，可以看到人体器官的逼真的三维图像，从而可以提高诊断的正确率。

目前在医学影像处理领域使用的比较流行的三维可视化方法可以分为两大类，一类是基于体素的(Voxel-based)，另外一类是基于三角面片的(Triangle-based)。

基于体素的方法是最早应用于医学影像领域的，它用体素作为最基本的单元。所谓体素，就是以原始数据集的每个网格点为中心的一个立方体，一个感兴趣的器官的表面可以用一套体素来表达，在绘制的时候将可见的体素表面绘制出来。由于这种方法直接将原始的网格点和待抽取的表面联系起来，并且体素的表达方法很简单有效，因此受到了很多研究者的重视。但是基于体素的方法必须使用纯软件实现，这使得它不能用来实时处理超大规模数据集。

基于三角面片的方法是从传统的图形学领域引进医学影像领域的，随着MC(Marching Cubes移动立方体)算法的巨大成功，它在医学影像领域也得到了越来越广泛的应用。但是，原始的MC算法存在三大严重问题：拓扑不一致、计算效率低以及输出的三角片太多，很大程度地限制了它在实际中的使用，并且它也不适应于超大规模数据集的实时可视化。

随着医学影像设备技术的不断进步，医学影像的空间解析度越来越高，所得到的切片张数也越来越多，给已有的三维可视化算法带来了严峻的挑战。特别是自从虚拟人体项目出现以来，美国的虚拟人体得到的数据量大于36GB，而中国的虚拟人体得到的数据量更是大于100GB，如此海量的数据使得三维的实时处理和显示变得更为困难，如果使用传统的可视化算法的话，将会得到上千万甚至上亿的体素或者三角面片，这在目前的情况下，即使是使用超级计算机也无法实时处理，从而开发新的能够高效处理海量、超大规模的医学影像数据集的方法显得尤为迫切。

发明内容

本发明的目的是充分利用现有的普通PC(个人电脑)的硬件能力，开发一种能够高效快速地处理超大规模医学影像数据集，从而能够处理目前越来越多的海量的医学影像数据，进而从数据中更好地得到有用的信息。

为了实现上述目的，本发明的技术解决方案是提供一种基于点重建的超大规模医学影像三维可视化方法，其包括：

(1)分割步骤，使用域值分割或者区域增长分割方法从医学影像的二维切片中将感兴趣的部分分割出来；

(2)基于点的重建步骤，又包括两个子步骤：单层表面跟踪和生成边界球分层；

a)单层表面跟踪步骤，使用单层表面跟踪的方法，将分割所得的感兴趣的表面提取出来，并以点云的形式表达；

b)生成边界球分层步骤，将上一步生成的点云存储到一棵树中，树中每个节点记录这个点的边界球信息；

(3)三维绘制步骤，对提取出来的基于点的三维表面模型，使用基于点的绘制方法来进行交互式地显示。

在上面所述方法中，第(2)步基于点的重建步骤，以点作为重建的基本单元，且每个点的半径是相等的。

在上面所述方法中，第(2)步中a)单层表面跟踪步骤，在进行单层表面跟踪时，对x坐标进行量化、编码，不对y、z坐标进行量化、编码计算。

在上面所述方法中，第(2)步中b)生成边界球分层步骤，在生成边界球分层时，使用八叉树组织节点信息。

在上面所述方法中，八叉树每个节点信息的量化编码是32位，分为两组，一组是x坐标组，一组是法向量坐标组，都为16位。

在上面所述方法中，x坐标组分为x坐标、树结构和法向量锥三个字段。x坐标被量化、编码成10位；在树结构中，编码成4位；法向量锥为2位。

在上面所述方法中，树结构使用3位去记录每个节点的子节点个数，使用1位去记录当前节点是否有子孙节点；法向量锥的2位，用来记录当前节点法向量锥的范围。

鉴于目前个人电脑(PC)的运算能力日益增强，并且适用于PC的主流显卡的绘制速度也在不断地增加，本发明的主要目标是在配有一块中档显卡(如NVidia GeForce 2 GTS)的普通的PC上(在我国普及率很高)实现与虚拟人体数据集相当规模的超大规模数据集的快速三维重建与实时显示。使用传统的重建方法无法处理如此海量的医学影像数据集，或者处理速度非常慢，且要在昂贵的工作站上进行。本发明具有高可信度、可应用性和可采纳性，并且在医学领域具有重要的应用价值。

附图说明

图1是基于点重建的超大规模医学影像三维可视化方法的构成图；

图2是基于点的重建步骤的方框图；

图3是表面体素的半径的选取的示意图；

图4是使用基于点的单层表面跟踪算法能正确处理的连接性示意图；

图5是使用基于点的单层表面跟踪算法能正确处理的连接性示意图；

图6是输出点云的数据结构的示意图；

图7是量化编码后节点的32bits的意义解释及其示意图；

图8是基于点的三维绘制的流程图；

图9是本方法的实验结果1：在虚拟人体数据集上重建出来的骨骼模型；

图10是本方法的实验结果2：在虚拟人体数据集上重建出来的皮肤模型。

具体实施方式

下面详细描述本发明的基于点重建的超大规模医学影像三维可视化方法。本实现方案由三个主要步骤组成，结构图可以参见图1。这三个步骤分别是：分割步骤、基于点的重建步骤和三维绘制步骤，下面逐一进行介绍。

步骤一：分割

这一步的目的是为基于点的重建算法做预处理，将目标物体从背景中分割出来，也称为二值化的过程。分割对于高质量的三维重建是至关重要的，因为它决定着最终显示出来的物体是否是我们感兴趣的器官。

分割方法有很多不同的种类，每一种类都适合于不同的来源影像。比如阈值分割对CT(计算机断层成像)比较有效，但是对于MRI(核磁共振成像)图像来讲，由于人体内部结构复杂、生物组织的蠕动和MRI(核磁共振)成像的特点，造成医学图像中目标物体不可避免的受到其它物体甚至是噪声的干扰，使得物体局部边缘特征模糊，用阈值分割就难以得到较好的效果。所以最好的方法就是将分割方法和三维重建方法结合起来，提供尽可能多的分割方法，针对不同的来源影像选用不同的分割方法，得到高精确度的分割结果以后，再应用表面重建方法。

在这里我们介绍两种比较实用的分割方法：域值法和区域增长法。阈值法的关键是阈值的选择，可以由用户选择区分背景与非背景的灰度阈值，也可用自动阈值法确定阈值。常见的自动阈值法有P-参数法，状态法，微分直方图法，判别分析法和可变阈值法。针对医学图像噪声多的特点，可以采用判别分析法。即在图像灰度值的直方图中，求得阈值t把灰度值的集合分成两组，使得两组得到最佳分离。最佳分离的标准是两组的平均值的方差和各组方差的比为最大。该方法在直方图中有两个波峰时，可作为状态法起作用；即使不存在波峰时也可求出阈值。设给定图像具有L级灰度值，阈值为k，k将图像的像素分成两组1，2。组1的像素数设为ω₁(k)，平均灰度值为M₁(k)，方差为σ₁(k)；组2的像素数设为ω₂(k)，平均灰度值为M₂(k)，方差为σ₂(k)。设全体像素的平均灰度值定为Mτ。则：

组内的方差 ${\mathrm{\σ}}_W^2={\mathrm{\ω}}_1{\mathrm{\σ}}_1^2+{\mathrm{\ω}}_2{\mathrm{\σ}}_2^2$

组间的方差 ${\mathrm{\σ}}_B^2={\mathrm{\ω}}_1{(M_1-\mathrm{M\τ})}^2+{\mathrm{\ω}}_2{(M_2-\mathrm{M\τ})}^2={\mathrm{\ω}}_1{\mathrm{\ω}}_2{(M_1-M_2)}^2$

最佳标准 ${\mathrm{\σ}}_B^2/{\mathrm{\σ}}_W^2$ 值为最大，即σ_B ²取最大值。

对于区域增长法，需要用户选择皮肤轮廓上的一个点作为种子点。区域生长的基本思想是将具有相似性质的象素集合起来构成区域，该方法需要先选取一个种子点，然后依次将种子象素周围的相似象素合并到种子象素所在的区域中。区域生长算法的研究重点一是特征度量和区域增长规则的设计，二是算法的高效性和准确性。我们使用对称区域增长算法，可以有效地弥补区域增长算法的两大弱点：对初始种子点的选择敏感，以及内存占用过多，而且对3D连接对象标记和删除空洞的算法效率高。

步骤二：基于点的重建

基于点的重建是本发明中最重要的一部分，是考虑到医学影像数据集的特点，比如：均匀采样，无需颜色和纹理信息等以后，结合图形学中基于点的绘制算法，提出的一个新的方法。本方法使用点作为表面重建的基本单元，可以省去传统的重建算法中非常耗时的插值运算。同时本发明的目标是在普通的PC(个人电脑)机上实现对海量数据集的快速重建和实时绘制，除了提出基于点的重建算法以求在算法层次上提高速度以外，还利用了现代CPU(中央处理单元)提供的多媒体指令集(如MMX和SSE)，以及新近的显卡中的GPU(图形处理单元)提供的编程能力，来进一步在实现层次上提高整个方法的速度。

基于点的重建的方框图如图2所示，总共包括两个部分：单层表面跟踪和生成边界球分层。

第一部分：单层表面跟踪

在单层表面跟踪部分，主要的目的是从输入的分割后的二值图像中快速地将感兴趣的器官的三维表面形状以点的形式提取出来。这里使用的方法是首先遍历整个数据集，找出代提取的表面所经过的体素，这些体素被叫做表面体素，其它的体素被叫做背景体素。每个表面体素使用它的中心点坐标和一个半径来表示，为了避免在绘制时出现空洞，半径的大小必须仔细选择。这里考虑到医学影像数据是在规则的长方体网格上采样的，因此半径的选择就可以直接取长方体网格的外接球的半径，图3显示了一个两维的例子，来说明如何选择表面体素的半径。图中深颜色的格子代表表面体素，浅颜色的格子代表背景体素。对于表面体素，使用它的外接圆的半径来作为本体素的半径，图中的圆圈代表每个表面体素的外接圆。

通过使用点作为重建的基本单元，可以消除掉非常耗时的线性插值运算(如Marching Cubes移动立方体算法)或者三线性插值运算(如DividingCubes剖分立方体算法)，并且因为只需要计算中心点的坐标，可以完全使用整数运算来代替浮点数运算，加上通过使用CPU(中央处理单元)的多媒体指令集，可以使重建速度得到成倍的提升。

除了运算速度以外，为了从海量的数据集中快速有效地抽取表面，还有很多的因素需要综合考虑。第一必须考虑内存消耗问题，因为海量的数据集不可能一次同时调入内存，所以必须考虑如何最有效地将数据集来分层处理；第二必须考虑数据集的遍历效率问题，因为有上亿个体素需要遍历，其中有一大部分属于背景体素，如何快速地过滤掉这些背景体素也是一个需要认真考虑的问题。

不幸的是，内存消耗和遍历效率是两个相互冲突的因素。如果使用一些空间分割技术如八叉树等来加快体素的遍历速度，则除了数据集以外，在内存里还要存储八叉树等辅助结构，这对于大规模数据集显然是不实际的；如果要节省内存，则必然要顺序地遍历体素，这将降低重建的速度。

这里使用基于点的单层表面跟踪，可以在消耗较少内存的情况下，快速地重建出以点云形式表达的感兴趣的表面。算法在运行时，切片数据按照从下到上的顺序被分片读入内存，每两片切片组成了一层。算法在一层中只在四个方向进行表面跟踪，处理完一层后再读进下一片切片数据。这样可以既利用表面跟踪的快速的遍历速度，也可以节省内存。

单层表面跟踪与六个方向的表面跟踪相比，可能会导致一个问题：抽取出来的表面不完全。因为缺少两个方向的自由度，在三维空间中是相互连接的曲面可能在一张切片中不连接，如图4所示。

为了解决这个问题，我们在第一层的处理中采用顺序扫描，同时将在上方向上有连接的体素加入到一个种子点集合中，在处理下一层切片时，就从这个种子点集合出发，在四个方向上进行跟踪，同样的，记录在上方向上有连接的体素。通过这样的种子点传播，不仅提高了运算速度，也可以部分解决上述问题。在图4中，用粗线表示的和用细线表示的部分在大部分切片里都不连通，但是通过这种方法，最终连在了一个表面里。但尽管这样，仍然损失了一个向下的自由度，当第一层中不包括所有的种子点时，算法仍然可能只搜索出部分表面，如图5所示，只有用粗线表示的部分被抽取出来。尽管这时可以再向下搜索，但一部分切片数据就要被重复读取，同时整个算法的复杂度也大为增加。这里使用了一个简单的方法，设置一个三角面片个数域值，如果某一层里抽取出来的三角面片个数少于这个域值，则在此层中再进行顺序扫描，得到此层中完整的种子点，然后再向上传播。实践证明，如果域值选择合适，这样作是比较有效的。

从上面的描述可以看出，使用单层表面跟踪可以很大程度上地节约内存，事实上，考虑到计算梯度的需要，我们在内存里面一次容纳四张切片就足够了。尽管这样，但是考虑到成亿个要处理的体素，还必须仔细考虑表面跟踪的输出数据结构怎么组织。在传统的表达方式中，提取出来的点云中的每个点都被添加进一个顶点列表中，并且一个顶点使用6个浮点数来记录，三个记录x、y、z坐标，三个记录法向量的三个分量。对于非常大量的点来说，这种表达方式非常耗费内存，尤其是在只有有限内存资源的普通的PC(个人电脑)机上，这种表达方式是不可行的。本方法使用了一种紧凑的数据结构来表达抽取出的表面上的点的信息，请参看图6，图中，左边的方块代表一个数据集，每一行代表一个扫描行，中间示意出了一个扫描行包括很多个点，最右边是一个点的数据结构。这种数据结构类似于基于扫描线的RLE(行程编码)编码。假设数据集的三个方向上的维数分别为：Ix，Iy和Iz，首先分配一个数组，其大小为Ix×Iz，每个元素是一个指针，代表了一个扫描行的信息。这个指针指向另外一个链表，链表中记录了每个表面体素所抽取出来的点的坐标信息和法向量信息。每个点的信息被定义在一个叫做PointList的结构中，如下面的伪代码所示：

struct PointList

{

unsigned short xPos；       //x坐标：16位；

unsigned short normal；     //编码后的法向量：16位；

}

为了后续的生成边界球分层能够更有效地进行，每个点在向链表里面插入的时候，必须按照x坐标升序排列。通过使用这样的数据结构，本方法可以仅仅记录每个点的x坐标，因为其它两个坐标y和z都可以在数组寻址的时候隐式地得到。在这里，因为x坐标是记录的中点的坐标，通过简单的放大10倍或者记录一个0.5的偏差，就可以完全使用整数来记录坐标，而不需要浮点数。并且对于目前来说最大的数据集，16位来表达坐标也已经足够用了。对于法向量来说，为了节省空间，必须对其进行编码和量化。本方法中法向量在一个单位立方体上被量化，每个面被分成100×100个点，总共6×100×100个不同的法向量，被量化成16位，这样的精度并不会造成视觉上图像质量的下降。

这里再从理论上分析一下使用这种紧凑的数据结构所带来的好处。通过使用上述的数据结构，对于一个具有M个点的点云数据来说，需要使用(Ix×Iz×4+M×4)字节，其中Ix×Iz×4是指针数组所额外耗费的内存空间。如果使用传统的表达方式，将要使用M×6×4字节。对于海量的医学数据集来说，M往往能达到上千万，甚至上亿，而Ix×Iz相对来说小得多，所以通常能成倍的节省内存空间。

第二部分：生成边界球分层

因为后续要使用基于点的绘制来得到交互式的显示速度，在这里必须将上面一步得到的点云数据组织成分层的数据结构。在传统的基于点的绘制当中，一般情况下是将点云组织成分层的边界球，放在一个四叉树中。当四叉树建立后，每个节点的属性(坐标、半径)被量化、编码成32位。为了处理超大规模的数据量，整棵四叉树作为一个连续的流被写入磁盘中，以备后续绘制时使用。

为了快速、有效地建立分层的边界球表示，本发明再次利用了医学影像的特有的优点。考虑到点云中每个点的半径是相等的，这里没有使用四叉树，而是使用了八叉树，这样就可以不必每个节点都存储半径信息。事实上，如果知道一个节点所在的层，就可以根据八叉树的性质和每一层节点的半径相等这个特性来计算出这个节点的半径。并且通过使用上面描述的PointList数据结构，也可以简化八叉树的生成过程。因为x坐标是按照升序排列的，并且指针数组允许随机访问它的元素，那么在构造八叉树的过程中将一个节点分裂成子节点时所需要的信息能很容易地从这些数据结构中得到。另外，本方法仅仅需要量化、编码x坐标，y、z坐标可以从数据结构中直接得到，这也节省了很大的计算量。

本方法将每个节点的信息量化并编码成32位，图7给出了具体每一位的含义，其中，下面的汉字代表每个字段名，上面相对应的数字代表此字段占多少位。

法向量是直接从PointList结构中的normal字段得来的，x坐标是从PointList结构中的xPos字段得来的，并且被量化、编码成10位。在树结构中，使用3位去记录每个节点的子节点个数，使用1位去记录当前节点是否有子孙节点。另外还有两位用来记录当前节点法向量锥的范围。

步骤三：三维绘制

通过上面的两个步骤，已经得到了感兴趣的器官的表面的以点云形式表达的数据。三维绘制的目的就是快速地将这些数据显示在屏幕上，并可以让使用者来旋转模型、改变颜色等等，为了达到交互式地显示的目的，必须采用有效的算法。

在绘制阶段，本方法使用了上面一个步骤计算出来的分层边界球结构，也就是八叉树结构来进行可视性检测和层次细节控制。这里从顶向下遍历八叉树，在每个节点，如果其边界球是可视的，并且它的投影大小低于某个设定的阈值，那么这个节点就马上被绘制出来。为了维护比较好的实时交互，在使用者旋转模型时，设置一个比较高的阈值，以便得到比较高的帧速率；而在系统空闲时，逐步降低阈值直至屏幕的一个像素大小，从而得到最精细的图像。

这里本方法再次利用了医学影像的特点，因为八叉树中同一层的节点的半径大小是相同的，那么它们的投影大小也是相同的。因此在每一帧的开始，算法预先计算出根节点的投影大小，然后可以在后续的节点中重复利用此信息，而不必要对每一个节点都进行一层投影计算。

图8给出了整个绘制部分的流程图，为了最大限度地提高绘制速度，本方法还使用了Intel CPU(英特尔中央处理器)提供的SIMD(单指令多数据)多媒体指令集来处理视锥裁剪、矩阵矢量计算等，这些都充分地利用了目前个人用PC(电脑)  的硬件能力，来得到最高的性能。

实施例

我们将此方法运用到我们自行设计实现的三维医学影像处理与分析系统中。我们研制开发的基于微机的三维医学影像处理与分析系统3DMed是在微机Windows XP/2000/NT/98环境下，采用面向对象的设计方法和软件工程规范，用C++语言实现的、面向医学领域的三维影像处理与分析系统。本系统具有丰富的图形图像处理与分析功能，不仅具有完善的二维图像处理分析功能，而且具有强大的三维处理与分析、网络传输与存储等功能。系统提供的功能包括数据输入、影像数据管理、二维处理、三维数据处理、切片重组、三维显示、手术模拟、虚拟内窥镜、PACS(影像归档与通信系统)与远程诊断等一系列功能。

下面说明利用基于点重建的超大规模医学影像三维可视化方法来处理海量医学影像数据集的具体实施过程。实验数据为美国国立医学图书馆网站上提供的虚拟人体数据集，这里使用的是CT(计算机断层成像)数据集，数据集规模为512×512×1876×8位。使用的微机配备了一个PIII 866的CPU，256MB内存，显卡为NVidia GeForce 2GTS，配备32MB显存。

具体操作步骤如下：

1首先通过数据接口读取数据。

2点击“高级构造”按钮，进入分割界面。

3系统提供了多种分割方法可供选择，有种子生长、腐蚀膨胀、模糊连接、域值、交互式分割等，此时可以选取一种分割方法对数据进行分割。因为处理的是CT数据，域值方法是比较有效的，用鼠标指定好低域值和高域值以后，系统就会将在此域值之内的物质分割出来。

4分割以后，点“3D显示”按钮，系统就会调用本申请中所描述的算法进行三维重建，并将重建后的真实感图形显示出来，允许用户使用鼠标进行交互式观察。对虚拟人体数据我们分别重建了皮肤和骨骼两个模型，结果如图9和图10所示。图9和图10并没有对原始数据进行任何前处理，主要目的只是为了演示本方法能对海量数据进行快速的三维重建。图片上的一些非骨骼和皮肤的杂质是原始图像中的噪音，可以通过精细的手工分割去掉，在这里不再详细阐述。

5我们记录了基于点的重建和绘制所消耗的时间，对于本次实验所处理的数据，重建时间为：骨骼30.39秒，皮肤83.88秒；绘制时间为：骨骼1.62秒，皮肤3.13秒。而使用传统的方法，对于如此海量的数据集无法处理，只能通过磁盘交换来进行处理。我们实现了一个基于传统的方法进行重建的算法，并且使用磁盘空间来作为临时交换空间，得到的重建时间都在500秒以上，绘制时间在60秒以上。

上述实验结果与利用基于点重建的超大规模医学影像三维可视化方法来处理海量的医学影像数据集的理论分析结论一致，具有高可信度、可应用性和可采纳性。

Claims (4)

1.一种基于点重建的超大规模医学影像三维可视化方法，其特征在于，包括下列步骤：

(1)分割步骤，使用域值分割或者区域增长分割方法从医学影像的二维切片中将感兴趣的部分分割出来；

(2)基于点的重建步骤，以点作为重建的基本单元，且每个点的半

径是相等的；该重建步骤又包括两个子步骤：单层表面跟踪和生成边界球分层；

a)单层表面跟踪步骤，使用单层表面跟踪的方法，将分割所得的感兴趣的表面提取出来，并以点云的形式表达；

b)生成边界球分层步骤，将上一步生成的点云存储到一棵树中，树中每个节点记录这个点的边界球信息；

(3)三维绘制步骤，对提取出来的基于点的三维表面模型，使用基于点的绘制方法来进行交互式地显示。

2.按权利要求1所述的方法，其特征在于，所述第(2)步中a)单层表面跟踪步骤，在进行单层表面跟踪时，对x坐标进行量化、编码，不对y、z坐标进行量化、编码计算。

3.按权利要求1所述的方法，其特征在于，所述第(2)步中b)生成边界球分层步骤，在生成边界球分层时，使用八叉树组织节点信息。

4.按权利要求3所述的方法，其特征在于，所述节点信息，其每个节点信息的量化编码是32位，分为两组，一组是x坐标组，一组是法向量坐标组，都为16位。

Images