Suppression Method for Cone-Beam CT Artifact Based on FDK Compensation and Dual-Source Weighting
-
摘要:
FDK算法具有结构简单、计算速度快、重建质量高等优点，至今仍是锥束CT解析重建的主流算法。然而，若重建点锥角过大，会出现密度下降的锥束伪影。为此学界提出了许多抑制伪影的方法，包括投影加权、束线重排与引入补偿项等，其中引入补偿项的方法计算效率高，并在锥角为十余度的情形取得了良好的伪影抑制效果，在诸多算法中具有一定的优势。然而，对于更大的锥角情形，由于数据不完备的固有缺陷，该算法的重建的效果并不十分理想，因此需要对重建算法与扫描几何作进一步优化，以适应实际的工业或临床需求。本文基于FDK重建的两个补偿项，将锥束CT几何拓展至z向分布的双光源情形，通过减小重建点平均锥角的方法进一步抑制锥束伪影，并提出一种双光源重建的锥角加权方法，将双光源的重建结果进行融合，得到最终重建图像。仿真实验结果表明，与单光源情形相比，本文的算法可以更加有效地提高重建图像的质量，从而有望在更大锥角的CBCT重建系统中得到应用。
Abstract:The Feldkamp-Davis-Kress (FDK) algorithm has the advantages of a simple structure, fast calculation speed, and high reconstruction quality, and it remains the mainstream algorithm for cone-beam computed tomography (CBCT) analytical reconstruction. However, if the cone angle for the reconstructed point is too large, the reconstruction exhibits artifacts characterized by an axial density drop. Numerous methods for suppressing this artifact have been proposed, including projection weighting, cone-beam rebinning, and introducing compensation terms. The compensation method has high computational efficiency and achieves good artifact suppression when the cone angle is not more than 16°. However, for cases with larger cone angles, the effectiveness of this algorithm is diminished. Thus, further optimization of the reconstruction algorithms and scanning geometry is required to accommodate industrial and clinical demands. Based on the two compensation terms for FDK reconstruction, this study extends cone-beam CT geometrically to the case of dual source in the z-direction distribution, further suppressing cone-beam artifacts by reducing the average cone angle of the reconstructed points. A cone-angle weighting method is proposed for dual-source reconstruction, to fuse the reconstructed results of the two sources and obtain the final reconstructed image. The simulation results show that the proposed algorithm can improve the quality of the reconstructed image more effectively than single-source case.
-
Keywords:
- FDK /
- cone-beam artifact /
- dual-source weighting
-
X射线计算机断层成像（computed tomography，CT）作为一种物质内部结构的无损探测手段，经历了从二维成像到三维成像的发展过程。其中，三维锥束CT具有射线利用率高、空间分辨率高等优点，在口腔CT等领域有着广阔的应用前景。
早在20世纪80年代，Feldkamp等[1]就针对圆轨道锥束CT提出了一类近似重建算法，被后人称之为FDK算法。FDK算法具有经典的滤波反投影结构，计算速度快、重建质量高，至今仍是大部分圆轨道锥束CT的主流重建算法。然而，由于圆轨道锥束扫描不满足Tuy[2]和Smith[3]数据充分条件，FDK算法对大锥角处的重建并不准确，会出现密度下降的伪影，这类伪影通常被称为锥束伪影（cone beam artifact）。学者们提出了相应的许多抑制此伪影的方法，主要包括投影加权、束线重排与引入补偿项等。其中，投影加权的方法主要包括HT-FDK[4]、CW-FDK[5]、Parker-FDK[6]、BPW-FDK[7]等，主要思想是在滤波前或滤波后对来自不同方向的投影加权，以降低大锥角射线重建带来的误差；束线重排的方法包括P-FDK[8]、T-FDK[9]、v-FDK[10]、C-FDK[11]等，这些算法的思想是通过将锥束数据重排为一系列平行的扇束，并改变
$ z $ 向坐标映射关系，以提高重建效率与精度，束线重排方法本质上是改变滤波路径，并且可以进一步设计滤波路径参数可调的重建算法[12]；引入补偿项的方法则是通过对Radon空间数据完备性进行分析，在FDK重建的基础上加上额外项，以修正FDK无法利用的数据产生的伪影，其中最具代表性的为Hu[13]和Zhu[14]提出的补偿项$ {f}_{H} $ 和$ {f}_{Z} $ ，他们的算法对锥角不超过16度的锥束伪影抑制效果较好，并且计算复杂度不高，在诸多算法中具有一定的优势，然而，对于更大的锥角情形，由于数据不完备的固有缺陷，该算法重建的效果并不十分理想，因此需要对重建算法与扫描几何作进一步优化，以适应实际的工业或临床需求。在扫描几何方面，学界曾提出过诸多新型扫描系统焦点轨道，如“圆弧+圆弧[15]”、“圆弧+直线[16]”、“圆弧+螺旋[17]”等，这些轨道可以使系统采集到更完备的投影数据，然而对应的重建算法也随之变得更加复杂，使得重建效率不高。近年来，各类静态CT的蓬勃发展推动着分布式X光源技术的逐渐成熟[18]。所谓分布式X光源，其靶点通常沿直线排列，通过控制系统使各靶点交替出束，采集来自不同靶点的投影信号。这种出束几何为进一步抑制锥束伪影、增大成像视野提供了新的思路：使分布式光源靶点沿z向分布，则各重建点都能找到与z向距离最近的两个靶点，对于这两个靶点来说，该重建点的平均锥角不会太大，Radon空间数据缺失问题不严重，因此加权重建后的锥束伪影也应更弱。
本文的研究建立在Hu等[13]和Zhu等[14]的工作基础之上，注意到分布式光源的在锥束伪影抑制方面的优势，以两个理想X光源代表分布式光源，提出一种基于双光源-面阵探测器几何的简洁高效的成像算法，并通过仿真实验验证算法对此类伪影的抑制效果。
1. FDK算法与补偿项
FDK算法专用于圆轨道锥束CT的解析重建，其扫描几何关系如图1所示：源探距离为
$ D $ 的光源$ S $ 与面阵探测器相对，沿半径为$ R $ 的圆形轨道旋转，$ \beta $ 为旋转角参数，$ x,\;y,\;z $ 为物空间坐标，$ {x}_{\beta },\;{y}_{\beta },\;{z}_{\beta } $ 为旋转角$ \beta $ 下的坐标，$ u $ 和$ v $ 为探测器坐标，$ {H}_{D} $ 为探测器高度，光源$ S $ 发出的锥束光线由平板探测器接收，投影数据为$ q(u,\;v,\;\beta ) $ 。为表示方便，这里设$ u $ 和$ v $ 为中心虚拟探测器坐标，重建的目标是由投影数据$ q(u,\;v,\;\beta ) $ 计算物空间衰减系数$ f\left(x,\;y,\;z\right) $ 。FDK算法给出的重建算法为
$$ \begin{aligned} {f}_{FDK}\left(x,\;y,\;z\right)=& \dfrac{1}{2}{\int }_{0}^{2\pi }\dfrac{{R}^{2}}{{\left(R-{y}_{\beta }\right)}^{2}}{\int }_{-\infty }^{\infty }\dfrac{R}{\sqrt{{R}^{2}+{u}^{2}+{v}^{2}}}\\& q\left(u,\;v,\;\beta \right)h\left({u'}-u\right)\mathrm{d}u\mathrm{d}\beta 。\\[-5pt]\end{aligned} $$ (1) 其中
$ h\left(u\right) $ 是斜坡滤波器，在本文的仿真模拟中为S-L滤波器，各符号的计算公式如下$$\begin{split} &{u'}=\frac{R}{R-{y}_{\beta }}{x}_{\beta },\;v=\frac{R}{R-{y}_{\beta }}z,\\& {x}_{\beta }=x\cos\beta +y\sin\beta ,\;{y}_{\beta }=-x\sin\beta +y\cos\beta。 \end{split} $$ (2) 由于FDK是由平面等距扇束的滤波反投影公式近似而来，在
$ z $ 较大的点重建误差大，Hu[13]证明这些误差源自于Radon空间的数据缺失，且FDK的重建结果大体上可以认为是把缺失部分看作零时的Radon逆变换结果。因此，衰减系数的重建结果$ f $ 可以分为3部分，分别为：①FDK重建项$ {f}_{FDK} $ ；②圆轨道锥束扫描中可以得到但没有被FDK公式利用到的Radon数据的逆变换$ {f}_{H} $ ；③圆轨道锥束扫描中缺失的Radon数据的逆变换$ {f}_{N} $ ，即$$ \begin{array}{c}f={f}_{FDK}+{f}_{H}+{f}_{N}。\end{array} $$ (3) Hu利用Radon逆变换与Grangeat公式[19]对
$ {f}_{H} $ 进行了数学推导，得出$$ \begin{split} {f}_{H}\left(x,\;y,\;z\right)=\;& -\dfrac{1}{2\pi }{\int }_{0}^{2\pi }\dfrac{z}{{\left(R-{y}_{\beta }\right)}^{2}}\Bigg[\dfrac{1}{2\pi }\dfrac{\partial }{\partial v}{\int }_{-\infty }^{\infty }\\& \dfrac{R}{\sqrt{{R}^{2}+{u}^{2}+{v}^{2}}}q\left(u,\;v,\;\beta \right)\mathrm{d}u\Bigg]\mathrm{d}\beta 。\end{split} $$ (4) 公式（4）中
$ v $ 和$ {y}_{\beta } $ 的定义与公式（2）中相同，在此基础之上，Yang[20]等利用Radon空间插值的方法得到了$ {f}_{N} $ 的一个近似项，然而，由于这个近似项滤波是移变的，庞大的计算量限制了其实际应用。Zhu等[14]在此基础之上，进一步利用平行束近似对缺失的Radon空间数据进行了估计，推导出$ {f}_{N} $ 的一个形式简单的近似项$ {f}_{Z} $ ，计算公式为$$ \begin{split} {f}_{Z}\left(x,\;y,\;z\right)=\;& -\dfrac{1}{4{\pi }^{2}}\dfrac{{z}^{2}+{R}^{2}}{{R}^{2}}\left(1-\dfrac{\sqrt{{R}^{2}-{z}^{2}}}{R}\right){\int }_{0}^{2\pi }\dfrac{{\partial }^{2}}{\partial {z}^{2}}\\& \left[{\int }_{-\infty }^{\infty }\dfrac{R}{\sqrt{{R}^{2}+{u}^{2}+{z}^{2}}}q\left(u,\;z,\;\beta \right)\mathrm{d}u\right]d\beta 。\end{split} $$ (5) 在
$ {f}_{Z} $ 表达式中需要对数据沿$ z $ 向存在求二阶导的操作，对应地，$ {f}_{H} $ 中需要求一阶导，这些求导会在补偿项的$ z $ 向引入高频噪声。由于锥角伪影主要是低频成分，因此在实际操作时，可以适当对补偿项进行中值滤波，以减少补偿重建结果的$ z $ 向噪声。可以看出，由于投影数据
$ q\left(u,\;v,\;\beta \right) $ 都经历了预加权、滤波、后加权与反投影等过程，三项的计算结构具有许多相似之处，也存在一些差异：①预加权因子相同，都是一个角度的余弦，这个角度是每个探测器单元对应的射线路径与中心射线路径的夹角；②滤波过程有区别，$ {f}_{FDK} $ 项使用斜坡滤波器对$ u $ 高通滤波，$ {f}_{H} $ 项先对变量$ u $ 积分、再对$ v $ 求一阶导，$ {f}_{Z} $ 项先做变量替换$ z=v $ 、接着对变量$ u $ 积分、再对$ z $ 求二阶导；③后加权因子不同，由于$ {f}_{FDK} $ 项和$ {f}_{H} $ 项反投影前要用到$ {y}_{\beta } $ ，其后加权因子与重建点和角度都有关，而$ {f}_{Z} $ 项的后加权因子则只与重建点的$ z $ 坐标有关，这使得这一加权的过程可以从反投影积分中提出，因此$ {f}_{Z} $ 项的计算复杂度相对前两者更低。通过引入$ {f}_{H} $ 与$ {f}_{Z} $ 两个补偿项，FDK算法的锥束伪影可以得到较好的抑制。2. 双光源加权方法
对于更大锥角的重建问题，双光源锥束CT有更大的优势，双光源的锥束CT扫描几何与单光源情形差别不大，区别在于将
$ S $ 换成了沿$ z $ 轴对称分布的$ {S}_{1} $ 和$ {S}_{2} $ ，其间距为$ H $ （图2）。从探测器投影数据重建三维图像的过程，可以看作是向Radon空间填充数据进行Radon逆变换的过程，而填充数据的过程可以借助Grangeat公式实现。图3展示了
$ R=10\;\mathrm{c}\mathrm{m}，H=2\;\mathrm{c}\mathrm{m} $ 时单光源与双光源CBCT的Radon空间数据分布情况，在不考虑投影截断的情况下，单光源圆轨道CBCT的Radon空间数据分布形状为图3（a）所示的圆环体，这一圆环体由图中单视角投影的球壳旋转而来，若要重建中心一球形区域，数据完备性要求Radon数据也应该填满此区域，单光源情形的Radon数据无法满足此区域，会有中心部分较大的数据缺失区（shadow zone），正是这部分数据缺失导致重建结果沿z向密度下降，产生了锥束伪影。对于本文提出的双光源锥束CT系统，Radon空间数据分布形状发生了改变，变为了图3（c）所示的旋转体。单光源CT存在一个中心平面，双光源CT系统存在两个中心平面，位于中心平面上的点$ P $ 重建所需的Radon数据是完备的，而偏离中心平面上的点$ Q $ 的Radon数据是缺失的，可以看出，光源$ {S}_{1} $ 投影数据未能填充的部分一定程度上通过$ {S}_{2} $ 的投影数据得到了弥补，从而改善了圆轨道锥束CT数据不完备的问题，抑制了锥束伪影。![]()
图 3 圆轨道CBCT的Radon空间数据分布Figure 3. Data distribution in Radon space of circular cone-beam computed tomography (CBCT)值得注意的是，根据上面的分析，即使双光源有望降低锥束伪影，但其测量数据仍然是不完备的，只是不完备部分区域有所减小，而数据缺失的区域和边界区域仍然需要使用两个补偿项进行估计。双光源的Radon空间数据存在重合部分，这部分自然需要加权叠加，而直接在Radon空间加权，并通过Radon逆变换进行双光源系统的三维重建过程复杂、重建时间长，且Radon空间直接插值也会引起额外的误差和伪影。因此本文不在Radon空间而在物重建空间进行加权处理，并在加权前分别对两个光源下单独的重建结果进行补偿，此方法在减少计算量的同时，也可以使重建结果更加稳定。
本文提出的双光源加权思路如下：首先定义重建点相对于两个光源的锥角
$ {\alpha }_{1} $ 和$ {\alpha }_{2} $ ，再设计以$ {\alpha }_{1} $ 和$ {\alpha }_{2} $ 为变量的权重函数$ {w}_{1} $ 与$ {w}_{2} $ ，最后将带补偿项$ {f}_{H} $ 和$ {f}_{Z} $ 的双光源各自解析重建结果以权重$ {w}_{1} $ 和$ {w}_{2} $ 叠加，得到最终的重建图像。对于重建点
$ P(x,\;y,\;z) $ ，如图4所示，连接$ {S}_{1} $ 与$ P $ 、$ {S}_{2} $ 与$ P $ 分别与探测器平面相交于$ {P}_{1 D} $ 和$ {P}_{2 D} $ ，在探测器直线$ u=0 $ 上$ {P}_{1 D} $ 和$ {P}_{2 D} $ 的正投影点分别为$ M $ 和$ N $ ，设$ {S}_{1} $ 与$ {S}_{2} $ 在探测器平面的正投影点分别为$ {S}_{1 D} $ 和$ {S}_{2 D} $ ，本文中定义$ P $ 相对于两个光源$ {S}_{1} $ 和$ {S}_{2} $ 的锥角分别为有向角$ {\alpha }_{1}=\angle {S}_{1 D}{S}_{1}M, {\alpha }_{2}=\angle {S}_{2 D}{S}_{2}N $ ，因此图中的$ {\alpha }_{1} < 0,\;{\alpha }_{2} > 0 $ 。两个锥角只与$ \beta $ 角度下的侧视图参数$ z $ 与旋转坐标$ {y}_{\beta } $ 有关，容易推导其计算式分别如下$$ \left\{\begin{array}{l} {\alpha }_{1}\left(x,\;y,\;z,\;\beta \right)={\mathrm{tan}}^{-1}\displaystyle\frac{z-H/2}{R-{y}_{\beta }}\\ \qquad\qquad\quad\;\;\; ={\mathrm{tan}}^{-1}\displaystyle\frac{z-H/2}{R+x\mathrm{s}\mathrm{i}\mathrm{n}\beta -y\mathrm{c}\mathrm{o}\mathrm{s}\beta }\\ {\alpha }_{2}\left(x,\;y,\;z,\;\beta \right)={\mathrm{tan}}^{-1}\displaystyle\frac{z+H/2}{R-{y}_{\beta }}\\ \qquad\qquad\quad\;\;\; ={\mathrm{tan}}^{-1}\displaystyle\frac{z+H/2}{R+x\mathrm{s}\mathrm{i}\mathrm{n}\beta -y\mathrm{c}\mathrm{o}\mathrm{s}\beta }\end{array}\right.。 $$ (6) 对于能够被射线照射到的重建点，
$ {\alpha }_{1} $ 和$ {\alpha }_{2} $ 的取值范围为$$ \begin{array}{c}\left\{\begin{array}{c}{\alpha }_{1}\in \left[-{\mathrm{tan}}^{-1}\displaystyle\frac{{H}_{D}+H}{2D},\;{\mathrm{tan}}^{-1}\displaystyle\frac{{H}_{D}-H}{2D}\right]\\ {\alpha }_{2}\in \left[-{\mathrm{tan}}^{-1}\displaystyle\frac{{H}_{D}-H}{2D},\;{\mathrm{tan}}^{-1}\displaystyle\frac{{H}_{D}+H}{2D}\right]\end{array}\right.。\end{array} $$ (7) 由于
$ {\alpha }_{1} $ 和$ {\alpha }_{2} $ 都是$ (x,\;y,\;z,\;\beta ) $ 的多元函数，二者并不独立。其中$ {\alpha }_{1} < 0 $ 且$ {\alpha }_{2} > 0 $ 的重建点高度位于两个光源之间，是需要将$ {S}_{1} $ 与$ {S}_{2} $ 重建结果加权叠加的部分；对于超出某一光源照射区域的重建点，只利用来自另一光源的投影数据，也即对应的权重因子为0。设$ {S}_{1} $ 与$ {S}_{2} $ 加权重建的权重因子分别为$ {w}_{1}({\alpha }_{1},\;{\alpha }_{2}) $ 和$ {w}_{2}({\alpha }_{1},\;{\alpha }_{2}) $ ，则在各区域的重建点$ {w}_{1} $ 与$ {w}_{2} $ 的分布如图5所示。作为权重因子，
$ {w}_{1} $ 和$ {w}_{2} $ 是$ ({\alpha }_{1},\;{\alpha }_{2}) $ 在对应定义域内的二元函数，需满足以下条件：（1）归一化：在公共定义域范围内，
$ {w}_{1}+ {w}_{2}=1 $ ；（2）缺失权重处理：若
$ {\alpha }_{1} $ 超出定义域，则$ {w}_{1}=0 $ ；若$ {\alpha }_{2} $ 超出定义域，则$ {w}_{2}=0 $ ；（3）锥角依赖性：在
$ {S}_{1} $ 与$ {S}_{2} $ 的公共照射区域，若$ {\alpha }_{2} $ 固定，则$ {|\alpha }_{1}| $ 越小，$ {w}_{1} $ 越大；若$ {\alpha }_{1} $ 固定，则$ \left|{\alpha }_{2}\right| $ 越小，$ {w}_{2} $ 越大；此外，由双光源几何关系，$ {\alpha }_{1} > 0 $ 时，$ {w}_{1}=1 $ ；$ {\alpha }_{2} < 0 $ 时，$ {w}_{2}=1 $ ；（4）连续性：除了有限点集以外，
$ {w}_{1} $ 与$ {w}_{2} $ 应是$ {\alpha }_{1} $ 和$ {\alpha }_{2} $ 的连续函数，特别地，在边界处连续且任意方向导数连续。在实际计算时，只需要考虑图5中绿色区域的权重因子即可，一个满足条件的权重计算公式为
$$ \begin{array}{c}\left\{\begin{array}{c}{w}_{1}({\alpha }_{1},\;{\alpha }_{2})=\displaystyle\frac{{\left[{\alpha }_{2}\left({\alpha }_{1}-{\alpha }_{1\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right]}^{2}}{{\left[{\alpha }_{2}\left({\alpha }_{1}-{\alpha }_{1\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right]}^{2}+{\left[{\alpha }_{1}\left({\alpha }_{2}-{\alpha }_{2\mathrm{m}\mathrm{a}\mathrm{x}}\right)\right]}^{2}}\\ {w}_{2}({\alpha }_{1},\;{\alpha }_{2})=\displaystyle\frac{{\left[{\alpha }_{1}\left({\alpha }_{2}-{\alpha }_{2\mathrm{m}\mathrm{a}\mathrm{x}}\right)\right]}^{2}}{{\left[{\alpha }_{2}\left({\alpha }_{1}-{\alpha }_{1\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right]}^{2}+{\left[{\alpha }_{1}\left({\alpha }_{2}-{\alpha }_{2\mathrm{m}\mathrm{a}\mathrm{x}}\right)\right]}^{2}}\end{array}\right.\end{array} $$ (8) 这里的
$ {\alpha }_{1\mathrm{m}\mathrm{i}\mathrm{n}} $ 和$ {\alpha }_{2\mathrm{m}\mathrm{a}\mathrm{x}} $ 分别是绿色区域中$ {\alpha }_{1} $ 的最小值和$ {\alpha }_{2} $ 的最大值，当$ {\alpha }_{1\mathrm{m}\mathrm{i}\mathrm{n}}=-10^\circ，{\alpha }_{2\mathrm{m}\mathrm{a}\mathrm{x}}=10^\circ $ 时，在所考虑区域的权重因子分布如图6所示，$ {w}_{1} $ 和$ {w}_{2} $ 在该区域外是0或1的分片常数，本文所提出的权重函数在边界较为平坦，保证了加权叠加时重建图像的平滑过渡。![]()
图 6 本文提出的权重因子$ {w}_{1} $ 与$ {w}_{2} $ Figure 6. Proposed weight factors$ {w}_{1} $ and$ {w}_{2} $ 具体的重建流程如图7所示，两个光源带补偿项的FDK重建先独立进行，再将二者加权叠加，根据公式（6），锥角
$ \alpha $ 是重建角度$ \beta $ 与重建点$ (x, y,\;z) $ 的函数，进一步，重建权重$ w $ 也是二者的函数，因此加权过程适合在反投影的过程中完成。由于单光源FDK可重建的区域是钻石形旋转体区域[8]，在靠近FOV边界处会出现反投影的截断，通过适当的加权，可以处理这些点处的不连续性，从而把两个光源的单独重建的公共区域平滑融合。加权之前，来自两个光源的重建结果在各自的中心水平面上是准确的，基于前面对权重因子的讨论，加权以后在这两个平面分别对应$ {w}_{1}=1 $ 与$ {w}_{2}=1 $ ，因此也是精确重建。在大于$ S1 $ 高度的区域，有$ {w}_{1}=1 $ ，加权结果与$ S1 $ 带补偿的FDK重建结果相等；在小于$ S2 $ 高度的区域，有$ {w}_{2}=1 $ ，加权结果与$ S2 $ 带补偿的FDK重建结果相等;在两个平面之间的区域，$ 0 < {w}_{1},\;{w}_{2} < 1 $ ，加权结果与介于二者各自重建结果之间，并可由FDK的重建误差进行估计。在本文仿真实验的系统中，我们将$ {S}_{1} $ 与$ {S}_{2} $ 对称放置，以便于计算与分析，由于重建误差仅与重建点与光源的相对位置有关，本文加权算法实际并不严格要求光源对称放置，并且很容易推广到$ z $ 方向多光源的系统之中。对于光源间距$ H $ ，我们将其设置到尽量保证$ {S}_{1} $ 与$ {S}_{2} $ 均匀照射，使平均的FDK误差较小。从满足性质（1）~（4）的约束条件而言，权重因子的选择显然不是唯一的，例如考虑下面更一般的表达形式
$$ \begin{array}{l}\left\{ \begin{array}{l}{w}_{1}\left({\alpha }_{1},\;{\alpha }_{2}\right)\displaystyle\frac{{f}_{1}\left({\alpha }_{1}\right){g}_{2}\left({\alpha }_{2}\right)}{{f}_{1}\left({\alpha }_{1}\right){g}_{2}\left({\alpha }_{2}\right)+{f}_{2}\left({\alpha }_{1}\right){g}_{1}\left({\alpha }_{2}\right)}\\ {w}_{2}\left({\alpha }_{1},\;{\alpha }_{2}\right)=\displaystyle\frac{{f}_{2}\left({\alpha }_{1}\right){g}_{1}\left({\alpha }_{2}\right)}{{f}_{1}\left({\alpha }_{1}\right){g}_{2}\left({\alpha }_{2}\right)+{f}_{2}\left({\alpha }_{1}\right){g}_{1}\left({\alpha }_{2}\right)}\end{array}\right. 。 \end{array} $$ (9) 其中各子函数满足
$ {f}_{1}({\alpha }_{1\mathrm{min}})=0,\;{f}_{2}({\alpha }_{2\mathrm{max}})=0, \; {g}_{1}\left(0\right)= 0,\;{g}_{2}\left(0\right)=0 $ ，当$ {f}_{1},\;{f}_{2},\;{g}_{1},\;{g}_{2} $ 具有较好的单调性质且在零点导数值为0时，所计算的权重因子总是符合要求的。根据前面的讨论，所有满足所述性质权重因子重建误差表现形式与本文的结果应当类似，而如何选择更为最合适的权重因子使得重建的绝对误差最小，可以通过参数最优化等方法完成，而对于通常的工程应用问题而言，往往需要在精度与计算复杂度之间实现平衡。本文所提出的权重因子表达式（8），正是当$ {f}_{1},\;{f}_{2},\;{g}_{1},\;{g}_{2} $ 均为二次函数时的特殊情形，具有形式简单、计算高效、结果平滑等优点，从而有一定的实用价值。3. 数值仿真
本文使用计算机仿真实验对提出的算法性能进行验证，所搭建双光源CT系统的参数如表1所示，实验结果还与单光源CT系统进行了对比，除了光源所在位置位于
$ z=0 $ 处以外，单光源CT系统的其他几何参数与双光源系统相同。单光源系统的最大锥角约为17度，更换为双光源后，最大锥角约为11度。实验使用的数字体模包括Shepp-Logan体模与Disfrise-Discs体模，这两个体模由一系列椭球构成，形状参数与Turbel[8]与Zhu等[14]的文献中取相同值，体模的直径均取60 mm。Shepp-Logan体模是人体头部的近似模型，所包含的椭球密度相近，适合验证成像算法对低对比度分辨率的重建能力；Disfrise-Discs体模包含$ z $ 方向上均匀分布的一系列密度相同的圆盘体，适合验证成像算法对$ z $ 向高频信息的恢复能力；两个体模被广泛应用于FDK相关算法的数值仿真实验之中。表 1 本文中的双光源CT系统参数Table 1. Parameters of the dual-source CT system in this study参数 值 旋转半径$ R $ 100 mm 源探距离$ D $ 200 mm 面阵探测器高度$ {H}_{D} $ 130 mm 双光源间距$ H $ 20 mm 探测器单元数$ {N}_{D}^{2} $ 256×256 重建体素数$ {N}_{V}^{3} $ 128×128×128 重建立方体区域边长$ W $ 60 mm 全扫描角度采样数$ {N}_{\theta } $ 400 图8是使用各种重建算法对对Shepp-Logan体模和Defrise-Discs体模的重建结果，对应的重建真值展示在了最左侧一列，并将本文的算法与FDK、T-FDK、Zhu的方法等进行了对比。对于直接FDK重建，可以观察到明显的锥束伪影，并且由于Defrise-Discs体模的Radon空间z向有更多的高频部分，而这部分数据在圆轨道锥束扫描几何中是缺失的，因此其锥束伪影也比Shepp-Logan体模更明显；T-FDK通过将滤波路径更改为平行束中心虚拟探测器平面上的一条水平线，使锥束伪影得到一定程度上的抑制，但从Defrixe-Discs体模的重建结果可以看到，抑制的效果并不十分理想；在Zhu的方法中，由于
$ {f}_{H} $ 项与$ {f}_{Z} $ 项的补偿，两个体模重建结果中的锥束伪影均得到更好的抑制，但在更大锥角处表现仍然欠佳；最右侧一列是带补偿的双光源FDK重建加权叠加的结果，$ {S}_{1} $ 和$ {S}_{2} $ 的高度分别为$ {h}_{1}=10\;\mathrm{m}\mathrm{m},\;{h}_{2}=-10\;\mathrm{m}\mathrm{m} $ ，因此其各自的FDK重建中心平面位于$ z={h}_{1} $ 和$ z={h}_{2} $ 处，在这两个中心平面上的重建是精确的，偏离中心平面越远，密度下降越显著。可以看到，本文的方法基于双光源$ z $ 向分布的优点，改变了锥束伪影的分布特性，相对于单光源各类成像算法的图像质量得到了进一步提高。![]()
图 8 Shepp-Logan体模和Defrise-Discs体模的重建结果，x-z方向，y=−7.7 mm，Shepp-Logan视窗：[0.98，1.05]，Defrise-Discs视窗：[0.7，1.1]Figure 8. Reconstruction result of Shepp-Logan phantom and Defrise-Discs phantom, x-z views. Shepp-Logan Display window: [0.98, 1.05], Defrise-Discs display window: [0.7, 1.1]为了更准确地量化重建质量，图9展示两个体模重建结果
$ z $ 方向上的profile曲线图，黑色虚线为真值，其余实线是各类上述算法的结果。对于单光源FDK重建，$ z=0 $ 处的重建结果是准确的，偏离此中心平面越远，密度下降越明显；单光源的T-FDK重建结果密度下降得到了一些缓和，引入$ {f}_{H} $ 与$ {f}_{Z} $ 后，Shepp-Logan体模上的重建已接近真值，但Defrise-Discs体模结果仍不理想。对于本文的双光源加权补偿重建，$ z={h}_{1} $ 和$ z={h}_{2} $ 处的重建结果都是准确的，$ {h}_{2} < z < {h}_{1} $ 区域在加权因子$ {w}_{1} $ 和$ {w}_{2} $ 的作用下使重建误差尽可能小，在此区域以外，$ z $ 向距离邻近光源越远，密度下降越显著。总的来说，双光源加权重建改变了锥束伪影的分布特性，使得重建的平均伪影程度得以降低；从图9中可以看出，双光源加权重建结果相对于单光源的各类算法，与真值的平均距离有了显著减小。![]()
图 9 Shepp-Logan体模与Defrise-Discs重建结果的profile曲线Figure 9. Profile curve of Shepp-Logan and Defrise-Discs phantoms4. 结论
本文基于FDK锥束伪影的两个补偿项
$ {f}_{H} $ 与$ {f}_{Z} $ ，提出了一种双光源锥束CT成像系统，并给出了该系统的图像重建方法：对于每一个重建点，首先分别利用每一个光源的投影数据计算出带补偿项的FDK重建结果，接着计算重建点相对于两个光源的锥角$ {\alpha }_{1} $ 、$ {\alpha }_{2} $ 和对应的权重$ {w}_{1} $ 、$ {w}_{2} $ ，根据此权重对双光源各自的重建结果加权叠加，得到最终的重建图像。仿真结果表明，相对于单光源的重建结果，本文的双光源锥束CT可以在两个补偿项的基础之上进一步减小锥束伪影，从而能对更大锥角、更深视野的区域进行三维成像。本文的方法还具有以下优点：①计算简单，加权在图像域的反投影过程完成，权重因子形式简单、计算复杂度低；②与经典解析算法适配性高，适用于单光源的具有滤波反投影结构的锥束CT重建算法，均可利用本文的加权策略推广至双光源情形；③可迁移性好，容易进一步拓展到其他领域，如分布式多光源锥束CT、能谱CT等，从而在工业与临床医学中更好地发挥作用。
权重因子
$ {w}_{1} $ 与$ {w}_{2} $ 的设计部分借鉴了扇束短扫描重建中Parker窗[21]的思路，本文给出了一个权重的经验公式，有一定的主观性，旨在指明双光源加权的可行性与在锥束伪影抑制方面的独特优势。如何优化权重因子与设计更合适的成像系统参数等，仍是需要进一步研究的内容。 -
![]()
图 3 圆轨道CBCT的Radon空间数据分布
Figure 3. Data distribution in Radon space of circular cone-beam computed tomography (CBCT)
![]()
图 6 本文提出的权重因子
$ {w}_{1} $ 与$ {w}_{2} $ Figure 6. Proposed weight factors
$ {w}_{1} $ and$ {w}_{2} $ ![]()
图 8 Shepp-Logan体模和Defrise-Discs体模的重建结果，x-z方向，y=−7.7 mm，Shepp-Logan视窗：[0.98，1.05]，Defrise-Discs视窗：[0.7，1.1]
Figure 8. Reconstruction result of Shepp-Logan phantom and Defrise-Discs phantom, x-z views. Shepp-Logan Display window: [0.98, 1.05], Defrise-Discs display window: [0.7, 1.1]
![]()
图 9 Shepp-Logan体模与Defrise-Discs重建结果的profile曲线
Figure 9. Profile curve of Shepp-Logan and Defrise-Discs phantoms
表 1 本文中的双光源CT系统参数
Table 1 Parameters of the dual-source CT system in this study
参数 值 旋转半径$ R $ 100 mm 源探距离$ D $ 200 mm 面阵探测器高度$ {H}_{D} $ 130 mm 双光源间距$ H $ 20 mm 探测器单元数$ {N}_{D}^{2} $ 256×256 重建体素数$ {N}_{V}^{3} $ 128×128×128 重建立方体区域边长$ W $ 60 mm 全扫描角度采样数$ {N}_{\theta } $ 400 
下载: 导出CSV
-
[1] FELDKAMP L A, DAVIS L C, KRESS J W. Practical cone-beam algorithm[J/OL]. Journal of the Optical Society of America. A, Optics, image science, and vision, 1984, 1(6): 612. DOI: 10.1364/JOSAA.1.000612.
[2] TUY H K. An inversion formula for cone-beam reconstruction[J/OL]. SIAM Journal on Applied Mathematics, 1983, 43(3): 546-552. DOI: 10.1137/0143035.
[3] SMITH B D. Image reconstruction from cone-beam projections: Necessary and sufficient conditions and reconstruction methods[J/OL]. IEEE Transactions on Medical Imaging, 1985, 4(1): 14-25. DOI: 10.1109/TMI.1985.4307689.
[4] GRASS M, KöHLER T, PROKSA R. Angular weighted hybrid cone-beam CT reconstruction for circular trajectories[J/OL]. Physics in Medicine & Biology, 2001, 46(6): 1595-1610. DOI: 10.1088/0031-9155/46/6/301.
[5] MORI S, ENDO M, KOMATSU S, et al. A combination-weighted Feldkamp-based reconstruction algorithm for cone-beam CT[J/OL]. Physics in Medicine & Biology, 2006, 51(16): 3953-3965. DOI: 10.1088/0031-9155/51/16/005.
[6] LIU Y, LIU H, WANG Y, et al. Half-scan cone-beam CT fluoroscopy with multiple X-ray sources[J/OL]. Medical physics (Lancaster), 2001, 28(7): 1466-1471. DOI: 10.1118/1.1381549.
[7] 杨宏成, 高欣, 张涛. 基于FDK反投影权重的锥束DSA重建算法[J/OL]. 江苏大学学报(自然科学版), 2013, 34(2): 190-195. DOI: 10.3969/j.issn.1671-7775.2013.02.012. [8] TURBELL H. Cone-Beam Reconstruction Using Filtered Backprojection[D]. ProQuest Dissertations & Theses, 2001.
[9] GRASS M, KöHLER T, PROKSA R. 3D cone-beam CT reconstruction for circular trajectories[J/OL]. Physics in Medicine & Biology, 2000, 45(2): 329-347. DOI: 10.1088/0031-9155/45/2/306.
[10] 王琛. 圆轨道锥束滤波反投影重建算法研究[D]. 中国工程物理研究院, 2020. DOI: 10.27498/d.cnki.gzgwy.2020.000073. [11] LI L, XING Y, CHEN Z, et al. A curve-filtered FDK (C-FDK) reconstruction algorithm for circular cone-beam CT[J]. Journal of X-Ray Science and Technology, 2011, 19(3): 355-371. DOI: 10.3233/XST-2011-0299.
[12] 穆子扬, 卢荣胜, 何攀, 等. 基于滤波路径变换的板状物X射线三维重建算法[J]. 光学学报, 2024, 44(9): 312-325. [13] HU H. An improved cone-beam reconstruction algorithm for the circular orbit[J/OL]. Scanning, 1996, 18(8): 572-581. DOI: 10.1002/sca.4950180807.
[14] ZHU L, STARMAN J, FAHRIG R. An efficient estimation method for reducing the axial intensity drop in circular cone-beam CT[J/OL]. International Journal of Biomedical Imaging, 2008, 2008(1): 242841-242841. DOI: 10.1155/2008/242841.
[15] ZAMBELLI J, NETT B E, LENG S, et al. Novel c-arm based cone-beam CT using a source trajectory of two concentric arcs[C/OL]//Proceedings of SPIE: volume 6510. 2007: 65101Q-65101Q-10. DOI: 10.1117/12.713813.
[16] ZAMYATIN A A, KATSEVICH A, CHIANG B S. Exact image reconstruction for a circle and line trajectory with a gantry tilt[J/OL]. Physics in Medicine & Biology, 2008, 53(23): N423-N435. DOI: 10.1088/0031-9155/53/23/N02.
[17] SCHOMBERG H, VAN DE HAAR P, BAATEN W. Cone-beam CT using a c-arm system as front end and a spherical spiral as source trajectory[C/OL]//Proceedings of SPIE: volume 7258. 2009: 72580F-72580F-12. DOI: 10.1117/12.811545.
[18] 张涛. 新型静态CT成像理论与重建算法研究[M]. 北京: 清华大学出版社, 2023. [19] Grangeat P. Mathematical framework of cone beam 3D reconstruction via the first derivative of the Radon transform[J]. Lecture Notes in Mathematics -Springer-verlag-, 1990, 1497. DOI: 10.1007/bfb0084509.
[20] YANG H Q, LI M H, KOIZUMI K, et al. FBP-type cone-beam reconstruction algorithm with radon space interpolation capabilities for axially truncated data from a circular Orbit[J]. Medical Imaging Technology, 2006, 24(3): 201-208.
[21] PARKER D L. Optimal short scan convolution reconstruction for fanbeam CT[J/OL]. Medical physics (Lancaster), 1982, 9(2): 254-257. DOI: 10.1118/1.595078.
计量
- 文章访问数: 251
- HTML全文浏览量: 22
- PDF下载量: 48
