CN111951346B - 一种联合运动估计与时空张量增强表示的4d-cbct重建方法 - Google Patents

一种联合运动估计与时空张量增强表示的4d-cbct重建方法 Download PDF

Info

Publication number
CN111951346B
CN111951346B CN202010828378.6A CN202010828378A CN111951346B CN 111951346 B CN111951346 B CN 111951346B CN 202010828378 A CN202010828378 A CN 202010828378A CN 111951346 B CN111951346 B CN 111951346B
Authority
CN
China
Prior art keywords
tensor
image
phase
cbct
reconstruction
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010828378.6A
Other languages
English (en)
Other versions
CN111951346A (zh
Inventor
刘进
亢艳芹
王勇
汪军
王族
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Anhui Polytechnic University
Original Assignee
Anhui Polytechnic University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Anhui Polytechnic University filed Critical Anhui Polytechnic University
Priority to CN202010828378.6A priority Critical patent/CN111951346B/zh
Publication of CN111951346A publication Critical patent/CN111951346A/zh
Application granted granted Critical
Publication of CN111951346B publication Critical patent/CN111951346B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/008—Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/246—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
    • G06T7/251—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments involving models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20081—Training; Learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30004—Biomedical image processing
    • G06T2207/30061—Lung
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/412—Dynamic
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/424—Iterative
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/428—Real-time
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Multimedia (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种联合运动估计与时空张量增强表示的4D‑CBCT重建方法，属于计算机断层成像技术领域。本发明中，首先利用呼吸信号对病人4D‑CBCT投影数据进行相位划分，并采用EM‑TV算法重建获得初始图像，利用初始图像估计初始的形变矢量场；然后，将时空张量增强表示引入到图像重建模型中，对形变矢量场处理后的重建图像进行约束，并获得更新重建图像；接下来，利用更新后的重建图像及投影数据建立对称的运动估计模型，以获得新的形变矢量场；最后，通过迭代交替更新重建图像及形变矢量场，获得高质量重建图像。本发明可有效减缓4D‑CBCT重建中的运动模糊及混叠伪影，保留更多细节，提高重建图像对比度，促进4D‑CBCT在图像引导放疗领域中的使用。

Description

一种联合运动估计与时空张量增强表示的4D-CBCT重建方法
技术领域
本发明涉及计算机断层成像技术领域，更具体地说，涉及一种联合运动估计与时空张量增强表示的4D-CBCT重建方法。
背景技术
随着医学技术的飞速发展，肿瘤放射治疗已经迈入精准时代，即精确定位、精确计划与精确放疗。而在放射治疗过程中，肿瘤与周围正常组织的摆位、位移的误差都会影响放疗精度，从而影响最终的治疗效果。在这过程中，图像质量保证尤其重要，是实现肿瘤精确定位跟踪的前提。为解决这一问题，人们将放疗设备和影像设备结合在一起，实现全自动肿瘤定位和跟踪，确定治疗靶区和重要组织结构的位置、运动轨迹，并在必要时进行姿态校正，这就是图像引导放疗(Image-Guided Radiation Therapy,IGRT)技术。IGRT技术中以锥形束CT (Cone Beam CT,CBCT)成像引导的放疗系统最具代表性，它是采用直线加速器机载影像系统架构，在放射治疗的同时实现CBCT扫描、重建、肿瘤定位及摆位误差计算等，为精准放疗提供依据。该技术可实现最大程度的肿瘤照射及射线利用率，能尽可能的避免对周围正常组织的损害，在提高放疗精度及治疗疗效方面起到重要作用。
当前，动态的4D-CBCT重建技术以逐步应用于人体呼吸运动、心脏与大血管搏动以及食管蠕动等状态成像目标的放疗过程中。但在基于时间序列的4D-CBCT图像重建中，由于扫描速度慢、投影角度少原因，会导致单个时间相位下投影数据严重缺乏，致使重建图像中含有严重的条状伪影，肿瘤边界难以确定，从而最终导致放疗计划中肿瘤靶区目标的定位不准，放疗剂量的分配存在较大误差。目前主流的放疗系统中，CBCT旋转一周的时间约为1～3 分钟，而扫描病人的呼吸周期约为3～5秒，这也就导致CBCT的照射时间分辨率远低于呼吸运动的时间分辨率，使得4D-CBCT重建图像质量下降。组织器官的运动会不可避免的造成肿瘤的漏照和/或正常组织的过多卷入与损伤，这是影响精准放疗实施的主要因素。针对4D-CBCT成像这一扫描模式的特殊性，近年来科研界和工业界开展了大量研究，旨在不过多增加扫描剂量的情况下，并通过充分利用X射线尽可能的提高图像成像质量，使其更好的服务于临床。当前解决该问题的主要技术之一就是基于运动补偿的重建。该类方法的基本策略是通过对非刚性的人体解剖组织图像配准，来估计不同时间相位下的重建图像之间的形变矢量场(Deformation Vector Field，DVF)，进而对组织器官运动信息进行补偿。
在4D-CBCT图像重建过程中，每个相位下的重建图像需要进行重投影，并对比扫描投影数据这一过程是依据呼吸信号进行分组的，而待配准图像与参考图像中均包含严重的噪声和条状伪影，因此单个相位图像的重建是一个病态问题。直接采用图像域中的形变矢量场作为重建中的运动补偿信息，难以保证重建后图像的拓扑结构准确。总的来说，该类方法重建的图像质量依然不高，临床应用程度不广，相关理论仍存在一些亟待解决的问题，如：如何提高重建过程中形变矢量场估计的准确性，如何充分利用不同相位之间的信息引入先验知识，如何在重建中恰当的将形变矢量场和先验知识相结合，并提高成像质量具有重要的临床应用价值。
公开于该背景技术部分的信息仅仅旨在增加对本发明的总体背景的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域一般技术人员所公知的现有技术。
发明内容
1.发明要解决的技术问题
本发明的目的在于克服现有技术中4D-CBCT重建方法存在的图像质量不高、组织细节丢失、对比度低等问题，拟提供一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，称之为联合运动估计与时空张量增强表示重建(Motion Estimation andSpatiotemporal Tensor Enhanced representation Reconstruction,简称MESTER)，在不改变现有的CBCT硬件成本条件下，通过在重建中结合运动估计与时空张量增强表示的方法，来抑制由呼吸运动和角度缺失引起的图像模糊和细节丢失现象，提高4D-CBCT重建图像质量，最终实现为放疗过程医生提供实时准确的影像信息，从而提高放射治疗精度和疗效。
2.技术方案
为达到上述目的，本发明提供的技术方案为：
本发明的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，包括以下步骤：
S1、获取初始重建图像及初始形变矢量场；
利用给定的呼吸信号对4D-CBCT投影数据进行相位划分，获得划分后的投影数据P＝{p1,p1,…pt,…pT}，其中T为总相位数，pt为t相位下投影数据；对划分后的投影数据进行EM-TV重建，获得初始的重建图像
Figure SMS_1
其中ft 0为t相位下的初始重建图像，上标0表示初始值；并利用初始重建图像及demons配准算法获得初始的形变矢量场
Figure SMS_2
其中t1≠t，表示t1相位到t相位下重建图像的形变矢量场；
S2、构建时空张量增强表示的4D-CBCT图像重建模型；
S3、构建4D-CBCT图像运动估计模型；
S4、求解重建模型及运动估计模型获得最终的重建图像；
对于重建模型及运动估计模型进行分别求解，并采用交替更新的方式对重建图像及形变矢量场进行迭代更新，获得最终的重建图像。
更进一步地，S2中的4D-CBCT图像重建模型为：
Figure SMS_3
Figure SMS_4
其中At为t相位下的投影矩阵，ft为t相位下的重建图像，pt为t相位下的投影数据，F表示所有相位下的重建图像，D为t1(t1≠t)相位到t相位下重建图像的形变矢量场集合，D·F 表示t相位下补偿后的重建图像集合，λ及μ均为正则化参数，
Figure SMS_5
为四维梯度算子，||·||0为 L0范数，Zc(·)和C分别表示第c个张量块组及张量块组总数，/>
Figure SMS_6
表示分解处理后的第c 个张量块组，λc为第c个张量块组的分解权重函数，Ui(1≤i≤5)表示张量分解后第i维度上的归一化向量，[[·]]表示张量分解过程。
更进一步地，S3中的4D-CBCT图像运动估计模型为：
Figure SMS_7
其中At为t相位下的投影矩阵，ft为t相位下的重建图像，pt为t相位下的投影数据，
Figure SMS_8
为t1相位图像/>
Figure SMS_9
在t相位下补偿后的重建图像，η为正则化参数，TV(·)为TV正则化约束项，用于减小噪声干扰，获得边界增强的形变矢量场。
更进一步地，S1中相位划分是以呼吸信号强度进行均匀划分的；且EM-TV重建最大迭代次数为150次。
更进一步地，S2中
Figure SMS_10
为补偿后的重建图像集合的梯度L0范数，具体表示为：
Figure SMS_11
其中
Figure SMS_12
与/>
Figure SMS_13
分别表示x,y,z和时间t方向上的梯度算子，θ(·)为计数算子，J 为总像素，D·F表示t相位下补偿后的重建图像集合。
更进一步地，S2中张量块组获得的步骤为：首先对运动补偿后的重建图像集合D·F进行分块，得到时空张量图像块集合；然后对这些张量图像块集合进行分组，最后获得分组后的张量块组，其中Zc(D·F)为D·F的第c个张量块组。
更进一步地，S2中张量图像块集合分组是按照K均值算法进行相似性分组。
更进一步地，S4中4D-CBCT图像重建模型采用正反向分裂算法求解，分别表示为：
Figure SMS_14
Figure SMS_15
Sk+1＝Sk+Gk+1-Fk+1 (3)
Figure SMS_16
其中At为t相位下的投影矩阵，ft为t相位下的重建图像，pt为t相位下的投影数据，F表示所有相位下的重建图像，则Fk为第k次迭代后有相位下的重建图像，D为t1(t1≠t)相位到t相位下重建图像的形变矢量场集合，D·F表示t相位下补偿后的重建图像集合，λ及μ均为正则化参数，
Figure SMS_17
为四维梯度算子，||·||0为L0范数，Zc(·)和C分别表示第c个张量块组及张量块组总数，/>
Figure SMS_18
表示分解处理后的第c个张量块组，λc为第c个张量块组的分解权重函数，Ui(1≤i≤5)表示张量分解后第i维度上的归一化向量，[[·]]表示张量分解过程；G及Gk分别为F的对偶变量及第k次迭代后的对偶变量，Sk为第k次迭代后的对偶辅助变量，上标 k表示第k次迭代结果，β为拉格朗日乘子参数；
执行过程中，依次求解公式(1)、公式(2)、公式(3)和公式(4)，其中公式(1)采用可分离抛物面替代算法求解，公式(2)采用L0范数梯度最小化算法求解，公式(4)采用张量CP分解方法求解。
更进一步地，S4中4D-CBCT图像运动估计模型采用交替方向乘子法求解。
更进一步地，S4中采用交替更新的方式对重建图像及形变矢量场进行迭代更新，直到达到最大迭代次数20次。
3.有益效果
采用本发明提供的技术方案，与现有技术相比，具有如下有益效果：
本发明的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，首先对给定的呼吸信号及4D-CBCT投影数据进行相位划分，获得划分后的投影数据并进行重建，得到初始重建图像，并由此估计初始形变矢量场；然后将初始重建图像及初始形变矢量场代入4D-CBCT图像重建模型及4D-CBCT图像运动估计模型中进行求解；最后通过交替更新重建图像及形变矢量场进行迭代重建。通过将运动估计模型和时空张量增强表示有机结合，可有效的改善了现有的重建方法在对呼吸运动重建时细节容易丢失及对比度下降的问题。实验结果验证了在常规呼吸状态下进行4D-CBCT重建，本发明的方法与运动补偿学习的重建方法 (Motion Compensation Learning Reconstruction,简称MCLR)相比，可够有效的抑制CBCT 图像中由投影角度缺失导致的条状伪影，减缓器官运动带来的组织模糊及细节丢失，同时软组织对比度更高，视觉效果更好。本发明的方法有望辅助医师实现肿瘤的精准定位和跟踪，为患者降低额外辐射，提高放疗收益，具有较高的应用和推广前景。
附图说明
图1为本发明中联合运动估计与时空张量增强表示的4D-CBCT重建方法流程图；
图2为本发明实施例中模拟的动态肺部数据及真实采集的肺部数据不划分相位下的横断面重建图像及使用的呼吸信号曲线(a：模拟数据；b：真实数据；c：呼吸信号曲线)；
图3为本发明实施例中模拟的动态肺部数据第5个相位下的横断面重建结果(a：参考图； b：FDK重建图像；c：MCLR重建图像；d：MESTER重建图像)；
图4为本发明实施例中模拟的动态肺部数据第5个相位下的横断面重建结果的局部放大图(a：参考图；b：FDK重建图像；c：MCLR重建图像；d：MESTER重建图像)；
图5为本发明实施例中模拟的动态肺部数据在本发明方法下所有相位的横断面重建结果 (a：参考图像；b：MESTER重建图像)；
图6为本发明实施例中真实采集的肺部数据第2个相位下的横断面重建结果(a：FDK 重建图；b：MCLR重建图像；c：MESTER重建图像)；
图7为本发明实施例中真实采集的肺部数据第2个相位下的冠状面重建结果(a：FDK 重建图；b：MCLR重建图像；c：MESTER重建图像)；
图8为本发明实施例中真实采集的肺部数据第2个相位下的冠状面重建结果的Profile曲线；
具体实施方式
为进一步了解本发明的内容，结合附图对本发明作详细描述。
在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。
下面结合实施例对本发明作进一步的描述。
实施例1
本实施例的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法流程图如图1 所示，具体步骤如下：
步骤S1、获取初始图像及初始形变矢量场；
利用给定的呼吸信号对4D-CBCT投影数据进行相位划分，可获得划分后的投影数据 P＝{p1,p1,…pt,…pT}，其中T为总相位数，pt为t相位下投影数据；对划分后的投影数据进行EM-TV重建，获得初始的重建图像
Figure SMS_19
其中ft 0为t相位下的初始重建图像，上标0表示初始值；并利用初始图像及demons配准算法获得初始的形变矢量场
Figure SMS_20
其中t1≠t，表示t1相位到t相位下重建图像的形变矢量场；
具体的，相位划分是以呼吸信号强度进行均匀划分的(如：呼吸信号最大强度为1，需重建10个相位图，则以呼吸信号强度间隔0.1为依据进行相位划分)，并对投影数据进行分组。分组后的投影数据按照EM-TV算法进行最大迭代次数为150次的迭代重建，图像更新过程可表示为：
Figure SMS_21
其中At和
Figure SMS_22
分别为t相位下的投影矩阵及反投影矩阵，ft为t相位下的重建图像，上标k 表示第k次迭代结果，pt为t相位下的投影数据，I为全1向量，γ为正则化参数，
Figure SMS_23
为对ft k求TV范数的导数。
步骤S2、构建时空张量增强表示的4D-CBCT图像重建模型；
4D-CBCT图像重建模型为：
Figure SMS_24
Figure SMS_25
其中At为t相位下的投影矩阵，ft为t相位下的重建图像，pt为t相位下的投影数据，F为所有相位下的重建图像，D为t1(t1≠t)相位到t相位下重建图像的形变矢量场集合，D·F表示t相位下补偿后的重建图像集合，λ及μ均为正则化参数，
Figure SMS_26
为四维梯度算子，||·||0为L0 范数，Zc(·)和C分别表示第c个张量块组及张量块组总数，/>
Figure SMS_27
表示分解处理后的第c个张量块组，λc为第c个张量块组的分解权重函数，Ui(1≤i≤5)表示张量分解后第i维度上的归一化向量，[[·]]表示张量分解过程；
具体的，4D-CBCT图像重建模型中
Figure SMS_28
为补偿后的重建图像集合的梯度L0范数，具体可表示为：/>
Figure SMS_29
其中/>
Figure SMS_30
Figure SMS_31
与/>
Figure SMS_32
分别表示x,y,z和时间t方向上的梯度算子，θ(·)为计数算子，J为总像素，D·F 表示t相位下补偿后的重建图像集合。张量块组获得的步骤为：首先对运动补偿后的重建图像 D·F按照一定大小及间隔进行分块(如：图块尺寸为16×16×16×4像素，分块间隔为 12×12×12×2像素)，获得时空张量图像块集合/>
Figure SMS_33
其中N为总时空张量图像块数目；然后对这含N个张量图像块的集合按照K均值算法进行相似性分组，共分为C组；最后获得分组后的张量块组Zc(D·F)。
步骤S3、构建4D-CBCT图像运动估计模型：
4D-CBCT图像运动估计模型为：
Figure SMS_34
其中At为t相位下的投影矩阵，ft为t相位下的重建图像，pt为t相位下的投影数据，
Figure SMS_35
为t1相位图像/>
Figure SMS_36
在t相位下补偿后的重建图像，η为正则化参数，TV(·)为TV正则化约束项，可减小噪声干扰，获得边界增强的形变矢量场。
步骤S4、求解重建模型及运动估计模型获得最终的重建图像；
对于重建模型及运动估计模型进行分别求解，并采用交替更新的方式对重建图像及形变矢量场进行迭代更新，获得高质量的最终的重建图像。
具体的，首先对4D-CBCT图像重建模型采用正反向分裂算法求解，获得如下子问题为：
Figure SMS_37
Figure SMS_38
Sk+1＝Sk+Gk+1-Fk+1 (3)
Figure SMS_39
其中At为t相位下的投影矩阵，ft为t相位下的重建图像，pt为t相位下的投影数据，F为所有相位下的重建图像，Fk为第k次迭代后的重建图像，D为t1(t1≠t)相位到t相位下重建图像的形变矢量场集合，D·F表示t相位下补偿后的重建图像集合，λ及μ均为正则化参数，
Figure SMS_40
为四维梯度算子，||·||0为L0范数，Zc(·)和C分别表示第c个张量块组及张量块组总数，/>
Figure SMS_41
表示分解处理后的第c个张量块组，λc为第c个张量块组的分解权重函数，Ui(1≤i≤5)表示张量分解后第i维度上的归一化向量，[[·]]表示张量分解过程；G及Gk分别为F的对偶变量及第k次迭代后的对偶变量，Sk为第k次迭代后的对偶辅助变量，上标k表示第k次迭代结果，β为拉格朗日乘子参数。在4D-CBCT图像重建模型求解过程中，依次求解上面的公式(1)、公式(2)、公式(3)和公式(4)获得重建图像，其中公式(1)采用可分离抛物面替代算法求解，公式(2)采用L0范数梯度最小化算法求解，公式(4)采用张量 CP分解方法进行求解。然后，进入4D-CBCT图像运动估计模型求解过程中，本实施例是采用交替方向乘子法对运动估计模型进行求解获得形变矢量场。最后，对重建图像及形变矢量场交替迭代更新，最大迭代20次后得到最终的重建图像。
效果评估准则
对模拟的动态肺部数据及真实采集的肺部数据进行重建，模拟扫描使用的参数为：探测器单元个数为1024×200，探测器单元尺寸为0.78×0.178mm2，射线源到物体中心和探测器中心的距离分别为50cm和100cm，扫描采集644个投影数据，其它参数采用默认值。真实数据扫描使用的参数为：探测器大小为504×504，探测器单元尺寸为0.5×0.5mm2，射线源到物体中心和探测器中心的距离分别为57cm和106.23cm，扫描采集644个投影数据，其它参数采用默认值。模拟数据采用与真实数据相同的呼吸信号。重建模型求解中正则化参数λ与μ分别为0.02和0.001，拉格朗日乘子β为600，C为100，G0＝F0，S0＝0。运动估计模型求解过程中正则化参数η为0.05。
模拟的动态肺部数据重建图像的显示窗宽为600HU(Housfield Units,HU)，窗位为50HU，真实采集的肺部数据重建图像的显示窗宽为0.025，窗位为0.0175。通过比较模拟的动肺部数据下FDK重建方法、MCLR重建方法和MESTER方法的重建结果来验证本发明的有效性(如图3和图4所示)；通过比较所有相位的横断面重建结果图像观察组织器官的运动过程及重建偏差(如图5所示)；通过对真实采集的肺部数据进行实验，进一步验证本发明方法在真实数据中的重建效果(如图6，图7和图8所示)。为客观有效的做出对比，采用视觉评估和量化评估两种方法进行评估。
视觉评估
通过观察比较图3、4、5、6、7中CBCT图像的结构细节、对比度、解剖组织纹理、噪声伪影强度等特点，可以看出本发明方法能获得更高质量的4D-CBCT重建图像。同时，从重建结果中可以看到MESTER重建方法较MCLR重建方法可以更好的保持图像细节，提高对比度。从局部放大中可以看到，直接FDK重建方法噪声伪影严重，无法区分组织细节； MCLR重建虽然能抑制部分伪影，但对比度有所下降，仍有部分噪声残留；而使用本发明方法MESTER重建后，图像质量有明显提高，组织细节更加清晰，效果更好，在模拟数据中更接近参考图。
量化评估
在使用视觉评估本发明方法的有效性的同时，引入两个量化指标对本发明方法的有效性进行评判，峰值信噪比(Peak Signal to Noise Ratio,PSNR)和结构相似度(Structural Similarity Index,SSIM)。其计算方法如下：
Figure SMS_42
Figure SMS_43
其中ft为t相位下的重建图像，rt为t相位下的CBCT参考图，N为图像像素总数；μft和μrt分别表示图ft和rt的均值；σft和σrt分别为对应的标准差，σfrt为对应的协方差；C1和 C2为两个常数，其中
Figure SMS_44
Figure SMS_45
为rt的最大值。计算出各相位下重建图像的PSNR和SSIM的均值±方差值，如表1所示。
表1
Figure SMS_46
Figure SMS_47
从表1中可以看出，在模拟的动态肺部投影数据重建中，FDK重建图像的量化指标最差， MCLR重建结果具有一定的提高，而采用本发明方法可以有更高的SSIM和PSNR(本发明 MESTER方法与MCLR方法结果相比，SSIM约高出0.02左右，PSNR约高出1dB左右)，且方差更小。从图8中可以看出，在真实采集的肺部数据重建中，不同组织边界像素值跳跃更加明显，组织边界更锐利。从上述实验可以看到，在相同投影角度扫描条件下本发明方法可获得高质量的4D-CBCT重建图像，且稳定性高，具有一定的应用前景。
以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均应属于本发明的保护范围。

Claims (9)

1.一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，其特征在于，包括以下步骤：
S1、获取初始重建图像及初始形变矢量场；
利用给定的呼吸信号对4D-CBCT投影数据进行相位划分，获得划分后的投影数据
Figure QLYQS_2
，其中T为总相位数，/>
Figure QLYQS_5
为t相位下投影数据；对划分后的投影数据进行EM-TV重建，获得初始的重建图像/>
Figure QLYQS_6
，其中/>
Figure QLYQS_3
为t相位下的初始重建图像，上标0表示初始值；并利用初始重建图像及demons配准算法获得初始的形变矢量场
Figure QLYQS_4
，其中/>
Figure QLYQS_7
，表示/>
Figure QLYQS_8
相位到/>
Figure QLYQS_1
相位下重建图像的形变矢量场；
S2、构建时空张量增强表示的4D-CBCT图像重建模型；
4D-CBCT图像重建模型为：
Figure QLYQS_9
其中
Figure QLYQS_19
为/>
Figure QLYQS_12
相位下的投影矩阵，/>
Figure QLYQS_15
为/>
Figure QLYQS_21
相位下的重建图像，/>
Figure QLYQS_25
为/>
Figure QLYQS_28
相位下的投影数据，/>
Figure QLYQS_31
表示所有相位下的重建图像，/>
Figure QLYQS_24
为/>
Figure QLYQS_26
（/>
Figure QLYQS_10
）相位到/>
Figure QLYQS_16
相位下重建图像的形变矢量场集合，
Figure QLYQS_13
表示/>
Figure QLYQS_18
相位下补偿后的重建图像集合，/>
Figure QLYQS_20
及/>
Figure QLYQS_22
均为正则化参数，/>
Figure QLYQS_27
为四维梯度算子，
Figure QLYQS_30
为L0范数，/>
Figure QLYQS_29
和C分别表示第c个张量块组及张量块组总数，/>
Figure QLYQS_32
表示分解处理后的第c个张量块组，/>
Figure QLYQS_11
为第c个张量块组的分解权重函数，/>
Figure QLYQS_14
（/>
Figure QLYQS_17
）表示张量分解后第i维度上的归一化向量，/>
Figure QLYQS_23
表示张量分解过程；
S3、构建4D-CBCT图像运动估计模型；
S4、求解重建模型及运动估计模型获得最终的重建图像；
对于重建模型及运动估计模型进行分别求解，并采用交替更新的方式对重建图像及形变矢量场进行迭代更新，获得最终的重建图像。
2.根据权利要求1所述的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，其特征在于：S3中的4D-CBCT图像运动估计模型为：
Figure QLYQS_33
其中
Figure QLYQS_35
为/>
Figure QLYQS_39
相位下的投影矩阵，/>
Figure QLYQS_42
为/>
Figure QLYQS_37
相位下的重建图像，/>
Figure QLYQS_40
为/>
Figure QLYQS_43
相位下的投影数据，
Figure QLYQS_45
为/>
Figure QLYQS_36
相位图像/>
Figure QLYQS_38
在/>
Figure QLYQS_41
相位下补偿后的重建图像，/>
Figure QLYQS_44
为正则化参数，/>
Figure QLYQS_34
为TV正则化约束项，用于减小噪声干扰，获得边界增强的形变矢量场。
3.根据权利要求1所述的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，其特征在于：S1中相位划分是以呼吸信号强度进行均匀划分的；且EM-TV重建最大迭代次数为150次。
4.根据权利要求1所述的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，其特征在于：S2中
Figure QLYQS_46
为补偿后的重建图像集合的梯度L0范数，具体表示为：
Figure QLYQS_47
其中
Figure QLYQS_48
，/>
Figure QLYQS_49
，/>
Figure QLYQS_50
与/>
Figure QLYQS_51
分别表示x, y, z和时间t方向上的梯度算子，/>
Figure QLYQS_52
为计数算子，J为总像素，/>
Figure QLYQS_53
表示/>
Figure QLYQS_54
相位下补偿后的重建图像集合。
5.根据权利要求1所述的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，其特征在于：S2中张量块组获得的步骤为：首先对运动补偿后的重建图像集合
Figure QLYQS_55
进行分块，得到时空张量图像块集合；然后对这些张量图像块集合进行分组，最后获得分组后的张量块组，其中/>
Figure QLYQS_56
为/>
Figure QLYQS_57
的第c个张量块组。
6.根据权利要求5所述的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，其特征在于：S2中张量图像块集合分组是按照K均值算法进行相似性分组。
7.根据权利要求2所述的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，其特征在于：S4中4D-CBCT图像重建模型采用正反向分裂算法求解，分别表示为：
Figure QLYQS_58
（1）
Figure QLYQS_59
（2）
Figure QLYQS_60
（3）
Figure QLYQS_61
（4）
其中
Figure QLYQS_86
为/>
Figure QLYQS_88
相位下的投影矩阵，/>
Figure QLYQS_90
为/>
Figure QLYQS_65
相位下的重建图像，/>
Figure QLYQS_77
为/>
Figure QLYQS_81
相位下的投影数据，/>
Figure QLYQS_85
表示所有相位下的重建图像，则/>
Figure QLYQS_64
为第k次迭代后有相位下的重建图像，/>
Figure QLYQS_67
为/>
Figure QLYQS_69
（/>
Figure QLYQS_71
）相位到/>
Figure QLYQS_73
相位下重建图像的形变矢量场集合，/>
Figure QLYQS_75
表示/>
Figure QLYQS_79
相位下补偿后的重建图像集合，/>
Figure QLYQS_83
及
Figure QLYQS_80
均为正则化参数，/>
Figure QLYQS_84
为四维梯度算子，/>
Figure QLYQS_87
为L0范数，/>
Figure QLYQS_89
和C分别表示第c个张量块组及张量块组总数，/>
Figure QLYQS_62
表示分解处理后的第c个张量块组，/>
Figure QLYQS_66
为第c个张量块组的分解权重函数，/>
Figure QLYQS_78
（/>
Figure QLYQS_82
）表示张量分解后第i维度上的归一化向量，/>
Figure QLYQS_63
表示张量分解过程；/>
Figure QLYQS_68
及
Figure QLYQS_70
分别为/>
Figure QLYQS_72
的对偶变量及第k次迭代后的对偶变量，/>
Figure QLYQS_74
为第k次迭代后的对偶辅助变量，上标k表示第k次迭代结果，/>
Figure QLYQS_76
为拉格朗日乘子参数；
执行过程中，依次求解上述公式（1）、公式（2）、公式（3）和公式（4），其中公式（1）采用可分离抛物面替代算法求解，公式（2）采用L0范数梯度最小化算法求解，公式（4）采用张量CP分解方法求解。
8.根据权利要求7所述的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，其特征在于：S4中4D-CBCT图像运动估计模型采用交替方向乘子法求解。
9.根据权利要求1-8任一项所述的一种联合运动估计与时空张量增强表示的4D-CBCT重建方法，其特征在于：S4中采用交替更新的方式对重建图像及形变矢量场进行迭代更新，直到达到最大迭代次数20次。
CN202010828378.6A 2020-08-18 2020-08-18 一种联合运动估计与时空张量增强表示的4d-cbct重建方法 Active CN111951346B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010828378.6A CN111951346B (zh) 2020-08-18 2020-08-18 一种联合运动估计与时空张量增强表示的4d-cbct重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010828378.6A CN111951346B (zh) 2020-08-18 2020-08-18 一种联合运动估计与时空张量增强表示的4d-cbct重建方法

Publications (2)

Publication Number Publication Date
CN111951346A CN111951346A (zh) 2020-11-17
CN111951346B true CN111951346B (zh) 2023-06-23

Family

ID=73342655

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010828378.6A Active CN111951346B (zh) 2020-08-18 2020-08-18 一种联合运动估计与时空张量增强表示的4d-cbct重建方法

Country Status (1)

Country Link
CN (1) CN111951346B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113077403B (zh) * 2021-04-22 2022-08-09 河北工业大学 基于局部数据块的张量增强技术的彩色图像重建方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105849778A (zh) * 2013-12-20 2016-08-10 皇家飞利浦有限公司 成像中的移动结构运动补偿
CN107025632A (zh) * 2017-04-13 2017-08-08 首都师范大学 一种图像超分辨率重建方法及系统
CN108898642A (zh) * 2018-06-01 2018-11-27 安徽工程大学 一种基于卷积神经网络的稀疏角度ct成像方法
CN108961237A (zh) * 2018-06-28 2018-12-07 安徽工程大学 一种基于卷积神经网络的低剂量ct图像分解方法

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2008141656A1 (en) * 2007-05-17 2008-11-27 Elekta Ab (Publ) On-line cone beam ct reconstruction
CN104268914B (zh) * 2014-10-24 2017-01-18 山东师范大学 一种4d‑ct不同时相序列图像的重建方法
US9990731B2 (en) * 2016-01-13 2018-06-05 Varian Medical Systems International Ag Systems and methods for evaluating motion tracking for radiation therapy
WO2019005180A1 (en) * 2017-06-26 2019-01-03 Elekta, Inc. METHOD FOR ENHANCING CT QUALITY QUALITY WITH CONIC BEAM THROUGH DEEP CONVOLUTIVE NEURAL NETWORK
CN109496327A (zh) * 2017-12-13 2019-03-19 上海联影医疗科技有限公司 用于诊断和治疗的系统和方法
US10672153B2 (en) * 2018-04-23 2020-06-02 Elekta Ab (Publ) Posterior image sampling using statistical learning model
CN110400357B (zh) * 2019-07-05 2020-04-14 北京航空航天大学 一种基于运动感知图像约束的4d-cbct重建方法
CN110390361B (zh) * 2019-07-25 2021-04-09 安徽工程大学 一种基于运动补偿学习的4d-cbct成像方法
CN110570489B (zh) * 2019-09-05 2021-07-30 重庆医科大学附属第一医院 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法
CN111275669B (zh) * 2020-01-13 2022-04-22 西安交通大学 一种先验信息引导的四维锥束ct图像重建算法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105849778A (zh) * 2013-12-20 2016-08-10 皇家飞利浦有限公司 成像中的移动结构运动补偿
CN107025632A (zh) * 2017-04-13 2017-08-08 首都师范大学 一种图像超分辨率重建方法及系统
CN108898642A (zh) * 2018-06-01 2018-11-27 安徽工程大学 一种基于卷积神经网络的稀疏角度ct成像方法
CN108961237A (zh) * 2018-06-28 2018-12-07 安徽工程大学 一种基于卷积神经网络的低剂量ct图像分解方法

Also Published As

Publication number Publication date
CN111951346A (zh) 2020-11-17

Similar Documents

Publication Publication Date Title
Li et al. Motion correction for improved target localization with on-board cone-beam computed tomography
CN107072595B (zh) 基于多模态成像的自适应重计划
Schreibmann et al. Image interpolation in 4D CT using a BSpline deformable registration model
Lu et al. Tomographic motion detection and correction directly in sinogram space
CN110390361B (zh) 一种基于运动补偿学习的4d-cbct成像方法
Brehm et al. Self‐adapting cyclic registration for motion‐compensated cone‐beam CT in image‐guided radiation therapy
Sarrut et al. Nonrigid registration method to assess reproducibility of breath-holding with ABC in lung cancer
JP4356831B2 (ja) 複数画像間の非剛体レジストレーション方法
CN109961419B (zh) 对pet活度分布图像进行衰减校正的校正信息获取方法
CN110570489B (zh) 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法
US20190046145A1 (en) Methods for performing digital subtraction angiography, hybrid imaging devices, computer programs, and electronically readable storage media
CN111951346B (zh) 一种联合运动估计与时空张量增强表示的4d-cbct重建方法
CN114387364A (zh) 用于pet图像重建的线性衰减系数获取方法及重建方法
KR101460908B1 (ko) 4차원 컴퓨터 단층촬영 영상의 폐종양 위치 추적 시스템 및 그 방법
CN116993848B (zh) Cbct图像重建方法、系统、计算机装置及存储介质
Chen et al. Motion-compensated mega-voltage cone beam CT using the deformation derived directly from 2D projection images
CN115439572A (zh) 一种衰减校正系数图像的获取方法、pet图像重建方法
CN115049752A (zh) 一种基于三维卷积神经网络的pet呼吸运动图像伪影配准校正方法
CN116322902A (zh) 图像配准系统和方法
He et al. Helical mode lung 4D-CT reconstruction using Bayesian model
Li et al. High-Quality 4D-CBCT Imaging from Single Routine Scan
CN118014909A (zh) 一种基于自监督先验图像引导的4d cbct成像方法
Zhang et al. A review on 4D cone‐beam CT (4D‐CBCT) in radiation therapy: Technical advances and clinical applications
Qi et al. Reconstruction with angular compensation in respiratory-gated cardiac SPECT
Jin et al. Non-local means denoising of SG-KS-4D-MRI using block matching 3D: Implications for pancreatic tumor registration and segmentation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant