多源Tomosynthesis迭代算法在动态模型中的应用

【作 者】彭家驹，赵俊

上海交通大学生物医学工程系,医学图像处理实验室，上海市，200240

Tomosynthesis又叫数字融合X射线成像技术，1972年由David G.Grant提出[1]。此后，各种重建算法相继提出，如G Lauritsch等人提出并证明了tomosynthesis的滤波反投影算法[2]，Beilei Wang等人提出了迭代算法[3]，James T.Dobbins等人提出了矩阵翻转法[4]，还有日本的平野圭藏等人提出了频域反卷积成像法（FDDI，Frequency Domain Deconvolution Imaging）[5]。这些算法在胸腺，乳腺及关节成像中都取得了很好的效果[1,6-8]。

Tomosynthesis的出现，很好地补充了传统CT和X光成像的不足。Tomosynthesis与传统的CT成像相比，不需要全方位的投影，具有辐射剂量低的优势，通过几次投影的数据就能得到任意断层的图像。

在提高采集速度方面，为了抑制运动所引起的伪迹，Jun Zhao等人曾提出了CT下的三源扫描算法[9,10]，参考这一算法，本文提出了一种多源的Tomosynthesis扫描方法，并在此基础上将多源重建算法应用到动态模型重建中。通过比较单源、双源与四源的重建结果，发现多源重建算法在动态重建方面的性能优于单源重建算法，而四源重建又比双源重建效果更好。

1 方法与材料

1.1 迭代tomosynthesis重建算法

在描述多源Tomosynthesis重建之前，先简单介绍单源扫描下的重建算法。一般情况下，由于Tomosynthesis投影较少(一周10-20个投影)，往往无法达到解析重建的要求。因此，为得到较清晰的图像，本文将采用迭代重建法(Algebraic Reconstruction Technique, ART)。迭代法的缺点是计算量大，因而重建耗时较长，但随着计算机运算速度的提高，这个缺点将逐渐被克服[11]。下面简要介绍经典的单源迭代重建算法[3]。

我们用如下的方程定义一次投影得到的数据：

其中，P是测量的投影集，F是待测物体体素点的密度矩阵，W是体素点的加权系数矩阵，n是测量过程中的随机噪声。这个方程可以扩展为一个线性方程组形式：

其中，M表示有M条X射线的线积分，即探测器探测单元的个数；N表示待重建区域被离散化以后体素点的个数。

迭代过程是一个不断反馈矫正的过程，对于第i条X射线上的第j个体素，我们用如下的方程式计算第k+1次迭代结果：

其中，λ是松弛系数，一般而言，λ是一个远小于1，大于0的数，调整λ的大小可以控制迭代过程的收敛速度。我们的目标是使代价函数J(f)=│P-W.F│2最小。

迭代过程可以简单的表述为以下4个步骤：

(1) 猜测函数的初始值，通常将其初始化为全0；

(2) 用给的投影计算公式(2)计算一条X射线的投影估测值；

(3) 计算实际测量的投影值和计算出的投影估测值之差并运用公式(3)更新结果，用松弛因子控制收敛速度；

(4) 重复2、3步，直到代价函数J(f)达到最小，或者两次结果间的差值小于给定的阈值，此时迭代结束。

1.2 多源Tomosynthesis扫描结构

图1 两种扫描结构Fig.1 Two scanning con fi guration

双源和四源tomosynthesis扫描方式如图1所示，平行光射线源在圆周上对称分布，每个X射线源阵列发出的X射线穿过被测对象后，分别被下方对应放置的接收器接收。当X射线源转动时，下方的检测器作同步反向转动。由于使用了多源-多探测器的采集方式，能成倍地提高采集速度。对于动态模型，因为采集速度提高，在一周的采集时间内物体运动的时间减少，运动幅度减小，因而由运动带来的伪迹被抑制。

1.3 动态模型定义

本文定义如下模型来评价重建算法在运动中的成像效果。模型由1段圆管和4个圆球组成（图2），其参数由表1给出。投影过程中，外层圆管静止不动，位于中央的4个小球以两种方式变化：(1)小球圆心绕着圆柱的中轴每秒转动3.6°；(2) 小球半径每秒膨胀1.05倍。在采集过程中，由于中央的小球连续运动，采集投影时与静止不动时相比会产生偏差，而基于多源的采集则能有效的抑制这种因运动带来的偏差。

图2 动态数字模型Fig.2 Dynamic digital phantom

表1 模型参数Tab.1 Model parameters

表1中x, y, z 分别表示圆管内外壁和圆球中心空间位置的归一化坐标，r表示归一化半径，h是圆管的归一化高度，ρ是圆管/圆球的吸收系数。

1.4 仿真实验参数设置

仿真实验中投影数据采用解析法求得，为了保证重建结果的可比性，单源，双源与四源都采用如表2所示的参数设置，值得注意的是，在检测器尺寸相同的情况下，双源、四源Tomosynthesis的采集速度分别是单源Tomosynthesis的2倍和4倍，对投影的采集时间缩短到单源情形下的1/2和1/4，虽然一次照射的辐射剂量变大了，但是病人接受辐射的总剂量并没有发生变化。

2 模型仿真结果及定量分析

表2 仿真参数Tab.2 Simulation parameters

2.1 不同扫描方式下的仿真重建结果

根据上述图像重建原理，采用解析法计算动态模型的投影值，再分别用单源，双源和四源迭代算法进行重建，松弛系数λ=0.025，迭代30次，对模型正中间z=0的断面显示，结果如图3和图4所示。图3中(a)、(b)、(c)分别对应了单源，双源和四源扫描下旋转模型的重建结果；图4中(a)、(b)、(c)分别对应了单源，双源和四源扫描下膨胀模型的重建结果。从重建结果可以看出，采用多源tomosynthesis能够有效地抑制运动伪迹。

图3 单源、双源、四源重建算法重建结果比较（转动模型）Fig.3 Reconstruction results of single-source, dual-source and four-source(Rotation phantom)

图4 单源、双源、四源重建算法重建结果比较（膨胀模型）Fig.4 Reconstruction results of single-source, dual-source and four-source(Expansion phantom)

2.2 对重建结果进行定量分析

为了进一步比较重建图像与原始图像的差异，我们选用下面的三个评价参数[11]进行定量比较：

(1) 归一化均方距离判据D

(2) 归一化平均绝对距离判据R

(3) 最坏情况距离判据E

其中，tij是原始值，rij是计算值，t是原始值的平均值，N是图像的边长。

对两种扫描方式下的重建结果，分别计算三种定量指标，不同迭代次数下的结果如图5和图6所示，图中纵坐标表示误差的大小，横坐标表示迭代的次数。从图中可以看到，随着迭代次数的增加，D，R，E的变化总体呈现下降趋势， 并且误差减小的速度和最终的收敛结果，是四源重建优于双源重建优于单源重建。

图5 旋转模型下不同扫描方式的定量对比Fig.5 Quantitative analysis of different scanning con fi guration in rotation phantom

图6 膨胀模型下不同扫描方式的定量对比Fig.6 Quantitative analysis of different scanning con fi guration in expansion phantom

3 结论

本文提出了一种新的多源Tomosynthesis扫描方法，在测量参数及仿真模型一样的情况下，比较单源，双源与四源的采集速度，迭代重建次数及定量误差。结果表明，多源扫描方式（尤其是四源扫描）在动态模型仿真中能加快收敛速度，减小误差，并有效地抑制动态伪迹，该方法对进一步改善心脏及呼吸运动成像具有指导意义。

[1] David G.Grant. tomosynthesis: A three-dimensional radiographic imaging Technique[J]. IEEE Transactions on Biomedical Engineering, 1972, 19(1): 20-28.

[2] Gunter Lauritsch, Wolfgang H. Harer. A theoretical framework for filtered backprojection in tomosynthesis[C]. Proc. SPIE, 1998,3338: 1127-1137.

[3] Beilei Wang, Kenneth Barner, Denny Lee. Algebraic tomosynthesis reconstruction[C]. Proc. of SPIE, 2004, 5370: 711-718.

[4] Richard J. Warp, Devon J. Godfrey, James T. Dobbins III.Applications of matrix inverse tomosynthesis[C]. Proc. SPIE, 2000,3977: 376-383.

[5] 平野圭藏. 频域反卷积成像(FDDI): 一种新的X线断层摄影法[D].日本:信州大学, 1993.

[6] Devon J. Godfrey, Amber Rader, James T. Dobbins III. Practical strategies for the clinical implementation of matrix inversion tomosynthesis[C]. Proc. SPIE, 2003, 5030: 379-390.

[7] Maravilla K R, Murry R C Jr, Dech M, et al. Clinical application of digital tomosynthesis: a preliminary report[J]. Am. J.Neuroradiology, 1983, 4(3): 277-280.

[8] Richard L. Webber, Hunter R. Underhill, Paul F. Hemler, et al.A Nonlinear algorithm for task-specific tomosynthetic image reconstruction[C]. Proc. SPIE, 1999, 3659: 258-265.

[9] Jun Zhao, Ming Jiang, Tiange Zhuang, et al. Minimum detection window and inter-helix PI-line with triple-source helical cone-beam scanning[J]. Journal of X-Ray Science and Technology, 2006, 14:95-107.

[10] Jun Zhao, Ming Jiang, Tiange Zhuang, et al. An exact reconstruction algorithm for triple-source helical cone-beam CT[J].Journal of X-Ray Science and Technology, 2006, 14: 191-206.

[11] 庄天戈. CT原理与算法[M]. 上海: 上海交通大学出版社, 1992