CN114387359A - 一种三维x射线低剂量成像方法及装置 - Google Patents
一种三维x射线低剂量成像方法及装置 Download PDFInfo
- Publication number
- CN114387359A CN114387359A CN202111453783.5A CN202111453783A CN114387359A CN 114387359 A CN114387359 A CN 114387359A CN 202111453783 A CN202111453783 A CN 202111453783A CN 114387359 A CN114387359 A CN 114387359A
- Authority
- CN
- China
- Prior art keywords
- projection
- dimensional
- image
- data
- integral image
- 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.)
- Pending
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 38
- 238000000034 method Methods 0.000 claims abstract description 38
- 238000011084 recovery Methods 0.000 claims abstract description 33
- 230000008707 rearrangement Effects 0.000 claims abstract description 29
- 238000007408 cone-beam computed tomography Methods 0.000 claims abstract description 20
- 238000004422 calculation algorithm Methods 0.000 claims description 25
- 230000006870 function Effects 0.000 claims description 25
- 230000014509 gene expression Effects 0.000 claims description 17
- 238000001914 filtration Methods 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 9
- 238000005457 optimization Methods 0.000 claims description 9
- 239000000126 substance Substances 0.000 claims description 9
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 8
- 238000009795 derivation Methods 0.000 claims description 4
- 230000009711 regulatory function Effects 0.000 claims description 4
- 230000001105 regulatory effect Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 2
- 238000002591 computed tomography Methods 0.000 abstract description 2
- 230000000694 effects Effects 0.000 description 7
- 230000008569 process Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000004590 computer program Methods 0.000 description 4
- 238000012805 post-processing Methods 0.000 description 4
- 210000001519 tissue Anatomy 0.000 description 4
- 230000005855 radiation Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 230000008901 benefit Effects 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 241000272443 Penelope Species 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 210000000481 breast Anatomy 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000002526 effect on cardiovascular system Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000013307 optical fiber Substances 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 210000004872 soft tissue Anatomy 0.000 description 1
- 238000013179 statistical model Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
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
-
- 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
- G06T15/00—3D [Three Dimensional] image rendering
- G06T15/005—General purpose rendering architectures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- 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/10004—Still image; Photographic image
- G06T2207/10012—Stereo images
-
- 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/10116—X-ray image
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Graphics (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
本发明公开了一种三维X射线低剂量成像方法及装置，方法包括采集空曝数据，估计三维X射线系统内参数；根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像；对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图；根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图；对所述理想投影弦图进行解重排，得到恢复的积分图像；对所述恢复的积分图像进行重建，得到目标CBCT断层图像。本发明噪声小且图像分辨率高，能够有效解决当前低剂量三维图像噪声大、图像分辨率差的技术痛点，可广泛应用于计算机断层成像技术领域。
Description
技术领域
本发明涉及计算机断层成像技术领域，尤其是一种三维X射线低剂量成像方法及装置。
背景技术
三维X射线成像，以锥形束计算机断层成像(Cone-beam Computed Tomography，CBCT)为代表，可极快地完成数据采集，且一次扫描即可三维成像，在牙科正畸、乳腺断层合成扫描、精准图像引导及4D-CBCT的心脑血管检查等临床领域有广泛的应用。为遵循低剂量原则(As Low As Reasonably Achievable,ALARA)，X射线低剂量成像方法是近年成像技术领域的研究热点。
低剂量成像主要通过降低X射线管电压或管电流、改进探测器对能量的响应、减少辐射时间及稀疏投影采样实现辐射剂量的降低。据相关研究，低剂量X射线对成像最关键的影响是图像的信噪比严重下降，部分或众多组织信息(特别是低对比度的软组织信息)被噪声淹没而无法满足诊断需求。为了解决此问题，许多专家学者已在CT领域开展了深入的研究。CT低剂量成像的解决策略主要分为对投影数据滤波、迭代重建及后处理去噪。其中投影滤波有非线性滤波器滤波，如常见的均值、中值滤波，以及基于统计模型的投影滤波，如泊松分布模型的投影去噪。但这类方法通常忽略射束能量多色性、探测器响应及邻域相关性等因素的影响，不能准确反应投影的噪声特性。迭代重建方法由于其天然的数学优化表达式而具有良好的抗噪声性能，当投影数据较少时仍旧可重建出较优的图像质量。但迭代重建的时间复杂度较高，目前大多数迭代重建算法还不能满足临床成像实时性的要求。后处理去噪方法也可显著降低重建图像的噪声，提高重建图像的信噪比。但大多数后处理方法都会不同程度地降低图像空间分辨率，造成组织结构边缘模糊，且后处理的滤波器设计会极大影响去噪的效果。
尽管对CT领域低剂量成像的研究众多，但对CBCT三维成像的低剂量研究却很少。一个重要原因CT探测器是线阵探测器，而CBCT的探测器是平板探测器，后者无法使用滤线栅对散射线进行滤除，使得本来就差的低剂量CBCT图像雪上加霜。三维低剂量成像的问题仍是亟待解决的技术痛点。
发明内容
有鉴于此，本发明实施例提供一种噪声小且图像分辨率高的，三维X射线低剂量成像方法及装置。
本发明的一方面提供了一种三维X射线低剂量成像方法，包括：
采集空曝数据，估计三维X射线系统内参数；
根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像；
对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图；
根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图；
对所述理想投影弦图进行解重排，得到恢复的积分图像；
对所述恢复的积分图像进行重建，得到目标CBCT断层图像。
可选地，所述采集空曝数据，估计三维X射线系统内参数，包括：
配置三维低剂量成像设备的几何参数，所述几何参数包括射线源到旋转中心距离、旋转中心到探测器的距离以及激光中心到成像中心的距离；
在预设角度下利用所述三位低剂量成像设备重复采集多张空曝数据；
根据所述空曝数据中投影序列，确定与投影相同大小的二维调节函数矩阵及尺度系数矩阵；
所述调节函数及尺度系数的表达式为：
其中，
为数组p(u,v,:)的均值；σ2(u,v)为数组p(u,v,:)的方差；g(u,v)为调节函数；an为尺度系数；N为尺度系数的阶数。
可选地，所述根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像，包括：
对所述原始投影数据进行余晖校准；
通过蒙特卡洛方法对所述余晖校准后的数据进行散射校准；
对所述散射校准后的数据进行硬化校准；
根据所述余晖校准、所述散射校准和所述硬化校准的结果，确定积分图像。
可选地，所述对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图这一步骤中，
所述数据重排的结果包括三维伪扇形束数据或三维伪平行束数据；
所述数据重排的表达式为：
其中，u表示积分图像中投影图像的行；v表示积分图像中投影图像的列；s表示积分图像中第s张投影图像；d为探测器像素大小；D为射线源到探测器的距离；p(s,v,u)代表积分图像数据；q(u′,v′,s′)代表重排后的三维伪扇束数据。
可选地，所述根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图，包括：
通过最大后验概率描述投影恢复模型；
对所述投影恢复模型进行求导，得到优化的迭代求解公式；
根据所述迭代求解公式，通过计算目标滤波窗与所述目标滤波窗下图像梯度的乘积，确定投影恢复模型；
根据滤波窗的加权均值和所述投影恢复模型，计算无噪声干扰的理想投影弦图；
其中，所述投影恢复模型的表达式为：
其中，
为待恢复的无损投影图像；qi为经过重排的投影弦图；
为噪声模型的一维表示；m为总像素数量；β为松弛因子；Ω为像素i的邻域；γij为纹理权重函数。
可选地，所述对所述理想投影弦图进行解重排，得到恢复的积分图像中，所述解重排的表达式为：
其中，p(u,v,s)代表解重排后恢复的积分图像；
代表无噪声干扰的理想投影弦图；u表示积分图像中投影图像的行；v表示积分图像中投影图像的列；s表示积分图像中第s张投影图像；d为探测器像素大小。
可选地，所述对所述恢复的积分图像进行重建，得到目标CBCT断层图像，包括：
采用经过锐化改造的迭代重建算法对所述恢复的积分图像进行重建；
根据所述重建图像确定目标函数；
其中，所述目标函数的表达式为：
其中，f表示迭代优化得到的重建图像；Pf表示对重建图像f进行正向投影；p为解重排的积分图像；μ为松弛调节因子；R(f)为正则化项；
所述重建图像采用GPU并行计算方式进行加速；所述正则化项通过梯度方差构造而来；所述迭代重建算法采用CGLS共轭梯度下降算法求解得到。
本发明实施例的另一方面还提供了一种三维X射线低剂量成像装置，包括：
第一模块，用于采集空曝数据，估计三维X射线系统内参数；
第二模块，用于根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像；
第三模块，用于对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图；
第四模块，用于根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图；
第五模块，用于对所述理想投影弦图进行解重排，得到恢复的积分图像；
第六模块，用于对所述恢复的积分图像进行重建，得到目标CBCT断层图像。
本发明实施例的另一方面还提供了一种电子设备，包括处理器以及存储器；
所述存储器用于存储程序；
所述处理器执行所述程序实现如前面所述的方法。
本发明实施例的另一方面还提供了一种计算机可读存储介质，所述存储介质存储有程序，所述程序被处理器执行实现如前面所述的方法。
本发明实施例还公开了一种计算机程序产品或计算机程序，该计算机程序产品或计算机程序包括计算机指令，该计算机指令存储在计算机可读存储介质中。计算机设备的处理器可以从计算机可读存储介质读取该计算机指令，处理器执行该计算机指令，使得该计算机设备执行前面的方法。
本发明的实施例首先采集空曝数据，估计三维X射线系统内参数；根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像；对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图；根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图；对所述理想投影弦图进行解重排，得到恢复的积分图像；对所述恢复的积分图像进行重建，得到目标CBCT断层图像。本发明噪声小且图像分辨率高，能够有效解决当前低剂量三维图像噪声大、图像分辨率差的技术痛点。
附图说明
为了更清楚地说明本申请实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的整体步骤流程图；
图2为本发明实施例提供的某一投影像素点序列的直方分布示意图；
图3为本发明实施例中投影重排前后的示意图；
图4为本发明实施例提供的投影恢复模型中像素均值的计算方法示意图；
图5为本发明实施例提供的低剂量重建方法处理前后图像效果对比图；
图6为本发明实施例提供的低剂量重建方法与其他方法重建图像对比图。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
针对现有技术存在的问题，本发明的一方面提供了一种三维X射线低剂量成像方法，包括：
采集空曝数据，估计三维X射线系统内参数；
根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像；
对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图；
根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图；
对所述理想投影弦图进行解重排，得到恢复的积分图像；
对所述恢复的积分图像进行重建，得到目标CBCT断层图像。
可选地，所述采集空曝数据，估计三维X射线系统内参数，包括：
配置三维低剂量成像设备的几何参数，所述几何参数包括射线源到旋转中心距离、旋转中心到探测器的距离以及激光中心到成像中心的距离；
在预设角度下利用所述三位低剂量成像设备重复采集多张空曝数据；
根据所述空曝数据中投影序列，确定与投影相同大小的二维调节函数矩阵及尺度系数矩阵；
所述调节函数及尺度系数的表达式为：
其中，
为数组p(u,v,:)的均值；σ2(u,v)为数组p(u,v,:)的方差；g(u,v)为调节函数；an为尺度系数；N为尺度系数的阶数。
可选地，所述根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像，包括：
对所述原始投影数据进行余晖校准；
通过蒙特卡洛方法对所述余晖校准后的数据进行散射校准；
对所述散射校准后的数据进行硬化校准；
根据所述余晖校准、所述散射校准和所述硬化校准的结果，确定积分图像。
可选地，所述对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图这一步骤中，
所述数据重排的结果包括三维伪扇形束数据或三维伪平行束数据；
所述数据重排的表达式为：
其中，u表示积分图像中投影图像的行；v表示积分图像中投影图像的列；s表示积分图像中第s张投影图像；d为探测器像素大小；D为射线源到探测器的距离；p(s,v,u)代表积分图像数据；q(u′,v′,s′)代表重排后的三维伪扇束数据。
可选地，所述根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图，包括：
通过最大后验概率描述投影恢复模型；
对所述投影恢复模型进行求导，得到优化的迭代求解公式；
根据所述迭代求解公式，通过计算目标滤波窗与所述目标滤波窗下图像梯度的乘积，确定投影恢复模型；
根据滤波窗的加权均值和所述投影恢复模型，计算无噪声干扰的理想投影弦图；
其中，所述投影恢复模型的表达式为：
其中，
为待恢复的无损投影图像；qi为经过重排的投影弦图；
为噪声模型的一维表示；m为总像素数量；β为松弛因子；Ω为像素i的邻域；γij为纹理权重函数。
可选地，所述对所述理想投影弦图进行解重排，得到恢复的积分图像中，所述解重排的表达式为：
其中，p(u,v,s)代表解重排后恢复的积分图像；
代表无噪声干扰的理想投影弦图；u表示积分图像中投影图像的行；v表示积分图像中投影图像的列；s表示积分图像中第s张投影图像；d为探测器像素大小。
可选地，所述对所述恢复的积分图像进行重建，得到目标CBCT断层图像，包括：
采用经过锐化改造的迭代重建算法对所述恢复的积分图像进行重建；
根据所述重建图像确定目标函数；
其中，所述目标函数的表达式为：
其中，f表示迭代优化得到的重建图像；Pf表示对重建图像f进行正向投影；p为解重排的积分图像；μ为松弛调节因子；R(f)为正则化项；
所述重建图像采用GPU并行计算方式进行加速；所述正则化项通过梯度方差构造而来；所述迭代重建算法采用CGLS共轭梯度下降算法求解得到。
本发明实施例的另一方面还提供了一种三维X射线低剂量成像装置，包括：
第一模块，用于采集空曝数据，估计三维X射线系统内参数；
第二模块，用于根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像；
第三模块，用于对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图；
第四模块，用于根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图；
第五模块，用于对所述理想投影弦图进行解重排，得到恢复的积分图像；
第六模块，用于对所述恢复的积分图像进行重建，得到目标CBCT断层图像。
本发明实施例的另一方面还提供了一种电子设备，包括处理器以及存储器；
所述存储器用于存储程序；
所述处理器执行所述程序实现如前面所述的方法。
本发明实施例的另一方面还提供了一种计算机可读存储介质，所述存储介质存储有程序，所述程序被处理器执行实现如前面所述的方法。
本发明实施例还公开了一种计算机程序产品或计算机程序，该计算机程序产品或计算机程序包括计算机指令，该计算机指令存储在计算机可读存储介质中。计算机设备的处理器可以从计算机可读存储介质读取该计算机指令，处理器执行该计算机指令，使得该计算机设备执行前面的方法。
下面结合说明书附图，对本发明的具体实现原理进行详细说明：
本发明的具体操作步骤如附图1所示。所述的发明共分为6个步骤，以下分别对此6个步骤进行详细解释说明：
101：采集空曝数据，估计系统内参数。本发明实施例所使用的三维低剂量成像设备的几何参数如下：射线源到旋转中心距离为100cm，旋转中心到探测器的距离为50cm，激光中心到成像中心的距离为55cm。为估计系统的内参数，需要在固定角度下重复采集若干张空曝数据。如图2所示，其中图2(a)为本实施例采集的投影序列示意图，此序列的每个投影像素可用符号I(u,v,s)表示，其中u表示投影图像的行，v表示投影图像的列，s表示第s张投影图像，s≥500。为准确估计系统内参数，本实施例中重复采集的投影序列设定为1000张。由于扫描的是空气，可以认为射线不会发生余晖、散射及硬化的物理现象，而仅受到噪声的干扰。因此空气的积分图像可用以下公式计算：
所有位于(u,v)位置的像素可提取构成1000个元素的数组p(u,v,:)，其中：表示所有张。根据射线物理知识可知，粒子穿透物体时具有随机性质，数组p(u,v,:)的数值分布符合非平稳高斯分布，如图2(b)所示。根据此物理特性，本实施例采用以下公式对噪声建模：
其中，
为数组p(u,v,:)的均值；σ2(u,v)为数组p(u,v,:)的方差；g(u,v)为调节函数；an为尺度系数；N为尺度系数的阶数。
本实施例中的N选择为2阶尺度。计算g(u,v)和an的过程可描述如下：对噪声模型公式两边取对数，可得：
对上式进行多项式拟合，即可得到logg(u,v)和an的值，亦即得到g(u,v)和an的值。遍历所有位置的像素，即可得到与投影相同大小的二维调节函数矩阵及尺度系数矩阵。
102：对病人原始投影数据余晖校准、散射校准及硬化校准，获取积分图像。三维X射线系统通常采用平板探测器，这类探测器对光子的捕获率较低，容易产生余晖，并且无法使用滤线栅对散射线进行滤除，采集的投影图像质量相比CT类的线阵探测器会差很多。此外，三维X射线系统的射束硬化现象也比CT设备严重得多。上述的余晖、散射及射束硬化可视为一种噪声，会造成101所建立的噪声模型偏移和干扰，准确性下降，因此有必要在应用投影恢复模型前对上述的偏移进行校准。校准的先后顺序为先执行余晖校准，再执行散射校准，最后对前两者修正的数据进行硬化校准。因余晖校准技术及硬化校准技术在本行业内已经非常成熟，本发明不再进行展开叙述。本实施例主要叙述基于蒙特卡洛模拟的散射校正技术。蒙特卡洛方法可精准模拟X射线光子从射线源出射后，穿透物体并发生相互作用(康普顿、瑞利和光电作用)，最终到达探测器并沉积能量的整个物理过程，是散射分布估计的金标准。本实施例所使用的蒙特卡洛模拟算法为自主研发的算法，并经过实验验证结果与开源蒙特卡洛代码包(如PENELOPE、EGSnrc、MCNP和GEANT4等)一致。算法需要输入病人的先验三维容积图像作为发生相互作用的物体，本实施例先验数据的来源为病人螺旋CT容积数据，若病人无螺旋CT容积数据，则选择为FDK解析重建的三维容积数据。假设S(u,v)为蒙特卡洛模拟得到的散射分布，Ip(u,v)为余晖校准的投影图像，则散射校准后的投影可描述为：
I′p(u,v)＝Ip(u,v)-λS(u,v)
其中，I′p(u,v)为校准后的投影，λ为放大系数，用来匹配蒙特卡洛模拟数据的幅度与真实数据的幅度。
103：对积分图像数据重排，得到噪声污染的低剂量投影弦图。为更高效恢复低剂量投影数据，需要对102得到的积分图像进行数据重排。数据重排公式可描述为：
其中，d为探测器像素大小，D为射线源到探测器的距离。数据重排前后的图像如图3(a)和(b)所示。以上投影数据重排为三维伪扇形束数据，期间涉及的插值算法选用线性插值。也可将投影数据重排为三维伪平行束数据，但会存在数据冗余丢失。
104：应用投影恢复模型，估计无噪声干扰的理想投影弦图。根据贝叶斯理论，从带噪图像恢复无损图像的过程可用最大后验概率去描述以下优化模型：
其中，
为待恢复的无损投影图像。qi为经过重排的投影弦图，此处用一维索引表示。
为101步骤所定义的噪声模型的一维表示。m为总像素数量，β为松弛因子。Ω为像素i的邻域，本实施例中选择为4邻域。γij为纹理权重函数，可根据图像的纹理进行自适应调节参数的值，实现去除噪声的同时尽可能保留投影的细节信息。
称为模型的数据保真项，
称为模型的惩罚项。对上述投影恢复模型进行求导，可得到上述优化的迭代求解公式：
其中，γij可拆解为一个滤波窗与该窗下的图像梯度的乘积，本实施例中所使用的一个滤波窗为：
图像梯度的计算方法为：
由于
的原始计算需要多次采集数据以获得数据均值，但实际诊断中不可能多次采集，因此采用滤波窗的加权均值
来替换101噪声模型中所定义的
的计算如图4所示，可描述为此均值不仅考虑当前滤波窗的影响，还考虑不同层的弦图之间的影响。用数学公式可表示为：
其中，ω(s′)＝[0.25 0.5 1 0.5 0.25]为不同层弦图的权重。
可尽可能的靠近理想均值
105：对理想投影弦图解重排，得到恢复的积分图像。经过步骤104后，低剂量采集的带噪弦图已被去除大部分随机噪声的干扰，可认为目前的弦图已趋近于理想投影弦图。对此弦图进行解重排，即可得到恢复的积分图像。解重排的数学公式可描述如下，与103步骤的重排操作互为逆运算：
其中，符号d与D的定义与步骤103一致。
106：采用改进加速的迭代重建算法对恢复的积分图像进行重建，得到噪声低、分辨率高的优质CBCT断层图像。由于迭代重建自身的良好抗噪性能，在图像重建过程中不会像解析重建一样引入重建算法导致的噪声。但由于多次迭代，重建图像会趋于模糊，组织结构的边缘会不同程度被弱化。为了尽可能提高重建图像的图像分辨率，本实施例采用经过锐化改造的迭代重建算法，算法的优化过程可用以下公式描述：
其中，f表示迭代优化得到的重建图像，Pf表示对重建图像f进行正向投影，p为解重排的积分图像。μ为松弛调节因子，R(f)为正则化项。通常R(f)由重建图像的全变差(TV)计算得到，但TV项仍不足以改善弱化的边缘。本实施例使用梯度方差来对R(f)进行构造，其数学表达式描述如下：
其中，
为G的均值。上述优化函数的求解采用CGLS共轭梯度下降算法，可比常规的GDLS梯度下降算法更快收敛。因CGLS算法已有众多文献公开，本发明不再赘述也不作为保护范围。为了进一步加速迭代重建的运算，本实施例采用GPU并行计算方式进行迭代重建，可在1分钟之内完成512x512x300的重建矩阵计算。经过投影恢复及改进的迭代重建算法，即可得到噪声低、分辨率高的优质CBCT断层图像。如附图5所示，(a)图为未使用本发明方案所得到的重建图像，(b)图为采用本发明方案所得到的重建图像，可明显看到本发明方法对图像效果的改进效果。附图6所示的(a)图为基于TV的迭代重建图像，(b)图为本发明方法的迭代重建图像，也可看出在平坦区域本发明方法的重建图像均匀性及锐利度更好，说明了本发明方法不仅可有效解决低剂量成像导致的噪声问题，还可尽可能地保持组织结构本身的纹理特性。
综上所述，相较于现有技术，本发明的三维X射线系统低剂量成像方法有以下优点：
(1)目前尚未有很好的CBCT低剂量成像解决方案，本发明提供的方法可有效解决当前低剂量三维图像噪声大、图像分辨率差的技术痛点。
(2)本发明提供的技术方案可充分考虑平板探测器的余晖效应、射线散射及硬化效应对投影分布模型造成的干扰，使得投影恢复模型更加精准。
(3)本发明提供的迭代重建算法为经过锐化改造及加速优化的并行迭代重建算法，图像的空间分辨率更高，同时重建速度较其他迭代算法更快。
在一些可选择的实施例中，在方框图中提到的功能/操作可以不按照操作示图提到的顺序发生。例如，取决于所涉及的功能/操作，连续示出的两个方框实际上可以被大体上同时地执行或所述方框有时能以相反顺序被执行。此外，在本发明的流程图中所呈现和描述的实施例以示例的方式被提供，目的在于提供对技术更全面的理解。所公开的方法不限于本文所呈现的操作和逻辑流程。可选择的实施例是可预期的，其中各种操作的顺序被改变以及其中被描述为较大操作的一部分的子操作被独立地执行。
此外，虽然在功能性模块的背景下描述了本发明，但应当理解的是，除非另有相反说明，所述的功能和/或特征中的一个或多个可以被集成在单个物理装置和/或软件模块中，或者一个或多个功能和/或特征可以在单独的物理装置或软件模块中被实现。还可以理解的是，有关每个模块的实际实现的详细讨论对于理解本发明是不必要的。更确切地说，考虑到在本文中公开的装置中各种功能模块的属性、功能和内部关系的情况下，在工程师的常规技术内将会了解该模块的实际实现。因此，本领域技术人员运用普通技术就能够在无需过度试验的情况下实现在权利要求书中所阐明的本发明。还可以理解的是，所公开的特定概念仅仅是说明性的，并不意在限制本发明的范围，本发明的范围由所附权利要求书及其等同方案的全部范围来决定。
所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
在流程图中表示或在此以其他方式描述的逻辑和/或步骤，例如，可以被认为是用于实现逻辑功能的可执行指令的定序列表，可以具体实现在任何计算机可读介质中，以供指令执行系统、装置或设备(如基于计算机的系统、包括处理器的系统或其他可以从指令执行系统、装置或设备取指令并执行指令的系统)使用，或结合这些指令执行系统、装置或设备而使用。就本说明书而言，“计算机可读介质”可以是任何可以包含、存储、通信、传播或传输程序以供指令执行系统、装置或设备或结合这些指令执行系统、装置或设备而使用的装置。
计算机可读介质的更具体的示例(非穷尽性列表)包括以下：具有一个或多个布线的电连接部(电子装置)，便携式计算机盘盒(磁装置)，随机存取存储器(RAM)，只读存储器(ROM)，可擦除可编辑只读存储器(EPROM或闪速存储器)，光纤装置，以及便携式光盘只读存储器(CDROM)。另外，计算机可读介质甚至可以是可在其上打印所述程序的纸或其他合适的介质，因为可以例如通过对纸或其他介质进行光学扫描，接着进行编辑、解译或必要时以其他合适方式进行处理来以电子方式获得所述程序，然后将其存储在计算机存储器中。
应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。
以上是对本发明的较佳实施进行了具体说明，但本发明并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可做出种种的等同变形或替换，这些等同的变形或替换均包含在本申请权利要求所限定的范围内。
Claims (10)
1.一种三维X射线低剂量成像方法，其特征在于，包括：
采集空曝数据，估计三维X射线系统内参数；
根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像；
对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图；
根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图；
对所述理想投影弦图进行解重排，得到恢复的积分图像；
对所述恢复的积分图像进行重建，得到目标CBCT断层图像。
2.根据权利要求1所述的一种三维X射线低剂量成像方法，其特征在于，所述采集空曝数据，估计三维X射线系统内参数，包括：
配置三维低剂量成像设备的几何参数，所述几何参数包括射线源到旋转中心距离、旋转中心到探测器的距离以及激光中心到成像中心的距离；
在预设角度下利用所述三位低剂量成像设备重复采集多张空曝数据；
根据所述空曝数据中投影序列，确定与投影相同大小的二维调节函数矩阵及尺度系数矩阵；
所述调节函数及尺度系数的表达式为：
其中，
为数组p(u,v,:)的均值；σ2(u,v)为数组p(u,v,:)的方差；g(u,v)为调节函数；an为尺度系数；N为尺度系数的阶数。
3.根据权利要求1所述的一种三维X射线低剂量成像方法，其特征在于，所述根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像，包括：
对所述原始投影数据进行余晖校准；
通过蒙特卡洛方法对所述余晖校准后的数据进行散射校准；
对所述散射校准后的数据进行硬化校准；
根据所述余晖校准、所述散射校准和所述硬化校准的结果，确定积分图像。
4.根据权利要求1所述的一种三维X射线低剂量成像方法，其特征在于，所述对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图这一步骤中，
所述数据重排的结果包括三维伪扇形束数据或三维伪平行束数据；
所述数据重排的表达式为：
其中，u表示积分图像中投影图像的行；v表示积分图像中投影图像的列；s表示积分图像中第s张投影图像；d为探测器像素大小；D为射线源到探测器的距离；p(s,v,u)代表积分图像数据；q(u′,v′,s′)代表重排后的三维伪扇束数据。
5.根据权利要求1所述的一种三维X射线低剂量成像方法，其特征在于，所述根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图，包括：
通过最大后验概率描述投影恢复模型；
对所述投影恢复模型进行求导，得到优化的迭代求解公式；
根据所述迭代求解公式，通过计算目标滤波窗与所述目标滤波窗下图像梯度的乘积，确定投影恢复模型；
根据滤波窗的加权均值和所述投影恢复模型，计算无噪声干扰的理想投影弦图；
其中，所述投影恢复模型的表达式为：
其中，
为待恢复的无损投影图像；qi为经过重排的投影弦图；
为噪声模型的一维表示；m为总像素数量；β为松弛因子；Ω为像素i的邻域；γij为纹理权重函数。
6.根据权利要求1所述的一种三维X射线低剂量成像方法，其特征在于，所述对所述理想投影弦图进行解重排，得到恢复的积分图像中，所述解重排的表达式为：
其中，p(u,v,s)代表解重排后恢复的积分图像；
代表无噪声干扰的理想投影弦图；u表示积分图像中投影图像的行；v表示积分图像中投影图像的列；s表示积分图像中第s张投影图像；d为探测器像素大小。
7.根据权利要求1所述的一种三维X射线低剂量成像方法，其特征在于，所述对所述恢复的积分图像进行重建，得到目标CBCT断层图像，包括：
采用经过锐化改造的迭代重建算法对所述恢复的积分图像进行重建；
根据所述重建图像确定目标函数；
其中，所述目标函数的表达式为：
其中，f表示迭代优化得到的重建图像；Pf表示对重建图像f进行正向投影；p为解重排的积分图像；μ为松弛调节因子；R(f)为正则化项；
所述重建图像采用GPU并行计算方式进行加速；所述正则化项通过梯度方差构造而来；所述迭代重建算法采用CGLS共轭梯度下降算法求解得到。
8.一种三维X射线低剂量成像装置，其特征在于，包括：
第一模块，用于采集空曝数据，估计三维X射线系统内参数；
第二模块，用于根据所述三维X射线系统参数构建的三维X射线系统采集病人的原始投影数据，并对病人的原始投影数据进行校准处理，得到积分图像；
第三模块，用于对所述积分图像进行数据重排，得到噪声污染的低剂量投影弦图；
第四模块，用于根据所述低剂量投影弦图，通过投影恢复模型估计无噪声干扰的理想投影弦图；
第五模块，用于对所述理想投影弦图进行解重排，得到恢复的积分图像；
第六模块，用于对所述恢复的积分图像进行重建，得到目标CBCT断层图像。
9.一种电子设备，其特征在于，包括处理器以及存储器；
所述存储器用于存储程序；
所述处理器执行所述程序实现如权利要求1至7中任一项所述的方法。
10.一种计算机可读存储介质，其特征在于，所述存储介质存储有程序，所述程序被处理器执行实现如权利要求1至7中任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111453783.5A CN114387359A (zh) | 2021-12-01 | 2021-12-01 | 一种三维x射线低剂量成像方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111453783.5A CN114387359A (zh) | 2021-12-01 | 2021-12-01 | 一种三维x射线低剂量成像方法及装置 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN114387359A true CN114387359A (zh) | 2022-04-22 |
Family
ID=81195641
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111453783.5A Pending CN114387359A (zh) | 2021-12-01 | 2021-12-01 | 一种三维x射线低剂量成像方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114387359A (zh) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115778412A (zh) * | 2023-02-09 | 2023-03-14 | 之江实验室 | X光声成像中造影剂剂量的优化方法、装置及存储介质 |
| CN116071450A (zh) * | 2023-02-07 | 2023-05-05 | 深圳扬奇医芯智能科技有限公司 | 鲁棒的低剂量ct成像算法、装置 |
| CN117671074A (zh) * | 2024-02-01 | 2024-03-08 | 赛诺威盛科技(北京)股份有限公司 | 图像重建方法、装置、电子设备及ct成像系统 |
-
2021
- 2021-12-01 CN CN202111453783.5A patent/CN114387359A/zh active Pending
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116071450A (zh) * | 2023-02-07 | 2023-05-05 | 深圳扬奇医芯智能科技有限公司 | 鲁棒的低剂量ct成像算法、装置 |
| CN116071450B (zh) * | 2023-02-07 | 2024-02-13 | 深圳扬奇医芯智能科技有限公司 | 鲁棒的低剂量ct成像算法、装置 |
| CN115778412A (zh) * | 2023-02-09 | 2023-03-14 | 之江实验室 | X光声成像中造影剂剂量的优化方法、装置及存储介质 |
| CN115778412B (zh) * | 2023-02-09 | 2023-04-28 | 之江实验室 | X光声成像中造影剂剂量的优化方法、装置及存储介质 |
| CN117671074A (zh) * | 2024-02-01 | 2024-03-08 | 赛诺威盛科技(北京)股份有限公司 | 图像重建方法、装置、电子设备及ct成像系统 |
| CN117671074B (zh) * | 2024-02-01 | 2024-05-14 | 赛诺威盛科技(北京)股份有限公司 | 图像重建方法、装置、电子设备及ct成像系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Tang et al. | Performance comparison between total variation (TV)-based compressed sensing and statistical iterative reconstruction algorithms | |
| US10274439B2 (en) | System and method for spectral x-ray imaging | |
| US9245320B2 (en) | Method and system for correcting artifacts in image reconstruction | |
| JP2024054204A (ja) | ニューラルネットワークの学習方法、プログラム、医用画像処理方法及び医用装置 | |
| US6507633B1 (en) | Method for statistically reconstructing a polyenergetic X-ray computed tomography image and image reconstructor apparatus utilizing the method | |
| Siltanen et al. | Statistical inversion for medical x-ray tomography with few radiographs: I. General theory | |
| US9159122B2 (en) | Image domain de-noising | |
| US8885903B2 (en) | Method and apparatus for statistical iterative reconstruction | |
| US20030076988A1 (en) | Noise treatment of low-dose computed tomography projections and images | |
| CN114387359A (zh) | 一种三维x射线低剂量成像方法及装置 | |
| US9261467B2 (en) | System and method of iterative image reconstruction for computed tomography | |
| WO2002067201A1 (en) | Statistically reconstructing an x-ray computed tomography image with beam hardening corrector | |
| EP3195265A2 (en) | Iterative image reconstruction with a sharpness driven regularization parameter | |
| US10813616B2 (en) | Variance reduction for monte carlo-based scatter estimation | |
| JP2016152916A (ja) | Ｘ線コンピュータ断層撮像装置及び医用画像処理装置 | |
| Yoon et al. | Simultaneous segmentation and reconstruction: A level set method approach for limited view computed tomography | |
| Xu et al. | Statistical iterative reconstruction to improve image quality for digital breast tomosynthesis | |
| Zhang et al. | PET image reconstruction using a cascading back-projection neural network | |
| WO2013146283A1 (ja) | 画像処理装置及び画像処理方法 | |
| Friot et al. | Iterative tomographic reconstruction with TV prior for low-dose CBCT dental imaging | |
| EP3404618B1 (en) | Poly-energetic reconstruction method for metal artifacts reduction | |
| US20220375038A1 (en) | Systems and methods for computed tomography image denoising with a bias-reducing loss function | |
| Bai et al. | PET image reconstruction: methodology and quantitative accuracy | |
| US20240127500A1 (en) | Projection-domain material decomposition for spectral imaging | |
| Tilley | High-quality computed tomography using advanced model-based iterative reconstruction |
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 |