CN111105476A - 一种基于Marching Cubes的CT图像三维重建方法 - Google Patents

Abstract

本发明是一种基于Marching Cubes的CT图像三维重建方法。所述方法为选取三维重建的CT序列图像，建立立方体素与等值面相交关系索引表，确定等值面与立方体相交的边；设定三维数据场中的等值面的阈值，确定三维图像的立方体素；在相邻两层图片中取一顶点，遍历所有立方体素，得到等值点集合；缩小所述CT序列图像中顶点的距离，确定等值面与立方体素每条边的交点，遍历所有长度为1的立方体素，得到缩小后的等值点集合；根据缩小后的等值点集合，重新构建三角面片，输出等值面图像。本发明可以在保证精度的条件下大大缩短CT序列图像的三维重建时间。

Description

一种基于Marching Cubes的CT图像三维重建方法

技术领域

本发明涉及CT图像三维重建技术领域，是一种基于Marching Cubes的CT图像三维重建方法。

背景技术

近年来，随着现代医学成像技术的发展，先后出现了计算机断层扫描(CT)，核磁共振成像(MRI)，超声(US)等医学成像技术。但这些医疗仪器只能提供人体内部组织或器官的二维图像，这无法满足对数据进行精确分析的要求。在目前的医疗诊断中，医生通过观察CT序列图片，在脑海中建立出三维模型，并凭借自身的经验做出相应的治疗方案。但是这样的方式很明显增加了医生诊断的难度，并且通过这种方式做出的诊断结果准确性很大程度上取决于医生的临床经验，难以达到准确的判断。

考虑到人工诊断的不足，借助计算机进行医学图像的分割以及三维重建工作，可以减少主观判断和临床经验不足对诊断带来的负面影响，极大的提高手术安全性。医学图像的三维重建给医生与患者都带来了很多方便。通过图像三维重建，可构造器官的三维模型，形象逼真地显示人体器官的立体模型，并能实现对重建模型的旋转、缩放等交互操作，从而使医生能够更充分地了解病灶的性质及周围组织的三维结构关系，达到更准确诊疗的目的。

目前，随着医疗设备的发展，医学影像数据集越来越大，海量数据给传统的三维重建技术带来了巨大的挑战，研究大量数据的三维医学影像的快速重建和绘制技术一直是一个非常有挑战性的问题。

三维重建的方法主要分为三维面绘制和三维体绘制两大类。由于体绘制运算量巨大，无法满足临床交互的需求，因此，面绘制仍是临床三维重建的主流，面绘制使由原始的三维医学图像数据先重建出中间几何单位，通过几何单元拼接来描述物体的三维结构。在医学图像三维重建中，经典的等值面重建算法移动立方体(Marching Cubes，MC)算法是一种比较常用的算法。但由于面MC算法中抽取等值面构建的表面几何模型所包含的三角面片数量较大，导致算法速度较慢。

发明内容

本发明为解决经典的等值面重建算法Marching Cubes在提取任意三维数据场的等值面时，需要逐个计算立方体与等值面的相交情况，得到所有的相交三角面，所以巨大的计算量会导致三维等值面显示延时的问题，本发明提供了一种基于Marching Cubes的CT图像三维重建方法，本发明提供了以下技术方案：

一种基于Marching Cubes的CT图像三维重建方法，包括以下步骤：

步骤1：选取三维重建的CT序列图像；

步骤2：建立立方体素与等值面相交关系拓扑信息表，确定等值面与立方体相交的边；

步骤3：设定三维数据场中的等值面的阈值，确定三维图像的立方体素；

步骤4：在相邻两层图片中取一顶点，所述CT序列图像包含8个顶点，确定等值面与立方体素每条边的交点即为等值点，遍历所有立方体素，得到等值点集合；

步骤5：缩小所述CT序列图像中顶点的距离，确定等值面与立方体素每条边的交点，遍历所有长度为1的立方体素，得到缩小后的等值点集合；

步骤6：根据等值点集合和缩小后的等值点集合，重新构建三角面片，输出等值面图像。

优选地，所述步骤2具体为：

步骤2.1：选取立方体素上的任意一个顶点，当所述任意一顶点位于待成像物体内部时，则设所述任意一顶点的状态值为1，表示此顶点位于等值面内或等值面上；当所述任意一顶点位于待成像物体外部时，则设置所述任意一顶点的状态值为0，表示此顶点位于表面外；

步骤2.2：按顺序将立方体素上的顶点的状态值读入一个8位二进制数，得到立方体素的索引值；

步骤2.3：利用旋转对称性与互补对称性，将立方体素所有情况归结为15种基本情况，构建含有15个表项的立方体素与等值面相交关系拓扑信息表，根据所述拓扑信息表确定等值面与立方体相交的边。

优选地，所述步骤4具体为：

步骤4.1：在相邻两层图片中相隔2个点取一个顶点，CT上下两个断层面共8个顶点，8个顶点构成一个长度为3、宽度高度为1的立方体素，假定函数值沿立方体素边界呈线性变化，等值面与立方体素每条边的交点为等值点，采用线性插值法确定等值点的坐标值；

步骤4.2：当等值点所在边两个端点A₁(x₁,y₁,z₁)和A₂(x₂,y₂,z₂)的灰度值分别为F₁(x₁,y₁,z₁)和F₂(x₂,y₂,z₂)时，确定等值点K的坐标(x,y,z)，通过下式表示所述等值点K的坐标(x,y,z)：

其中，x，y和z为等值点K的坐标，c为等值面阈值；

步骤4.3：遍历所有长度为3、宽度高度为1立方体素，得到等值点集合。

优选地，所述步骤5具体为：

步骤5.1：缩小顶点间的距离，CT上下两个断层面共8个顶点，8个顶点组成一个长度为1的立方体素，根据点A₁与点A₃计算得到的等值点K1，根据点A₃与点A₄计算得到的等值点K2，根据点A₄与点A₂计算得到的等值点K3；

步骤5.2：确定等值点K在x轴上的坐标，通过下式表示等值点K在x轴上的坐标：

根据等值点K在x轴上的坐标，确定K1点在x轴上的坐标，通过下式表示K1点在x轴上的坐标：

其中，K1_x为点K1在x轴上的坐标值，F(x₃,y₃,z₃)为坐标为(x₃,y₃,z₃)的点A₃的灰度值，F(x₄,y₄,z₄)为坐标为(x₄,y₄,z₄)的点A₄的灰度值；

确定K3点在x轴上的坐标，通过下式表示K3点在x轴上的坐标：

其中，K3_x为点K3在x轴上的坐标值；

步骤5.3：在求取K2点的坐标时，采用棱边的中点来代替线性插值点，通过下式表示K2在x轴上的坐标：

其中，K2_x为点K2在x轴上的坐标值；

步骤5.4：重复步骤5.1到步骤5.3，计算出长度为1的立方体素的全部等值点的集合。

优选地，所述步骤6具体为：

步骤6.1：根据等值点集合和缩小后的等值点的集合，比较每个等值点与等值面阈值c的大小，确定该顶点是0或者1，把这8个顶点构成的01串组成一个8位的索引值；根据立方体素与等值面相交关系拓扑信息表，确定立方体素内三角面的构型；

步骤6.2：计算三角面片顶点所在边的两个端点的梯度,再通过线性插值得到三角面片顶点的法向量，采用中心差分法确定立方体素顶点的梯度，通过下式表示三角面片顶点法向量在坐标系的分量：

其中，g_x(x,y,z)为三角面片顶点法向量在x轴方向的分量，g_y(x,y,z)为三角面片顶点法向量在y轴方向的分量，g_z(x,y,z)为三角面片顶点法向量在z轴方向的分量，Δx为两个端点在x轴方向上的增量，Δy为两个端点在y轴方向的增量，Δz为两个端点在z轴方向上的增量；

步骤6.3：根据计算出的三角面片顶点的法向量，使用绘图函数绘制各三角面片，选择Gouraud光照模型进行光照计算，生成平滑的等值面图像。

优选地，选取三维重建的CT序列图像，选取格式为.bmp文件，像素大小为512×512，所述CT序列图像的灰度值范围为0到255之间的

本发明具有以下有益效果：

本发明通过设定所建立的等值面的灰度值进行三维重建，然后只需要将CT序列图像作为三维重建模型的输入，即可得到三维重建结果。由于最终等值点集计算过程中利用到之前的交点集计算结果，使得本发明可以在保证精度的条件下大大缩短CT序列图像的三维重建时间。

附图说明

图1是基于Marching Cubes的CT图像三维重建方法流程图；

图2是立方体素的示意图；

图3是同一个方向，以x轴为例上的4个顶点的示意图；

图4是利用CT序列进行三维重建的示意图；

图5为15种立方体素与等值面相交关系情况示意图。

具体实施方式

以下结合具体实施例，对本发明进行了详细说明。

具体实施例一：

根据图1和图2所示，本发明提供一种基于Marching Cubes的CT图像三维重建方法，包括以下步骤：

步骤1：选取三维重建的CT序列图像，选取格式为.bmp文件，像素大小为512×512，所述CT序列图像的灰度值范围为0到255之间；

步骤2：建立立方体素与等值面相交关系拓扑信息表，确定等值面与立方体相交的边；

对于立方体素上的任意一个顶点，如果它位于待成像物体内部，则设该顶点状态值为1，表示此顶点位于等值面内或等值面上；

如果该顶点位于待成像物体外部，则设置其状态值为0，表示此顶点位于表面外；

按同一顺序把个顶点的状态值(0或1)分别读入一个8位二进制数，得到该立方体素的索引值；

因为每个立方体素有8个顶点，每个顶点有两种状态，即等值面内和等值面外，因而立方体素与等值面的关系有2⁸＝256种情况；

利用旋转对称性与互补对称性，将立方体素所有256种情况归结为15种基本情况，具体情况如图5所示，构建含有15个表项的立方体素与等值面相交关系拓扑信息表。

步骤3：设定三维数据场中的等值面的阈值，确定三维图像的立方体素；

步骤4：在相邻两层图片中取一顶点，所述CT序列图像包含8个顶点，确定等值面与立方体素每条边的交点即为等值点，遍历所有立方体素，得到等值点集合；

当立方体素很小时，假定函数值沿立方体素边界呈线性变化，基于这个假设，等值面与立方体素每条边的交点，即等值点的坐标值可以采用线性插值法求出；

假设等值点所在边两个端点A₁(x₁,y₁,z₁)和A₂(x₂,y₂,z₂)的灰度值分别为F₁(x₁,y₁,z₁)和F₂(x₂,y₂,z₂)，那么在该边上的三角面片顶点K的坐标(x,y,z)，计算公式如下：

遍历所有长度为3、宽度高度为1立方体素，得到等值点集合。

步骤5：缩小所述CT序列图像中顶点的距离，确定等值面与立方体素每条边的交点，遍历所有长度为1的立方体素，得到缩小后的等值点集合；

缩小顶点间的距离，CT上下两个断层面共8个顶点，这8个顶点组成一个长度为1的立方体素，图3为同一个方向(以x轴方向为例)上的4个顶点；，根据点A₁与点A₃计算得到的等值点K1，根据点A₃与点A₄计算得到的等值点K2，根据点A₄与点A₂计算得到的等值点K3；

计算K1点的坐标，由于K在x轴上的坐标如下：

根据上述结果，K1在x轴上的坐标计算公式如下：

其中，K1_x为点K1在x轴上的坐标值，F(x₃,y₃,z₃)为坐标为(x₃,y₃,z₃)的点A₃的灰度值，F(x₄,y₄,z₄)为坐标为(x₄,y₄,z₄)的点A₄的灰度值；

同理，K3在x轴上的坐标计算公式如下：

在求取K2点的坐标时，本文采用棱边的中点来代替线性插值点，K2在x轴上的坐标计算公式如下：

重复上述步骤，计算出全部等值点的坐标集合。

步骤6：根据等值点集合和缩小后的等值点集合，重新构建三角面片，输出等值面图像。

步骤6.1：根据等值点集合和缩小后的等值点的集合，比较每个等值点与等值面阈值c的大小，确定该顶点是0或者1，把这8个顶点构成的01串组成一个8位的索引值；根据立方体素与等值面相交关系拓扑信息表，确定立方体素内三角面的构型；

那么在该边上的三角面片顶点坐标(x,y,z)，为了求得三角形顶点的法向量，先计算顶点所在边的两个端点的梯度,再通过线性插值得到三角形顶点法向量。立方体素顶点的梯度用中心差分法计算，公式如下：

其中，g_x(x,y,z)为三角面片顶点法向量在x轴方向的分量，g_y(x,y,z)为三角面片顶点法向量在y轴方向的分量，g_z(x,y,z)为三角面片顶点法向量在z轴方向的分量，Δx为两个端点在x轴方向上的增量，Δy为两个端点在y轴方向的增量，Δz为两个端点在z轴方向上的增量；

根据计算出的三角面片顶点的法向量，使用绘图函数绘制各三角面片，选择Gouraud光照模型进行光照计算，生成平滑的等值面图像，如图4所示利用CT序列进行三维重建的示意图。

以上所述仅是一种基于Marching Cubes的CT图像三维重建方法的优选实施方式，一种基于Marching Cubes的CT图像三维重建方法的保护范围并不仅局限于上述实施例，凡属于该思路下的技术方案均属于本发明的保护范围。应当指出，对于本领域的技术人员来说，在不脱离本发明原理前提下的若干改进和变化，这些改进和变化也应视为本发明的保护范围。

Claims (6)

1.一种基于Marching Cubes的CT图像三维重建方法，其特征是：包括以下步骤：

步骤1：选取三维重建的CT序列图像；

步骤2：建立立方体素与等值面相交关系拓扑信息表，确定等值面与立方体相交的边；

步骤3：设定三维数据场中的等值面的阈值，确定三维图像的立方体素；

步骤4：在相邻两层图片中取一顶点，所述CT序列图像包含8个顶点，确定等值面与立方体素每条边的交点即为等值点，遍历所有立方体素，得到等值点集合；

步骤5：缩小所述CT序列图像中顶点的距离，确定等值面与立方体素每条边的交点，遍历所有长度为1的立方体素，得到缩小后的等值点集合；

步骤6：根据等值点集合和缩小后的等值点集合，重新构建三角面片，输出等值面图像。

2.根据权利要求1所述的一种基于Marching Cubes的CT图像三维重建方法，其特征是：所述步骤2具体为：

步骤2.1：选取立方体素上的任意一个顶点，当所述任意一顶点位于待成像物体内部时，则设所述任意一顶点的状态值为1，表示此顶点位于等值面内或等值面上；当所述任意一顶点位于待成像物体外部时，则设置所述任意一顶点的状态值为0，表示此顶点位于表面外；

步骤2.2：按顺序将立方体素上的顶点的状态值读入一个8位二进制数，得到立方体素的索引值；

步骤2.3：利用旋转对称性与互补对称性，将立方体素所有情况归结为15种基本情况，构建含有15个表项的立方体素与等值面相交关系拓扑信息表，根据所述拓扑信息表确定等值面与立方体相交的边。

3.根据权利要求1所述的一种基于Marching Cubes的CT图像三维重建方法，其特征是：所述步骤4具体为：

步骤4.1：在相邻两层图片中相隔2个点取一个顶点，CT上下两个断层面共8个顶点，8个顶点构成一个长度为3、宽度高度为1的立方体素，假定函数值沿立方体素边界呈线性变化，等值面与立方体素每条边的交点为等值点，采用线性插值法确定等值点的坐标值；

步骤4.2：当等值点所在边两个端点A₁(x₁,y₁,z₁)和A₂(x₂,y₂,z₂)的灰度值分别为F₁(x₁,y₁,z₁)和F₂(x₂,y₂,z₂)时，确定等值点K的坐标(x,y,z)，通过下式表示所述等值点K的坐标(x,y,z)：

其中，x，y和z为等值点K的坐标，c为等值面阈值；

步骤4.3：遍历所有长度为3、宽度高度为1立方体素，得到等值点集合。

4.根据权利要求1所述的一种基于Marching Cubes的CT图像三维重建方法，其特征是：所述步骤5具体为：

步骤5.1：缩小顶点间的距离，CT上下两个断层面共8个顶点，8个顶点组成一个长度为1的立方体素，根据点A₁与点A₃计算得到的等值点K1，根据点A₃与点A₄计算得到的等值点K2，根据点A₄与点A₂计算得到的等值点K3；

步骤5.2：确定等值点K在x轴上的坐标，通过下式表示等值点K在x轴上的坐标：

根据等值点K在x轴上的坐标，确定K1点在x轴上的坐标，通过下式表示K1点在x轴上的坐标：

其中，K1_x为点K1在x轴上的坐标值，F(x₃,y₃,z₃)为坐标为(x₃,y₃,z₃)的点A₃的灰度值，F(x₄,y₄,z₄)为坐标为(x₄,y₄,z₄)的点A₄的灰度值；

确定K3点在x轴上的坐标，通过下式表示K3点在x轴上的坐标：

其中，K3_x为点K3在x轴上的坐标值；

步骤5.3：在求取K2点的坐标时，采用棱边的中点来代替线性插值点，通过下式表示K2在x轴上的坐标：

其中，K2_x为点K2在x轴上的坐标值；

步骤5.4：重复步骤5.1到步骤5.3，计算出长度为1的立方体素的全部等值点的集合。

5.根据权利要求1所述的一种基于Marching Cubes的CT图像三维重建方法，其特征是：

所述步骤6具体为：

步骤6.1：根据等值点集合和缩小后的等值点的集合，比较每个等值点与等值面阈值c的大小，确定该顶点是0或者1，把这8个顶点构成的01串组成一个8位的索引值；根据立方体素与等值面相交关系拓扑信息表，确定立方体素内三角面的构型；

步骤6.2：计算三角面片顶点所在边的两个端点的梯度,再通过线性插值得到三角面片顶点的法向量，采用中心差分法确定立方体素顶点的梯度，通过下式表示三角面片顶点法向量在坐标系的分量：

其中，g_x(x,y,z)为三角面片顶点法向量在x轴方向的分量，g_y(x,y,z)为三角面片顶点法向量在y轴方向的分量，g_z(x,y,z)为三角面片顶点法向量在z轴方向的分量，Δx为两个端点在x轴方向上的增量，Δy为两个端点在y轴方向的增量，Δz为两个端点在z轴方向上的增量；

步骤6.3：根据计算出的三角面片顶点的法向量，使用绘图函数绘制各三角面片，选择Gouraud光照模型进行光照计算，生成平滑的等值面图像。

6.根据权利要求1所述的一种基于Marching Cubes的CT图像三维重建方法，其特征是：选取三维重建的CT序列图像，选取格式为.bmp文件，像素大小为512×512，所述CT序列图像的灰度值范围为0到255之间。

Images