CN110570489B - 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法 - Google Patents

一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法 Download PDF

Info

Publication number
CN110570489B
CN110570489B CN201910837231.0A CN201910837231A CN110570489B CN 110570489 B CN110570489 B CN 110570489B CN 201910837231 A CN201910837231 A CN 201910837231A CN 110570489 B CN110570489 B CN 110570489B
Authority
CN
China
Prior art keywords
phase
image
dvf
cbct
quality
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
CN201910837231.0A
Other languages
English (en)
Other versions
CN110570489A (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.)
Affiliated Hospital of Jiangsu University
First Affiliated Hospital of Chongqing Medical University
Original Assignee
Affiliated Hospital of Jiangsu University
First Affiliated Hospital of Chongqing Medical 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 Affiliated Hospital of Jiangsu University, First Affiliated Hospital of Chongqing Medical University filed Critical Affiliated Hospital of Jiangsu University
Priority to CN201910837231.0A priority Critical patent/CN110570489B/zh
Publication of CN110570489A publication Critical patent/CN110570489A/zh
Priority to PCT/CN2020/096443 priority patent/WO2021042807A1/zh
Priority to US17/276,938 priority patent/US11386591B2/en
Application granted granted Critical
Publication of CN110570489B publication Critical patent/CN110570489B/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/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
    • 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/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
    • A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
    • A61B6/5264—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise due to motion
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00—Radiation therapy
    • A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A61N5/1048—Monitoring, verifying, controlling systems and methods
    • A61N5/1049—Monitoring, verifying, controlling systems and methods for verifying the position of the patient with respect to the radiation beam
    • 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
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00—Radiation therapy
    • A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A61N5/1048—Monitoring, verifying, controlling systems and methods
    • A61N5/1049—Monitoring, verifying, controlling systems and methods for verifying the position of the patient with respect to the radiation beam
    • A61N2005/1052—Monitoring, verifying, controlling systems and methods for verifying the position of the patient with respect to the radiation beam using positron emission tomography [PET] single photon emission computer tomography [SPECT] imaging
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2210/00—Indexing scheme for image generation or computer graphics
    • G06T2210/41—Medical
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Biomedical Technology (AREA)
  • Medical Informatics (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Radiology & Medical Imaging (AREA)
  • Pathology (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Biophysics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Optics & Photonics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Algebra (AREA)
  • Image Processing (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明公开了一种设计基于双边滤波的运动补偿高质量4D‑CBCT图像重建方法，包括如下步骤：1)建立基于瞬时代数重建技术(SART)的3D‑CBCT运动补偿重建模型，重建出参考相位(即0％相位)的高质量CBCT图像；2)构建迭代优化过程，估计0％相位与其他4D相位图像之间的4D图像变形场(Deformable Vector Field,DVF)模型，进而求解出包含运动器官表面之间反向滑动运动的精确4D‑DVF解；3)将高质量0％相位图像按照4D‑DVF最优解依次变形，获取最终高质量4D‑CBCT图像序列。该方法能在不改变目前常规放疗线性加速器硬件结构的基础上，实现高质量4D‑CBCT图像的精确重建；从而进一步提高肺癌立体定向放射治疗(Stereotactic Body Radiation Therapy,SBRT)前图像引导定位的精确性。

Description

一种设计基于双边滤波的运动补偿高质量4D-CBCT图像重建 方法
技术领域
本发明应用于针对癌变运动器官(如肺癌)在进行立体定向放射治疗(SBRT)过程中的精确图像引导放射治疗。
背景技术
针对肺癌的SBRT治疗,即体部立体定向放射治疗，目前越来越广泛的被应用在图像引导放射治疗过程中。其治疗效果已经被证实优于传统的调强放射治疗(IntensityModulated Radiation Therapy,IMRT)。但是SBRT治疗过程中采用了基于小射野(射野尺寸小于5x5厘米)的非均整滤波器射束模式(Flattening Filter Free,FFF)进行出束。其射束中心剂量率和强度相比于加了均整滤波器的调强射束来说，均提高许多，如图1所示。因此对SBRT治疗中的FFF射束进行精确定位，保证射束在精确打击癌变病灶的同时最大限度保护周围正常组织，进而实现SBRT治疗的安全性来说，意义重大。
这就要求医生在对病人进行SBRT治疗前，再次通过图像引导手段，精确确认治疗时刻病人肺癌区域的呼吸运动情况是否和SBRT治疗计划中的放疗靶区运动范围相吻合后再出束治疗。目前集成在线性加速器上的锥束CT(Cone Beam CT,CBCT)只能进行3D成像，其缺乏对运动器官的4D动态成像功能。
目前集成在线性加速器机架上的CBCT系统，只能进行传统3D-CBCT成像功能。因此CBCT目前只能用于放疗前摆位验证。尽管目前国际主流放疗加速器厂商已经提出基于解析法的4D-CBCT成像功能模块，但其成像效果欠佳，图象质量有待进一步改进。虽然目前最近推出的加速器CBCT已经提供了4D-CBCT成像功能，但其采用的成像方法伪影较多，图像质量较差。对于一些需要精确定位引导SBRT治疗的肺癌场景，不能起到良好的监控引导作用。
发明内容
针对现有技术中存在的上述不足之处，本发明提出了一种设计基于双边滤波的运动补偿高质量4D-CBCT图象重建方法，双边滤波操作实施在4D-DVF迭代优化过程中。
为了解决上述技术问题，本发明采用了如下技术方案：
一种设计基于双边滤波的运动补偿高质量4D-CBCT(4D-Cone Beam ComputedTomography)图像重建方法，该方法能够实现运动器官表面反向滑动运动的精确建模重建；该方法包括如下步骤：
1)建立基于瞬时代数重建技术的3D-CBCT运动补偿重建模型，重建出参考相位的高质量CBCT图像，参考相位为0％相位；
2)构建迭代优化过程，估计0％相位与其他4D相位图像之间的4D图像变形场模型，进而求解出包含运动器官表面之间反向滑动运动的精确4D-DVF(4D-Deformable VectorField)解；
3)将高质量0％相位图像按照4D-DVF最优解依次变形，获取最终高质量4D-CBCT图像序列。
作为本发明的一种优选方案，在步骤1)的建模中，将0％相位与其他4D-CBCT相位的4D-DVF优化解带入SART(simultaneous algebraic reconstruction technique)重建过程，进而得到改进后的运动补偿SART算法，可以获得基于整个4D投影数据的高质量0％相位初始图像；运动补偿SART的数学模型描述如下：
令
Figure BDA0002192566430000021
表示对数压缩后相位t的4D-CBCT投影线积分,且
Figure BDA0002192566430000022
表示相位t图像的衰减系数；运动补偿SART重建的数学模型为:
Figure BDA0002192566430000023
k是迭代次数；j是相位0％的图像体素索引；n是相位t图像的体素索引；ain是投影DRR射线i与其穿过图像体素n的穿透长度；
Figure BDA0002192566430000031
表示将相位t图像变形到相位0的反向DVF；DVF的引入即体现了运动补偿的思想；％相位初始图像
Figure BDA0002192566430000032
通过全变差最小化重建获得。
作为本发明的另一种优选方案，对每一个4D相位的DVF求解，该算法通过配准1)：该相位的测量投影数据；和2)：变形到该相位的0％参考相位图像的模拟正向投影数据；来产生该求解相位的DVF；在投影配准优化DVF的过程中采用的能量函数如下：
Figure BDA0002192566430000033
Figure BDA0002192566430000034
s.t.v0→t(x+vt→0)+vt→0＝vt→0(x+v0→t)+v0→t＝0 (2)
其中，f1和f2表示对称的反向连续能量函数；其上标0表示0％相位，t表示其他相位t；A是对变形后的0％相位图像求正向投影的矩阵；v0→t表示正向变形场DVF；vt→0表示反向变形场DVF；
Figure BDA0002192566430000035
和
Figure BDA0002192566430000036
为数据保真项；
Figure BDA0002192566430000037
和
Figure BDA0002192566430000038
是相对应的归一化项；DVF的反向连续性由上面方程中第三行保证；β控制着平滑项的平滑程度；为了考虑到运动器官表面的反向滑动运动建模，
Figure BDA0002192566430000039
被定义为：
Figure BDA00021925664300000310
其中：
Figure BDA00021925664300000311
表示DVF中每一个体素周围领域的差分；
Figure BDA00021925664300000312
Figure BDA00021925664300000313
Figure BDA0002192566430000041
Gx是空间域的高斯子滤波核，其方差是
Figure BDA0002192566430000042
Gμ是像素强度域的高斯子核，方差是
Figure BDA0002192566430000043
Gv是DVF域方差为
Figure BDA0002192566430000044
的子核；'xj'是子核中心的体素，'yj'是N(xj)范围内的相邻体素；
Figure BDA0002192566430000045
是在3x3x3的领域范围内计算；经过推导，得到
Figure BDA0002192566430000046
的梯度如下：
Figure BDA0002192566430000047
作为本发明的一种改进方案，在得到高质量0％相位CBCT图像和最优4D-DVF解基础上，可以通过变形该0％相位图像获得最终高质量4D-CBCT图像序列；该过程的数学描述如下：
Figure BDA0002192566430000048
Figure BDA0002192566430000049
表示相位t的CBCT图像可以通过用4D-DVFs，即
Figure BDA00021925664300000410
来变形
Figure BDA00021925664300000411
而获得。
本发明提供了一套能够提供高质量图像的运动补偿4D-CBCT图像重建方法。该方法利用所有4D相位投影，通过精确运动补偿重建构建一个初始参考相位的CBCT高质量图像。之后，对该参考相位和其他各个4D相位之间的DVF进行优化求解，得到最优4D-DVF。最终，通过用4D-DVF变形高质量初始参考相位图象，得到高质量4D-CBCT图像序列。在对4D-DVF进行估计的过程中，采用对DVF进行双边滤波，不仅可以获得整个胸腔内的整体均匀运动信息，还可以获得运动器官表面发生的局部非均匀反向滑动运动信息，这进一步提高了4D-DVF估计的精确性，从而获得精确的4D-CBCT重建图像序列。
与现有技术相比，本发明具有如下技术效果：
1、本发明的方法不仅可以对运动器官的整体均匀运动进行估计，还能对运动器官表面局部非均匀运动进行精确估计。
2、无需对加速器CBCT成像硬件系统和临床扫描协议进行任何改动，完全通过本发明的方法即能获得高质量4D-CBCT图象序列。通过本发明的方法获得的高质量4D-CBCT图象可以为肺癌SBRT放疗提供亚毫米级误差的肿瘤运动4D轨迹。大大加强SBRT肺癌治疗的安全性和肿瘤追踪的精确性。该方法弥补了目前放疗加速器CBCT无法进行高质量4D-CBCT成像的短板，可以作为肺癌SBRT治疗中一个有效的图象引导工具。
附图说明
图1为现有技术中的FFF射束和均整射束切面曲线剂量率对比图；
图2为基于双边滤波的高质量4D-CBCT运动补偿技术路线框图；
图3为4D NCAT数字模体重建结果对比图；
图4为NCAT体模ROI区域4D运动轨迹(z轴方向)重建结果对比图；
图5为病人临床数据算法验证结果图；
图6为第二例临床数据重建结果对比图。
具体实施方式
一种设计基于双边滤波的运动补偿高质量4D-CBCT(4D-Cone Beam ComputedTomography)图像重建方法，该方法包括如下步骤：
1、4D-CBCT投影数据采集
首先进行4D-CBCT投影数据采集。采用目前成熟的基于照相机的呼吸运动检测系统，通过即可获取病人治疗时随时间的的呼吸运动幅度和相位波形。随着CBCT投影数据的采集，每一帧采集到的投影数据都会被监测到的呼吸信号打上相应的时间标签。在完成整个投影数据采集后，带有呼吸信号时间和幅度标签的CBCT投影数据就可以按照呼吸幅度或者相位信号进行4D相位划分，从而获得4D-CBCT投影数据。
该步骤属于数据预采集过程。具体实施中，采用瓦里安(Varian)的RespirationPositioning Management(RPM)系统来获取4D-CBCT投影数据。该系统通过在病人腹部放置一个红外光线反光球底座，同时在治疗室屋顶安装红外照相机用于监控反光球运动幅度。病人随时间的呼吸幅度和相位信号就可被红外相机记录到。记录到的呼吸信号分配给CBCT系统采集到的含病人呼吸运动的投影数据，相当于给每一个采集到的投影数据打时间标签。之后对这些投影数据按照呼吸运动曲线的时间标签进行4D分组即可得到4D-CBCT投影数据。
2、4D-CBCT图像初始重建和4D-DVF初始产生
利用全变差最小化(Total Variation minimization,TV)和Demons配准算法分别进行4D-CBCT初始重建和4D-DVF(4D-Deformable Vector Field)的初始产生。
采集到的4D-CBCT投影首先通过全变差最小化(Total Variation minimization,TV)重建，得到4D-CBCT初始图像序列。采用TV重建方法的原因是该方法能有效去除重建图像中的噪声，同时对图像边界有良好的保留功能。在投影经过4D相位划分后，每个相位的投影数目约在50到60个之间。利用如此稀疏的投影进行初始4D-CBCT重建，传统解析法会产生强烈的条纹伪影，严重恶化图像质量。而TV方法能避免条纹伪影对图像的污染，从而保证后续4D-CBCT图像质量的进一步提升。
得到初始4D-CBCT图像后，采用Demons配准算法来产生0％参考相位和其他各个4D相位之间的4D-DVF初始值。Demons配准算法是一种经典的基于图像像素值相似度的可变形图像配准算法，其能产生两幅配准图像之间的三维图像变形场DVF。利用Demons配准算法，得到4D-DVF初始值，用该初始值触发4D-DVF的进一步优化求解，有利于快速求得最优解。
3、4D-DVF迭代优化求解
0％参考相位图像与某一个其他4D相位之间的DVF，通过配准该目标相位的测量投影数据和变形后的0％相位之正向投影数据而产生。其过程如公式(2)所示：pt表示目标求解相位下测量到的投影数据。μ0(x+v0→t)表示DVF像素v0→t，将0％相位图像μ0从参考相位0％变形到目标相位t。而变形后的0％相位图像前面加上正投矩阵A后就得到该变形后图像的正投数据。
Figure BDA0002192566430000061
表示目标相位t的测量数据和变形到目标相位t后的0％相位图像之正投数据之间差值的l2范数，即差平方。这构成了DVF优化目标函数的保真项。对DVF优化目标函数的惩罚项
Figure BDA0002192566430000071
设计，兼顾了让DVF自然平滑和保持运动器官局部非均匀滑动运动的提取两方面。具体的平滑项设计如方程(3)所示：该平滑项主要包含基于双边滤波的三个高斯滤波核。分别是基于空间域的高斯核Gx，基于图像像素值域的高斯核Gμ，和基于图像变形场域的高斯核
Figure BDA0002192566430000072
在优化过程中，通过设置不同的高斯核相对应的方差值，如Gx对应的方差
Figure BDA0002192566430000073
Gμ对应的方差
Figure BDA0002192566430000074
和
Figure BDA0002192566430000075
对应的方差
Figure BDA0002192566430000076
就能调节整个双边滤波的效果。最终让DVF既能估计出运动器官的整体均匀运动矢量场，又能抓取到运动器官表面发生的局部非均匀反向滑动运动。
如前述所述，通过配准目标相位投影数据和变形后的参考相位正投数据来调整两相位之间的DVF。在此过程中，需要针对双边滤波核中的不同方差赋值。具体实施中，首先通过对理想情况下的4D人体呼吸和心血管运动数字模体,简称NCAT模体(4D NURBS-basedCardiac-Torso(NCAT)Phantom)进行算法验证。发现当σx＝3毫米时，能得到最理想的结果。过大或者过小的σx要么会过平滑DVF，要么不足够抓取局部空间特征。对于图像像素值域方差σμ，为了让图像像素从胸壁处自然过渡到肺部组织，令σμ等于肺部组织处像素平均值和邻近胸壁处像素平均值之差。针对NCAT模体，σμ＝0.03mm-1满足此条件并且重建结果理想。最后也是最关键的，就是如何确定合适的DVF域方差σv。理论上，σv应该比在σx空间范围内的任何两点之间的DVF差值都大，除非当σx包含胸壁像素值突变区域。在胸壁像素突变区域，σv应该比此区域的DVF像素差值小。这样才能滤除较小的像素运动量同时保留胸壁处较大的像素运动。为了避免过度分割，设定只有较大的锐利运动不连续性运动可以被保留(如运动器官表面的反向滑动运动)。通过对临床病人数据的观察，认为10毫米的移动量可以相对认为是比较大的滑动运动阈值。在此假设下，发现当σv＝2.5mm时，能获得基于NCAT模体数据较好的重建结果。
在基于NCAT模体验证算法性能基础上，进一步使用病人临床数据对算法进行了验证。由于临床数据是在真实环境下采集的，其受到一系列因素的干扰包括噪声、散射、射束硬化等。适合初步临床数据特点的相应滤波核方差值与前述适合NCAT数字模体的房差之略有区别。发现σx＝3mm，σμ＝0.02mm-1和σv＝2mm比较适合临床数据，并且能取得较好的临床重建结果。
由上述方差值进行DVF的迭代优化估计，对目标函数的优化中，采用非线性共轭梯度算子，因此需要计算目标函数的梯度。目标函数中惩罚项的梯度推导如方程(4)所示。
4、0％参考相位的不断迭代更新重建
在每一次4D-DVF迭代更新后，都要对0％相位参考图像进行更新重建，从而不断提高参考相位图像的图像质量。重建过程如方程(1)所描述。其中的DVF像素
Figure BDA0002192566430000081
随着每次4D-DVF的更新而更新。
具体实施中，4D-DVF的优化估计和0％参考相位图像更新是交互式进行的。即对4D-DVF先进行n次的迭代优化(n可以设置为10～20次)，并保存这时刻的4D-DVF中间值。然后用此时的4D-DVF中间值带入运动补偿SART算法，重建出该n次4D-DVG迭代优化后的0％相位图像。然后再进行下一个n次的4D-DVF迭代优化后，实施第二次运动补偿SART(simultaneous algebraic reconstruction technique)重建获取第二次更新后的0％相位图像。每一次4D-DVF迭代优化后，其各个相位的目标函数值均被记录下来。在第m个运动补偿SART重建迭代后，如果各个相位的4D-DVF优化目标函数走势均达到收敛状态，则停止整个4D-DVF优化和运动补偿SART重建。并保存此时已经达到最优状态的4D-DVF和经历了m次运动补偿SART重建的0％相位高质量CBCT图像。
5、最终高质量4D-CBCT图像序列生成
随着4D-DVF优化过程中目标函数的不断减小并最终趋向稳定，其得到的最优解即被用来变形经过不断迭代更新重建后的0％参考相位图像。最终得到高质量4D-CBCT图像序列。整个过程如方程(5)所示。
上述步骤的技术路线框图如附图2所示。
实施例1
首先在4D NCAT数字模体上验证该算法的有效性。该4D模体呼吸周期为4s。横膈膜沿上内侧(SI)的最大运动为20mm，胸部前后(AP)的最大运动为12mm。模体有10个相位，每个相位模拟出20个投影，用于DVF估计和4D-CBCT重建。模体图像尺寸为256x256x150，体素大小为1x1x1 mm3。投影尺寸为300x240x20，投影体素尺寸为1x1x1 mm3。
图3所示为在有和没有滑动运动约束的情况下获得的40％相位重建图像。3(a)显示了在不进行双边滤波的情况下获得的重建40％相位的矢状面视图；图3(b)显示了用双边滤波重建的相同矢状面切片；3(c)显示了模体金标准参考图片；3(d)、3(e)和3(f)是感兴趣区域(Region Of Interest,ROI)，其中滑动运动发生在心脏边缘和静脉部位。ROI显示心脏边缘边界(由(d)、(e)和(f)中的箭头指示)通过双边过滤得到了良好的获取。此外，静脉(用黄点线标记)通过双边过滤更精确地被重建(例如静脉长度更好地被重建)。图3(g)、(h)和(i)显示了带和不带双边滤波的肋骨位置重建差异。在图3(h)中，肋骨顶部边缘1和2与图3(i)中的金标准肋骨参考位置相匹配，而图3(g)中没有进行双边过滤。
为了检测运动器官边界处由于滑动运动影像而产生的运动误差，在上图中心脏边缘处(参见图3(f)中标记位置的点线)，提取该处沿z轴方向的4D运动轨迹。每一幅4D图像中，相应检测到的点线位置用于绘制运动轨迹。图4表明，采用基于双边滤波的滑动考虑提取的轨迹与4D轨迹金标准的匹配程度要比不考虑此因素的轨迹好得多。取轨道的每一个运动量进行均方根误差(Root Mean Square Error,RMSE)和最大运动误差MaxE计算。RMSE的定义如下：
Figure BDA0002192566430000091
Figure BDA0002192566430000101
表示重建图像中第phth个相位图像中的特征结构位置；
Figure BDA0002192566430000102
表示金标准图像中相应相位中同一个特征结构的位置；ph表示相位索引；此处4D相位划分为10个。MaxE表示在10个相位中最大的那个误差值。
结果发现基于双边滤波的结果，4D轨迹的RMSE/MaxE分别为0.796mm和1.02mm。与不带双边滤波的原始重建结果相比，该轨迹的RMSE/MaxE分别为2.704mm和4.08mm。
在NCAT模体数据基础上，进行了初步病人临床数据实验。采集了两个肺癌患者的CBCT投影数据，对基于双侧滤波的4D-CBCT滑动补偿重建算法进行了初步试验。为了找到临床数据的金标准，进行量化比较，提交了IRB临床伦理批件。每个病人数据实在全扇扫模式(Full Fan Mode)下进行的投影数据采集。扫描时间为4到6分钟。获得大约2000个投影。然后利用采集过程中相应的RPM信号，将获得的投影分为10个相位。以这种方式，每个阶段的平均投影数约为200个，并应用TV最小化重建技术重建高质量4D-CBCT，可作为临床结果量化的参考金标准图像。这种情况在临床常规扫描协议下是不允许的。这样做的目的只是为了采集到足够多的投影数据，保证即使在4D相位划分后，每个相位仍然有足够多投影数目重建出高质量4D-CBCT临床数据图像，作为参考金标准进行算法临床效果评价验证。
考虑到CBCT在临床实际1分钟扫描协议下，几下旋转360度最多只能采集到约670个投影数据。经过4D相位划分后，则每个4D相位实际上只能获得约60个投影。为了模拟真实的1分钟CBCT扫描场景，对采集到的4D全投影进行了欠采样操作。即把每个相位的平均投影欠采样至约为40。然后采用几种不同的重建方法得到一系列平行对照组4D-CBCT重加结果和提出的基于双边滤波的运动补偿重建结果进行对比。所采用的重建方法包括：传统FDK重建、TV最小化重建、以及在算法中不考虑滑动运动建模的重建。
相应的结果如图5所示。图5(a)、(b)显示了无/有基于双侧滤波的滑动运动约束的第一名患者重建结果的矢状面视图。图5(c)显示了由电视从全采样投影重建的患者参考图像。箭头所示为一个紧贴胸壁的肿瘤，肿瘤和胸壁之间存在一个小洞。与5(c)中的参考值相比，5(a)中的小空腔明显模糊。与对照组比较，5(b)组肿瘤边界较5(a)组明显。这表明基于双边滤波的重建比不采用双边滤波的重建获得了更多的参考可比结果。为了进行综合比较，图5(d)和图5(e)分别列出了通过TV和FDK的相应重建切片。这两个结果表明，TV重建图像中的结构已经过平滑；同时，噪声污染严重恶化了FDK图像中的图像内容。
另一方面，图5(f)、5(g)和5(h)显示了该患者的肺肋笼边界，因此可以比较其肋骨运动。点线表示肋骨上边缘位置。箭头表明，在不进行双边滤波(如图5(f))的重建中，与双边滤波重建结果(图5(g))和病人参考图像图5(h)中的参考肋骨位置相比，图5(f)中肋骨顶部边缘略微向上偏移约2mm。因此，双边滤波还是较好的矫正了肋骨在重建图像中的真实位置。
同时还用另一个病人数据验证了该算法。图6显示了第二例病人临床数据重建的冠状面视图结果。图6(a)是未经双边滤波的原始运动补偿重建图像。图6(b)显示了双边滤波重建结果。图6(c)是通过TV最小化从全投影数据重建的病人临床参考图像。图6(d)、6(e)和6(f)分别是6(a)、6(b)和6(c)中相同位置的缩放ROI(见6(a)中的ROI框)。它清楚地表明，与图6(f)相比，6(e)中的肺心边界(用箭头标记)比6(d)中的更清晰。这进一步证明了双边滤波的优点，与原始运动补偿结果相比，可以获得更精确的重建结果。
最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (3)

1.一种设计基于双边滤波的运动补偿高质量4D-CBCT图像重建方法，其特征在于，该方法能够实现运动器官表面反向滑动运动的精确建模重建；该方法包括如下步骤：
1)建立基于瞬时代数重建技术SART的3D-CBCT运动补偿重建模型，重建出参考相位的高质量CBCT图像，参考相位为0％相位；
2)构建迭代优化过程，估计0％相位与其他4D相位图像之间的4D图像变形场模型，进而求解出包含运动器官表面之间反向滑动运动的精确4D-DVF(4D-Deformable VectorField)解；
3)将高质量0％相位图像按照4D-DVF最优解依次变形，获取最终高质量4D-CBCT图像序列；
对每一个4D相位的DVF求解，求解算法通过配准：1)该相位的测量投影数据；和2)变形到该相位的0％参考相位图像的模拟正向投影数据；来产生求解相位的DVF；在投影配准优化DVF的过程中采用的能量函数如下：
Figure FDA0003107411850000011
Figure FDA0003107411850000012
s.t.v0→t(x+vt→0)+vt→0＝vt→0(x+v0→t)+v0→t＝0 (2)
其中，f1和f2表示对称的反向连续能量函数；其上标0表示0％相位,t表示其他相位t；A是对变形后的0％相位图像求正向投影的矩阵；v0→t表示正向变形场DVF；vt→0表示反向变形场DVF；
Figure FDA0003107411850000013
和
Figure FDA0003107411850000014
为数据保真项，其中x代表图像的每一个体素，p0代表参考相位即0％相位的测量投影数据，pt代表相位t的测量投影数据，μ0(x+v0→t)表示经过正向变形后的参考相位图像 ，μt(x+vt→0)表示经过反向变形的相位图像；其前面加上正向投影变换矩阵A后则表示将变形后的图像 转换为该相位各个投影角度下的正投数据；
Figure FDA0003107411850000015
和
Figure FDA0003107411850000016
是相对应的归一化项；β控制着平滑项的平滑程度；DVF的反向连续性由上面方程中第三行保证；为了考虑到运动器官表面的反向滑动运动建模，
Figure FDA0003107411850000021
被定义为:
Figure FDA0003107411850000022
其中：
Figure FDA0003107411850000023
表示DVF中每一个体素周围领域的差分，j表示图像 中体素的纬度索引，i表示DVF沿x,y,z三个方向的索引；三个子滤波核的具体形式如下：
Figure FDA0003107411850000024
Figure FDA0003107411850000025
Figure FDA0003107411850000026
其中Gx是空间域的高斯子滤波核，其方差是
Figure FDA0003107411850000027
Gμ是像素强度域的高斯子核，方差是
Figure FDA0003107411850000028
Gv是DVF域方差为
Figure FDA0003107411850000029
的子核；'xj'是子核中心的体素，'yj'是N(xj)范围内的相邻体素；
Figure FDA00031074118500000210
表示
Figure FDA00031074118500000211
针对v的偏导数，其是在3x3x3的领域范围内计算；经过推导，得到：
Figure FDA00031074118500000212
2.根据权利要求1所述的一种设计基于双边滤波的运动补偿高质量4D-CBCT图像重建方法，其特征在于，在步骤1)的建模中，将0％相位与其他4D-CBCT相位的4D-DVF优化解带入SART重建过程，进而得到改进后的运动补偿SART算法，可以获得基于整个4D投影数据的高质量0％相位初始图像；运动补偿SART的数学模型描述如下：
令
Figure FDA00031074118500000213
表示对数压缩后相位t的4D-CBCT投影线积分,且
Figure FDA0003107411850000031
表示相位t图像的衰减系数；运动补偿SART重建的数学模型为:
Figure FDA0003107411850000032
k是迭代次数；j是相位0％的图像体素索引；n是相位t图像的体素索引；ain是正向投影模拟射线i与其穿过图像体素n的穿透长度；
Figure FDA0003107411850000033
表示将相位t图像变形到相位0的反向DVF；DVF的引入即体现了运动补偿；0％相位初始图像
Figure FDA0003107411850000034
通过全变差最小化重建获得。
3.根据权利要求1所述一种设计基于双边滤波的运动补偿高质量4D-CBCT图像重建方法，其特征在于，在得到高质量0％相位CBCT图像和最优4D-DVF解基础上，可以通过变形该0％相位图像获得最终高质量4D-CBCT图像序列；该过程的数学描述如下：
Figure FDA0003107411850000035
Figure FDA0003107411850000036
表示相位t的CBCT图像可以通过用4D-DVF，即
Figure FDA0003107411850000037
来变形
Figure FDA0003107411850000038
而获得。
CN201910837231.0A 2019-09-05 2019-09-05 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法 Active CN110570489B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN201910837231.0A CN110570489B (zh) 2019-09-05 2019-09-05 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法
PCT/CN2020/096443 WO2021042807A1 (zh) 2019-09-05 2020-06-17 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法
US17/276,938 US11386591B2 (en) 2019-09-05 2020-06-17 Reconstruction method for motion compensated high quality 4D-CBCT image based on bilateral filtering

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910837231.0A CN110570489B (zh) 2019-09-05 2019-09-05 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法

Publications (2)

Publication Number Publication Date
CN110570489A CN110570489A (zh) 2019-12-13
CN110570489B true CN110570489B (zh) 2021-07-30

Family

ID=68778025

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910837231.0A Active CN110570489B (zh) 2019-09-05 2019-09-05 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法

Country Status (3)

Country Link
US (1) US11386591B2 (zh)
CN (1) CN110570489B (zh)
WO (1) WO2021042807A1 (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110570489B (zh) 2019-09-05 2021-07-30 重庆医科大学附属第一医院 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法
CN111583296B (zh) * 2020-04-28 2022-09-09 山东大学 基于卡尔曼滤波和4d ct图像配准的肺部呼吸运动估计方法
CN111932475B (zh) * 2020-07-31 2023-07-04 东软医疗系统股份有限公司 滤波方法、装置、ct设备及ct系统
CN111951346B (zh) * 2020-08-18 2023-06-23 安徽工程大学 一种联合运动估计与时空张量增强表示的4d-cbct重建方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101743568A (zh) * 2007-05-17 2010-06-16 伊利克塔股份有限公司 在线锥形束ct重建
CN104361568A (zh) * 2014-09-18 2015-02-18 南方医科大学 基于配准的肺4d-ct图像呼气过程中间相位图像的重建方法
CN107206251A (zh) * 2014-12-10 2017-09-26 医科达有限公司 用于构建四维图像信息的磁共振投影

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101428005B1 (ko) * 2012-10-29 2014-08-07 한국과학기술원 소수의 저선량 ｃｔ 영상을 이용하여 ｐｅｔ 영상을 움직임 보상 및 감쇠 보정하는 방법
CN107886478B (zh) 2017-09-22 2020-10-30 深圳先进技术研究院 一种ct图像重建方法及系统、终端及可读存储介质
CN110570489B (zh) 2019-09-05 2021-07-30 重庆医科大学附属第一医院 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101743568A (zh) * 2007-05-17 2010-06-16 伊利克塔股份有限公司 在线锥形束ct重建
CN104361568A (zh) * 2014-09-18 2015-02-18 南方医科大学 基于配准的肺4d-ct图像呼气过程中间相位图像的重建方法
CN107206251A (zh) * 2014-12-10 2017-09-26 医科达有限公司 用于构建四维图像信息的磁共振投影

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
DeepOrganNet: On-the-Fly Reconstruction and Visualization of 3D/ 4D Lung Models from Single-View Projections by Deep Deformation Network;Yifan Wang,等;《https://arxiv.org/search/?query=DeepOrganNet%3A+On-the-Fly+Reconstruction+and+Visualization+of+3D+%2F+4D+Lung+Models+from+Single-View+Projections+by+Deep+Deformation+Network&searchtype=all&source=header》;20190722;全文 *
Simultaneous 4D-CBCT reconstruction with sliding motion constraint;Dang Jun,等;《MEDICAL PHYSICS》;20161031;第5453-5458页 *
基于运动补偿的压缩感知4D-CBCT优质重建;杨轩，等;《南方医科大学学报》;20160629;全文 *
运动配准先验图像的4D-CBCT优质重建;陈美玲，等;《南方医科大学学报》;20190314;全文 *

Also Published As

Publication number Publication date
CN110570489A (zh) 2019-12-13
US20210248790A1 (en) 2021-08-12
US11386591B2 (en) 2022-07-12
WO2021042807A1 (zh) 2021-03-11

Similar Documents

Publication Publication Date Title
CN110570489B (zh) 一种设计基于双边滤波的运动补偿高质量4d-cbct图像重建方法
Xie et al. Scatter artifacts removal using learning-based method for CBCT in IGRT system
US8923590B2 (en) Method and system for 3D cardiac motion estimation from single scan of C-arm angiography
Chou et al. 2D/3D image registration using regression learning
US10524756B2 (en) Methods and systems for image artifacts reduction
CN112435307B (zh) 一种深度神经网络辅助的四维锥束ct影像重建方法
CA2583831A1 (en) Method and apparatus for metal artifact reduction in computed tomography
Zhang et al. A biomechanical modeling guided CBCT estimation technique
CN110390361B (zh) 一种基于运动补偿学习的4d-cbct成像方法
Qi et al. 4-D reconstruction with respiratory correction for gated myocardial perfusion SPECT
CN111275669B (zh) 一种先验信息引导的四维锥束ct图像重建算法
Manjeshwar et al. Motion compensated image reconstruction of respiratory gated PET/CT
Dang et al. A pilot evaluation of a 4-dimensional cone-beam computed tomographic scheme based on simultaneous motion estimation and image reconstruction
Lewis et al. Mitigation of motion artifacts in CBCT of lung tumors based on tracked tumor motion during CBCT acquisition
Qi et al. A quantitative study of motion estimation methods on 4D cardiac gated SPECT reconstruction
Zhao et al. Local metric learning in 2D/3D deformable registration with application in the abdomen
CN105976412B (zh) 一种基于离线字典稀疏正则化的低管电流强度扫描的ct图像重建方法
Star‐Lack et al. A modified McKinnon‐Bates (MKB) algorithm for improved 4D cone‐beam computed tomography (CBCT) of the lung
McClelland Estimating internal respiratory motion from respiratory surrogate signals using correspondence models
CN111951346B (zh) 一种联合运动估计与时空张量增强表示的4d-cbct重建方法
Hirai et al. Regression model-based real-time markerless tumor tracking with fluoroscopic images for hepatocellular carcinoma
Dang et al. Fully automatic sliding motion compensated and simultaneous 4D-CBCT via bilateral filtering
Turco et al. Impact of CT-based attenuation correction on the registration between dual-gated cardiac PET and high-resolution CT
Dang et al. Deformation vector fields (DVF)-driven image reconstruction for 4D-CBCT
Bögel et al. Diaphragm tracking in cardiac C-Arm projection data

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