








本申请基于申请号为201510163516.2、申请日为2015年4月8日的中国专利申请提出,并要求该中国专利申请的优先权,该中国专利申请的全部内容在此引入本申请作为参考。The present application is based on a Chinese patent application filed on Apr. 8, 2015, the entire disclosure of which is hereby incorporated by reference.
本公开涉及CT成像技术,具体而言涉及能谱CT成像系统及其数据采集方法,以及重建能谱CT图像的方法。The present disclosure relates to CT imaging techniques, and in particular to an energy spectrum CT imaging system and data acquisition method thereof, and a method of reconstructing an energy spectrum CT image.
X射线CT成像技术在医疗、安检、工业无损检测等领域有着非常广泛的应用。能谱CT是近年来受到普遍关注和迅速发展的研究方向。目前基于光子计数探测器的能谱CT系统是能谱CT的主要实现方式之一。国际上有多个研究单位和公司致力于研究和生产适用于X光CT的光子计数探测器。X-ray CT imaging technology has a wide range of applications in medical, security, and industrial non-destructive testing. Energy spectrum CT is a research direction that has received widespread attention and rapid development in recent years. At present, the energy spectrum CT system based on photon counting detector is one of the main implementation methods of energy spectrum CT. There are several research units and companies in the world dedicated to the research and production of photon counting detectors for X-ray CT.
光子计数探测器对接收到的光子进行处理,按照能量阈值把接收到的光子归到相应的计数器,从而使得可以选择能量窗采集多个能量下的X射线穿过物体的信号数据。通过对应的能谱CT图像重建算法可以得到物体的三维物理信息,其中包括线性衰减系数、电子密度、等效原子序数、材料分布等等。因此,能谱CT相在功能上能够比传统CT提供被成像体更多层面的信息。The photon counting detector processes the received photons, and the received photons are returned to the corresponding counters according to the energy threshold, so that the energy window can be selected to collect the signal data of the X-rays passing through the object under multiple energies. The three-dimensional physical information of the object can be obtained by the corresponding energy spectrum CT image reconstruction algorithm, including linear attenuation coefficient, electron density, equivalent atomic number, material distribution and the like. Therefore, the spectral CT phase is functionally capable of providing more levels of information than the conventional CT.
另一方面,由于计数方式可以通过阈值方式消除电子学噪声,进一步提高数据的可靠性和信噪比。因此,能谱CT的图像可以具有更高的图像质量,或者在同等图像质量情况下降低剂量。On the other hand, since the counting method can eliminate the electronic noise by the threshold method, the reliability and the signal-to-noise ratio of the data are further improved. Therefore, an image of the energy spectrum CT can have a higher image quality or a lower dose in the case of equivalent image quality.
光子计数探测器的性能是影响能谱CT成像效果的关键因素之一。目前应用于X光CT的像素化光子计数面阵探测器的像素分辨率主要在50-100微米之间,光子计数率在106-108/mm2之间。一方面,为了获得高质量的重建图像,需要采集足够多的光子以降低量子噪声的影响;另一方面,每个像素计数器的计数速度受电子学的硬件信号处理速度限制。由此,虽然减小像素大小有助于获得更高的每平方毫米计数率,但是减小像素使得像素之间的电荷共享问题严重化。而且,像素数目的增加对于探测器的读出速率要求也显著提高。The performance of photon counting detectors is one of the key factors affecting the imaging performance of energy spectrum CT. The pixel resolution of the pixelated photon counting area array detector currently applied to X-ray CT is mainly between 50-100 microns, and the photon counting rate is between 106 -108 /mm2 . On the one hand, in order to obtain high quality reconstructed images, it is necessary to collect enough photons to reduce the influence of quantum noise; on the other hand, the counting speed of each pixel counter is limited by the electronic hardware signal processing speed. Thus, while reducing the pixel size helps to achieve a higher count rate per square millimeter, reducing the pixel complicates the charge sharing problem between pixels. Moreover, the increase in the number of pixels also significantly increases the read rate requirements of the detector.
因此,需要对能谱CT进行改进或开发信的能谱CT系统。Therefore, there is a need for an energy spectrum CT system that improves or develops a spectral CT.
在所述背景技术部分公开的上述信息仅用于加强对本公开的背景的理解,因此它可以包括不构成对本领域普通技术人员已知的现有技术的信息。The above information disclosed in this Background section is only for enhancement of understanding of the background of the present disclosure, and thus it may include information that does not constitute a prior art known to those of ordinary skill in the art.
发明内容Summary of the invention
本申请公开一种能谱CT成像系统及其数据采集方法,以及重建能谱CT图像的方法,能够克服上述问题中的一个或多个。The present application discloses an energy spectrum CT imaging system and a data acquisition method thereof, and a method for reconstructing an energy spectrum CT image.One or more of the above problems can be overcome.
本公开的其他特性和优点将通过下面的详细描述变得显然,或部分地通过本公开的实践而习得。Other features and advantages of the present disclosure will be apparent from the following detailed description.
根据本公开的一个方面,提供一种能谱CT成像系统的数据采集方法,所述能谱CT成像系统包括具有多个像素的光子探测器,所述数据采集方法包括在多个投影角度从所述光子探测器对穿过被成像物的光子进行数据采集,对于每一投影角度,该数据采集方法包括:按照多个随机分布的像素地址编码从所述光子探测器获取数据。According to an aspect of the present disclosure, a data acquisition method of a spectral CT imaging system is provided, the energy spectrum CT imaging system including a photon detector having a plurality of pixels, the data acquisition method including a plurality of projection angles The photon detector performs data acquisition on photons passing through the imaged object. For each projection angle, the data acquisition method includes: acquiring data from the photon detector according to a plurality of randomly distributed pixel address codes.
根据本公开的一些实施例,所述像素地址编码利用设定的像素地址分布的概率密度函数而生成。According to some embodiments of the present disclosure, the pixel address encoding is generated using a probability density function of a set pixel address distribution.
根据本公开的一些实施例,对于第一个投影角度,所述像素地址编码利用设定的像素地址分布的概率密度函数而生成;对于除第一个投影角度外的其余投影角度,所述像素地址编码或者利用为前一投影角度设定的像素地址分布的概率密度函数而生成、或者利用重新设定的像素地址分布的概率密度函数而生成、或者与用于前一投影角度的所述像素地址编码相同。According to some embodiments of the present disclosure, for a first projection angle, the pixel address encoding is generated using a probability density function of the set pixel address distribution; for the remaining projection angles other than the first projection angle, the pixels The address code is generated either by a probability density function of the pixel address distribution set for the previous projection angle, or by a probability density function of the reset pixel address distribution, or with the pixel for the previous projection angle The address code is the same.
根据本公开的一些实施例,所述按照多个随机分布的像素地址编码从所述光子探测器采集数据包括:从自所述光子探测器获得的全部数据中按照多个随机分布的像素地址编码选取数据。According to some embodiments of the present disclosure, the collecting data from the photon detector according to a plurality of randomly distributed pixel address encodings comprises: encoding, according to a plurality of randomly distributed pixel addresses, from all data obtained from the photon detector Select the data.
根据本公开的一些实施例,所述按照多个随机分布的像素地址编码从所述光子探测器采集数据包括:从所述光子探测器的全部像素中按照多个随机分布的像素地址编码实时选择像素来获取数据。According to some embodiments of the present disclosure, the acquiring data from the photon detector according to a plurality of randomly distributed pixel address encodings comprises: real-time selection from a plurality of randomly distributed pixel address encodings from all pixels of the photon detector Pixels to get the data.
根据本公开的另一发明,提供一种能谱CT成像系统的数据采集方法,所述能谱CT成像系统包括具有多个像素的光子探测器,所述数据采集方法包括在多个投影角度从所述光子探测器对穿过被成像物的光子进行数据采集,对于每一投影角度,该数据采集方法包括:从按照多个随机分布的像素地址编码设置的所述多个像素采集数据。根据本公开的一些实施例,所述光子探测器为面阵探测器或环形探测器。According to another invention of the present disclosure, there is provided a data acquisition method for a spectral CT imaging system, the energy spectrum CT imaging system comprising a photon detector having a plurality of pixels, the data acquisition method comprising The photon detector performs data acquisition on photons passing through the imaged object. For each projection angle, the data acquisition method includes: collecting data from the plurality of pixels set according to a plurality of randomly distributed pixel address encodings. According to some embodiments of the present disclosure, the photon detector is an area array detector or a ring detector.
根据本公开的另一方面,提供一种重建能谱CT图像的方法,包括:设置数据采集的多个能量窗;利用如前所述的任一数据采集方法从所述多个能量窗进行数据采集;利用所述像素地址编码生成对应的系统矩阵;将采集到的数据在投影域进行能谱信息分解,得到分解系数投影;通过迭代方法,利用所述分解系数投影进行分解系数的空间分布重建,得到分解系数;根据重建的分解系数进行合成,得到单能量的衰减系数和/或估算被成像物的电子密度和/或等效原子序数分布图。In accordance with another aspect of the present disclosure, a method of reconstructing an energy spectrum CT image is provided, comprising: setting a plurality of energy windows for data acquisition; and performing data from the plurality of energy windows using any of the data acquisition methods described above Collecting; generating a corresponding system matrix by using the pixel address encoding; decomposing the collected data in the projection domain to obtain a decomposition coefficient projection; and using the iterative method, using the decomposition coefficient projection to perform spatial distribution reconstruction of the decomposition coefficient Obtaining a decomposition coefficient; synthesizing according to the reconstructed decomposition coefficient, obtaining a single energy attenuation coefficient and/or estimating an electron density and/or an equivalent atomic number distribution map of the imaged object.
根据本公开的一些实施例,所述系统矩阵的维度是M×N,其中N是被重建的图像域像素数目,M是投影角度数与所述像素地址编码的数量之积。According to some embodiments of the present disclosure, the dimension of the system matrix is M×N, where N is the number of reconstructed image domain pixels, and M is the product of the number of projection angles and the number of pixel address codes.
根据本公开的一些实施例,根据如下公式求分解系数投影:According to some embodiments of the present disclosure, the decomposition coefficient projection is obtained according to the following formula:
其中τ为整数,对应第τ个分解项,wk(E)为归一化的第k个能量窗下的X光能谱分布,pk为对应于光子探测器上采集的数据的M维向量,φτ(E)为第τ个分解基函数,Aτ为第τ个分解系数投影,k=1,2,...K,K为能量窗数目,E是光子能量。Where τ is an integer corresponding to the τth decomposition term, wk (E) is the X-ray energy spectrum distribution under the normalized k-th energy window, and pk is the M-dimensional corresponding to the data collected on the photon detector Vector, φτ (E) is the τth decomposition basis function, Aτ is the τth decomposition coefficient projection, k=1, 2, ... K, K is the number of energy windows, and E is the photon energy.
根据本公开的一些实施例,根据如下公式通过迭代进行分解系数的空间分布重建:According to some embodiments of the present disclosure, the spatially distributed reconstruction of the decomposition coefficients is iteratively performed according to the following formula:
其中ε为根据数据噪声确定的阈值,Ф(i,j,z)为关于成像视野的像素位置(i,j,z)的先验函数,aτ为分解系数,H是所述系统矩阵。Where ε is the threshold determined from the data noise, Ф(i, j, z) is the a priori function with respect to the pixel position (i, j, z) of the imaging field of view, aτ is the decomposition coefficient, and H is the system matrix.
根据本公开的一些实施例,根据重建结果通过下式计算单能量的衰减系数:According to some embodiments of the present disclosure, the attenuation coefficient of the single energy is calculated by the following formula according to the reconstruction result:
根据本公开的一些实施例,在两个能量窗的情况下,K=2,τ∈[1,2],φτ(E)分别取为光电效应系数、康普顿散射系数。According to some embodiments of the present disclosure, in the case of two energy windows, K=2, τ∈[1, 2], φτ(E) are taken as the photoelectric effect coefficient and the Compton scattering coefficient, respectively.
根据本公开的一些实施例,等效原子序数分布图、电子密度分布图分别通过下式估算:According to some embodiments of the present disclosure, the equivalent atomic number distribution map and the electron density distribution map are respectively estimated by the following formula:
ρe=2a2ρe =2a2
其中a1、a2对应于τ分别为1和2时的aτ值。Where a1 and a2 correspond to aτ values when τ is 1 and 2, respectively.
根据本公开的另一方面,提供一种重建能谱CT图像的方法,包括:设置数据采集的多个能量窗;利用如前所述的任一数据采集方法从所述多个能量窗进行数据采集;利用所述像素地址编码生成对应的系统矩阵;通过迭代方法,对每个能量窗下采集到的数据进行衰减系数的空间分布重建,得到衰减系数;利用衰减系数的重建结果,对被重建的图像域像素逐个进行能谱信息分解,得到分解系数;估算被成像物的电子密度和/或等效原子序数分布图。In accordance with another aspect of the present disclosure, a method of reconstructing an energy spectrum CT image is provided, comprising: setting a plurality of energy windows for data acquisition; and performing data from the plurality of energy windows using any of the data acquisition methods described above Collecting; using the pixel address coding to generate a corresponding system matrix; by using an iterative method, spatially reconstructing the attenuation coefficient of each data collected under each energy window to obtain an attenuation coefficient; using the reconstruction result of the attenuation coefficient, the reconstruction is performed The image domain pixels are decomposed one by one to obtain the decomposition coefficient; the electron density of the imaged object and/or the equivalent atomic number distribution map are estimated.
根据本公开的一些实施例,所述系统矩阵的维度是M×N,其中N是被重建的图像域像素数目,M是投影角度数与所述像素地址编码的数量之积。According to some embodiments of the present disclosure, the dimension of the system matrix is M×N, where N is the number of reconstructed image domain pixels, and M is the product of the number of projection angles and the number of pixel address codes.
根据本公开的一些实施例,根据如下公式通过迭代进行衰减系数重建:According to some embodiments of the present disclosure, the attenuation coefficient reconstruction is performed by iteration according to the following formula:
其中k表示第k个能量窗,gk为对采集到的数据进行归一化和负对数处理得到的投影值,μk为衰减系数,pk为对应于光子探测器上采集的数据的M维向量,ε为根据数据噪声确定的阈值,Ф(i,j,z)为关于成像视野的像素位置(i,j,z)的先验函数,H是所述系统矩阵。Where k is the kth energy window, gk is the projection value obtained by normalizing and negative logarithmically processing the acquired data, μk is the attenuation coefficient, and pk is the data corresponding to the data collected on the photon detector The M-dimensional vector, ε is a threshold determined according to data noise, Ф(i, j, z) is a priori function with respect to the pixel position (i, j, z) of the imaging field of view, and H is the system matrix.
根据本公开的一些实施例,根据重建结果μk,通过求解如下线性方程组对被重建的图像域像素逐个进行能谱信息分解而得到分解系数:According to some embodiments of the present disclosure, according to the reconstruction result μk , the decomposition coefficients are obtained by performing spectral information decomposition on the reconstructed image domain pixels one by one by solving the following linear equations:
其中τ为整数,对应第τ个分解项,k=1,...,K,K为能量窗数目,φτ(E)为第τ个分解基函数,aτ(i,j,z)为第τ个分解系数,Ek表示第k个能量窗的等效光子能量。Where τ is an integer corresponding to the τth decomposition term, k=1,...,K,K is the number of energy windows, φτ (E) is the τth decomposition basis function, aτ (i,j,z) For the τth decomposition coefficient, Ek represents the equivalent photon energy of the kth energy window.
根据本公开的一些实施例,在两个能量窗的情况下,K=2,τ∈[1,2],φτ(E)分别取为光电效应系数、康普顿散射系数。According to some embodiments of the present disclosure, in the case of two energy windows, K = 2, τ ∈ [1, 2], φτ (E) are taken as the photoelectric effect coefficient and the Compton scattering coefficient, respectively.
根据本公开的一些实施例,等效原子序数分布图、电子密度分布图分别通过下式估算:According to some embodiments of the present disclosure, the equivalent atomic number distribution map and the electron density distribution map are respectively estimated by the following formula:
ρe=2a2ρe =2a2
其中a1、a2对应于τ分别为1和2时的aτ值。Where a1 and a2 correspond to aτ values when τ is 1 and 2, respectively.
根据本公开的另一方面,提供一种能谱CT成像系统,包括:射线发生装置,包括射线源;光子探测器,包括多个像素;及数据采集系统,从所述光子探测器对穿过被成像物的光子进行数据采集,其特征在于,所述数据采集系统按照多个随机分布的像素地址编码从所述光子探测器获取数据。According to another aspect of the present disclosure, there is provided an energy spectrum CT imaging system comprising: a radiation generating device comprising a radiation source; a photon detector comprising a plurality of pixels; and a data acquisition system passing through the photon detector pair The photon of the imaged object is subjected to data acquisition, wherein the data acquisition system acquires data from the photon detector according to a plurality of randomly distributed pixel address codes.
根据本公开的一些实施例,所述像素地址编码利用设定的像素地址分布的概率密度函数而生成。According to some embodiments of the present disclosure, the pixel address encoding is generated using a probability density function of a set pixel address distribution.
根据本公开的一些实施例,所述数据采集系统从所述光子探测器输出的全部数据中按照多个随机分布的像素地址编码选取数据。According to some embodiments of the present disclosure, the data acquisition system encodes selected data from a plurality of randomly distributed pixel address codes from all data output by the photon detector.
根据本公开的一些实施例,所述数据采集系统包括电子学系统,所述电子学系统配置为从所述光子探测器的全部像素中按照多个随机分布的像素地址编码实时选择像素来读取数据。According to some embodiments of the present disclosure, the data acquisition system includes an electronics system configured to encode real-time selected pixels from a plurality of randomly distributed pixel address encodings from all pixels of the photon detector data.
根据本公开的一些实施例,能谱CT成像系统还包括数据处理系统,用于根据前述的任一项重建能谱CT图像的方法利用所述数据采集系统获取的数据重建能谱CT图像。根据本公开的一些实施例,所述光子探测器为面阵探测器或环形探测器。According to some embodiments of the present disclosure, the energy spectrum CT imaging system further includes a data processing system for reconstructing the energy spectrum CT image using the data acquired by the data acquisition system in accordance with any of the foregoing methods of reconstructing an energy spectrum CT image. According to some embodiments of the present disclosure, the photon detector is an area array detector or a ring detector.
根据本公开的另一方面,提供一种能谱CT成像系统,包括:射线发生装置,包括射线源;光子探测器,包括多个像素;及数据采集系统,从所述光子探测器对穿过被成像物的光子进行数据采集。按照多个随机分布的像素地址编码设置所述光子探测器的所述多个像素。According to another aspect of the present disclosure, there is provided an energy spectrum CT imaging system comprising: a radiation generating device comprising a radiation source; a photon detector comprising a plurality of pixels; and a data acquisition system passing through the photon detector pair Imaged objectPhotons are used for data acquisition. The plurality of pixels of the photon detector are set according to a plurality of randomly distributed pixel address encodings.
根据本公开的一些实施例,所述像素地址编码利用设定的像素地址分布的概率密度函数而生成。According to some embodiments of the present disclosure, the pixel address encoding is generated using a probability density function of a set pixel address distribution.
根据本公开的一些实施例,能谱CT成像系统还包括数据处理系统,用于根据前述任一重建能谱CT图像的方法利用所述数据采集系统获取的数据重建能谱CT图像。根据本公开的一些实施例,所述光子探测器为面阵探测器或环形探测器。According to some embodiments of the present disclosure, the energy spectrum CT imaging system further includes a data processing system for reconstructing the energy spectrum CT image using the data acquired by the data acquisition system according to any of the methods of reconstructing the energy spectrum CT image. According to some embodiments of the present disclosure, the photon detector is an area array detector or a ring detector.
根据本公开的一些实施例的能谱CT成像系统及其数据采集方法以及重建能谱CT图像的方法,能够有效控制电荷共享效应,并提高CT系统处理速度,以及降低系统成本。The energy spectrum CT imaging system and the data acquisition method thereof and the method of reconstructing the energy spectrum CT image according to some embodiments of the present disclosure can effectively control the charge sharing effect, improve the processing speed of the CT system, and reduce the system cost.
通过参照附图详细描述其示例实施例,本公开的上述和其它特征及优点将变得更加明显。The above and other features and advantages of the present disclosure will become more apparent from the detailed description of exemplary embodiments.
图1示意性示出根据本公开一些示例实施例的能谱CT成像系统的系统结构示意图;FIG. 1 is a schematic block diagram showing a system configuration of an energy spectrum CT imaging system according to some example embodiments of the present disclosure; FIG.
图2示意性示出根据本公开一些示例实施例的数据采集流程图;FIG. 2 schematically illustrates a data collection flow diagram in accordance with some example embodiments of the present disclosure; FIG.
图3a和图3b示意性示出根据本公开一些示例实施例的随机分布或随机选择的探测器像素,其中图3a为随机均匀分布模式,图3b为按照特定概率密度例如中间概率大、边缘概率小的随机分布模式;3a and 3b schematically illustrate randomly distributed or randomly selected detector pixels, wherein FIG. 3a is a random uniform distribution pattern and FIG. 3b is a specific probability density, eg, an intermediate probability, edge probability, according to some example embodiments of the present disclosure. Small random distribution pattern;
图4a和图4b示意性示出根据本公开一些示例实施例的投影数据比较,其中图4a示出普通能谱CT采集到的一个能量窗的投影数据,图4b示出根据随机分布的像素地址编码采集前者十分之一数据得到的一个能量窗下的投影数据;4a and 4b schematically illustrate projection data comparisons according to some example embodiments of the present disclosure, wherein FIG. 4a shows projection data of one energy window acquired by ordinary energy spectrum CT, and FIG. 4b shows pixel addresses according to random distribution. Coding the projection data under an energy window obtained by collecting one tenth of the data of the former;
图5a和图5b示出根据本公开一些示例实施例进行能谱CT重建的结果示例,其中图5a是电子密度分布图,图5b是等效原子序数分布图;及5a and 5b illustrate examples of results of performing spectral CT reconstruction according to some example embodiments of the present disclosure, wherein FIG. 5a is an electron density distribution map, and FIG. 5b is an equivalent atomic number distribution map;
图6示出根据本公开一实施例的能谱CT结构示意图,其包括具有预定义的随机分布像素的探测器。6 shows a schematic diagram of an energy spectrum CT structure including a detector having predefined randomly distributed pixels, in accordance with an embodiment of the present disclosure.
现在将参考附图更全面地描述示例实施例。然而,示例实施例能够以多种形式实施,且不应被理解为限于在此阐述的实施例;相反,提供这些实施例使得本公开将全面和完整,并将示例实施例的构思全面地传达给本领域的技术人员。在图中相同的附图标记表示相同或类似的部分,因而将省略对它们的重复描述。Example embodiments will now be described more fully with reference to the accompanying drawings. However, the exemplary embodiments can be embodied in a variety of forms and should not be construed as being limited to the embodiments set forth herein. To those skilled in the art. The same reference numerals in the drawings denote the same or similar parts, and the repeated description thereof will be omitted.
此外,所描述的特征、结构或特性可以以任何合适的方式结合在一个或更多实施例中。在下面的描述中,提供许多具体细节从而给出对本公开的实施例的充分理解。然而,本领域技术人员将意识到,可以实践本公开的技术方案而没有所述特定细节中的一个或更多,或者可以采用其它的方法、组元、材料、装置、步骤等。在其它情况下,不详细示出或描述公知结构、方法、装置、实现、材料或者操作以避免模糊本公开的各方面。Furthermore, the described features, structures, or characteristics may be combined in any suitable manner in one or more embodiments. In the following description, numerous specific details are set forth However, one skilled in the art will appreciate that the technical solution of the present disclosure may be practiced without one or more of the specific details, or other methods, components, materials, devices, steps, etc. may be employed. In other cases, it is not shown or depicted in detail.Well-known structures, methods, devices, implementations, materials, or operations are omitted to avoid obscuring aspects of the present disclosure.
附图中所示的方框图仅仅是功能实体,不一定必须与物理上独立的实体相对应。即,可以采用软件形式来实现这些功能实体,或在一个或多个硬件模块或集成电路中实现这些功能实体,或在不同网络和/或处理器装置和/或微控制器装置中实现这些功能实体。The block diagrams shown in the figures are merely functional entities and do not necessarily have to correspond to physically separate entities. That is, these functional entities may be implemented in software, or implemented in one or more hardware modules or integrated circuits, or implemented in different networks and/or processor devices and/or microcontroller devices. entity.
附图中所示的流程图仅是示例性说明,不是必须包括所有的步骤。例如,有的步骤还可以分解,而有的步骤可以合并或部分合并,因此实际执行的顺序有可能根据实际情况改变。The flowcharts shown in the figures are merely illustrative and not necessarily all of the steps. For example, some steps may be decomposed, and some steps may be combined or partially merged, so the actual execution order may vary depending on the actual situation.
图1示意性示出根据本公开一些示例实施例的能谱CT成像系统100的系统结构示意图。FIG. 1 schematically illustrates a system structure diagram of a spectral
如图1所示,能谱CT成像系统100可包括射线发生装置105、机械运动系统110、光子探测器120、以及数据采集系统130。根据本公开的系统可以通过圆轨道扫描实现,也可以通过螺旋轨迹扫描实现,可用于三维能谱CT成像。As shown in FIG. 1, the energy spectrum
射线发生装置105可包括射线源,用于发射例如X射线。The
机械运动系统110用于使被成像物与所述射线源发生相对运动。机械运动系统110可例如包括机械运动装置和对应的控制系统(未示出)。可以是被成像物体115运动而射线源和/或探测器120保持静止(图1所示方式),也可以是射线源和/或探测器运动,而物体保持静止。一般在医疗领域中避免转动病人,可通过转动源和/或探测器实现。在工业无损检测中,转动和平移物体的方式比较常见。对于CT成像,起作用的是相对运动,所以两种方式等效。The
光子探测器120可包括多个像素1202(如图3a所示),用于接收并处理穿过物体的光子。
数据采集系统130从所述光子探测器120对穿过被成像物115的光子进行数据采集。根据本公开的发明构思,数据采集系统130按照多个随机分布的像素地址编码从光子探测器120获取数据,如下面将要详细描述的。The
根据本公开的一些实施例,光子探测器120包括电子学系统125。数据采集系统130从电子学系统125输出的全部数据中按照多个随机分布的像素地址编码选取数据。According to some embodiments of the present disclosure,
根据本公开的另一些实施例,数据采集系统130包括电子学系统125。电子学系统125配置为从光子探测器120的全部像素中按照多个随机分布的像素地址编码实时选择像素来读取数据。According to further embodiments of the present disclosure,
根据本公开的另一些实施例,如图3a、图3b和图6所示,可按照多个随机分布的像素地址编码设置光子探测器120的多个像素1205。另外,如图3a、图3b和图6所示,光子探测器120可为例如面阵探测器或环形探测器。According to further embodiments of the present disclosure, as shown in Figures 3a, 3b, and 6, a plurality of
能谱CT成像系统100还可包括主控制器140及数据处理装置135。主控制器140负责能谱CT系统运行过程的主控制,包括机械转动、电气控制、安全连锁控制等。数据处理装置135对由数据采集系统130获得的数据进行处理,获得能谱CT任意能量下的衰减系数图像。另外,也可由此计算等效原子序数和电子密度分布图。这些图可通过断层或者三维可视化方式在显示器上显示。主控制器140及数据处理装置135可以是单个PC,也可以是工作站或计算机集群。The energy spectrum
图2示意性示出根据本公开一些示例实施例的数据采集流程图。下面参照图2、图3a和图3b描述根据本公开的从光子探测器120获取数据的方法。FIG. 2 schematically illustrates a data collection flow diagram in accordance with some example embodiments of the present disclosure. A method of acquiring data from
根据本公开的发明构思的数据采集方法,对于每一投影角度,该方法包括:按照多个随机分布的像素地址编码从所述光子探测器获取数据。According to a data acquisition method of the inventive concept of the present disclosure, for each projection angle, the method includes acquiring data from the photon detector according to a plurality of randomly distributed pixel address codes.
根据本公开的一些实施例,光子探测器120包括电子学系统125,该方法包括从自所述光子探测器获得的全部数据中按照多个随机分布的像素地址编码选取数据。In accordance with some embodiments of the present disclosure,
根据本公开的另一些实施例,数据采集系统130包括电子学系统125。电子学系统125配置为从光子探测器120的全部像素中按照多个随机分布的像素地址编码实时选择像素来读取数据。该方法包括从所述光子探测器的全部像素中按照多个随机分布的像素地址编码实时选择像素来获取数据。According to further embodiments of the present disclosure,
参照图2,在各个投影角度,可以按照一定的概率产生Ns个随机分布的像素地址编码,按照这些地址编码从探测器读取数值。这Ns个探测器像素的地址可以通过参数设置来控制其在整个探测器上的分布模式,即设置像素地址的概率密度函数。Referring to FIG. 2, at each projection angle, Ns randomly distributed pixel address codes can be generated with a certain probability, and values are read from the detector according to these address codes. The addresses of the Ns detector pixels can be parameterized to control their distribution pattern over the entire detector, ie, the probability density function of setting the pixel address.
随机分布的被读取像素的具体示例如图3a和3b所示。虚线网格表示实际的探测器单元矩阵,深灰色小方格表示在某一角度数据采集过程中选择的数据输出像素。图3a为按照均匀分布随机采样。图3b为按照特定概率密度,例如中间概率大、边缘概率小的方式进行的随机采样。Specific examples of randomly distributed read pixels are shown in Figures 3a and 3b. The dashed grid represents the actual detector cell matrix, and the dark gray small squares represent the data output pixels selected during a certain angle of data acquisition. Figure 3a shows random sampling in a uniform distribution. Figure 3b shows random sampling in a manner that has a specific probability density, such as a large intermediate probability and a small edge probability.
记探测器的总像素单元数目为Ndet,数据采集过程中每个视角读取的数据数目为Ns≤Ndet。根据本公开的技术方案,Ns可远小于Ndet,例如小于Ndet的50%,或更进一步地,小于Ndet的20%。The total number of pixel units of the detector is Ndet , and the number of data read by each view in the data acquisition process is Ns ≤ Ndet . According to the technical solution of the present disclosure, Ns may be much smaller than Ndet , for example, less than 50% of Ndet or, more completely, less than 20% of Ndet .
如图2所示,根据本公开一些实施例的方法包括如下处理。As shown in FIG. 2, a method in accordance with some embodiments of the present disclosure includes the following processing.
数据采集开始205之后,判断是否重新设置地址分布概率210。After the
如果为是,则设置地址分布的概率密度函数215,并进行生成地址编码225,然后按照地址编码读取探测器单元数据230。If so, the
如果为否,则判断是否更新地址编码220。如果为是,则转到225,否则转到230。例如,对于除第一个投影角度外的其余投影角度,像素地址编码可以利用为前一投影角度设定的像素地址分布的概率密度函数而生成(转到225);或者,与用于前一投影角度的所述像素地址编码相同(转到230)。If not, it is judged whether or not the
在按照地址编码读取探测器单元数据230之后,判断是否完成所有角度扫描。如果为否,则对下一个角度判断是否重新设置地址分布概率210。如果为是,则完成数据采集240。After the
下面说明根据本公开的利用通过前述方法采集的数据进行能谱CT图像重建的方法。A method of performing spectral CT image reconstruction using data acquired by the foregoing method according to the present disclosure will be described below.
设采集数据的多个能量窗(或谱)wk(E),k=1,2,...K。K为能量窗数目,K=2时为双能成像系统。各个wk(E)覆盖的能量范围可以有重叠部分也可以完全分开。多能CT的数据分别为用pk表示:Let multiple energy windows (or spectra) wk (E), k = 1, 2, ... K of the acquired data. K is the number of energy windows, and K=2 is a dual-energy imaging system. The energy range covered by each wk (E) may have overlapping portions or may be completely separated. The data of the multi-energy CT is denoted by pk :
-ln∫wk(E)exp(-Hμ(E))dE=pk (1)-ln∫wk (E)exp(-Hμ(E))dE=pk (1)
其中,E是光子能量。wk(E)为归一化的第k个能量窗下的X光能谱分布,可以是多种方式产生,例如设置光子计数探测器的能窗阈值,或者在光源处使用不同的滤光片等。μ(E)是物体的衰减系数,H是M×N维投影矩阵,pk为M维向量,即从光子计数探测器的随机分布的像素上采集M条射线的数据。在本公开中,这样的光子计数探测器可称为随机地址编码的光子计数探测器Coded-PCD。Where E is the photon energy. wk (E) is the X-ray energy spectrum distribution under the normalized kth energy window, which can be generated in various ways, such as setting the energy window threshold of the photon counting detector, or using different filtering at the light source. Film and so on. μ(E) is the attenuation coefficient of the object, H is the M×N-dimensional projection matrix, and pk is the M-dimensional vector, that is, the data of the M rays is collected from the randomly distributed pixels of the photon counting detector. In the present disclosure, such a photon counting detector may be referred to as a random address encoded photon counting detector Coded-PCD.
能谱CT的重建从框架上可以有三种方式(参见Y.Xing,et al.,″A Reconstruction Method for Dual High-Energy CT With MeV X-Rays,″Nuclear Science,IEEE Transactions on,vol.58,pp.537-546,2011):(1)预处理方式:首先把采集到的数据在投影域进行能谱信息分解,然后通过迭代进行空间信息重建,再根据重建的空间信息进行合成;(2)后处理方式:首先对于每一个能谱下的数据进行空间信息重建,然后在图像域进行能谱信息分解;(3)综合处理方式:建立综合的数据模型,通过迭代重建,同时获得能谱相关和空间位置相关信息。The reconstruction of the energy spectrum CT can be done in three ways from the framework (see Y.Xing, et al., "A Reconstruction Method for Dual High-Energy CT With MeV X-Rays," Nuclear Science, IEEE Transactions on, vol. Pp.537-546, 2011): (1) Preprocessing method: Firstly, the collected data is decomposed in the projection domain, and then the spatial information is reconstructed by iteration, and then synthesized according to the reconstructed spatial information; (2) Post-processing method: firstly, spatial information reconstruction is performed for each data under the energy spectrum, and then spectral information decomposition is performed in the image domain; (3) Comprehensive processing method: establishing a comprehensive data model, and obtaining energy spectrum through iterative reconstruction Relevant and spatial location related information.
根据本公开的能谱CT重建可以在预处理和后处理方法框架下完成。其中的能谱信息分解方式以及合成可以使用领域内现有的技术完成(参见G.Zhang,et al.,″A practical reconstruction method for dual energy computed tomography,″Journal of X-ray Science and Technology,vol.16,pp.67-88 2008以及Y.Xing,et al.,″A General Adaptive Decomposition method for Multi-Energy Spectral CT,″in Nuclear Science Symposium and Medical Imaging Conference(NSS/MIC),2011 IEEE,Seul,Korea,2013,pp.M12-15)。Energy spectrum CT reconstruction in accordance with the present disclosure can be accomplished within the framework of pre- and post-processing methods. The decomposition of the energy spectrum information and the synthesis can be accomplished using existing techniques in the field (see G. Zhang, et al., "A practical reconstruction method for dual energy computed tomography," Journal of X-ray Science and Technology, vol .16, pp. 67-88 2008 and Y. Xing, et al., "A General Adaptive Decomposition method for Multi-Energy Spectral CT," in Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC), 2011 IEEE, Seul , Korea, 2013, pp. M12-15).
空间信息重建方法可结合光子探测器的数据分布情况针对性地如下进行处理。The spatial information reconstruction method can be processed in accordance with the data distribution of the photon detector in a targeted manner as follows.
记空间信息分布的投影数据为y,待重建的空间信息为x,根据光子探测器的数据采集方式,有:The projection data of the spatial information distribution is y, and the spatial information to be reconstructed is x. According to the data collection method of the photon detector, there are:
y=Hx。 (2)y=Hx. (2)
此处H的维度是M×N,M=投影角度数×Ns,N是被重建的图像域像素数目。(2)式异于常规CT的空间信息重建之处在于系统矩阵H对应的是随机分布的探测器单元位置,而且可以是稀疏采样,即M<<N。即,每个角度下数据采集的探测器单元数目可以远小于总的探测器单元数目Ndet。当M<<N时,(2)是一个不定方程组求解问题。在本公开中,通过稀疏条件限制求解:Here, the dimension of H is M×N, M=the number of projection angles×Ns , and N is the number of pixels of the image domain to be reconstructed. (2) The spatial information reconstruction of the conventional CT is different in that the system matrix H corresponds to the randomly distributed detector unit position, and may be sparse sampling, that is, M<<N. That is, the number of detector units for data acquisition at each angle can be much smaller than the total number of detector units Ndet . When M<<N, (2) is an indefinite equation solving problem. In the present disclosure, the solution is limited by sparse conditions:
ε为由数据噪声情况确定的常数,(i,j,z)为被重建像素的三维位置坐标,Ф(i,j,z)是与成像视野的像素位置(i,j,z)有关的先验函数,可以预先设定或者通过探测器数据采样分布的概率密度函数计算得到(参见E.Y.Sidky,et al.,″Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,″Journal of X-Ray Science and Technology,vol.14,pp.119-139,2006)。ε is a constant determined by the data noise condition, (i, j, z) is the three-dimensional position coordinate of the reconstructed pixel, and Ф(i, j, z) is related to the pixel position (i, j, z) of the imaging field of view. The prior function can be pre-set or calculated from the probability density function of the detector data sample distribution (see EYSidky, et al., "Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT," Journal of X-Ray Science and Technology, vol. 14, pp. 119-139, 2006).
|·|为1阶范数:|x|=∑|xi|。|·| is a first-order norm: |x|=∑|xi |.
根据拉格朗日乘子法,(3)可以通过下式求解:According to the Lagrange multiplier method, (3) can be solved by:
此处λ是调节稀疏条件作用和数据保真性的因子,可以根据经验设置为一个常数。当Ns与Ndet接近时可以提高λ的取值。公式(3)和(4)的数值求解可以通过使用本领域公知的数值优化方法如ART-TV(参见Y.Sidky,et al.,″Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,″Journal of X-Ray Science and Technology,vol.14,pp.119-139,2006)或者ADMM(参见J.Yang and Y.Zhang,″Altemating Direction Algorithms for l1-problems in Compressive Sensing,″SIAM J.Sci.Comput.,vol.33,pp.250-278,2011)的方法迭代完成。Here λ is a factor that adjusts the sparse condition action and data fidelity and can be set to a constant based on experience. The value of λ can be increased when Ns is close to Ndet . The numerical solution of equations (3) and (4) can be performed by using numerical optimization methods well known in the art such as ART-TV (see Y. Sidky, et al., "Accurate image reconstruction from few-views and limited-angle data in divergent -beam CT, "Journal of X-Ray Science and Technology, vol. 14, pp. 119-139, 2006) or ADMM (see J. Yang and Y. Zhang, "Altemating Direction Algorithms for l1-problems in Compressive Sensing, The method of "SIAM J. Sci. Comput., vol. 33, pp. 250-278, 2011) is iteratively completed.
下面说明根据本公开一些实施例的利用通过前述方法采集的数据进行能谱CT图像重建的方法。A method of performing spectral CT image reconstruction using data acquired by the foregoing method in accordance with some embodiments of the present disclosure is described below.
根据一些实施例,采用图1所示的能谱CT系统架构,使用普通X光机,采集2个能量窗的数据,一圈采集Nθ个角度,每个角度使用同样参数的高斯分布的概率密度函数,但生成不同的样本,即不同的地址编码,进行探测器读数,且Ns=Ndet/l0。According to some embodiments, using the energy spectrum CT system architecture shown in FIG. 1, using a common X-ray machine, collecting data of two energy windows, acquiring Nθ angles in one circle, and using a probability of Gaussian distribution of the same parameter for each angle Density function, but generates different samples, ie different address codes, for detector readings, and Ns =Ndet /l0.
图4a和图4b示意性示出根据本公开一些示例实施例的投影数据比较,其中图4a示出普通能谱CT采集到的一个能量窗的投影数据,图4b示出根据随机分布的像素地址编码采集前者十分之一数据得到的一个能量窗下的投影数据。图5a和图5b示出根据本公开一些示例实施例进行能谱CT重建的结果示例,其中图5a是电子密度分布图,图5b是等效原子序数分布图。4a and 4b schematically illustrate projection data comparisons according to some example embodiments of the present disclosure, wherein FIG. 4a shows projection data of one energy window acquired by ordinary energy spectrum CT, and FIG. 4b shows pixel addresses according to random distribution. The coded projection data under an energy window obtained by collecting one tenth of the data. 5a and 5b illustrate examples of results of performing spectral CT reconstruction in accordance with some example embodiments of the present disclosure, wherein FIG. 5a is an electron density distribution map and FIG. 5b is an equivalent atomic number distribution map.
下面分别以预处理能谱CT重建的方法框架、及后处理能谱CT重建的方法框架为例进行说明。The following is an example of a method framework for pre-processing energy spectrum CT reconstruction and a method framework for post-processing energy spectrum CT reconstruction.
预处理方式进行能谱CT重建Pre-processing method for energy spectrum CT reconstruction
1)按照概率密度函数生成地址编码,按照地址编码进行数据采集。记录地址编码,生成对应的系统矩阵H。但本公开不限于此,例如,也可以根据记录的地址编码,在重建过程中实时计算系统矩阵的元素值。系统矩阵H的生成为本领域所公知,此处不再赘述。1) The address code is generated according to the probability density function, and the data is collected according to the address code. The address code is recorded to generate a corresponding system matrix H. However, the present disclosure is not limited thereto, and for example, it may be reconstructed based on the recorded address code.The element values of the system matrix are calculated in real time. The generation of the system matrix H is well known in the art and will not be described here.
2)计算Ф(i,j,z)。假设定义各个投影角度θ下探测器数据读取的采样点分布的概率密度函数为pdf(u,v,θ),u,v表示探测器位置,θ表示投影角度,可以设(但是不限于)2) Calculate Ф(i,j,z). Suppose that the probability density function of the sampling point distribution for reading the detector data under each projection angle θ is defined as pdf(u, v, θ), u, v is the detector position, and θ is the projection angle, which can be set (but not limited to)
Ф(i,j,z)=[H0T×pdf(u,v,θ)]-p,1/3<p<3,Ф(i,j,z)=[H0T ×pdf(u,v,θ)]-p , 1/3<p<3,
其中H0为采集全部探测器数据情况下的系统矩阵。此步骤通过一次反投影和结果求幂次运算完成。特殊情况下,例如无任何先验知识且均匀分布随机采样的情况下,可以设Where H0 is the system matrix in the case of collecting all detector data. This step is done by a back projection and a result power exponentiation. In special cases, for example, without any prior knowledge and evenly distributed random sampling, it can be set
Ф(i,j,z)=1。Ф(i,j,z)=1.
3)从2个能量窗的所有探测器数据进行解析得到分解系数投影。记得到的分解系数投影为Aτ。3) Analyze the decomposition coefficients from all the detector data of the two energy windows. The resulting decomposition coefficient is projected as Aτ .
用公式表示为根据下式求Aτ:Formulated to find Aτ according to the following formula:
Aτ={A1,τ,A2,τ,…,AL,τ}。Aτ ={A1,τ ,A2,τ ,...,AL,τ }.
τ为整数,对应第τ个分解项,在2个能量窗的情况下τ∈[1,2]。φτ(E)为第τ个分解基函数,Aτ为第τ个分解系数投影。该步骤可以使用本领域内现有的方法完成,例如双效应分解或者基材料分解。此处可选择φ1,2(E)为光电效应系数、康普顿散射系数,但本公开不限于此。τ is an integer corresponding to the τth decomposition term, and τ ∈ [1, 2] in the case of two energy windows. φτ (E) is the τth decomposition basis function, and Aτ is the τth decomposition coefficient projection. This step can be accomplished using methods available in the art, such as double effect decomposition or base material decomposition. Here, φ1,2 (E) may be selected as the photoelectric effect coefficient and the Compton scattering coefficient, but the present disclosure is not limited thereto.
4)使用3)得到的结果进行分解系数重建,得到分解系数图像,即对所有Aτ,使用如下方式重建获得aτ:4) Using the result obtained in 3), the decomposition coefficient is reconstructed to obtain an image of the decomposition coefficient, that is, for all Aτ , reconstruction is performed as follows to obtain aτ :
此处ε为一个小的阈值。此处给出一种具体实现过程示例,但不排斥其它迭代方法实现这一步骤。因为对每个aτ均使用下面的过程得到,所以在下面的叙述中省略了下标τ。Here ε is a small threshold. An example of a specific implementation process is given here, but other iterative methods are not excluded to implement this step. Since the following procedure is used for each aτ , the subscript τ is omitted in the following description.
a.设定迭代初值为a0;a. Set the initial value of the iteration to a0 ;
b.进行保真项更新,即计算b. Perform fidelity item update, ie calculation
c.非负性约束更新(这一步根据能谱分解方法而定,例如基材料分解,也可以省略)c. Non-negative constraint update (this step depends on the energy spectrum decomposition method, such as base material decomposition, can also be omitted)
d.先验约束更新:d. A priori constraint update:
i)
ii)进行Q次全变分最小化迭代(此处可以使用梯度下降法,但不限于这个方法),Q可为自行选择的整数,例如在5至100之间。Ii) Perform Q total variation minimization iterations (the gradient descent method can be used here, but not limited to this method), and Q can be a self-selected integer, for example between 5 and 100.
e.令再进行b-e步,直到满足收敛条件停止迭代。E. order Then perform the be step until the convergence condition is met to stop the iteration.
5)根据上述重建结果aτ,τ=1,...,Γ,合成计算单能量的μ(E),5) synthesizing the single energy μ(E) according to the above reconstruction result aτ , τ=1,...,Γ,
也可由μ(E)、aτ,τ=1,...,Γ估计被成像物的电子密度、等效原子序数分布图。在使用光电效应系数、康普顿散射系数作为φ(E)的情况下,The electron density and the equivalent atomic number distribution map of the imaged object can also be estimated from μ(E), aτ , τ=1, . In the case where the photoelectric effect coefficient and the Compton scattering coefficient are used as φ(E),
ρe=2a2。ρe = 2a2 .
后处理方式进行能谱CT重建Post-processing method for energy spectrum CT reconstruction
1)按照概率密度函数生成地址编码,按照地址编码进行数据采集。记录地址编码,生成对应的系统矩阵H。但本公开不限于此,例如,也可以根据记录的地址编码,在重建过程中实时计算系统矩阵的元素值。系统矩阵H的生成为本领域所公知,此处不再赘述。1) The address code is generated according to the probability density function, and the data is collected according to the address code. The address code is recorded to generate a corresponding system matrix H. However, the present disclosure is not limited thereto, and for example, the element values of the system matrix may be calculated in real time during the reconstruction process according to the recorded address encoding. The generation of the system matrix H is well known in the art and will not be described here.
2)计算Ф(i,j,z)。例如定义各个投影角度θ下探测器数据读取的采样点分布的概率密度函数为pdf(u,v,θ),可以设(但是不限于)2) Calculate Ф(i,j,z). For example, the probability density function of the sampling point distribution for reading the detector data read under each projection angle θ is pdf(u, v, θ), which can be set (but not limited to)
Ф(i,j,z)=[H0T×pdf(u,v,θ)]-p,1/3<p<3,Ф(i,j,z)=[H0T ×pdf(u,v,θ)]-p , 1/3<p<3,
其中H0为采集全部探测器数据情况下的系统矩阵。此步骤通过一次反投影和结果求幂次运算完成。特殊情况下,例如无任何先验知识且均匀分布随机采样的情况下,可以设Where H0 is the system matrix in the case of collecting all detector data. This step is done by a back projection and a result power exponentiation. In special cases, for example, without any prior knowledge and evenly distributed random sampling, it can be set
Ф(i,j,z)=1。Ф(i,j,z)=1.
3)从2个能量窗的所有探测器数据进行两个线衰减系数的空间分布重建。记对采集到的数据进行归一化和负对数处理得到的投影值为:3) Spatial distribution reconstruction of the two line attenuation coefficients from all detector data of the two energy windows. The projection values obtained by normalizing and negative logarithmically processing the collected data are:
对所有gk,使用如下方式重建获得μk:For all gk , use the following method to reconstruct to obtain μk :
其中ε为一个小的阈值。此处给出一种具体实现过程示例,但不排斥其它迭代方法实现这一步骤。因为对每个μk均使用下面的过程得到,所以在下面的叙述中省略了下标k。Where ε is a small threshold. An example of a specific implementation process is given here, but other iterative methods are not excluded to implement this step. Since the following procedure is used for each μk , the subscript k is omitted in the following description.
a.设定迭代初值为μ0;a. Set the initial value of the iteration to μ0 ;
b.进行保真项更新,即计算b. Perform fidelity item update, ie calculation
c.非负性约束更新c. Non-negative constraint update
d.先验约束更新:d. A priori constraint update:
i)
ii)进行Q次全变分最小化迭代(此处可以使用梯度下降法,但不限于这个方法),Q可为自行选择的整数,例如在5至100之间。Ii) Perform Q total variation minimization iterations (the gradient descent method can be used here, but not limited to this method), and Q can be a self-selected integer, for example between 5 and 100.
e.令再进行b-e步,直到满足收敛条件停止迭代。E. order Then perform the be step until the convergence condition is met to stop the iteration.
4)根据上述重建结果μk,k=1,...,K,对被重建的图像域像素逐个像素点进行能谱信息分解得到分解系数aτ(i,j,z),即求解线性方程组:4) According to the above reconstruction result μk , k=1, . . . , K, the spectral information is decomposed pixel by pixel from the reconstructed image domain pixel to obtain a decomposition coefficient aτ (i, j, z), that is, to solve the linearity equation set:
τ为整数,对应第τ个分解项。Ek表示第k个能量窗的等效光子能量。对所有像素点完成求解后,得到aτ。φτ(E)为第τ个分解基函数。在使用(但不限于)光电效应系数、康普顿散射系数作为φτ(E)的情况下,可以得到等效原子序数和电子密度分布图:τ is an integer corresponding to the τth decomposition term. Ek represents the equivalent photon energy of the kth energy window. After solving all the pixels, aτ is obtained. Φτ (E) is the τth decomposition basis function. In the case of using (but not limited to) the photoelectric effect coefficient and the Compton scattering coefficient as φτ (E), an equivalent atomic number and electron density distribution map can be obtained:
ρe=2a2。ρe = 2a2 .
图6示出根据本公开一实施例的能谱CT结构示意图,其包括具有预定义的随机分布像素的探测器。6 shows a schematic diagram of an energy spectrum CT structure including a detector having predefined randomly distributed pixels, in accordance with an embodiment of the present disclosure.
根据本公开的发明构思,可以采用具有预定义的随机分布像素的探测器。例如,参照图3a和3b,虚线网格处可不设置探测器单元(或像素),仅在深灰色小方格处设置探测器单元(或像素)。即,可按照多个随机分布的像素地址编码设置光子探测器的多个像素。像素地址编码可利用设定的像素地址分布的概率密度函数而生成,如前所述。In accordance with inventive concepts of the present disclosure, a detector having predefined randomly distributed pixels can be employed. For example, referring to Figures 3a and 3b, the detector unit (or pixel) may not be provided at the dotted grid, and the detector unit (or pixel) is only provided at the dark gray square. That is, a plurality of pixels of the photon detector can be set in accordance with a plurality of randomly distributed pixel address encodings. Pixel address coding can be generated using a probability density function of the set pixel address distribution, as previously described.
参照图6,该示例中没有采用一块完整的面阵光子计数探测器进行CT扫描实时编码采样的方式,而是采用另一种方式实现对投影数据的编码采集。如图6所示的CT系统结构,光子计数探测器以预先确定好的采样模式,例如利用设定的像素地址分布的概率密度函数而生成的采样模式,分布在一个环形探测器605的支架上。此外还有一与该环形探测器支架同轴的环形X射线源610。其中,环形X射线源610可以使用多个普通X射线源沿圆环排列的方式,也可以采用第五代电子束CT的环形靶组成。这种系统结构没有旋转结构,可以使用更高的速率采集投影数据,克服物体运动的模糊,提高CT重建图像的时间分辨率。Referring to FIG. 6, in this example, a complete area array photon counting detector is not used for CT scanning real-time encoding sampling, but another method is adopted for encoding and collecting projection data. As shown in the CT system structure of FIG. 6, the photon counting detector is distributed on a support of a
对于这种结构,其数据处理方式基本与上述方式相同,只有在计算Ф(i,j,z)时,由于探测器环上光子计数像素已经是预先定义好的,所以可以提前计算好存储下来,在重建时直接使用。For this structure, the data processing method is basically the same as the above method. Only when calculating Ф(i, j, z), since the photon counting pixels on the detector ring are already predefined, it can be calculated and stored in advance. , used directly when rebuilding.
通过以上的详细描述,本领域的技术人员易于理解,根据本公开实施例的系统和方法具有以下优点中的一个或多个。From the above detailed description, those skilled in the art will readily appreciate that systems and methods in accordance with embodiments of the present disclosure have one or more of the following advantages.
本公开通过一种新的探测器设计方式和对应的能谱CT重建方法,提出了一种新的能谱CT成像系统。通过低采样率随机编码方式读取光子计数探测器的像素数据,降低探测器的电荷共享效应和探测器的数据读取速率需求。The present disclosure proposes a new energy spectrum CT imaging system through a new detector design method and corresponding energy spectrum CT reconstruction method. The pixel data of the photon counting detector is read by low sampling rate random encoding, which reduces the charge sharing effect of the detector and the data reading rate requirement of the detector.
可以通过优化的能谱CT重建方法,获得常规系统固有的分辨率。The resolution inherent in conventional systems can be obtained by an optimized energy spectrum CT reconstruction method.
可以同时实现高分辨率、小数据量和大成像视野。High resolution, small data volume and large imaging field of view can be achieved simultaneously.
可以通过调整探测器数据采样分布的概率函数,灵活实现扫描视野内差异分辨率成像。The differential resolution imaging within the scanning field of view can be flexibly implemented by adjusting the probability function of the detector data sample distribution.
通过以上的实施例的描述,本领域的技术人员易于理解,本公开实施例可以通过硬件实现,也可以通过软件结合必要的硬件的方式来实现。因此,本公开实施例的技术方案可以以软件产品的形式体现出来,该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM,U盘,移动硬盘等)中,包括若干指令用以使得一台计算设备(可以是个人计算机、服务器、移动终端、或者网络设备等)执行根据本公开实施例的方法。Through the description of the above embodiments, those skilled in the art can easily understand that the embodiments of the present disclosure can be implemented by hardware or by software in combination with necessary hardware. Therefore, the technical solution of the embodiments of the present disclosure may be embodied in the form of a software product, which may be stored in a non-volatile storage medium (which may be a CD-ROM, a USB flash drive, a mobile hard disk, etc.), including a plurality of instructions. A method for causing a computing device (which may be a personal computer, server, mobile terminal, or network device, etc.) to perform a method in accordance with an embodiment of the present disclosure.
本领域技术人员可以理解,附图只是示例实施例的示意图,附图中的模块或流程并不一定是实施本公开所必须的,因此不能用于限制本公开的保护范围。It will be understood by those skilled in the art that the drawings are only a schematic diagram of the exemplary embodiments, and the modules or processes in the drawings are not necessarily required to implement the disclosure, and therefore are not intended to limit the scope of the disclosure.
本领域技术人员可以理解上述各模块可以按照实施例的描述分布于装置中,也可以进行相应变化位于不同于本实施例的一个或多个装置中。上述实施例的模块可以合并为一个模块,也可以进一步拆分成多个子模块。It will be understood by those skilled in the art that the above various modules may be distributed in the device according to the description of the embodiments, or the corresponding changes may be located in one or more devices different from the embodiment. The modules of the above embodiments may be combined into one module, or may be further split into multiple sub-modules.
以上具体地示出和描述了本公开的示例性实施例。应该理解,本公开不限于所公开的实施例,相反,本公开意图涵盖包含在所附权利要求的精神和范围内的各种修改和等效布置。The exemplary embodiments of the present disclosure have been particularly shown and described above. It is to be understood that the invention is not to be construed as being limited
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510163516.2 | 2015-04-08 | ||
| CN201510163516.2ACN106153647B (en) | 2015-04-08 | 2015-04-08 | Energy spectral CT imaging system and method for data acquisition and reconstruction of energy spectral CT image |
| Publication Number | Publication Date |
|---|---|
| WO2016161844A1true WO2016161844A1 (en) | 2016-10-13 |
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| PCT/CN2016/073918CeasedWO2016161844A1 (en) | 2015-04-08 | 2016-02-17 | Energy-spectrum ct imaging system, data acquisition method and energy-spectrum ct image reconstructing method |
| Country | Link |
|---|---|
| CN (1) | CN106153647B (en) |
| WO (1) | WO2016161844A1 (en) |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112529980A (en)* | 2020-12-14 | 2021-03-19 | 重庆师范大学 | Multi-target finite angle CT image reconstruction method based on maximum minimization |
| CN112666194A (en)* | 2020-12-22 | 2021-04-16 | 上海培云教育科技有限公司 | Virtual digital DR image generation method and DR virtual simulation instrument |
| CN113358674A (en)* | 2021-04-01 | 2021-09-07 | 西安交通大学 | Neutron resonance CT imaging system and method designed for reinforced concrete member |
| CN115541633A (en)* | 2022-09-20 | 2022-12-30 | 明峰医疗系统股份有限公司 | Energy spectrum CT acquisition device and method |
| CN116125517A (en)* | 2023-02-13 | 2023-05-16 | 成都理工大学 | Single-pixel radiation imaging method and system based on rotation measurement |
| CN116258673A (en)* | 2022-12-28 | 2023-06-13 | 上海联影医疗科技股份有限公司 | Image reconstruction method, system, electronic equipment and storage medium of spectral CT |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106559676A (en)* | 2016-11-25 | 2017-04-05 | 中北大学 | A kind of ciphered compressed of image and decryption decompression method |
| CN109924998A (en)* | 2019-03-22 | 2019-06-25 | 上海联影医疗科技有限公司 | Medical imaging procedure and photon counting power spectrum CT imaging device |
| CN110051387B (en)* | 2019-04-11 | 2020-05-19 | 华中科技大学 | Ray theory-based ultrasonic CT image reconstruction method and system |
| CN111289544A (en)* | 2020-02-25 | 2020-06-16 | 沈阳先进医疗设备技术孵化中心有限公司 | CT (computed tomography) equipment and parameter configuration method of detector array of CT equipment |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040101087A1 (en)* | 2002-11-27 | 2004-05-27 | Jiang Hsieh | Methods and apparatus for generating CT scout images |
| US20050271293A1 (en)* | 2004-06-04 | 2005-12-08 | Zhengrong Ying | Method of and system for destreaking the photoelectric image in multi-energy computed tomography |
| CN102695074A (en)* | 2012-06-08 | 2012-09-26 | 青岛海信电器股份有限公司 | Three-dimensional image display method and device |
| CN203149136U (en)* | 2012-12-31 | 2013-08-21 | 清华大学 | Multi-energy CT imaging system |
| CN103472074A (en)* | 2013-06-19 | 2013-12-25 | 清华大学 | CT imaging system and method |
| CN103559699A (en)* | 2013-11-18 | 2014-02-05 | 首都师范大学 | Multi-energy-spectrum CT image reconstruction method based on projection estimation |
| CN103900931A (en)* | 2012-12-26 | 2014-07-02 | 首都师范大学 | Multi-energy-spectrum CT imaging method and imaging system |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040101087A1 (en)* | 2002-11-27 | 2004-05-27 | Jiang Hsieh | Methods and apparatus for generating CT scout images |
| US20050271293A1 (en)* | 2004-06-04 | 2005-12-08 | Zhengrong Ying | Method of and system for destreaking the photoelectric image in multi-energy computed tomography |
| CN102695074A (en)* | 2012-06-08 | 2012-09-26 | 青岛海信电器股份有限公司 | Three-dimensional image display method and device |
| CN103900931A (en)* | 2012-12-26 | 2014-07-02 | 首都师范大学 | Multi-energy-spectrum CT imaging method and imaging system |
| CN203149136U (en)* | 2012-12-31 | 2013-08-21 | 清华大学 | Multi-energy CT imaging system |
| CN103472074A (en)* | 2013-06-19 | 2013-12-25 | 清华大学 | CT imaging system and method |
| CN103559699A (en)* | 2013-11-18 | 2014-02-05 | 首都师范大学 | Multi-energy-spectrum CT image reconstruction method based on projection estimation |
| Title |
|---|
| HAO, JIA ET AL.: "Multi-energy X-ray Imaging Technique and Its Application in Computed Tomography", CT THEORY AND APPLICATIONS, vol. 20, no. 1, 31 March 2011 (2011-03-31)* |
| SHEN, LE ET AL., STRUCTURAL PRIOR ENHANCED COMPRESSED SENSING FOR CT RECONSTRUCTION WITH INCOMPLETE DATA, 31 December 2013 (2013-12-31)* |
| XING, YUXIANG ET AL., A GENERAL ADAPTIVE DECOMPOSITION METHOD FOR MULTI-ENERGY SPECTRAL CT, 31 December 2013 (2013-12-31)* |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112529980A (en)* | 2020-12-14 | 2021-03-19 | 重庆师范大学 | Multi-target finite angle CT image reconstruction method based on maximum minimization |
| CN112529980B (en)* | 2020-12-14 | 2022-11-08 | 重庆师范大学 | A Multi-target Finite-Angle CT Image Reconstruction Method Based on Minimization |
| CN112666194A (en)* | 2020-12-22 | 2021-04-16 | 上海培云教育科技有限公司 | Virtual digital DR image generation method and DR virtual simulation instrument |
| CN112666194B (en)* | 2020-12-22 | 2022-12-20 | 上海培云教育科技有限公司 | Virtual digital DR image generation method and DR virtual simulation instrument |
| CN113358674A (en)* | 2021-04-01 | 2021-09-07 | 西安交通大学 | Neutron resonance CT imaging system and method designed for reinforced concrete member |
| CN113358674B (en)* | 2021-04-01 | 2024-05-24 | 西安交通大学 | Neutron resonance CT imaging system and method designed for reinforced concrete member |
| CN115541633A (en)* | 2022-09-20 | 2022-12-30 | 明峰医疗系统股份有限公司 | Energy spectrum CT acquisition device and method |
| CN116258673A (en)* | 2022-12-28 | 2023-06-13 | 上海联影医疗科技股份有限公司 | Image reconstruction method, system, electronic equipment and storage medium of spectral CT |
| CN116125517A (en)* | 2023-02-13 | 2023-05-16 | 成都理工大学 | Single-pixel radiation imaging method and system based on rotation measurement |
| Publication number | Publication date |
|---|---|
| CN106153647A (en) | 2016-11-23 |
| CN106153647B (en) | 2021-04-13 |
| Publication | Publication Date | Title |
|---|---|---|
| WO2016161844A1 (en) | Energy-spectrum ct imaging system, data acquisition method and energy-spectrum ct image reconstructing method | |
| CN104240270B (en) | CT imaging methods and system | |
| CN103472074B (en) | Ct imaging system and method | |
| Wu et al. | Image-domain material decomposition for spectral CT using a generalized dictionary learning | |
| CN106530366B (en) | Power spectrum CT image rebuilding method and power spectrum CT imaging system | |
| JP2016536032A (en) | Concatenated reconstruction of electron density images | |
| JP6615509B2 (en) | Reconstruction apparatus, X-ray computed tomography apparatus and reconstruction method | |
| JP2015226753A (en) | X-ray CT apparatus and control method thereof | |
| US20240374225A1 (en) | Systems and methods for energy bin downsampling | |
| US20210262947A1 (en) | Method and device for acquiring tomographic image data by oversampling, and control program | |
| Zhang et al. | Spectral CT image-domain material decomposition via sparsity residual prior and dictionary learning | |
| Li et al. | Multienergy cone-beam computed tomography reconstruction with a spatial spectral nonlocal means algorithm | |
| Hsieh et al. | Improving pulse detection in multibin photon-counting detectors | |
| WO2014107651A1 (en) | System and method for ultra-high resolution tomorgraphic imaging | |
| Zhang et al. | Reconstruction method for DECT with one half-scan plus a second limited-angle scan using prior knowledge of complementary support set (Pri-CSS) | |
| CN113506355A (en) | Scattering correction method, device, imaging system and computer readable storage medium | |
| Zhang et al. | Multi-energy CT reconstruction using tensor nonlocal similarity and spatial sparsity regularization | |
| Jumanazarov et al. | Significance of the spectral correction of photon counting detector response in material classification from spectral x-ray CT | |
| JP6462397B2 (en) | X-ray computed tomography apparatus and image reconstruction method | |
| Zhang et al. | CT image reconstruction algorithms: A comprehensive survey | |
| US9770221B1 (en) | Imaging using multiple energy levels | |
| WO2024173939A1 (en) | Customized cadmium zinc telluride (czt) detection system and methods | |
| Zhu et al. | Photon allocation strategy in region-of-interest tomographic imaging | |
| CN118297882A (en) | Photon counting type energy spectrum computed tomography multi-material decomposition method and device | |
| CN114916950B (en) | High-spatial-resolution energy spectrum CT image reconstruction method based on multilayer flat panel detector |
| Date | Code | Title | Description |
|---|---|---|---|
| 121 | Ep: the epo has been informed by wipo that ep was designated in this application | Ref document number:16776019 Country of ref document:EP Kind code of ref document:A1 | |
| NENP | Non-entry into the national phase | Ref country code:DE | |
| 122 | Ep: pct application non-entry in european phase | Ref document number:16776019 Country of ref document:EP Kind code of ref document:A1 |