在图像引导放射治疗（IGRT）中，三维锥束CT（CBCT）被广泛地应用于在射束照射前检查病人的位置。然而，3D-CBCT不能捕捉靶区运动，也不能反映放射治疗过程中的呼吸运动。目前，体部立体定向放射治疗(Stereotactic Body Radiotherapy, SBRT）由于其较传统的调强放射治疗（Intensity Modulated Radiation Therapy, IMRT）具有更好的治疗效果，在肺癌治疗中得到了广泛的应用。目前SBRT射束照射前仅行3D-CBCT检查病人位置。医生不能确认病人呼吸是否与定位扫描的4D-CT相匹配。这不利于GTV靶区的准确追踪。为了弥补不足，在本研究使用双边滤波开发了一种全自动滑动运动补偿4D-CBCT重建算法，该算法在运动优化过程中对DVF进行双边滤波。双边滤波以前曾用于估计4D-CT^[1]滑动运动，但在这里将它用于4D-CBCT，这是一个更具挑战性的应用场景。这是因为从4D-CBCT几何投影成像中精确地估计4D-DVF，特别是对于滑动运动提取，比4D-CT更困难。这不仅是因为获得的CBCT投影被散射线严重的污染^[2]，而且由于传统的1 min临床扫描方案，每个阶段的可用投影数是相当有限的。通过将每个目标相位的测量投影与变形相位0%的数字重建射线照相（Digital Reconstructed Radiography, DRR）相匹配来估计4D-DVF。同时，将双边滤波引入到4D-DVF估计过程中，用于滑动运动建模。该优化过程采用非线性共轭梯度优化器。

结果表明，基于双边滤波的运动建模和重建能够更好地进行滑动运动建模。对于算法验证，作者应用了一种基于心-躯干（NURBS-based Cardiac-Torso phantom, NCAT）模体的非均匀旋转B样条。随后，4个经临床伦理批件批准的患者数据被用于初步的实验性临床验证。

1 材料与方法： 1.1 基于双边滤波的4D-CBCT图像同步重建算法

首先对原同步运动补偿重建算法进行了简要的回顾。算法有2个步骤：①利用运动补偿SART（同时代数重建技术，mSART），利用所有获得的投影重建高质量的相位0%图像。在方程1）中对运动补偿SART进行了数学描述。②通过将每个相位的测量投影与变形相位0%的DRR（数字重建射线照相）匹配来估计4D-DVF%。这2个步骤是以交错的方式执行的，以允许收敛的能量函数曲线。损失函数被设计成对称形式，以确保逆一致的DVF解，见方程2）。一旦得到4D-DVF优化解，将用于变形最终迭代重建的高质量相位0%图像，得到最终的高质量4D-CBCT[见方程6）]。从数学上讲，上述步骤可以表示如下：

设 ${\mathrm{p}}^{\mathrm{t}}=({\mathrm{p}}_{1}^{\mathrm{t}},{\mathrm{p}}_{2}^{\mathrm{t}},\dots,{\mathrm{p}}_{\mathrm{I}}^{\mathrm{t}})$ 表示相位t的对数变换4D-CBCT投影（即线积分）， $ \;{\mathrm{\mu }}^{\mathrm{t}}=({\mathrm{\mu }}_{1}^{\mathrm{t}},{\mathrm{\mu }}_{2}^{\mathrm{t}},\dots,{\mathrm{\mu }}_{\mathrm{J}}^{\mathrm{t}}) $ 表示相位t图像的衰减系数，修改后的mSART由（1）给出：

$ {\mathrm{\mu }}_{\mathrm{j}}^{0,(\mathrm{k}+1)}={\mathrm{\mu }}_{\mathrm{j}}^{0,\left(\mathrm{k}\right)}+\frac{{\displaystyle\sum }_{\mathrm{t},\mathrm{n}}{\mathrm{d}}_{\mathrm{j}\mathrm{n}}^{\mathrm{t}\to 0}{\displaystyle\sum }_{\mathrm{i}}\left[{\mathrm{a}}_{\mathrm{i}\mathrm{n}}\frac{{\mathrm{p}}_{\mathrm{i}}^{\mathrm{t}}-{\displaystyle\sum }_{\mathrm{n}}{\mathrm{a}}_{\mathrm{i}\mathrm{n}}{\mathrm{\mu }}_{\mathrm{n}}^{\mathrm{t},\left(\mathrm{k}\right)}}{{\displaystyle\sum }_{\mathrm{n}=1}^{\mathrm{J}}{\mathrm{a}}_{\mathrm{i}\mathrm{n}}}\right]}{{\displaystyle\sum }_{\mathrm{i}}{\mathrm{a}}_{\mathrm{i}\mathrm{n}}} $

(1)

其中k是迭代步骤， $ \mathrm{j} $ 是相位0%的体素指数， $ \mathrm{n} $ 是相位t的体素索引。 $ {\mathrm{a}}_{\mathrm{i}\mathrm{n}} $ 是投影射线i与体素n的交点长度，这是通过射线追踪技术得到的。 $ {\mathrm{d}}_{\mathrm{j}\mathrm{n}}^{\mathrm{t}\to 0} $ 表示将相位t变形为相位0的反向DVF。首先通过TV最小化^[3]重建初始图像 $\; {\mathrm{\mu }}_{\mathrm{j}}^{0,\left(0\right)} $ 。之后不断优化DVF和 $\; {\mathrm{\mu }}_{\mathrm{j}}^{0,\left(0\right)} $ ，并通过变形高质量 $\; {\mathrm{\mu }}_{\mathrm{j}}^{0,\left(0\right)} $ 依次得到其他相位的高质量CBCT图像。

$ {\mathrm{\mu }}_{\mathrm{n}}^{\mathrm{t},\left(\mathrm{k}\right)}={\sum }_{\mathrm{j}}{\mathrm{d}}_{\mathrm{j}\mathrm{n}}^{0\to \mathrm{t}}{\mathrm{\mu }}_{\mathrm{j}}^{0,\left(\mathrm{k}\right)} $

(2)

为优化DVF，通过配准：1）目标相位的测量投影和2）变形到目标相位的参考相位图像 $ \;{\mathrm{\mu }}_{\mathrm{j}}^{0,\left(0\right)} $ 的正投获得。相应的损失函数如下^[4-5]：

$ \begin{array}{l} {\mathrm{f}}_{1}\left({\mathrm{v}}^{0\to \mathrm{t}}\right)={‖{\mathrm{p}}^{\mathrm{t}}-\mathrm{A}{\mathrm{\mu }}^{0}\left(\mathrm{x}+{\mathrm{v}}^{0\to \mathrm{t}}\right)‖}_{{\mathrm{l}}_{2}}^{2}+\mathrm{\beta }\mathrm{\phi }\left({\mathrm{v}}^{0\to \mathrm{t}}\right)\\ {\mathrm{f}}_{2}\left({\mathrm{v}}^{\mathrm{t}\to 0}\right)={‖{\mathrm{p}}^{0}-\mathrm{A}{\mathrm{\mu }}^{\mathrm{t}}\left(\mathrm{x}+{\mathrm{v}}^{\mathrm{t}\to 0}\right)‖}_{{\mathrm{l}}_{2}}^{2}+\mathrm{\beta }\mathrm{\phi }\left({\mathrm{v}}^{\mathrm{t}\to 0}\right) \\ \mathrm{s}.\mathrm{t}.{\mathrm{v}}^{0\to \mathrm{t}}\left(\mathrm{x}+{\mathrm{v}}^{\mathrm{t}\to 0}\right)+{\mathrm{v}}^{\mathrm{t}\to 0}={\mathrm{v}}^{\mathrm{t}\to 0}\left(\mathrm{x}+{\mathrm{v}}^{0\to \mathrm{t}}\right)+{\mathrm{v}}^{0\to \mathrm{t}}=0 \end{array} $

(3)

$ {\mathrm{f}}_{1} $ 和 $ {\mathrm{f}}_{2} $ 表示对称能量函数。0代表相位0%，t代表任何其他相位t。 $ \mathrm{A} $ 是投影矩阵。 $\; {\mathrm{\mu }}^{0}\left(\mathrm{x}+{\mathrm{v}}^{0\to \mathrm{t}}\right) $ 表示图像 $ \;{\mathrm{\mu }}^{0} $ 在从相位0变形到相位t的变形场 $ {\mathrm{v}}^{0\to \mathrm{t}} $ 表示每个体素的前向DVF元素，而 $ {\mathrm{v}}^{\mathrm{t}\to 0} $ 表示每个体素从相位0到t的DVF元素。 $ {‖{\mathrm{p}}^{\mathrm{t}}-\mathrm{A}{\mathrm{\mu }}^{0}\left(\mathrm{x}+{\mathrm{v}}^{0\to \mathrm{t}}\right)‖}_{{\mathrm{l}}_{2}}^{2} $ 和 $ {‖{\mathrm{p}}^{0}-\mathrm{A}{\mathrm{\mu }}^{\mathrm{t}}\left(\mathrm{x}+{\mathrm{v}}^{\mathrm{t}\to 0}\right)‖}_{{\mathrm{l}}_{2}}^{2} $ 是对称的损失函数数据保真度项。 $ \mathrm{\phi }\left({\mathrm{v}}^{0\to \mathrm{t}}\right) $ 、 $ \mathrm{\phi }\left({\mathrm{v}}^{\mathrm{t}\to 0}\right) $ 是相应的正则化项。 $ \;\mathrm{\beta } $ 控制着数据保真度项与平滑正则化项之间的平衡。假设肺仅具有简单的各向同性运动，则 $ \mathrm{\phi }\left(\mathrm{v}\right) $ 设计如下：

$ \mathrm{\phi }\left(\mathrm{v}\right)=\sum\nolimits_{\mathrm{v}\in {\mathrm{R}}^{3}}{\sum }_{\mathrm{i}=1}^{3}{\sum }_{\mathrm{j}=1}^{3}{\left(\frac{\partial {\mathrm{v}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right)}^{2} $

(4)

其中 $\dfrac{\partial {\mathrm{v}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{j}}}$ 表示每个DVF分量沿3个方向的相邻体素之间的差异。索引“i”表示沿x、y和z方向的DVF分量。索引“j”表示3个笛卡尔坐标之一；“ $ {\mathrm{v}}_{\mathrm{i}} $ ”表示沿每个直角坐标方向的DVF元素；“ $ {\mathrm{x}}_{\mathrm{j}} $ ”表示沿每个直角坐标方向的每个图像体素。如果进一步考虑器官滑动，则重新设计基于双边滤波的正则化项：

$\begin{split}& \mathrm{\phi }\left(\mathrm{v}\right)=\sum\nolimits_{\mathrm{v}\in {\mathrm{R}}^{3}}{\sum }_{\mathrm{i}=1}^{3}{\sum }_{\mathrm{j}=1,{,\mathrm{y}}_{\mathrm{j}}\in \mathrm{N}\left({\mathrm{x}}_{\mathrm{j}}\right)}^{3}\\& \left({\mathrm{G}}_{\mathrm{x}}\left({\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{\sigma }}_{\mathrm{x}}^{2}\right){\bullet \mathrm{G}}_{ \mathrm{\mu }} \left({\mathrm{\mu }}^{\mathrm{t}}\left({\mathrm{x}}_{\mathrm{j}}\right),{\mathrm{\mu }}^{\mathrm{t}}\left({\mathrm{y}}_{\mathrm{j}}\right),{\mathrm{\sigma }}_{\mathrm{\mu }}^{2}\right)\bullet \right.\\& \left.{\mathrm{G}}_{{\mathrm{v}}_{\mathrm{i}}}\left({\mathrm{v}}_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{j}}\right),{\mathrm{v}}_{\mathrm{i}}\left({\mathrm{y}}_{\mathrm{j}}\right),{\mathrm{\sigma }}_{\mathrm{v}}^{2}\right){\left(\frac{\partial {\mathrm{v}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right)}^{2}\right) \end{split}$

(5)

其中

$ {\mathrm{G}}_{\mathrm{x}}({\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{\sigma }}_{\mathrm{x}}^{2})=\mathrm{e}\mathrm{x}\mathrm{p}\Bigg(-\frac{{({\mathrm{x}}_{\mathrm{j}}-{\mathrm{y}}_{\mathrm{j}})}^{\mathrm{T}}({\mathrm{x}}_{\mathrm{j}}-{\mathrm{y}}_{\mathrm{j}})}{2{\mathrm{\sigma }}_{\mathrm{x}}^{2}}\Bigg) $

$ {\mathrm{G}}_{\mathrm{\mu }}({\mathrm{\mu }}^{\mathrm{t}}({\mathrm{x}}_{\mathrm{j}}),{\mathrm{\mu }}^{\mathrm{t}}({\mathrm{y}}_{\mathrm{j}}),{\mathrm{\sigma }}_{\mathrm{\mu }}^{2})=\mathrm{e}\mathrm{x}\mathrm{p}\Bigg(-\frac{{‖{\mathrm{\mu }}^{\mathrm{t}}({\mathrm{x}}_{\mathrm{j}})-{\mathrm{\mu }}^{\mathrm{t}}({\mathrm{y}}_{\mathrm{j}})‖}^{2}}{2{\mathrm{\sigma }}_{\mathrm{\mu }}^{2}}\Bigg) $

$\begin{split}& {\mathrm{G}}_{{\mathrm{v}}_{\mathrm{i}}}({\mathrm{v}}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{j}}),{\mathrm{v}}_{\mathrm{i}}({\mathrm{y}}_{\mathrm{j}}),{\mathrm{\sigma }}_{\mathrm{v}}^{2})=\\&\mathrm{e}\mathrm{x}\mathrm{p}\Bigg(-\frac{{({\mathrm{v}}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{j}})-{\mathrm{v}}_{\mathrm{i}}({\mathrm{y}}_{\mathrm{j}}))}^{\mathrm{T}}({\mathrm{v}}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{j}})-{\mathrm{v}}_{\mathrm{i}}({\mathrm{y}}_{\mathrm{j}}))}{2{\mathrm{\sigma }}_{\mathrm{v}}^{2}}\Bigg)\end{split} $

$ {\mathrm{G}}_{\mathrm{x}} $ 是空间域上的高斯核，方差 $ {\mathrm{\sigma }}_{\mathrm{x}}^{2} $ ； $ {\mathrm{G}}_{\mathrm{\mu }} $ 是另一个基于像素值域的高斯核，方差 $ {\mathrm{\sigma }}_{\mathrm{\mu }}^{2} $ ； $ {\mathrm{G}}_{{\mathrm{v}}_{\mathrm{i}}} $ 是具有方差 $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ 的DVF域高斯核。索引i、j、 $ {\mathrm{v}}_{\mathrm{i}} $ 和 $ {\mathrm{x}}_{\mathrm{j}} $ 与前述公式中含义相同。 $ {\mathrm{X}}_{\mathrm{j}} $ 是每个双边滤波器核中的中心体素。 $ {\mathrm{Y}}_{\mathrm{j}} $ 表示邻域体素包围 $ {\mathrm{x}}_{\mathrm{j}} $ ， $ {\mathrm{y}}_{\mathrm{j}} $ 总个数为 $ \mathrm{N} $ 。梯度算子是在围绕每个感兴趣体素的3 × 3 × 3邻域内计算。 $ {\left.\nabla \mathrm{\phi }\left(\mathrm{v}\right)\right|}_{\mathrm{v}} $ 采用非线性共轭梯度算子优化。 $ \mathrm{\phi }\left(\mathrm{v}\right) $ 的梯度计算公式如下：

$ \begin{split} \nabla \mathrm{\phi }(\mathrm{v})=&\sum\limits_{\mathrm{v}\in {\mathrm{R}}^{3}}{\sum }_{\mathrm{i}=1}^{3}{\sum }_{\mathrm{j}=1{,\mathrm{y}}_{\mathrm{j}}\in \mathrm{N}({\mathrm{x}}_{\mathrm{j}})}^{3}\\& ({\mathrm{G}}_{\mathrm{x}}({\mathrm{x}}_{\mathrm{j}},{\mathrm{y}}_{\mathrm{j}},{\mathrm{\sigma }}_{\mathrm{x}}^{2}){\bullet \mathrm{G}}_{\mathrm{\mu }}({\mathrm{\mu }}^{\mathrm{t}}({\mathrm{x}}_{\mathrm{j}}),{\mathrm{\mu }}^{\mathrm{t}}({\mathrm{y}}_{\mathrm{j}}),{\mathrm{\sigma }}_{\mathrm{\mu }}^{2})\bullet \\& {\mathrm{G}}_{\mathrm{v}}({\mathrm{v}}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{j}}),{\mathrm{v}}_{\mathrm{i}}({\mathrm{y}}_{\mathrm{j}}),{\mathrm{\sigma }}_{\mathrm{v}}^{2})({\mathrm{v}}_{\mathrm{i}}({\mathrm{x}}_{\mathrm{j}})-{\mathrm{v}}_{\mathrm{i}}({\mathrm{y}}_{\mathrm{j}}))) \end{split} $

(6)

为了加速能量函数的收敛，需要生成初始的4D-DVF来启动优化过程。测量的CBCT投影最初被分类为4D，用于初始的4D-CBCT重建。采用Demons配准算法得到各个相位和0%相位的4D-DVF初值。算法工作流图如图1。

1.2 算法验证实验设计

本文使用NCAT体模来评估基于双边滤波的滑动运动补偿4D-CBCT重建性能。模拟4D-NCAT的10个呼吸阶段，呼吸周期为4 s，最大膈肌运动沿上下方向（SI）为20 mm，最大胸前后（AP）运动为12 mm。利用10个相位的投影，每个相位20个投影角度数据，进行DVF估计和4D-CBCT重建。数字模体图像大小为256 mm × 256 mm × 150 mm，体素大小为（1 × 1 × 1） mm³。投影尺寸为300 mm × 240 mm × 20 mm，投影体素尺寸为（1 × 1 ）mm²。为了进行运动跟踪比较，从冠状视图切片中的心脏边缘提取沿z方向的4D-NCAT运动轨迹进行定量评价。

1.3 滤波子核方差的选择

$ {\mathrm{\sigma }}_{\mathrm{x}}^{2} $ 、 $ {\mathrm{\sigma }}_{\mathrm{\mu }}^{2} $ 和 $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ 的选择标准将极大地影响最终的DVF解决方案。 $ {\mathrm{\sigma }}_{\mathrm{x}}^{2} $ 、 $ {\mathrm{\sigma }}_{\mathrm{\mu }}^{2} $ 的确定相对容易。过大或过小的 $ {\mathrm{\sigma }}_{\mathrm{x}}^{2} $ （空间平滑性）将过度平滑图像内容，或阻止其充分捕获局部稀疏特征。合理的 $ {\mathrm{\sigma }}_{\mathrm{x}}^{2} $ 不应小于2 mm。通过试验，确定了 $ {\mathrm{\sigma }}_{\mathrm{x}}^{2} $ = 3 mm能得到重建图像的最佳结果。 $ {\mathrm{\sigma }}_{\mathrm{\mu }}^{2} $ 控制肺表面与胸壁界面之间的图像强度域平滑度。可将其设置为肺组织与周围胸腔组织区域像素差，以保留从胸腔边界到肺内部的自然过渡。对NCAT数据， $ {\mathrm{\sigma }}_{\mathrm{\mu }}^{2} $ = 0.03 mm⁻¹得到最好结果。关于 $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ 的确定， $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ 应大于2个点之间的DVF差，同时其点间距小于 $ {\mathrm{\sigma }}_{\mathrm{x}}^{2} $ 。然而，在胸膜腔区域， $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ 应小于两者之间的DVF像素差。为了避免运动过度分割，本文设置只有尖锐的不连续性（例如，大的滑动运动）才能被捕获。在本文以前的工作中，相关结果比较了原始同时重建方法和地面真相的结果，发现心脏边缘位置的滑动运动估计误差约为7.5 mm。因此，可认为10 mm是一个相对较大的滑动运动。在此假设下，本文测试了几个 $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ 值，确定 $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ = 2.5 mm对NCAT数据给出了合理的结果。

1.4 量化评价标准 1.4.1 肿瘤运动准确性

从重建图像和模体金标准中提取肿瘤运动轨迹。分析各自结果中提取的肿瘤位置的均方根误差（RMSE）。

$ \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{9}\times {\sum }_{\mathrm{p}\mathrm{h}=1}^{9}{\left({\mathrm{P}\mathrm{o}\mathrm{s}}_{\mathrm{p}\mathrm{h}}^{\mathrm{R}}-{\mathrm{P}\mathrm{o}\mathrm{s}}_{\mathrm{p}\mathrm{h}}^{\mathrm{T}}\right)}^{2}} $

(7)

其中 $ {\mathrm{P}\mathrm{o}\mathrm{s}}_{\mathrm{p}\mathrm{h}}^{\mathrm{R}} $ 表示第ph个相位的估计图像特征点位置； $ {\mathrm{P}\mathrm{o}\mathrm{s}}_{\mathrm{p}\mathrm{h}}^{\mathrm{T}} $ 表示相应的金标准结果。MaxE定义为从所有9个相位图像提取的肿瘤位置的最大误差。

1.4.2 Dice系数

在最后的4D重建完成后，使用Dice系数来测量分段肺边界轮廓，以查看滑动运动补偿结果是否与真实参考相比具有更多的轮廓相似性。分割是通过ITK快照软件工具进行的。设A是由有无滑动运动补偿的结果得到的等高线面积，B是由真实基准得到的等高线。Dice系数如下：

$ \mathrm{s}=\frac{2\left|\mathrm{A}\cap \mathrm{B}\right|}{\left|\mathrm{A}\right|+\left|\mathrm{B}\right|} $

(8)

在本文研究中，使用器官轮廓内的体素数作为精确区域面积的替代。

1.4.3 相对重建错误

利用滑动建模重建的4D-CBCT与地面真相/参考之间的相对误差（RE）通过定义来量化图像重建精度

$ \mathrm{R}\mathrm{E}=\sqrt{\frac{\displaystyle\sum {\left({\mathrm{u}}_{\mathrm{R}}\left(\mathrm{x}\right)-{\mathrm{u}}_{\mathrm{T}}\left(\mathrm{x}\right)\right)}^{2}}{\displaystyle\sum {\left({\mathrm{u}}_{\mathrm{R}}\left(\mathrm{x}\right)\right)}^{2}}}\times 100\mathrm{\%} $

(9)

$ {\mathrm{u}}_{\mathrm{T}}\left(\mathrm{x}\right) $ 代表模体金标准， $ {\mathrm{u}}_{\mathrm{R}}\left(\mathrm{x}\right) $ 是重建图像。

2 结　果 2.1 NCAT模体结果 2.2 NCAT模体运动轨迹结果

图2显示了从原始重建中获得的40%相位重建图像、基于分割的滑动运动补偿重建和基于双边滤波的滑动运动补偿重建。图2a显示了从原始重建中获得的重建40%相位的矢状面视图；图2b显示了从基于分割的重建中重建的同层面矢状面；图2c显示了从双边滤波重建中重建的矢状面；图2d显示了模体金标准；白色箭头标记肋骨，这在基于双边滤波重建图像和金标准中清晰可见。在基于分割的重建中，肋骨也部分可见。但在原始重建图像中几乎看不到。白色箭头清楚地标记了肋骨的比较。图2e～h是感兴趣区域（Region on Interest，ROI），其中滑动运动存在于心脏边缘和静脉部位。经过滑动运动校正，图像中静脉（黄色虚线所示）被准确重建（即静脉长度得到了更好的重建，白色箭头示）。见双边滤波和分割的滑动运动补偿重建结果。图2i～l显示了原始重建、基于分割的重建、双边滤波重建和金标准之间的肋骨位置重建差异。在图2j、图2k，肋骨顶部边缘1和2与图2l的金标准位置匹配完好，相反没有考虑滑动运动校正的原始重建结果肋骨位置与金标准之间有较明显差异。

沿z方向的4D NCAT运动轨迹是从冠状视图切片中的心脏边缘提取的（参考图2f中的虚线位置）。通过建立每个相位的均匀阈值，从ROI图像中检测虚线位置。检测到的虚线位置用于绘制运动轨迹。图3中的运动轨迹分别来自于原始重建、基于分割的滑动补偿重建、基于双边滤波的滑动运动补偿重建和金标准。该图表明，双边滤波重建和分割重建结果中提取的轨迹与金标准更加接近。本文考虑每个轨迹的运动幅度进行RMSE计算，并确定基于双边滤波的结果的轨迹RMSE和MaxE分别为0.796 mm和1.02 mm。同时，基于分割的RMSE/MaxE与基于双边滤波的结果相当接近。

2.3 Dice系数

本文从原始同时重建、基于分割和双边滤波重建中提取4D肺边界（通过ITK-SNAP软件）。将每个不同运动建模方案的NCAT幻影实验中提取的4D Dice系数汇总。基于分割的Dice系数和基于双边滤波的Dice系数都始终大于原始的同步重建。结果表明，与原始重建相比，基于分割和双边滤波的运动估计可以更准确地分割肺边界。

2.4 参数 $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ 灵敏度分析

相对重建误差与 $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ 之间的对应关系如图4。此图表明，用 $ {\mathrm{\sigma }}_{\mathrm{v}}^{2} $ = 2.5 mm可以得到最小相对重建误差。

3 讨　论

本文讨论了一种基于双边滤波的滑动运动补偿4D-CBCT图像重建方法。相对于本文提出的方法，之前的很多工作已经报道了一系列将4D-CBCT用于精确的在体运动追踪工作。归纳后大体分为以下几类：第一类是通过增加机架旋转角度或降低机架旋转速度来增加每个机架角度下所获得的投影数^[6-7]，但延长了成像时间，增加了成像的剂量^[8-9]。第二类4D-CBCT是基于全变差最小化（Total Variation Minimization, TV Minimization）的算法^[10]。但它过度平滑了微小结构，使低对比度区域的图像质量进一步变劣。第三类是基于完全数据初始化的重建，如McKinnon-Bates（MKB）算法^[11]和基于先验图像约束压缩感知（Prior Image Constrained Compressed Sensing, PICCS）的算法^[12]。然而，残差运动将从初始重建到最终图像的伪影。第四类是低秩模型^[13]和基于张量框架^[14]的重建。然而，低秩方法不能完全实现时间分化，这2种方法都缺乏临床支持的可行性检验结果。近年来，基于可变形向量场（DVF）的4D-CBCT图像重建算法在高质量的4D-CBCT重建^[15-17]中显示出优势。大多数这些方法假设肺内各处均沿发生各向同性运动，却忽略了肺的非均匀局部运动（即滑动运动）。但滑动运动广泛存在于运动器官之间的交界面，如肺和胸壁的交界面。一些研究试图通过肺边界分割来模拟滑动运动，但由于其对肺边界半自动分割的硬性要求，使其临床转化应用受到了阻碍。本文在前述工作基础上，进一步挖掘基于DVF的精细运动补偿重建方法来获得高质量4D-CBCT。并取得了满意的初步结果。

为了客观评价平行的集中运动补偿方法，本文进一步讨论基于双边滤波核肺部表面分割方法各自得到的结果比较。基于双边滤波和基于分割的相对重建误差分别是是7.3%和7.4%。然而，结果表明在某些层面还是能发现明显的重建差异。在图2k和2l中，肋骨的形状在双边滤波方法下得到了更准确的重建。图5为分别采用双边滤波和肺表面分割得到的4D-CBCT图像重建结果，显示了基于双边滤波，肋骨位置均得到了修正（肋骨1和2均能够和金标准结果相匹配）。与金标准相比，这2种算法也对静脉长度（用黄色圆圈表示）进行了校正。双边滤波得到的静脉重建明显优于通过肺部分割得到的相应结果。

图6为双边滤波滑动运动补偿和不考虑滑动运动补偿各自得到的DVF与相应图像叠加的结果对比。图6A可见，考虑了滑动运动补偿的DVF在肺和胸腔的交界处流向发生了突变。这种DVF流量差别直接导致肋骨重建图像中的位置差别。而图6B中则没有这种情况。说明双边滤波有效的保留了发生在肺表面和胸腔的滑动运动。

关于计算时间，双边滤波重建的收敛性类似于原始重建，在theDVF中总共有200次迭代。通过检测损失函数的收敛，发现200此迭代足以达到良好的收敛性。对重建一副尺寸为200 mm × 200 mm × 150 mm的图像，该算法一次DVF优化迭代的计算时间为18 s。算法通过一块中低端GPU卡加速（英伟达GeforceGTX 980）。在加速过程中算法仅对正向投影部分进行了加速。因此后期重建速度可以进一步加快。

可能的策略包括完整的GPU实现和多个GPU联合使用。同时近期有不少基于深度学习的4D-CBCT成像研究工作发表^[18-19]，基本思路是通过建立一个前期基于4D-CT的运动预测数据集，来提取各个相位DVF之间的特征标签，实现对在线治疗时刻仅基于有限个投影情况下的呼吸运动预测。但相关工作用语建立呼吸运动预测的数据无法真实反映病人治疗时刻的运动模式，因此今后将继续拓展并建立能真实反映病人治疗时刻呼吸运动数据集的构建，并开展相应的深度学习4D-CBCT预测研究。

综上，本文提出了一种基于双边滤波的全自动滑动运动补偿4D-CBCT重建方案。数字NCAT模体实验和初步临床验证表明，该方案是一种有效的同步高质量4D-CBCT图像重建算法。实验还表明，基于双边滤波的算法在4D-CBCT重建中优于基于分割的滑动运动建模算法。该算法是一种前瞻性的4D-CBCT工具，用于图像引导放射治疗中的临床转化。