CN105809723A - Cbct重建方法及系统 - Google Patents

Abstract

本发明提供了一种CBCT重建方法，包括如下步骤：在全局坐标系的Z轴方向上，设定相邻两个重建切片之间的间隔等于射线源的相邻两个位置点之间的间隔的非负整数倍；利用射线源获取K帧投影数据；若投影数据对当前重建切片为部分贡献时，根据当前投影数据对当前重建切片的贡献范围边界建立临界数组，根据临界数组和当前投影数据针对每帧投影数据进行反投影运算。本发明还提供了一种CBCT重建系统，包括预设模块、数据获取模块、反投影重建模块。本发明还提供了一种CBCT系统，包括射线源、探测器、上位机和FPGA加速电路。本发明的CBCT重建方法及系统，简化了系统的计算量，提高了反投影运算的效率，降低了系统的复杂程度，易于FPGA实现。

Description

CBCT重建方法及系统

技术领域

本发明涉及成像技术领域，特别是涉及一种CBCT重建方法及系统。

背景技术

CBCT(conebeamcomputedtomography)即锥束CT，其射线束为锥形，相比于传统的平行束、扇形束CT(ComputedTomography)，具有扫描速度快、图像分辨率高、辐射利用率高、各向同性空间分辨力等优点，也因此广泛用在口腔医疗和放射治疗等领域。

CBCT重建算法按精确程度可分为近似重建算法和精确重建算法。其中近似重建算法以FDK(Feldkamp-Davis-Kress)算法为代表，该重建算法的轨迹是圆形，数学形式上简单，实现起来容易，在锥角小于5°的情况下，能够很好的抑制伪影。但随着角度的增加，FDK算法中的加权相对误差越来越大，由此引入的噪声将导致重建图像的质量下降，且扫描轨迹只适合圆形。

精确重建算法以Katsevich算法为代表，Katsevich能够兼顾重建质量和重建速度，其分为求导、长度校正、前向插值、DHT滤波、后向插值和反投影六大步，反投影部分是耗时比较多且与重建质量密切相关的步骤，且传统上又分为基于PI线(PI-line)和基于Tam窗(Tam-DanielssonWindow)两种反投影方法。

基于PI线的原理是每个待重建像素点对应着一条PI线，对每个重建点在PI线区间内做投影数据积分运算。当投影数据过来时，需要逐个像素点判断此投影帧的位置是否在该像素点的PI区间内，若是，则进行反投影累加，若否则略过。但求解PI线区间需要求解非线性方程，运算量大，且每个像素点都需要开辟两个数据空间用于存储PI线区间的两个端点，大大增加了内存，电路结构较复杂且成本高。

基于Tam窗需要先计算重建像素点的在当前射线源照射下的投影坐标，然后判断该投影坐标是否在Tam窗内，若是，则取投影数据进行累加，若否则略过该射线源处的投影。并且在多切片重建时，由于重建像素点规模大，在射线源的位置更新时需要重新计算所有像素点是否在更新后的射线源锥束覆盖的范围内，这样就存在很多无效的投影点坐标计算，计算量较大。

此外，由于数据量大的问题，单纯的用一个CPU(CentralProcessingUnit)计算重建速度慢，效率低。目前多采用计算机集群、GPU(GraphicsProcessingUnit)、FPGA(FieldProgrammableGateArray)进行加速。但是，采用计算机集群对重建算法进行加速，成本高且难以突破通用处理器自身的瓶颈，不适合临床应用。GPU虽然可以通过大量运算单元(流处理器)对数据进行并行加速处理，对算法并行化加速相当有效，但GPU对于本身不能高线程并行化的工作带来的帮助不大。同时，针对Katsevich算法的FPGA实现，都是在现有的PI线和Tam窗理论上进行，没有提出新的方案。

发明内容

鉴于现有技术的现状，本发明的目的在于提供一种CBCT重建方法及系统，减小了系统的计算量，降低了系统的复杂程度，易于通过FPGA实现。

为实现上述目的，本发明采用如下技术方案：

一种CBCT重建方法，包括如下步骤：

在全局坐标系的Z轴方向上，设定相邻两个重建切片之间的间隔等于射线源相邻两个位置点之间的间隔的非负整数倍；

利用所述射线源获取K帧投影数据，其中，K为正整数；

若所述投影数据对当前重建切片为部分贡献时，针对每帧所述投影数据，按照以下步骤进行反投影过程：

获取当前投影数据对所述重建切片的贡献范围边界，并根据所述贡献范围边界建立临界数组；

基于射线源相邻两个位置点之间的间隔建立相邻两个重建切片之间的贡献范围边界的坐标旋转矩阵；

利用所述贡献范围边界的坐标旋转矩阵，由已重建切片的临界数组迭代获得与所述已重建切片相邻的所述重建切片的临界数组；

根据所述重建切片的临界数组和当前投影数据，将当前投影数据反投影累加至上一次的重建数据中，生成更新的重建数据；

根据所述更新的重建数据，获得所述重建切片的重建图像。

在本发明的一个实施例中，每帧所述投影数据对应获得一个长度为M的临界数组；其中，M为重建切片的行分辨率，所述临界数据的索引表示所述当前重建切片的行，所述临界数组内的元素值表示相应行的边界点列坐标。

在本发明的一个实施例中，所述反投影过程还包括：

获取当前投影数据的帧序号；

根据所述当前投影数据的帧序号，判断所述当前投影数据对所述重建切片的贡献度，生成用以表示所述当前投影数据类别的贡献系数，其中，所述投影数据类别包括无贡献数据、部分贡献数据以及完全贡献数据；

若所述贡献系数表示所述投影数据为无贡献数据，则略过所述当前投影数据；

若所述贡献系数表示所述投影数据为完全贡献数据，则逐行对所述重建切片的所有重建点进行反投影累加；

若所述贡献系数表示所述投影数据为部分贡献数据，则获取每帧所述投影数据对所述重建切片的贡献范围边界，根据所述贡献范围边界逐行对所述重建切片的重建点进行反投影累加。

在本发明的一个实施例中，若所述贡献系数表示当前投影数据为完全贡献数据，则逐行对所述重建点进行反投影累加的步骤包括：

针对每一个重建切片，逐行获取所述重建切片的所有重建点的上一次的重建数据；

针对每一个所述重建点，计算反投影系数；

根据所述反投影系数，将所述当前投影数据累加至所述上一次的重建数据中，生成更新的重建数据。

在本发明的一个实施例中，若所述贡献系数表示为部分贡献数据，则根据所述贡献范围边界逐行对所述重建切片的重建点进行反投影累加的步骤还包括：

获取当前重建切片的序号以及所述当前重建切片的当前重建行的行号；

根据所述当前投影数据的帧序号、所述当前重建切片的序号以及所述当前重建行的行号生成所述临界数组的查找地址；

根据所述查找地址从所述临界数组中获得与所述当前重建行的行号对应的边界值；

判断重建切片的个数，若所述重建切片的数量为一个，则根据所述边界值逐行判断所述重建切片的当前重建行是否需要所述当前投影数据，若是，则获取所述当前重建行的上一次的重建数据，将所述当前投影数据累加至所述上一次的重建数据中，生成所述更新的重建数据；若否，则直接跳过所述当前重建行，进行下一行的反投影过程，直至遍历所述当前重建切片的所有行之后，继续下一帧投影数据的反投影过程；

若所述重建切片的数量为两个以上时，则执行如下步骤：

根据所述贡献范围边界的坐标旋转矩阵，针对所述边界值进行坐标旋转变换，获得转换边界值；

根据所述转换边界值逐行判断所述当前重建切片的当前重建行是否需要所述当前投影数据，若是，则获取所述当前重建行的上一次的重建数据，将所述当前投影数据累加至所述上一次的重建数据中，生成所述更新的重建数据；若否，则直接跳过所述当前重建行，进行下一行的反投影过程，直至遍历所述当前重建切片的所有行之后继续下一重建切片的重建，当完成当前投影数据对所有重建切片的反投影运算后，继续下一帧投影数据的反投影过程。

在本发明的一个实施例中，所述反投影过程还包括如下步骤：

在所述重建切片确定的范围内遍历所有重建点；

根据所述参考重建点的坐标计算重建所述重建切片所需射线源的最大位置和最小位置；

根据所述所需射线源的最大位置和最小位置判断是否需要略过所述当前投影数据。

在本发明的一个实施例中，基于所述射线源的相邻两个位置点之间的间隔建立相邻两个重建切片之间的贡献范围边界的坐标旋转矩阵的步骤包括：

获取所述K帧投影数据对应的所述射线源的K个位置点；

根据射线源的K个位置点获得相邻两个位置点之间的角度差Δs，其中，Δs＝2π/(K-1)；

根据转换关系式A＝YB获得相邻的下一个重建切片的转换边界值，其中，A表示所述射线源在第n(n＝1,2……K)个位置点处产生的锥束覆盖对第i个重建切片的贡献范围边界上的一个点的xy平面坐标，B表示所述射线源在第n+1个位置点处产生的锥束覆盖对第i+1个重建切片的贡献范围边界上的一个点的xy平面坐标，Y表示坐标转换矩阵，其中，

$Y=\left\{\begin{array}{cc}cos\frac{2\mathrm{\π}}{(K-1)}& -sin\frac{2\mathrm{\π}}{(K-1)}\\ sin\frac{2\mathrm{\π}}{(K-1)}& \cos \frac{2\mathrm{\π}}{(K-1)}\end{array}\right\}.$

在本发明的一个实施例中，所述反投影过程之前还包括对每帧所述投影数据依次进行求导、长度校正、前向插值、滤波、后向插值的过程。

基于重建方法，本发明还提供了一种基于同一发明构思的CBCT重建系统，包括：

预设模块，在全局坐标系的Z轴方向上，设定相邻两个重建切片之间的间隔等于射线源的相邻两个位置点之间的间隔的非负整数倍；

数据获取模块，用于利用所述射线源获取K帧投影数据，其中，K为正整数；

反投影重建模块，用于当投影数据对当前重建切片为部分贡献时，针对每帧所述投影数据，进行反投影过程，包括：

边界生成单元，用于获取当前投影数据对所述重建切片的贡献范围边界，并根据所述贡献范围边界建立临界数组；

转换矩阵单元，用于基于所述射线源的相邻两个位置点之间的间隔建立相邻两个重建切片之间的贡献范围边界的坐标旋转矩阵；

旋转单元，用于利用所述贡献范围边界的坐标旋转矩阵，由已重建切片的临界数组迭代获得与所述已重建切片相邻的所述重建切片的临界数组；

反投影单元，用于根据所述重建切片的临界数组和当前投影数据，将当前投影数据反投影累加至上一次的重建数据中，生成更新的重建数据；以及

重建单元，用于根据所述更新的重建数据，获得所述重建切片的重建图像。

此外，本发明还提供了一种CBCT重建系统，包括：

射线源，用于发射锥束射线；

与所述射线源相对设置的探测器，用于接收锥束射线穿过被测物体后的射线，并将光信号转化为数字信号，获取K帧投影数据，其中，K为正整数；

FPGA加速电路，用于对K帧投影数据进行反投影处理，生成更新的重建数据；以及

上位机，用于根据所述更新的重建数据，获得所述重建切片的重建图像并显示，其中，所述FPGA加速电路包括：

存储模块，用于保存反投影过程中所述重建切片的重建数据，以及设定的相邻两个重建切片之间的间隔与所述射线源的相邻两个位置点之间的间隔的倍数关系；

反投影模块，用于基于贡献范围边界的坐标旋转矩阵，进行反投影叠加生成更新的重建数据，所述反投影模块包括：

边界查询单元，用于获取当前投影数据对所述重建切片的贡献范围边界，并根据所述贡献范围边界建立临界数组，利用所述贡献范围边界的坐标旋转矩阵，由已重建切片的临界数组迭代获得与所述已重建切片相邻的所述重建切片的临界数组；以及

反投影运算单元，用于根据所述重建切片的临界数组和所述当前投影数据，将当前投影数据反投影累加至上一次的重建数据中，生成所述更新的重建数据。

本发明的有益效果是：

本发明的CBCT重建方法及系统，先判断当前投影数据对当前重建切片的贡献度，继而通过贡献范围边界来筛选重建像素点进行反投影，从而简化了系统的计算量，并通过将系统参数设定为相邻两个重建切片之间的间隔等于相邻射线源之间的间隔的非负整数倍，使得在进行多切片重建时，能够利用坐标旋转矩阵，由已重建切片的临界数组迭代获得与已重建切片相邻的重建切片的临界数组，从而简化了多切片重建时的计算过程，提高了反投影运算的效率，降低了系统的复杂程度，易于FPGA实现。

附图说明

图1为本发明的CBCT重建方法一实施例的流程图；

图2为本发明的CBCT重建方法中反投影运算一实施例的流程图；

图3为本发明的CBCT重建系统一实施例的原理框图；

图4为本发明的CBCT重建系统另一实施例的原理框图；

图5为本发明的CBCT重建系统中单切片重建一实施例的原理框图；

图6为本发明的CBCT重建系统中多切片重建一实施例的原理框图；

图7为本发明中投影数据的贡献度以及投影数据对当前重建切片的贡献范围边界的示意图。

具体实施方式

为了使本发明的技术方案更加清楚，以下结合附图，对本发明的CBCT重建方法及系统作进一步详细的说明。应当理解，此处所描述的具体实施例仅用以解释本发明并不用于限定本发明。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

图1所示的CBCT重建方法，可以应用于一包括射线源(未示出)、探测器(未示出)、上位机和FPGA加速电路的CBCT系统中，其中，射线源与探测器相对设置，上位机和FPGA加速电路通过PCI-E总线连接，如图4所示。

如图1所示，本发明提供了一种CBCT重建方法，可应用于单切片的重建以及多切片的重建，包括如下步骤：

S100、设定系统参数，即在全局坐标系的Z轴方向上，设定相邻两个重建切片之间的间隔等于射线源的相邻两个位置点之间的间隔的非负整数倍。由于在CT扫描成像的过程中，射线源做螺旋轨迹运动，射线源每运动一定的弧度，就停下来照射一次，获得一帧投影数据。假设射线源在的2π旋转区间，停下来K次，就产生了K个离散的位置点，获得K帧投影数据。

若相邻两个重建切片在Z轴方向上的间隔为Δx₃，相邻两个射线源之间的间隔为h_s，则Δx₃＝m*h_s，其中m可以取0、1、2等非负整数。

本实施例中，设定重建切片的数量为i，其中i取正整数。当m＝0时，说明重建切片的数量为单个，即重建切片的数量i＝1。当m为正整数时，说明重建切片的数量i≥2，重建切片为多切片。通过设定系统参数，使得相邻两个射线源对相邻两个重建切片的贡献度之间存在一个对应关系，根据该对应关系可以根据重建切片i的相关数据获得重建切片i+1的相关数据，从而可以根据该对应关系进行迭代，加速多个重建切片的重建过程，节省了计算射线源对所有重建切片贡献度的时间，减少了系统的运算量。

例如，当重建切片的数量为两个以上时，射线源的相邻两个位置点分别标记为s_n和s_n+1，则射线源的相邻两个位置点s_n和s_n+1之间的角度差为Δs(单位为弧度)，。假设射线源螺旋轨迹的螺距为h，在射线源的2π旋转区间内取K个离散的位置点，相应的产生K帧投影数据，相邻两个射线源之间在Z轴上的高度差h_s，即：h_s＝hΔs/2π，K-1＝2π/Δs。若相邻两个重建切片的Z轴距离Δx₃取h_s的整数倍，则相邻两个重建切片之间的间隔与相邻两个射线源之间的间隔存在一个很好的对应关系。

此处取特殊值假定两者相等，即h_s＝Δx₃，则射线源的相邻两个位置点之间的高度增加h_s，相邻两个重建切片的坐标值增加Δx₃，则s_n对重建切片i的贡献度和s_n+1对重建切片i+1的贡献度一样(s_n、s_n+1分别表示第n个和第n+1个位置点)，因为s_n和重建切片i、s_n+1和重建切片i+1的相对位置是一样的，则重建切片i+1时可以利用重建切片i的相关数据，以此类推，多切片重建时重建后边的切片可以用到前面的数据，节省了计算射线源对所有切片贡献程度所需要的时间。

S200、利用射线源获取K帧投影数据，其中，K为正整数。本实施例中，在单切片重建时，射线源以该重建切片为中心上下各π弧度形成的2π弧度范围内产生K帧投影数据，根据上述K帧投影数据即可获得该重建切片的重建图像。当重建切片两个以上的单切片的形成的多切片时，射线源在2π的螺旋区间内选取K个离散的位置点，产生K帧投影数据，其中，射线源的K个离散的位置点与K帧投影数据一一对应设置。

S300、针对每帧投影数据依次进行求导、长度校正、前向插值、滤波、后向插值的过程；本实施例中，上位机用于对投影数据进行求导、长度校正、前向插值等预处理，并将预处理后的投影数据通过PCI-E总线传送至FPGA加速电路。FPGA加速电路中设置有PCI-E接口模块，通过PCI-E接口模块接收上位机预处理后的投影数据，之后对预处理后的投影数据进行滤波、后向插值以及反投影运算，以获得重建数据。当然，在其他实施例中，也可以通过上位机对投影数据进行求导、长度校正、前向插值、滤波、后向插值等预处理，通过FPGA加速电路对预处理后的投影数据进行反投影运算。

若所述投影数据对当前重建切片为部分贡献时，针对每帧投影数据，按照以下步骤进行反投影过程：

S410、获取当前投影数据对当前重建切片的贡献范围边界，并根据贡献范围边界建立临界数组；其中，每帧投影数据对应获得一个长度为M的临界数组，其中，M为重建切片的行分辨率；所述临界数组的索引表示当前重建切片的行，所述临界数组内的元素值表示相应行的边界点列坐标。此处的临界数组的索引可以理解为临界数组的查找地址。

如图7所示，圆形区域表示重建切片的感兴趣区域，曲线部分表示贡献范围边界，该贡献范围边界左侧的区域内的重建点在射线源的锥束覆盖范围内，即贡献边界左侧的区域内的重建点需要该帧投影数据。该贡献范围边界右侧的区域内的重建点不在射线源的锥束覆盖范围内，即贡献范围边界右侧的区域内的重建点不需要该帧投影数据。

具体地，每帧投影数据对重建切片的贡献范围边界的推导过程如下：

任意射线源y(s)在参考坐标系中的坐标表示为：

$\left\{\begin{array}{c}u=-Rcos\left(s\right)sin\left(s_0\right)+Rsin\left(s\right)cos\left(s_0\right)\\ v=\frac{h}{2\mathrm{\π}}s-\frac{h}{2\mathrm{\π}}{s}_0\\ w=R\cos \left(s\right)\cos \left(s_0\right)+R\sin \left(s\right)\sin \left(s_0\right)+R\end{array}\right.---\left(1\right)$

通过射线源y(s)与射线源y(s0)的空间直线的参数方程表示为：

$\left\{\begin{array}{c}u=Rt(-cos(s\left)sin\right(s_0)+sin(s\left)cos\right(s_0\left)\right)\\ v=\frac{ht}{2\mathrm{\π}}(s-s_0)\\ w=Rt\left(cos\right(s\left)\cos \right(s_0)+sin(s\left)sin\right(s_0)-1)+2R\end{array}\right.---\left(2\right)$

当前重建切片在全局坐标系中的平面方程z在参考坐标系中的参数方程表示：

$v=z-\frac{h}{2\mathrm{\π}}{s}_0---\left(3\right)$

则联立方程(2)、(3)可以获得上述空间直线与当前重建切片的交点的集合，具体表示为：

$\left\{\begin{array}{c}u=\frac{R(2\mathrm{\π}z-{\mathrm{hs}}_0)}{h(s-s_0)}sin(s-s_0)\\ w=\frac{R(2\mathrm{\π}z-{\mathrm{hs}}_0)}{h(s-s_0)}\left(cos\right(s-s_0)-1)+2R\end{array}\right.---\left(4\right)$

其中，R表示射线源的螺旋线轨迹半径，s表示参数方程的参变量，s0表示当前投影数据所对应的射线源的位置点，h表示螺旋轨迹的螺距。

S420、基于所述射线源的相邻两个位置点之间的间隔建立相邻两个重建切片之间的贡献范围边界的坐标旋转矩阵；

S430、利用贡献范围边界的坐标旋转矩阵，由已重建切片的临界数组迭代获得与该已重建切片相邻的重建切片的临界数组；本实施例中，当重建切片的数量为1时，则坐标旋转矩阵中各个元素值可以均取1，以使得通过坐标旋转后的坐标值保持不变。当然，为了简化系统的运算量，在重建切片的数量为1时，可以不设置坐标转换矩阵。

当重建切片的数量为一个以上时，由于s_n对切片i的贡献度和s_n+1对切片i+1的贡献度一样，此处的贡献度包括投影范围的大小与形状，因此，切片i的边界值与切片i+1的边界值之间存在一个旋转关系，通过坐标旋转矩阵对切片i的边界值进行坐标旋转变换，即可获得切片i+1的边界值(即此处的转换边界值)，因此，通过切片i的贡献范围边界即可迭代出切片i+1的贡献范围边界，从而使得边界贡献范围的计算简单，减少了系统的运算量。

S440、根据重建切片的临界数组和当前投影数据，计算反投影系数，将当前投影数据反投影累加至上一次的重建数据中，生成更新的重建数据；

S450、根据更新的重建数据，获得重建切片的重建图像。本实施例中，FPGA加速电路将生成的更新的重建数据通过PCI-E总线反馈至上位机，上位机根据更新的重建数据进行图像重建并显示重建图像。

作为进一步的改进，本发明的反投影过程在获取当前投影数据对当前重建切片的贡献范围边界的步骤之前还包括如下步骤：

S501、获取当前投影数据及其帧序号，即识别当前投影数据是第几帧投影数据，如，当前投影数据的帧序号可以是1、2、3或n等等。

S502、根据当前投影数据的帧序号判断当前投影数据对重建切片的贡献度，生成用以表示当前投影数据类别的贡献系数，其中，当前投影数据类别包括无贡献数据、部分贡献数据以及完全贡献数据；例如，贡献系数可以分别设定为0、1、#，若贡献系数为0，则表示当前投影数据为无贡献数据；若贡献系数为1，则表示当前投影数据为完全贡献数据；若贡献系数为#，则表示当前投影数据为部分贡献数据。

若贡献系数表示当前投影数据为无贡献数据，则直接略过当前投影数据，进行下一帧投影数据的反投影运算；这样，不仅减小了系统的运算量，同时降低了系统的时间复杂程度和空间复杂程度，使得系统结构简单，易于FPGA实现。

若贡献系数表示当前投影数据为完全贡献数据，则计算反投影系数，根据当前投影数据以及反投影系数逐行对当前重建切片的所有重建点进行反投影累加，生成更新的重建数据。

具体地，若贡献系数表示当前投影数据为完全贡献数据，则逐行对重建点进行反投影累加的步骤包括：

逐行获取当前重建切片的所有重建点的上一次的重建数据；即FPGA加速电路的反投影重建模块首先从存储模块或外部内存模块中逐行读取当前重建切片的相应行的上一次的重建数据。

针对每一个重建点，计算反投影系数；

根据反投影系数将当前投影数据累加至上一次的重建数据中，生成更新的重建数据。

若贡献系数表示当前投影数据为部分贡献数据，则执行步骤S410，即首先确定当前投影数据对当前重建切片的贡献范围边界，根据当前投影数据对当前重建切片的贡献范围边界，逐行对当前重建切片的重建点进行反投影累加，其具体实现方法参见下文中的描述。

作为一种可实施方式，如图2所示，若贡献系数表示为部分贡献数据，则根据贡献范围边界逐行对重建切片的所有重建点进行反投影累加的步骤还包括：

S503、获取当前重建切片的序号；当重建切片的数量为i＝1时，则当前重建切片的序号可以始终为1，即当前重建切片的序号值保持不变。当重建切片的数量为两个以上时，则当前重建切片的序号可以依次表示为1、2、3等等，以根据当前重建切片的序号辨识当前是对第几个重建切片进行重建。

S504、获取当前重建切片的当前重建行的行号，以便于根据当前重建行的行号逐行对当前重建切片的重建点进行反投影运算。

S505、根据当前投影数据的帧序号、当前重建切片的序号以及当前重建行的行号生成所述临界数组的查找地址；

S506、根据查找地址从临界数组中获得与当前重建切片的当前重建行的行号对应的边界值；

S507、判断重建切片的个数，若重建切片为单切片，即当重建切片的数量为1个时，则执行步骤S508，根据当前重建行对应的边界值逐行判断该重建切片的当前重建行是否需要当前投影数据，若是，则执行步骤S511，获取当前重建行的上一次的重建数据，将当前投影数据累加至当前重建行的上一次的重建数据中，生成当前重建行的更新的重建数据；若否，则直接跳过当前重建行，进行下一行的反投影过程，直至遍历当前重建切片的所有行之后继续下一帧投影数据的反投影过程。

若所述重建切片的数量为两个以上时，则执行如下步骤：

S509、根据贡献范围边界的坐标旋转矩阵，针对当前重建切片的当前重建行对应的边界值进行坐标旋转变换，获得当前重建切片的当前重建行对应的转换边界值；

S510、根据当前重建切片的当前重建行对应的转换边界值逐行判断当前重建切片的当前重建行是否需要当前投影数据，若是，则执行步骤S511，获取当前重建切片的当前重建行的上一次的重建数据，将当前投影数据累加至当前重建切片的当前重建行的上一次的重建数据中，生成更新的重建数据。若否，则直接跳过所述当前重建行，进行下一行的反投影过程，直至遍历所述当前重建切片的所有行之后继续下一重建切片的重建，当完成当前投影数据对所有重建切片的反投影运算后，继续下一帧投影数据的反投影过程。

具体地，当重建切片为单个的重建切片时，即重建切片的数量i＝1时，若贡献系数表示当前投影数据为部分贡献数据，则逐行对当前重建切片的重建点进行反投影累加的步骤包括：

获取当前重建切片的当前重建行的行号，在本实施例中，逐行对当前重建切片进行反投影重建，因此，首先应当确认当前是对当前重建切片的第几行重建点进行反投影运算，即获取当前重建行的行号。

根据当前投影数据的帧序号、当前重建行的行号生成临界数组的查找地址；

根据查找地址从临界数组中获得与当前重建行的行号对应的边界值；

根据边界值判逐行断当前重建行的重建点是否需要当前投影数据，若是，则首先从存储模块中获取当前重建行的上一次的重建数据，计算反投影系数，根据反投影系数将当前投影数据累加至上一次的重建数据中，生成更新的重建数据，完成对当前重建行的反投影运算；若否，则直接跳过当前重建行，进行该重建切片的下一行的反投影过程，直至遍历当前重建切片的所有行之后继续下一帧投影数据的反投影过程，以减少系统的运算量，提高运算效率。

作为另一种可实施方式，当重建切片的数量为两个以上时，即重建切片为多切片时，在单切片重建的基础上实现多切片的重建。具体地，若贡献系数表示当前投影数据为部分贡献数据，则逐行对重建切片的重建点进行反投影运算的步骤包括：

设定系统参数，即在全局坐标系的Z轴方向上，设定相邻两个重建切片之间的间隔为等于射线源的相邻两个位置点之间间隔的整数倍，相邻两个重建切片之间的间隔可以为射线源的相邻两个位置点之间的间隔的1倍、2倍或3倍等等。

获取当前重建切片的序号以及当前重建切片的当前重建行的行号；

根据当前投影数据的帧序号、当前重建切片的序号以及当前重建行的行号生成临界数组的查找地址；

根据查找地址从临界数组中获得与当前重建行的行号对应的边界值；

根据贡献范围边界的坐标旋转矩阵，针对边界值进行坐标旋转变换，获得当前重建切片的当前重建行对应的转换边界值；由于s_n对重建切片i的贡献度和s_n+1对重建切片i+1的贡献度一样，此处的贡献度包括投影范围的大小与形状，因此，切片i的边界值与重建切片i+1的边界值之间存在一个旋转关系，通过对重建切片i的边界值进行坐标旋转变换，即可获得重建切片i+1的边界值(即此处的转换边界值)，因此，通过重建切片i的贡献范围边界即可迭代出重建切片i+1的贡献范围边界，从而使得边界贡献范围的计算简单，减少了系统的运算量。

根据当前重建切片的当前重建行对应的转换边界值逐行判断当前重建切片的当前重建行是否需要当前投影数据，若是，则从外部内存模块中获取当前重建行的上一次的重建数据，计算反投影系数，并根据反投影系数将当前投影数据累加至上一次的重建数据中，生成更新的重建数据；若否，则直接跳过当前重建切片的当前重建行，进行当前重建切片的下一行的反投影过程，直至遍历当前重建切片的所有行之后继续下一重建切片的重建，当完成当前投影数据对所有重建切片的反投影运算后，继续下一帧投影数据的反投影过程。

进一步的，基于所述射线源的相邻两个位置点之间的间隔建立相邻两个重建切片之间的贡献范围边界的坐标旋转矩阵的步骤包括：

获取K帧投影数据对应的射线源的K个位置点；

根据所述射线源的K个位置点获得相邻两个位置点之间的角度差Δs，其中，Δs＝2π/(K-1)；

根据转换关系式A＝YB获得相邻的下一个重建切片的转换边界值，其中，A表示所述射线源在第n(n＝1,2……K)个位置点处产生的锥束覆盖对第i个重建切片的贡献范围边界上的一个点的xy平面坐标，B表示所述射线源在第n+1个位置点处产生的锥束覆盖对第i+1个重建切片的贡献范围边界上的一个点的xy平面坐标，Y表示坐标转换矩阵，其中，

$Y=\left\{\begin{array}{cc}cos\frac{2\mathrm{\π}}{(K-1)}& -sin\frac{2\mathrm{\π}}{(K-1)}\\ sin\frac{2\mathrm{\π}}{(K-1)}& \cos \frac{2\mathrm{\π}}{(K-1)}\end{array}\right\}.$

在本发明的另一个实施例中，所述反投影过程还包括如下步骤：

在重建切片确定的范围内遍历所有重建点；

根据参考重建点的坐标计算重建所述重建切片所需射线源的最大位置和最小位置；

根据所需射线源的最大位置和最小位置判断是否需要略过当前投影数据。

例如，重建切片的切片数为P，分辨率为M*N，则重建切片确定的范围为P*M*N，设定参考重建点x的坐标为x＝(x₁,x₂,x₃)，则x₃存在最小值x_min和最大值x_max，根据锥束覆盖原理，当x₃处于最小值x_min时，存在以下关系式：

$\frac{R-x_{1}cos\left(s_0\right)-x_{2}sin\left(s_0\right)}{2R}\stackrel{^bot}{v}+\frac{{\mathrm{hs}}_0}{2\mathrm{\π}}\≤x_{min}---\left(5\right)$

当x₃处于最大值x_max时，存在以下关系式：

$x_{max}\≤\frac{R-x_{1}cos\left(s_0\right)-x_{2}sin\left(s_0\right)}{2R}\stackrel{^top}{v}+\frac{{\mathrm{hs}}_0}{2\mathrm{\π}}---\left(6\right)$

其中，是对应的Tam窗的下边界的值，是对应的Tam窗的上边界的值，R表示射线源的螺旋线轨迹半径；s0表示当前投影数据所对应的射线源的位置点；h表示螺旋轨迹的螺距。

根据上述公式(5)遍历x₁和x₂可以得到所需射线源的最小位置，根据上述公式(6)遍历x₁和x₂可以得到所需射线源的最大位置，根据所需射线源的最大位置和最小位置判断是否需要略过当前投影数据。

以上各个实施例在具体说明中仅只针对相应步骤的实现方式进行了阐述，然后在逻辑不相矛盾的情况下，上述各个实施例是可以相互组合的而形成新的技术方案的，而该新的技术方案依然在本具体实施方式的公开范围内。

如图3所示，基于上述CBCT重建算法，本发明还提供了一种CBCT重建系统，包括预设模块001、数据获取模块002、反投影重建模块003，其中，预设模块001用于在全局坐标系的Z轴方向上，设定相邻两个重建切片之间的间隔等于射线源的相邻两个位置点之间的间隔的非负整数倍；数据获取模块002用于利用所述射线源获取K帧投影数据，其中，K为正整数；反投影重建模块003用于当投影数据对当前重建切片为部分贡献时，针对每帧所述投影数据，进行反投影过程。

具体地，反投影重建模块003包括边界生成单元031、转换矩阵单元032、旋转单元033、反投影运算单元034以及重建单元035。其中，边界生成单元031用于获取当前投影数据对所述重建切片的贡献范围边界，并根据所述贡献范围边界建立临界数组；转换矩阵单元032用于基于所述射线源的相邻两个位置点之间的间隔建立相邻两个重建切片之间的贡献范围边界的坐标旋转矩阵；旋转单元033用于利用贡献范围边界的坐标旋转矩阵，由已重建切片的临界数组迭代获得与该已重建切片相邻的重建切片的临界数组；反投影运算单元034用于根据重建切片的临界数组和当前投影数据，计算反投影系数，将当前投影数据反投影累加至上一次的重建数据中，生成更新的重建数据；重建单元035用于根据更新的重建数据，获得重建切片的重建图像。

本实施例中，预设模块001、数据获取模块002、边界生成单元031、转换矩阵单元032、旋转单元033、反投影运算单元034以及重建单元035，与上述方法中的步骤S100、步骤S200、以及步骤S410、步骤S420、步骤S430、步骤S440以及步骤S450一一对应，用于实现上述步骤。

在本发明的一个实施例中，反投影重建模块003还包括：

用于获取当前投影数据的帧序号的单元；

用于根据所述当前投影数据的帧序号，判断所述当前投影数据对重建切片的贡献度，生成用以表示所述当前投影数据类别的贡献系数的单元，其中，所述投影数据类别包括无贡献数据、部分贡献数据以及完全贡献数据；

若所述贡献系数表示所述投影数据为无贡献数据，则略过当前投影数据；

若所述贡献系数表示所述投影数据为完全贡献数据，则逐行对当前重建切片的所有重建点进行反投影累加；

若所述贡献系数表示所述投影数据为部分贡献数据，则获取每帧所述投影数据对当前重建切片的贡献范围边界，根据所述贡献范围边界逐行对当前重建切片的重建点进行反投影累加。

在本发明的一个实施例中，若所述贡献系数表示为部分贡献数据，则反投影运算单元还包括：

用于获取所述当前重建切片的序号以及所述当前重建切片的当前重建行的行号的单元；

用于根据所述当前投影数据的帧序号、所述当前重建切片的序号以及所述当前重建行的行号生成所述临界数组的查找地址的单元；

用于根据所述查找地址从所述临界数组中获得与所述当前重建行的行号对应的边界值的单元；

用于判断重建切片的个数的单元，若所述重建切片为单切片，则调用第一贡献判断单元，用于根据所述边界值逐行判断该重建切片的当前重建行是否需要所述当前投影数据，若是，则获取所述当前重建行的上一次的重建数据，将所述当前投影数据累加至所述上一次的重建数据中，生成所述更新的重建数据；若否，则直接跳过所述当前重建行，进行下一行的反投影过程，直至遍历所述当前重建切片的所有行之后继续下一帧投影数据的反投影过程；

若所述重建切片的数量为两个以上时，则调用边界值转换单元，用于根据所述贡献范围边界的坐标旋转矩阵，针对所述边界值进行坐标旋转变换，获得转换边界值；

之后调用第二贡献判断单元，用于根据所述转换边界值逐行判断当前重建切片的当前重建行是否需要所述当前投影数据，若是，则获取所述当前重建行的上一次的重建数据，将所述当前投影数据累加至所述上一次的重建数据中，生成所述更新的重建数据，；若否，则直接跳过所述当前重建行，进行下一行的反投影过程，直至遍历所述当前重建切片的所有行之后继续下一重建切片的重建，当完成当前投影数据对所有重建切片的反投影运算后，继续下一帧投影数据的反投影过程。

在本发明的一个实施例中，所述反投影重建模块还包括对每帧所述投影数据依次进行求导、长度校正、前向插值、滤波、后向插值的单元。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品承载在一个非易失性计算机可读存储载体(如ROM、磁碟、光盘，服务器存储空间)中，包括若干指令用以使得一台终端设备(可以是手机，计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法流程和系统架构。

此外，如图4所示，本发明还提供了一种CBCT重建系统，包括：射线源(未示出)、与射线源相对设置的探测器(未示出)、FPGA加速电路200以及上位机100。其中，射线源用于发射锥束射线，射线源做螺旋轨迹运动，每旋转一定角度后，停下来照射一次，获得一帧投影数据。本实施例中，射线源在2π旋转区间内取K个离散的位置点，相应的产生K帧投影数据。探测器用于接收锥束射线穿过被测物体后的射线，并将光信号转化为数字信号，获取K帧投影数据，其中，K为正整数。FPGA加速电路200用于对K帧投影数据进行反投影处理，生成更新的重建数据；上位机100用于根据更新的重建数据，获得重建切片的重建图像并显示，本实施例中，上位机是PC或者计算机集群。

具体地，上位机100对探测器获得的K帧投影数据依次进行求导、长度校正以及前向插值预处理，并将预处理后的所述投影数据传送至FPGA加速电路200。FPGA加速电路200对上位机预处理后的投影数据进行滤波、后向插值以及反投影，生成更新的重建数据，并将更新的重建数据传送至上位机100，用以在上位机100上显示重建图像。本实施例中，FPGA加速电路200中设置有PCI-E接口模块210，上位机通过PCI-E总线连接至PCI-E接口模块210，从而实现上位机与FPGA加速电路的电连接。

其中，FPGA加速电路200包括PCI-E接口模块210、DHT滤波模块220、后向插值模块230、反投影模块240以及存储模块250。PCI-E接口模块210、DHT滤波模块220以及后向插值模块230依次电连接，DHT滤波模块220用于接收上位机100预处理后的投影数据，并对上述投影数据进行DHT滤波处理，之后，后向插值模块230对滤波后的投影数据进行后向插值处理。后向插值模块230连接反投影模块240以及存储模块250，同时反投影模块240以及存储模块250均连接至PCI-E接口模块210。

存储模块250用于保存反投影过程中所述重建切片的重建数据，即存储模块250用于存储上一次的重建数据以及反投影运算后获得的更新的重建数据。当然，存储模块250还用于存储设定的相邻两个重建切片之间的间隔与所述射线源的相邻两个位置点之间的间隔的倍数关系。本实施例中，存储模块250可以是DDR3SDRAM或DDR2SDRAM。

反投影模块240用于基于贡献范围边界的坐标旋转矩阵，进行反投影叠加生成更新的重建数据，反投影模块240包括边界查询单元以及反投影运算单元245。其中，边界查询单元用于获取当前投影数据对所述重建切片的贡献范围边界，并根据所述贡献范围边界建立临界数组，利用贡献范围边界的坐标旋转矩阵，由已重建切片的临界数组迭代获得与已重建切片相邻的重建切片的临界数组。反投影运算单元245用于根据重建切片的临界数组和当前投影数据，将当前投影数据反投影累加至上一次的重建数据中，生成更新的重建数据，并将更新的重建数据传送至存储模块250。

在本发明的一个实施例中，反投影模块240还包括缓存单元242、控制单元241、第一判断单元243，缓存单元242连接PCI-E接口模块210，缓存单元242用于接收上位机100处理后的投影数据，缓存生成一帧投影数据。控制单元241和反投影运算单元245连接缓存单元241，控制单元241用于当缓存满一帧数据时，控制缓存单元241向反投影运算单元245传送当前投影数据，并获取当前投影数据的帧序号。第一判断单元243串联在控制单元241与反投影运算单元245之间，第一判断单元241根据控制单元传送的当前投影数据的帧序号，判断当前投影数据对重建切片的贡献度，生成用以表示所述当前投影数据类别的贡献系数。当贡献系数表示所述投影数据为无贡献数据，则反投影运算单元245略过所述当前投影数据；当贡献系数表示投影数据为完全贡献数据，则控制单元241控制反投影运算单元245逐行对重建切片的所有重建点进行反投影累加；当贡献系数表示所述投影数据为部分贡献数据，则控制单元241控制反投影运算单元245获取当前投影数据对重建切片的贡献范围边界，根据贡献范围边界逐行对重建切片的重建点进行反投影累加。

本实施例中，边界查询单元包括第一地址生成单元246、第一边界查询表单元247以及转换单元248，反投影运算单元245、第一地址生成单元246、第一边界查询表单元247以及转换单元248依次循环连接。其中，第一地址生成单元246用于从反投影运算单元245获取当前重建切片的序号、当前重建切片的当前重建行的行号以及当前投影数据的帧序号，根据当前投影数据的帧序号、当前重建切片的序号以及当前重建行的行号生成临界数组的查找地址。第一边界查询表单元247用于根据查找地址从临界数组中获得与所述当前重建切片的当前重建行的行号对应的边界值。转换单元248用于根据坐标旋转矩阵，针对边界值进行坐标旋转变换，获得当前重建切片的当前重建行对应的转换边界值，并将转换边界值反馈至反投影运算单元245。

当重建切片的数量为一个，可以设定转换单元248中转换矩阵的各个元素均为1。当然，也可以不设置转换单元248，此时，反投影运算单元245、第一地址生成单元246、第一边界查询表单元247依次循环电连接，如图5所示。第一边界查询表单元247将该重建切片的当前重建行对应的边界值反馈至反投影运算单元245。反投影运算单元245用于根据边界值判断该重建切片的当前重建行是否需要当前投影数据，若是，则获取当前重建行的上一次的重建数据，将当前投影数据累加至所述上一次的重建数据中，生成所述更新的重建数据；若否，则直接跳过所述当前重建行，进行下一行的反投影过程，直至遍历所述当前重建切片的所有行之后继续下一帧投影数据的反投影过程。

下面结合图4和图5说明针对单一重建切片进行反投影运算的工作过程：

首选，经过上位机100预处理后的投影数据经PCI-E总线以及PCI-E接口模块210进入DHT滤波模块220，经过DHT滤波模块220的滤波以及后向插值模块230的后向插值运算后进入缓存单元242中进行缓存。当满一帧投影数据时，缓存单元242输出full信号给控制单元241，控制单元241记录当前投影数据的帧序号，并将当前投影数据的帧序号传送至第一判断单元243，同时控制单元241向缓存单元242反馈ctrl信号，使得缓存单元242将当前投影数据传送至反投影运算单元245。第一判断单元243根据当前投影数据的帧序号判断当前投影数据是无贡献数据、完全贡献数据还是部分贡献数据，生成用以表示当前投影数据类别的贡献系数，并将当前投影数据的帧序号以及贡献系数传送至反投影运算单元245。

反投影运算单元245根据贡献系数执行不同的反投影运算过程，具体过程如下：

若贡献系数表示当前投影数据为无贡献数据，则反投影运算单元245直接略过当前投影数据，不进行反投影操作，继续下一帧投影数据的反投影运算。

若贡献系数表示当前投影数据为完全贡献数据，则反投影运算单元245逐行对该重建切片上的所有重建点进行反投影运算。具体地，反投影运算单元245首先逐行从存储模块250中获取上一次的重建数据，之后，针对每一个重建点，将当前投影数据累加值上一次的重建数据中，生成更新的重建数据，之后，反投影运算单元245将更新的重建数据重新写回存储模块250中，当遍历该重建切片的所有行后，反投影运算单元向控制单元发送finish信号，控制单元根据接收的finish信号控制缓存信号输入下一帧投影数据。

若贡献系数表示当前投影数据为部分贡献数据，则反投影运算单元245将当前投影数据的帧序号、当前重建切片的当前重建行的行号传送至边界查询单元的第一地址生成单元246。第一地址生成单元246根据当前投影数据的帧序号以及当前重建的行号生成临界数组的查找地址，之后将查找地址传送至第一边界查询表单元247，第一边界查询表单元247根据上述查找地址获得当前重建行对应边界值，并将边界值反馈给反投影运算单元245。反投影运算单元245根据边界值该重建切片的当前重建行判断是否需要当前投影数据，若是，则反投影运算单元从存储模块250中获取当前重建行的上一次的重建数据，根据上一次的重建数据和当前投影数据计算反投影系数，根据反投影系数将当前投影数据累加至上一次的重建数据中，生成更新的重建数据，完成对当前重建行的反投影运算。若否，则表示当前重建切片不需要该帧投影数据，直接跳过所述当前重建行，进行下一行的反投影过程，以减少系统的运算量，提高运算效率。直至遍历所述当前重建切片的所有行，反投影运算单元输出finish信号给控制单元，之后继续下一帧投影数据的反投影运算。

作为另一种可实施方式，如图6所示，当重建切片的数量为两个以上时，第一判断单元243可以集成于反投影运算单元245内，反投影运算单元245、第一地址生成单元246、第一边界查询表单元247、转换单元248依次循环电连接。转换单元248将当前重建切片的当前重建行对应的转换边界值反馈至反投影运算单元245。此时，反投影运算单元245用于根据转换边界值判断当前重建切片的当前重建行是否需要当前投影数据，若是，则获取所述当前重建行的上一次的重建数据，将所述当前投影数据累加至所述上一次的重建数据中，生成所述更新的重建数据。若否，则直接跳过当前重建行，进行下一行的反投影过程，直至遍历当前重建切片的所有行之后继续下一重建切片的重建，当完成当前投影数据对所有重建切片的反投影运算后，继续下一帧投影数据的反投影过程。

作为进一步的扩展，该CBCT重建系统还包括外部内存模块260，外部内存模块260连接存储模块250，用于存储多个重建切片的上一次的重建数据以及更新的重建数据，以进一步扩展系统的内存，易于FPGA的实现。此时，存储模块250用于存储当前重建切片的上一次的重建数据和更新的重建数据。

下面结合图4和图6说明针对多个重建切片进行反投影运算的工作过程：

首选，经过上位机100预处理后的投影数据经PCI-E总线以及PCI-E接口模块210进入DHT滤波模块220，经过DHT滤波模块220的滤波以及后向插值模块230的后向插值运算后进入缓存单元242中进行缓存。当满一帧投影数据时，缓存单元242输出full信号给控制单元241，控制单元241记录当前投影数据的帧序号，并将当前投影数据的帧序号传送至反投影运算单元245，同时控制单元向缓存单元242反馈ctrl信号，使得缓存单元242将当前投影数据传送至反投影运算单元245。反投影运算单元245内置的第二判断单元根据当前投影数据的帧序号与当前重建切片的序号判断当前投影数据对当前重建切片是无贡献数据、完全贡献数据还是部分贡献数据，生成用以表示当前投影数据类别的贡献系数。

反投影运算单元245根据其生成的贡献系数执行不同的反投影运算过程，具体过程如下：

若贡献系数表示当前投影数据对当前重建切片为无贡献数据，则反投影运算单元245直接略过当前投影数据，不对当前重建切片进行反投影操作，继续判断当前投影数据对下一重建切片的贡献度，直至遍历所有的重建切片。

若贡献系数表示当前投影数据为完全贡献数据，则反投影运算单元245逐行对该重建切片上的所有重建点进行反投影运算。具体地，反投影运算单元245首先逐行从外部内存模块260中读取当前重建切片的上一次的重建数据，并将当前重建切片的上一次的重建数据存储于存储模块250中，之后，针对当前重建切片的每一个重建点，计算反投影系数，并根据反投影系数将当前投影数据累加至当前重建切片的上一次的重建数据中，生成当前重建切片的更新的重建数据。之后，反投影运算单元245将当前重建切片的更新的重建数据重新写入外部内存模块260中，当遍历当前重建切片的所有行后，反投影运算单元245进行下一个重建切片的反投影重建过程，直至遍历所有的重建切片后，反投影运算单元245向控制单元241发送finish信号，控制单元241根据接收的finish信号控制缓存信号输入下一帧投影数据。

若贡献系数表示当前投影数据为部分贡献数据，则反投影运算单元245将当前投影数据的帧序号、当前重建切片的序号、当前重建切片的当前重建行的行号传送至边界查询单元的第一地址生成单元246。第一地址生成单元246根据当前投影数据的帧序号以及当前重建的行号生成临界数组的查找地址，之后将查找地址传送至第一边界查询表单元247，第一边界查询表单元247根据上述查找地址获得当前重建行对应边界值，并将边界值传送至转换单元248，转换单元248利用坐标转换矩阵生成转换边界值，并将该转换边界值反馈给反投影运算单元245。

反投影运算单元245根据该转换边界值判断当前重建切片的当前重建行是否需要当前投影数据，若是，则反投影运算单元245从外部内存模块260中读取当前重建切片的上一次的重建数据，并将当前重建切片的上一次的重建数据存储于存储模块250中。之后，逐行计算反投影系数，根据反投影系数将当前投影数据累加至上一次的重建数据中，生成更新的重建数据，完成对当前重建行的反投影运算，并将当前重建行的更新的重建数据写入存储模块250中。若否，则直接跳过当前重建行，进行下一行的反投影过程，以减少系统的运算量，提高运算效率。直至遍历所述当前重建切片的所有行之后继续下一重建切片的重建，当完成当前投影数据对所有重建切片的反投影运算后，投影运算单元245输出finish信号给控制单元241，继续下一帧投影数据的反投影过程。

本发明的CBCT重建方法及系统，先判断当前投影数据对当前重建切片的贡献度，继而通过贡献范围边界来筛选重建像素点进行反投影，从而简化了系统的计算量，并通过将系统参数设定为相邻两个重建切片之间的间隔等于相邻射线源之间的间隔的非负整数倍，使得在进行多切片重建时，能够利用坐标旋转矩阵，由已重建切片的临界数组迭代获得与已重建切片相邻的重建切片的临界数组，从而简化了多切片重建时的计算过程，提高了反投影运算的效率，降低了系统的复杂程度，易于FPGA实现。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种CBCT重建方法，其特征在于，包括如下步骤：

在全局坐标系的Z轴方向上，设定相邻两个重建切片之间的间隔等于射线源的相邻两个位置点之间的间隔的非负整数倍；

利用所述射线源获取K帧投影数据，其中，K为正整数；

若所述投影数据对当前重建切片为部分贡献时，针对每帧所述投影数据，按照以下步骤进行反投影过程：

获取当前投影数据对所述重建切片的贡献范围边界，并根据所述贡献范围边界建立临界数组；

基于所述射线源的相邻两个位置点之间的间隔建立相邻两个重建切片之间的贡献范围边界的坐标旋转矩阵；

利用所述贡献范围边界的坐标旋转矩阵，由已重建切片的临界数组迭代获得与所述已重建切片相邻的所述重建切片的临界数组；

根据所述重建切片的临界数组和所述当前投影数据，将当前投影数据反投影累加至上一次的重建数据中，生成更新的重建数据；

根据所述更新的重建数据，获得所述重建切片的重建图像。

2.根据权利要求1所述的CBCT重建算法，其特征在于，每帧所述投影数据对应获得一个长度为M的临界数组，其中，M为重建切片的行分辨率；所述临界数组的索引表示所述当前重建切片的行，所述临界数组内的元素值表示相应行的边界点列坐标。

3.根据权利要求1所述的CBCT重建算法，其特征在于，所述反投影过程还包括：

获取当前投影数据的帧序号；

根据所述当前投影数据的帧序号，判断所述当前投影数据对所述重建切片的贡献度，生成用以表示所述当前投影数据类别的贡献系数，其中，所述投影数据类别包括无贡献数据、部分贡献数据以及完全贡献数据；

若所述贡献系数表示所述投影数据为无贡献数据，则略过所述当前投影数据；

若所述贡献系数表示所述投影数据为完全贡献数据，则逐行对所述重建切片的所有重建点进行反投影累加；

若所述贡献系数表示所述投影数据为部分贡献数据，则获取所述当前投影数据对所述重建切片的贡献范围边界，根据所述贡献范围边界逐行对所述重建切片的重建点进行反投影累加。

4.根据权利要求3所述的CBCT重建算法，其特征在于，若所述贡献系数表示当前投影数据为完全贡献数据，则逐行对所述重建点进行反投影累加的步骤包括：

针对每一个重建切片，逐行获取所述重建切片的所有重建点的上一次的重建数据；

针对每一个所述重建点，计算反投影系数；

根据所述反投影系数，将所述当前投影数据累加至所述上一次的重建数据中，生成更新的重建数据。

5.根据权利要求3所述的CBCT重建算法，其特征在于，若所述贡献系数表示为部分贡献数据，则根据所述贡献范围边界逐行对所述重建切片的重建点进行反投影累加的步骤还包括：

获取当前重建切片的序号以及所述当前重建切片的当前重建行的行号；

根据所述当前投影数据的帧序号、所述当前重建切片的序号以及所述当前重建行的行号生成所述临界数组的查找地址；

根据所述查找地址从所述临界数组中获得与所述当前重建切片的当前重建行的行号对应的边界值；

判断所述重建切片的个数，若所述重建切片的数量为一个，则根据所述边界值逐行判断所述重建切片的当前重建行是否需要所述当前投影数据，若是，则获取所述当前重建行的上一次的重建数据，将所述当前投影数据累加至所述上一次的重建数据中，生成所述更新的重建数据；若否，则直接跳过所述当前重建行，进行下一行的反投影过程，直至遍历所述当前重建切片的所有行之后继续下一帧投影数据的反投影过程；

若所述重建切片的数量为两个以上时，则执行如下步骤：

根据所述贡献范围边界的坐标旋转矩阵，针对所述边界值进行坐标旋转变换，获得所述当前重建切片的当前重建行对应的转换边界值；

根据所述转换边界值逐行判断所述当前重建切片的当前重建行是否需要所述当前投影数据，若是，则获取所述当前重建行的上一次的重建数据，将所述当前投影数据累加至所述上一次的重建数据中，生成所述更新的重建数据；若否，则直接跳过所述当前重建行，进行下一行的反投影过程，直至遍历所述当前重建切片的所有行之后继续下一重建切片的重建，当完成当前投影数据对所有重建切片的反投影运算后，继续下一帧投影数据的反投影过程。

6.根据权利要求1所述的CBCT重建算法，其特征在于，所述反投影过程还包括如下步骤：

在所述重建切片确定的范围内遍历所有重建点；

根据所述参考重建点的坐标计算重建所述重建切片所需射线源的最大位置和最小位置；

根据所述所需射线源的最大位置和最小位置判断是否需要略过所述当前投影数据。

7.根据权利要求1所述的CBCT重建算法，其特征在于，基于所述射线源的相邻两个位置点之间的间隔建立相邻两个重建切片之间的贡献范围边界的坐标旋转矩阵的步骤包括：

获取所述K帧投影数据对应的所述射线源的K个位置点；

根据所述射线源的K个位置点获得相邻两个位置点之间的角度差Δs，其中，Δs＝2π/(K-1)；

根据转换关系式A＝YB获得相邻的下一个重建切片的转换边界值，其中，A表示所述射线源在第n(n＝1,2……K)个位置点处产生的锥束覆盖对第i个重建切片的贡献范围边界上的一个点的xy平面坐标，B表示所述射线源在第n+1个位置点处产生的锥束覆盖对第i+1个重建切片的贡献范围边界上的一个点的xy平面坐标，Y表示坐标转换矩阵，其中，

$Y=\left\{\begin{array}{cc}cos\frac{2\mathrm{\π}}{(K-1)}& -sin\frac{2\mathrm{\π}}{(K-1)}\\ sin\frac{2\mathrm{\π}}{(K-1)}& \cos \frac{2\mathrm{\π}}{(K-1)}\end{array}\right\}.$

8.根据权利要求1所述的CBCT重建算法，其特征在于，所述反投影过程之前还包括对每帧所述投影数据依次进行求导、长度校正、前向插值、滤波、后向插值的过程。

9.一种CBCT重建系统，其特征在于，包括：

预设模块，在全局坐标系的Z轴方向上，设定相邻两个重建切片之间的间隔等于射线源的相邻两个位置点之间的间隔的非负整数倍；

数据获取模块，用于利用所述射线源获取K帧投影数据，其中，K为正整数；

反投影重建模块，用于当投影数据对当前重建切片为部分贡献时，针对每帧所述投影数据，进行反投影过程，包括：

边界生成单元，用于获取当前投影数据对所述重建切片的贡献范围边界，并根据所述贡献范围边界建立临界数组；

转换矩阵单元，用于基于所述射线源的相邻两个位置点之间的间隔建立相邻两个重建切片之间的贡献范围边界的坐标旋转矩阵；

旋转单元，用于利用所述贡献范围边界的坐标旋转矩阵，由已重建切片的临界数组迭代获得与所述已重建切片相邻的所述重建切片的临界数组；

反投影单元，用于根据所述重建切片的临界数组和当前投影数据，将当前投影数据反投影累加至上一次的重建数据中，生成更新的重建数据；以及

重建单元，用于根据所述更新的重建数据，获得所述重建切片的重建图像。

10.一种CBCT重建系统，其特征在于，包括：

射线源，用于发射锥束射线；

与所述射线源相对设置的探测器，用于接收锥束射线穿过被测物体后的射线，并将光信号转化为数字信号，获取K帧投影数据，其中，K为正整数；

FPGA加速电路，用于对K帧投影数据进行反投影处理，生成更新的重建数据；以及

上位机，用于根据所述更新的重建数据，获得所述重建切片的重建图像并显示，其中，所述FPGA加速电路包括：

存储模块，用于保存反投影过程中所述重建切片的重建数据，以及设定的相邻两个重建切片之间的间隔与所述射线源的相邻两个位置点之间的间隔的倍数关系；

反投影模块，用于基于贡献范围边界的坐标旋转矩阵，进行反投影叠加生成更新的重建数据，所述反投影模块包括：

边界查询单元，用于获取当前投影数据对所述重建切片的贡献范围边界，并根据所述贡献范围边界建立临界数组，利用所述贡献范围边界的坐标旋转矩阵，由已重建切片的临界数组迭代获得与所述已重建切片相邻的所述重建切片的临界数组；以及

反投影运算单元，用于根据所述重建切片的临界数组和所述当前投影数据，将当前投影数据反投影累加至上一次的重建数据中，生成所述更新的重建数据。