CN109509259B

CN109509259B - 一种医学图像三维重建等值面网格提取方法

Info

Publication number: CN109509259B
Application number: CN201811233339.0A
Authority: CN
Inventors: 王沫楠; 段正华
Original assignee: Harbin University of Science and Technology
Current assignee: Harbin University of Science and Technology
Priority date: 2018-10-22
Filing date: 2018-10-22
Publication date: 2020-11-20
Anticipated expiration: 2038-10-22
Also published as: CN109509259A

Abstract

本发明属于医学图像技术领域，提供一种医学体图像三维重建等值面网格提取方法。本发明按顺序依次求解等值面与立方体的各个交点，然后将交点组成三角形网格，并用相邻两个交点一组去确定共用该两点的相邻立方体是否要处理，用队列记录要处理的立方体的信息。整个过程处理的立方体有两种来源：一种为中间一层的所有立方体，另一种为通过队列取出的数据构造的立方体。当队列为空时，生成指定阈值的等值面网格。

Description

一种医学图像三维重建等值面网格提取方法

技术领域

本发明属于医学图像领域，提供一种医学图像三维重建等值面网格提取方法。

背景技术

通过医学影像技术，医生可以得到很多病人体内的特征。但传统方法是医生直接通过CT/MRI二维断层图像去想象三维来分析和诊断，这需要医生有丰富的经验；同时依据传统方法很难准确地把握组织的生理病理的性质、空间几何特征及其周围组织之间的空间关系。通过计算机辅助，从很多二维断层组成的体数据中提取等值面网格，最后将等值面网格以三维模型的形式呈现给医生，可以极大减轻医生的工作量和病症判断的准确性。

网格是三维几何物体表示的基本形式，网格生成技术是医学图像三维重建的最重要的一环。等值面网格提取是通过设定一个标准灰度值，称为阈值，将灰度值等于阈值的点提取出来，由这些点组成的网格。当今，医学图像网格生成最为成熟是MC算法，它是由Lorensen Cline(“Marching cubes A high resolution3D surface constructionalgorithm，”in Computer Graphics(Proc.SIGGRAPH)，vol.21，1987，pp.163-169.)提出的，称为Marching Cubes(简称MC)，之后又有很多方法对该算法进行了改进。

先分别从相邻两层数据中各提取四个点构成立方体构成六面体(图1)，然后对六面体的各个顶点灰度值与阈值进行比较，判断两个点组成的边是否与等值面有交点，最后将这些交点组合成网格。MC算法通过遍历所有的形成的六面体最后生成等值面，其中很大部分时间浪费在了空体素(立方体和等值面没有相交)的计算上，还会产生一些与主体网格不相连的小网格。

发明内容

本方案提供一种基于一种医学图像三维重建等值面网格生成方案，以解决现有方法等值面提取产生大量多余网格的问题，去除一些不必要的计算，消除部分不需要的小网格。

本发明采用的技术方案为以下步骤：

A、通过相邻两层医学断层图像各4个相邻点构造立方体(如图1所示，一个M×N二维数据，相邻两层数据构成(M-1)×(N-1)个立方体；层数为L时，构成(L-1)×(M-1)×(N-1)个立方体，每个立方体会对应一个三维坐标aa、bb、cc，分别对应笛卡尔坐标系的x、y、z三个方向该立方体在所有立方体中的排列序数。立方体的8个顶点灰度值以及空间坐标由立方体的三维坐标计算得到，灰度值通过坐标算出对应二维图像中对应的像素求得，空间坐标先通过坐标算出其中基准点的坐标，再用二维图像像素间对应的层内距离和两个相邻断层的层间距离求得其他顶点的坐标，其中所有的层内、层间距离数据在通用的Dicom格式医学图片中有记录)。从等值面的角度出发，等值面与立方体的面的相交分为两种情况(如图2所示)：一种情况是相交线与相交立方体面上的点组成的最少边多边形为一个三角形，另一种情况是相交线与相交立方体面的点形成一个四边形；以等值面一侧的立方体顶点中，选择一个与相交线组成三角形的顶点为选定点(如图2左图标记为4的点)。

B、遍历立方体的各个顶点，将顶点的灰度值和阈值(等值面的灰度值)进行比较找出所有等值面对应的选定点。先比较立方体顶点灰度值与阈值的大小，将顶点分为大于阈值和不大于阈值的两部分，以数量size(size为点的数量)较小的那部分点为基准(后面称这部分点为满足要求的点，当size等于4时，满足要求的点为大于阈值的点)，size等于0不再对立方体进行后续处理(等值面和立方体不相交)，所有的选定点都在满足要求的点中寻找；将立方体的各个顶点与阈值进行比较，找出立方体中的选定点，具体分为两种情况：

(1)通过下标遍历立方体各个顶点时，找到一个顶点为满足要求的点，和它通过边相连的三个立方体顶点都不是满足要求的点，则该点为选定点(如图4的7号点)，然后size减1，若size等于0时，说明找到所有的选定点；

(2)通过下标遍历各个顶点时，找到一个顶点为满足要求的点(称为当前点)，但和它通过边相连的顶点刚好只有两个不是满足要求的点，具体分为两种情况：一种情况为，size小于4，则当前点为选定点，然后size减2，若size等于0，说明找到所有的选定点，不等于0继续找是否还有(1)所述的选定点；一种情况为，size等于4，再分为两种情况：①满足要求的4个点为两条对角边上的四个点(如图3，两条对角边上的四个点7、3、5和1号点都是满足要求的点)，则当前点和它同平面的满足要求的点都是选定点，此时找到了所有的选定点；②满足要求的4个点不为两条对角边上的四个点，则只有当前点为选定点，然后size减2，若size不等于0找是否还有(1)所述的选定点。

C、求等值面与立方体的交点。等值面与立方体的交点，最后落实到边与等值面相交，一条边若和等值面相交，说明边上的两个顶点的灰度值，一个大于阈值，一个不大于阈值，假设c为阈值，v₁、v₂、c₁、c₂分别为两个顶点的空间坐标和灰度值，则等值面与边的交点坐标为

分两种情况：

一种情况为在立方体中找到了选定点。从选定点出发，依次求出等值面与立方体各边的交点：以选定点出发，在等值面一侧，满足要求的点是相连的，并且满足要求的点最多只有4个(如图5的7、3、0和2号点是相连的)。基于这些特性，先找出选定点通过边相连的两个不是满足要求的点(这两个点是立方体的顶点)，假设与选定点通过边相连的第三个立方体的顶点为基准点，选定点与最开始找到的两个立方体的顶点所在的面命名为SS，面SS相对的面命名为S_S；假设面SS上选定点右侧的立方体顶点为上右点，面SS上选定点左侧的立方体顶点为上左点，面SS上选定点对角的立方体顶点为上对角点，面S_S上基准点右侧的立方体顶点为下右点，面S_S上基准点左侧的立方体顶点为下左点，面S_S上基准点对角的立方体顶点为下对角点，先插值求出选定点和上右点之间的边上、选定点和上左点之间的边上与等值面的交点(如图4，7号点为选定点，分别在7和6号点的边上、7和4号点的边上求等值面的交点)，接下来分情况求剩下的各个交点：(1)、基准点不是满足要求的点：求选定点与基准点之间的边上与等值面的交点(如图4，7和3号点的边上求等值面的交点)；(2)、基准点是满足要求的点，下左点和下右点都是满足要求的点：分别求上左点和下左点之间的边上、下左点和下对角点之间的边上、下右点和下对角点之间的边上、上右点和下右点之间的边上与等值面的交点(如图5，7号点为选定点，分别求0和4之间、0和1之间、1和2之间、2和6之间与等值面的交点)；(3)、基准点是满足要求的点，下左点与下右点中只有下左点是满足要求的点：先求出上左点和下左点的边上与等值面的交点，再分两种情况：一种情况为，下对角点是满足要求的点，分别求上对角点和下对角点之间的边上、下右点和下对角点之间的边上、下右点和基准点之间的边上与等值面的交点；另一种情况为，下对角点不是满足要求的点，分别求下左点和下对角点之间的边上、下右点和基准点之间的边上与等值面的交点(如图6中，7号点为选定点，a图对应第一种情况，求0和4之间、1和5之间、1和2之间、3和2之间与等值面之间的交点；b图对应另一种情况，求0和4、0和1之间、2和3之间与等值面的交点)；(4)、基准点是满足要求的点，下左点与下右点中只有下右点是满足要求的点，先求出下左点和基准点之间的边上与等值面的交点，再分两种情况：一种情况为，下对角点是满足要求的点，分别求下左点和下对角点之间的边上、上对角点和下对角点之间的边上、上右点和下右点之间的边上与等值面的交点；另一种情况为，下对角点不是满足要求的点，分别求下右点和下对角点之间的边上、上右点和下右点之间的边上与等值面的交点(如图7，7号点为选定点，a图对应第一种情况，求0和3之间、0和1之间、1和5之间、2和6之间与等值面之间的交点；b图对应另一种情况，求0和3之间、1和2之间、2和6之间与等值面之间的交点)；(5)、基准点满足要求时，下左点和下右点都不是满足要求的点，分别求下左点和基准点之间的边上、下右点和基准点之间的边上与等值面的交点(如图8中，7号点为选定点，求3和0之间、3和2之间与等值面的交点)。

另一种情况为在立方体中没有找不到选定。，它是因为等值面从中间将立方体一分为二了，交线没有和立方体的点形成三角形，并且此时size等于4(如图9)。此种情况只需要找到立方体那个面的点灰度值都大于阈值，那个面的点灰度都不大于阈值(即找到等值面在那个方向平分立方体)，然后对两个面上相连的点之间依次(顺时针或逆时针)求与等值面的交点(如图9，依次求0和4之间、1和5之间、2和3之间、3和7之间与等值面的交点)。

D、用一个标记数组记录所有立方体的状态，防止同一个立方体被反复计算多次，立方体通过自己的坐标访问数组中相对应的数据：数据为0表示当前立方体没有处理，为1表示当前立方体的信息已经放入了增长队列，为2表示当前立方体已经处理；

E、将所有交点按交点求解顺序三个点一组形成三角形，把形成的三角形放入三角链表中。当所有的交点都找到后，按照顺序3个点一起组成个三角形，然后将三角形放入三角链表中进行渲染。将所有交点按顺序三个点一起组成三角形时，把交点按求解的顺序从0依次递增设定相应的标号。按交点数量分类：交点数为3时，直接三个点形成三角形；交点数为4时，标号为(0，1，2)和(2，3，0)的交点分别组成两个三角形；交点数为5时，标号为(0，1，2)、(2，3，4)和(4，0，2)的交点分别组成三个三角形；交点数为6时，标号为(0，1，2)、(2，3，4)、(4，5，0)、(0，2，4)的交点分别组成四个三角形(如图10所示)。

F、将交点相邻两个为一组，判断共用该个交点相邻立方体是否被处理：首先，先找到两个相邻点所在的平面，令垂直该平面的坐标方向为方向i(相邻两点坐标在i方向的数值相等，点和立方体的坐标都有三个方向，对应笛卡尔坐标系的x，y，z)，然后和相邻第三点比较各自坐标在i方向上的数值大小，若两相邻点在i方向的坐标数值小于第三点在i方向的坐标数值，则所有相邻立方体的坐标在i方向上的数值为本立方体坐标在i方向数字减一，否则所求相邻立方体的坐标在i方向上的数值为本立方体坐标在i方向数字加一，相邻立方体坐标在其他两个方向数值等于本立方体在这两个方向的数值(如图11，以相邻点b和c为例，b和c的坐标在x方向相等，但b和c的x坐标小于a点x的坐标，说明相邻立方体坐标x方向等于本立方体x坐标减一，其他两个方向相等，若本立方体坐标为(aa,bb,cc)，则相邻立方体坐标为(aa-1,bb,cc))；接下来根据所求坐标在标记数组中找出相邻立方体对应的数值，如果对应标记数组中的值为0，先将相应数值置为1，再将所求相邻立方体的坐标生成队列成员的对象，将对象压入队列中；对应标记数组中的值不为0，则不再进行处理。

G、在立方体构造过程中，按构造数据来源，立方体可分为两种：一种为种子立方体，另一种为通过队列的数据构造的立方体。医学体数据是有很多二维图片组成的，本发明取中间两层图片数据组成的立方体为种子立方体，对中间一层所有的立方体依次遍历是否与等值面相交，相交生成网格，并把三角形放入三角链表中，并计算要求增长的另一个立方体的数据，把数据生成对象放入队列中；另一种为从队列中取出的数据，由数据构造的立方体，要先判断立方体对应标记数组的数值是否为2，只有数值不为2才进行对该立方体进行处理，整个过程一直持续到队列中数据为空为止。当一个立方体处理完后，将其状态标记为已处理。

附图说明

图1为从二维图片构成立方体的示意图；

图2为等值面与立方体的面相交的情况；

图3为选定点选择过程中对应的一种特殊情况，可能存在拓扑二义性，需要用从队列传来的函数指针消除二义性；

图4～图8为以选定点出发，等值面与立方体边的相交情况；

图9为立方体中找不到选定点时，等值面与立方体边的相交情况；

图10为将相交点组合三角形；

图11为将相交点两个相邻点一组，判断判断共享该两个交点的另一个立方体是否被处理，没有处理将其数据放入队列，等待处理；

图12为具体程序实现时，以选定点在立方体的坐标出发，其它顶点在立方体中对应的坐标。

具体实施方式

为使本发明的特点更加清晰，结合具体实例，以VTK为平台进行详细说明。

将医学断层dicom格式数据导入程序后，先经过去噪等图像预处理过程，再将图像数据转发为数据流，便于根据下标构建立方体，为本发明算法使用。在VTK平台中，仿照MC算法的实现过程，在它的RequsetData函数(RequsetData为VTK管线中的函数)中有dims[3](表示数据在图像三个空间尺度的大，对应像素的个数)、spacing[3](表示空像在三个空间尺度实际空间大小，对应像素间物体的实际尺寸)、origin[3]和extent[6](origin[3]和extent[6]为基准点相对坐标系原点坐标的偏差系数)用来求解立方体顶点的灰度值和空间坐标。先通过立方体顶点在数据中的坐标找到所有顶点的灰度值，将所有顶点的灰度值和阈值比较，根据发明内容中的判断依据得知等值面与立方体相交时，在交点求解之前，需要根据中心差分计算立方体各个点的梯度，以便后续网格的渲染，在之后的交点求解过程中，也需要插值计算出交点出的梯度。此外，需要根据RequestData提供的信息算出立方体中(0，0，0)点的坐标，其坐标在三个方向数值分别为：pts[0][0][0][0]＝origin[0]+(x+extent[0])*spacing[0]、pts[0][0][0][1]＝origin[1]+(y+extent[2])*spacing[1]、pts[0][0][0][2]＝origin[2]+(z+extent[4])*spacing[2]，其中x,y,z代表立方体在所有立方体中的序数下标，其他顶点的坐标根据空间位置加上像素间的距离spacing。之后就根据灰度值大小求解交点的空间坐标和梯度信息，并将将交点组合成三角形放入三角形链表；在交点求解过程中的具体代码实现时，选定点和不满足要求的两个相连点所组成的平面SS会是x、y和z三个方向的任意平面，需要根据三个方向对应三个方向的代码(按在图12立方体中顶点下标的命名规则，需要对应i、j和k三个方向的代码，右左下方向按i、j和k的循序排列就行：如在图12中，选定点下标为(i,j,k)，不满足要求的两个相邻点坐标若为(i，(j+1)％2，k)和((i+1)％2,j,k)，选定点的右点下标为((i+1)％2,j,k)，选定点左点的下标为(i，(j+1)％2，k)，第三个坐标值为(k+1)％2的为基准点(发明内容交点求解过程中所述)，其它点的下标通过(i+1)％2、(j+1)％2、(k+1)％2求得)，其它过程如发明内容所示。在等值面生成后，就可以用来等值面网格显示三维重建的模型了。

Claims

1.一种医学图像三维重建等值面网格提取方法，其特征在于，包括以下步骤：

A、通过相邻两层医学断层图像各4个相邻点构造立方体，等值面与立方体面相交有两种情况：一种情况是相交线与相交立方体面上的顶点组成的最少边多边形为一个三角形，另一种情况是相交线与相交立方体面的点形成一个四边形；以等值面一侧的立方体顶点中，选择一个与相交线组成三角形的顶点为选定点；

B、遍历立方体的各个顶点，将顶点的灰度值和阈值进行比较找出所有等值面对应的选定点，阈值为等值面的灰度值，选定点寻找过程为：先比较立方体顶点灰度值与阈值的大小，将顶点分为大于阈值和不大于阈值的两部分，以数量size较小的那部分点为基准，size为点的数量，并称为基准的点为满足要求的点，当size等于4时，满足要求的点为大于阈值的点，size等于0时，等值面和立方体不相交，不对立方体进行后续处理，所有的选定点都在满足要求的点中寻找；将立方体的各个顶点与阈值进行比较，找出立方体中的选定点，具体分为两种情况：

(1)通过下标遍历立方体各个顶点时，找到一个顶点为满足要求的点，和它通过边相连的三个立方体顶点都不是满足要求的点，则该点为选定点；然后size减1，若size等于0时，说明找到所有的选定点；

(2)通过下标遍历各个顶点时，找到一个顶点为满足要求的点，称为当前点，但和它通过边相连的顶点刚好只有两个不是满足要求的点，具体分为两种情况：一种情况为，size小于4，则当前点为选定点，然后size减2，若size等于0，说明找到所有的选定点，不等于0继续找是否还有(1)所述的选定点；另一种情况为，size等于4，再分为两种情况：①满足要求的4个点为两条对角边上的四个点，则当前点和它同平面的满足要求的点都是选定点，此时找到了所有的选定点；②满足要求的4个点不为两条对角边上的四个点，则只有当前点为选定点，然后size减2，若size不等于0找是否还有(1)所述的选定点；

C、从选定点出发，依次求出等值面与立方体各边的交点；当在步骤B的过程中找不到选定点时，需要找出立方体大于阈值和不大于阈值各4个顶点组成的两个对面，然后对两个对面之间的相连边依次顺时针或逆时针求与等值面的交点；等值面与立方体边的交点求解的过程为：找出选定点通过边相连的两个不是满足要求的点，假设与选定点通过边相连的第三个立方体的顶点为基准点，选定点和其通过边相连且不满足要求的两个立方体的顶点所在的面命名为SS，面SS相对的面命名为S_S；假设面SS上选定点右侧的立方体顶点为上右点，面SS上选定点左侧的立方体顶点为上左点，面SS上选定点对角的立方体顶点为上对角点，面S_S上基准点右侧的立方体顶点为下右点，面S_S上基准点左侧的立方体顶点为下左点，面S_S上基准点对角的立方体顶点为下对角点，先插值求出选定点和上右点之间的边上、选定点和上左点之间的边上与等值面的交点，接下来分情况求剩下的各个交点：

(1)、基准点不是满足要求的点：求选定点与基准点之间的边上与等值面的交点；

(2)、基准点是满足要求的点，下左点和下右点都是满足要求的点：分别求上左点和下左点之间的边上、下左点和下对角点之间的边上、下右点和下对角点之间的边上、上右点和下右点之间的边上与等值面的交点；

(3)、基准点是满足要求的点，下左点与下右点中只有下左点是满足要求的点：先求出上左点和下左点的边上与等值面的交点，再分两种情况：一种情况为，下对角点是满足要求的点，分别求上对角点和下对角点之间的边上、下右点和下对角点之间的边上、下右点和基准点之间的边上与等值面的交点；另一种情况为，下对角点不是满足要求的点，分别求下左点和下对角点之间的边上、下右点和基准点之间的边上与等值面的交点；

(4)、基准点是满足要求的点，下左点与下右点中只有下右点是满足要求的点，先求出下左点和基准点之间的边上与等值面的交点，再分两种情况：一种情况为，下对角点是满足要求的点，分别求下左点和下对角点之间的边上、上对角点和下对角点之间的边上、上右点和下右点之间的边上与等值面的交点；另一种情况为，下对角点不是满足要求的点，分别求下右点和下对角点之间的边上、上右点和下右点之间的边上与等值面的交点；

(5)、基准点满足要求时，下左点和下右点都不是满足要求的点，分别求下左点和基准点之间的边上、下右点和基准点之间的边上与等值面的交点；

D、用一个标记数组记录所有立方体的状态，防止同一个立方体被反复计算多次，立方体通过自己的坐标访问数组中相对应的数据：数据为0表示当前立方体没有处理，为1表示当前立方体的信息已经放入了增长队列，为2表示当前立方体已经处理；用一个增长队列存储将要处理立方体的坐标；

E、将所有交点按交点求解顺序三个点一组形成三角形，把形成的三角形放入三角链表中；

F、将所有交点按顺序两个相邻点一组，找出共用该两点相邻立方体的状态，步骤D所述的状态，若相邻立方体没有处理，需要将相邻立方体的坐标和本立方体的选定点灰度值放入增长队列中，并将相邻立方体的状态标记为1；

G、立方体处理顺序为：先从中间一层遍历所有立方体进行处理，然后通过增长队列取出的数据构造其他层的立方体进行处理，直到队列为空；当一个立方体处理完后，将其状态标记为已处理。

2.根据权利要求1所述的一种医学图像三维重建等值面网格提取方法，其特征在于，在所述步骤E中，将所有交点按顺序三个点一起组成三角形时，把交点按求解的顺序从0依次递增设定相应的标号，按交点数量分类：交点数为3时，直接三个点形成三角形；交点数为4时，标号为(0，1，2)和(2，3，0)的交点分别组成两个三角形；交点数为5时，标号为(0，1，2)、(2，3，4)和(4，0，2)的交点分别组成三个三角形；交点数为6时，标号为(0，1，2)、(2，3，4)、(4，5，0)、(0，2，4)的交点分别组成四个三角形。

3.根据权利要求1所述的一种医学图像三维重建等值面网格提取方法，其特征在于，在所述步骤F中，共用相邻两点的相邻立方体的状态查找过程为：首先，先找到两个相邻点所在的平面，令垂直该平面的坐标方向为方向i，相邻两点坐标在i方向的数值相等，点和立方体的坐标都有三个方向，对应笛卡尔坐标系的x，y，z，然后和相邻第三点比较各自坐标在i方向上的数值大小，若两相邻点在i方向的坐标数值小于第三点在i方向的坐标数值，则所有相邻立方体的坐标在i方向上的数值为本立方体坐标在i方向数字减一，否则所求相邻立方体的坐标在i方向上的数值为本立方体坐标在i方向数字加一；所求相邻立方体坐标在其他两个方向数值等于本立方体在这两个方向的数值，最后根据所求坐标代入标记数组中判断相邻立方体的状态。