CN111915695B - 一种基于方程正交化修正的能谱ct多基材料快速迭代分解方法 - Google Patents

一种基于方程正交化修正的能谱ct多基材料快速迭代分解方法 Download PDF

Info

Publication number
CN111915695B
CN111915695B CN202010777292.5A CN202010777292A CN111915695B CN 111915695 B CN111915695 B CN 111915695B CN 202010777292 A CN202010777292 A CN 202010777292A CN 111915695 B CN111915695 B CN 111915695B
Authority
CN
China
Prior art keywords
projection
energy
base material
projection data
energy spectrum
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
CN202010777292.5A
Other languages
English (en)
Other versions
CN111915695A (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.)
Capital Normal University
Original Assignee
Capital Normal 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 Capital Normal University filed Critical Capital Normal University
Priority to CN202010777292.5A priority Critical patent/CN111915695B/zh
Publication of CN111915695A publication Critical patent/CN111915695A/zh
Application granted granted Critical
Publication of CN111915695B publication Critical patent/CN111915695B/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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4023—Scaling of whole images or parts thereof, e.g. expanding or contracting based on decimating pixels or lines of pixels; based on inserting pixels or lines of pixels
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/77—Retouching; Inpainting; Scratch removal
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/90—Determination of colour characteristics
    • 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]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

本发明公开了一种基于方程正交化修正的能谱CT多基材料快速迭代分解方法，该方法包括：步骤1，采用不同的X射线能谱对被测物体扫描，获得被测物体的多色投影数据；步骤2，为被测物体的各基材料密度图像赋初值；步骤3，对各基材料密度图像的估计值进行正投影，并根据X射线能谱信息和物质的质量衰减系数信息，获得多色投影估计值；步骤4，计算多色投影估计值与多色投影数据之间的误差，利用正交修正求解各基材料投影残差；步骤5，获得各基材料的图像残差，并更新被测物体的各基材料密度图像的估计值；步骤6，重复步骤3至5，直到满足终止条件。本发明能够由采集到的多个能谱的多色投影数据重建出被测物体的多个基材料密度图像。

Description

一种基于方程正交化修正的能谱CT多基材料快速迭代分解 方法
技术领域
本发明涉及X射线CT成像技术领域，特别是关于一种基于方程正交化修正的能谱CT多基材料快速迭代分解方法。
背景技术
X射线计算机断层成像技术(X ray Computed Tomography，简称X射线CT)可以在不破坏或损伤物体的情况下呈现物体的内部细节，已广泛应用于医学、生物、工业、材料、古化石、航天等众多领域。传统CT成像理论假设X射线由单一能量的光子组成，忽略了X射线的多色性，因此利用传统单能量CT重建算法重建实采数据时会产生射束硬化伪影，如杯状伪影和条状伪影，严重影响成像质量。
能谱CT使用多个能量下的X射线能谱扫描被测物体，能够测得比传统单能谱CT更多的被测物体信息。利用这些信息可以重建被测物体的等效原子序数和电子密度或是基材料的密度图像，具有更好的物质区分能力，因此有广泛的应用前景。
能谱CT多色投影数据获取方式主要有以下几种：全扫描模式、慢电压切换的扫描模式和基于光子计数型探测器的扫描模式。全扫描模式采用同一套射线源设备进行多次扫描，间隔使用不同能量，单圈扫描过程中电压不变。此种扫描模式可在常规CT系统平台上实现，简单方便，且单圈扫描时可以灵活设置电压，但需要进行多次扫描。慢电压切换的扫描模式是指单圈扫描过程中，射线源电压随着角度呈规律变换(一般假设该变换为线性的)。在一个变换周期内，射线源电压可以从最低到最高能量或从最高到最低能量完成一次完整的过渡，因此单圈扫描一次即可获得多个能谱下的真实多色投影数据。但是其获得的多色投影数据是不一致的，需要进行配准校正。基于光子计数型探测器的扫描模式采用特殊的光子计数型探测器，其可以同时探测不同的能谱，获得的多色投影数据是几何一致的。但是由于技术不完备，光子计数型探测器存在电荷共享和计数率不足等缺点。
现有的能谱CT重建方法大致可分为三类，标定方法、深度学习方法与迭代方法。标定方法通过扫描标定模体，建立基材料投影与多色投影或基材料图像与多色投影直接重建获得的图像之间的映射关系，然后利用映射关系进行求解，包括直接标定方法与查找表方法。标定方法不需要获取X射线的能谱与物质的质量衰减系数，且一般采用低次多项式进行标定、运算速度快。但是其需要额外扫描标定模体，且获得的映射关系依赖于标定模体的形状与材质。深度学习方法是近年来的研究热点，在有完备集的情况下可以获得高质量的图像重建结果。然而很多情况下，我们无法获得充足的训练样本，并且深度学习方法训练周期较长。迭代法是求解能谱CT问题最常用的方法，其利用数值方法或者优化方法构造迭代结构，通过对图像重建结果逐步修正，可以获得高精度的被测物体的各基材料密度图像信息。由于能谱CT重建问题的非线性性以及病态性，所以现有的能谱CT重建算法无法快速重建出高质量的基材料图像。
发明内容
本发明的目的在于提供一种基于方程正交化修正的能谱CT多基材料快速迭代分解方法，适用于多种常用的能谱CT扫描模式，可以由采集到的多个能谱的多色投影数据重建出被测物体的多个基材料密度图像。
为实现上述目的，本发明提供一种基于方程正交化修正的能谱CT多基材料快速迭代分解方法，该方法包括以下步骤：
步骤1，采用多个不同的X射线能谱对含有多种基材料的被测物体进行X射线扫描，获得被测物体在各个能谱下的真实多色投影数据；
步骤2，根据步骤1获得的多色投影数据，为被测物体的各基材料密度图像赋初值，作为各基材料密度图像的估计值；
步骤3，对各基材料密度图像的估计值进行正投影，获得各基材料的投影估计值，并根据X射线能谱信息和物质的质量衰减系数信息，获得各个能谱下的多色投影估计值；
步骤4，计算步骤3获得的多色投影估计值与步骤1获得的真实多色投影数据之间的误差，并利用正交修正求解各基材料投影残差；
步骤5，将步骤3获得的各基材料的投影残差进行反投影操作，获得各基材料的图像残差，并更新被测物体的各基材料密度图像的估计值；
步骤6，重复步骤3至步骤5，直到满足终止条件。
进一步地，步骤4具体包括：
步骤4.1，计算预设能量下的未插值的真实多色投影数据与由其获得的多色投影估计值之间的误差，并按照下式(3)更新基材料的投影估计值：
Figure BDA0002618915130000031
Figure BDA0002618915130000032
Figure BDA0002618915130000033
Figure BDA0002618915130000034
式中，
Figure BDA0002618915130000035
表示第k次迭代过程中第i次的计算结果，
Figure BDA0002618915130000036
分别表示第k次迭代得到的第m、t个基材料的投影估计值，pi为步骤1获得的多色投影数据，
Figure BDA0002618915130000037
表示第k次迭代获得的第i个能谱对应的多色投影估计值，E表示能量变量，Si，n(E)为第i个归一化能量谱的离散形式，θm(E)、θt(E)分别表示第m、t个基材料的质量衰减系数的离散形式，h＝1,2,…,M,i＝1,2,…,I，I表示能谱的个数，M表示被测物体的待分解基材料的个数；
步骤4.2，根据下式(4)，利用
Figure BDA0002618915130000038
更新正交修正矩阵：
Figure BDA0002618915130000039
式中，Pi称为正交修正矩阵，初值P1设为单位阵，α为正交修正矩阵的参数；
步骤4.3，计算第i且i≥2个能量下的多色投影估计值与相应的真实多色投影数据之间的误差，按照下式(5)更新基材料的投影估计值：
Figure BDA00026189151300000310
式中，
Figure BDA00026189151300000311
为下降方向，
Figure BDA00026189151300000312
为下降步长，lr为松弛因子；
步骤4.4，重复步骤4.2-4.3，直到遍历完所有能量的多色投影数据，根据下式(6)，获得第k+1迭代的各基材料的投影残差：
Figure BDA0002618915130000041
进一步地，步骤1中，当获得的多色投影数据几何不一致时，则对多色投影数据进行插值，获得几何一致的插值投影数据；
步骤4.3中，在几何不一致的情况下，计算多色投影估计值与插值获得的多色投影数据之间的误差。
进一步地，采用下式(2)获得步骤3中的多色投影估计值：
Figure BDA0002618915130000042
进一步地，步骤3中，采用式(1)获取各基材料的投影估计值：
Figure BDA0002618915130000043
式中，L表示X射线路径，通常情况下，L由射线源和探测器的位置确定；RL表示沿X射线路径L的投影算子，
Figure BDA0002618915130000044
表示第k次迭代获得的第m个基材料密度图像的估计值。
进一步地，步骤5中，根据式(7)获得被测物体的各基材料密度图像的估计值：
Figure BDA0002618915130000045
式中，Δf(k+1)表示各基材料的图像残差，
Figure BDA0002618915130000046
为沿X射线路径的重建算子，λ为图像更新时的松弛因子。
进一步地，步骤2中，为被测物体的各基材料密度图像赋初值的方法包括：以0作为初始图像，即初始图像中的每个像素均为0。
进一步地，步骤1的扫描方式为全扫描模式、慢电压切换的扫描模式、或者基于光子计数型探测器的扫描模式。
进一步地，慢电压切换的扫描模式是指在单圈扫描过程中，射线源电压随着角度呈规律变换，在一个变换周期内，射线源电压从最低到最高能量、或者从最高到最低能量完成一次完整的过渡，以通过一次单圈扫描获得多个能谱下的真实多色投影数据。
进一步地，步骤3中的将各基材料密度图像的估计值进行正投影得到各基材料的投影估计值、步骤4所述由多色投影估计值与真实多色投影数据之间的误差利用正交修正求解得到各基材料的投影残差、以及步骤5所述将各基材料的投影残差进行反投影操作得到各基材料的图像残差，均采用并行方法计算，并基于硬件并行计算平台加速实现。
本发明由于采取以上技术方案，其具有以下优点：
本发明方法适用于多种不同的能谱CT扫描方式，可根据获得的不同能谱下的多色投影数据进行多基材料分解重建，且对能谱CT扫描数据无几何一致的要求。与现有方法相比，本发明方法具有重建图像质量高、收敛速度快、抗噪性能强的优点。
附图说明
图1为本发明实施例提供的基于方程正交化修正能谱CT多基材料快速迭代分解方法的流程图；
图2为本发明实施例提供的用作测试模体的FORBILD胸腔模型图像；
图3a为本发明实施例采用的慢电压切换的扫描模式示意图；
图3b为不同切换速度下射线源的能量随着角度的规律变换示意图；
图4为本发明实施例所用的X射线源分别在140kV和80kV管电压下发出的高能谱和低能谱示意图；
图5a为本发明实施例的高能谱下采集到的测试模体添加初始光子数为106的泊松噪声的多色投影数据图像；
图5b为本发明实施例的低能谱下采集到的测试模体添加初始光子数为106的泊松噪声的多色投影数据图像；
图6a为采用本发明方法的骨基材料密度图像重建的结果图像；
图6b为采用本发明方法的水基材料密度图像重建的结果图像；
图6c为本发明实施例所重建的关于能量为70keV的光子的线性衰减系数图像；
图7为本发明方法与已公开的相关专利文献方法的对比示意图；
图8a为本发明实施例用作测试模体的牙模型中骨基材料图像；
图8b为本发明实施例用作测试模体的牙模型中金基材料图像放大图；
图8c为本发明实施例用作测试模体的牙模型中水基材料图像；
图9为本发明实施例所用的X射线源在不同管电压下发出的能谱示意图；
图10a为图9中的40kV管电压下采集到的测试模体添加初始光子数为106的泊松噪声的多色投影数据图像；
图10b为图9中的85kV管电压下采集到的测试模体添加初始光子数为106的泊松噪声的多色投影数据图像；
图10c为本发明图9中的140kV管电压下采集到的测试模体添加初始光子数为106的泊松噪声的多色投影数据图像；
图11a为采用本发明方法的骨基材料密度图像重建的结果图像；
图11b为采用本发明方法的金基材料密度图像重建结果的放大图像；
图11c为采用本发明方法的水基材料密度图像重建的结果图像；
图11d为本发明实施例重建的关于能量为70keV的光子的线性衰减系数图像。
具体实施方式
下面结合附图和实施例对本发明进行详细的描述。
参阅图1，为本发明提出的一种基于方程正交化修正的能谱CT多基材料快速迭代分解方法的流程图。由图1可以看出，本发明方法包括以下步骤：
步骤1：采用多个不同的X射线能谱对含有多种基材料的被测物体进行X射线扫描，获得被测物体在各个能谱下的真实多色投影数据。当获得的多色投影数据几何不一致时，需要对多色投影数据进行插值，获得几何一致的插值投影数据。
步骤2：根据步骤1获得的多色投影数据，为被测物体的各基材料密度图像赋初值，作为各基材料密度图像的估计值。
步骤3：对各基材料密度图像的估计值进行正投影，获得各基材料的投影估计值，并根据X射线能谱信息和物质的质量衰减系数信息，获得各个能谱下的多色投影估计值。
步骤4：计算步骤3获得的多色投影估计值与步骤1获得的真实多色投影数据之间的误差，并用方程进行正交化修正，快速求解获得各基材料的投影残差。
步骤5：将步骤3获得的各基材料的投影残差进行反投影操作，获得各基材料的图像残差，并更新被测物体的各基材料密度图像的估计值。
步骤6：重复步骤3至步骤5，直到满足终止条件。
在一个实施例中，步骤1中，扫描方式可以采用慢电压切换下的扫描模式，也可以采用全扫描模式，还可以采用基于光子计数探测器的扫描模式，甚至采用其它常用的能谱CT扫描模式。全扫描模式是指采用同一套射线源设备、单圈扫描时设置不同的能量、扫描过程中电压不变、多次扫描得到多个能谱下的多色投影数据。基于光子计数型探测器的扫描模式是采用特殊的光子计数型探测器，其可以同时探测不同的能谱，一次扫描即可得到多个能谱下的多色投影数据。其他常用的能谱CT扫描模式还包括快速kVp切换的扫描模式、基于三明治型探测器的扫描模式等。这些扫描模式同样可以利用本发明方法快速迭代，获得各基材料密度图像的估计值。下文主要以慢电压切换的能谱CT扫描模式为例进行说明，其它扫描模式与此类似，在此不再一一展开描述。
慢电压切换的能谱CT扫描模式是指在单圈扫描过程中，射线源电压随着角度呈规律变换，一般假设该变换呈线性。在一个变换周期内，射线源电压可以从最低到最高能量、或者从最高到最低能量完成一次完整的过渡，因此，单圈扫描一次便可以获得多个能谱下的真实多色投影数据。
慢电压切换的能谱CT扫描模式下，使用多个能谱进行X射线扫描是指射线源在不同电压下发出不同能量的X射线光子，即通过调节射线源电压，使用不同电压扫描被测物体。根据多能谱CT系统中所用探测器不同，获得的X射线多色投影数据维度也将不同。当探测器为包含W个探测器单元的线阵性探测器时，在圆轨迹扇束扫描模式下，在V个扫描角度共采集到V×W个X射线投影数据，其中，每个角度的投影包括N个X射线投影数据，将投影数据按照角度顺序进行排序，组合获得一幅二维X射线多色投影图。当探测器为包含W×H个探测器单元的面阵性探测器时，在圆轨迹扇束扫描模式下，在V个扫描角度共采集到V×W×H个X射线投影数据，其中，每个角度的投影包括W×H个X射线投影数据，将投影数据按照角度顺序进行排序，组合获得一幅三维X射线多色投影图。
在一个实施例中，步骤2中的为被测物体的各基材料密度图像赋初值的方法包括：以0作为初始图像，即初始图像中的每个像素均为0。
在一个实施例中，以
Figure BDA0002618915130000071
表示第k次迭代获得的第m个基材料的投影估计值、
Figure BDA0002618915130000072
表示第k次迭代获得的第m个基材料密度图像的估计值、M表示待分解基材料的个数，则步骤3中的对各基材料密度图像的估计值进行正投影，获得各基材料的投影估计值的公式表示为式(1)：
Figure BDA0002618915130000081
其中，L表示X射线路径，通常情况下，L由射线源和探测器的位置确定；RL表示沿X射线路径L的投影算子，其根据探测器个数沿路径采样点个数和路径确定，为已知矩阵。当然，除了式(1)提供的获得各基材料的投影估计值的方法之外，还可以采用现有公知的其它方法获得。
在一个实施例中，以
Figure BDA0002618915130000082
表示第k次迭代获得的第i个能谱对应的多色投影估计值，I表示能谱的个数，则步骤3中的由各基材料的投影估计值获得各个能谱下的多色投影估计值的过程是非线性的，其计算公式表示为下式(2)：
Figure BDA0002618915130000083
式中，E为能量变量，Si，n(E)为第i个归一化能量谱的离散形式，
Figure BDA0002618915130000084
为第m个基材料的质量衰减系数的离散形式。需要说明的是，除了式(2)提供的获得各个能谱下的多色投影估计值的方法之外，还可以采用现有公知的其它方法获得。
在一个实施例中，以pi(i＝1，2，…，I)表示由步骤1扫描或插值获得的第i个能谱对应的多色投影数据，步骤4具体包括：
步骤4.1：将预设能量(假设为第一个能量)下的未插值的真实多色投影数据视为较为准确的数据，计算该能量下的多色投影估计值与真实多色投影数据之间的误差，并按照下式(3)更新基材料的投影估计值：
Figure BDA0002618915130000085
其中：
Figure BDA0002618915130000086
Figure BDA0002618915130000087
Figure BDA0002618915130000091
式中，
Figure BDA0002618915130000092
表示第k次迭代过程中第i次的计算结果，
Figure BDA0002618915130000093
分别表示第k次迭代得到的第m个基材料的投影估计值，
Figure BDA0002618915130000094
表示第k次迭代得到的第t个基材料的投影估计值，pi为步骤1获得的多色投影数据，
Figure BDA0002618915130000095
表示第k次迭代获得的第i个能谱对应的多色投影估计值，E表示能量变量，Si，n(E)为第i个归一化能量谱的离散形式，n为离散能谱序列的索引，θm(E)表示第m个基材料的质量衰减系数的离散形式，θt(E)表示第t个基材料的质量衰减系数的离散形式，h＝1,2,…,M,i＝1,2,…,I，I表示能谱的个数，M表示被测物体的待分解基材料的个数。
步骤4.2：根据下式(4)，利用
Figure BDA0002618915130000096
更新正交修正矩阵：
Figure BDA0002618915130000097
式中，Pi称为正交修正矩阵，初值P1设为单位阵，α为正交修正矩阵参数，一般设为0或者一个极小值。
步骤4.3：计算第i(i≥2)个能量下的多色投影估计值与相应的真实多色投影数据之间的误差(在几何不一致的情况下，则计算多色投影估计值与插值获得的多色投影数据之间的误差)，按照下式(5)更新基材料的投影估计值：
Figure BDA0002618915130000098
式中，
Figure BDA0002618915130000099
为下降方向，
Figure BDA00026189151300000910
为下降步长，lr为松弛因子，pi为步骤1获得的真实多色投影数据或插值获得的多色投影数据。
步骤4.4：重复步骤4.2-4.3，直到遍历完所有能量的多色投影数据，根据下式(6)，获得第k+1迭代的各基材料的投影残差：
Figure BDA00026189151300000911
需要说明的是，除了步骤4.1至步骤4.4提供的获得各基材料的投影残差的方法之外，还可以采用现有公知的其它方法获得。
在一个实施例中，根据式(7)，步骤5中的将各基材料的投影残差进行反投影操作，获得各基材料的图像残差，进而更新被测物体的各基材料密度图像的估计值：
Figure BDA0002618915130000101
式中，
Figure BDA0002618915130000102
为与上文中的RL对应的沿X射线路径L的重建算子，λ为图像更新时的松弛因子λ∈[0，1]。
实质上，除了采用式(7)提供的更新被测物体的各基材料密度图像的估计值的方法之外，还可以采用现有公知的其它方法获得.
在一个实施例中，步骤6中的终止条件，包括但不限于最大迭代次数，获得的被测物体的各基材料密度图像的估计值是否符合主观判断或所需要求，两次迭代获得的各基材料密度图像估计值的差小于设定的阈值，例如0.001。
在一个实施例中，无论是步骤3中的将各基材料密度图像的估计值进行正投影获得各基材料的投影估计值，还是步骤4中的由多色投影估计值与真实多色投影数据之间的误差利用正交修正求解获得各基材料的投影残差，或是步骤5中的将各基材料的投影残差进行反投影操作获得各基材料的图像残差，均可采用并行方法计算，并基于硬件并行计算平台加速实现。
下面通过一个具体实施例来说明本发明方法的具体实现过程。
图2为本具体实施例中用作测试模体的FORBILD胸腔模型图像。本实施例中采用慢电压切换的CT扫描方式，即在单圈扫描过程中，射线源电压随着角度呈规律变换，扫描过程如图3a所示，A、B表示不同电压处的数据采样点。图3b为射线源设定最高能量140kVp与最低能量80kVp时，不同切换速度下射线源kVp随着角度的规律变换。图3b中的圆点为变换周期为0.5个角度时数据采样点，其连接起来的折线表示射线源电压的变换周期为0.5个角度时射线源kVp随着角度变换的关系；五角星点为变换周期为2个角度时数据采样点，其连接起来的折线表示射线源电压的变换周期为2个角度时射线源kVp随着角度变换的关系。在本实施例中，射线源电压的变换周期为0.5个角度，即在0.5个角度内，射线源电压可以从最低到最高能量或从最高到最低能量完成一次过渡，分别在140kVp和80kVp处进行采样。使用开源软件Spectrum GUI，仿真了GE Maxiray 125球管在80kVp和140kVp下的能谱，其中后者在射线源前添加了1mm的铝滤波器，并对能谱做了归一化处理，如图4所示。选用的基材料分别为水和骨，密度分别为1.0g/cm3和1.8g/cm3，线性衰减系数从美国国家标准技术研究院(NIST)网站获得，并根据能谱的取值对其做了相应的插值。
能谱CT系统的扫描参数为：射线源到探测器的距离为880mm，射线源到转台中心的距离为490mm，探测器采用线阵型，由512个0.2mm的探测单元组成，采样角度为720个，所获得的多色投影数据尺寸为720×512，且两个能量下的多色投影数据是几何不一致的。对投影数据添加初始光子数为106的泊松噪声，采集到的测试模体添加泊松噪声的高、低能多色投影数据如图4所示。
具体实施步骤如下：
1)为被测物体的各基材料密度图像赋初值
Figure BDA0002618915130000111
2)假设经过了k(≥0)次迭代的更新，获得的被测物体的各基材料密度图像为
Figure BDA0002618915130000112
和
Figure BDA0002618915130000113
在第k+1次迭代过程中将被测物体的各基材料密度图像进行正投影操作，获得两种基材料的投影估计值
Figure BDA0002618915130000114
和
Figure BDA0002618915130000115
进而利用X射线能谱信息与物质的质量衰减系数信息生成多色投影估计值
Figure BDA0002618915130000116
3)初始化正交修正矩阵P为单位阵，利用未插值的真实多色投影数据p1，按照式(3)更新基材料的投影估计值；
4)利用3)中所用的
Figure BDA0002618915130000117
按照式(4)更新正交修正矩阵。
5)利用真实多色投影数据p2，插值获得与真实多色投影数据p1几何一致的插值多色投影数据
Figure BDA0002618915130000118
按照式(5)更新基材料的投影估计值。
6)按照式(6)获得第k+1迭代的各基材料的投影残差。
7)对6)中获得的各基材料的投影残差进行反投影操作，获得各基材料的图像残差，进而按照式(7)更新被测物体的各基材料密度图像的估计值。
8)对下一次迭代，重复步骤2)-7)，直到满足终止条件。
选取α＝1E-10、lr＝0.01时进行30次迭代获得的基材料图像，其中图6a是骨基材料密度图像(灰度窗为[0，1.92])，图6b是水基材料密度图像(灰度窗为[0.3，1.25])，图6c是关于能量为70keV的光子的线性衰减系数图像(灰度窗为[0，0.06])。这里线性衰减系数图像是通过将重建的水基图像与水对指定能量的X射线光子质量衰减系数的乘积以及重建的骨基图像与骨对指定能量的X射线光子质量衰减系数的乘积求和获得的。由图5a和图5b结果可看出：本发明方法适用于慢电压切换下的能谱CT多基材料分解重建，可以由采集到的多个能谱的多色投影数据重建出被测物体的多个基材料密度图像，且对能谱CT扫描数据无几何一致的要求。与现有方法相比，本发明方法可以降低硬件花费、减少扫描剂量，具有重建图像质量高、收敛速度快、抗噪性能强的优点。
以下通过一个具体实施例来说明本发明与已公开的相关专利文献的不同。
使用开源软件Spectrum GUI生成，仿真了GE Maxiray 125球管在140kVp下的能谱，在射线源前添加了1mm的铝滤波器，并对能谱做了归一化处理。为了清楚显示不同算法的迭代路径，本发明将140kVp的能谱分为两部分：能谱1为140kV p能谱中1-70kVp的部分，能谱2为140kV p能谱中71-140kVp的部分。选用的基材料分别为水和骨，密度分别为1.0g/cm3和1.8g/cm3，线性衰减系数从美国国家标准技术研究院(NIST)网站获得，并根据能谱的取值对其做了相应的插值。将水基材的投影设为4，骨基材的投影设为1，分别获得能谱1和能谱2的多色投影值。图7为利用本发明方法与已公开的相关专利文献方法求解上述问题的迭代轨迹图。其中，两条黑色加粗虚线表示能谱1和能谱2的多色投影曲线，其利用本发明方法获得的。大括号H和J对应的虚线为选取参数α＝1E-10、lr＝0.5时本发明方法求解上述问题的迭代轨迹图。括号J和之字形折线对应的实线为利用E-ART方法求解上述问题的迭代轨迹图。括号I对应的实线为利用I-FBP方法求解上述问题的迭代轨迹图。可以看出，本发明方法与上述两种已公开的相关专利文献方法的不同之处在于本发明方法在迭代至信赖度较高的解集(即多色投影曲线)上后，每一次迭代会沿着上次迭代的正交方向上步进。在当前解集的正交方向上寻找最优解，以保证每次迭代的解可以落在当前解集上，因此获得的解具有较高的准确性；并且在正交方向上求解时解耦了变量之间的耦合信息，可以加快收敛速度。当lr＝1时，获得的解与I-FBP方法获得的解相同，但是由于噪声等因素的存在，此时获得的解往往准确度不高；本发明方法可以通过调节参数lr，在保证解的高精确度的同时加速解的迭代收敛。
最后通过一个具体实施例来说明本发明可以对三个或更多个能谱下的扫描数据进行，获得三个或更多个基材料的密度图像。
本具体实施例中用作测试模体的模型图像原本是一个牙齿模型，将一颗牙分割出来当做金属材质假牙进行模拟，获得图8a为骨基材料图像，图8b为金骨基材料放大图像，图8c为水基材料图像。本实施例中采用慢电压切换的CT扫描方式，射线源设定最高能量140kVp、最低能量40kVp，射线源电压的变换周期为1个角度，即在1个角度内，射线源电压可以从最低到最高能量或从最高到最低能量完成一次过渡，分别在40kVp、85kVp和140kVp处进行采样。使用开源软件Spectrum GUI，仿真了GE Maxiray 125球管在40kVp、85kVp和140kVp下的能谱，其中在利用140kVp电压进行扫描时在射线源前添加了1mm的铝滤波器，并对能谱做了归一化处理，如图9所示。选用的基材料分别为水、骨和金，密度分别为1.0g/cm3、1.8g/cm3和19.3g/cm3，线性衰减系数从美国国家标准技术研究院(NIST)网站获得，并根据能谱的取值对其做了相应的插值。
能谱CT系统的扫描参数为：射线源到探测器的距离为880mm，射线源到转台中心的距离为490mm，探测器采用线阵型，由512个0.2mm的探测单元组成，采样角度为720个，所获得的多色投影数据尺寸为720×512，且三个能量下的多色投影数据是几何不一致的。对投影数据添加初始光子数为106的泊松噪声，采集到的测试模体添加泊松噪声的多色投影数据如图10所示。
具体实施步骤如下：
1)为被测物体的各基材料密度图像赋初值
Figure BDA0002618915130000131
2)假设经过了k(≥0)次迭代的更新，获得的被测物体的各基材料密度图像为
Figure BDA0002618915130000132
和
Figure BDA0002618915130000133
在第k+1次迭代过程中将被测物体的各基材料密度图像进行正投影操作，获得基材料的投影估计值
Figure BDA0002618915130000134
和
Figure BDA0002618915130000135
进而利用X射线能谱信息与物质的质量衰减系数信息生成多色投影估计值
Figure BDA0002618915130000136
3)初始化正交修正矩阵P为单位阵，利用未插值的真实多色投影数据p1，按照(3)式更新基材料的投影估计值；
4)利用3)中所用的
Figure BDA0002618915130000137
按照式(4)更新正交修正矩阵。
5)利用真实多色投影数据p2，插值获得与真实多色投影数据p1几何一致的插值多色投影数据
Figure BDA0002618915130000138
按照式(5)更新基材料的投影估计值。
6)根据
Figure BDA0002618915130000139
按照式(4)更新正交修正矩阵。
7)利用真实多色投影数据p3，插值获得与真实多色投影数据p1几何一致的插值多色投影数据
Figure BDA00026189151300001310
按照式(5)更新基材料的投影估计值。
8)按照式(6)获得第k+1迭代的各基材料的投影残差。
9)对8)中获得的各基材料的投影残差进行反投影操作，获得各基材料的图像残差，进而按照式(7)更新被测物体的各基材料密度图像的估计值。
10)对下一次迭代，重复步骤2)-9)，直到满足终止条件。
选取α＝1E-10、lr＝0.005时进行30次迭代获得的基材料图像如图11所示，其中图11a是骨基材料密度图像(灰度窗为[0.2，1])，图11b是金基材料密度图像的放大图(灰度窗为[0，1])，图11c是水基材料密度图像(灰度窗为[0，1])，图11d是关于能量为70keV的光子的线性衰减系数图像(灰度窗为[0，0.06])。由图11结果可看出，本发明方法可推广至三个或更多个能谱下的扫描数据，获得三个或更多个基材料的密度图像，且对能谱CT扫描数据无几何一致的要求，具有重建图像质量高、收敛速度快、抗噪性能强的优点。
本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。
本领域普通技术人员可以理解：本实施例采用慢电压切换的方式获得所需数据，但是本发明保护范围并不限于此扫描模式。
本领域普通技术人员可以理解：本实施例获得的两组数据是几何不一致的，但是本发明可以对几何一致的数据进行，获得相应的基材料密度图像。
本领域普通技术人员可以理解：本发明可以对三个或更多个能谱下的扫描数据进行，获得三个或更多个基材料的密度图像。
最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (9)

1.一种基于方程正交化修正的能谱CT多基材料快速迭代分解方法，其特征在于，包括以下步骤：
步骤1，采用多个不同的X射线能谱对含有多种基材料的被测物体进行X射线扫描，获得被测物体在各个能谱下的真实多色投影数据；
步骤2，根据步骤1获得的多色投影数据，为被测物体的各基材料密度图像赋初值，作为各基材料密度图像的估计值；
步骤3，对各基材料密度图像的估计值进行正投影，获得各基材料的投影估计值，并根据X射线能谱信息和物质的质量衰减系数信息，获得各个能谱下的多色投影估计值；
步骤4，计算步骤3获得的多色投影估计值与步骤1获得的真实多色投影数据之间的误差，并利用正交修正求解各基材料投影残差；
步骤5，将步骤4获得的各基材料的投影残差进行反投影操作，获得各基材料的图像残差，并更新被测物体的各基材料密度图像的估计值；
步骤6，重复步骤3至步骤5，直到满足终止条件；
步骤4具体包括：
步骤4.1，计算预设能量下的未插值的真实多色投影数据与由其获得的多色投影估计值之间的误差，并按照下式(3)更新基材料的投影估计值：
Figure FDA0003631509900000011
Figure FDA0003631509900000012
Figure FDA0003631509900000013
Figure FDA0003631509900000014
式中，
Figure FDA0003631509900000015
表示第k次迭代过程中第i次的计算结果，
Figure FDA0003631509900000016
ft (k)分别表示第k次迭代得到的第m、t个基材料的投影估计值，pi为步骤1获得的多色投影数据，
Figure FDA0003631509900000027
表示第k次迭代获得的第i个能谱对应的多色投影估计值，E表示能量变量，Si，n(E)为第i个归一化能量谱的离散形式，θm(E)、θt(E)分别表示第m、t个基材料的质量衰减系数的离散形式，h＝1,2,…,M,i＝1,2,…,I，I表示能谱的个数，M表示被测物体的待分解基材料的个数，N表示X射线投影数据的个数；
步骤4.2，根据下式(4)，利用
Figure FDA0003631509900000021
更新正交修正矩阵：
Figure FDA0003631509900000022
式中，Pi为正交修正矩阵，Pi-1为第i-1个能谱对应的正交修正矩阵，初值P1设为单位阵，α为正交修正矩阵的参数；
步骤4.3，计算第i且i≥2个能量下的多色投影估计值与相应的真实多色投影数据之间的误差，按照下式(5)更新基材料的投影估计值：
Figure FDA0003631509900000023
式中，
Figure FDA0003631509900000024
为下降方向，
Figure FDA0003631509900000025
为下降步长，lr为松弛因子；
步骤4.4，重复步骤4.2-4.3，直到遍历完所有能量的多色投影数据，根据下式(6)，获得第k+1迭代的各基材料的投影残差：
Figure FDA0003631509900000026
2.根据权利要求1所述的基于方程正交化修正的能谱CT多基材料快速迭代分解方法，其特征在于，步骤1中，当获得的多色投影数据几何不一致时，则对多色投影数据进行插值，获得几何一致的插值投影数据；
步骤4.3中，在几何不一致的情况下，计算多色投影估计值与插值获得的多色投影数据之间的误差。
3.如权利要求1至2中任一项所述的基于方程正交化修正的能谱CT多基材料快速迭代分解方法，其特征在于，采用下式(2)获得步骤3中的多色投影估计值：
Figure FDA0003631509900000031
4.根据权利要求3所述的基于方程正交化修正的能谱CT多基材料快速迭代分解方法，其特征在于，步骤3中，采用式(1)获取各基材料的投影估计值：
Figure FDA0003631509900000032
式中，L表示X射线路径，L由射线源和探测器的位置确定；RL表示沿X射线路径L的投影算子，
Figure FDA0003631509900000033
表示第k次迭代获得的第m个基材料密度图像的估计值。
5.根据权利要求4所述的基于方程正交化修正的能谱CT多基材料快速迭代分解方法，其特征在于，步骤5中，根据式(7)获得被测物体的各基材料密度图像的估计值：
Figure FDA0003631509900000034
式中，Δf(k+1)表示各基材料的图像残差，
Figure FDA0003631509900000035
为沿X射线路径的重建算子，λ为图像更新时的松弛因子。
6.根据权利要求1所述的基于方程正交化修正的能谱CT多基材料快速迭代分解方法，其特征在于，步骤2中，为被测物体的各基材料密度图像赋初值的方法包括：以0作为初始图像，即初始图像中的每个像素均为0。
7.根据权利要求1所述的基于方程正交化修正的能谱CT多基材料快速迭代分解方法，其特征在于，步骤1的扫描方式为全扫描模式、慢电压切换的扫描模式、或者基于光子计数型探测器的扫描模式。
8.根据权利要求7所述的基于方程正交化修正的能谱CT多基材料快速迭代分解方法，其特征在于，慢电压切换的扫描模式是指在单圈扫描过程中，射线源电压随着角度呈规律变换，在一个变换周期内，射线源电压从最低到最高能量、或者从最高到最低能量完成一次完整的过渡，以通过一次单圈扫描获得多个能谱下的真实多色投影数据。
9.根据权利要求1所述的基于方程正交化修正的能谱CT多基材料快速迭代分解方法，其特征在于，步骤3中的将各基材料密度图像的估计值进行正投影得到各基材料的投影估计值、步骤4由多色投影估计值与真实多色投影数据之间的误差利用正交修正求解得到各基材料的投影残差、以及步骤5将各基材料的投影残差进行反投影操作得到各基材料的图像残差，均采用并行方法计算，并基于硬件并行计算平台加速实现。
CN202010777292.5A 2020-08-05 2020-08-05 一种基于方程正交化修正的能谱ct多基材料快速迭代分解方法 Active CN111915695B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010777292.5A CN111915695B (zh) 2020-08-05 2020-08-05 一种基于方程正交化修正的能谱ct多基材料快速迭代分解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010777292.5A CN111915695B (zh) 2020-08-05 2020-08-05 一种基于方程正交化修正的能谱ct多基材料快速迭代分解方法

Publications (2)

Publication Number Publication Date
CN111915695A CN111915695A (zh) 2020-11-10
CN111915695B true CN111915695B (zh) 2022-07-08

Family

ID=73287185

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010777292.5A Active CN111915695B (zh) 2020-08-05 2020-08-05 一种基于方程正交化修正的能谱ct多基材料快速迭代分解方法

Country Status (1)

Country Link
CN (1) CN111915695B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114916950B (zh) * 2022-07-21 2022-11-01 中国科学院深圳先进技术研究院 基于多层平板探测器的高空间分辨能谱ct图像重建方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103559729A (zh) * 2013-11-18 2014-02-05 首都师范大学 一种双能谱ct图像迭代重建方法
CN103900931A (zh) * 2012-12-26 2014-07-02 首都师范大学 一种多能谱ct成像方法及成像系统
CN108010098A (zh) * 2017-12-04 2018-05-08 首都师范大学 一种双能谱ct基材料图像迭代重建方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9269168B2 (en) * 2013-03-15 2016-02-23 Carestream Health, Inc. Volume image reconstruction using data from multiple energy spectra

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103900931A (zh) * 2012-12-26 2014-07-02 首都师范大学 一种多能谱ct成像方法及成像系统
CN103559729A (zh) * 2013-11-18 2014-02-05 首都师范大学 一种双能谱ct图像迭代重建方法
CN108010098A (zh) * 2017-12-04 2018-05-08 首都师范大学 一种双能谱ct基材料图像迭代重建方法

Also Published As

Publication number Publication date
CN111915695A (zh) 2020-11-10

Similar Documents

Publication Publication Date Title
US11864939B2 (en) Apparatus and method that uses deep learning to correct computed tomography (CT) with sinogram completion of projection data
CN110189389B (zh) 基于深度学习的双能谱ct投影域基材料分解方法及装置
CN108010098B (zh) 一种双能谱ct基材料图像迭代重建方法
Fessler Hybrid Poisson/polynomial objective functions for tomographic image reconstruction from transmission scans
US7298812B2 (en) Image-based material decomposition
US7391844B2 (en) Method and apparatus for correcting for beam hardening in CT images
EP2727083B1 (en) Iterative image reconstruction
CN111968060B (zh) 一种基于倾斜投影修正技术的多能谱ct快速迭代重建方法
CN110533734B (zh) 基于传统单能ct的多能谱分段稀疏扫描迭代重建方法
US9683948B2 (en) Systems and methods for iterative multi-material correction of image data
CN108010099A (zh) 一种x射线多能谱ct有限角扫描和图像迭代重建方法
Johnston et al. Temporal and spectral imaging with micro‐CT
JP2019188135A (ja) Ａｌｖａｒｅｚ−Ｍａｃｏｖｓｋｉ減衰モデルを使用した断層像再構成におけるＸ線ビームハードニング補正
JP2020512900A (ja) コンピュータ断層撮影における単一エネルギーデータを用いたアーチファクト低減の方法
CN111915695B (zh) 一种基于方程正交化修正的能谱ct多基材料快速迭代分解方法
EP3671646A1 (en) X-ray computed tomography (ct) system and method
Gao et al. Beam hardening correction for middle-energy industrial computerized tomography
CN112304987B (zh) 基于光子计数能谱ct的含能材料等效原子序数测量方法
CN110706299A (zh) 一种用于双能ct的物质分解成像方法
CN111899312B (zh) 一种迭代补偿的有限角度ct投影重建方法
Nguyen et al. BeadNet: a network for automated spherical marker detection in radiographs for geometry calibration
CN110264536B (zh) 一种在平行束超分重建中计算高低分辨率投影关系的方法
CN110731788B (zh) 一种基于双能ct扫描仪对基物质进行快速分解的方法
CN107292846B (zh) 一种圆轨道下不完全ct投影数据的恢复方法
Sun et al. Simultaneous correction of motion and metal artifacts in head CT scanning

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