CN102279970A - 基于gpu的螺旋锥束ct重建方法 - Google Patents

Abstract

公开了一种基于GPU的螺旋锥束CT重建方法，其中，在GPU中执行以下步骤：对投影数据进行预处理；对预处理后的数据进行三维反投影，得到重建的图像数据。其中预处理的步骤可包括求导、重排、滤波等操作。通过GPU实现上述方法，能大大提高重建速度，从而使在实际应用中进行实时重建成为可能。

Description

基于GPU的螺旋锥束CT重建方法

技术领域

本发明的实施例涉及计算机断层成像(CT：ComputerTomography)技术，具体涉及一种基于图形处理器(GPU：GraphicProcessor Unit)的螺旋锥束CT重建方法。

背景技术

在CT应用领域，Katsevich型重建方法作为螺旋CT精确重建算法中的一种，其重建效果已被广泛认可。但由于其计算量巨大，限制了在实际系统中的应用。

通用GPU计算技术是近几年兴起的快速计算技术由不断增长的市场对实时三维图形处理的需求。图形处理器(Graphic ProcessingUnit，GPU)以惊人的速度在发展。从固定的渲染管道发展到可编程的非固定渲染管道的专用处理器，发展到现在的可用于通用计算目的的通用并行处理器GPGPU(General Purpose GPU)，其应用的范围也逐渐从最初的加速图形几何变换、光照和纹理等，发展到可编程的顶点和片段控制，最终成为了可满足常见计算需求的通用计算设备。

GPU的浮点计算能力也以远超CPU的速度在增长，到08年底NVIDIA GT200核心的显卡达到了约1T FLOPs，而同时期的Intel CPUHarpertown浮点运算能力约65G FLOPs。(1T＝1000G)NVIDIDA即将推出的GT300构架显卡预计将达到2.7TFLOP/s，远超CPU的计算能力。

设计目的不同是造成CPU和GPU之间的差异主要原因：GPU为了完成复杂的图形渲染工作，需要强大的并行计算能力，因此将更多晶体管用于计算单元；CPU为了更为通用，则不断的增加缓存和流程控制单元。CPU主要用于不规则的数据结构，其存取模式也就变得不可预测，因此只有不断增加缓存来解决存储器延迟问题。而GPU使用单程序多数据执行模式(SPMD)、规则的数据结构，使得存储可预测，从而可以有更多的空间用于增加计算单元。另外，GPU通过提升计算通信比例来隐藏存储器延时。GPU中同时运行大量线程，等待内存访问的线程马上被切换，而GPU中线程切换只需要一个时钟周期，相对的CPU切换线程则要几百个时钟周期。

从图1中可以看到，GPU以更多的算术逻辑单元(Arithmetic LogicUnit，ALU)用于计算，相应的牺牲了缓存(Cache)和程序流控制能力。因此，如果能充分的利用GPU的优点，合理的避开其缺点，将大幅提高程序的运行效率，提高单机运算能力。

发明内容

本发明的目的是提出一种基于图形处理器的螺旋锥束CT重建方法。

在本发明的一个方面，提出了基于GPU的螺旋锥束CT重建方法，其中的反投影步骤在GPU中执行时，为每个待重建像素分配单独线程来执行，重建所需预处理后的投影数据存储在GPU的纹理内存中，其中线程分配过程包括以下步骤：根据GPU的特性设置每个线程块(Block)的尺寸；根据待重建图像的尺寸设置线程块(Block)的个数；按照线程块(Block)设置执行内核程序。

根据本发明的实施例，对投影数据的预处理步骤包括：对投影数据求导；对求导后的数据进行加权；将加权的数据沿着k线方向重排；对重采样的数据进行滤波；对滤波后的数据从k线方向重排到u方向。

根据本发明的实施例，使用stream技术，从而隐藏数据读取时间。

根据本发明的实施例，使用通过修改函数库中函数的精度来提高运算速度。

通过GPU实现上述方法，能大大提高重建速度，从而使在实际应用中进行实时重建成为可能。

附图说明

图1示出了CPU与GPU结构对比图。

图2A是平面探测器螺旋扫描示意图。

图2B是描述Tam窗的示意图。

图2C是描述k线的示意图。

图3是描述Katsevich算法推导坐标示意图。

图4是根据本发明实施例的基于GPU来实现Katsevich重建算法的流程图。

图5示出了根据本发明实施例的图像像素与线程对应示意图。

图6示出了根据本发明实施例的基于GPU的滤波过程的流程图。

图7示出了根据本发明实施例的三维体数据重建线程分配示意图。

图8示出了作为比较用的基于CPU的执行加权过程的计算流程。

图9示出了根据本发明实施例的GPU执行加权过程流程图。

图10示出了根据本发明实施例的基于GPU的MicroCT系统数据处理流程图。

图11示出了根据本发明实施例的使用stream技术时间对比。

具体实施方式

根据本发明的实施例，提出了在GPU中实现Katsevich算法，以便充分利用GPU和Katsevich算法的特点，使得在单机上实时进行精确的CT重建称为可能。

Katsevich算法是第一个具有滤波反投影形式，且使用一维移不变滤波核的精确螺旋锥束重建算法。

螺旋锥束CT系统一般由固定在可旋转机架上的光源S的和一个面阵列探测器组成，一般探测器为平面探测器或是柱面探测器，但算法原理并不局限于此，因此本实施例中所使用的公式均以平面探测器为例。本领域的技术人员很容易将平面探测器扩展到其他类型探测器，并能成功实施本案。

将待成像物体放置在光源和探测器之间，如图2A所示。

以固定在物体上的坐标系(x，y，z)为基准坐标系，并假设物体函数的支撑被限制在一个半径为R₀的圆柱内。

$f(x,y,z)=0{,x}^2+y^2>R_0^2---2-1$

光源和探测器则以z轴为旋转轴围绕物体旋转，光源S到探测器中心的连线垂直于探测器平面，且距离为D。在光源S和探测器旋转的同时，被扫描物体匀速沿z轴方向运动，因此从物体上看，光源的运动轨迹为螺旋线。设光源到旋转轴的距离为R，光源每旋转一周物体移动的距离为h，则光源的位置可表示为：

$\stackrel{\→}{S}=(R*\cos \mathrm{\β},R*\sin \mathrm{\β},\frac{\mathrm{\β}}{2\mathrm{\π}}h+z_0)---2-2$

其中β为光源相对于起始位置的旋转角，z₀为光源的起始位置。另外，在本发明的实施例中将螺距定义为h。为了表述该算法，首先引入了κ平面，该平面按如下方式定义。

令κ平面与螺旋线相交于三个点

其中

该平面的法线为：

Katsevich证明了对于给定的点

如果该点在κ平面内，并且光源

与

的连线与探测器的交点落于Tam窗内，则κ平面被唯一确定，即

完全确定。该Tam窗由能够在径向可以完全包围被扫描物体、z方向刚好覆盖两条螺旋线的半径为R的柱面探测器在平面探测器的投影组成，如图2B所示。κ平面的物理含义可理解为，该平面为与螺旋线有三个交点，并且包含光源和待重建点的平面，Katsevich证明了如果光源

与待重建点

的连线与探测器的焦点落于Tam窗内，则该平面是唯一的。定义κ平面与探测器平面的交线为κ线，如图2C所示。

定义光源到重建点的单位向量为：

$\widehat{\mathrm{\ω}}(\mathrm{\β},\stackrel{\→}{r})=\frac{\stackrel{\→}{r}-\stackrel{\→}{S}\left(\mathrm{\β}\right)}{\left|\right|\stackrel{\→}{r}-\stackrel{\→}{S}\left(\mathrm{\β}\right)\left|\right|}---2-29$

另外定义

为包含穿过

的方向为

的射线的κ平面的法向量，则可定义κ平面内另外一个单位向量：

$\hat{e}(\mathrm{\β},\stackrel{\→}{r})=\widehat{\mathrm{\ω}}(\mathrm{\β},\stackrel{\→}{r})\×\hat{m}(\mathrm{\β},\widehat{\mathrm{\ω}})---2-30$

经过以上定义后，Katsevich算法的重建公式可如下表述：

$f(x,y,z)=-\frac{1}{2{\mathrm{\π}}^2}{\∫}_{I_{\mathrm{PI}}\left(\stackrel{\→}{r}\right)}\frac{1}{\left|\right|\stackrel{\→}{r}-\stackrel{\→}{S}\left|\right|}\mathrm{PV}{\∫}_0^{2\mathrm{\π}}\frac{d}{d{\mathrm{\β}}^{\′}}D\left[f(\stackrel{\→}{S}\left({\mathrm{\β}}^{\′}\right),\widehat{\mathrm{\θ}}(\mathrm{\β},\stackrel{\→}{r},\mathrm{\γ}))\right]{|}_{{\mathrm{\β}}^{\′}=\mathrm{\β}}\frac{\mathrm{d\γ}}{\sin \mathrm{\γ}}\mathrm{d\β}---2-31$

其中

而D(f)定义为：

$D\left[f(\stackrel{\→}{S},\widehat{\mathrm{\θ}})\right]={\∫}_0^{\∞}f(\stackrel{\→}{S}+t\widehat{\mathrm{\θ}})\mathrm{dt}---2-32$

易知D(f)的物理含义为κ平面内，经过光源角度为

的投影值。

探测器平面内的κ线可由如下公式计算：

图3是描述Katsevich算法推导坐标示意图。

定义探测器坐标系(u，v，w)：

${\hat{e}}_u\left(\mathrm{\β}\right)={[-\sin \mathrm{\β},\cos \mathrm{\β},0]}^T$

${\hat{e}}_v\left(\mathrm{\β}\right)={\left[\mathrm{0,0,1}\right]}^T---2-34$

${\hat{e}}_w\left(\mathrm{\β}\right)={[-\cos \mathrm{\β},-\sin \mathrm{\β},0]}^T$

将式2-32用单位向量表示，则探测器的投影可以表述为：

$D\left[f(\stackrel{\→}{S},\widehat{\mathrm{\ω}}(\mathrm{\β},\stackrel{\→}{r}))\right]=P(\mathrm{\β},u^{*},v^{*})---2-35$

(u^*，v^*)可以通过如下方式计算：

$u^{*}(\mathrm{\&beta;},\stackrel{\&RightArrow;}{r})=D\frac{<\widehat{\mathrm{\&omega;}},{\hat{e}}_u>}{<\widehat{\mathrm{\&omega;}},{\hat{e}}_w>}=\frac{D}{w^{*}(\mathrm{\&beta;},\stackrel{\&RightArrow;}{r})}(-x\sin \mathrm{\&beta;}+y\cos \mathrm{\&beta;})$

$v^{*}=D\frac{<\widehat{\mathrm{\&omega;}},{\hat{e}}_v>}{<\widehat{\mathrm{\&omega;}},{\hat{e}}_w>}=\frac{D}{w^{*}(\mathrm{\&beta;},\stackrel{\&RightArrow;}{r})}(z-\frac{\mathrm{\&beta;}}{2\mathrm{\&pi;}}h)---2-36$

$w^{*}(\mathrm{\&beta;},\stackrel{\&RightArrow;}{r})=R+<\stackrel{\&RightArrow;}{r},{\hat{e}}_w>=R-x\cos \mathrm{\&beta;}-y\sin \mathrm{\&beta;}$

从形式上，式2-31可以分为三步来计算：

1、对投影数据进行求导：

$g_1(\mathrm{\&beta;},u,v)=\frac{d}{\mathrm{d\&beta;}}P(\mathrm{\&beta;},u,v)$

2-37

$g_1(\mathrm{\&beta;},u,v)=(\frac{\&PartialD;}{\&PartialD;\mathrm{\&beta;}}+\frac{D^2+u^2}{D}\frac{\&PartialD;}{\&PartialD;u}+\frac{\mathrm{uv}}{D}\frac{\&PartialD;}{\&PartialD;v})P(\mathrm{\&beta;},u,v)$

2、对求导之后的数据进行Hilbert滤波，滤波的方向沿着κ线方向，可拆分为4步：

a.对投影数据g₁(β，u，v)加权：

$g_2(\mathrm{\&beta;},u,v)=g_1(\mathrm{\&beta;},u,v)\frac{D}{\sqrt{D^2+u^2+v^2}}---2-38$

b.对投影数据g₂(β，u，v)沿κ线方向重采样(式2-33)：

c.对投影数据做Hilbert滤波：

d.对投影数据沿u方向进行重采样：

此处

仍需满足式2-33。如图2A和2B所示，u方向指平面探测器的水平方向。

3、滤波之后即可进行三维反投影：

$f(x,y,z)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{{\mathrm{\&beta;}}_1}^{{\mathrm{\&beta;}}_2}\frac{1}{w^{*}(\mathrm{\&beta;},\stackrel{\&RightArrow;}{r})}{g}_5(\mathrm{\&beta;},u^{*},v^{*})\mathrm{d\&beta;}---2-42$

图4是根据本发明实施例的基于GPU来实现Katsevich重建算法的流程图。

根据本发明的实施例，反投影步骤在GPU中执行时，为每个待重建像素分配单独线程来执行。重建所需预处理后的投影数据存储在GPU的纹理内存中。线程分配过程包括：根据GPU的特性设置每个线程块(Block)的尺寸；根据待重建图像的尺寸设置线程块(Block)的个数；按照线程块(Block)设置执行内核程序。

如图4所示，在步骤S40，GPU读取各个投影角度下的投影数据。在步骤S41，GPU将投影读入显存中。然后，在步骤S42，按照式2-37和2-38对投影数据做求导运算和加权运算。

以每一帧微分后的投影数据为处理对象，设虚拟探测器尺寸为M×N，则可将整个图像分为X×Y块(Block)，每个块中包含m×n个线程(Thread)，而每个线程负责计算对应的像素P(u，v)。在CUDA(Compute Unified Device Architecture，一种计算机编程语言，可用于GPU编程)中可以通过内置变量来计算对应(u，v)：

blockDim：block的尺寸(m×n)

blockIdx：block在网格中的索引(x，y)

threadIdx：为线程在block内的索引(i，j)

则线程与像素的对应关系可按下式计算：

u＝blockDim.x^*blockIdx.x+threadIdx.x；

                                         3-1

v＝blockDim.y^*blockIdx.y+threadIdx.y；

而在OPENCL(Open Computing Language，一种计算机编程语言，可用于GPU编程)中则可以通过内置函数get_global_id(k)；直接得到，其中k为维度：

u＝get_global_id(0)；

                          3-2

v＝get_global_id(1)；

图5示出了根据本发明实施例的图像像素与线程对应示意图。

在步骤S43，GPU通过式2-39对微分后的数据进行重排到k线方向，重排时需要根据式2-33计算出v_k。另外该步骤也可以通过预先计算插值位置及系数查找表的方式实现。一般v_k不会落在整数像素上，为了得到更好的重建效果，需要做线性插值。为此，我们引入二维纹理协助完成该操作。首先将微分后投影数据拷贝到纹理内存，绑定为二维纹理，然后将计算得出的物理坐标转换为像素坐标，然后就可以通过使用纹理拾取函数来快速实现该操作。

在步骤S44，GPU执行滤波运算。滤波属于卷积运算，可以通过快速傅立叶变换(FFT)实现。图6示出了根据本发明实施例的基于GPU的滤波过程的流程图。

根据Nyquist采样定理可知，采样的频率需要是原信号的两倍才能将其复原。因此在进行FFT之前，需要将每行投影数据补0。使用FFT函数库进行滤波的流程如图6所示。

在步骤S60，在进行FFT之前，需要将每行投影数据补0。

在步骤S61，对投影数据进行FFT，将其转换成频域数据。

在步骤S62，对频域将频域数据域频域滤波核相乘，以进行频域滤波。

在步骤S63，将频域滤波的数据进行IFFT，将其转换到空间域，完成滤波。

在步骤S45，GPU执行逆向重排，将滤波后的投影数据从k线方向重排到u方向时，是式2-33的逆向操作。由于该式无法直接逆向求解，所以可采用计算查找表的方式实现求解。同步骤S43类似，在步骤S46，将滤波后的投影数据拷贝到纹理内存，协助完成插值操作，提高存取速度。

在步骤S47，GPU执行反投影过程。反投影过程主要是根据2-36公式计算出反投影地址。得到该地址后将该点的投影值加权后累加到待重建点f(x，y，z)上即可。由上一章中我们知道反投影计算的复杂度为O(N⁴)，而每一次反投影还对应着一次读取，如果使用双线性插值则是四次。因此我们将投影数据存储到有缓存的纹理内存中实现双线性插值。

由于反投影时需要计算的点数是三维尺寸的X×Y×Z，因此线程管理稍微复杂些，见图7。可以将整个重建图像理解为宽度为X，高度为Y×Z的二维图像再进行块的划分。而在程序内部根据高度方向的索引来计算z坐标，易知：

x＝blockDim.x^*blockIdx.x+threadIdx.x

y_temp＝blockDim.y^*blockIdx.y+threadIdx.y

$z=\mathrm{floor}\left(\frac{y_{\mathrm{temp}}}{Y}\right)---3-3$

y＝y_temp-z^*Y；

其中floor表示向下取整运算。

在步骤S48，GPU判断是否已经针对所有的投影角度执行了上述的重建过程，如果否的话，流程转到步骤S41，继续上述过程。如果已经针对所有的投影角度执行了上述的重建过程，在步骤S49，将重建后的图像数据从显存读到内存中。

下面以步骤2-38为例，说明CPU与GPU在计算过程中的不同。

$g_2(\mathrm{\&beta;},u,v)=g_1(\mathrm{\&beta;},u,v)\frac{D}{\sqrt{D^2+u^2+v^2}}---2-38$

对微分后的投影数据g₁(β，u，v)，其中(0≤u＜M，0≤v＜N)，在CPU中求解该方程流程如下：

在步骤S80，首先进行初始化工作，并将探测器行方向索引v置为0，即从第0行数据开始处理。

在步骤S81，将探测器列方向索引u置为0，即从第0列数据开始处理。

在步骤S82，需要根据当前的u，v值，利用式2-38中右侧部分来计算加权系数

在步骤S83，使用步骤S82中计算所得加权系数，对g₁(β，u，v)进行修正，得到g₂(β，u，v)完成对一个像素单元的加权。探测器数据列方向索引u指向下一个单元。

在步骤S84，判断探测器数据列方向索引u是否小于探测器列数M，如果是则循环步骤S82，S83；如果不是则执行步骤S85。

在步骤S85，将行方向索引v指向下一行。

在步骤S86，判断探测器数据行方向索引v是否小于探测器行数N，如果是则循环步骤S81，S82，S83，S84，S85；如果不是则执行步骤S87。

在步骤S87，所有像素计算完毕。

可见在，上述基于CPU的过程中，必须进行两层嵌套循环。

但是，在GPU中，以每一帧微分后的投影数据为处理对象，将整个图像分为X×Y块(Block)，每个块中包含m×n个线程(Thread)，而每个线程负责计算对应的像素g₂(β，u，v)。图9示出了根据本发明实施例的GPU执行加权过程流程图。

在步骤S90，进行初始化。在步骤S91，按照如上所述的方法分配线程参数。在步骤S92，在线程内得到u，v，并计算加权系数，然后对投影数据进行修正。最后，在步骤S93，完成整个加权过程。

从图9中可以看出，GPU中的执行过程省去了两层循环过程，将整幅图的加权过程一步完成，充分利用GPU多线程并行的能力，从而大幅度提高计算速度。

图10示出了根据本发明实施例的基于GPU的MicroCT系统数据处理流程图。

根据本发明实施例的CT系统由若干子系统构成，比如扫描控制系统，电器系统，机械系统，数据采集系统和数据处理与可视化系统等。其中数据处理系统负责数据的校正、预处理、重建、后处理和显示等，其处理流程如图10所示。

在通用GPU(GPGPU：General Purpose GPU)计算技术出现以前，只有步骤S105由显卡来完成，而现在S100-S105均可由GPU来实现，并且由于CT系统本身数据适于并行处理，性能比传统的使用CPU有较大提高。另外由于整条流程都在显卡内部完成，则可以通过计算的时间隐藏掉与系统内存的IO时间，从而进一步提高性能。

根据本发明的实施例，重建软件支持两种重建模式：在线重建和离线重建。在线重建是指在整个系统的工作流程过程中，每采集一帧投影数据，就处理一帧。离线重建是指已经有全部经过校正和预处理的投影数据，直接选择进行重建。两种重建方式区别在于，在线重建由于系统采集较慢，约2.5fps，即400ms。对于GPU来说，单个投影角度的重建时间远小于采集时间，相当于有大部分时间在等待投影数据。因此无需对IO进行优化。而对于离线重建来讲，由于数据直接从硬盘读取，无需等待采集系统，因此如何优化IO就成为关键。

根据本发明的实施例，在在线重建过程中使用与系统相关的流技术(Stream)来进行优化IO。

设IO操作和计算有先后关系，并且有相同的量级。在CUDA中，IO API和内核程序都支持同步和异步执行。同步是指在函数执行时，线程将阻塞，必须等待该函数执行完毕再执行下一个函数。而异步函数则调用后直接返回，可继续执行下一个函数。在CUDA中，我们可以使用stream来实现批量的数据存储和内核执行。

假设有N个float数据需要读取并由kernel处理，使用stream处理过程的代码如下：

//给stream分配内存

cudaStream_t^＊streams＝

         (cudaStream_t^＊)malloc(nStreams^＊sizeof(cudaStream_t))；

//创建stream

for(int i＝0；i＜nStreams；i++){

    cudaStreamCreate(&(streams[i]))；

}

//每次读取的字节数

int size＝N^＊sizeof(float)/nStreams；

//分nStreams次从host拷贝数据，每次size个字节

for(int i＝0；i＜nStreams；i++){

    offset＝i^＊N/nStreams ；

    //d为设备内存指针，h为主机内存指针，dir表示从Host-＞Device

    cudaMemcpyAsync(d+offset，h+offset，size，dir，stream[i])；

}

for(i＝0；i＜nStreams；i++){

    offset＝i^＊N/nStreams；

    //将处理数据的内核放入对应的stream

    kernel<<<N/(nThreads^＊nStreams)，nThreads，0，stream[i]>>>

    (d+offset)；

}

通过使用stream，可以将数据读取的时间部分隐藏(见图11)，从而进一步提高程序执行效率。在本系统中，由于计算量较大，可以通过stream技术将IO时间全部覆盖。

根据本发明的实施例，可以按照探测器A/D的精度来设定或修改重建过程中的使用的库函数的精度，使其与探测器A/D的精度匹配，从而避免冗余计算过程。例如，在该系统中，使用的探测器A/D转换位数为12位，即输出值范围为[0，4095]。这意味着，最高的有效精度为

因此我们可以通过降低运算精度来提高计算速度。在CUDA中，提供了内置的低精度数学函数，例如：

t＝__sinf(x)                    4-12

用来计算x的正弦值，当x∈[-π，π]时，最大误差为2^-21.41。其他同类型低精度数学函数的运算精度，都远高于系统探测器的精度，因此完全可以使用低精度函数库。本实施例中用到的低精度数学函数及其误差见表1。以PI算法的重排步骤为例，通过使用低精度函数库可以将速度提高15％。

表1CUDA低精度函数误差

虽然以上以Katsevich算法为例描述了CT重建算法在GPU中的实现，但是本领域的技术人员也可以将如PI Original算法，OS-EM算法，BPF，FDK类，ART类算法在GPU中实现。

通过GPU实现上述重建算法，能将重建速度提高50～100倍，从而使该算法在实际应用中进行实时重建成为可能。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可轻易想到的变换或替换，都应涵盖在本发明的包含范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (4)

1.一种基于GPU的螺旋锥束CT重建方法，其中的反投影步骤在GPU中执行时，为每个待重建像素分配单独线程来执行，重建所需预处理后的投影数据存储在GPU的纹理内存中，其中线程分配过程包括以下步骤：

根据GPU的特性设置每个线程块(Block)的尺寸；

根据待重建图像的尺寸设置线程块(Block)的个数；

按照线程块(Block)设置执行内核程序。

2.根据权利要求1所述的方法，其中对投影数据的预处理步骤包括：

对投影数据求导；

对求导后的数据进行加权；

将加权的数据沿着k线方向重排；

对重采样的数据进行滤波；

对滤波后的数据从k线方向重排到u方向。

3.根据权利要求1所述的方法，其中使用stream技术，从而隐藏数据读取时间。

4.根据权利要求1所述的方法，其中使用通过修改函数库中函数的精度来提高运算速度。

Images