基于二维X射线图像序列滤波反投影的分块三维重建方法

摘要

本发明公开了一种基于二维X射线图像序列滤波反投影的分块三维重建方法，包含如下步骤：步骤一，根据初始分块数目，建立立方体分块；步骤二，对各个立方体分块进行局部重建获得分块重建结果；步骤三，组合各个分块的重建结果；步骤四，根据分块的重建结果计算所有分块数目集合，根据步骤一至步骤三进行重建并反馈其时间复杂度，直至分块数目集合中各元素均计算完毕，获得反馈的时间复杂度集合；从反馈的时间复杂度集合中选取具有最小时间复杂度的分块数目，将其作为最优分块数目，重新执行步骤一至步骤三，获得最优分块重建结果，完成三维重建。本发明能够使多线程编程技术在三维重建领域得到充分的运用，从而极大提高了重建速度。

主权项

1.一种基于二维X射线图像序列滤波反投影的分块三维重建方法，其特征在于，包含如下步骤：步骤一，根据初始分块数目，建立立方体分块：包括建立离散坐标系和确定各分块坐标区；步骤二，对各个立方体分块进行局部重建获得分块重建结果：包括待重建数据规范化、待重建数据预处理以及基于二维X射线图像的分块反投影重建；步骤三，组合各个分块的重建结果；步骤四，根据分块的重建结果计算所有分块数目集合，根据步骤一至步骤三进行重建并反馈其时间复杂度，直至分块数目集合中各元素均计算完毕，获得反馈的时间复杂度集合；从反馈的时间复杂度集合中选取具有最小时间复杂度的分块数目，将其作为最优分块数目，重新执行步骤一至步骤三，获得最优分块重建结果，完成三维重建；步骤一包括如下步骤：初始化分块数目N₀=4；建立立方体分块结构：先构建一个重建立方体，以重建立方体其中一个顶点为坐标原点，与原点相连的三边分别作为x、y、z三个坐标轴建立右手坐标系，用坐标点（x，y，z）表示任一分块中的一个体素，则任意体素坐标点(x，y，z)满足：，其中参数XDIM、YDIM、ZDIM分别表示x、y、z轴的最大值，表示自然数集合；确定各分块坐标区：保持y、z轴不变，对x轴进行分割，通过如下公式调整参数XDIM使其为分块数目N₀的整数倍：$\mathrm{XDIM}=N_0\left[\frac{1}{N_0}\mathrm{XDIM}\right];$把x轴[0，XDIM]的范围平均分割为N₀个坐标范围，即：$[0,\frac{1}{N_0}\mathrm{XDIM}-1],[\frac{1}{N_0}\mathrm{XDIM},\frac{2}{N_0}\mathrm{CDIM}-1]···[\frac{N_0-1}{N_0}\mathrm{XDIM},\mathrm{XDIM}-1];$N₀个x轴坐标范围分别结合y、z轴构成N₀个立方体分块；步骤二中基于二维X射线图像的分块反投影重建包括如下步骤：建立用于重建的二维X射线图像序列中每幅图像的像素和重建立方体中体素坐标的关系：图像序号用proj_num表示，θ表示相邻图像之间的拍摄夹角，依次把每幅图像数据存入一维序列，体素坐标为(x，y，z)，当前要进行反投影的二维X射线图像的拍摄角度为φ=proj_num·θ，所述体素在序号为proj_num的二维X射线图像中对应的坐标(x′_p，y′_p)为：$(x_p^{\prime },y_p^{\prime })=(\frac{D(y\cos \phi -x\sin \phi )}{d+x\cos \phi +y\sin \phi },\frac{D·z}{d+x\cos \phi +y\sin \phi })$其中，D为成像装置中X射线发射源到探测器的距离，d为成像装置中X射线发射源到待重建物体的几何中心的距离；令P(i′，j′)为一幅二维X射线图像中任意坐标点(i′，j′)的像素值，用差值法获得对应的经滤波后的修正像素值用三个维度的相对位置(i，j，k)表示一个立方体数据中的体素，以立方体数据的几何中心为原点，分别以i、j、k三个维度为x，y，z轴建立右手直角坐标系；对于每一个基于i轴分割开来的分块区域，其i轴范围分别为$[0,\frac{1}{N_0}\mathrm{XDIM}-1],[\frac{1}{N_0}\mathrm{XDIM},\frac{2}{N_0}\mathrm{CDIM}-1]···[\frac{N_0-1}{N_0}\mathrm{XDIM},\mathrm{XDIM}-1],$结合j、k轴表示为N₀个长方体的形式，以三元组的形式表示为：其中V表示体素值，l为分块序号；对于每一个分块内体素(i，j，k)，二维X射线图像序列在分块V_l[i][j][k]上的反投影如下式所示：$V_l\left[i\right]\left[j\right]\left[k\right]=\underset{\mathrm{proj}\_\mathrm{num}}{\overset{\mathrm{NUM}}{\Sigma }}{\hat{P}}_0{(i^{\prime },j^{\prime })}_{\mathrm{proj}\_\mathrm{num};}$其中表示在序号为proj_num的二维X射线图像中第i′行第j′个像素经滤波后的修正像素值，NUM表示二维X射线图像的最大序号，二维X射线图像总数为NUM+1，对于分块中的每个体素都用进行反投影运算，获得最后的分块重建结果步骤三中组合分块重建结果步骤如下：利用步骤二中获得的分块重建结果进行拼装组合，描述为：其中当(i，j，k)遍历立方体中所有的体素时，V[i][j][k]即为重建结果；步骤四中根据反馈的时间复杂度集合计算最优分块数目包括如下步骤：步骤41，定义综合贡献因子：分块数目N₀由x轴最大值XDIM约束，令参数α=XDIMmodN₀；定义综合贡献因子η：$\eta =\lambda ·(1-\frac{t}{t_0})+(\lambda -1)·\frac{\alpha }{\mathrm{XDIM}}$λ为表征精度和性能的参数，λ∈[0，1]；t₀和t分别表示分块数目N₀取初始值4时和取当前的数值时的重建时间；步骤42，确定分块数目：先选取所有满足参数α=0的分块，并统计其数目为m，生成集合对应的重建时间记为然后从集合中选取重建时间最少且分块数目也最少的分块数记为，对应重建时间为t_k，定义最优分块数目集合计算最优分块数目集合中各个分块数目所对应的综合贡献因子η_i：${\eta }_i=\lambda ·(1-\frac{t_{(i+N_k)}}{t_{N_k}})+(\lambda -1)·\frac{{\alpha }_{(i+N_k)}}{\mathrm{XDIM}},i=\mathrm{0,1},···,(N_{(k+1)}-N_k),$其中，和分别表示分块数目为(i+N_k)和N_k时的重建时间，${\alpha }_{(i+N_k)}=\mathrm{XDIM}\mathrm{mod}(i+N_k);$求出最优分块数目集合中每个分块数对应的综合贡献因子η_i后，取该值最小且分块数目最少的分块数为最优分块数目N₀^*。