CN103617648A - 一种锥形束ct重建方法和系统 - Google Patents

Abstract

本发明公开了一种锥形束CT重建系统，包括：计算机图形计算卡和计算机主机；计算机主机中包括预先计算权重系数模块；计算机图形计算卡包括显卡显存、并行运算模块；预先计算权重系数模块用于预先计算出锥形束FDK重建所需的权重系数，这些权重系数加载至计算机主机内存中，并将其传输至计算机图形计算卡GPU的显卡显存中；三维体数据重建过程中，计算机图形计算卡先将该投影图像对应的系数传输至GPU的显卡显存中，再用并行运算模块对该投影图像中的像素进行加权、滤波，最后进行反投影操作。此外，本发明还公开了一种锥形束CT重建方法。本发明能显著提高锥形束CT图像重建速度，且可整合至CT扫描过程，实现边采集图像，边重建图像。

Description

一种锥形束CT重建方法和系统

技术领域

本发明属于CT成像技术领域，具体涉及一种基于计算机图形计算卡（GPU）的CT图像重建加速方法，尤其涉及一种锥形束CT重建方法和系统。

背景技术

CT成像技术在现代医疗诊断领域中发挥着重要的作用。CT断层图像重建是CT成像系统的重要组成部分之一，其质量的优略和速度的快慢是衡量CT系统性能的重要指标。1971年第一台CT原型机进行一次CT扫描后大约需要2个半小时进行图像重建过程，然后可以获得一幅CT断层图像。至今，经过40年的发展，由于重建算法的改进和计算硬件的提升，大型螺旋CT的图像重建速度已经达到了每秒40张图像。这些CT成像系统通常使用专用集成电路（ASIC）或者可编程逻辑阵列（FPGA）电路加速图像重建过程。由于该类硬件是针对具体算法进行逻辑功能按需配置，基于该类硬件的CT系统的图像重建速度可以达到每秒40张图像。但其硬件针对具体问题设定，不利于软件升级等缺点限制了其在医疗CT系统中的广泛应用。

目前广泛使用的锥形束CT重建算法是FDK算法，以平板探测器的锥形束CT为例，其主要包括三个步骤：

1.投影图像加权：

$R_{\mathrm{\β}}^{\′}(x,y)=\frac{D_{\mathrm{SO}}}{\sqrt{D_{\mathrm{SO}}^2+x^2+y^2}}{R}_{\mathrm{\β}}(x,y),$

其中β描述该投影图像R对应的CT成像系统的旋转角度，(x，y)是位于旋转中心的虚拟探测器平面的像素位置，D_SO是射线源到虚拟探测器平面的距离。

2.投影图像滤波:

$Q_{\mathrm{\β}}(x,y)=R_{\mathrm{\β}}^{\′}(x,y)*h\left(x\right),$

其中*表示卷积操作，h(x)是滤波核函数。

3.将滤波后的投影图像加权反投影至三维体数据

$f(i,j,k)={\∫}_0^{2\mathrm{\π}}\frac{D_{\mathrm{SO}}^2}{{(D_{\mathrm{SO}}-s)}^2}\·Q_{\mathrm{\β}}(\frac{D_{\mathrm{SO}}x}{D_{\mathrm{SO}}-s},\frac{D_{\mathrm{SO}}y}{D_{\mathrm{SO}}-s})\mathrm{d\β},$

其中s是待重建像素相对于虚拟探测器平面的距离坐标。

由于对每个投影进行步骤1-3是相互独立的，所以可以很直接地对CT系统采集的多幅投影图像进行并行化处理。

通用型图形处理器（GPGPU）具有单指令多数据流（SIMD）运算特点，可以很好的完成FDK重建算法的并行运算。目前基于图像计算卡（GPU）的FDK的图像重建软硬件模块已经成功应用于现有的锥形束CT成像系统。

经相关专利检索发现，现在已经有多个专利描述使用GPU技术加速CT图像重建，如申请号为200810113846.0的中国专利、申请号为201010208314.2的中国专利、申请号为201210010806.X的中国专利和申请号为201310277674.1的中国专利。这些专利主要都是利用GPU的并行计算特点来加速CT图像的重建。除此，目前市场上已经有一些厂商使用GPU加速技术来提高锥形束CT的图像重建速度，如NewTom VGi、Carestream CS90003D和Planmeca Promax3D。根据他们产品的说明书，他们的CT图像重建时间大概在30秒到1分钟。现有的相关的公开专利文献均是使用GPU的SIMD运算并行化CT图像重建过程以达到加速图像重建效果。但GPU的每个计算单元只拥有很简单的计算逻辑硬件，且GPU的核心计算频率不足CPU的一半，使用GPU进行繁杂的计算不利于进一步加速CT图像的重建速度。

发明内容

本发明要解决的技术问题在于针对现有技术中存在的不足提出一种更快的基于GPU的锥形束CT重建方法及系统。该方法相对已有方法可以显著提高锥形束CT图像重建速度，而且可以整合至CT扫描过程，实现一边采集图像，一边重建图像，不需要严格限制单幅图像的处理时间要小于数据采集时间间隔。即使在CT重建时间大于CT扫描时间的情况下也能正确地完成图像重建过程。因为在图像采集过程中已经开始进行CT断层图像重建过程，这样大大减少了用户的等待时间，达到完成CT扫描，也就基本上完成CT图像重建的效果，实现实时CT图像重建。为此，本发明还提供一种锥形束CT重建系统。

本发明是通过以下技术方案实现的：本发明首先根据背景技术中介绍的FDK重建算法预先计算出锥形束FDK重建所需的权重系数，

h(x)和然后在进行锥形束FDK重建过程前，将这些权重系数加载至主机内存中，根据所占空间大小将和h(x)传输至GPU显存中。为GPU并行化计算创建操作流，将FDK重建步骤及数据传输步骤依次记录在操作流中，待GPU顺序执行。在对每个投影角度图像的处理过程中，GPU依次把每幅投影图像及其对应的系数

传输至GPU显存中，并使用已预先传输至GPU中的

和h(x)对投影图像进行加权、滤波操作。最后将滤波后图像映射至纹理显存（texture）中，使用该幅图像对应的系数

进行反投影操作。

本发明包括以下步骤：

第一步，预先计算出锥形束FDK重建所需的权重系数，

h(x)和

其中,h(x)是滤波核函数，(x，y)是位于旋转中心的虚拟探测器平面的像素位置，D_SO是射线源到虚拟探测器平面的距离，s是待重建像素相对于虚拟探测器平面的距离坐标，s与投影角度β相关。

所述的预先计算，是指在使用GPU进行CT并行重建前将相关的权重系数采用CPU预先计算策略离线计算好。

第二步，在进行锥形束FDK重建过程前，将第一步预先计算好的这些权重系数加载至主机内存中，根据所占空间大小将

和h(x)传输至GPU显存中。

所述的在进行锥形束FDK重建过程前，将这些权重系数加载至主机内存中，是指将上一步预先计算好的权重系数加载至主机内存中，获得其内存访问地址入口。

所述的根据所占空间大小将

和h(x)传输至GPU显存中，是指在GPU显存中申请同样大小的空间，将空间占用小的

和h(x)系数直接存放至GPU显存中，供GPU重建过程调用。

第三步，为GPU并行化计算创建操作流。

所述的为GPU并行化计算创建操作流，是指在GPU重建前声明一个流（stream），在GPU进行并行化重建过程中，将第四步中所述的数据传输、图像加权滤波和第五步的反投影等操作加载至这个流中。

第四步，在三维体数据重建过程中，GPU依次对每幅投影图像进行处理：先将该投影图像对应的系数

传输至GPU显存中，然后使用已预先传输至GPU中的和h(x)对投影图像进行加权、滤波操作，最后进行反投影操作。

所述的在三维体数据重建过程中，GPU依次把每幅投影图像及其对应的系数传输至GPU显存中，是指把投影图像写入GPU的CUDA数组，并且只将与该投影图像相关的系数

传输至GPU显存，并步骤记录在第三步创建的操作流中。GPU并行化的基本计算单元是投影图像中的像素。

所述的使用已预先传输至GPU中的

和h(x)对投影图像进行加权、滤波操作，是指使用GPU的多个线程并行地对该投影内的多个图像像素进行加权、滤波，其中使用已经预先计算好，并在第二步中传输至GPU显存中的系数。该操作记录在第三步创建的操作流中。

所述反投影操作中，可以将滤波后图像映射至纹理显存（texture）中，使用该幅图像对应的系数

进行反投影操作。所述的将滤波后图像映射至纹理显存（texture）中，使用该幅图像对应的系数进行反投影操作，是指将滤波后的投影图像与GPU的纹理坐标系进行绑定，使用纹理拾取函数访问纹理存储器，加速锥形束CT重建的反投影过程。该操作记录在第三步创建的流中。

此外，本发明还提供一种锥形束CT重建系统，该系统包括：计算机图形计算卡和计算机主机；所述计算机主机中包括预先计算权重系数模块；所述计算机图形计算卡包括显卡显存、并行运算模块；

所述预先计算权重系数模块用于预先计算出锥形束FDK重建所需的权重系数，这些权重系数加载至所述计算机主机内存中，并将

和h(x)传输至计算机图形计算卡GPU的显卡显存中；在三维体数据重建过程中，计算机图形计算卡先将该投影图像对应的系数

传输至GPU的显卡显存中，然后采用并行运算模块对该投影图像中的像素进行加权、滤波，最后进行反投影操作。

作为优选的技术方案，所述计算机图形计算卡中还可以包括：纹理存储器，在所述反投影操作中使用纹理拾取函数访问纹理存储器，以加速反投影过程。

由于采用了上述技术方案，本发明的有益效果包括：

1.采用CPU预先计算策略，将计算好的参数通过数据流的方式传递给GPU直接使用，减少了GPU的计算操作，进一步加速了GPU的CT图像重建过程。

2.将GPU并行化进行锥形束CT重建的数据传输（主机内存到GPU显存）和GPU内数据操作的过程记录在流中，然后再进行执行是因为在CT成像过程中，CT投影数据的采集速度可能快于使用GPU对每张投影图像的处理速度，那么将投影数据从主机内存到GPU显存的传输和后续的针对该幅图像的一系列操作记录在流中可以等到GPU执行完前一幅图像的操作后再执行当前图像的操作。这样，就可以将CT数据采集与CT图像重建整合在一起，不论CT重建的时间是否大于数据采集时间都可以缩短了总的CT成像时间，否则需要等待完成所有数据采集过程再进行图像重建。

附图说明

图1是本发明锥形束CT重建方法的流程示意图；

图2是本发明CT数据采集过程和CT图像重建过程的示意图；

图3是本发明的锥形束CT重建中GPU操作流中的执行序列示意图；

图4是本发明实施例的CT重建结果示意图，其中，图4A是本发明重建结果的CT断层图；图4B是将一系列断层图像组合成的三维图。

图5是本发明锥形束CT重建系统的结构框图。

具体实施方法

下图结合附图对本发明的实施例做进一步详细说明。

实施例：使用锥形束CT对羊头进行扫描，使用本发明所述方法将CT数据采集与CT图像重建过程整合在一起，实现实时锥形束CT重建。本实施例数据采集系统使用Varian的2050D平板探测器；CT重建使用显卡是Nvidia TITAN。

请参阅图1，依据本发明方法的实施过程包括：

1.根据该成像系统的参数（D_SO，探测器分辨率等）预先计算出锥形束FDK重建所需的权重系数（图1中步骤101），并存储至主机内存中（图1中步骤106）。将内存占用少的系数，

和h(x)保存至主机内存1（图1中步骤102）；将内存占用多的系数，如保存至主机内存2（图1中步骤103）。

2.在显卡显存（图1中步骤107）中申请空间，显卡显存0用于存储单张投影图像数据，显卡显存1用于存储

和h(x)参数，显卡显存2用于存储一个角度下的

参数，显卡显存3用于存储三维体数据（三维数据体是指以三维空间坐标（x,y,z）为位置函数的数据集合体）。

3.在CT成像开始前，将占用内存较少的主机内存1存储至显卡显存1中（图1中步骤109）。

4.创建GPU执行流（stream）。

5.在CT数据采集过程中（图2中步骤201），将CT投影图像序列（图1中步骤104和步骤105）依次保存在主机内存0（图2中步骤202，图1中的步骤105），每张投影图像拥有一个内存地址。

6.在CT数据采集过程中同时进行CT断层重建。每采集一幅投影图像，就将该图像传输至GPU显存中的一个固定的CUDA数组（图1中步骤108，图2中步骤203）的操作加入步骤4创建的流中，并且将存放在主机中的该幅图像对应的

系数传输至显卡显存2中的步骤（图1中步骤110）也记录在步骤4创建的流中。

7.在步骤4创建的流中记录对该张投影图像进行的加权、滤波、加权反投影等操作（图1中步骤111、112）。反投影步骤中使用纹理拾取函数访问纹理存储器，以加速反投影过程。

8.在数据采集过程中，对于投影数据的一系列操作都保存在步骤4创建的流中（图2中步骤204），流中的这一系列操作是依照顺序执行的，与主机CPU的运算无关（图3中的步骤301）。在步骤4创建的流中依次进行着针对每张投影图像的数据传输（图3中步骤303）和并行计算过程（图3中步骤302），直到完成流中的所有操作（图3中步骤304）。

9.几乎在完成CT数据采集的同时完成CT图像重建过程，重建出CT断层图像（见图4A），并可以将一系列断层图像组合成三维图像显示（见图4B）。

如图5所示，本发明一种锥形束CT重建系统，该系统包括：计算机图形计算卡和计算机主机；所述计算机主机中包括预先计算权重系数模块；所述计算机图形计算卡包括显卡显存、并行运算模块；

所述预先计算权重系数模块用于预先计算出锥形束FDK重建所需的权重系数，这些权重系数加载至所述计算机主机内存中，并将和h(x)传输至计算机图形计算卡GPU的显卡显存中；在三维体数据重建过程中，计算机图形计算卡先将该投影图像对应的系数传输至GPU的显卡显存中，然后采用并行运算模块对该投影图像中的像素进行加权、滤波，最后进行反投影操作。所述计算机图形计算卡中还可以包括：纹理存储器，在所述反投影操作中使用纹理拾取函数访问纹理存储器，以加速反投影过程。

Claims (9)

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

第一步，预先计算出锥形束FDK重建所需的权重系数，包括：h(x)和

其中，h(x)是滤波核函数，(x，y)是位于旋转中心的虚拟探测器平面的像素位置，D_SO是射线源到虚拟探测器平面的距离，s是待重建像素相对于虚拟探测器平面的距离坐标；

第二步，在进行锥形束FDK重建过程前，将第一步预先计算好的这些权重系数加载至主机内存中，并将和h(x)传输至GPU显存中；

第三步，为GPU并行化三维体数据重建创建操作流；

第四步，在三维体数据重建过程中，GPU依次对每幅投影图像进行处理：先将该投影图像对应的系数

传输至GPU显存中，然后使用已预先传输至GPU中的和h(x)对投影图像进行加权、滤波操作，最后进行反投影操作。

2.如权利要求1所述的方法，其特征在于，第一步中，所述的预先计算，是指在使用GPU进行CT并行重建前将相关的权重系数离线计算好。

3.如权利要求1所述的方法，其特征在于，第二步中，所述的在进行锥形束FDK重建过程前，将第一步预先计算好的这些权重系数加载至主机内存中，是指将第一步预先计算好的权重系数加载至主机内存中，获得其内存访问地址入口；所述的将

和h(x)传输至GPU显存中是指在GPU显存中申请同样大小的空间，将

和h(x)存放至GPU显存中，供GPU重建过程调用。

4.如权利要求1所述的方法，其特征在于，第三步中，所述的为GPU并行化三维体数据重建创建操作流，是指在GPU重建前声明一个流，在GPU进行并行化重建过程中，将第四步中所述的数据传输和操作加载至这个流中。

5.如权利要求1所述的方法，其特征在于，第四步中，所述的先将该投影图像对应的系数

传输至GPU显存中，是指：把投影图像写入GPU的CUDA数组，并且只将与该投影图像相关的系数

传输至GPU显存，并将步骤记录在第三步创建的操作流中，GPU并行化的基本计算单元是投影图像中的像素。

6.如权利要求1所述的方法，其特征在于，第四步中，所述的使用已预先传输至GPU中的和h(x)对投影图像进行加权、滤波操作，是指使用GPU的多个线程并行地对该投影内的多个图像像素进行加权、滤波，其中使用已经预先计算好，并在第二步中传输至GPU显存中的系数，该操作记录在第三步创建的流中。

7.如权利要求1所述的方法，其特征在于，第四步中，所述反投影操作中，将滤波后的投影图像映射至纹理显存中，并进行反投影操作，具体指：将滤波后的投影图像与GPU的纹理坐标系进行绑定，然后使用纹理拾取函数访问纹理存储器，加速锥形束CT重建的反投影过程，该操作记录在第三步创建的流中。

8.一种锥形束CT重建系统，其特征在于，该系统包括：计算机图形计算卡和计算机主机；所述计算机主机中包括预先计算权重系数模块；所述计算机图形计算卡包括显卡显存、并行运算模块；

所述预先计算权重系数模块用于预先计算出锥形束FDK重建所需的权重系数，这些权重系数加载至所述计算机主机内存中，并将

和h(x)传输至计算机图形计算卡GPU的显卡显存中；在三维体数据重建过程中，计算机图形计算卡先将该投影图像对应的系数

传输至GPU的显卡显存中，然后采用并行运算模块对该投影图像中的像素进行加权、滤波，最后进行反投影操作。

9.如权利要求8所述的系统，其特征在于，所述计算机图形计算卡中还包括：纹理存储器，在所述反投影操作中使用纹理拾取函数访问纹理存储器，以加速反投影过程。

Images