CN107784684A - 一种锥束ct三维重建方法及系统 - Google Patents

Abstract

本发明提供的锥束CT三维重建方法和系统，采用GPU加速的方式，同时利用几何对称性减少计算量，对每一条像素与射线源之间的连线建立一个线程，计算该连线在三维体数据中每一个体素内的长度，进而通过迭代方法得到x最优解，从而达到高效进行三维重建的目的，且使得重建图像质量得到保证的同时大大改善计算效率。

Description

一种锥束CT三维重建方法及系统

技术领域

本发明涉及X射线CT技术领域重建算法技术，尤其是涉及一种锥束CT三维重建方法及系统。

背景技术

相比于传统的CT技术，平板探测器由于有更高的图像分辨率和更大的FOV，因此应用平板探测器技术的锥束CT在成像质量与成像效率上都更有优势。因此如何在锥束CT模型下快速、准确的进行三维重建，成为十分重要的问题。

专利CN104899903A，CN102609978A，CN102279970A提供了基于FDK算法的三维重建。作为滤波反投影累算法，虽然计算速度较快，对系统资源要求低，但重建图像质量差，存在伪影，影响可视化效果和医生诊断。

有文献提出代数迭代算法ART，以解决滤波反投影算法的图像质量问题。ART算法需要计算投影矩阵，目前通常是采用Siddon等提出的算法，该算法是逐点计算方法，即根据射线上某一点和待重建体素间的相对位置关系判断射线上下一点的位置，计算量较大，计算效率较低，尤其是在待重建体数据为三维的情况下缺陷更为明显，限制了ART类算法在临床上的应用。虽然有一些改进算法，但是多数是采用近似假设，如仅判断射线是否穿越体素，穿越记为1，否则记为0的01近似法，会影响投影矩阵的准确度。

以FDK算法为代表的滤波反投影算法，虽然计算速度较快，对系统资源要求低，但重建图像质量差，存在伪影，影响可视化效果和医生诊断。ART算法需要计算投影矩阵，目前通常是采用Siddon等提出的算法，该算法是逐点计算方法，即根据射线上某一点和待重建体素间的相对位置关系判断射线上下一点的位置，计算量较大，计算效率较低，尤其是在待重建体数据为三维的情况下缺陷更为明显，限制了ART类算法在临床上的应用。虽然有一些改进算法，但是多数是采用近似假设，如仅判断射线是否穿越体素，穿越记为1，否则记为0的01近似法，会影响投影矩阵的准确度。

发明内容

有鉴如此，有必要针对现在技术存在的缺陷，提供一种能够快速进行锥束CT三维重建方法。

为实现上述目的，本发明采用下述技术方案：

一种锥束CT三维重建方法，包括下述步骤：

步骤S110：采集投影数据；

步骤S120：构建第一目标函数，F＝|(Ax-b)|²+λR(x)，其中，A为投影矩阵，x为待重建三维体数据，b为实际采集投影图像，b＝[b₀,b₁...b_i...b_N]^T，R为约束项，λ为调整系数，对每一帧图像b_i，计算对应帧图像的射线源所处三维位置坐标；

步骤S130：构建第二目标函数A^T(Ax_n-b)，其中，T_i为每一帧图像i对应的分量，并根据计算x_n+1；

步骤S140：重复步骤S120至S130，直到达到终止条件。

在一些实施例中，步骤S110中，所述投影数据以unsigned short类型数据保存在N*H*W三维数组中，其中N为投影数据帧数，H为投影数据高度，W为投影数据宽度。

在一些实施例中，在完成步骤S120后进行步骤S130前，还包括下述步骤：

步骤S121：计算对每一帧图像b_iX轴左半轴的每一个像素点分配一个线程，计算该的像素点在空间中的坐标，求得连接该坐标与该帧图像对应射线源坐标的直线方程；

步骤S122：判断该直线与待重建的三维体数据中位于X轴左半轴的每一个体素是否有相交；

若有相交，则计算该直线在该体素内部分的长度值，并记录该体素在三维体数据中的三维位置(x,y,z)，并计算该位置关于X轴的轴对称位置，对应的长度值保持不变；

若无相交，记录该体素在三维体数据中的三维位置(x,y,z)，设长度值为0，且该体素关于X轴的轴对称位置的对应长度也设置为0。

步骤S123：当所有线程执行完成后，即得到对应当前帧图像的所有像素点b_ij的投影矩阵A_ij，并计算对应分量M为图像所含像素数。

另外，本发明还提供了一种锥束CT三维重建系统，包括：

投影数据采集单元，用于采集投影数据；

第一计算单元，构建第一目标函数，F＝|(Ax-b)|²+λR(x)，其中，A为投影矩阵，x为待重建三维体数据，b为实际采集投影图像，b＝[b₀,b₁...b_i...b_N]^T，R为约束项，λ为调整系数，对每一帧图像b_i，计算对应帧图像的射线源所处三维位置坐标；

第二计算单元，用于构建第二目标函数A^T(Ax_n-b)，其中，T_i为每一帧图像i对应的分量，并根据计算x_n+1；

迭代单元，重复所述第一计算单元及第二计算单元，直到达到终止条件。

在一些实施例中，还包括位于所述第一计算单元和第二计算单元之间的线程分配单元，包括：

第一计算模块：对每一帧图像b_iX轴左半轴的每一个像素点分配一个线程，计算该像素点在空间中的坐标，求得连接该坐标与该帧图像对应射线源坐标的直线方程；

判断模块：判断该直线与待重建的三维体数据中位于X轴左半轴的每一个体素是否有相交；

若有相交，则计算该直线在该体素内部分的长度值，并记录该体素在三维体数据中的三维位置(x,y,z)，并计算该位置关于X轴的轴对称位置，对应的长度值保持不变；

若无相交，记录该体素在三维体数据中的三维位置(x,y,z)，设长度值为0，且该体素关于X轴的轴对称位置的对应长度也设置为0。

第二计算模块：当所有线程执行完成后，即得到对应当前帧图像的所有像素点b_ij的投影矩阵A_ij，并计算对应分量M为图像所含像素数。

本发明采用上述技术方案的优点是：

本发明提供的锥束CT三维重建方法和系统，采用GPU加速的方式，同时利用几何对称性减少计算量，对每一条像素与射线源之间的连线建立一个线程，计算该连线在三维体数据中每一个体素内的长度，进而通过迭代方法得到x最优解，从而达到高效进行三维重建的目的，且使得重建图像质量得到保证的同时大大改善计算效率。

附图说明

图1为本发明实施例提供的锥束CT三维重建方法的步骤流程图。

图2为本发明实施例设定的三维体数据坐标系和投影数据坐标系。

图3为本发明实施例提供的锥束CT三维重建系统的结构示意图。

图4为本发明实施例提供的锥束CT三维重建系统的线程分配单元的结构示意图。

具体实施方式

请参阅图1，本发明实施例提供的一种锥束CT三维重建方法，包括下步骤：

步骤S110：采集投影数据；

优选地，投影数据以unsigned short类型数据保存在N*H*W三维数组中。其中N为投影数据帧数，H为投影数据高度，W为投影数据宽度。

请参阅图2，设定三维体数据坐标系和投影数据坐标系，为后续方便叙述，假设光源绕X轴进行旋转。

步骤S120：构建第一目标函数，F＝|(Ax-b)|²+λR(x)，其中，A为投影矩阵，x为待重建三维体数据，b为实际采集投影图像，b＝[b₀,b₁...b_i...b_N]^T，R为约束项，λ为调整系数，对每一帧图像b_i，计算对应帧图像的射线源所处三维位置坐标；

可以理解，ART类算法的基本方法如下：令目标函数F＝|(Ax-b)|²+λR(x)。其中A为投影矩阵，x为待重建三维体数据，b为投影图像，R为约束项，λ为权重。迭代寻找x，令目标函数F取最小值。

可以理解，基于GPU并行加速技术，对ART类算法进行并行加速，从而达到高效进行三维重建的目的。

步骤S130：构建第二目标函数A^T(Ax_n-b)，其中，T_i为每一帧图像i对应的分量，并根据计算x_n+1；

优选地，在进行步骤S130后还进行下述步骤：

步骤S121：对每一帧图像b_iX轴左半轴的每一个像素点分配一个线程，计算该像素点在空间中的坐标，求得连接该坐标与该帧图像对应射线源坐标的直线方程；

步骤S122：判断该直线与待重建的三维体数据中位于X轴左半轴的每一个体素是否有相交；

若有相交，则计算该直线在该体素内部分的长度值，并记录该体素在三维体数据中的三维位置(x,y,z)，并计算该位置关于X轴的轴对称位置，对应的长度值保持不变；

若无相交，记录该体素在三维体数据中的三维位置(x,y,z)，设长度值为0，且该体素关于X轴的轴对称位置的对应长度也设置为0。

步骤S123：当所有线程执行完成后，即得到对应当前帧图像的所有像素点b_ij的投影矩阵A_ij，并计算对应分量M为图像所含像素数。

可以理解，本申请采用最速下降法，即其中γ为迭代步长，可以有效降低计算量，提高重建速度。

步骤S140：重复步骤S120至S130，直到达到终止条件。

请参阅图3，为本发明实施例提供的锥束CT三维重建系统的结构示意图，包括：

投影数据采集单元110，用于采集投影数据；

第一计算单元120，构建第一目标函数，F＝|(Ax-b)|²+λR(x)，其中，A为投影矩阵，x为待重建三维体数据，b为实际采集投影图像，b＝[b₀,b₁...b_i...b_N]^T，R为约束项，λ为调整系数，对每一帧图像b_i，计算对应帧图像的射线源所处三维位置坐标；

第二计算单元130，用于构建第二目标函数A^T(Ax_n-b)，其中，T_i为每一帧图像i对应的分量，并根据计算x_n+1；

迭代单元140，重复所述第一计算单元及第二计算单元，直到达到终止条件。

请参阅图4，为本发明实施例提供的锥束CT三维重建系统的线程分配单元150的结构示意图。

优选地，线程分配单元150位于所述第一计算单元120和第二计算单元130之间的，包括：

第一计算模块151：对每一帧图像b_iX轴左半轴的每一个像素点分配一个线程，计算该像素点在空间中的坐标，求得连接该坐标与该帧图像对应射线源坐标的直线方程；

判断模块152：判断该直线与待重建的三维体数据中位于X轴左半轴的每一个体素是否有相交；

若有相交，则计算该直线在该体素内部分的长度值，并记录该体素在三维体数据中的三维位置(x,y,z)，并计算该位置关于X轴的轴对称位置，对应的长度值保持不变；

若无相交，记录该体素在三维体数据中的三维位置(x,y,z)，设长度值为0，且该体素关于X轴的轴对称位置对应的长度也设置为0。

第二计算模块153：当所有线程执行完成后，即得到对应当前帧图像的所有像素点b_ij的投影矩阵A_ij，并计算对应分量M为图像所含像素数。

本发明提供的锥束CT三维重建方法和系统，采用GPU加速的方式，同时利用几何对称性减少计算量，对每一条像素与射线源之间的连线建立一个线程，计算该连线在三维体数据中每一个体素内的长度，进而通过迭代方法得到x最优解，从而达到高效进行三维重建的目的，且使得重建图像质量得到保证的同时大大改善计算效率。

当然本发明的锥束CT投影矩阵的快速计算方法还可具有多种变换及改型，并不局限于上述实施方式的具体结构。总之，本发明的保护范围应包括那些对于本领域普通技术人员来说显而易见的变换或替代以及改型。

Claims (5)

1.一种锥束CT三维重建方法，其特征在于，包括下述步骤：

步骤S110：采集投影数据；

步骤S120：构建第一目标函数，F＝|(Ax-b)|²+λR(x)，其中，A为投影矩阵，x为待重建三维体数据，b为实际采集投影图像，b＝[b₀,b₁...b_i...b_N]^T，R为约束项，λ为调整系数，对每一帧图像b_i，计算对应帧图像的射线源所处三维位置坐标；

步骤S130：构建第二目标函数A^T(Ax_n-b)＝∑T_i，其中，T_i为每一帧图像i对应的分量，并根据计算x_n+1；

步骤S140：重复步骤S120至S130，直到达到终止条件。

2.根据权利要求1所述的锥束CT三维重建方法，其特征在于，步骤S110中，所述投影数据以unsigned short类型数据保存在N*H*W三维数组中，其中N为投影数据帧数，H为投影数据高度，W为投影数据宽度。

3.根据权利要求1所述的锥束CT三维重建方法，其特征在于，在完成步骤S120后进行步骤S130前，还包括下述步骤：

步骤S121：对每一帧图像b_iX轴左半轴的每一个像素点分配一个线程，计算该像素点在空间中的坐标，求得连接该坐标与该帧图像对应射线源坐标的直线方程；

步骤S122：判断该直线与待重建的三维体数据中位于X轴左半轴的每一个体素是否有相交；

若有相交，则计算该直线在该体素内部分的长度值，并记录该体素在三维体数据中的三维位置(x,y,z)，并计算该位置关于X轴的轴对称位置，对应的长度值保持不变；

若无相交，记录该体素在三维体数据中的三维位置(x,y,z)，设长度值为0，且该体素关于X轴的轴对称位置的对应长度也设置为0；

步骤S123：当所有线程执行完成后，即得到对应当前帧图像的所有像素点b_ij的投影矩阵Ai_j，并计算对应分量M为图像所含像素数。

4.一种锥束CT三维重建系统，其特征在于，包括：

投影数据采集单元，用于采集投影数据；

第一计算单元，构建第一目标函数，F＝|(Ax-b)|²+λR(x)，其中，A为投影矩阵，x为待重建三维体数据，b为实际采集投影图像，b＝[b₀,b₁...b_i...b_N]^T，R为约束项，λ为调整系数，对每一帧图像b_i，计算对应帧图像的射线源所处三维位置坐标；

第二计算单元，用于构建第二目标函数A^T(Ax_n-b)，其中，T_i为每一帧图像i对应的分量，并根据计算x_n+1；

迭代单元，重复所述第一计算单元及第二计算单元，直到达到终止条件。

5.根据权利要求4所述的锥束CT三维重建系统，其特征在于，还包括位于所述第一计算单元和第二计算单元之间的线程分配单元，包括：

第一计算模块：对每一帧图像b_iX轴左半轴的每一个像素点分配一个线程，计算该像素点在空间中的坐标，求得连接该坐标与该帧图像对应射线源坐标的直线方程；

判断模块：判断该直线与待重建的三维体数据中位于X轴左半轴的每一个体素是否有相交；

若有相交，则计算该直线在该体素内部分的长度值，并记录该体素在三维体数据中的三维位置(x,y,z)，并计算该位置关于X轴的轴对称位置，对应的长度值保持不变；

若无相交，记录该体素在三维体数据中的三维位置(x,y,z)，设长度值为0，且该体素关于X轴的轴对称位置的对应长度也设置为0；

第二计算模块：当所有线程执行完成后，即得到对应当前帧图像的所有像素点b_ij的投影矩阵A_ij，并计算对应分量M为图像所含像素数。