CN102110310B - 利用图形处理器实现三维反投影的方法 - Google Patents

Abstract

本发明涉及一种利用图形处理器实现三维反投影的方法，包括以下步骤：将CT扫描器采集的投影数据保存在主机内存中，同时计算待重建的CT图像即目标图像中各像素点的Hermite系数及螺旋权的采样值；将主机内存中的投影数据及各像素点的Hermite系数及螺旋权的采样值拷贝到图形处理器的显存中，将显存中的数据进行纹理绑定；在GPU中进行三维反投影，得到CT重建图像；将上述得到CT重建图像拷贝到主机内存。应用本发明方法在GPU中实现三维反投影，借助GPU三维数组存储及三维纹理绑定技术，并且以float4数据类型做为二维参数数组的基本元素进行取数，避免了因参数数组过多导致GPU缓存利用率低而带来的取数慢的问题，使得运算效率进一步提高，获得了较为满意的结果。

Description

利用图形处理器实现三维反投影的方法

技术领域

本发明涉及医学图像重建技术，具体的说是一种利用图形处理器实现三维反投影的方法。

背景技术

CT图像三维重建速度一直是重建算法领域人们热议的话题，目前制约三维重建速度的主要瓶颈是其中的反投影部分，如何从软、硬件方面提高此部分的运算速度一直是人们关注的焦点。

目前，主流计算机中的处理器主要是中央处理器(CPU)和图形处理器(GPU)。传统上，GPU只负责图形渲染，而大部分的处理都交给了CPU。早在2000年就已经有人开始把部分CPU工作放在GPU上执行，以此来加速完成某些具有计算密度高、逻辑分支简单的大规模数据并行任务。但此方法受软、硬件资源的限制，无法大面积推广使用，这种情况随着NVDIA公司于2007年正式发布的CUDA(Computer Unified Deviece Architecture计算统一设备架构)发生了改变，CUDA是一种使用类C语言(支持现有C语言的基础，进行了部分扩展)进行通用计算的开发环境和软件体系，它提供了更丰富的硬件资源，同时使用类C语言的开发易于人们掌握。与CPU相比，CUDA借助GPU的卓越并行运算能力可以显著提高具有并行性特征的算法的运算速度。

美国US 2007/0014486A1专利文献中描述了使用GPU进行反投影处理的方法，该方法利用投影射线驱动的方式实现反投影。该专利利用二维反投影原理，其主要思想是：使用一个视图(View)下的所有射线数据，借助GPU纹理与顶点着色器的自动插值功能，从所采用的视图的所在角度，一次性投射到整幅待建目标图像区中，然后在下一个View所在角度再次投射其下所有通道的值到整幅目标图像区中，重复迭代该过程，直到完成所有View的投射过程。

上述专利方法不适应更高层数的医学CT图像的重建。

文章《Accelerating Backproj ections via CUDA Architecture》(HaiquanYang，Meihua Li，Kazuhito Koizumi，and Hiroyuki Kudo 2007.07.9 9thInternaltional Meeting on Fully Three-Dimensional Image Reconstruction inRadiology and Nuclear Medicine)公开了FDK三维反投影方法，是一种仅基于基本的三维反投影原理的GPU实现，图像重建速度慢。

发明内容

针对现有技术中的方法不适应更高层数的医学CT图像的重建或图像重建速度慢等不足之处，本发明要解决的技术问题是提供一种可以显著提高反投影的运算速度、适应更高层数的利用图形处理器实现三维反投影的方法。

为解决上述技术问题，本发明采用的技术方案是：

本发明一种利用图形处理器实现三维反投影的方法包括以下步骤：

将CT扫描器采集的投影数据保存在主机内存中，同时计算待重建的CT图像即目标图像中各像素点的Hermite系数及螺旋权的采样值；

将主机内存中的投影数据及各像素点的Hermite系数及螺旋权的采样值拷贝到图形处理器的显存中，将显存中的数据进行纹理绑定；

在GPU中进行三维反投影，得到CT重建图像；

将上述得到CT重建图像拷贝到主机内存。

所述计算Hermite系数的公式为：

$C^0=\frac{y_0{x}_1^2{x}_2^2}{4{\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})x_1^2{x}_2^2}{4{\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1{x}_0^2{x}_2^2}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}{x}_0^2{x}_2^2}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_2{x}_0^2{x}_1^2}{4{\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}+{3y}_2/\mathrm{\Δx})x_1^2{x}_0^2}{4{\left(\mathrm{\Δx}\right)}^4};$

$C^1=\frac{y_0(-{2x}_1{x}_2^2-{2x}_2{x}_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})(-{2x}_1{x}_2^2-{2x}_2{x}_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})x_1^2{x}_2^2}{{4\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1(-{2x}_0{x}_2^2-{2x}_2{x}_0^2)}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}(-{2x}_0{x}_2^2-{2x}_2{x}_0^2)}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1^{\′}{x}_0^2{x}_2^2}{{\left(\mathrm{\Δx}\right)}^4};$

$+\frac{y_2(-{2x}_0{x}_1^2-{2x}_1{x}_0^2)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}-{3y}_2/\mathrm{\Δx})(-{2x}_0{x}_1^2-{2x}_1{x}_0^2)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}+{3y}_2/\mathrm{\Δx})x_1^2{x}_0^2}{{4\left(\mathrm{\Δx}\right)}^4}$

$C^2=\frac{y_0(x_2^2+{4x}_1{x}_2+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})(x_2^2+{4x}_1{x}_2+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})(-{2x}_1{x}_2^2-{2x}_2{x}_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1(x_0^2+{4x}_0{x}_2+x_2^2)}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}(x_0^2+{4x}_0{x}_2+x_2^2)}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1^{\′}(-{2x}_0{x}_2^2-{2x}_2{x}_0^2)}{{\left(\mathrm{\Δx}\right)}^4};$

$+\frac{y_2(x_0^2+{4x}_0{x}_1+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}-{3y}_2/\mathrm{\Δx})(x_0^2+{4x}_1{x}_0+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}-{3y}_2/\mathrm{\Δx})(-{2x}_1{x}_0^2-{2x}_0{x}_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}$

$C^3=\frac{y_0(-{2x}_2-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})(-{2x}_2-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})(x_2^2+{4x}_1{x}_2+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1(-{2x}_0-{2x}_2)}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}(-{2x}_0-{2x}_2)}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1^{\′}(x_0^2+{4x}_0{x}_2+x_2^2)}{{\left(\mathrm{\Δx}\right)}^4};$

$=\frac{y_2(-{2x}_0-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}-{3y}_2/\mathrm{\Δx})(-{2x}_0-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}-{3y}_2/\mathrm{\Δx})(x_0^2+{4x}_1{x}_0+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}$

$C^4=\frac{y_0}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})(-{2x}_2-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1^{\′}(-2x_0-{2x}_2)}{{\left(\mathrm{\Δx}\right)}^4};$

$+\frac{y_2}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}-{3y}_2/\mathrm{\Δx})}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}-{3y}_2/\mathrm{\Δx})(-{2x}_0-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}$

$C^5=\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})}{4{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1}{{\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}-{3y}_2/\mathrm{\Δx})}{4{\left(\mathrm{\Δx}\right)}^4};$

其中，C₀～C₅为索引得到的Hermite系数，x_i为采样投影角，y_i为相应的层采样值，y’₁为层采样值差分，Δx为相邻采样的步距即投影角，i为采样序号，i＝0、1、2。

所述相应的层采样值y_i通过以下公式得到：

$y=\frac{R*(\mathrm{ReconZ}-(\mathrm{\θ}+\mathrm{\γ})*h)}{\sqrt{R^2-t^2}+(-Y\cos \mathrm{\θ}-X\sin \mathrm{\θ})}+\mathrm{MidSlice}$

其中R为CT扫描机架旋转半径，ReconZ为目标图像的z位置，θ为投影视图的角度，γ为投影视图的扇角，h为单位角度的进床长度值，t为投影视图平行束数据中旋转中心到某个射线通道的垂直距离，X为待重建图像世界坐标的横坐标，Y为待重建图像世界坐标的纵坐标，MidSlice为CT扫描器检测器的中心层值。

所述在GPU中进行三维反投影包括以下步骤：

将待重建的CT图像即目标图像按纵向、横向分解成若干个相等大小的子块；

对子块中各个目标像素点通过借助GPU并行计算能力，分别并行地采用反投影算法，即从数据集中索引视图位置、相应层以及射线通道的数据，并对上述索引到的数据与螺旋权值相乘后进行累加，所有目标图像中的像素点都完成上述累加过程后得到CT重建图像。

所述索引相应层为：

利用GPU纹理绑定功能对显存中的Hermite系数表中的各元素进行快速索引，然后根据索引得到的系数，利用秦九韶算法计算Hermite多项式的值，即层的索引值，秦九韶算法计算公式如下：

y＝((((C₅*x+C₄)*x+C₃)*x+C₂)*x+C₁)*x+C₀

其中C₀～C₅为索引得到的Hermite系数，x为某一投影视图的角度值。

所述螺旋权值为GPU显存中权值采样值的线性插值，该线性插值通过GPU纹理提供的硬件线性插值功能实现索引和计算。

所述累加时，其它计算参数以float4为单元从二维数组中取数。

本发明具有以下有益效果及优点：

1.三维图像重建速度快。应用本发明方法在GPU中实现三维反投影，其中三维反投影过程中的螺旋权值计算和层索引，可以借助GPU三维数组存储及三维纹理绑定技术，做到硬件实现计算螺旋权值时的线性插值过程和层索引的快速取数，达到加速功能；除上述螺旋权值和层索引参数数组外，其它建像参数都存储一个二维数组中，并且以float4数据类型做为二维参数数组的基本元素进行取数，减少了参数数组的声明，进而避免了因参数数组过多导致GPU缓存利用率低而带来的取数慢的问题，使得运算效率进一步提高，获得了较为满意的结果。

2.三维图像重建质量无损失。本发明利用Hermite多项式进行层索引的逼近，从而避免了精确但复杂的层索引计算。

附图说明

图1为本发明方法流程图；

图2为本发明方法中重建图像某点对应的所有投影视图的层索引趋势示意图；

图3为本发明方法中目标图像及其中一个图像子块示意图；

图4为本发明方法中某个图像子块对应的所有投影视图的通道索引块示意图。

具体实施方式

如图1所示，本发明利用图形处理器实现三维反投影的方法包括以下步骤：

将CT扫描器采集的投影数据保存在主机内存中，同时计算待重建的CT图像即目标图像中各像素点的Hermite系数及螺旋权的采样值；

将主机内存中的投影数据及各像素点的Hermite系数及螺旋权的采样值拷贝到图形处理器的显存中，将显存中的数据进行纹理绑定；

在GPU中进行三维反投影，得到CT重建图像；

将上述得到CT重建图像拷贝到主机内存。

计算hermite系数的过程如下：

对待重建图像点(X，Y)，起始结束投影分别为VS_X，Y和VE_X，Y：对于k＝0，1，…N-1，执行如下步骤：

A.计算投影视图的采样值：

θ_k＝(VS_XY+k·step_X，Y)·AglPerView

θ′_k＝(VS_XY+(k+1)·step_X，Y)·AglPerView

B.计算投影视图在θ_k处的层

$y\left({\mathrm{\θ}}_k\right)=\frac{R*(\mathrm{ReconZ}-({\mathrm{\θ}}_k+\mathrm{\γ})*h)}{\sqrt{R^2-t^2}+(-Y\cos {\mathrm{\θ}}_k-X\sin {\mathrm{\θ}}_k)}+\mathrm{MidSlice}$

C.计算投影视图在θ′_k处的层

$y\left({\mathrm{\θ}}_k^{\′}\right)=\frac{R*(\mathrm{ReconZ}-({\mathrm{\θ}}_k^{\′}+\mathrm{\γ})*h)}{\sqrt{R^2-t^2}+(-Y\cos {\mathrm{\θ}}_k-X\sin {\mathrm{\θ}}_k)}+\mathrm{MidSlice}$

D.计算差分

DY(θ_k)＝Y(θ′_k)-Y(θ_k)

上述各式中，AglPerView为每个投影的角度值，MidSlice为CT扫描器检测器的中心层值，R为CT扫描机架旋转半径，ReconZ为目标图像的z位置，θ为投影视图的角度，γ为投影视图的扇角，h为单位角度的进床长度值，t为投影视图平行束数据中旋转中心到某个射线通道的垂直距离，X为待重建图像世界坐标的横坐标，Y为待重建图像世界坐标的纵坐标，k为投影视图序号，k＝0，1，…N-1，N为总投影视图的个数，step_X，Y为投影视图采样间距。

根据计算出的各投影视图下的层y(θ_k)和层差分Dy(θ_k)，代入计算Hermite系数的公式进行计算，得到Hermite系数C_k ⁰，C_k ¹，…，C_k ⁵。

所述计算Hermite系数的公式为：

$C^0=\frac{y_0{x}_1^2{x}_2^2}{4{\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})x_1^2{x}_2^2}{4{\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1{x}_0^2{x}_2^2}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}{x}_0^2{x}_2^2}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_2{x}_0^2{x}_1^2}{4{\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}+{3y}_2/\mathrm{\Δx})x_1^2{x}_0^2}{4{\left(\mathrm{\Δx}\right)}^4};$

$C^1=\frac{y_0(-{2x}_1{x}_2^2-{2x}_2{x}_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})(-{2x}_1{x}_2^2-{2x}_2{x}_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})x_1^2{x}_2^2}{{4\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1(-{2x}_0{x}_2^2-{2x}_2{x}_0^2)}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}(-{2x}_0{x}_2^2-{2x}_2{x}_0^2)}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1^{\′}{x}_0^2{x}_2^2}{{\left(\mathrm{\Δx}\right)}^4};$

$+\frac{y_2(-{2x}_0{x}_1^2-{2x}_1{x}_0^2)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}-{3y}_2/\mathrm{\Δx})(-{2x}_0{x}_1^2-{2x}_1{x}_0^2)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}+{3y}_2/\mathrm{\Δx})x_1^2{x}_0^2}{{4\left(\mathrm{\Δx}\right)}^4}$

$C^2=\frac{y_0(x_2^2+{4x}_1{x}_2+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})(x_2^2+{4x}_1{x}_2+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})(-{2x}_1{x}_2^2-{2x}_2{x}_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1(x_0^2+{4x}_0{x}_2+x_2^2)}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}(x_0^2+{4x}_0{x}_2+x_2^2)}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1^{\′}(-{2x}_0{x}_2^2-{2x}_2{x}_0^2)}{{\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_2(x_0^2+{4x}_0{x}_1+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}-{3y}_2/\mathrm{\Δx})(x_0^2+{4x}_1{x}_0+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}-{3y}_2/\mathrm{\Δx})(-{2x}_1{x}_0^2-{2x}_0{x}_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}$

；

$C^3=\frac{y_0(-{2x}_2-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})(-{2x}_2-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})(x_2^2+{4x}_1{x}_2+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1(-{2x}_0-{2x}_2)}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}(-{2x}_0-{2x}_2)}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1^{\′}(x_0^2+{4x}_0{x}_2+x_2^2)}{{\left(\mathrm{\Δx}\right)}^4};$

$=\frac{y_2(-{2x}_0-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}-{3y}_2/\mathrm{\Δx})(-{2x}_0-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}-{3y}_2/\mathrm{\Δx})(x_0^2+{4x}_1{x}_0+x_1^2)}{{4\left(\mathrm{\Δx}\right)}^4}$

$C^4=\frac{y_0}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_0(y_0^{\′}+{3y}_0/\mathrm{\Δx})}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})(-{2x}_2-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}$

$+\frac{y_1}{{\left(\mathrm{\Δx}\right)}^4}-\frac{x_1{y}_1^{\′}}{{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1^{\′}(-2x_0-{2x}_2)}{{\left(\mathrm{\Δx}\right)}^4};$

$+\frac{y_2}{{4\left(\mathrm{\Δx}\right)}^4}-\frac{x_2(y_2^{\′}-{3y}_2/\mathrm{\Δx})}{{4\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}-{3y}_2/\mathrm{\Δx})(-{2x}_0-{2x}_1)}{{4\left(\mathrm{\Δx}\right)}^4}$

$C^5=\frac{(y_0^{\′}+{3y}_0/\mathrm{\Δx})}{4{\left(\mathrm{\Δx}\right)}^4}+\frac{y_1}{{\left(\mathrm{\Δx}\right)}^4}+\frac{(y_2^{\′}-{3y}_2/\mathrm{\Δx})}{4{\left(\mathrm{\Δx}\right)}^4};$

其中，C₀～C₅为索引得到的Hermite系数，x_i为采样投影角，y_i为相应的层采样值，y′_i为层采样值差分，Δx为相邻采样的步距即投影角，i为采样序号，i＝0、1、2。

在反投影过程中，根据(X，Y)点计算出来的hermite系数，计算多项式的值即可，过程如下：

对i＝0，1，…VE_X，Y-VS_X，Y，执行如下过程：

A.计算系数索引：

B.计算多项式的值：

$y\left({\mathrm{\θ}}_i\right)=C_k^0+C_k^1{\mathrm{\θ}}_i+C_k^2{\left({\mathrm{\θ}}_i\right)}^2+C_k^3{\left({\mathrm{\θ}}_i\right)}^3+c_k^4{\left({\mathrm{\θ}}_i\right)}^4+C_k^5{\left({\mathrm{\θ}}_i\right)}^5$

其中B步多项式的计算可以用秦九韶算法迭代优化为复杂度O(n)，n为多项式次数。

如图2所示，该曲线为待重建图像中的某一点随着投影视图的变化，其选择的层的变化近似规律，可以用Hermite多项式来近似逼近。

所述在GPU中进行三维反投影包括以下步骤：

如图3所示，将待重建的CT图像即目标图像按纵向、横向分解成若干个相等大小的子块；

如图4所示，对子块中各个目标像素点分别并行地采用反投影算法，即从数据集中索引相应层、射线通道以及视图位置的数据，并对上述索引到的数据与螺旋权值相乘后进行累加，得到CT重建图像。

由于三维数组的元素个数比较多，且GPU缓存相对较小，在累加取数时，容易造成取数慢的问题，因此本实施例采用float4做为二维维数组的基本元素进行取数，进而提高了数组的读取速度，使得运算效率进一步提高，获得了较为满意的结果。

所述索引相应层为：

利用GPU纹理绑定功能对显存中的Hermite系数表中的各元素进行快速索引，然后根据索引得到的系数，利用秦九韶算法计算Hermite多项式的值，即层的索引值，秦九韶算法计算公式如下：

y＝((((C₅*x+C₄)*x+C₃)*x+C₂)*x+C₁)*x+C₀

其中C0～C5为索引得到的Hermite系数，x为某一投影视图的角度值。

在三维反投影过程中，不仅需要对投影射线通道进行索引，而且也需要对这些投影射线通道所在层进行索引，随着投影视图的不同，投影射线通道所在的层也是发生变化的，即重建中的某一像素值，它需要对穿过该点的各个view下某一层中某个投影射线通道值累加。关于层的计算过程如下：

(1)利用GPU纹理绑定功能实现对Hermite系数表的快速索引。然后根据这些系数，利用计算多项式值的秦九韶算法进行计算，过程如下：

秦九韶算法：

把一个n次多项式

f(x)＝a[n]x^n+a[n-1]x^(n-1)+......+a[1]x+a[0]

改写成如下形式

f(x)＝(......((a[n]x+a[n-1])x+a[n-2])x+......+a[1])x+a[0]；

对螺旋权值的逼近可以简单在主机中对各象素点用到的view螺旋权值进行采样，然后存储下来，在GPU的反投影过程中作硬件上的线性插值即可；由于有效减少了乘法运算次数，使得运行效率明显提升。

Claims (5)

1.一种利用图形处理器实现三维反投影的方法，其特征在于包括以下步骤：

将CT扫描器采集的投影数据保存在主机内存中，同时计算待重建的CT图像即目标图像中各像素点的Hermite系数及螺旋权的采样值；将主机内存中的投影数据及各像素点的Hermite系数及螺旋权的采样值拷贝到图形处理器的显存中，将显存中的数据进行纹理绑定；在GPU中进行三维反投影，得到CT重建图像；将上述得到CT重建图像拷贝到主机内存；

所述在GPU中进行三维反投影包括以下步骤：

将待重建的CT图像即目标图像按纵向、横向分解成若干个相等大小的子块；

对子块中各个目标像素点通过借助GPU并行计算能力，分别并行地采用反投影算法，即从数据集中索引视图位置、相应层以及射线通道的数据，并对上述索引到的数据与螺旋权值相乘后进行累加，所有目标图像中的像素点都完成上述累加过程后得到CT重建图像；

所述索引相应层为：

利用GPU纹理绑定功能对显存中的Hermite系数表中的各元素进行快速索引，然后根据索引得到的系数，利用秦九韶算法计算Hermite多项式的值，即层的索引值，秦九韶算法计算公式如下：

y＝((((C₅*x+C₄)*x+C₃)*x+C₂)*x+C₁)*x+C₀

其中C₀～C₅为索引得到的Hermite系数，x为某一投影视图的角度值。

2.按权利要求1所述的利用图形处理器实现三维反投影的方法，其特征在于所述计算Hermite系数的公式为：

；

；

；

；

其中，C₀～C₅为索引得到的Hermite系数，x_i为采样投影角，y_i为相应的层采样值，y′_i为层采样值差分，Δx为相邻采样的步距即投影角，i为采样序号，i＝0、1、2。

3.按权利要求2所述的利用图形处理器实现三维反投影的方法，其特征在于：所述相应的层采样值y_i通过以下公式得到：

其中R为CT扫描机架旋转半径，ReconZ为目标图像的z位置，θ为投影视图的角度，γ为投影视图的扇角，h为单位角度的进床长度值，t为投影视图平行束数据中旋转中心到某个射线通道的垂直距离，X为待重建图像世界坐标的横坐标，Y为待重建图像世界坐标的纵坐标，MidSlice为CT扫描器检测器的中心层值。

4.按权利要求1所述的利用图形处理器实现三维反投影的方法，其特征在于所述螺旋权值为GPU显存中螺旋权的采样值的线性插值，该线性插值通过GPU纹理提供的硬件线性插值功能实现索引和计算。

5.按权利要求1所述的利用图形处理器实现三维反投影的方法，其特征在于：所述累加时，除上述螺旋权值和层索引参数数组外，其他建像参数都存储于一个二维数组中，并且以float4数据类型作为二维参数数组的基本元素进行取数。 

Images