增加了并行计算某个像素点反投影 的线程数目，但需要使用原子加操作进行写入全局存储器。
[0039] 进一步地，其中所述步骤S5的具体特征在于：
[0040] S51 :根据重建的两种分解图像，在同一内核函数中按照像素级并行进行计算等效 原子序数图像和电子密度图像。
[0041] 本发明所述一种基于CUDA架构的GPU并行加速双能谱CT解析重建方法的有益效 果在于：加快速度、提升效率。双能CT重建加速方法与采用现有单一 CT加速重建方法重复 计算两个能谱相比，在配准的前提下，在加权步骤避免了对原始投影和分解后投影重复计 算权值，在反投影步骤节省了重复计算四幅重建图像的反投影地址的时间，同时也减少了 计算中间过程的数据存取时间；在滤波步骤中又能充分利用复数的实部和虚部进行双能同 步滤波。在CT重建步骤中，反投影是最耗时的，约占总时间的85%，采用本发明的双能重建 方法，重建图像由一类增加到两类时，反投影耗时仅增加约5. 5%的时间，重建图像由两类 增加到四类时，反投影耗时仅增加约8%的时间。
【附图说明】
[0042] 图1是本发明所述双能谱CT解析重建方法的流程示意图。
【具体实施方式】
[0043] 为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例 中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是 本发明一部分实施例，而不是全部的实施例。
[0044] 如以基于基材料进行投影分解的双能扇束等角CT重建为例，阐述本发明所述基 于CUDA架构的GPU并行加速双能谱CT解析重建方法的具体实施例。
[0045] 物质的线性衰减系数的基材料分解模型如下：
[0046] y(E) =b:y! (E) +b2y2 (E) (1)
[0047] yi (E)、y2 (E)分别为两种基材料的线性衰减系数。bpb#对应两种基材料的分 解系数，对于某一固定的物质，bi、b2是两个常数。该分解公式表示任何一种物质的线性衰 减系数都可由两种基材料的线性衰减系数线性叠加而成。根据这个分解模型，基于基材料 进行投影分解的双能CT预处理重建算法的核心即为
[0049] SH(E)、SJE)分别为高低能谱，PH、P^>别为高低能投影值。B i、B2S分解后的投影 值，即V匕的线积分投影值。根据CT原理，利用滤波反投影重建算法，便可由B i、B2计算 出bpb2，由此可以计算物质的等效原子序数和电子密度信息。
[0051] Pe=b:Pel+b2Pe2 (4)
[0052] n = 3~4，ZpZ2分别为两种基材料的原子序数；P el、P e2分别为两种基材料的电 子密度。
[0053] 首先依据CT系统能谱信息生成合理物理意义范围内的各种&、B2数值对应的1\、 P H，建立投影分解查找表。投影分解时根据实际得到的1\、?11去匹配查找表中与之吻合的高 低能投影，同时获取对应的&、B 2值。
[0054] 基于参考文献"针对X射线双能CT成像的正弦图快速配准方法"，快速配准通过计 算两个正弦图fjx，y)与fjx，y)的一维投影及0〇与C(>')的最大归一化互相关系数来实 现。
[0060] A、B、C、D均为常数，所以最大的归一化互相关系数即为&(.)0.(.V)的卷积的最 大值。在GPU加速中，通过内核函数实现并行计算一维投影，并为计算FFT而进行补零操 作。假设采样数据的宽度和高度为width，height，为加快处理速度，可以仅对有效宽度为 validwidth(validwidth小于width)，高度height的采样数据进行二值化并计算一维投 影。
[0062] 将补零后的一维投影通过CUFFT函数库的FFT变换到频域cufftSafeCall(cuff tExecC2C(ProjCufftPlan, (cufftComplex*)d_lDProj, (cufftComplex*)d_lDProj,CUFFT_ FORWARD))；
[0063] 为加快频域计算时访问数据的速度，可以将一维投影绑定到纹理。
[0064] 通过内核函数实现频域并行计算点积。然后通过逆FFT变换获得卷积结果，卷积 结果的最大值所在位置，即为配准需要移动的视角数量。
[0065] 将离线生成投影分解查找表，系统初始化时装入显存并绑定到纹理。
[0066] cudaMemcpyToArray(d-LUTArrayH，0, 0，（void*)h-mLUTH，Size，cudaMemcpyHostT oDevice)；
[0067] cudaMemcpyToArray(d-LUTArrayL，0, 0，（void*)h-mLUTL，Size，cudaMemcpyHostT oDevice)；
[0068] cudaBindTextureToArray(LUTTexH,d_LUTArrayH,channelDesc)；
[0069] cudaBindTextureToArray(LUTTexL,d_LUTArrayL,channelDesc)；
[0070] 对配准后的双能采样数据，采用GPU内核函数并行计算双能投影，再根据双能投 影值通过查找表进行投影分解，并同时进行四类投影值的加权操作，加权值可以在该核函 数中直接计算，也可以事先计算好，存储在查找表中直接读取。
[0071] uintx= -umul24 (blockldx.x,blockDim.x) +threadldx.x；
[0072] uinty= -umul24 (blockldx.y,blockDim.y) +threadldx.y；
[0073] if(x<projectWidth&&y<projectHeight)
[0074] {
[0075] uintindex= -umul24(y,projectWidth)+x；
[0076] ......
[0077] d_outputH[index] =projectValueHXweight；
[0078] d_outputL[index] =projectValueLXweight；
[0079] d_decomoutputH[index] =tex2D(LUTTexH,projectValueH+0. 5,projectValueL +0. 5)Xweight；
[0080] d_decomoutputL[index] =tex2D(LUTTexL,projectValueH+0. 5,projectValueL +0. 5)Xweight；
[0081] 初始化时将滤波核补零位后用FFT变换到频域，并存储到一维纹理存储器中cu fftExecC2C(CufftPlan，（cufftComplex*)d-filter，（cufftComplex*)d-filter，CUFFT- FORWARD)
[0082] cudaBindTexture(0,FilterTex，d-filter)
[0083] 利用复数的实部和虚部同时进行高低能投影的FFT变换或者分解后的双投影的 FFT变换，处理一批一维离散傅里叶变换，并将变换后的复数值存储到复数纹理以便计算频 域点积时加快读取速度。
[0084] cufftExecC2C(ProjCufftPlan， （cufftComplex*)d-signal， （cufftComplex*)d- signal，CUFFT-FORWARD)
[0085] cudaMemcpyToArray(d-ProjArray，0, 0,d-signal，Size，cudaMemcpyDeviceToDev ice)
[0086] cudaBindTextureToArray(ProjTexH，d-ProjArray，channelDesc)
[0087] 使用内核函数对频域的双能投影和分解后投影同步进行与滤波核的点积计算。
[0088] ComplexTexMulAndScale<<<gridSize,blockSize>>> (d_signal，d_ Decomsignal,size,views,scale)；
[0089] 将滤波后的结果数据绑定到纹理存储器，以便反投影重建时加快读取速度 cudaMemcpyToArray(d_filteredProjArray, 0, 0,d_Signal,size,cudaMemcpyDeviceToDev ice)
[0090]cudaBindTextureToArray(filteredProjTex,d_ filteredProjArray,channelDesc)
[0091] 使用像素驱动的反投影方式，采用像素级并行在同一内核函数完成四种类型重建 图像的同一像素位置的全部视角views的反投影计算