CN111134709A - 一种多能量ct基材料分解方法 - Google Patents
一种多能量ct基材料分解方法 Download PDFInfo
- Publication number
- CN111134709A CN111134709A CN202010054618.1A CN202010054618A CN111134709A CN 111134709 A CN111134709 A CN 111134709A CN 202010054618 A CN202010054618 A CN 202010054618A CN 111134709 A CN111134709 A CN 111134709A
- Authority
- CN
- China
- Prior art keywords
- energy spectrum
- intensity information
- transmission intensity
- ray transmission
- scattering
- 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.)
- Granted
Links
- 239000000463 material Substances 0.000 title claims abstract description 103
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 68
- 238000000034 method Methods 0.000 title claims abstract description 68
- 238000001228 spectrum Methods 0.000 claims abstract description 144
- 230000005540 biological transmission Effects 0.000 claims abstract description 121
- 230000009467 reduction Effects 0.000 claims abstract description 25
- 238000012545 processing Methods 0.000 claims abstract description 23
- 238000003384 imaging method Methods 0.000 claims abstract description 13
- 230000008030 elimination Effects 0.000 claims abstract description 6
- 238000003379 elimination reaction Methods 0.000 claims abstract description 6
- 238000005457 optimization Methods 0.000 claims description 15
- 238000013507 mapping Methods 0.000 claims description 9
- 230000002829 reductive effect Effects 0.000 claims description 4
- 230000005855 radiation Effects 0.000 abstract description 30
- 238000002591 computed tomography Methods 0.000 description 50
- 239000000126 substance Substances 0.000 description 11
- 230000008569 process Effects 0.000 description 8
- 230000003595 spectral effect Effects 0.000 description 6
- 238000012937 correction Methods 0.000 description 5
- ZOKXTWBITQBERF-UHFFFAOYSA-N Molybdenum Chemical compound [Mo] ZOKXTWBITQBERF-UHFFFAOYSA-N 0.000 description 4
- 238000013170 computed tomography imaging Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 4
- 239000010931 gold Substances 0.000 description 4
- 229910052737 gold Inorganic materials 0.000 description 4
- 229910052750 molybdenum Inorganic materials 0.000 description 4
- 239000011733 molybdenum Substances 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000010521 absorption reaction Methods 0.000 description 2
- 230000000670 limiting effect Effects 0.000 description 2
- 238000002203 pretreatment Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000001131 transforming effect Effects 0.000 description 2
- 239000003570 air Substances 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000004888 barrier function Effects 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000002405 diagnostic procedure Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 239000002355 dual-layer Substances 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 230000005764 inhibitory process Effects 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000011946 reduction process Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/032—Transmission computed tomography [CT]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5205—Devices using data or image processing specially adapted for radiation diagnosis involving processing of raw data to produce diagnostic data
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis 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
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- Pathology (AREA)
- Heart & Thoracic Surgery (AREA)
- High Energy & Nuclear Physics (AREA)
- Physics & Mathematics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Veterinary Medicine (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Pulmonology (AREA)
- Theoretical Computer Science (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
本发明公开了一种多能量CT基材料分解方法，属于辐射成像领域，包括：获取多能谱成像系统的多能能谱信息和扫描物体时探测器探测到的射线透射强度信息；对多能能谱进行降维处理，得到新的多能能谱；对射线透射强度信息进行降维处理消除散射项，得到新的射线透射强度信息；根据新的多能能谱和新的射线透射强度信息进行基材料分解。该方法可以消除散射对多能能谱和投影数据的影响，根据消除散射后的数据进行基材料分解。
Description
技术领域
本发明涉及辐射成像技术领域，特别涉及一种多能量CT基材料分解方法。
背景技术
X射线计算机断层成像(Computed Tomography，CT)是一种目前广泛使用的医学诊断手段，其原理是根据不同物质对X射线的吸收能力不同，通过对探测器上接收到的衰减后X射线信号进行三维图像重建得到人体的断层或三维图像，从而更好地协助医生发现病灶。
射线散射一直是CT成像领域的一大难题。传统的CT重建理论将原始探测器采集到的信号认为是射线直射过程中穿过物质衰减系数的线积分，没有考虑射线散射的影响，然而在实际的CT扫描时必然有散射光子到达探测器上被收集。散射信号的出现会导致原始线积分模型出现明显偏差，大大降低重建图像质量并导致图像中出现散射伪影。
随着CT技术朝着快速成像的方向发展，现代CT系统逐渐趋于增加每个机架旋转覆盖的体积。先进的CT系统通常使用发散的锥形束射线源和大面积平板探测器，由于X射线的物理特性，这导致散射信号被进一步增强。
因此，在CT成像领域中的核心问题之一就是如何去除或估计出射线散射的影响。目前主要的射线散射修正方法可以分为两大类：前处理方法和后处理方法。前处理方法是通过阻止散射到达探测器，降低散射占比，前处理方法主要有：栅格阻挡法、增大成像物体与探测器距离的方法。后处理方法是基于对投影数据中散射贡献的估计，后处理方法主要有：直接测量法、基于软件模拟(蒙卡模拟等)、基于硬件的分解方法等。其中，前处理校正方法可以直接借助阻挡物抑制散射，但其抑制能力有限，无法达到让人满意的水平。而后处理方法虽然已经有很多种对于散射的估计方法，但目前仍然没有一种可靠的、最优的散射估计方法。
多能量CT成像是目前CT领域的研究热点之一，主要包括双能CT和基于光子计数探测器的能谱CT。相比于传统CT，多能量CT利用物质在不同X射线能量下吸收能力的差异，能提供比常规CT更多的影像信息，具有提高图像质量、抑制射束硬化伪影、降低辐射剂量等优势。目前主流的临床双能CT扫描仪包括重复进行两次不同能谱的扫描、双源双探测器、快速千伏切换技术、双层探测器技术等，这些技术都各有其优势与劣势。近年来，多能量CT基材料物质分解方法研究主要集中于图像域方法。
发明内容
本发明旨在至少在一定程度上解决相关技术中的技术问题之一。
为此，本发明的目的在于提出一种多能量CT基材料分解方法，该方法消除散射对投影数据的影响，可进行投影域或图像域的基材料物质分解、稀疏数据图像重建等工作。
为达到上述目的，本发明一方面实施例提出了一种多能量CT基材料分解方法：包括以下步骤：
S1，获取多能谱成像系统的多能能谱信息和扫描待成像物体时探测器探测到的射线透射强度信息；
S2，对所述多能能谱进行降维处理，得到新的多能能谱；
S3，对所述射线透射强度信息进行降维处理，得到新的射线透射强度信息；
S4，根据所述新的多能能谱和所述新的射线透射强度信息对所述待成像物体进行基材料分解。
本发明实施例的多能量CT基材料分解方法，计算多能能谱和射线透射强度信息并进行优化，对射线透射强度信息和多能能谱降维消除散射，根据消除散射后的射线透射强度信息和多能能谱直接进行基材料分解，在多能数据不足的情况下，根据消除散射后的射线透射强度信息和多能能谱进行散射估计，根据散射估计对已有的透射数据进行去散射处理，得到去散射的透射数据，进而生成去散射的投影数据，再进行基材料分解。由此，可以进行投影域或图像域的基材料物质分解、稀疏数据图像重建等工作。
另外，根据本发明上述实施例的多能量CT基材料分解方法还可以具有以下附加的技术特征：
进一步地，在本发明的一个实施例中，还包括：在步骤S4之后还包括：
S5，根据所述新的多能能谱、所述新的射线透射强度信息以及消除散射方程计算基材料分解系数；
S6，根据计算出的所述基材料分解系数与所述多能能谱和所述射线透射强度信息之间的映射关系反解出散射值，根据所述散射值和散射值在邻域分布的相似性得到整体散射估计；
S7，根据所述整体散射估计对所述射线透射强度信息或优化后的射线透射强度信息进行去散射处理，得到去散射后的射线透射强度信息，根据所述去散射后的射线透射强度信息进行基材料分解。
进一步地，在本发明的一个实施例中，在对所述多能能谱和所述射线透射强度信息进行降维处理之前，包括：对所述多能能谱和所述射线透射强度信息进行插值优化处理，以使所述多能能谱和所述射线透射强度信息趋于一致。
进一步地，在本发明的一个实施例中，以所述多能能谱或所述射线透射强度信息作为自变量构建降维函数，通过所述降维函数消除一维的散射，将N维的所述多能能谱和所述射线透射强度信息降为不包含散射信息的N-1维信息。
进一步地，在本发明的一个实施例中，所述降维函数为所述多能能谱两两相减和所述射线透射强度信息两两相减；或
优化后的多能能谱两两相减和优化后的射线透射强度信息两两相减。
进一步地，在本发明的一个实施例中，所述S5包括：
所述消除散射方程参数包括所述新的射线透射强度信息、所述新的多能能谱、多个基材料参数，根据所述新的射线透射强度信息、所述新的多能能谱和部分基材料参数求解所述基材料分解系数，其中，所述多个基材料参数包括所述基材料分解系数。
进一步地，在本发明的一个实施例中，所述映射关系为探测器上接收到的射线透射强度信息与所述多能能谱、基材料衰减系数、等效厚度、对应能谱、散射值的关系式。
进一步地，在本发明的一个实施例中，所述S5中，计算所述基材料分解系数的方法包括但不限于迭代求解、查找表法和多项式拟合法。
进一步地，在本发明的一个实施例中，通过N能能谱和N能透射数据对小于等于N-1种基材料分解，其中，N大于等于3。
进一步地，在本发明的一个实施例中，在所述S6中，所述根据所述散射值和散射值在邻域分布的相似性得到整体散射估计，包括：
利用包括但不限于插值或加权的方式估计出整体的散射分布，完成散射估计。
本发明附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
本发明上述的和/或附加的方面和优点从下面结合附图对实施例的描述中将变得明显和容易理解，其中：
图1为根据本发明一个实施例的一种多能量CT基材料分解方法流程图；
图2为根据本发明一个实施例的多能量CT系统结构示意图；
图3为根据本发明一个实施例的多能量CT系统成像平面图；
图4为根据本发明一个实施例的两种条状板重叠得到的能谱调制器示意图；
图5为根据本发明一个实施例的含散射校正下的三能量CT基材料物质分解方法流程图。
具体实施方式
下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
下面参照附图描述根据本发明实施例提出的一种多能量CT基材料分解方法。
图1为根据本发明一个实施例的一种多能量CT基材料分解方法流程图。
如图1所示，该多能量CT基材料分解方法包括以下步骤：
S1，获取多能谱成像系统的多能能谱信息和扫描待成像物体时探测器探测到的射线透射强度信息。
在本发明的一个实施例中，作为一种实现方式，可以通过下述方式获得多能能谱信息和射线透射强度信息，首先计算未放置待成像物体且射线源焦点位于不同位置时，射线穿过能谱调制器不同子单元后探测器上的多能能谱，再计算放置待成像物体且射线源焦点位于不同位置时，射线穿过能谱调制器不同子单元后探测器上的射线透射强度信息。
可以理解的是，本发明实施例的方法在进行基材料分解时，可以基于图2所示的多能量CT系统，也可以应用于其它的成像系统，以多能量CT系统为例介绍本发明实施例的基材料分解方法。
如图2所示，为本发明实施例中的多能量CT系统结构示意图，包括：射线源，能谱调制器、探测器以及机械/电气控制与数据传输/处理单元。图中的旋转箭头表示成像时射线源、能谱调制器和探测器绕基材料旋转。
首先在多能量CT系统中，不放置待成像物体，计算射线源焦点位于不同位置时射线穿过能谱调制器不同子单元后探测器上的多能能谱，然后在多能量CT系统放置待成像物体，射线源在不同焦点位置时发出的射线穿过能谱调制器和待成像物体后，采集探测器上的射线透射强度信息。
如图2所示，射线源为具有变焦点功能的射线源，本实施例中使用的多能量数据是通过变换射线源焦点位置，穿过一块设计好的网格状能谱调制器，使得探测器上得到分别穿过不同衰减的子单元的不同能谱射线。且两种能量(图3中实线、虚线)的透射数据仅相隔一个探测器单元，在物体上的投影线距离很小，可视为基本一致。在CT扫描旋转时，由于旋转间隔较小，同一探测器单元上旋转角度连续的不同能量的透射数据也可视作同一角度下的多能量数据。
如图3所示，在多能量CT系统中，由于射线源焦点移动距离很小，且被扫描物体、能谱调制器、探测器均不移动，因此在不同能谱下的散射分布基本一致，某方向移动时的平面示意图如图3所示，因此有Is1≈Is2≈Is3≈Is＝const。Isi为在能谱Si(E)下的散射强度。
在本发明的实施例中，根据CT投影理论，多能量下的能谱{Si(E),i＝1,2,3}与投影到探测器上的射线透射强度信息{Pi,i＝1,2,3}关系可以表示为：
其中，μ1(E)和μ2(E)分别为两种基材料对应的衰减值，L1和L2为两种基材料等效厚度积分(即加权系数投影值)，Isi为在能谱Si(E)下的散射强度。
步骤S2，对多能能谱进行降维处理，得到新的多能能谱。
步骤S3，对射线透射强度信息进行降维处理，得到新的射线透射强度信息。
步骤S4，根据新的多能能谱和新的射线透射强度信息对待成像物体进行基材料分解。
可以理解的是，根据成像扫描视场中有物体和没有物体时的射线透射数据，通过取负对数等常规数据处理，得到投影数据。然后根据新的多能能谱和新的投影数据进行基材料分解。
进一步地，在本发明的一个实施例中，在对多能能谱和射线透射强度信息进行降维处理之前，包括：对多能能谱和射线透射强度信息进行插值优化处理，以使多能能谱和射线透射强度信息趋于一致。
可以理解的是，在步骤S1得到的多能能谱和射线透射强度信息为N维的信息，需要求解一维的散射和N-1维的基材料分解系数，需要对N维多能能谱和射线透射强度信息进行降维操作，去除散射信息，留下只包含基材料分解系数的信息。
对N维信息进行降维具有多种方式，在本发明的实施例中，通过降维函数对N维信息进行降维，将N维的所述多能能谱和所述射线透射强度信息降为N-1维，其中，作为一种实现方式，降维函数以多能能谱或射线透射强度信息作为自变量。
其中，作为一种可能实现的方式，降维函数的具体形式为：降维函数为多能能谱两两相减和射线透射强度信息两两相减；或
优化后的多能能谱两两相减和优化后的射线透射强度信息两两相减。
具体地，降维函数的形式可以为：
Sij＝gij(Si(E)-Sj(E))或Sij＝Si(E)-Sj(E)
Pij＝fij(Pi-Pj)或Pij＝Pi-Pj
其中，Sij为新的多能能谱，Si(E)和Sj(E)为多能能谱，gij表示对多能能谱进行优化，Pij为新的射线透射强度信息，Pi和Pj为射线透射强度信息，fij表示对射线透射强度信息进行优化。
可以理解的是，在降维函数中，可以利用优化前的多能能谱或射线透射强度信息，也可以利用优化后的多能能谱或者射线透射强度信息。在降维处理之前进行如插值或其他优化处理可以使得多能能谱数据更加趋于一致。
具体地，以在X、Z方向变动射线源焦点位置为例，在探测器上得到三种能谱的透射数据后，如图3所示，原有透射射线实际未对准，可以对其进行处理，如进行插值或其他优化处理使得两种能谱的透射数据对应的射线更加趋于一致，然后两两相减，可以消去散射项，得到新的透射数据结果如下：
Pij＝fij(Pi-Pj)＝∫[gij(Si(E)-Sj(E))]exp(-μ1(E)L1-μ2(E)L2)dE
其中，fij(Pi-Pj)是对X进行根据投影数据位置的插值，或者是其他对Pi或Pj优化处理的结果，gij(Si(E)-Sj(E))是相应的对能谱的处理。此处也可以对Pi、Pj、Si(E)和
Sj(E)中的部分或全部进行优化处理，得到效果与对Pj、Sj(E)处理效果相同。fij(Pj)、gij(Sj(E))存在一个特例，该特例为不进行优化处理的情况，如下式：
总的来说，fij(Pi-Pj)、gij(Si(E)-Sj(E))、fij(Pj)、gij(Sj(E))、fij(Pj)、gij(Sj(E))是对于透射数据、能谱的优化处理，使得两种能谱和相应透射数据对应的射线更加趋于一致。如图3所示，原有投影射线实际未对准，且存在其他导致射线偏差的因素，优化后，数据更加趋于一致。
本发明实施例的中使用的能谱调制器包括多个单元，多个单元包括多个子单元，每个子单元具有不同的衰减，由于三种能谱是由射线穿过不同衰减的子单元后得到，由X射线特性可知，通常最后得到的高能能谱的数据衰减最多，低能能谱的数据衰减最少，即P1>P2>P3，其中1,2,3指代高能、中能、低能的透射数据。因此可以将Pij视作新的去除了散射部分的透射数据，gij(Si(E)-Sj(E))视作新的能谱，用已有的三种能谱数据，可以得到新的去除散射的两种能谱的CT透射数据。然后利用去除了射线散射部分的两种新能谱的CT透射数据，生成新的投影数据，可以直接进行常规的双能CT基材料物质分解。
在进行双能CT的基材料物质分解时，通常需要两种能谱有足够大的能量差异，即两种能谱的平均能量差。因此，存在N种能量的数据不足，在数据不足的情况下，还可以包括：
S5，根据新的多能能谱、新的射线透射强度信息以及消除散射方程计算基材料分解系数；
S6，根据计算出的基材料分解系数与多能能谱和射线透射强度信息之间的映射关系反解出散射值，根据散射值和散射值在邻域分布的相似性得到整体散射估计；
S7，根据整体散射估计对射线透射强度信息或优化后的射线透射强度信息进行去散射处理，得到去散射后的射线透射强度信息，根据去散射后的射线透射强度信息进行基材料分解。
可以理解的是，在进行双能CT的基材料物质分解时，通常需要两种能谱有足够大的能量差异，即两种能谱的平均能量差。在用图2所示多能量CT系统时，可能无法令所有探测器像素处新的能谱gij(Si(E)-Sj(E))都有足够大的能量差异。因此，也可以只利用部分相减后新的能谱gij(Si(E)-Sj(E))的能量差异足够大的数据，利用这部分新的去除散射的去除散射的两种能谱的CT透射数据进行基材料物质分解后，重新根据公式
计算散射，利用散射在邻域分布的相似性，得到整体的散射估计，然后进行含散射校正的基材料物质分解或CT重建。
根据数值模拟可以得到当所用系统中的能谱调制器为特定不同衰减的子单元时，相减后进而得到新的能谱gij(Si(E)-Sj(E))的平均能量。如图4所示，在此以用两种条状板重叠得到的能谱调制器为例，其中，无颜色区域代表镂空区域，有颜色区域为特定材料。分别取两种条状板材料为钼、金，钼厚度为0.3mm，金厚度为0.0253mm，则对应的四个子单元分别为：空气，钼，金，钼+金。通过数值计算得出垂直穿过能谱调制器不同子单元得到的其中三种能谱按上述方法相减后，两种新的能谱的平均能量差最高为20.83keV，具备进行双能物质分解的可行性。
进一步地，消除散射方程参数包括新的射线透射强度信息、新的多能能谱、多个基材料参数，根据新的射线透射强度信息、新的多能能谱和部分基材料参数求解基材料分解系数，其中，多个基材料参数包括基材料分解系数。
通过上述步骤得到了新的多能能谱Sij，新的射线透射强度信息Pij，再根据消除散射方程消除散射并得到基材料分解系数。
其中，消除散射方程可以为上述提到的：
Pij＝fij(Pi-Pj)＝∫[gij(Si(E)-Sj(E))]exp(-μ1(E)L1-μ2(E)L2)dE
其中，Pij为新的射线透射强度信息，i、j为正整数，fij为对射线透射强度信息进行优化，gij为对多能能谱进行优化，L1和L2为两种基材料等效厚度积分，μ1(E)和μ2(E)分别为两种基材料对应的衰减值。在上式中，基材料分解系数为L1和L2。
进一步地，在本发明的一个实施例中，计算基材料分解系数的方法包括但不限于迭代求解、查找表法和多项式拟合法。
进一步地，在本发明的一个实施例中，基材料分解系数与多能能谱和射线透射强度信息之间的映射关系为探测器上接收到的射线透射强度信息与多能能谱、基材料衰减系数、等效厚度、对应能谱、散射值的关系式。
其中，映射关系可以为：
其中，Si(E)为多能能谱，Pi为射线透射强度信息，μ1(E)和μ2(E)分别为两种基材料对应的衰减值，L1和L2为两种基材料等效厚度积分，Isi为在能谱Si(E)下的散射值。
在步骤S7中，计算出整体的散射估计后，根据散射估计对原有的射线透射强度信息或者优化后的射线透射强度信息直接进行去散射处理，根据去散射后的射线透射强度信息进行基材料分解。
下面以三能量情况下的基材料分解方法对本发明实施例的分解过程进行说明，通过上述推导过程，消去射线透射强度信息中的散射项，得到如下方程：
P12＝f12(P1-P2)＝∫[g12(S1(E)-S2(E))]exp(-μ1(E)L1-μ2(E)L2)dE
P13＝f13(P1-P3)＝∫[g13(S1(E)-S3(E))]exp(-μ1(E)L1-μ2(E)L2)dE
需要说明的是，此处不一定为1-2,1-3，也可用其他组合方式，例如1-3,2-3等实现类似效果，1,2,3指代不同的能谱。
在方程中，P1,P2,P3,S1(E),S2(E),S3(E),μ1(E),μ2(E)均已知，需求解L1,L2，由两个方程求解两个未知量，可采用多种方法进行求解，例如迭代求解、查找表法、多项式拟合法等都有充分的研究，可直接进行基材料物质分解或图像重建。也可只利用部分多能数据，利用求得的L1，L2，可根据映射公式重新计算散射，映射公式中Pi，Si(E)可以是fij(x)、gij(x)进行优化处理后的结果，也可以直接用未优化结果或部分用未优化结果进行计算。利用散射在邻域分布的相似性，得到整体的散射估计，然后结合部分双能数据或单能数据进行含散射校正的基材料物质分解或稀疏数据重建。
上述具体实施例为三能量的两种基材料的分解方法，在实际应用中，对能谱调制器的结构进行优化后，变动射线源的焦点位置，可以实现在N能量的情况下，最多进行N-1种基材料分解(N≥3)，也可以对散射分布进行更加精确的计算，比如在四能量下可以分别用三种能量数据进行两种基材料的分解，然后对不同组合求出的散射进行平均或其他加权方式得到更为精确的散射分布。
如图5所示，展示了含散射校正下的三能量CT基材料物质分解流程。具体为：
S0：在本发明所用的多能CT成像系统中(示意图如图2所示)，多能谱数据是通过变换射线源焦点位置，使得探测器上得到分别穿过不同衰减的子单元的射线而获得，先建立不含能谱调制器时射线源焦点处于不同位置时的有效能谱模型。
S1：根据测量的能谱调制器的不同子单元厚度，结合几何关系、能谱调制器的不同子单元的线衰减系数、射线源焦点位置等信息，通过数值计算求得含能谱调制器的CT系统有效能谱。
S2：按之前所述，对三种能谱进行相减，在相减前可以先对能谱进行优化，如进行插值及其它优化方法使不同能谱的投影射线更加趋于一致。
S3：如图2所示系统中，在CT扫描时变动射线源焦点位置，得到不同能谱的透射数据。此数据可以是部分为三种能量，部分为两种能量、部分为一种能量。
S4：由于源变动距离很小，且被成像物体、能谱调制器、探测器均不移动，旋转基本不影响散射分布，因此在焦点变换前后，光子散射分布具有一定相似性。以X、Z方向变动射线源位置为例，在探测器上得到三种能谱的透射数据后，对其进行优化处理，如插值使透射数据位置更为趋于一致，然后两两相减，可以消去散射项。
S5：利用新的去除散射的两种能谱的CT透射数据，按上述推导，散射项消去。对新的双能能谱和双能透射数据，可以直接进行常规的双能CT的基材料物质分解，还可以根据新的双能能谱和双能透射数据计算基材料系数。
在这一步骤中，根据成像扫描视场中有物体和没有物体时的射线透射数据，通过取负对数等常规数据处理，得到投影数据。然后根据新的多能能谱和新的投影数据进行基材料分解。
S6：根据计算得到的基材料系数，反推出探测器单元上的散射大小。利用散射项在像素邻域基本不变或相近的特性，通过不完整的散射分布，利用插值、加权等方式估计出整体的散射分布，完成散射估计。
S7：根据散射估计，对已有的数据进行去散射处理，得到去散射的透射数据。根据原有的双能数据或单能数据，可进行投影域或图像域的基材料物质分解、稀疏数据图像重建等工作。
直接通过低能量下重建图像，会含有大量散射造成的伪影，通过图5所示的流程采用一种迭代求解方法进行基材料物质分解，可以有效消除散射对分解结果的影响。
根据本发明实施例提出的多能量CT基材料分解方法，计算多能能谱和射线透射强度信息并进行优化，对射线透射强度信息和多能能谱降维消除散射，根据消除散射后的射线透射强度信息和多能能谱直接进行基材料分解，在多能数据不足的情况下，根据消除散射后的射线透射强度信息和多能能谱进行散射估计，根据散射估计对已有的透射数据进行去散射处理，进而得到去散射的投影数据，再进行基材料分解。由此，可以进行投影域或图像域的基材料物质分解、稀疏数据图像重建等工作。
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
Claims (10)
1.一种多能量CT基材料分解方法，其特征在于，包括以下步骤：
S1，获取多能谱成像系统的多能能谱信息和扫描待成像物体时探测器探测到的射线透射强度信息；
S2，对所述多能能谱进行降维处理，得到新的多能能谱；
S3，对所述射线透射强度信息进行降维处理，得到新的射线透射强度信息；
S4，根据所述新的多能能谱和所述新的射线透射强度信息对所述待成像物体进行基材料分解。
2.根据权利要求1所述的方法，其特征在于，在步骤S4之后还包括：
S5，根据所述新的多能能谱、所述新的射线透射强度信息以及消除散射方程计算基材料分解系数；
S6，根据计算出的所述基材料分解系数与所述多能能谱和所述射线透射强度信息之间的映射关系反解出散射值，根据所述散射值和散射值在邻域分布的相似性得到整体散射估计；
S7，根据所述整体散射估计对所述射线透射强度信息或优化后的射线透射强度信息进行去散射处理，得到去散射后的射线透射强度信息，根据所述去散射后的射线透射强度信息进行基材料分解。
3.根据权利要求1所述的方法，其特征在于，在对所述多能能谱和所述射线透射强度信息进行降维处理之前，包括：对所述多能能谱和所述射线透射强度信息进行插值优化处理，以使所述多能能谱和所述射线透射强度信息趋于一致。
4.根据权利要求1所述的方法，其特征在于，以所述多能能谱或所述射线透射强度信息作为自变量构建降维函数，通过所述降维函数消除一维的散射信息，将N维的所述多能能谱和所述射线透射强度信息降为不包含散射信息的N-1维信息。
5.根据权利要求4所述的方法，其特征在于，
所述降维函数为所述多能能谱两两相减和所述射线透射强度信息两两相减；或
优化后的多能能谱两两相减和优化后的射线透射强度信息两两相减。
6.根据权利要求2所述的方法，其特征在于，所述S5包括：
所述消除散射方程参数包括所述新的射线透射强度信息、所述新的多能能谱、多个基材料参数，根据所述新的射线透射强度信息、所述新的多能能谱和部分基材料参数求解所述基材料分解系数，其中，所述多个基材料参数包括所述基材料分解系数。
7.根据权利要求2所述的方法，其特征在于，所述映射关系为探测器上接收到的射线透射强度信息与所述多能能谱、基材料衰减系数、等效厚度、对应能谱、散射值的关系式。
8.根据权利要求2所述的方法，其特征在于，所述S5中，计算所述基材料分解系数的方法包括但不限于迭代求解、查找表法和多项式拟合法。
9.根据权利要求1所述的方法，其特征在于，通过N能能谱和N能透射数据对小于等于N-1种基材料分解，其中，N大于等于3。
10.根据权利要求2所述的方法，其特征在于，在所述S6中，所述根据所述散射值和散射值在邻域分布的相似性得到整体散射估计，包括：
利用包括但不限于插值或加权的方式估计出整体的散射分布，完成散射估计。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010054618.1A CN111134709B (zh) | 2020-01-17 | 2020-01-17 | 一种多能量ct基材料分解方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010054618.1A CN111134709B (zh) | 2020-01-17 | 2020-01-17 | 一种多能量ct基材料分解方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111134709A true CN111134709A (zh) | 2020-05-12 |
| CN111134709B CN111134709B (zh) | 2021-09-14 |
Family
ID=70525891
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010054618.1A Active CN111134709B (zh) | 2020-01-17 | 2020-01-17 | 一种多能量ct基材料分解方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111134709B (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114549360A (zh) * | 2022-02-25 | 2022-05-27 | 清华大学 | 多能量ct去散射降噪物质分解方法、装置、设备及介质 |
Citations (21)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101647706A (zh) * | 2008-08-13 | 2010-02-17 | 清华大学 | 高能双能ct系统的图象重建方法 |
| CN102095663A (zh) * | 2007-10-05 | 2011-06-15 | 清华大学 | 液态物品检查方法和设备 |
| CN102313897A (zh) * | 2010-06-29 | 2012-01-11 | 成都理工大学 | 一种放射性能谱识别方法 |
| CN102608135A (zh) * | 2007-10-05 | 2012-07-25 | 清华大学 | 在危险品检查系统中确定ct扫描位置的方法和设备 |
| CN103364419A (zh) * | 2012-04-04 | 2013-10-23 | 通用电气公司 | 用于多能计算机断层扫描成像的系统和方法 |
| CN105246410A (zh) * | 2013-05-30 | 2016-01-13 | 株式会社东芝 | 光子计数x射线计算机断层摄影装置、以及光子计数x射线计算机断层摄影方法 |
| CN106251383A (zh) * | 2016-08-31 | 2016-12-21 | 东北大学 | 一种能谱ct基物质正弦图的估计方法 |
| CN106473761A (zh) * | 2016-10-14 | 2017-03-08 | 山东大学 | 一种医学双能ct 电子密度图像的重建与数值校准方法 |
| CN106796733A (zh) * | 2014-10-13 | 2017-05-31 | 皇家飞利浦有限公司 | 用于谱成像的噪声降低 |
| CN107019518A (zh) * | 2016-02-01 | 2017-08-08 | 通用电气公司 | 用于计算机断层扫描中的散射校正的信号处理方法及成像系统 |
| US20180035965A1 (en) * | 2014-06-20 | 2018-02-08 | Marquette University | Material decomposition of multi-spectral x-ray projections using neural networks |
| CN108010098A (zh) * | 2017-12-04 | 2018-05-08 | 首都师范大学 | 一种双能谱ct基材料图像迭代重建方法 |
| CN108280860A (zh) * | 2017-01-03 | 2018-07-13 | 西门子医疗有限公司 | 减少多能量ct成像中的光谱通道数目 |
| WO2019006310A1 (en) * | 2017-06-29 | 2019-01-03 | Cuadros Angela | PIXELIZED K-ENCODED OPENING SYSTEM FOR X-RAY IMAGING WITH COMPRESSION SPECTRUM |
| CN109875593A (zh) * | 2019-02-02 | 2019-06-14 | 中北大学 | Ct成像方法、存储介质和装置 |
| CN109959669A (zh) * | 2017-12-26 | 2019-07-02 | 清华大学 | 基材料分解方法和装置及计算机可读存储介质 |
| CN110175957A (zh) * | 2019-04-11 | 2019-08-27 | 清华大学 | 多能量ct基材料物质分解方法 |
| CN110189389A (zh) * | 2019-05-31 | 2019-08-30 | 首都师范大学 | 基于深度学习的双能谱ct投影域基材料分解方法及装置 |
| CN110179486A (zh) * | 2019-04-11 | 2019-08-30 | 清华大学 | 多能量ct成像系统及其应用 |
| US20190313993A1 (en) * | 2018-04-12 | 2019-10-17 | Canon Medical Systems Corporation | Method and apparatus for computed tomography (ct) and material decomposition with pile-up correction calibrated using a real pulse pileup effect and detector response |
| CN110353716A (zh) * | 2019-07-25 | 2019-10-22 | 沈阳开普医疗影像技术有限公司 | 双能ct短程扫描协议中的数字域分解方法 |
-
2020
- 2020-01-17 CN CN202010054618.1A patent/CN111134709B/zh active Active
Patent Citations (21)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102095663A (zh) * | 2007-10-05 | 2011-06-15 | 清华大学 | 液态物品检查方法和设备 |
| CN102608135A (zh) * | 2007-10-05 | 2012-07-25 | 清华大学 | 在危险品检查系统中确定ct扫描位置的方法和设备 |
| CN101647706A (zh) * | 2008-08-13 | 2010-02-17 | 清华大学 | 高能双能ct系统的图象重建方法 |
| CN102313897A (zh) * | 2010-06-29 | 2012-01-11 | 成都理工大学 | 一种放射性能谱识别方法 |
| CN103364419A (zh) * | 2012-04-04 | 2013-10-23 | 通用电气公司 | 用于多能计算机断层扫描成像的系统和方法 |
| CN105246410A (zh) * | 2013-05-30 | 2016-01-13 | 株式会社东芝 | 光子计数x射线计算机断层摄影装置、以及光子计数x射线计算机断层摄影方法 |
| US20180035965A1 (en) * | 2014-06-20 | 2018-02-08 | Marquette University | Material decomposition of multi-spectral x-ray projections using neural networks |
| CN106796733A (zh) * | 2014-10-13 | 2017-05-31 | 皇家飞利浦有限公司 | 用于谱成像的噪声降低 |
| CN107019518A (zh) * | 2016-02-01 | 2017-08-08 | 通用电气公司 | 用于计算机断层扫描中的散射校正的信号处理方法及成像系统 |
| CN106251383A (zh) * | 2016-08-31 | 2016-12-21 | 东北大学 | 一种能谱ct基物质正弦图的估计方法 |
| CN106473761A (zh) * | 2016-10-14 | 2017-03-08 | 山东大学 | 一种医学双能ct 电子密度图像的重建与数值校准方法 |
| CN108280860A (zh) * | 2017-01-03 | 2018-07-13 | 西门子医疗有限公司 | 减少多能量ct成像中的光谱通道数目 |
| WO2019006310A1 (en) * | 2017-06-29 | 2019-01-03 | Cuadros Angela | PIXELIZED K-ENCODED OPENING SYSTEM FOR X-RAY IMAGING WITH COMPRESSION SPECTRUM |
| CN108010098A (zh) * | 2017-12-04 | 2018-05-08 | 首都师范大学 | 一种双能谱ct基材料图像迭代重建方法 |
| CN109959669A (zh) * | 2017-12-26 | 2019-07-02 | 清华大学 | 基材料分解方法和装置及计算机可读存储介质 |
| US20190313993A1 (en) * | 2018-04-12 | 2019-10-17 | Canon Medical Systems Corporation | Method and apparatus for computed tomography (ct) and material decomposition with pile-up correction calibrated using a real pulse pileup effect and detector response |
| CN109875593A (zh) * | 2019-02-02 | 2019-06-14 | 中北大学 | Ct成像方法、存储介质和装置 |
| CN110175957A (zh) * | 2019-04-11 | 2019-08-27 | 清华大学 | 多能量ct基材料物质分解方法 |
| CN110179486A (zh) * | 2019-04-11 | 2019-08-30 | 清华大学 | 多能量ct成像系统及其应用 |
| CN110189389A (zh) * | 2019-05-31 | 2019-08-30 | 首都师范大学 | 基于深度学习的双能谱ct投影域基材料分解方法及装置 |
| CN110353716A (zh) * | 2019-07-25 | 2019-10-22 | 沈阳开普医疗影像技术有限公司 | 双能ct短程扫描协议中的数字域分解方法 |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114549360A (zh) * | 2022-02-25 | 2022-05-27 | 清华大学 | 多能量ct去散射降噪物质分解方法、装置、设备及介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111134709B (zh) | 2021-09-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhu et al. | Scatter correction method for X-ray CT using primary modulation: Theory and preliminary results | |
| Sun et al. | Improved scatter correction using adaptive scatter kernel superposition | |
| JP7139119B2 (ja) | 医用画像生成装置及び医用画像生成方法 | |
| US7272429B2 (en) | Methods and apparatus for facilitating a reduction in artifacts | |
| US6639964B2 (en) | Device and method for forming a computed X-ray tomogram with scatter correction | |
| US6891918B2 (en) | Methods and apparatus for acquiring perfusion data | |
| US8055039B2 (en) | System and method to obtain noise mitigated monochromatic representation for varying energy level | |
| Gao et al. | Scatter correction method for x‐ray CT using primary modulation: phantom studies | |
| CN110175957B (zh) | 多能量ct基材料物质分解方法 | |
| CN1635850A (zh) | X射线散射校正 | |
| Byl et al. | Photon‐counting normalized metal artifact reduction (NMAR) in diagnostic CT | |
| JP2009201885A (ja) | Ｘ線ｃｔ装置 | |
| CN110636796B (zh) | X射线暗场成像中的射束硬化校正 | |
| Niu et al. | Overview of x-ray scatter in cone-beam computed tomography and its correction methods | |
| EP3230953A1 (en) | Methods and systems for spectral ct imaging | |
| US9943279B2 (en) | Methods and systems for task-based data generation and weighting for CT spectral imaging | |
| JP7370989B2 (ja) | スペクトルボリューム画像データを生成するように構成された非スペクトルコンピュータ断層撮影（ｃｔ）スキャナ | |
| CN110084864A (zh) | 一种基于能谱ct的电子密度图像重建方法 | |
| Li et al. | Implementation of dual‐energy technique for virtual monochromatic and linearly mixed CBCTs | |
| CN110179486A (zh) | 多能量ct成像系统及其应用 | |
| Kyriakou et al. | Intensity distribution and impact of scatter for dual-source CT | |
| CN111134709B (zh) | 一种多能量ct基材料分解方法 | |
| Tang et al. | Investigation into the influence of x-ray scatter on the imaging performance of an x-ray flat-panel imager-based cone-beam volume CT | |
| US10383589B2 (en) | Direct monochromatic image generation for spectral computed tomography | |
| CN114549360A (zh) | 多能量ct去散射降噪物质分解方法、装置、设备及介质 |
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 |