CN109146987B - 一种基于gpu的快速锥束计算机断层成像重建方法 - Google Patents

Abstract

本发明属于用于放射诊断的仪器技术领域，公开了一种基于GPU的快速锥束计算机断层成像重建方法，包括：获取多个角度的物体测量数据；获取计算机断层成像系统的各个参数；对测量数据做预处理；对预处理后的数据做加权滤波，基于GPU的重排滤波算法；对滤波过后的数据做反投影，基于GPU的反投影使用纹理映射理论以及FDK算法的本身对称性实现并行加速，完成计算机的断层成像。本发明基于滤波反投影的重建方法，利用系统参数等先验信息，采用X射线照射，有限角度测量，利用GPU策略，实现了在不改变原有的图像质量基础上有效缩小了重建的时间，在计算机断层三维锥束重建算法等领域有重要的应用价值。

Description

一种基于GPU的快速锥束计算机断层成像重建方法

技术领域

本发明属于用于放射诊断的仪器技术领域，尤其涉及一种基于GPU的快速锥束计算机断层成像重建方法。

背景技术

目前，业内常用的现有技术是这样的：作为二十一世纪影响人类发展的十大技术之一，计算机断层成像技术(Computed Tomography，即CT)是一种在各个角度下利用X射线穿透物体得到的投影信息反求物体密度分布的成像技术。该技术的主要理论是数学和核物理学，在此基础上又集成了自动控制技术、计算机技术和探测器技术等多个学科，且其已广泛应用于各个重要领域，如生物研究、医疗诊断、公共安全、工业检测，国防建设等，尤其是作为获取人体内部结构信息的最佳手段在临床诊断中应用。其中，锥形束CT技术，则以扫描速度快、空间分辨率高等特点已成为目前CT研究工作的重点和热点。锥束CT(Cone-BeamComputed Tomography，CBCT)重建和可视化一直是三维重建软件中最关键的部分。扫描仪对敏感组织以及病变体进行扫描并重建出三维模型，为了使医学诊断能更准确和快速的放射治疗和病变定位，利用可视化技术为医生提供需要直接可见的三维模型部位或者整体。近年来，重建算法的加速以及其可视化已逐渐成为研究热点。锥形束CT扫描获得的数据量庞大，算法相对复杂，重建时间较长，无法满足实际应用中实时性的要求。近半个世纪以来，随着图形处理器(GPU，GraphicsProcessing Unit)的迅速发展，其强大的并行运算能力对CT重建算法的加速产生了革命性的意义。尤其是，CUDA架构的提出，促使GPU的可编程性增加，方便了开发者快速掌握其应用，实时性的要求是重建的必须因素，但现阶段的加速算法速度并不是很理想，因此实现锥形束CT重建算法的加速是很有必要的。

综上所述，现有技术存在的问题是：锥形束CT扫描获得的数据量庞大，算法相对复杂，重建时间较长，无法满足实际应用中实时性的要求。

解决上述技术问题的难度和意义：锥束CT重建算法的复杂度高，根据算法以及硬件对其进行加速，可以使观察人员实时观测物体的内部结构。

发明内容

针对现有技术存在的问题，本发明提供了一种基于GPU的快速锥束计算机断层成像重建方法。

本发明是这样实现的，一种基于GPU的快速锥束计算机断层成像重建方法，所述基于GPU的快速锥束计算机断层成像重建方法包括：获取多个角度的物体测量数据；获取计算机断层成像系统的各个参数；对测量数据做预处理；对预处理后的数据做加权滤波，基于GPU的重排滤波算法；对滤波过后的数据做反投影，基于GPU的反投影使用纹理映射理论实现并行加速，最终完成计算机的断层成像；

所述对预处理后的数据做加权滤波，基于GPU的重排滤波算法；对滤波过后的数据做反投影具体包括：

(1)对像素和源点之间的角度差和距离进行修正，对处理过后原始数据进行加权，使用并行算法对每个角度的测量数据同时处理：

g'(u,v,θ)＝cos(λ)g(u,v,θ)；

其中g(u,v,θ)为测量数据,

R(u,v)为射线源到旋转中心的距离，(u,v)为探测器上物体的位置；

(2)数据在频域空间进行滤波，对加权后的数据进行滤波，使用CUFFT库，并且进行数据重排，复数和虚数分别占一行数据；

g”(u,v,θ)＝g'(u,v,θ)h(u)；

其中h(u)为滤波核函数；

(3)数据进行加权反投影对滤波过后的数据按照X射线源的投影方向做加权反投影，重建最终的结果就是把所有角度的射线贡献进行求和。

进一步，所述基于GPU的快速锥束计算机断层成像重建方法进一步包括：

1)反投影时的权值计算

只需要计算0～90°的权值，其中R为射线源到旋转中心的距离，U²(x,y,θ)为射线源到物体的距离；

2)将θ，θ+90°，θ+180°和θ+270°的四个旋转对称视图中的投影图像作为一组进行排列，并打包放入GPU的2D纹理的四个颜色通道，每个通道存放一个投影图像；

3)根据纹理映射矩阵，由顶点着色器计算每一层切片四个顶点的纹理映射坐标：

U(x,y,θ)＝w(x,y,θ)；

其中射线源到探测器距离为D，射线源到旋转中心的距离为R，射线源S绕着中心轴z进行旋转，扇形角为α，锥角为β，从射线源出发过点A的射线与中心射线的夹角为λ，射线源到物体中心与虚拟坐标系x轴之间的夹角为θ；

4)在光栅化器中，矩形内部点的纹理映射坐标是由GPU自动完成插值计算；

5)在片元着色器中，根据纹理映射坐标取出相应的投影数据，完成加权累加；

其中U(x,y,θ)＝R+xcosθ+ysinθ；

其中射线源到旋转中心的距离为R，射线源S绕着中心轴z进行旋转，投影角度为θ；

6)写入帧缓存，四个角度同时进行渲染，M为体积切片数，每一层切片都由一个纹理内存来保存，每个纹理中都有一个单独的四个颜色通道，并且每个通道用于从一个投影视图保存反投影值，利用坐标FDK算法的垂直对称性同时计算1层和M层的反投影数据。

进一步，所述获取多个角度的物体测量数据利用计算机断层成像系统，依次采集测量，其中测量数据的获取为有限个角度，为绕物体一圈，隔1度获取一次，总共为360组。

进一步，所述获取计算机断层成像系统的各个参数，对测量数据做预处理具体包括：

(1)使用计算机断层成像系统来进行重建时，需要获取系统的各个参数：射线源到中心物体的距离、射线源到探测器之间的距离，探测器的分辨率、探测器的偏移度；

(2)对数据做预处理时，测量断层成像系统做没有物体时的数据，对和(1)中的测量数据做对数运算。

本发明的另一目的在于提供一种应用所述基于GPU的快速锥束计算机断层成像重建方法的计算机断层成像系统。

综上所述，本发明的优点及积极效果为：本发明采用的基于GPU锥束重建，FDK算法的三个步骤(加权、滤波、反投影)分别使用计算设备编程架构进行了加速。加权时，并行处理每个角度的数据；滤波时，使用计算设备编程架构中的CUFFT库，重排加权后的数据，减少无效的运算和存储空间的使用；反投影时，利用FDK算法本身的对称性(垂直对称、旋转对称)以及计算设备编程架构的四纹理通道提高了算法的运算速度。在不改变原有的图像质量基础上有效缩小了重建的时间。

CT系统参数	小鼠一	小鼠二
			探测器大小	506×516	1944×1535
旋转角度	360	360
			探测器分辨率	0.1	0.074
探测器偏移	14.7193×51.91535	271.4419×-19.3603
			射线源到探测器距离	482.2066	632.1605
射线源到物体中心距离	405.7135	418.7630

重建尺寸大小	CPU重建时间	原GPU重建时间	GPU重建时间
				512×512×512(小鼠一)	1575s	17.58s	9.73s
512×512×512(小鼠二)	6358s	64.32s	27.45s

附图说明

图1是本发明实施例提供的基于GPU的快速锥束计算机断层成像重建方法流程图。

图2是本发明实施例提供的基于GPU的快速锥束计算机断层成像重建方法实现流程图。

图3是本发明实施例提供的用于仿真实验的数字鼠模型截面显示示意图。

图4是本发明实施例提供的用于仿真实验的数字鼠模型冠面显示示意图。

图5是本发明实施例提供的矢面表示示意图。

图6是本发明实施例提供的重建算法获得的重建结果；

图中：(a)图显示为506×516×360的探测数据重建为512×512×512的三维展示；(b)图显示为1944×1535×360的探测数据重建为512×512×512的三维展示。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明结合CUDA架构的存储结构以及FDK算法自身的特性，实现了基于GPU的FDK算法，对重建进行了加速本发明提出了基于GPU策略的有限投影角度的计算机断层成像重建方法；在不改变原始图像的重建质量上，为了实现实时重建，采用基于GPU的重建，对锥束计算机断层成像的每个步骤分别进行了加速。

下面结合附图对本发明的应用原理作详细的描述。

如图1所示，本发明实施例提供的基于GPU的快速锥束计算机断层成像重建方法包括以下步骤：

S101：获取多个角度的物体测量数据；

S102：获取计算机断层成像系统的各个参数；

S103：对测量数据做预处理；

S104：对预处理后的数据做加权滤波，基于GPU的重排滤波算法；

S105：对滤波过后的数据做反投影，基于GPU的反投影使用纹理映射理论实现并行加速，最终完成计算机的断层成像。

下面结合附图对本发明的应用原理作进一步的描述。

如图2所示，本发明实施例提供的基于GPU的快速锥束计算机断层成像重建方法具体包括以下步骤：

(1)获取有限角度的断层数据：

1a)X射线源对固定在电控旋转台上的重建目标进行有限角度的计算机断层成像；

计算机断层成像，X射线束对人体某部一定厚度的层面进行扫描，由探测器接收透过该层面的X射线，转变为可见光后，由光电转换变为电信号，再经模拟/数字转换器(analog/digital converter)转为数字，输入计算机处理。有限个角度计算机断层扫描，用电脑控制电控旋转台等间隔旋转一定角度，一般是间隔1°，在平板探测器上得到投影的图像，一般转一个角度激发一次，采集一次测量数据。

(2)获得系统的物理信息：

2a)系统的物理信息以及对数据进行预处理

对计算机断层成像系统的数据进行采集，主要有；探测器的大小、旋转角度、探测器的分辨率、探测器的偏移量、射线源到探测器的距离以及射线源到物体中心的距离。

2b)数据预处理

对数据做预处理时，测量断层成像系统做没有物体时的数据，对其和(1)中的测量数据做对数运算。

(3)基于反投影重建理论，使用FDK算法进行锥束CT的重建，重建分成三步，加权、滤波、反投影，最终实现三维物体的重建。将上述中整个重建过程放在GPU中进行，对加权、滤波、反投影三个步骤分别做加速；

3a)加权时，对像素和源点之间的角度差和距离进行修正，对处理过后原始数据进行加权，只有加权因子影响着加权的整个过程，可以使用并行算法对每个角度的测量数据同时处理：

g'(u,v,θ)＝cos(λ)g(u,v,θ)；

其中g(u,v,θ)为测量数据,

R(u,v)为射线源到旋转中心的距离，(u,v)为探测器上物体的位置。

3b)滤波时，可以使用CUFFT库，并且进行数据重排，复数和虚数分别占一行数据，再在滤波完成后进行还原，实现滤波时间的减半；

g”(u,v,θ)＝g'(u,v,θ)h(u)；

其中h(u)为滤波核函数，一般采用斜坡滤波。

3c)反投影时，按照X射线源的投影方向做加权反投影，重建最终的结果就是把所有角度的射线贡献进行求和。

1)反投影时的权值计算

只需要计算0～90°的权值，其中R为射线源到旋转中心的距离，U²(x,y,θ)为射线源到物体的距离。

2)将θ，θ+90°，θ+180°和θ+270°的四个旋转对称视图中的投影图像作为一组进行排列，并将其打包放入GPU的2D纹理的四个颜色通道(红色/绿色/蓝色/alpha或RGBA)，每个通道存放一个投影图像。

3)根据纹理映射矩阵，由顶点着色器计算每一层切片四个顶点的纹理映射坐标：

U(x,y,θ)＝w(x,y,θ)；

其中射线源到探测器距离为D，射线源到旋转中心的距离为R，射线源S绕着中心轴z进行旋转，扇形角为α，锥角为β，从射线源出发过点A的射线与中心射线的夹角为λ，射线源到物体中心与虚拟坐标系x轴之间的夹角为θ。

4)在光栅化器中，矩形内部点的纹理映射坐标是由GPU自动完成插值计算。

5)在片元着色器中，根据纹理映射坐标取出相应的投影数据，完成加权累加：

其中U(x,y,θ)＝R+xcosθ+ysinθ；

其中射线源到旋转中心的距离为R，射线源S绕着中心轴z进行旋转，投影角度为θ。

6)写入帧缓存，这一步骤可以使用循环渲染到纹理(render-to-texture)技术。把四个角度同时进行渲染，M为体积切片数，每一层切片都由一个纹理内存来保存，每个纹理中都有一个单独的四个颜色通道，并且每个通道用于从一个投影视图保存反投影值，可以利用坐标FDK算法的垂直对称性同时计算1层和M层的反投影数据。

上述步骤(2)-(3)为本发明的基于GPU的快速锥束计算机断层成像重建过程。

(7)显示结果，将步骤(6)的重建结果最终进行展示，用MITK软件进行显示。

下面结合仿真对本发明的应用效果做详细的描述。

图3-图5用于仿真实验的二组数字鼠模型的三种不同视图。其中图3是截面显示，图4是冠面显示，图5是矢面表示。

图6是基于本发明的重建结果，投影角度为360个，其中(a)图显示为506×516×360的探测数据重建为512×512×512的三维展示；(b)图显示为1944×1535×360的探测数据重建为512×512×512的三维展示。(c)图显示为在小鼠二的重建数据中的第200层切片，对比其在149行的GPU与CPU所对应的灰度图。

用本发明所示的重建方法对小鼠一进行GPU重建时所需时间为9.73，前者是后者时间的160倍，当随着原始数据的增大，小鼠二在CPU和GPU下的重建时间分别是6358和27.45，速度大概提升到230倍，这也说明随着数据量的增大，GPU的处理速度也会越来越快，这说明在不改变原有的图像质量基础上有效缩小了重建的时间。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于GPU的快速锥束计算机断层成像重建方法，其特征在于，所述基于GPU的快速锥束计算机断层成像重建方法包括：获取多个角度的物体测量数据；获取计算机断层成像系统的各个参数；对测量数据做预处理；对预处理后的数据做加权滤波，基于GPU的重排滤波算法；对滤波过后的数据做反投影，基于GPU的反投影使用纹理映射理论实现并行加速，完成计算机的断层成像；

所述对预处理后的数据做加权滤波，基于GPU的重排滤波算法；对滤波过后的数据做反投影具体包括：

(1)对像素和源点之间的角度差和距离进行修正，对处理过后原始数据进行加权，使用并行算法对每个角度的测量数据同时处理：

g'(u,v,θ)＝cos(λ)g(u,v,θ)；

其中g(u,v,θ)为测量数据,

R(u,v)为射线源到旋转中心的距离，(u,v)为探测器上物体的位置；

(2)数据在频域空间进行滤波，对加权后的数据进行滤波，使用CUFFT库，并且进行数据重排，复数和虚数分别占一行数据；

g”(u,v,θ)＝g'(u,v,θ)h(u)；

其中h(u)为滤波核函数；

(3)数据进行加权反投影对滤波过后的数据按照X射线源的投影方向做加权反投影，重建最终的结果就是把所有角度的射线贡献进行求和；

所述基于GPU的快速锥束计算机断层成像重建方法进一步包括：

1)反投影时的权值计算

只需要计算0～90°的权值，其中R为射线源到旋转中心的距离，U²(x,y,θ)为射线源到物体的距离；

2)将θ，θ+90°，θ+180°和θ+270°的四个旋转对称视图中的投影图像作为一组进行排列，并打包放入GPU的2D纹理的四个颜色通道，每个通道存放一个投影图像；

3)根据纹理映射矩阵，由顶点着色器计算每一层切片四个顶点的纹理映射坐标：

U(x,y,θ)＝w(x,y,θ)；

其中射线源到探测器距离为D，射线源到旋转中心的距离为R，射线源S绕着中心轴z进行旋转，扇形角为α，锥角为β，从射线源出发过点A的射线与中心射线的夹角为λ，射线源到物体中心与虚拟坐标系x轴之间的夹角为θ；

4)在光栅化器中，矩形内部点的纹理映射坐标是由GPU自动完成插值计算；

5)在片元着色器中，根据纹理映射坐标取出相应的投影数据，完成加权累加；

其中U(x,y,θ)＝R+xcosθ+ysinθ；

其中射线源到旋转中心的距离为R，射线源S绕着中心轴z进行旋转，投影角度为θ；

6)写入帧缓存，四个角度同时进行渲染，M为体积切片数，每一层切片都由一个纹理内存来保存，每个纹理中都有一个单独的四个颜色通道，并且每个通道用于从一个投影视图保存反投影值，利用坐标FDK算法的垂直对称性同时计算1层和M层的反投影数据。

2.如权利要求1所述的基于GPU的快速锥束计算机断层成像重建方法，其特征在于，所述获取多个角度的物体测量数据利用计算机断层成像系统，依次采集测量，其中测量数据的获取为有限个角度，为绕物体一圈，隔1度获取一次，总共为360组。

3.如权利要求1所述的基于GPU的快速锥束计算机断层成像重建方法，其特征在于，所述获取计算机断层成像系统的各个参数，对测量数据做预处理具体包括：

(1)使用计算机断层成像系统来进行重建时，需要获取系统的各个参数：射线源到中心物体的距离、射线源到探测器之间的距离，探测器的分辨率、探测器的偏移度；

(2)对数据做预处理时，测量断层成像系统做没有物体时的数据，对和(1)中的测量数据做对数运算。

4.一种应用权利要求1～3任意一项所述基于GPU的快速锥束计算机断层成像重建方法的计算机断层成像系统。

Images