CT图像中血管的三维模型重建方法及系统与流程

本发明涉及医学成像技术领域，特别涉及一种对ct(电子计算机断层扫描)图片进行血管的三维模型重建方法及系统。

背景技术：

目前的ct技术已经能使人体体内病变区域实现2d(平面图形)可视化，然而从2d的ct图像中很难想象出病变组织的三维结构。而三维cad(计算机辅助设计)模型的构建能够协助医生对病情进行更加准确的判断。

现有技术中对于ct图像的三维重构，常常采取人为标记的方法：依靠有一定专业经验的人员对于ct图像进行人为标记，并用相关软件进行进一步分割处理。这种方法明显的缺陷在于额外的人工标记成本，人工标记时间和人工本身可能会带来的误差。其中后两点可能会造成医疗诊断的推迟和影响医疗诊断的精准，从而耽误病人进一步的诊疗，造成较为严重的后果。

近年来，基于像素对比度差实现图像自动分割的技术迅猛发展，但其适用于分割的对象与其他周围的对象有明显的色彩和对比度差别的图像。对于模糊的ct图像而言，分割对象(比如说血管)可能会很模糊，如此一来基于像素对比度差的方法分割图片效果会非常差，甚至会分割出多个不需要的对象。

可见，利用现有技术的图像自动分割技术对分割对象较模糊的ct图像进行分割，精确度很低，达不到后续进行三维模型重构的质量要求。

技术实现要素：

本发明要解决的技术问题是为了克服现有技术的图像自动分割技术对分割对象较模糊的ct图像进行分割，精确度很低的缺陷，提供一种ct图像中血管的三维模型重建方法及系统。

本发明是通过下述技术方案来解决上述技术问题：

一种ct图像中血管的三维模型重建方法，所述三维模型重建方法包括：

基于神经网络建立图像分割模型；所述图像分割模型用于对ct图像中的血管进行分割，并输出所述血管的二维分割图；

将待建模的ct图像输入至所述图像分割模型，得到多张二维分割图；

对所述二维分割图进行三维重建。

较佳地，基于神经网络建立图像分割模型的步骤，具体包括：

获取已标识的ct图像作为训练样本；

将所述训练样本输入神经网络模型中，根据所述神经网络模型的损失函数对所述神经网络模型进行训练，得到所述图像分割模型。

较佳地，所述损失函数为：

si＝∑i,jmpredict(i,j)×mtrue(i,j)；

sum(mpredict)＝∑i,j|mpredict(i,j)|；

sum(mtrue)＝∑i,j|mtrue(i,j)|；

其中，l表示所述损失函数；mtrue(i,j)表示训练样本的二维分割图中(i,j)位置的像素值；mpredict(i,j)表示所述神经网络模型输出的二维分割图中(i,j)位置的像素值；η表示平滑参数。

较佳地，对所述二维分割图进行三维重建的步骤，具体包括：

将所述二维分割图转化为三维二进制的analyze格式的三维图像；

基于推进立方体算法对所述analyze格式的三维图像进行三维图像重建；

基于nurbs对重建的三维图像进行拟合。

较佳地，对重建的三维图像进行拟合的步骤之前，还包括：

对所述重建的三维图像进行光滑处理。

一种ct图像中血管的三维模型重建系统，所述三维模型重建系统包括：

建模模块，用于基于神经网络建立图像分割模型；所述图像分割模型用于对ct图像中的血管进行分割，并输出所述血管的二维分割图；

三维重建模块，用于将待建模的ct图像输入至所述图像分割模型，得到多张二维分割图，并对所述二维分割图进行三维重建。

较佳地，所述建模模块具体包括：

数据获取单元，用于获取已标识的ct图像作为训练样本；

模型训练单元，用于将所述训练样本输入神经网络模型中，根据所述神经网络模型的损失函数对所述神经网络模型进行训练，得到所述图像分割模型。

较佳地，所述损失函数为：

si＝∑i,jmpredict(i,j)×mtrue(i,j)；

sum(mpredict)＝∑i,j|mpredict(i,j)|；

sum(mtrue)＝∑i,j|mtrue(i,j)|；

其中，l表示所述损失函数；mtrue(i,j)表示训练样本的二维分割图中(i,j)位置的像素值；mpredict(i,j)表示所述神经网络模型输出的二维分割图中(i,j)位置的像素值；η表示平滑参数。

较佳地，所述三维重建模块具体包括：

格式转化单元，用于将所述二维分割图转化为三维二进制的analyze格式的三维图像；

三维重建单元，用于基于推进立方体算法对所述analyze格式的三维图像进行三维图像重建；

拟合单元，用于基于nurbs对重建的三维图像进行拟合。

较佳地，所述三维重建模块还包括：

光滑处理单元，用于对所述重建的三维图像进行光滑处理后，调用所述拟合单元。

本发明的积极进步效果在于：本发明构建的图像分割模型能快速、准确的对ct图像中的血管进行分割，并输出血管的二维分割图，进而精确地构建血管的三维模型。

附图说明

图1为本发明实施例1的ct图像中血管的三维模型重建方法的第一流程图。

图2为本发明实施例1的ct图像中血管的三维模型重建方法的第二流程图。

图3为本发明实施例2的ct图像中血管的三维模型重建系统的模块示意图。

具体实施方式

下面通过实施例的方式进一步说明本发明，但并不因此将本发明限制在所述的实施例范围之中。

实施例1

本实施例提供一种ct图像中血管的三维模型重建方法，如图1所示，三维模型重建方法包括以下步骤：

步骤101、基于神经网络建立图像分割模型。

其中，图像分割模型用于对ct图像进行分割，并输出ct图像中血管的二维分割图。本实施例中，采用全卷积神经网络作为针对ct图像进行自动化血管分割的基本结构。

参见图2，以下对图像分割模型的建立过程进行说明：

步骤101-1、获取已标识的ct图像作为训练样本。

步骤101-1中获取大量的相关的高精度人工标记过的ct图像，一部分作为神经网络模型的训练样本，用于神经网络模型的训练；一部分作为测试样本，用于测试每次迭代训练后的神经网络模型的精准度。

步骤101-2、将训练样本输入神经网络模型中，根据神经网络模型的损失函数对神经网络模型进行训练，得到图像分割模型。

本实施例中，针对血管分割的具体情形，设计专属的卷积神经网络损失函数，损失函数l可以但不限于为：

si＝∑i,jmpredict(i,j)×mtrue(i,j)；

sum(mpredict)＝∑i,j|mpredict(i,j)|；

sum(mtrue)＝∑i,j|mtrue(i,j)|；

其中，mtrue表示训练样本ct图片的真实分割图(二维分割图)数据，真实分割图是经过训练的标记员对原ct图像用相应的标记软件进行标记所获得的；mpredict表示模型训练过程中神经网络模型输出的预测分割图(二维分割图)数据；mtrue(i,j)表示在真实分割图中(i,j)位置的像素值；mpredict(i,j)表示在预测分割图中(i,j)位置的像素值；si表示真实分割图与预测分割图的分割交集(segmentationintersection)；η表示平滑参数。

需要说明的是，η可根据实际需求自行设置，本实施例中，出于数值计算考略，将“平滑参数”设为1。对于神经网络模型，训练模型的架构可根据实际需求(比如说不同的病例ct图片、精确度等)自行选择，相应的模型参数也不同。

以一张训练样本(2d的ct图片，大小为w×l，w表示ct图片的宽度，t表示ct图片的长度)为例，“真实分割图”中分割区域(也即血管所在区域)的像素点为白色，分割区域以外的像素点为黑色。将“真实分割图”表示为矩阵mtrue∈r^w×l，则mtrue(i,j)＝1表示真实分割图中分割区域的像素点，也即白色的像素点；mtrue(i,j)＝0表示真实分割图中分割区域以外的像素点，也即黑色的像素点。同样的，将神经网络模型输出的“预测分割图”表示为mpredict∈r^w×l，则mpredict(i,j)＝1表示预测分割图中分割区域的像素点，也即白色的像素点；mtrue(i,j)＝0表示预测分割图中分割区域以外的像素点，也即黑色的像素点。

步骤101-2的模型训练过程中，参照整个训练和测试过程中观测到的精准度(根据已标识的测试样本计算模型的精准度)，来对卷积神经网络中的结构进行调整，当损失函数收敛到一个局部最优点且在测试数据上观测到满意的表现则停止迭代，从而提高训练的精度。

步骤102、将待建模的ct图像输入至图像分割模型，得到多张血管的二维分割图。

步骤103、对二维分割图进行三维重构。

本实施例中，步骤103具体包括：

步骤103-1、将每张二维分割图转化为三维二进制的analyze格式(一种医学影像格式)的三维图像。

具体的，采用itk(一种影像分析扩展软件工具)将每个二维分割图转化为三维二进制的analyze格式的三维图像。

其中，生成的三维图像包含多个不相互连通的分割出来的区域，利用itk计算出其中体积最大的连通区域并去除分割出来的不与最大分割区域连通的其它像素，也即去除图像中孤立的像素，只保留图像中要进行三维几何重构的部分，即为三维图像。

步骤103-2、基于推进立方体算法(marchingcubes)对analyze格式的三维图像进行三维图像重建。

本实施例中，推进立方体算法主要分成两部分：

第一部分是利用分割和获取(divide-and-conquer)的方法选找图像中分割对象的表面位置。每一个立方体确定对象表面是否如何与其相割并挪动到下一个立方体位置。每一个立方体的顶点根据它们在分割对象内还是外分别赋予1和0。一个立方体有6个顶点，对6个顶点的不同赋值组合代表了对象表面与立方体的不同的相割状态。除去对称性，一共有14种不同的相割模式。对每种相割方式，算法通过对顶点的像素值进行插值可以确定哪一个边与表面相交。

推进立方体算法的第二部分是计算由于对象与立方体形成的三角形的单位法向量。法向量可以用来确定表面的内和外。为了确定三角形的法向量，算法采取对立方体顶点的梯度向量进行插值的方法。每个顶点的梯度向量通过中心差分法来得到：

其中，d(i,j,k)表示三维图像中(i,j,k)位置的像素值。δx、δy和δz分别表示立方体的边长。gx(i,j,k)、gy(i,j,k)和gz(i,j,k)分别表示顶点处x方向、y方向和z方向的梯度值。算法通过在相交边对边的两个顶点进行插值获得相交点的梯度向量，并通过立方体最后在图像的全域推进获得整个分割对象的表面，得到三维网格的结点坐标信息。

步骤103-3、对重建后的三维图像做光滑化处理。

本实施例中，为了使生成的分割对象表面更加光滑或者说为了减少表面大曲率的变化，本实施例中，采用两步连续的高斯光滑化步骤实现光滑化处理算法。

在高斯光滑化处理中，每个三维网格的三角单元的顶点的位置被其邻近的三角形单元的顶点重新计算。一个顶点的坐标向量vm邻近的三角形单元的顶点坐标向量vn，指的是所有和vm共同在一条边上或者同在一个三角形单元上的三角形单元顶点，并且不包含顶点vm。如果将vm临近的三角形单元的顶点定义为sm＝{vn:n＝1,2,3,…,p}，其中p是邻近的顶点的总数，一个平均向量δvm被计算为：

其中，wmo是个权衡参数，其值被设定为1≤o≤p；顶点vm的坐标被更新为：

vm′＝vm+λδvm；

其中，延展参数λ的取值在0到1区间内。

本实施例中，第一步高斯光滑化处理采用一个正的延展参数λ。第二步高斯光滑化处理采用一个负的延展参数μ。并且两个参数满足关系式0<λ<-μ<1。假设第一次高斯光滑化处理得到的顶点坐标为vm′，第二次高斯光滑化处理为：

vm″＝vm′+μδvm′；

其中，δvm′是基于第一次光滑化处理后的新顶点坐标得到的平均向量，vm″是基于第一次光滑化处理后的新顶点坐标得到的第二次光滑化处理的新顶点坐标。通过对光滑处理之后的表面继续采用以上的连续两步高斯光滑化处理，所得到的表面能够足够光滑。

步骤103-4、基于nurbs对经过光滑化处理的三维图像进行拟合。

以下对图像拟合的过程作进一步说明：

(1)点云在三维空间上的主方向分析

为了方便表述，将经过光滑化处理的三维图像的三维点云集合vm″重新定义点云x{x1,x2,…,xr}∈r³，并定义在开始nurbscad三维重构时pca(principalcomponentanalysis)分解得到的前3个主方向为点云坐标系e{e1,e2,e3}。其求解过程简介如下：

构建伴随矩阵其中对伴随矩阵进行特征值分解，即求解∑e＝λe。取e的前3个特征向量为点云初始坐标系e{e1,e2,e3}。

(2)点云在三维空间的刚体转动

由于血管重构的点云拍摄时条件不一，在开始nurbscad三维重构时，空间取向，位置均不确定。为了确保曲面在进入曲面重构阶段时，都具有相同的坐标系，必须先将点云刚体转动至统一的全局坐标系。本实施例采用的方法是svd(singular-valuedecomposition)。设全局坐标系为e{e1,e2,e3}，则相似伴随矩阵构建如下：其中h可认为是一个奇异值矩阵。则其可被分解为h＝u∑v^t。通过分解相似伴随阵h，转动矩阵r可由uv^t求得。刚体位移t则可通过求得。于是，转动点云x{x1,x2,…,xr}∈r³为x＝rx+t。

(3)二维点云的nurbs曲线拟合

本专利采用分层式nurbs曲面的拟合。每个分层拟合为环形nurbs曲线。曲线方程为：

其中，pi为控制点坐标向量，ωi和ωj为其相应的权重值。基函数n由以下递归函数求得：

其中，u＝{u0,…,um}是节点向量，是任意选取的一组非负值、均匀且单调递增的节点向量；p为曲线阶数，p的数值可根据实际需求自行设置，例如本实施例中p＝2。每层的拟合曲线使用一套相同的节点向量。通过构建最小二乘法方程组n^tnp＝n^tx，可以求解得到控制点pi。

(4)张量乘积构建nurbs曲面

由于每层的每条曲线使用的是相同的节点向量，对于相同的的节点，可以沿层数方向建立一套新的节点方向v＝{v0,…,vm}。nurbs曲面则由两个方向的的nurbs曲线通过张量乘积算得。其曲面方程为：

其中，pi,j为在新构建v方向，第j层，第i个控制点的值；n，m分别为u和v方向的曲线阶数。n，m的数值可根据实际需求自行设置，例如本实施例中n＝2，m＝2。

本实施例中，构建的图像分割模型能快速、准确的对ct图像中的血管进行分割，并输出血管的二维分割图，进而实现精确地构建血管的三维模型。

实施例2

如图3所示，本实施例的ct图像中血管的三维模型重建系统包括：建模模块1和三维重建模块2。

建模模块1用于基于神经网络建立图像分割模型。其中，图像分割模型用于对ct图像中的血管进行分割，并输出ct图像中血管的二维分割图。本实施例中，采用全卷积神经网络作为针对ct图像进行自动化血管分割的基本结构。

本实施例中，建模模块1具体包括：数据获取单元11和模型训练单元12。

数据获取单元11用于获取已标识的ct图像作为训练样本。具体的，数据获取单元11获取大量的相关的高精度人工标记过的ct图像，一部分作为神经网络模型的训练样本，用于神经网络模型的训练；一部分作为测试样本，用于测试每次迭代训练后的神经网络模型的精准度。

模型训练单元12用于将训练样本输入神经网络模型中，根据神经网络模型的损失函数对神经网络模型进行训练，得到图像分割模型。

在模型训练的过程中，模型训练单元12参照整个训练和测试过程中观测到的精准度(根据已标识的测试样本计算模型的精准度)，对卷积神经网络中的结构进行调整，当损失函数收敛到一个局部最优点且在测试数据上观测到满意的表现则停止迭代，从而提高训练的精度。

本实施例中，针对血管分割的具体情形，设计专属的卷积神经网络损失函数，损失函数l可以但不限于为：

si＝∑i,jmpredict(i,j)×mtrue(i,j)；

sum(mpredict)＝∑i,j|mpredict(i,j)|；

sum(mtrue)＝∑i,j|mtrue(i,j)|；

其中，mtrue表示训练样本ct图片的真实分割图(二维分割图)数据，真实分割图是经过训练的标记员对原ct图像用相应的标记软件进行标记所获得的；mpredict表示模型训练过程中神经网络模型输出的预测分割图(二维分割图)数据；mtrue(i,j)表示在真实分割图中(i,j)位置的像素值；mpredict(i,j)表示在预测分割图中(i,j)位置的像素值；si表示真实分割图与预测分割图的分割交集(segmentationintersection)；η表示平滑参数。

需要说明的是，η可根据实际需求自行设置，本实施例中，出于数值计算考略，将“平滑参数”设为1。对于神经网络模型，训练模型的架构可根据实际需求(比如说不同的病例ct图片、精确度等)自行选择，相应的模型参数也不同。

以一张训练样本(2d的ct图片，大小为w×l，w表示ct图片的宽度，t表示ct图片的长度)为例，“真实分割图”中分割区域(也即血管所在区域)的像素点为白色，分割区域以外的像素点为黑色。将“真实分割图”表示为矩阵mtrue∈r^w×l，则mtrue(i,j)＝1表示真实分割图中分割区域的像素点，也即白色的像素点；mtrue(i,j)＝0表示真实分割图中分割区域以外的像素点，也即黑色的像素点。同样的，将神经网络模型输出的“预测分割图”表示为mpredict∈r^w×l，则mpredict(i,j)＝1表示预测分割图中分割区域的像素点，也即白色的像素点；mtrue(i,j)＝0表示预测分割图中分割区域以外的像素点，也即黑色的像素点。

三维重建模块2用于将待建模的ct图像输入至图像分割模型，得到多张二维分割图，并对二维分割图进行三维重建。

本实施例中，三维重建模块2具体包括：格式转化单元21、三维重建单元22、光滑处理单元23和拟合单元24。

格式转化单元21用于将二维分割图转化为三维二进制的analyze格式的三维图像，并输出至三维重建单元22。具体的，格式转化单元21采用itk(一种影像分析扩展软件工具)将每个二维分割图转化为三维二进制的analyze格式的三维图像。

其中，生成的三维图像包含多个不相互连通的分割出来的区域，利用itk计算出其中体积最大的连通区域并去除分割出来的不与最大分割区域连通的其它像素，也即去除图像中孤立的像素，只保留图像中要进行三维几何重构的部分，即为三维图像。

三维重建单元22用于基于推进立方体算法对analyze格式的三维图像进行三维图像重建，并输出至光滑处理单元23。

本实施例中，推进立方体算法主要分成两部分：

第一部分是利用分割和获取(divide-and-conquer)的方法选找图像中分割对象的表面位置。每一个立方体确定对象表面是否如何与其相割并挪动到下一个立方体位置。每一个立方体的顶点根据它们在分割对象内还是外分别赋予1和0。一个立方体有6个顶点，对6个顶点的不同赋值组合代表了对象表面与立方体的不同的相割状态。除去对称性，一共有14种不同的相割模式。对每种相割方式，算法通过对顶点的像素值进行插值可以确定哪一个边与表面相交。

推进立方体算法的第二部分是计算由于对象与立方体形成的三角形的单位法向量。法向量可以用来确定表面的内和外。为了确定三角形的法向量，算法采取对立方体顶点的梯度向量进行插值的方法。每个顶点的梯度向量通过中心差分法来得到：

其中，d(i,j,k)表示三维图像中(i,j,k)位置的像素值。δx、δy和δz分别表示立方体的边长。gx(i,j,k)、gy(i,j,k)和gz(i,j,k)分别表示顶点处x方向、y方向和z方向的梯度值。算法通过在相交边对边的两个顶点进行插值获得相交点的梯度向量，并通过立方体最后在图像的全域推进获得整个分割对象的表面，得到三维网格的结点坐标信息。

光滑处理单元23用于对重建的三维图像进行光滑处理后，调用拟合单元24。

本实施例中，为了使生成的分割对象表面更加光滑或者说为了减少表面大曲率的变化，本实施例中，采用两步连续的高斯光滑化步骤实现光滑化处理算法。

在高斯光滑化处理中，每个三维网格的三角单元的顶点的位置被其邻近的三角形单元的顶点重新计算。一个顶点的坐标向量vm邻近的三角形单元的顶点坐标向量vn，指的是所有和vm共同在一条边上或者同在一个三角形单元上的三角形单元顶点，并且不包含顶点vm。如果将vm临近的三角形单元的顶点定义为sm＝{vn:n＝1,2,3,…,p}，其中p是邻近的顶点的总数，一个平均向量δvm被计算为：

δvm＝∑vo∈smwmo(vo-vm)；

其中，wmo是个权衡参数，其值被设定为1≤o≤p；顶点vm的坐标被更新为：

vm′＝vm+λδvm；

其中，延展参数λ的取值在0到1区间内。

本实施例中，第一步高斯光滑化处理采用一个正的延展参数λ。第二步高斯光滑化处理采用一个负的延展参数μ。并且两个参数满足关系式0<λ<-μ<1。假设第一次高斯光滑化处理得到的顶点坐标为vm′，第二次高斯光滑化处理为：

vm″＝vm′+μδvm′；

其中，δvm′是基于第一次光滑化处理后的新顶点坐标得到的平均向量，vm″是基于第一次光滑化处理后的新顶点坐标得到的第二次光滑化处理的新顶点坐标。通过对光滑处理之后的表面继续采用以上的连续两步高斯光滑化处理，所得到的表面能够足够光滑。

拟合单元24用于基于nurbs对重建的三维图像进行拟合。

以下对拟合单元24进行图像拟合的工作原理作进一步说明：

(1)点云在三维空间上的主方向分析

为了方便表述，将经过光滑化处理的三维图像的三维点云集合vm″重新定义点云x{x1,x2,…,xr}∈r³，并定义在开始nurbscad三维重构时pca(principalcomponentanalysis)分解得到的前3个主方向为点云坐标系e{e1,e2,e3}。其求解过程简介如下：

构建伴随矩阵其中对伴随矩阵进行特征值分解，即求解∑e＝λe。取e的前3个特征向量为点云初始坐标系e{e1,e2,e3}。

(2)点云在三维空间的刚体转动

由于血管重构的点云拍摄时条件不一，在开始nurbscad三维重构时，空间取向，位置均不确定。为了确保曲面在进入曲面重构阶段时，都具有相同的坐标系，必须先将点云刚体转动至统一的全局坐标系。本实施例采用的方法是svd(singular-valuedecomposition)。设全局坐标系为e{e1,e2,e3}，则相似伴随矩阵构建如下：其中h可认为是一个奇异值矩阵。则其可被分解为h＝u∑v^t。通过分解相似伴随阵h，转动矩阵r可由uv^t求得。刚体位移t则可通过求得。于是，转动点云x{x1,x2,…,xr}∈r³为x＝rx+t。

(3)二维点云的nurbs曲线拟合

本专利采用分层式nurbs曲面的拟合。每个分层拟合为环形nurbs曲线。曲线方程为：

其中，pi为控制点坐标向量，ωi和ωj为其相应的权重值。基函数n由以下递归函数求得：

其中，u＝{u0,…,um}是节点向量，是任意选取的一组非负值、均匀且单调递增的节点向量；p为曲线阶数，p的数值可根据实际需求自行设置，例如本实施例中p＝2。每层的拟合曲线使用一套相同的节点向量。通过构建最小二乘法方程组n^tnp＝n^tx，可以求解得到控制点pi。

(4)张量乘积构建nurbs曲面

由于每层的每条曲线使用的是相同的节点向量，对于相同的的节点，可以沿层数方向建立一套新的节点方向v＝{v0,…,vm}。nurbs曲面则由两个方向的的nurbs曲线通过张量乘积算得。其曲面方程为：

其中，pi,j为在新构建v方向，第j层，第i个控制点的值；n，m分别为u和v方向的曲线阶数。n，m的数值可根据实际需求自行设置，例如本实施例中n＝2，m＝2。

本实施例中，构建的图像分割模型能快速、准确的对ct图像中的血管进行分割，并输出血管的二维分割图，进而实现精确地构建血管的三维模型。

虽然以上描述了本发明的具体实施方式，但是本领域的技术人员应当理解，这仅是举例说明，本发明的保护范围是由所附权利要求书限定的。本领域的技术人员在不背离本发明的原理和实质的前提下，可以对这些实施方式做出多种变更或修改，但这些变更和修改均落入本发明的保护范围。