基于衰减补偿的CBCT短扫描实时重建方法

【作 者】杨坪坪，冯汉升,，许继伟，杨洋，宋云涛,

1 中国科学技术大学，合肥市，230026

2 中科离子医学技术装备有限公司，合肥市，230088

3 中国科学院等离子体物理研究所，合肥市，230031

0 引言

锥形束CT（cone beam computed tomography,CBCT）[1-2]因锥形束射线利用率高，采集效率高，能有效减少辐射剂量，提高成像效率，提高层间分辨率而应用于口腔、头部、胸腔以及盆腔成像中。在临床应用中，每日就诊人数较多，每次患者扫描进行图像重建的时间有限，因此对CBCT系统的重建算法运行速度要求较高。基于短扫描的重建方式所需投影帧数更少，重建时间减短，且扫描剂量较少，因此成为当前研究热点。由于数据冗余，在扫描过程中没有必要让探测器做360o全扫描，扫描范围可以小于360o，这种扫描方式叫做短扫描。基于短扫描加权的FDK（Feldkamp）图像重建算法[3]既能保证重建速度，扫描剂量也更少，因此商业应用潜力更加巨大。1982年Parker首先提出了圆轨迹短扫描扇形束概念[4]，指出适当地对短扫描原始数据进行加权，重建得到的图像质量基本上等价于使用全扫描投影数据进行图像重建得到的图像质量，并提出了一种二维加权方式，获得了较好的短扫描重建图像质量。同时由于锥形束成像采用圆扫描轨迹，投影数据不完全，因此FDK重建图像存在强度衰减现象[5-6]，研究普遍采用反投影加权进行强度衰减补偿，其中倒高斯加权和倒余弦加权公式重建效果最佳，但都各有缺点，倒高斯加权重建水平向衰减补偿以及倒余弦加权重建轴向衰减补偿不佳。根据已有的衰减补偿算法的不足，本研究提出新的衰减补偿算法，在提高倒高斯加权轴向衰减补偿的同时提高水平向衰减补偿，有效地提高了重建图像质量；并针对基于三维重建开源库RTK库（Reconstruction Toolkit，医学图像重建库）实现的短扫描加权FDK算法进行流程改进，克服了基于RTK库一次读入所有投影图像的缺点，实时性地每次读入一帧投影图像进行重建，同步并行执行患者采集和三维重建过程，大大减少从图像开始采集到生成图像重建之间的总耗时，同时应用新的衰减补偿算法改善图像质量。

1 FDK算法原理

三维图像重建的定义是从物体的投影数据得到物体的内部断层图像的过程，FDK算法是其中应用最广泛的算法，需要利用投影图像数据与重建几何位置信息，实施加权、滤波、反投影三个步骤即可完成重建。

1.1 锥形束图像重建FDK算法

FDK算法是由Feldkamps等[7-8]在1984年提出，针对锥形束投影数据进行近似重建。该算法可以由扇形束的FBP（滤波反投影）算法推导得到。

针对圆轨迹的FDK算法公式如下:

式中，βmin为初始视角，βmax为最终视角，R为物体中心到射线源距离，Z为探测器平面上探测器像素对应探测器坐标系的距离。p(t,l,β)⊕h(t)行滤波后，经加权β)，即可反投影得到重建图像。

其中，p(t,l,β)为投影数据，t是探测器上的行坐标，l为探测器上的列坐标，h(t)为一维斜坡滤波核，参数几何位置关系见图1。

图1 锥形束示意图Fig.1 Schematic diagram of cone beam

1.2 短扫描加权FDK算法

根据1982年Parker首先提出的圆轨迹扇形束短扫描的概念以及加权概念，短扫描重建需要在加权基础上乘上额外的Parker加权因子。

针对第i帧的Parker加权因子wi(α,β)公式如下：

其中：α是投影数据位置相对于重建物体中心角度，β是X射线源对应角度。△是扇形束角度，其参数几何位置关系，如图2所示。

图2 Parker加权因子中各参数对应几何关系Fig.2 Geometric relationship in the Parker weighting factor

2 衰减补偿算法改进

2.1 算法介绍及改进

FDK重建域中，中心断层物体能被精确重建，非中心断层物体的重建精度随与中心断层距离迅速下降，其形状貌似山顶函数，称为强度衰减现象，如图3所示。

图3 均匀密度球模型横截面512层与FDK重建断层图像横截面对比Fig.3 Comparison of the cross section of the uniform density sphere model (512 layers) and the cross section of the FDK reconstruction central tomographic image

由此，学者们提出许多函数进行强度衰减补偿，其中最有效的是倒高斯加权函数[9]和倒余弦加权函数[5]。其一般形式如下:

倒高斯加权函数：

其中：z∈[-zmax,zmax]

倒余弦加权函数：

其中，x、y、z分别为重建点坐标，k、c1、c2分别为倒高斯加权函数和倒余弦加权函数的加权因子，D为射线源到重建中心点距离。根据文献[9]表明，倒高斯函数最优加权因子k=20.5；倒余弦加权函数最优加权因子c1=1.4，c2=1。倒高斯函数相比较倒余弦函数有更好的轴向衰减补偿效果；而倒余弦函数因为考虑了重建点到原点距离r因素，因此有更好的水平向衰减补偿效果。

由上，本研究提出新的加权公式，兼顾倒高斯加权函数与倒余弦加权函数的优势，在倒高斯函数的基础上加上距离r因素，称为改进倒高斯加权函数，如下：

2.2 改进算法比较

实验采用均匀球密度模型，利用Matlab进行仿真重建。重建角度从0o到209o，每隔1o取一张投影，锥角取28.72o。重建体素512×512×512个，体素尺寸0.08 mm×0.08 mm×0.08 mm，探测器像素512×512个，探测器尺寸51.2 mm×51.2 mm，源点到重建物体中心距离80 mm，到探测器中心距离100 mm。在医学成像的仿真实验中，常用归一化均方误差评价重建图像质量，定义为:

其中：xi为物体模型真实强度值，为重建图像像素值。

共设置5组数据，即原始球数据、FDK重建数据、倒高斯加权数据、倒余弦加权数据、改进倒高斯加权数据，其中倒高斯加权、倒余弦加权均采用最优加权因子，改进倒高斯加权采用加权因子k1=20.5，k2=1.4。对比如图4～图7所示。

图4 X轴第256层（均匀球模型)、X轴第244层（SheepLogan[10]模型）模型重建横截面图像对比（窗口[0.9 1.25]）Fig.4 X-axis layer 256 (uniform sphere model),X-axis layer 244 (SheepLogan[10] model) model reconstruction cross-sectional image comparison (window [0.9 1.25])

图5 轴向x=256,y=256扫描线强度对比Fig.5 Axial x=256,y=256 scan line intensity comparison

图6 水平向y=256,z=128扫描线强度对比Fig.6 Scanning line intensity comparison of horizontal y=256,z=128

图7 水平向x=256, z=128扫描线强度对比Fig.7 Scanning line intensity comparison of horizontal x=256,z=128

根据扫描线强度对比图及归一化均方误差表（见表1）比较可知，经过衰减强度补偿的FDK算法图像质量都有明显的提升，其中改进倒高斯加权重建轴向衰减补偿优势更大，比倒余弦重建轴向衰减补偿误差减小了11.019%，比倒高斯重建水平向128层衰减补偿误差减小了1.291%，轴向重建质量几乎趋近于原始图像，且水平向衰减补偿较倒高斯加权也有较大的提升，接近于倒余弦加权重建，同时竖直面方向上改进倒高斯加权重建衰减补偿误差更小，而水平面上重建衰减补偿误差接近于倒余弦加权重建。因此，改进倒高斯加权重建是一种更有效的衰减补偿重建方法。

表1 不同扫描位置下的归一化均方误差比较Tab.1 Comparison of normalized mean square error at different scanning positions

3 短扫描FDK算法实施流程及实验结果对比

利用改进倒高斯加权公式，本研究针对原有基于RTK的短扫描加权FDK重建方法进行了改进。首先，原有方法在开始重建之前即读入所有投影数据，无法满足实时性，改进方法通过同步患者采集和重建过程，实现了实时重建；其次，在改进方法反投影的过程中进行改进倒高斯加权，改善图像质量。基于以上两点，做了5组性能对比实验，分别是本地Parker加权重建（原有方法）、未加权实时重建、倒高斯加权实时重建、倒余弦加权实时重建、改进倒高斯加权实时重建（改进方法），滤波方式采用Hamming滤波，截止频率为0.5。其中本地Parker加权重建过程模拟基于RTK一次读入所有投影数据方法。

3.1 算法实施流程图

具体实施流程如图8所示。

图8 算法实施流程Fig.8 Algorithm implementation flowchart

3.1.1 本地Parker加权重建流程

（1）等待探测器采集投影图像时间26.4 s；

（2）一次性加载投影数据和几何位置信息；

（3）经过加权、滤波、Parker乘积遍历加权、反投影进行FDK重建；

（4）重建结果输出，计算重建时间。

3.1.2 未加权实时重建流程

（1）一次性加载几何位置信息；

（2）模拟患者采集过程，每间隔70 ms将投影数据发送给重建线程；

（3）经过加权、滤波、反投影进行FDK重建；

（4）重建结果输出，计算重建时间。

3.1.3 衰减补偿实时加权实时重建流程

（1）一次性加载几何位置信息；

（2）模拟患者采集，每间隔70 ms将投影数据发送给重建线程；

（3）经过加权、滤波、Parker乘积遍历加权、衰减补偿函数加权反投影进行FDK重建；

（4）重建结果输出，计算时间。

3.2 实验对比

本次测试以CatPhan500性能测试体模（见图9）代替患者，同时结合投影图像实例（见图10）进一步说明重建对比结果。

图9 CatPhan500性能测试体模Fig.9 CatPhan500 performance test body mask

3.2.1 测试采用的设备性能

所使用的实验平台为Intel Core i5-8400 CPU，内存8 GB，GPU为NVIDIA GeForce GTX 1050 Ti，设备上可用的全局内存总量为4 096 MB。

3.2.2 测试所用图像重建参数

射线源到重建物体中心距离为1 967 mm，射线源到探测器中心距离2 967 mm，锥角8.32o。投影图像角度从0.5o到189o，间隔0.5o取一帧，其中一帧投影图像见图10，每帧像素个数为1 440×1 440，像素大小为0.3 mm×0.3 mm，每个像素占2个字节。重建图像像素个数为512×150×512，像素大小为0.5 mm×2 mm×0.5 mm，每个像素占2个字节。测试数据采用短扫描数据，重建结果经过射野外视场消除、线性变换从图像像素值转换为CT值。

图10 投影图像示例Fig.10 Example of projected image

3.2.3 重建结果对比

重建结果对比如图11所示，可以看出，相比较未加权实时重建结果，加权重建图像变得更加清晰，图像伪影明显变少，成像质量提高了许多，这是由于重建方式为短扫描重建，且有冗余数据，当不进行Parker加权时，有冗余数据的探测器角度会多次反投影，而没有冗余数据的探测器角度只会进行一次反投影，造成重建图像质量明显不一。此外，图层2在不同重建方法中都可以看到比较明显的伪影，初步考虑为由于散射造成的边缘模糊现象。锥形束CT成像中，X射线一次与物体接触的体积更大，散射发生概率也会更大，因此可以考虑散射矫正进一步改善重建图像质量。

图11 短扫描未加权实时重建、本地Parker加权重建、倒高斯加权实时重建、倒余弦加权实时重建、改进倒高斯加权实时重建结果对比 Fig.11 Results comparison of short scan unweighted real-time reconstruction,local Parker weighted reconstruction,Gaussian weighted real-time reconstruction,inverse cosine weighted real-time reconstruction,improved inverse Gaussian weighted real-time reconstruction

从图11及表2可以看出，5种重建方式中实时重建耗费时间减小接近一半，更好地满足临床上对时效性的要求。本地Parker加权重建方法采用一次性加载所有投影图像数据，缺点是不能实时重建，无法满足实际需求。未加权实时重建方法重建过程中不计算Parker权重，所以重建时间最短。改进倒高斯加权实时重建实时同步计算Parker加权因子并作用于投影图像数据，只需要一次遍历，且在反投影过程中实施改进倒高斯加权，Parker加权实时重建方法重建图像质量更佳，时间更短。通过模拟患者扫描时间（26.4 s）和改进倒高斯加权实时重建方法重建时间（27.932 s）可以看出，改进倒高斯加权实时重建方法在患者采集完成后1～2 s内即可完成整个重建过程，达到了实时重建、改进重建图像质量的目的。

表2 五种重建方法耗费时间对比Tab.2 Time consumption comparison of five reconstruction methods

4 结论

通过同步患者采集和三维图像重建过程，在载入每帧投影数据的同时进行重建，且在反投影过程中进行新的衰减补偿加权，克服了基于RTK库的短扫描加权FDK算法一次载入所有投影图像无法满足实时重建的问题，同时能有效改善图像质量，在患者采集完成后1～2 s内完成重建，实现了短扫描Parker加权FDK算法实时重建。针对基于短扫描的Parker加权FDK算法实时重建问题进行了研究，但仍有许多不足的地方，需要进一步考虑三维图像质量的改进方法，对图像进行散射矫正以及设计滤波函数改善图像质量，更好满足临床需求。