锥形束CT解析算法进展的研究

    锥形束CT解析算法进展的研究 (pdf)

    ［摘要］  近年来锥形束CT解析法重构有了突破性进展，螺旋CT的非移变滤波反投影（FBP）算法首先由Katsevich提出，并得到不断完善。随后，这一重构系统被推广至普适轨道，同时，在此基础上又衍生出反投影滤波（BPF）的新思路。本文提出了锥形束CT争析算法发展中的关键问题，并作了深入剖析，对比了FBP与BPF算法的优缺点，指出了未来研究发展的要点。

    ［关键词］  锥形束CT；解析算法；滤波反投影算法；反投影滤波算法；Katsevich类算法

      锥形束CT的解析算法一直是三维体积CT领域的重要课题。锥形束重构属于弱病态问题［1］，数值计算方面的困难重重。理论上的公式虽然严格完备，却难以应用于实际设备，所以当前的CT设备仍采用2.5维的Z轴堆叠的空间重构。

    真正意义下的三维体积重构研究在近年有了突破性进展。2002年Katsevich提出了基于螺旋轨道的移不变滤波反投影（FBP）算法［2～4］，锥形束重构研究由此进入新阶段。Katsevich类的重构系统从数值仿真到系统实现的研究工作广泛展开，文献［5］基于实际探测器几何形态详细地讨论了Katsevich法重构系统的实现。随后，为改进重构精度Katsevich提出了3PI算法［6］。同样是源于Katsevich类算法，Pan小组引入Hilbert变换（HT）提出重构的新思路，即反投影滤波（FBP）算法［7，8］。相比FBP算法，BPF在横向截断投影数据情形下仍能获取更好的重构效果，因而在感兴趣区域重构方面有着广阔的应用前景。另一方面，螺旋轨道情形下的重构公式与一些定理也被推广到普适轨道的通用系统［9～11］。新轨道的开拓与基于新轨道重构算法实现也是当下重要的研究内容［12，13］。

    本文提出了锥形束CT解析算法发展中的若干关键问题，地比了FBP与BPF算法的优缺点，指出了未来发展的研究要点。

    1  Katsevich类FBP算法

    1.1  锥形束重构公式的困难  在Katsevich之前，锥形束FBP重构算法的主要困难是：锥形束变换和三维Radon变换不对等［14］，而通过等价关系变换过程［15］中，存在一个非一一映射的变换，从而变换后的重构公式表达为滤波反投影的形式时，滤波的过程是移变的。其中锥形束变换（或称Xray变换）表达为：

    P为锥形束投影算子，Λ为实轴上某区间，S2为三维实空间中单位球。

    2002年，Katsevich首先在螺旋轨道的情形下给出了非移变的FBP重构算法公式。这一公式的建立是基于螺旋轨道特性，并在滤波方向上有了很大的改进。

    1.2  Katsevich螺旋CT非移变FBP算法

    1.2.1  螺旋CT FBP初始公式  设螺旋轨道C：=R3：x=Rcosλ,y=Rsinλ,z=λ（h/2π）,λΛ,Katsevich的公式利用了螺旋轨道一个重要性质：对于轨道内的任一重构点x-，必然存在唯一的连接轨道上两源点的线段。该线段称为PIline记为LPI（x-），其对应的轨道曲线为CPI（x-）。给定重构点（x-），对CPI（x-）上各源点（λ，x-）的投影数据g（λ,θ）做滤波得到gF后，就可以再通过反投影重构出该点的密度f（x-）,即：

    关键问题就是对于每个源点（λ，x-）如何滤波得到相应的gF。Fatsevich最先提出的滤波算法［2］包含了两个方向上的滤波，记（λ，x-）指向x-的单位矢量为x（λ,x-），则：

    其中ek（λ，x-）=1，2表示与滤波有关的两方方向，e1（λ，x-）与（λ，x-）处的切线方向及x（λ,x-）两个矢量方向垂直，e2（λ，x-）与CPI（x-）相切，且与x（λ，x-）共面。其重构的最终结果是这两部分反投影的均值。

    1.2.2  改进的螺旋CT FBP算法  随后Katsevich滤波方向选择上做了改进，仅采用一个方向滤波，从而有效地减少了计算量［3］（图1）。对于给定点x-和某源点（λ0，x-），（λ*，x-）和［（λ0+λ*）/2，x-］确定的平面过点x-，该平面记为Kplane，法线方向记为u（λ0,λ*），这一滤波方向由u（λ0，λ*）矢量决定的。改进的滤波法就是用e（λ,x-）:=χ（λ,x-）×u（λ,x-）替代（3）中的ek（λ,x-），即：

    若采用平板探测器，那么（4）中的表达式cosγχ（λ,x-）+sinγe（λ，x-）指出了滤波就是沿Kplane和投影面的交线。在实际算法实现时，先将投影所得的矩阵上的数据沿滤波方向重排为行向量，再对各个行向量作滤波。特别地，（4）式体现了这一滤波具有Hilbert变换的形态。

    1.2.3  普适轨道的FBP算法  为改进精确的重构效果，Katsevich又提出了3PI算法，拓展了扫描角度，但是这一拓展使得上述情形中的u（λ0，λ*）并不唯一了，必须引入权因子处理多个μi（λ0,λ*）,0＜i＜M的积分。类同于3PI算法中的权因子的思路，Katsevich提出了普适轨道的滤波反投影的重构公式，所谓的PIline演化为Mline，CPI（x-）上的反投影积分对应为CMline（x-）。       其中（λ，x,e）相应的权因子，μm对应于k（λ,x）的不连续点，而cm（λ,x）是不连续点处的阶跃值。然而（5）还不是万能公式，在应用于具体的轨道时，仍需处理诸如“如何确定Kplane”、“如何设置权因子”等问题。

    图1  Katsevich螺旋锥形束CT滤波反投影算法示意图

    图2  螺旋锥形束CT反投影滤波算法示意图

    2  基于HT的BPF算法

    2.1  BPF算法  改进后的Katsevich算法启发了新算法的产生。考虑交换（5）中两个积分的顺序，即先做反投影，再做滤波，Pan小组指出了这样的变换是成立的，并由此提出了反投影滤波（BPF）的新思路（图2）。（x-）由gb（x-′）做Hilbert变换得到，其中K（x-,x-′）为Hilbert变换核，e（x-）为关于x-的PIline的方向；gb（x-′）是投影g（x-,λ）的加权反投影积分，即：

    这说明BPF算法中的Hilbert变换对应的滤波是沿PIline方向的，重构是在PIline上实现的。那么，若仅需重构物体的某局部，可以只计算与该区域相交的各条PIline线段上的各点的密度即可。

    2.2  BPF算法与FBP算法的性能比较  若投影数据集完备，这两种算法是一致的，不过实际系统中的投影数据集在横向与纵向常常是截断的，两种算法在处理纵向截断数据方面都具有良好的鲁棒性，但在横向截断时存在差异。

    FBP算法的滤波沿Kplane与投影面的交线，也就是说某个点x-的重构涉及Kplane上所有的其他点。若探测器在横向较为狭窄，就会导致横向数据截断，在重构中形成伪影。另一方面，BPF法的公式中蕴含了优越的局部特性（图2），gb（x-′）仅与CPI（x-）上各源点的x-点处投影的结果有关，Hilbert变换沿LPI（x-）方向，也是仅依赖于x-。这一局部特性使得BPF在横向及纵向截断投影数据情形下人都能获取更好的重构效果，在感兴趣区域（ROI）重构方面有着广阔的应用前景。

    3  进一步研究要点

    BPF算法的提出与完善给ROI方向的研究开拓了广阔的领域，其发展必定同时关注最小数据集重构与高分辨率的冗余扫描重构。前者有助于最大可能降低放射线源剂量及减小探测器的面积。

    另一方面，探求与ROI相匹配的重构轨道也是今后研究的要点之一。比如在脏器官的CT扫描中，若采用更灵活的源点轨道，使其几何特征与脏器官在体内不同的位置与形态相匹配，有希望提高重构系统的性能。

    4  结论

    近年来锥形束CT解析法重构有了突破性进展，Katsevich提出了螺旋轨道CT的非移变滤波反投影（FBP）公式及其改进形式。随后，这一重构系统推广至普适轨道。在此基础上Pan变换了反投影与滤波的积分顺序，基于Hilbert变换建立了反投影滤波（BPF）算法。两者在滤波方式上存在着较大的不同，FBP算法相对完善，数值计算精度较高，而BPF所具有的局部特性使其在ROI领域中有着旷阔的应用前景。基于锥形束CT的ROI研究将成为今后该领域的研究要点，包括最小数据集重构、冗余数据处理、自适应轨道等方面内容。

    ［参考文献］

    1  Adel Faridani.Introduction to the mathematics of computed tomography.Inverse Problems,2003,47:1-45.

    2  Katsevich A.Theoretically exact filtered backprojectiontype inversion algorithm for spiral CT.SIAM Appl Math,2002,162（6）:2010-2026.

    3  Katsevich A.An improved exact filtered backprojection algorithm for spiral computed tomography.Advances in Applied Mathematics,2004,32:681-697.

    4  Katsevich A.Analysis of an exact inversion algorithm for spiral conebeam CT.Phys Med Biol,2002,47:2583-2598.

    5  Noo F,Pack J,Heuscher D.Exact helical reconstruction using native conebeam geometries.Phys Med Biol,2003,48:3787-3818.

    6  Katsevich A.3PI algorithms for spiral CT.Adv Appl Math,2006,36:213-250.

    7  Yu Zou,Xiaochuan Pan.Exact image reconstruction on PIlines from minimum data in helical conebeam CT.Phys Med Biol,2004,49:941-959.

    8  Pack JD,Noo F,Clackdoyle R.Conebeam reconstruction using the backprojection of locally filtered projections.IEEE Trans Med Imag,2005,24（1）:70-84.

    9  Katsevich A.A general scheme for constructing inversion alforithms for cone beam CT.Int J Math Math Sci,2003,21:1305-1321.

    10  Pack JD,Noo F.Conebeam reconstruction using 1D filtering along the projection of Mlines.Inverse Problem,2005,21:1105-1120.

    11  Yangbo Ye.A general exact reconstruction for conebeam CT via backprojectionfiltration.IEEE Trans Med Imag,2005,24（9）:1190-1198.

    12  Alexander Katsevich.Image reconstruction for the circle and line trajectory.Phys Med Biol,2005,50:2249-2265.

    13  Haiquan Yang,Meihua Li,Kazuhito Koizumi,et al.Exact cone beam reconstruction for a saddle trajectory.Phys Med Biol,2006,51:1157-1172.

    14  Smith BD.Conebeam tomography:recent advances and a tutorial review.Optical Engineering,1990,29:524-534.

    15  Grangeat P.Mathematical framework of conebeam reconstruction via the first derivative of the Radon transform.Herman GT,Louis AK,Natterer F.Lecture Notes in Mathematics,1991,1497:66-97.

    *基金项目：广东自然科学基金（编号：05006593）

    作者单位： 510641 广东广州，华南理工大学电子与信息学院

　 （编辑：丁剑辉）