技术领域
[0001] 本
发明涉及计算成像技术领域，尤其涉及一种基于运动感知图像约束的 4D-CBCT重建方法。
背景技术
[0002] CBCT(Cone beam computed tomography)或3D-CBCT重建是一种三维断 层图像重建技术，而四维锥束
断层图像(4D-CBCT)重建是CBCT重建技术的 拓展，利用围绕成像目标旋转采集的锥束X光透视图像重建出不同时刻下的多 组(通常10组左右)3D-CBCT。
[0003] 为了得到比较理想的重建效果，现有的4D-CBCT成像设备采用延长采集 时间的方式(比如通过旋转多周或者降低旋转速度)来获取足够的透视图像进 行4D-CBCT重建。比如，医科达公司的Versa HD
加速器采集二千多幅透视图 像进行
三维重建。采用这种方式，虽然能够重建出实际可用的4D-CBCT图像， 但会大大增加X光
辐射剂量，限制其在实际中的应用，并且，相比于3D-CBCT 图像，其图像
质量仍然较差。
[0004] 由于上述采集策略的效果较差，目前4D-CBCT的研究热点集中在
算法方 面，
力求能够借助先进的算法利用较少的透视图像重建出高质量的4D-CBCT 图像。其中，基于时域非局部均值图像约束的重建方法是一种有效方法。借助 相邻时相3D-CBCT之间具有一定相似性这一先验知识，该方法通过在
迭代重 建过程中施加时域非局部均值图像约束，能够一定程度上提升4D-CBCT重建 的质量。然而，该方法没有考虑相邻时相之间的组织运动。对于成像目标运动 速度比较快的时相，相邻时
相图像间相同
位置处的组织形态差异比较大，影响 上述方法的效果。虽然可以通过扩大搜索邻域范围来解决这个问题，但会显著 增加计算量。
发明内容
[0005] 有鉴于此，本发明提供了一种基于运动感知图像约束的4D-CBCT重建方 法，用以解决对于成像目标运动速度比较快的时相，相邻时相图像间相同位置 处的组织形态差异比较大从而影响重建效果的问题。
[0006] 因此，本发明提供了一种基于运动感知图像约束的4D-CBCT重建方法， 包括如下步骤：
[0007] S1：利用分相后的X光透视图像，计算出各个时相下的初始重建图像 {fi0，i＝1,2,…,Np}和各个时相重建图像之间的形变场 作为初值；
[0008] S2：构建基于运动感知时域非局部均值图像约束的4D-CBCT重建优化数 学模型：
[0009]
[0010] 其中，{fi，i＝1,2,…,Np}为各个时相下待重建的3D-CBCT；Np为待重建的 时相个数；Pi为每个时相X光透视图像yi对应的成像系统矩阵； 为 数据一致性项，描述重建图像与透视图像之间的一致性；u为数据一致性项的 权重；为运动感知时域非局部均值图像约束项，用于提升 相邻时相3D-CBCT之间的相似性；Ti-1和Ti+1分别表示对相邻两个时相图像fi-1和 fi+1进行形变变换，从图像fi指向相邻两个时相图像fi-1和fi+1的形变变换的参数 分别Di-1,i和Di+1,i；S(Ti-1(fi-1,Di-1,i),fi)+S(Ti+1(fi+1,Di+1,i),fi)为
图像配准相似度测量项， S(Ti-1(fi-1,Di-1,i),fi)描述由相邻时相图像fi-1经形变变换后得到的Ti-1(fi-1,Di-1,i)与当 前时相图像fi之间的相似度测量，S(Ti+1(fi+1,Di+1,i),fi)描述由相邻时相图像fi+1经 形变变换后得到的Ti+1(fi+1,Di+1,i)与当前时相图像fi之间的相似度测量；
[0011] S3：利用交替方向乘子法算法将4D-CBCT重建优化数学模型分解成以下 三个子问题：
[0012]
[0013]
[0014] P(3):{Di,i-1,Di,i+1}＝argminS(Ti-1(fi-1,Di-1,i),fi)+S(Ti+1(fi+1,Di+1,i),fi) (4)[0015] 其中，gi为辅助变量，表示重建过程中针对各个时相的辅助重建结果； 为第k步迭代得到的辅助重建结果；k为迭代次数；
[0016] S4：利用迭代优化方法，对三个子问题P(1)、P(2)和P(3)进行求解。
[0017] 在一种可能的实现方式中，在本发明提供的上述4D-CBCT重建方法中， 步骤S4，利用迭代优化方法，对三个子问题P(1)、P(2)和P(3)进行求解，具体 包括如下步骤：
[0018] S41：利用迭代优化方法求解子问题P(1)，通过一步或多步迭代分别更新 [0019]
[0020] 其中，λ为更新步长；fik-1为第k-1步迭代得到的重建图像；PiT为Pi的转 置，PiT(Pifik-1-yi)表示对投影图像的差异进行反投影，得到三维误差图；
[0021] S42：利用当前形变参数Di,i-1和Di,i+1，对fi-1和fi+1进行形变变换，得到形变 后的图像 和
[0022] S43：根据 和 求解子问题P(2)，通过一步迭代得到更新后的fik：
[0023]
[0024] 其中， 和fik(x)分别为第k步迭代各时相下的辅助重建图像和重建图 像，x和y分别为当前时相和相邻时 相图像中的
像素位置； 和 为运动感知时域非局部均值图像约束的更新项，系数 和
根据如下高斯函数来计算：
[0025]
[0026]
[0027] 其中，Z为归一化系数，使得 h为第一预设值，控制高斯 函数的宽度；
[0028] S44：通过线性整流函数对fik的像素值进行处理，确保fik的非负性；
[0029] S4 5 ：计 算当 前重 建图 像与 上一 步迭 代重 建图 像之 间 变化 量 [0030] S46：判断delta是否小于预设的
阈值；若是，则执行步骤S47；若否，则 执行步骤S48；
[0031] S47：结束迭代优化，得到最终的重建图像；
[0032] S48：利用迭代优化方法，通过m步迭代更新形变参数Di,i-1和Di,i+1，其中， m为第二预设值或者根据迭代情况自动选择的值，并返回步骤S41，继续迭代 优化。
[0033] 本发明提供的上述4D-CBCT重建方法，采用运动感知时域非局部均值图 像约束，在现有时域非局部均值图像约束的
基础上考虑不同重建时相3D-CT 之间的运动信息，重建出更高质量的4D-CBCT。通过在重建过程中同步对不 同时相下的3D-CBCT进行配准，即采用边重建边配准的机制，在重建过程中 实时估计不同时相图像之间的空间位置对应关系，在优化求解时域非局部均值 图像约束时，对估计得到的运动信息进行补偿，利用形变变换后、组织形态更 一致的相邻时相下3D-CBCT图像对特定时相下3D-CBCT图像施加非局部均 值图像约束，可以避免时相间组织形态差异引起的误差，从而可以更合理、有 效的利用相邻时相图像之间的相似性，提升重建图像的质量。
附图说明
[0034] 图1为本发明
实施例提供的基于运动感知图像约束的4D-CBCT重建方法 的
流程图之一；
[0035] 图2为本发明实施例提供的基于运动感知图像约束的4D-CBCT重建方法 的流程图之二。
具体实施方式
[0036] 下面将结合本
申请实施方式中的附图，对本申请实施方式中的技术方案进 行清楚、完整的描述，显然，所描述的实施方式仅仅是作为例示，并非用于限 制本申请。
[0037] 本发明实施例提供的一种基于运动感知图像约束的4D-CBCT重建方法， 如图1所示，包括如下步骤：
[0038] S1：利用分相后的X光透视图像，计算出各个时相下的初始重建图像 {fi0，i＝1,2,…,Np}和各个时相重建图像之间的形变场 作为初值；
[0039] 具体地，收集成像设备采集的X光透视
图像序列，根据X光透视图像分 析，将收集的X光透视图像分成Np个时相；根据每个时相下包含的X光透视 图像，利用现有重建方法重建出每个时相下的3D-CBCT{fi0，i＝1,2,…,Np}，{fi0} 作为初值用于后续的
迭代重建；利用形变配准方法(比如基于B样条形变模型 和图像互信息相似性准则的形变配准方法)，对{fi0，i＝1,2,…,Np}进行配准，得 到每个时相与相邻两个时相之间的形变变换作为初值用 于后续的迭代配准；
[0040] S2：构建基于运动感知时域非局部均值图像约束的4D-CBCT重建优化数 学模型：
[0041]
[0042] 其中，{fi，i＝1,2,…,Np}为各个时相下待重建的3D-CBCT；Np为待重建的 时相个数；Pi为每个时相X光透视图像yi对应的成像系统矩阵； 为 数据一致性项，描述重建图像与透视图像之间的一致性；u为数据一致性项的 权重；为运动感知时域非局部均值图像约束项，用于提升 相邻时相3D-CBCT之间的相似性；Ti-1和Ti+1分别表示对相邻两个时相图像fi-1和 fi+1进行形变变换，从图像fi指向相邻两个时相图像fi-1和fi+1的形变变换的参数 分别Di-1,i和Di+1,i；S(Ti-1(fi-1,Di-1,i),fi)+S(Ti+1(fi+1,Di+1,i),fi)为图像配准相似度测量项， S(Ti-1(fi-1,Di-1,i),fi)描述由相邻时相图像fi-1经形变变换后得到的Ti-1(fi-1,Di-1,i)与当 前时相图像fi之间的相似度测量，S(Ti+1(fi+1,Di+1,i),fi)描述由相邻时相图像fi+1经 形变变换后得到的Ti+1(fi+1,Di+1,i)与当前时相图像fi之间的相似度测量；
[0043] S3：利用交替方向乘子法算法将4D-CBCT重建优化数学模型分解成以下 三个子问题：
[0044]
[0045]
[0046] P(3):{Di,i-1,Di,i+1}＝argminS(Ti-1(fi-1,Di-1,i),fi)+S(Ti+1(fi+1,Di+1,i),fi) (4)[0047] 其中，gi为辅助变量，表示重建过程中针对各个时相的辅助重建结果； 为第k步迭代得到的辅助重建结果；k为迭代次数；
[0048] S4：利用迭代优化方法，对三个子问题P(1)、P(2)和P(3)进行求解。
[0049] 本发明实施例提供的上述4D-CBCT重建方法，采用运动感知时域非局部 均值图像约束，在现有时域非局部均值图像约束的基础上考虑不同重建时相 3D-CT之间的运动信息，重建出更高质量的4D-CBCT。通过在重建过程中同 步对不同时相下的3D-CBCT进行配准，即采用边重建边配准的机制，在重建 过程中实时估计不同时相图像之间的空间位置对应关系，在优化求解时域非局 部均值图像约束时，对估计得到的运动信息进行补偿，利用形变变换后、组织 形态更一致的相邻时相下3D-CBCT图像对特定时相下3D-CBCT图像施加非 局部均值图像约束，可以避免时相间组织形态差异引起的误差，从而可以更合 理、有效的利用相邻时相图像之间的相似性，提升重建图像的质量。
[0050] 在具体实施时，在本发明实施例提供的上述4D-CBCT重建方法中，步骤 S4，利用迭代优化方法，对三个子问题P(1)、P(2)和P(3)进行求解，如图2所 示，具体包括如下步骤：
[0051] S41：利用迭代优化方法求解子问题P(1)，通过一步或多步迭代分别更新 [0052]
[0053] 其中，λ为更新步长；fik-1为第k-1步迭代得到的重建图像；PiT为Pi的转 置，PiT(Pifik-1-yi)表示对投影图像的差异进行反投影，得到三维误差图；
[0054] 具体地，可以利用梯度下降优化方法求解子问题P(1)，λ可以根据Wolfe条 件等非精确线搜索方法来确定；
[0055] S42：利用当前形变参数Di,i-1和Di,i+1，对fi-1和fi+1进行形变变换，得到形变 后的图像 和
[0056] S43：根据 和 求解子问题P(2)，通过一步迭代得到更新后的fik：
[0057]
[0058] 其中， 和fik(x)分别为当前迭代(即第k步迭代)各时相下的辅助重 建图像和重建图像，x和y分别为当前时相和相邻时相图像中的像素位置； 和为运动感知时域非局部均值图像约束 的更新项，系数 和
根据如下高斯函数来计算：
[0059]
[0060]
[0061] 其中，Z为归一化系数，使得 h为第一预设值，控制高斯 函数的宽度；
[0062] 具体地，迭代公式(6)可以看成利用相邻两个时相的重建图像( 和 ) 对进行加权的时域非局部均值滤波，即fik的一个像素的灰度值为 在x点 处的灰度值、和 中x邻域Px内所有像素点灰度值的加权和；
[0063] S44：通过线性整流函数对fik的像素值进行处理，确保fik的非负性；
[0064] S4 5 ：计 算当 前重 建图 像与 上一 步迭 代重 建图 像之 间 变化 量 [0065] S46：判断delta是否小于预设的阈值；若是，则执行步骤S47；若否，则 执行步骤S48；
[0066] S47：结束迭代优化，得到最终的重建图像；
[0067] S48：利用迭代优化方法，通过m步迭代更新形变参数Di,i-1和Di,i+1，其中， m为第二预设值或者根据迭代情况自动选择的值，并返回步骤S41，继续迭代 优化；
[0068] 具体地，可以利用随机梯度优化方法进行迭代优化。
[0069] 本发明实施例提供的上述4D-CBCT重建方法，采用运动感知时域非局部 均值图像约束，在现有时域非局部均值图像约束的基础上考虑不同重建时相 3D-CT之间的运动信息，重建出更高质量的4D-CBCT。通过在重建过程中同 步对不同时相下的3D-CBCT进行配准，即采用边重建边配准的机制，在重建 过程中实时估计不同时相图像之间的空间位置对应关系，在求解时域非局部均 值图像约束时，对估计得到的运动信息进行补偿，利用形变变换后、组织形态 更一致的相邻时相下3D-CBCT图像对特定时相下3D-CBCT图像施加非局部 均值图像约束，可以避免时相间组织形态差异引起的误差，从而可以更合理、 有效的利用相邻时相图像之间的相似性，提升重建图像的质量。
[0070] 显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发 明的精神和范围。这样，倘若本发明的这些
修改和变型属于本发明
权利要求及 其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
