肺癌以其高速增长的发病率和死亡率成为影响人类健康的头号杀手^[1-2], 图像引导放疗(IGRT)技术是肺癌的常规治疗手段之一, 用于各个时期的肺癌治疗。锥形束CT (CBCT)是目前最常用于图像引导放射治疗的成像技术^[2-3]。但在胸腹部CBCT的X线容积成像中, 呼吸运动所导致的模糊影像被噪声和伪影污染, 无法实现精准定位功能^[3]。4D-CBCT实现了平静呼吸运动状态下对肿瘤的实时定位, 消除了运动组织的成像模糊^[4-5], 可以通过对靶区的运动轨迹追踪实现放疗的实时监测^[6-7], 但在最基础的4D-CBCT算法中, 每一个呼吸相位的重建数据大大减少^[8-9], 降采样所产生的条状伪影和噪声严重影响了重建图像质量。为了克服降采样CBCT重建所产生的条状伪影和大量噪声, 各种4D-CBCT高分辨率重建算法被相继提出^[9-15]。先验图像约束的压缩感知方法(PICCS)将全采样重建图像作为先验图像^[16], 结合CS理论对4D-CBCT进行重建^{[12, 17]}, 提高了重建图像分辨率, 但全采样先验图的运动伪影在重建过程中也被带入到了待重建图像, 导致重建图像存在轻微的运动模糊。

为了进一步消除先验图像运动伪影对重建图像的影响, 本文提出了运动配准先验图像的4D-CBCT优质重建算法, 充分利用单一相位的位置信息, 并针对不同相位的图像重建更加准确地构建其对应的先验图像, 经过运动配准算法处理之后的先验图像可以消除运动导致的图像模糊, 在增加图像重建信息的同时防止了运动伪影的带入。

1 方法 1.1 基于先验图像约束的压缩感知算法(PICCS)的CBCT重建算法

2008年, Chen等^[18]提出了先验图像压缩感知PICCS算法, 该算法基于CS重建理论基础, 采用FBP重建算法所得到的降噪图像作为先验信息对即将重建的图像进行约束^[19-20]。PICCS作为最初的CS算法的延伸, 引入了先验图像, 先验图像I_p由全部投影数据直接重建所得, 包含了投影数据的完整信息, 先验图像的引入代表着重建信息量大幅提高, 但全部投影数据忽略运动信息直接重建的先验图像中存在大量的运动伪影和噪声。在相位为i的CBCT图像重建过程中, PICCS算法同样通过目标函数的最小化获取最优解, 算法目标函数如式(1)表示, 分为两部分, 第一部分将(I_i-I_p)引入最优化条件中, 旨在利用先验图像的信息, 第二部分I_i则调节了先验图像中伪影噪声对目标图像的影响, 两部分采用同一稀疏变化基Ψ进行稀疏化处理后分别为‖Ψ(I_i-I_p)‖₁和‖ΨI_i‖₁, 并引入权重因子α和(1-α)对约束条件的两个组成部分进行权重调节。

$ \begin{array}{l} \min \left( {\alpha {{\left\| {\Psi \left( {{I_i} - {I_p}} \right)} \right\|}_1} + \left( {1 - \alpha } \right){{\left\| {\Psi {I_i}} \right\|}_1}} \right), \\ \;\;\;\;\;\;\;\;\;s.t.\;\;\;\;\;{\Phi _i}{I_i} = {P_i}\;\;\;\;\left( {i = 1, 2, \cdots N} \right) \end{array} $

(1)

式中0 ≤ α ≤ 1, 代表先验图像对目标重建图像的影响因子, 一般依赖于先验图像的质量进行选择, 当α=1时, 待重建图像和先验图像无限接近, 当α=0时, PICCS算法退化为原始的CS算法, 式中的第一项作为保真项将投影数据信息用于待重建图像, 当对呼吸运动导致的形变组织进行采样时, 投影数据直接重建的先验图像I_p会存在一定的运动伪影, 第二项则可以作为补偿项, 降低运动伪影和噪声对待重建图像的影响, 在实际重建过程中, 一般通过α的调节来实现图像分辨率和伪影的折中。

1.2 运动配准先验图像的PICCS算法在CBCT重建中的应用

采用运动矫正后的重建图像作为先验图像, 即以重建时相i为基准, 对其余相位投影数据的重建图像进行运动矫正, 采用原i相位重建图像和其余各相位矫正后的图像进行叠加平均处理后的图像I_pi作为相位i的先验图像, 将传统的PICCS最优化问题转换为下式(2)的求解:

$ \begin{array}{l} \min \left( {\alpha {{\left\| {\Psi \left( {{I_i} - {I_{pi}}} \right)} \right\|}_1} + \left( {1 - \alpha } \right){{\left\| {\Psi {I_i}} \right\|}_1}} \right), \\ \;s.t.\;\;\;\;\;{\Phi _i}{I_i} = {P_i}\;\;\;\;\left( {i = 1, 2, \cdots N} \right) \end{array} $

(2)

其中

$ {I_{pi}} = \frac{1}{N}\sum\limits_{k = 1}^N {{M_{i, k}}} {H_{pk}} $

(3)

公式(3)中i表征一个常量, 代表待重建相位编号, k则为变量, 代表了所有相位中的任一编号, 公式(3)表征相位i的先验图像I_pi和每个相位的降采样重建图像H_pk之间的相关特性, M_{i, k}表示相位k的图像到相位i图像的运动补偿算子, 可以通过两个相位图像间的变形场得到, 由此可推论相位i到相位k运动补偿算子为M_{i, k}^-1, 当k=i时, M_{i, k}为单位矩阵, N为一个运动周期包含的相位数。在本文方法中, 对相位i进行重建时, 任何除了相位i以外的FDK降采样重建图像均需要做运动矫正后才能用于构建先验图像, 结合相位图像之间的相关性, 公式(2)可以等效为下式:

$ \begin{array}{l} \min \left( {\alpha {{\left\| {\Psi \left( {{I_i} - \frac{1}{N}\sum\limits_{k = 1}^N {{M_{i, k}}} {H_{pk}}} \right)} \right\|}_1} + \left( {1 - \alpha } \right){{\left\| {\Psi {I_i}} \right\|}_1}} \right), \\ \;s.t.\;\;\;\;\;{\Phi _i}{I_i} = {P_i}\;\;\;\;\left( {i = 1, 2, \cdots N} \right) \end{array} $

(4)

式(4)表示改进先验图像的PICCS算法, 式中I_i为相位i的待重建图像, k表示以一定时序排列的全部相位编号。本文在改进先验图像的PICCS算法实现中, 稀疏变化采用各向同性全变分构建正则化, TV正则化项可以通过以下公式计算:

$ {\left\| {{I_i}} \right\|_{TV}} = \sum\limits_j {\sqrt {\left[ {{D_x}{I_i}} \right]_j^2 + \left[ {{D_y}{I_i}} \right]_j^2 + \left[ {{D_z}{I_i}} \right]_j^2} } $

(5)

公式(5)中的D_x、D_y、D_z分别表示时域内的沿着x轴(横向)、y轴(纵向)、z轴方向的前向差分算子^[5], 其定义如下式(6):

$ \begin{array}{l} {D_x}{I_i} = vec({I_i}(x + {\rm{ }}1, y, z){\rm{ }} - {I_i}(x, y, z))\\ {D_y}{I_i} = vec({I_i}(x, y + {\rm{ }}1, z){\rm{ }} - {I_i}(x, y, z))\\ {D_z}{I_i} = vec({I_i}(x, y, z + {\rm{ }}1){\rm{ }} - {I_i}(x, y, z)) \end{array} $

(6)

改进先验图像的PICCS算法具体实现步骤如下:

(1) 获取每个相位的投影数据(降采样投影数据), 并根据成像几何分别构建系统矩阵Φ_k;

(2) 采用FDK算法通过每个相位的投影数据P_k进行单一相位的图像重建, 并利用相位之间的相关性, 通过变形场计算相位i与其余各个相位图像的运动补偿算子M_{i, k};

(3) 采用各向同性全变分正则化进行稀疏转化, 并确定比例权重因子α, 以此确立改进先验图像的PICCS算法代价函数(公式(4));

(4) 根据步骤3中确定的代价函数, 即以公式(4)中的$\min \left( {\alpha {{\left\| {\Psi \left( {{I_i} - \frac{1}{N}\sum\limits_{k = 1}^N {{M_{i, k}}} {H_{pk}}} \right)} \right\|}_1} + \left( {1 - \alpha } \right){{\left\| {\Psi {I_i}} \right\|}_1}} \right)$为约束条件, 结合正则化过程, 利用迭代算法通过求解目标函数Φ_iI_i=P_i (i=1, 2, ⋯N)优化重建图像。

1.3 实验设置

第1个实验采用肺部仿真数据, 机架旋转一周获取300个角度的投影数据, 数据采集期间组织进行有规律的模拟呼吸周期运动, X线源到探测器距离为750 mm, 旋转中心到探测器距离为250 mm, 原始投影数据大小为420×270, 像素大小为1 mm×1 mm, 在4D-CBCT重建过程中, 我们将呼吸周期分为10个相位, 因此每个相位对应30个投影数据。

第2个实验采用扫描病人数据, 整套数据采集了620个角度的投影, 并以呼吸阶段为标准将投影数据不等量地分配到10个相位组, 投影数据大小为512×512, 像素大小为0.8 mm×0.8 mm, X线源和探测器距离1536 mm, X线源和旋转中心距离1000 mm, 重建图像像素大小为1 mm×1 mm。

2 结果 2.1 肺部仿真数据结果

文章在实现改进的PICCS算法之前, 我们用传统的FDK算法分别对全部肺部仿真投影数据进行重建, 重建结果为256×256×150的3D图像, 对于三维图像的展示, 文章采用三个维度的中间切面图呈现, 所有256× 256×150的3D图像在文中均由3个维度上分别为大小256×256的二维横断面和另两个维度方向上大小为256×150的纵断面三幅图展示。

图 1展示了所有数据直接用FDK重建的结果。图 1 (A)、(B)、(C)为3D重建结果的3个截面效果图, 图 1中几乎没有条形伪影的污染, 但运动伪影使组织边界模糊。

传统的PICCS算法将全部投影数据的重建图像作为先验图像, 使运动伪影带入了待重建图像中, 文章将所有相位的投影数据分别进行FDK重建, 再以目标相位为基准, 将其余相位的图像进行运动矫正后, 用公式(3)得到目标相位的先验图像, 图 2展示了改进的PICCS算法第1个相位的先验图像的3个维度切面图, 图像组织轮廓清晰, 直观地表明了改进的先验图像和传统的PICCS算法先验图像相比, 可以去除运动伪影对待重建图像的影响, 虽然图中存在少量的条形伪影, 但图像边缘清晰度优于直接三维重建的FDK图像。

为了展示文章方法的优越性, 文章将传统PICCS算法和改进先验图像的PICCS算法进行了对比实验, 分别引入了对应的先验图像, 采用公式(4)进行迭代重建, 先验图像影响因子α=0.5, 迭代步数为80, 其余各参数均保持一致, 为了方便直观对比, 图 3呈现了先验图像和重建结果的冠状切面图, 图 3 (A)为传统PICCS算法的先验图像, 从图中可以看出由于投影数据的不一致导致的图像运动模糊, 图 3 (B)展示的传统PICCS重建结果也存在着运动伪影, 图 3 (C)比3 (A)图像组织边界清晰, 消除了由于投影数据不一致导致的运动伪影, 图 3(D)展现的组织清晰度远高于传统PICCS重建图像。

PICCS算法利用先验图像的信息提升重建图像质量, 文章方法相对传统的PICCS算法, 在保证先验图像信息的前提下进一步消除了运动伪影对重建结果的影响, 对运动的组织实现内部器官的精确定位。

为了对重建图像进行定量评估, 本文选取图 4所示的两种算法重建图像中红色线所代表的相同位置像素为感兴趣区域, 分别对其进行梯度计算后采用高斯拟合, 得到拟合曲线的FWHM(Full Width at Half Maximum), 从图 4中的拟合曲线和FWHM值可以看出, 改进先验图像的PICCS算法和传统的PICCS算法相比, 图像失真现象明显得到了改善, 图像分辨率提高。

2.2 病人数据结果

除了上文的肺部仿真实验, 本文还利用了临床病人投影数据进行重建算法对比, 图 5为传统PICCS算法和改进先验图像的PICCS算法分别采用的先验图像及其对应的重建结果在同一位置纵断面上的效果图, 除了先验图像不同, 两种结果所采用的重建算法参数均相等。仿真实验和病人数据实验均表明:相比传统PICCS重建算法, 配准后先验图像的PICCS算法所重建的CBCT图像分辨率上升, 细节对比度和图像质量均有所提高。

3 讨论

CT重建图像常因患者心跳、胃肠蠕动等生理性运动干扰, 使被扫层平面与实际扫描平面存在一定的偏离, 在重建的图像上产生伪影^[21]。CT图像的伪影受患者运动的严重程度及其方向密切相关, 通常呈现条状伪影或阴影。伪影的存在将影响临床医生及物理师对病变的范围和程度估计, 从而影响肿瘤靶区、临床靶区和计划靶区的勾画, 导致放疗计划的范围过大, 为患者带来不必要的损伤或范围过小, 导致放疗不彻底。消除或降低CT图像中的伪影是任何CT检测都要面临的问题。临床上可以通过指导受检者采用合适的呼吸方式配合完成检查、注射镇静剂、提高扫描速度、调整检查部位与球管的相对位置或使用特殊的重建技术(如运动伪影校正算法等)等使运动伪影最小化^[22-24]。对于图 1中所有数据直接采取FDK重建后得到图像, 虽然几乎消除了图像条形伪影的污染, 但由于投影数据不一致导致的运动伪影使组织边界模糊, 明显影响图像的清晰度。

4D-CBCT重建技术可以实现平静呼吸下组织动态成像, 在疾病监测和图像引导放疗等多个领域有着重大的临床应用意义^[5-7], 如何提高4D-CBCT重建图像的精确度是研究的重点^[9-15]。Sidky等^[11]提出了一种自适应梯度下降凸集投影算法(ASD-POCS), 抑制了4DCBCT每个相位稀疏采样所导致的重建图像中的条形伪影, 但此算法也导致了图像出现过平滑现象。Chen等^[18]和Chen等^[19]提出的PICCS算法在基于CS理论的所有4D-CBCT重建算法中表现上佳, 其原理为将全采样重建图像作为先验图像, 结合CS理论优化迭代公式^[17], 对4D-CBCT进行重建, 重建过程中通过先验图像的相关参数调节实现高分辨率重建图像, 而PICCS算法所有相位重建中均采用同一个先验图, 从文中图 1可以看出利用全采样数据重建的先验图像包含了运动伪影, 该现象会直接导致在4D-CBCT不同相位图像重建过程中的运动模糊被带入到待重建图像, 通过图 3 (B)和5 (B)观察到采用传统的PICCS算法的先验图像所得到的重建图像存在轻微运动模糊。

PICCS重建算法以CS理论为基础, 同样利用了数据的稀疏特性结合图像重建, 但PICCS算法引入了先验图像, 提高了单一相位CBCT图像重建的采样数据信息, 并在重建过程中通过参数调节降低先验图像运动伪影和噪声对重建结果的影响, 当先验图像不存在严重运动伪影时PICCS算法可以得到很好的重建结果, 但此方法高度依赖于先验图像质量, 当先验图像伪影和噪声太过严重时, 为了保证重建图像质量, 参数α应选择较小值, 当α的值过小, PICCS算法几乎与CS算法效果等同, 对图像质量的提升并没有实质的意义, 针对以上不足, 本文提出新的方法对先验图像进行重建, 提升质量的先验图像可以将所有采样信息用于单相位CBCT图像重建, 同时也避免了给待重建图像带来运动伪影和噪声。

为了克服传统PICCS算法先验图像的运动伪影对重建结果的影响^[25], 本研究充分利用各呼吸相位之间的相关性, 首先将所有时相的降采样重建结果以某一呼吸时相为基准对其余相位的重建图像进行运动矫正, 然后计算原相位图像和矫正后图像叠加平均所得的图像并以此作为此呼吸相位重建先验图, 最后引入每个相位对应的先验图像利用PICCS算法中的数据保真模型进行迭代重建。本研究所使用方法在传统PICCS算法基础上针对先验图像进行改进, 按照每一个相位的位置信息构建其对应的先验图像, 经过运动配准算法处理先验图像。图 3和图 5中所展示的图像重建结果直观地表明了本方法可以消除投影数据不匹配导致的重建图像运动模糊, 而图 4所展示的评估指标则进一步说明了本研究方法所重建的图像相比PICCS算法所得图像有着更优的图像分辨率, 且本算法的数据保真模型减少了重建图像失真。文章用仿真实验结果和病人数据实验结果对比展示了通过运动配准先验图像改进的PICCS算法的优越性, 文中不同数据所得到的实验结果均表明了利用运动配准先验图像所得到重建结果相比传统的PICCS重建结果组织边界清晰, 运动伪影减少, 图像分辨率提高。本文所提出的改进先验图像的PICCS算法进一步消除了运动伪影, 所得的重建图像质量优于传统的PICCS算法, 在4D-CBCT优质重建研究方向上迈进了一步, 为4D-CBCT图像重建质量的进一步提升奠定了基础。