基于CUDA架构的GPU并行加速双能谱CT重建方法

摘要

本发明公开一种基于CUDA架构的GPU并行加速双能谱CT解析重建方法，结合双能CT重建的特点，提供基于CUDA架构利用GPU的并行能力，进行双能CT快速解析重建同时获得四幅重建图像的方法，即：后处理重建的高低能线衰减系数图像和预处理重建的等效原子序数和电子密度图像。本发明所述方法的优点在于加快速度、提升效率，双能CT重建加速方法与现有单一CT加速重建方法相比，在配准的前提下，在加权步骤避免对原始投影和分解后投影重复计算权值，在反投影步骤节省重复计算四幅重建图像的反投影地址的时间，并减少计算中间过程的数据存取时间；在滤波步骤中又能充分利用复数的实部和虚部进行双能同步滤波；重建图像增加时，并不增加更多反投影耗时。

说明书

技术领域

本发明涉及X射线CT技术领域双能CT解析重建算法，尤其涉及一种基于 CUDA架构的GPU并行加速双能谱CT解析重建方法。

背景技术

对于CT加速领域，采用GPU加速已经成为重要的方法，而具体的加速方法 是提高CT产品的处理速度的关键环节，既是本领域技术人员研究的核心技术， 也直接决定产品的市场竞争力。

与单能CT重建不同，双能CT重建并不是前者的简单重复执行，而具有其特 殊性。双能CT解析重建算法分为后处理重建和预处理重建两种算法。其中，后 处理重建和单能CT解析重建方法一致，是对投影数据直接重建，获得被检物的 高低能线衰减系数信息，重建的形状精度较好，但由于不能从根本上消除射术硬 化效应，重建的被检物材料精度难以保证；而预处理重建与后处理算法相比，预 处理重建算法要对双能投影数据进行分解，预处理重建算法分为“基效应分解模 型”和“基材料分解模型”两种模型，两种模型是统一的，对分解后的基材料(或 基效应)投影数据进行重建，再根据基材料的重建图像计算等效原子序数和电子 密度图像。这种预处理算法理论上可以得到不受能谱影响的被检物的等效原子序 数和电子密度信息，重建的材料精度较高。双能CT重建想要同时获得上述四种 重建图像信息，时间因素是算法能否实际应用的关键。

如现有技术中，公开号为：101283913，名称：CT图像重建的GPU加速方法， 该发明专利申请公开了一种基于GPU硬件加速的CT图像重建方法，该发明软件 部分通常包括基于GPU的CT数据预处理模块、基于GPU的CT数据滤波模块、基 于GPU的CT正投影模块和基于GPU的CT图像重建反投影模块，该发明利用GPU 硬件实现了对CT图像重建算法的加速,重建部分和数据处理部分均在GPU上完 成，通过对大数据处理采用分块处理方法，对每块待重建体重建时只需要用到部 分投影数据,减少了数据传输,从而可提高整个体重建速度。

又如现有技术中，公开号为：101596113，名称：一种CT并行重建系统及成 像方法，该发明专利申请公开了一种CT并行重建系统及其成像方法,该发明所述 系统包括前端采样器、中心节点、以及与所述中心节点连接的多个子节点,每个 所述子节点均安装有图形处理器。该发明所述成像方法包括如下步骤:1)将重建 图像划分为多个BLOCK区域；2)中心节点接收前端采样器所采集的各角度下的原 始投影数据；3)中心节点为各节点分配计算任务,每个分配了计算任务的子节点 计算一个或多个BLOCK区域的重建值；4)各子节点完成所分配的重建计算；5)中 心节点组合出完整的重建图像。该发明采取划分重建图像子区域的方法充分挖掘 CT扫描及重建的并行特性,在数据采集的同时进行GPU的并行重建,重建时间由 现有技术的分钟级降低为毫秒级,达到准实时的被测对象的断层图像重建显示。

由此可知，现有技术基本上都是利用GPU并行实现单一的CT图像重建过程， 以达到加速图像重建的目的，由于CT重建本身计算量大，耗时高，算法速度一 直是CT相关算法能否在实际中使用的关键影响因素，对于双能CT解析重建算法 中预处理重建算法的特殊性，现有技术中并没有给出有针对性的完善的并行加速 解决方法，基于上述现有技术中的不足，这也成为双能CT解析重建算法发展过 程中亟待解决的问题。

发明内容

为了弥补现有技术的缺陷，本发明提供一种基于CUDA架构的GPU并行加速 双能谱CT解析重建方法，结合双能CT重建的特点，提供一种基于CUDA架构利 用GPU的并行能力，进行双能CT快速解析重建同时获得四幅重建图像的方法， 即：后处理重建的高低能线衰减系数图像和预处理重建的等效原子序数和电子密 度图像。

为了实现上述目的，本发明所述采取的技术解决方案是：本发明所述一种基 于CUDA架构的GPU并行加速双能谱CT解析重建方法，所述方法包括以下基本步 骤：

S1：高低能采样数据快速配准；

S2：并行计算原始投影数据及投影分解和投影加权；

S3：双能投影或者分解投影同步快速卷积滤波；

S4：同时并行反投影重建几类图像；

进一步的，也可以包括:

S5：由分解重建图像计算等效原子序数图像和电子密度图像。

根据不同应用，上述步骤也可以调整为S1对原始投影数据进行快速配准； S2双能投影分解和投影加权。

进一步的，其中所述步骤S1的具体特征为：

S11：所述高低能采样数据由双能谱CT系统采样获得，其双能谱为通过快速 切换射线源产生的真实双能谱或者通过使用双层夹心探测器获得的伪双能谱；

S12：所述双能谱CT系统的扫描模式为圆周扫描或者螺旋扫描；双能均采用 顺时针扫描或者逆时针扫描；

S13：所述双能谱CT系统的探测器为单排探测器或者多排探测器或者面阵探 测器；

S14：所述的CT系统相关参数在GPU设备端定义为常量参数，并在初始化时 赋值，以便加快GPU内核函数对其访问的速度；

S15：所述的配准是为保证双能CT预处理投影分解重建算法的分解精度，先 对双能采样数据进行视角配准，其特征在于：

S16：所述的双能快速配准使用计算最大互相关系数方法，为加速将计算双 能二维采样数据的互相关系数简化为对采样数据进行一维投影，再计算双能一维 投影的互相关系数；

S17：采用GPU并行方式对高低能采样进行局部探测器单元数据二值化，并 对二值化后的数据进行行方向的一维投影；

S18：通过双能一维投影的卷积实现计算互相关系数加速，使用CUFFT对双 能采样的一维投影变换到频域进行，获得互相关系数数组，寻找数组中最大互相 关系数所在位置即为配准需要调整的行数。

进一步地，其中所述步骤S2的具体特征在于：

S21：计算原始双能投影为加快读取数据的速度，将采样数据和空采明场数 据都存储在纹理存储器中；

S22：基于查找表实现双能投影分解，初始化时将预先离线生成的“投影分 解查找表”存储到显存的二维纹理存储器以便提高读取速度，该查找表纹理采用 非归一化的浮点型拾取坐标，将高低能投影值作为横纵坐标查找对应的分解投影 值，并利用纹理存储器硬件实现的线性模式滤波功能对读取的浮点型返回值进行 插值，实现快速高精度双能投影分解；

S23：在同一内核函数中按照像素级并行先计算高低能原始投影数据，再利 用高低能原始投影数据到查找表纹理中查找分解投影数据，并对原始投影和分解 投影的同时进行加权操作。根据具体应用场景不同，加权值可以离线计算后，初 始化时导入显存纹理中，在加权操作时读取，或者在该内核函数中实现同步计算 该像素点的加权值。在同一内核函数中同时实现上述计算的优点在于节省存取中 间结果的时间，加快处理速度。

进一步地，其中所述步骤S3的具体特征在于：

S31：滤波过程中使用CUDA架构提供的基于GPU的快速傅里叶变换函数库 CUFFT；

S32：初始化时将滤波核用FFT变换到频域，并存储到一维纹理存储器中， 以便在频域运算时加快读取速度；

S33：同步进行双能和双分解后投影的FFT变换，利用复数的实部和虚部同 时进行高低能投影的FFT变换或者分解后的双投影的FFT变换；

S34：利用FFT同时处理一批一维离散傅里叶变换，并将变换后的复数值存 储到复数纹理以便计算频域点积时加快读取速度；

S35：使用同一内核函数同步对双能和双分解后频域投影进行与滤波核的点 积计算，将滤波后的结果数据绑定到纹理存储器，以便反投影重建时加快读取速 度。

进一步地，其中所述步骤S4的具体特征在于：

S41：采用像素驱动的反投影重建，不存在对存储重建图像的全局显存的 “写冲突”；

S42：同一内核函数中按照像素级并行进行几类图像反投影重建，同一线程 完成四种类型重建图像的同一像素位置的反投影计算，即同时完成高低能线衰减 系数图像和两种分解图像中各自的同一像素位置的重建；

S43：重建图像存储在全局存储器中，内核函数内部使用寄存器变量临时存 储计算用的变量，计算过程中不访问全局存储器，每个内核函数完成X个均匀分 布视角的反投影计算，计算完X个视角后再给全局存储器的相应像素位置赋值， 以便节省反投影重建的访存时间，根据具体应用环境不同，该X＝N/M。N为扫描 一圈采集的全部视角；M为用来并行处理视角的线程数量，且能够被N整除，M 等于1时，对重建的像素点不存在“写冲突”，可以不使用原子加操作进行写 入全局存储器，当M大于1时，增加了并行计算某个像素点反投影的线程数目， 但需要使用原子加操作进行写入全局存储器。

进一步地，其中所述步骤S5的具体特征在于：

S51：根据重建的两种分解图像，在同一内核函数中按照像素级并行进行计 算等效原子序数图像和电子密度图像。

本发明所述一种基于CUDA架构的GPU并行加速双能谱CT解析重建方法的有 益效果在于：加快速度、提升效率。双能CT重建加速方法与采用现有单一CT 加速重建方法重复计算两个能谱相比，在配准的前提下，在加权步骤避免了对原 始投影和分解后投影重复计算权值，在反投影步骤节省了重复计算四幅重建图像 的反投影地址的时间，同时也减少了计算中间过程的数据存取时间；在滤波步骤 中又能充分利用复数的实部和虚部进行双能同步滤波。在CT重建步骤中，反投 影是最耗时的，约占总时间的85％，采用本发明的双能重建方法，重建图像由一 类增加到两类时，反投影耗时仅增加约5.5％的时间，重建图像由两类增加到四 类时，反投影耗时仅增加约8％的时间。

附图说明

图1是本发明所述双能谱CT解析重建方法的流程示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实 施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所 描述的实施例是本发明一部分实施例，而不是全部的实施例。

如以基于基材料进行投影分解的双能扇束等角CT重建为例，阐述本发明所 述基于CUDA架构的GPU并行加速双能谱CT解析重建方法的具体实施例。

物质的线性衰减系数的基材料分解模型如下：

μ(E)＝b₁μ₁(E)+b₂μ₂(E)   (1)

μ₁(E)、μ₂(E)分别为两种基材料的线性衰减系数。b₁、b₂为对应两种基材料的分 解系数，对于某一固定的物质，b₁、b₂是两个常数。该分解公式表示任何一种物 质的线性衰减系数都可由两种基材料的线性衰减系数线性叠加而成。根据这个分 解模型，基于基材料进行投影分解的双能CT预处理重建算法的核心即为

$\left(\begin{array}{c}{P}_L=-\ln \left[{\int S}_L\left(E\right)\exp \right[-B_1{\mu }_1\left(E\right)-B_2{\mu }_2\left(E\right)\left]\mathrm{dE}\right]+\ln {\int S}_L\left(E\right)\mathrm{dE}\\ P_H=-\ln \left[{\int S}_H\left(E\right)\exp \right[-B_1{\mu }_1\left(E\right)-B_2{\mu }_2\left(E\right)\left]\mathrm{dE}\right]+\ln {\int S}_H\left(E\right)\mathrm{dE}\end{array}\right)---\left(2\right)$

S_H(E)、S_L(E)分别为高低能谱，P_H、P_L分别为高低能投影值。B₁、B₂为分解后 的投影值，即b₁、b₂的线积分投影值。根据CT原理，利用滤波反投影重建算法， 便可由B₁、B₂计算出b₁、b₂，由此可以计算物质的等效原子序数和电子密度信 息。

$Z_{\mathrm{eff}}={\left[\frac{b_1{\rho }_{e1}{Z_1}^n+b_2{\rho }_{e2}{Z_2}^n}{b_1{\rho }_{e1}+b_2{\rho }_{e2}}\right]}^{\frac{1}{n}}---\left(3\right)$

ρ_e＝b₁ρ_e1+b₂ρ_e2   (4)

n＝3～4，Z₁、Z₂分别为两种基材料的原子序数；ρ_e1、ρ_e2分别为两种基材料的 电子密度。

首先依据CT系统能谱信息生成合理物理意义范围内的各种B₁、B₂数值对应的 P_L、P_H，建立投影分解查找表。投影分解时根据实际得到的P_L、P_H去匹配查找 表中与之吻合的高低能投影，同时获取对应的B₁、B₂值。

基于参考文献“针对X射线双能CT成像的正弦图快速配准方法”，快速配准通过 计算两个正弦图f₀(x,y)与f₁(x,y)的一维投影与的最大归一化互相 关系数来实现。

$\rho (P_{f_0},P_{f_1})=\frac{\underset{y=0}{\overset{n-1}{\Sigma }}{P}_{f_0}\left(y\right)P_{f_1}\left(y\right)-A-B-C}{D}---\left(5\right)$

$A=\underset{y=0}{\overset{n-1}{\Sigma }}{P}_{f_0}\left(y\right)\overline{P_{f_1}}$

$B=\underset{y=0}{\overset{n-1}{\Sigma }}{P}_{f_1}\left(y\right)\overline{P_{f_0}}$

$C=n\overline{P_{f_0}}\overline{P_{f_1}}$

$D=\sqrt{\underset{y=0}{\overset{n-1}{\Sigma }}{(P_{f_0}\left(y\right)-\overline{P_{f_0}})}^2\underset{y=0}{\overset{n-1}{\Sigma }}{(P_{f_1}\left(y\right)-\overline{P_{f_1}})}^2}$

A、B、C、D均为常数，所以最大的归一化互相关系数即为与的卷积的最大值。在GPU加速中，通过内核函数实现并行计算一维投影，并为计 算FFT而进行补零操作。假设采样数据的宽度和高度为width,height，为加快 处理速度，可以仅对有效宽度为validwidth(validwidth小于width)，高度 height的采样数据进行二值化并计算一维投影。

将补零后的一维投影通过CUFFT函数库的FFT变换到频域 cufftSafeCall(cufftExecC2C(ProjCufftPlan,(cufftComplex*)d_1DProj, (cufftComplex*)d_1DProj,CUFFT_FORWARD))；

为加快频域计算时访问数据的速度，可以将一维投影绑定到纹理。

通过内核函数实现频域并行计算点积。然后通过逆FFT变换获得卷积结果， 卷积结果的最大值所在位置，即为配准需要移动的视角数量。

将离线生成投影分解查找表，系统初始化时装入显存并绑定到纹理。

cudaMemcpyToArray(d_LUTArrayH,0,0,(void*)h_mLUTH,Size, cudaMemcpyHostToDevice)；

cudaMemcpyToArray(d_LUTArrayL,0,0,(void*)h_mLUTL,Size, cudaMemcpyHostToDevice)；

cudaBindTextureToArray(LUTTexH,d_LUTArrayH,channelDesc)；

cudaBindTextureToArray(LUTTexL,d_LUTArrayL,channelDesc)；

对配准后的双能采样数据，采用GPU内核函数并行计算双能投影，再根据双 能投影值通过查找表进行投影分解，并同时进行四类投影值的加权操作，加权值 可以在该核函数中直接计算，也可以事先计算好，存储在查找表中直接读取。

uint x＝__umul24(blockIdx.x,blockDim.x)+threadIdx.x；

uint y＝__umul24(blockIdx.y,blockDim.y)+threadIdx.y；

if(x<projectWidth&&y<projectHeight)

{

uint index＝__umul24(y,projectWidth)+x；

……

d_outputH[index]＝projectValueH×weight；

d_outputL[index]＝projectValueL×weight；

d_decomoutputH[index]＝tex2D(LUTTexH,projectValueH+0.5,projectValueL +0.5)×weight；

d_decomoutputL[index]＝tex2D(LUTTexL,projectValueH+0.5,projectValueL +0.5)×weight；

初始化时将滤波核补零位后用FFT变换到频域，并存储到一维纹理存储器中 cufftExecC2C(CufftPlan,(cufftComplex*)d_filter,(cufftComplex *)d_filter,CUFFT_FORWARD)

cudaBindTexture(0,FilterTex,d_filter)

利用复数的实部和虚部同时进行高低能投影的FFT变换或者分解后的双投 影的FFT变换，处理一批一维离散傅里叶变换，并将变换后的复数值存储到复数 纹理以便计算频域点积时加快读取速度。

cufftExecC2C(ProjCufftPlan,(cufftComplex*)d_signal,(cufftComplex *)d_signal,CUFFT_FORWARD)

cudaMemcpyToArray(d_ProjArray,0,0,d_signal, Size,cudaMemcpyDeviceToDevice)

cudaBindTextureToArray(ProjTexH,d_ProjArray,channelDesc)

使用内核函数对频域的双能投影和分解后投影同步进行与滤波核的点积计 算。

ComplexTexMulAndScale<<<gridSize,blockSize>>>(d_signal,d_Decomsignal, size,views,scale)；

将滤波后的结果数据绑定到纹理存储器，以便反投影重建时加快读取速度 cudaMemcpyToArray(d_filteredProjArray,0,0,d_Signal, size,cudaMemcpyDeviceToDevice)

cudaBindTextureToArray(filteredProjTex,d_filteredProjArray, channelDesc)

使用像素驱动的反投影方式，采用像素级并行在同一内核函数完成四种类型 重建图像的同一像素位置的全部视角views的反投影计算，假设某个视角下穿过 像素点的射线投射到探测器的水平和垂直方向的位置分别为xPos、yPos。

uint x＝__umul24(blockIdx.x,blockDim.x)+threadIdx.x；

uint y＝__umul24(blockIdx.y,blockDim.y)+threadIdx.y；

if(x<imageWidth&&y<imageHeight)

{

uint index＝__umul24(y,imageWidth)+x；

for(int i＝0；i<views；i++){

xPos＝……

yPos＝……

PixelH+＝tex2D(TexH,xPos+0.5,yPos+0.5)×scale

PixelL+＝tex2D(TexL,xPos+0.5,yPos+0.5)×scale

DecomPixelH+＝tex2D(filteredProjTexH,xPos+0.5,yPos+0.5)×scale

DecomPixelL+＝tex2D(filteredProjTexL,xPos+0.5,yPos+0.5)×scale

PixelH、PixelL、DecomPixelH、DecomPixelL即分别为四种类型重建图像的(x， y)像素的重建值。再根据公式(3)、(4)和DecomPixelH、DecomPixelL在另一核 函数中计算获得等效原子序数图像和电子密度图像。

最后应说明的是：以上实施例仅说明本发明的技术方案，而非对其限制；尽 管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理 解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技 术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本 发明各实施例技术方案的精神和范围。