A method, device, equipment and storage medium for spectral CT image reconstruction

Abstract

The present invention discloses a spectral CT image reconstruction method, device, equipment and storage medium, belonging to the field of image processing. The method includes: determining a target information extraction matrix based on a mass attenuation coefficient and a spectral data combination; determining a residual data combination between a projection data combination corresponding to a current base image combination and an actual projection data combination in the spectral dimension; if the residual data combination does not meet the iteration termination condition, determining an intermediate image combination based on the residual data combination; for each base material, extracting corresponding image information from each intermediate image using a target information extraction matrix to generate a target image corresponding to the current base material; updating the current base image combination based on all target images, returning to the step of determining the projection data combination, until the residual data combination meets the iteration termination condition, and using the current base image combination as the target base image combination. The embodiment of the present invention can quickly and accurately reconstruct a target base image combination based on projection data with inconsistent scanning geometry.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image processing, and in particular to a spectral CT image reconstruction method, device, equipment and storage medium. Background Art

[0002] Image reconstruction, as the name implies, is to reconstruct the projection data obtained by the imaging device to obtain an image, where the pixels of the image usually correspond to certain physical quantities of the scanned sample (target object). Compared with traditional CT (Computed Tomography), spectral CT is a more advanced imaging method. It collects projection data under multiple X-ray energy spectra and then reconstructs a spectral image that can reflect richer physical information inside the sample. Due to the nonlinearity of the imaging principle and the scanning geometry inconsistency caused by the imaging device, spectral image reconstruction needs to overcome the problem of scanning geometry inconsistency.

[0003] The existing technology can use the method of interpolating the projection data in advance to convert the projection data with inconsistent scanning geometry into the detection projection data with consistent scanning geometry, and then reconstruct the image of the detection projection data, but this often introduces errors, resulting in low quality of the reconstructed energy spectrum image.

[0004] In summary, it is difficult for existing technologies to quickly and accurately reconstruct high-quality energy spectrum images for projection data with inconsistent scanning geometry. Summary of the invention

[0005] The present invention provides a spectral CT image reconstruction method, device, equipment and storage medium to solve the problem that the prior art cannot directly reconstruct high-quality spectral images using projection data with inconsistent scanning geometry.

[0006] According to one aspect of the present invention, a method for reconstructing a spectral CT image is provided, comprising:

[0007] Obtaining a combination of energy spectrum data corresponding to Q energy spectra and a combination of actual projection data, and a combination of mass attenuation coefficients corresponding to D basis materials and a combination of a current basis image, where Q is equal to D, and D is greater than or equal to 2;

[0008] Determine a target information extraction matrix according to the mass attenuation coefficient combination and the energy spectrum data combination;

[0009] During the iterative image reconstruction process, determining the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension;

[0010] If the residual data combination does not meet the iteration termination condition, determining an intermediate image combination based on the residual data combination, and for each base material in the D base materials, extracting corresponding image information from each intermediate image of the intermediate image combination using the target information extraction matrix to generate a target image corresponding to the current base material;

[0011] The current base image combination is updated according to all target images, and the step of determining the projection data combination corresponding to the current base image combination in the energy spectrum dimension is returned until the residual data combination meets the iteration termination condition, and the current base image combination is used as the target base image combination.

[0012] According to another aspect of the present invention, there is provided a spectral CT image reconstruction device, comprising:

[0013] An acquisition module, used to acquire a combination of energy spectrum data corresponding to Q energy spectra and a combination of actual projection data, and a combination of mass attenuation coefficients corresponding to D basis materials and a combination of a current basis image, wherein Q is equal to D, and D is greater than or equal to 2;

[0014] A matrix determination module, used to determine a target information extraction matrix according to the mass attenuation coefficient combination and the energy spectrum data combination;

[0015] A residual data determination module, used to determine the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension during the image iterative reconstruction process;

[0016] an information extraction module, configured to determine an intermediate image combination based on the residual data combination if the residual data combination does not meet the iteration termination condition, and extract corresponding image information from each intermediate image of the intermediate image combination using the target information extraction matrix for each base material in the D base materials, so as to generate a target image corresponding to the current base material;

[0017] An iterative return module is used to update the current base image combination according to all target images, and return to the step of determining the projection data combination corresponding to the current base image combination in the energy spectrum dimension, until the residual data combination meets the iteration termination condition, and the current base image combination is used as the target base image combination.

[0018] According to another aspect of the present invention, an electronic device is provided, the electronic device comprising:

[0019] at least one processor; and

[0020] a memory communicatively connected to the at least one processor; wherein,

[0021] The memory stores a computer program executable by the at least one processor, and the computer program is executed by the at least one processor so that the at least one processor can perform the energy spectrum CT image reconstruction method described in any embodiment of the present invention.

[0022] According to another aspect of the present invention, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores computer instructions, and the computer instructions are used to enable a processor to implement the energy spectrum CT image reconstruction method described in any embodiment of the present invention when executed.

[0023] The technical solution of the energy spectrum CT image reconstruction method provided by the embodiment of the present invention determines the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension, so as to realize the determination of the residual data of each energy spectrum based on the actual projection data of each energy spectrum, which is helpful to improve the speed and accuracy of determining the residual data combination; by using the target information extraction matrix to extract the corresponding image information from the intermediate image combination, all the image information is integrated to generate the target image corresponding to the current base material, so as to improve the speed and accuracy of extracting image information from the intermediate image, thereby improving the speed and accuracy of target image generation, and the speed and accuracy of updating the current base image combination; moreover, in the image iteration process, the residual data combination is determined based on the actual projection data, and the corresponding image information is extracted from each intermediate information using the target information extraction matrix to generate the target image corresponding to the current base material, and all the target images are used to update the current base image combination. The coordinated use of technical features such as the above makes it possible to quickly and accurately reconstruct the energy spectrum data and the actual projection data with inconsistent scanning geometry without the need for image reconstruction based on the energy spectrum data and the actual projection data with consistent scanning geometry to obtain a high-quality target base image combination.

[0024] It should be understood that the contents described in this section are not intended to identify the key or important features of the embodiments of the present invention, nor are they intended to limit the scope of the present invention. Other features of the present invention will become easily understood through the following description. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the following briefly introduces the drawings required for use in the description of the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0026] Figure 1 is a flow chart of a spectral CT image reconstruction method provided in an embodiment of the present invention;

[0027] Figure 2 is another flow chart of a spectral CT image reconstruction method provided according to an embodiment of the present invention;

[0028] Figure 3 is another flow chart of a spectral CT image reconstruction method provided according to an embodiment of the present invention;

[0029] Figure 4 is an image combination result determined based on different algorithms according to an embodiment of the present invention;

[0030] Figure 5 is a curve showing a change in the relative error of a base image as a function of the number of iteration steps determined based on different algorithms according to an embodiment of the present invention;

[0031] Figure 6 is a curve showing a change in the relative error of projection data as a function of the number of iteration steps determined based on different algorithms according to an embodiment of the present invention;

[0032] Figure 7 is another image combination result determined based on different algorithms according to an embodiment of the present invention;

[0033] Figure 8 is a curve showing a change in relative error of another base image as a function of the number of iteration steps determined based on different algorithms according to an embodiment of the present invention;

[0034] Fig. 9 is another curve of a relative error of projection data determined based on different algorithms according to an embodiment of the present invention as a function of the number of iteration steps;

[0035] Fig.10 This is a schematic diagram of the structure of a spectral CT image reconstruction device provided according to an embodiment of the present invention:

[0036] Fig.11 It is a schematic diagram of the structure of an electronic device for implementing the energy spectrum CT image reconstruction method according to an embodiment of the present invention. DETAILED DESCRIPTION

[0037] In order to enable those skilled in the art to better understand the scheme of the present invention, the technical scheme in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work should fall within the scope of protection of the present invention.

[0038] It should be noted that the terms "first", "second", etc. in the specification and claims of the present invention and the above-mentioned drawings are used to distinguish similar objects, and are not necessarily used to describe a specific order or sequence. It should be understood that the data used in this way can be interchanged where appropriate, so that the embodiments of the present invention described herein can be implemented in an order other than those illustrated or described herein. In addition, the terms "including" and "having" and any variations thereof are intended to cover non-exclusive inclusions, for example, a process, method, system, product or device that includes a series of steps or units is not necessarily limited to those steps or units that are clearly listed, but may include other steps or units that are not clearly listed or inherent to these processes, methods, products or devices.

[0039] The imaging principle of CT is to use the absorption rate of the target material to X-rays to perform tomographic imaging on it, so as to obtain the internal structural information of the target material. Physically, this information is expressed by the distribution of the linear attenuation coefficient in space, which corresponds to the numerical value of the CT image pixel (or CT value). The linear attenuation coefficient of the sample is not only related to the spatial position, but also to the photon energy to which the sample is exposed. In order to reconstruct the linear attenuation coefficient of the sample, spectral CT uses the known energy spectrum and relevant information of the base material to remodel the data model, which is specifically divided into the following two aspects: On the one hand, spectral CT decomposes the sample into (regards it as) a combination of multiple base materials. If we divide the energy range under consideration into M energy intervals, then the linear attenuation coefficient of the sample in the mth energy interval, the pixel value in the ith pixel block is decomposed into the following form accordingly:

[0040]

[0041] Among them, b dm represents the mass attenuation coefficient of the dth basis material in the mth energy interval, f di Represents the value of the d-th basis material in the ith pixel block / voxel block. For example, assuming 50keV is in the 5th energy interval, μ 5i represents the linear attenuation coefficient of the sample in the ith pixel block / voxel block at 50keV. At this time, we also call is a virtual monoenergetic image at 50keV. In addition, if We call f d The constructed vector is the dth basis image. In spectral CT, according to the decomposition of the basis material, solving the unknown linear attenuation coefficient (or virtual monoenergetic image) is equivalent to solving the basis image.

[0042] On the other hand, the rays emitted by an X-ray tube usually contain photons of multiple energies. The distribution of photons of different energies in an X-ray is called the energy spectrum of the X-ray. Traditional CT is equivalent to scanning a sample using only single-energy X-rays, that is, assuming that all X-rays contain photons of only one energy. Spectral CT uses X-rays of different energy spectra to scan samples and can obtain multiple sets of observation data. For example, dual-energy CT uses two X-rays with different energy spectra (referred to as low-energy spectrum and high-energy spectrum) to scan samples.

[0043] The embodiment of the present invention uses represents the energy spectrum q, then It represents the proportion of photons in the mth energy interval to all photons. In summary, without considering noise and other physical factors, the imaging process of spectral CT can be modeled as the following data model:

[0044]

[0045] Where q is the energy spectrum identifier, j = 1,…,J [q] is the X-ray identification under the corresponding energy spectrum, It is the data model of the X-ray identified as j in the energy spectrum identified as q. The length (or area, or other related values) of the X-ray j of the energy spectrum denoted as q passing through the pixel block / voxel block i is determined by the geometric parameters corresponding to the X-ray j in the scanning device. The scanning geometry corresponding to the scanning device here can be a two-dimensional parallel beam, a fan beam, a three-dimensional parallel beam, a cone beam, a spiral cone beam, etc. is the discrete X-ray transformation vector, i.e., the scanning geometry information corresponding to the energy spectrum q (the energy spectrum marked as q), represents the column vector of all basis images concatenated. In addition, we call is the base image f d Projection data under X-ray j of energy spectrum q.

[0046] In the problem of spectral CT image reconstruction, the actual projection data (detection data), the X-ray energy spectrum distribution used, and the X-ray scanning geometry information are often known to reconstruct the base image. represents the actual projection data corresponding to X-ray j under energy spectrum q, then mathematically, the energy spectrum CT image reconstruction problem is equivalent to solving the following nonlinear equations:

[0047]

[0048] in, is a data model, and the data model is used to determine the projection data corresponding to the X-ray identified as j under the energy spectrum identified as q. It can be understood that K [q]It is used to determine the projection data corresponding to all X-rays under the energy spectrum marked as q.

[0049] For dual-energy CT, existing X-ray scanning and data acquisition devices mainly include dual-source dual-detector CT, fast low / high voltage switching CT, multi-layer detector CT, etc. Among these devices, only the scanning geometry of low-energy projection data and high-energy projection data acquired by multi-layer detector CT is consistent. Therefore, in the field of spectral CT imaging, the problem of inconsistent scanning geometry is ubiquitous. It is understandable that if an algorithm can reconstruct images for projection data with inconsistent scanning geometry, then it must also be able to reconstruct images for projection data with consistent scanning geometry. Therefore, the method described in this embodiment can reconstruct images for projection data acquired by various forms of data acquisition devices.

[0050] This embodiment rearranges the energy spectrum CT image reconstruction problem and combines the related equations of the same energy spectrum rays, thus obtaining the following Q nonlinear equation groups:

[0051]

[0052] Where q = 1, ..., Q corresponds to Q energy spectra. The discrete X-ray transformation of energy spectrum q is The matrix formed is A [q] ,Right now It includes the scanning geometry information of all actual projection data under the energy spectrum q. It should be noted that after the matrix row transformation, if A [i] =A [j] , where i and j are different and are respectively greater than or equal to 1 and less than or equal to Q, then this situation belongs to consistent scanning geometry; otherwise, it belongs to inconsistent scanning geometry. The scanning geometry in this embodiment can be a scanning geometry that appears in scanning modes such as a two-dimensional parallel beam, a fan beam geometry, and a three-dimensional parallel beam, a cone beam, or a spiral cone beam.

[0053] Figure 1 The flowchart of the energy spectrum CT image reconstruction method provided by the embodiment of the present invention is applicable to the situation where a high-quality base image is reconstructed quickly and accurately based on the projection data with consistent or inconsistent scanning geometry. The method can be executed by an energy spectrum CT image reconstruction device, which can be implemented in the form of hardware and / or software, and the energy spectrum CT image reconstruction device can be configured in the processor of an electronic device. Figure 1 As shown, the method includes:

[0054] S110, obtaining a combination of energy spectrum data corresponding to Q energy spectra and a combination of actual projection data, and a combination of mass attenuation coefficients corresponding to D basis materials and a combination of a current basis image, wherein Q is equal to D, and D is greater than or equal to 2.

[0055] The actual projection data combination is the projection data combination actually collected by the detector, which includes the actual projection data corresponding to Q energy spectra {g [q]}; The energy spectrum data combination includes Q energy spectrum data, and the energy spectrum data corresponding to each energy spectrum can be expressed as The mass attenuation coefficient of each base material in the mass attenuation coefficient combination can be expressed as {b dm}, where q is the energy spectrum identifier, d is the base material identifier, and m is the energy interval identifier. The energy spectrum data combination can be expressed as:

[0056]

[0057] Where M is the total number of energy intervals.

[0058] Similarly, the mass attenuation coefficient combination can be expressed as:

[0059]

[0060] In one embodiment, the energy spectrum data combination includes energy spectrum data corresponding to at least two energy spectra, such as high energy spectrum data and low energy spectrum data, or high energy spectrum data, medium energy spectrum data and low energy spectrum data.

[0061] The actual projection data corresponding to the at least two energy spectra in the actual projection data combination may be acquired based on a synchronous acquisition method or an asynchronous acquisition method.

[0062] It should be noted that the energy spectrum data corresponding to each energy spectrum in the energy spectrum data combination in the embodiment of the present invention is the energy spectrum data participating in the reconstruction of the energy spectrum CT image. The number of energy spectra corresponding to the actual energy spectrum data combination actually collected is greater than or equal to the number of base materials. When reconstructing the energy spectrum CT image, the energy spectrum data corresponding to Q energy spectra can be selected from the actual energy spectrum data combination based on the random selection method, and the energy spectrum data corresponding to the Q energy spectra can be used as the energy spectrum data combination; the energy spectrum data corresponding to the Q energy spectra that meet the energy spectrum constraint conditions in the actual projection data combination can also be used as the energy spectrum data combination; at least two energy spectrum data combinations can also be determined based on the actual energy spectrum data combination, and then the energy spectrum data combination that meets the set conditions can be selected from the at least two energy spectrum data combinations. Among them, the energy spectrum constraint conditions can be the maximum energy spectrum interval condition, the minimum energy spectrum interval condition, etc.

[0063] The current base image combination is the base image combination corresponding to the current iteration round in the image iterative reconstruction process. The current base image combination includes the current base images corresponding to each base material. Exemplarily, for the energy spectrum data combination obtained by high and low energy spectrum methods and the actual projection data combination, usually only two base images can be reconstructed, such as a water-based image and an iodine-based image, or a calcium-based image and a water-based image, or a water-based image and a bone-based image.

[0064] In the first iteration, the current base image combination is the initial base image combination. In one embodiment, the initial base image is a given base image, for example, the pixel value of each pixel is zero, and the resolution is a set resolution.

[0065] S120, determining a target information extraction matrix according to a combination of mass attenuation coefficients and a combination of energy spectrum data.

[0066] In one embodiment, the product of the mass attenuation coefficient combination and the energy spectrum data combination is used as an intermediate information extraction matrix; and the inverse matrix of the intermediate information extraction matrix is ​​used as a target information extraction matrix.

[0067] Specifically, define the mass attenuation coefficient combination as B, and the energy spectrum data combination as S, then the intermediate information extraction matrix is: The target information extraction matrix is: It can be understood that the target information extraction matrix includes sub-matrices corresponding to each energy spectrum in the Q energy spectra. The sub-matrices corresponding to each energy spectrum are configured to extract corresponding image information from the intermediate image corresponding to the corresponding energy spectrum. The target information extraction matrix in this embodiment can ensure the convergence of this iterative reconstruction algorithm, and the better convergence can make the method described in this embodiment have better universality, which can be used to perform iterative image reconstruction of energy spectrum data and actual projection data collected in different scenes.

[0068] S130. During the iterative image reconstruction process, determine the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension.

[0069] The projection data combination includes projection data corresponding to each energy spectrum, and the projection data corresponding to each energy spectrum is the projection data of the current base image combination under the corresponding energy spectrum determined based on the algorithm.

[0070] For each energy spectrum, the data model corresponding to the current energy spectrum is used to determine the projection data corresponding to the current base image combination under the current energy spectrum. The projection data corresponding to the intermediate image combination under all energy spectra are aggregated into a projection data combination; the difference between the projection data combination and the projection data of the corresponding energy spectrum in the actual projection data combination is determined to obtain the residual data corresponding to each energy spectrum; the residual data corresponding to all energy spectra are aggregated into a residual data combination.

[0071] The residual data corresponding to each current base image can reflect the data deviation of its projection data compared with the corresponding actual projection data, that is, reflect the image quality of each current base image. Specifically, for any current base image, if the norm of its corresponding residual data is large, then its image quality is low; if the norm of its corresponding residual data is small, then its image quality is high.

[0072] This embodiment uses actual projection data as a benchmark and residual data as a basis to judge the image quality of each current base image, thereby ensuring the accuracy of image quality judgment and thus ensuring the accuracy of the iterative loop; at the same time, determining the residual data in the energy spectrum dimension can increase the speed of determining the residual data.

[0073] S140. If the residual data combination does not meet the iteration termination condition, an intermediate image combination is determined based on the residual data combination, and for D base materials, a target information extraction matrix is ​​used to extract corresponding image information from each intermediate image of the intermediate image combination to generate a target image corresponding to the current base material.

[0074] In one embodiment, the iteration termination condition is that the number of iterations of the algorithm reaches the maximum number of iterations, or the norm of the residual data corresponding to each base material in the residual data combination is less than a set threshold.

[0075] Specifically, if the number of iterations corresponding to the residual data combination reaches the maximum number of iterations, or the norm of the residual data corresponding to each basis material in the residual data combination is less than the set threshold, it is determined that the iteration termination condition is met; otherwise, it is determined that the iteration termination condition is not met.

[0076] The maximum number of iterations and the set threshold can be determined based on experience and image quality requirements, such as being set to 20 times, 60 times, 1E-3, 1E-4, etc.

[0077] In one embodiment, image reconstruction is performed on the residual data corresponding to each energy spectrum in the residual data combination to obtain an intermediate image corresponding to each energy spectrum; and all the intermediate images are aggregated to obtain an intermediate image combination. The image reconstruction method in this embodiment is a single-energy CT image reconstruction method, which can specifically be a filtered back-projection algorithm, an algebraic reconstruction method, a reconstruction algorithm based on an optimization model, or a deep learning image reconstruction algorithm.

[0078] In one embodiment, according to the correspondence between the pre-created scanning geometry identifier and the single-energy CT image reconstruction algorithm, the single-energy CT image reconstruction algorithm corresponding to the scanning geometry identifier of each residual data is determined, and the image is reconstructed on the corresponding residual data based on the single-energy CT image reconstruction algorithm to obtain the corresponding intermediate image. Exemplarily, the scanning geometry identified as a three-dimensional cone beam corresponds to the FDK (Feldkamp-Davis-Kress) algorithm. This embodiment can improve the accuracy of each intermediate image.

[0079] Each residual data in the residual data combination can be expressed as Res [q] , where q is the energy spectrum identifier. The residual data combination can be expressed as: The intermediate image combination can be expressed as: Among them, A [q]-1 Indicates the single-energy CT reconstruction algorithm corresponding to the residual data marked as q.

[0080] After the intermediate image combination is determined, the target information extraction matrix is ​​used to extract the corresponding image information from each intermediate image in the intermediate image combination based on block matrix multiplication to obtain the target image combination, where the target image combination can be expressed as: in, is the symbol for block matrix multiplication.

[0081] Specifically, the number of sub-matrices included in the target information matrix is ​​the same as the number of energy spectra, so when using block matrix multiplication, each sub-matrix of the target information matrix is ​​configured to extract image information corresponding to an intermediate image, and then synthesize the image information extracted from each intermediate image to generate a target image corresponding to the current base material. The target images corresponding to each base material are summarized to obtain a target image combination. Therefore, the introduction of the target information extraction matrix achieves the technical effect of synthesizing the intermediate images corresponding to each energy spectrum to generate the target image corresponding to the current base material, and improves the accuracy of the target image determination.

[0082] In one embodiment, determining at least two energy spectrum data combinations based on the actual energy spectrum data combination, and selecting the energy spectrum data combination that meets the set conditions to determine the target image can be achieved by the following steps: determining at least two target information extraction matrices according to the mass attenuation coefficient combination and each energy spectrum data combination in the at least two energy spectrum data combinations; in the image iterative reconstruction process, determining the projection data combination corresponding to the current base image combination in the energy spectrum dimension, and the residual data combination between the projection data combination and the actual projection data combination; if the residual data combination does not meet the iteration termination condition, determining the intermediate image combination based on the residual data combination; for each base material in the D base materials The target intermediate image combinations corresponding to the target information extraction matrices are selected from the intermediate image combinations, wherein the target intermediate image combinations having a corresponding relationship correspond to the same energy spectrum identification combination as the target information extraction matrix; for each target intermediate image combination in all target intermediate image combinations, the corresponding target information extraction matrix is ​​used to extract corresponding image information from each intermediate image of the current target intermediate image combination to generate a target image to be evaluated corresponding to the current base material; the evaluation data of each target image to be evaluated corresponding to the current base material is determined, and the target image to be evaluated corresponding to the evaluation data corresponding to the highest image quality is used as the target image corresponding to the current base material. In this embodiment, for each target intermediate image combination in at least two target intermediate image combinations, the corresponding target information extraction matrix is ​​used to extract corresponding image information from each intermediate image of the current target intermediate image combination, and the target image to be evaluated corresponding to the current base material is determined based on the extracted image information; among the at least two target images to be evaluated, the target image to be evaluated corresponding to the evaluation data corresponding to the highest image quality is determined, and the target image to be evaluated is used as the target image corresponding to the current base material, thereby improving the accuracy of the target image determined in each iteration round, thereby reducing the iteration rounds in the image reconstruction process.

[0083] Specifically, in one embodiment, energy spectrum data of Q energy spectra are randomly selected from the actual energy spectrum data combination, and the energy spectrum data of the Q energy spectra are used as an energy spectrum data combination; and the step is repeated to obtain at least two energy spectrum data combinations. This embodiment can simply and quickly determine the required number of energy spectrum data combinations.

[0084] In one embodiment, a combination determination method is adopted to determine all energy spectrum data groups corresponding to the actual energy spectrum data combination. Specifically, determine The corresponding energy spectrum data combination mode, and the energy spectrum data combination corresponding to various energy spectrum data combination modes. This embodiment can quickly and accurately determine all energy spectrum data combinations corresponding to the actual energy spectrum data combination.

[0085] It should be noted that, in the actual energy spectrum data combination and each energy spectrum data combination, the arrangement order of the energy spectrum data of each energy spectrum is the same as the size order of the corresponding energy spectrum.

[0086] According to the combination of mass attenuation coefficients and each energy spectrum data combination, the corresponding target information extraction matrix is ​​determined. The target information extraction matrix is ​​used to extract image information of the intermediate image corresponding to the energy spectrum corresponding to the sub-matrix.

[0087] Exemplarily, the energy spectrum data combination includes energy spectrum data corresponding to the first energy spectrum, energy spectrum data corresponding to the third energy spectrum, and energy spectrum data corresponding to the fourth energy spectrum. Correspondingly, the target information extraction matrix includes a sub-matrix corresponding to the first energy spectrum, a sub-matrix corresponding to the third energy spectrum, and a sub-matrix corresponding to the fourth energy spectrum; then the target information extraction matrix is ​​used to extract image information of the intermediate image corresponding to the first energy spectrum, information of the intermediate image corresponding to the third energy spectrum, and image information of the intermediate image corresponding to the fourth energy spectrum. If a target information extraction matrix includes a first sub-matrix and a third sub-matrix, wherein the first sub-matrix corresponds to the first energy spectrum, and the third sub-matrix corresponds to the third energy spectrum; then the target intermediate image combination corresponding to the target information extraction matrix includes a first intermediate image and a third intermediate image, wherein the first intermediate image corresponds to the first energy spectrum, and the third intermediate image corresponds to the third energy spectrum.

[0088] S150, update the current base image combination according to all target images, and return to the step of determining the projection data combination corresponding to the current base image combination in the energy spectrum dimension, until the residual data combination meets the iteration termination condition, and use the current base image combination as the target base image combination.

[0089] In one embodiment, if each intermediate image in the intermediate image combination is reconstructed by each residual data in the residual data combination, the sum of the corresponding pixel values ​​of each current base image and the corresponding target image in the current base image combination is determined to obtain an image merging result; all image merging results are used as the updated current base image combination, and the updated current base image combination is used as the current base image combination in the next image iteration round, as follows:

[0090]

[0091] Among them, k is the current iteration round, k+1 is the next iteration round, and D is the total number of base materials.

[0092] If the number of iterations corresponding to the residual data combination is less than the maximum number of iterations, and the norm of each residual data in the residual data combination is greater than or equal to the set threshold, then return to the step of determining the projection data combination corresponding to the current basis image combination in the energy spectrum dimension, and continue the iterative calculation; if the number of iterations corresponding to the residual data combination is greater than or equal to the maximum number of iterations, or the norm of each residual data in the residual data combination is less than the set threshold, then the current basis image combination is output as the target basis image combination.

[0093] The technical solution of the energy spectrum CT image reconstruction method provided by the embodiment of the present invention determines the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension, realizes the determination of the residual data of each energy spectrum based on the actual projection data of each energy spectrum, and helps to improve the speed and accuracy of determining the residual data combination; by using the target information extraction matrix to extract the corresponding image information from the intermediate image combination, all the image information is integrated to generate the target image corresponding to the current base material, the speed and accuracy of extracting image information from the intermediate image is improved, thereby improving the speed and accuracy of target image generation, as well as the speed and accuracy of updating the current base image combination; moreover, in the image iteration process, the residual data combination is determined based on the actual projection data, and the corresponding image information is extracted from each intermediate information using the target information extraction matrix to generate the target image corresponding to the current base material, and all the target images are used to update the current base image combination. The coordinated use of technical features such as the above makes it possible to quickly and accurately reconstruct the energy spectrum data and the actual projection data with inconsistent scanning geometry to obtain a high-quality target base image combination without the need for image reconstruction based on the energy spectrum data and the actual projection data with consistent scanning geometry.

[0094] Figure 2 This is a flow chart of a spectral CT image reconstruction method provided in an embodiment of the present invention. This embodiment is used to refine the determination step of the target information extraction matrix in the above embodiment. Figure 2 As shown, the method includes:

[0095] S210. Obtain the energy spectrum data combination and the actual projection data combination corresponding to Q energy spectra, the mass attenuation coefficient combination and the current basis image combination corresponding to D basis materials, and the constant matrix corresponding to each energy spectrum, Q is equal to D, D is greater than or equal to 2, and the constant matrix includes matrix blocks corresponding to each energy spectrum in the Q energy spectra, each matrix block is a non-zero matrix, and each element is greater than or equal to 0.

[0096] S2201. For each energy spectrum data in the energy spectrum data combination, determine the scale coefficient of the current energy spectrum data according to the mass attenuation coefficient corresponding to each base material in the mass attenuation coefficient combination and the matrix block corresponding to the current energy spectrum.

[0097] In one embodiment, for the current energy spectrum, for each base material, the product of the mass attenuation coefficient corresponding to the current base material and the matrix block corresponding to the current energy spectrum is determined, and the product is used as the first data; the sum of the first data corresponding to all base materials is used as the second data; the exponential result with the set data as the base and the second data as the power is calculated, and the exponential result is used as the scale coefficient. The set data can be data such as e or 10.

[0098] S2202: reorganize each energy spectrum data according to each energy spectrum data in the energy spectrum data combination and a scale coefficient corresponding to each energy spectrum data.

[0099] The product of each energy spectrum data in the energy spectrum data combination and the corresponding scale coefficient is determined, and the product is used as the corresponding energy spectrum data after reorganization.

[0100] The scale factor is used to increase the difference between energy spectrum data corresponding to different energy spectra. The scale factor can be set to be greater than 1 or less than 1.

[0101] By processing the energy spectrum data corresponding to each energy spectrum in the energy spectrum data combination through the scale coefficient, the difference between the energy spectrum data corresponding to different energy spectra can be improved, thereby improving the accuracy and speed of the corresponding image information extracted by the target information extraction matrix from each intermediate image, as well as the accuracy and determination speed of the target base image combination.

[0102] S2203, normalizing each reorganized energy spectrum data to obtain a normalized reorganized energy spectrum data combination.

[0103] Determine the sum of the energy spectrum data corresponding to all the reorganized energy spectra to obtain the total energy spectrum data; determine the ratio of the energy spectrum data corresponding to each reorganized energy spectrum to the total energy spectrum data, and use the ratio as the initial information extraction matrix corresponding to the current energy spectrum. It can be understood that the normalized reorganized energy spectrum data combination includes Q initial information extraction matrices.

[0104] The initial information extraction matrix is ​​defined as: It can be determined by the following formula:

[0105]

[0106] in, represents the energy spectrum data of energy spectrum q in the mth energy interval, b dm represents the mass attenuation coefficient of the dth basis material in the mth energy interval, is the matrix block corresponding to the dth basis material and energy spectrum q in the initial information extraction matrix; D is the number of basis materials, and M is the number of energy intervals.

[0107] S2204, determining an intermediate information extraction matrix according to the mass attenuation coefficient combination and the normalized reshaped energy spectrum data combination, and using the inverse matrix of the intermediate information extraction matrix as the target information extraction matrix.

[0108] The product of the initial information extraction matrix and the mass attenuation coefficient is used as the intermediate information extraction matrix, which achieves the purpose of simply and quickly determining the intermediate information extraction matrix. The inverse matrix of the intermediate information extraction matrix is ​​used as the target information extraction matrix, so that the intermediate information extraction matrix and the target information extraction matrix are inverse matrices of each other.

[0109] It can be understood that different constant matrices correspond to different target information extraction matrices; for the same intermediate image, different target information extraction matrices extract different image information and generate different target images, and different target images correspond to different current base image combination update results, so different target information extraction matrices will correspond to different iteration speeds and different convergence speeds.

[0110] S230. During the iterative image reconstruction process, determine the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension.

[0111] S240. If the residual data combination does not meet the iteration termination condition, an intermediate image combination is determined based on the residual data combination, and for each of the D base materials, a target information extraction matrix is ​​used to extract corresponding image information from each intermediate image of the intermediate image combination to generate a target image corresponding to the current base material.

[0112] S250, updating the current base image combination according to all target images, and returning to the step of determining the projection data combination corresponding to the current base image combination in the energy spectrum dimension, until the residual data combination meets the iteration termination condition, and taking the current base image combination as the target base image combination.

[0113] The technical solution provided by the embodiment of the present invention realizes determining different target information extraction matrices by selecting different constant matrices; for the same intermediate image, different target information extraction matrices extract different image information, resulting in different update results of the current base image combination, so different target information extraction matrices will correspond to different iteration speeds and different convergence speeds. Compared with the target information extraction matrix used in the aforementioned embodiment, the target information extraction matrix described in this embodiment has a higher iteration speed.

[0114] Figure 3 This is a flow chart of a spectral CT image reconstruction method provided in an embodiment of the present invention. This embodiment provides another implementation method for determining an intermediate image combination based on a residual data combination. Figure 3 As shown, the method includes:

[0115] S310, obtaining a combination of energy spectrum data corresponding to Q energy spectra and a combination of actual projection data, and a combination of mass attenuation coefficients corresponding to D basis materials and a combination of a current basis image, wherein Q is equal to D, and D is greater than or equal to 2.

[0116] S320, determining a target information extraction matrix according to a combination of mass attenuation coefficients and a combination of energy spectrum data.

[0117] S330. During the iterative image reconstruction process, determine the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension.

[0118] S3401. Taking the sum of the residual data combination and the current iteration projection data combination as the updated current iteration projection data combination.

[0119] The iterative projection data is an update object of this embodiment, and this embodiment aims to update the current base image combination by updating the iterative projection data.

[0120] The iterative projection data combination in the first iteration round includes initial iterative projection data corresponding to each energy spectrum, and each initial iterative projection data is determined based on a nonlinear Taylor expansion of the data model at an initial basis image combination of the corresponding energy spectrum.

[0121] Since the Taylor expansion is actually an approximation, it is to approximate (or fit) the original function by superposition of power functions at a certain specified point. The nonlinear Taylor expansion is a multi-order approximation with high accuracy. Therefore, this embodiment uses the nonlinear Taylor expansion of the data model at the initial basis image under the corresponding energy spectrum to determine the initial iterative projection data under the corresponding energy spectrum, which can improve the accuracy of the initial iterative projection data under the corresponding energy spectrum. It can be understood that the initial iterative projection data is an approximate data for reconstructing the actual projection data of the corresponding initial basis image.

[0122] In one embodiment, the nonlinear Taylor expansion includes projection data items, actual projection data items and intermediate projection data items; the result of the projection data items is the projection data corresponding to the initial basis image of the corresponding energy spectrum; the result of the actual projection data items is the actual projection data corresponding to the corresponding energy spectrum in the actual projection data combination; the result of the intermediate projection data items is determined based on the target data, and the target data includes the scanning geometry information of the corresponding energy spectrum, the initial basis images of each basis material, and the submatrix of each basis material corresponding to the corresponding energy spectrum in the intermediate information extraction matrix. See the following formula for details:

[0123]

[0124] Among them, (Λ0 ) [q] is the initial iterative projection data corresponding to the energy spectrum q, K [q] is the data model of the basis image under the energy spectrum q, f 0 is the initial basis image combination under the energy spectrum q, g [q] is the actual projection data corresponding to the energy spectrum q, A is a matrix block corresponding to the energy spectrum q and the basis material marked as d in the intermediate information extraction matrix; [q] is the scanning geometry matrix corresponding to the energy spectrum q, is an initial base image corresponding to the base material identified as d. The initial base image is a given base image, such as a set resolution image in which all pixel values ​​are zero.

[0125] Since the above nonlinear Taylor expansion includes actual projection data items and intermediate projection data items, the initial iterative projection data includes scanning geometry information corresponding to each energy spectrum.

[0126] In one embodiment, the nonlinear Taylor expansion further includes other data items to better approximate the actual projection data used to reconstruct the initial basis image combination under the corresponding energy spectrum.

[0127] S3402 . Reconstruct the images of the respective current iteration projection data in the updated current iteration projection data combination to obtain an intermediate image combination.

[0128] The image reconstruction algorithm described in the above embodiment can be used to reconstruct the image of each iterative projection data in the updated current iterative projection data combination, which will not be described in detail in this embodiment.

[0129] Since the updated current iteration projection data combination includes the current iteration projection data of the current iteration round, the intermediate image includes image information corresponding to the current iteration projection data and image information corresponding to the residual data.

[0130] S3403. For each base material among the D base materials, use the target information extraction matrix to extract corresponding image information from each intermediate image of the intermediate image combination to generate a target image corresponding to the current base material.

[0131] In one embodiment, the target image corresponding to each base material is determined by the following formula:

[0132]

[0133] Among them, (Λ k ) [q] is the current iteration projection data corresponding to the energy spectrum q in the kth iteration round, is the intermediate image corresponding to the energy spectrum marked as 1, This is the middle image corresponding to the energy spectrum marked as Q.

[0134] S3501. Combining results of all target images as updated current base image combination.

[0135]

[0136] After the target images corresponding to the base materials are determined, all the target images are aggregated into a combination result, and the combination result is used as the current base image combination.

[0137] In one embodiment, if the current iteration round is the first iteration round, the initial base image combination is updated based on the target information extraction matrix and the initial iteration projection data, and the updated initial base image combination is used as the current base image combination of the current iteration round; if the current iteration round is not the first iteration round, the current base image combination updated in the previous iteration round is used as the current base image combination of the current iteration round. Compared with directly determining the projection data combination corresponding to the current base image combination, and the residual data of the projection data combination compared to the actual projection data combination, this embodiment can reduce the number of iterations of image iterative reconstruction and reduce the time of image iterative reconstruction.

[0138] S3502, return to the step of determining the projection data combination corresponding to the current basis image combination in the energy spectrum dimension, until the residual data combination meets the iteration termination condition, and use the current basis image combination as the target basis image combination.

[0139] In the embodiment of the present invention, the sum of the current iterative projection data combination and the residual data combination is used as the updated current iterative projection data, thereby achieving the effect of using the residual data to correct the current iterative projection data; the process of using the corresponding image information extracted from each intermediate image to generate the corresponding target image includes: comprehensively extracting the image information from the intermediate images corresponding to each energy spectrum, and determining the target image corresponding to the corresponding base material; in this way, after all the target images are combined as the updated current base image, the technical effect of correcting the current base image combination is achieved. In summary, in this embodiment, one projection data correction and one image correction are performed in one iteration round, which improves the information correction speed, helps to reduce the number of iterations in the image iterative reconstruction process, and improves the image reconstruction speed.

[0140] To facilitate the understanding of the technical solution, this embodiment simulates a noise-free dual-energy CT image reconstruction experiment to illustrate that the energy spectrum CT image reconstruction method described in the above embodiment has higher accuracy and faster image reconstruction speed.

[0141] This embodiment includes a low energy spectrum of 80 kVp and a high energy spectrum of 140 kVp, wherein a 1 mm wide copper filter is added to the X-rays of the high energy spectrum. Assume that the torso phantom is made of water-based materials and bone-based materials, distributed in an area of ​​[-5 cm, 5 cm] × [-5 cm, 5 cm], and the pixel size is 128 × 128. Figure 4 As shown in the first row of . Parallel beam geometry is used for both low-energy spectrum X-rays and high-energy ones. Specifically, there are 384 scanning angles of low-energy spectrum X-rays, which are evenly distributed from 0 to 180 degrees, and there are 384 detectors, which are evenly distributed in the range of [-7.05cm, 7.05cm]. The detector for collecting the actual projection data of the high-energy spectrum is the same as the detector for collecting the actual projection data of the low-energy spectrum, but its 384 scanning angles are evenly distributed between 0.2344 and 180.2344 degrees. Since the scanning angle distributions corresponding to the actual projection data of the high-energy spectrum and the actual projection data of the low-energy spectrum are different, the scanning geometries of the two are inconsistent.

[0142] The energy spectrum data, set mass attenuation coefficient, scanning geometry information and real image are substituted into the imaging model to simulate and generate noise-free dual energy spectrum data. Among them, the real image is the target base image corresponding to the scanning conditions in the simulation experiment. The scanning conditions include energy spectrum energy, scanning geometry, energy spectrum data, mass attenuation coefficient, etc.

[0143] Specifically, in the following formula represents the base image of the water-based material at the k-th iteration, Represents the base image of the bone-based material at the k-th iteration.

[0144] Initialization algorithm, the initial base image is set to zero value image, that is The actual projection data corresponding to each energy spectrum {g [q]}、Energy spectrum data And the mass attenuation coefficient matrix B: = (b dm ). Let k←0, let the constant matrix Then the matrix W(0) = S, calculate and store the target information extraction matrix

[0145] Set each current base image to a zero-value image, calculate the projection data corresponding to each current base image, and the residual data between each projection data and the corresponding actual projection data, i.e., Res [1] and Res [2] :

[0146]

[0147] If the maximum number of iterations is reached or the norm of at least one residual data in the residual data combination is greater than the set threshold, the algorithm is terminated and the current base image combination is output as the target base image combination; otherwise, the filtered back projection algorithm is used to reconstruct the residual data to obtain two intermediate images, which are A and B respectively. [1]-1 (Res [1] )and Next, for each base material, the target information extraction matrix is ​​used to extract the corresponding image information from each intermediate image, and the target image corresponding to the corresponding base material is generated based on the corresponding image information extracted from each intermediate image. The sum of the pixel values ​​of each target image and the corresponding pixel in the corresponding current base image is calculated to generate the corresponding result image, and all the result images are summarized to obtain the updated current base image combination. The specific formula is as follows:

[0148]

[0149] The number of iteration steps is increased by 1, that is, k←k+1, and the iteration is continued, jumping to the above step of calculating the projection data corresponding to each current base image and the residual data between each projection data and the corresponding actual projection data.

[0150] The method described in the embodiment of the present invention, the existing iterative filtered back projection algorithm, the existing non-convex primal dual algorithm, the data domain interpolation method, and the nonlinear Kaczmarz method are all iterated for 60 steps to plot the relative error of the base image (see Figure 5 ) and the curve of the relative error of the projection data changing with the number of iterations (see Figure 6 ). Relative error of base image Defined as:

[0151]

[0152] Among them, f * is the real image, f k is the image determined in the kth iteration round, such as each basis image, the virtual monoenergetic image under the corresponding energy, ∥·∥ F Represents the Frobenius norm of a matrix.

[0153] Relative error of projection data Defined as:

[0154]

[0155] Among them, K [q] (f k ) is the projection data corresponding to the current basis image under the energy spectrum q in the kth iteration round, g [q]is the actual projection data corresponding to the energy spectrum q in the actual projection data combination, K [q] (f k )-g [q] It is the residual data corresponding to the current basis image under the energy spectrum q.

[0156] From the attached Figure 4 With attached Figure 5 It can be seen that in the image iterative reconstruction process of the method described in the embodiment of the present invention, the relative error of the base image and the relative error of the projection data both decrease with the increase of the number of iteration steps; in the image iteration process of the iterative filtering back projection algorithm, the relative error of the base image and the relative error of the projection data both decrease first and then increase sharply, and the algorithm does not converge in the end; in the image iteration process of the non-convex primal dual algorithm and the nonlinear Kaczmarz method, the relative error of the base image and the relative error of the projection data both decrease with the increase of the number of iteration steps, but the error reduction speed is much smaller than the error reduction speed corresponding to the method described in this embodiment. Although the data domain interpolation method converges quickly, the relative error of its final result is larger than the relative error of the method described in the embodiment of the present invention.

[0157] It was found through experiments that when the image reconstruction was performed based on the iterative filtering back projection algorithm, it still did not converge after the 60th iteration, so the desired image combination result was not obtained after the 60th iteration. Figure 4 As shown, the first row is the aforementioned real image, and the second, third, fourth, and fifth rows represent the image combination determined by this embodiment and the image combination results determined based on the non-convex primal dual algorithm, the data domain interpolation method, and the nonlinear Kaczmarz method; the image combination results include, from left to right, the water-based image, the bone-based image, the 60keV virtual monoenergetic image, and the 100keV virtual monoenergetic image. Figure 4 It can be seen that the image combination result determined by the embodiment of the present invention is almost the same as the real image combination result, and there is a significant deviation between the image combination result determined based on the non-convex primal dual algorithm, the data domain interpolation method, and the nonlinear Kaczmarz method and the real image combination result.

[0158] In order to better demonstrate the difference between the image combination result determined by the method described in the embodiment of the present invention and the image combination determined by the non-convex primal-dual algorithm, evaluation indicators such as normalized root mean square error, structural similarity, structural difference, and peak signal-to-noise ratio are introduced, among which the structural difference is defined as: the difference between 1 and the structural similarity, that is, 1-SSIM.

[0159] Table 1 shows the evaluation index display results provided by the embodiment of the present invention. Each cell in Table 1, from top to bottom, is normalized root mean square error data, structural difference data and peak signal-to-noise ratio data.

[0160] Among them, the normalized root mean square error is defined as: u * is the image to be compared, and the image in the image combination result determined in this embodiment or the image in the image combination result determined based on the non-convex primal-dual algorithm can be substituted here. For u * The corresponding real image.

[0161] The structural similarity is defined as: Among them, x i for u * The i-th window in y i is the i-th window in u, N is the total number of windows, and the window can be set to 8×8, 16×16, etc. u is the image to be compared. In addition, ssim in the above formula is defined as

[0162]

[0163] where μ x and μ y represent the means of x and y respectively, and Represent the variance of x and y respectively, where c1=c2=0.05, x and y are the image to be compared and the corresponding real image respectively.

[0164] The peak signal-to-noise ratio is defined as: Among them, u is the image to be compared, u * is the real image corresponding to u.

[0165] Table 1 Evaluation index data display table

[0166]

[0167] It can be seen from the above table that under the same experimental conditions, for the noise-free spectral CT image reconstruction problem of the same scale, the accuracy of the image index data determined by the embodiment of the present invention is much higher than the accuracy of the image index data determined based on the non-convex primal dual algorithm, the data domain interpolation method, and the nonlinear Kaczmarz method.

[0168] It is found through image iteration timing that each iteration of the method described in the embodiment of the present invention only takes 0.08 seconds, while each iteration of the iterative filtering back projection algorithm and the non-convex primal dual algorithm takes 15.81 seconds and 0.22 seconds respectively, and each iteration of the data domain interpolation method and the nonlinear Kaczmarz method takes 16.63 seconds and 0.66 seconds respectively. This shows that the energy spectrum CT image reconstruction method described in the embodiment of the present invention has a good balance between image quality and reconstruction efficiency, and can quickly and accurately reconstruct a high-quality target base image combination.

[0169] The experiment of simulating the reconstruction of dual-energy CT images containing noise is used to illustrate that the energy spectrum CT image reconstruction method described in the above embodiment can take into account both image accuracy and reconstruction efficiency.

[0170] The CT image of the torso phantom is used as the real image. At least two energy spectra include 80kVp (low energy) spectrum and 140kVp (high energy) spectrum. A 1mm wide copper filter is added to the X-ray of the high energy spectrum. Assuming that the torso phantom is composed of water-based materials and bone-based materials, distributed in an area of ​​[-5cm, 5cm]×[-5cm, 5cm], with a pixel size of 512×512, the real image combination corresponding to this scanning condition is as follows: Figure 7 The first row in . Both low-energy spectrum X-ray and high-energy spectrum X-ray adopt parallel beam geometry. Specifically, there are 1024 scanning angles for low-energy spectrum X-ray, which are evenly distributed from 0 to 180 degrees, and 1536 detectors are evenly distributed in the range of [-7.05cm, 7.05cm]. The detector corresponding to high-energy spectrum X-ray is the same as that corresponding to low-energy spectrum X-ray, but its 1024 scanning angles are evenly distributed between 0.08789 and 180.08789 degrees. Therefore, the scanning geometry of the actual projection data collected by the detector corresponding to the high-energy spectrum is different from that of the actual standard projection data collected by the detector corresponding to the low-energy spectrum. The actual projection data, mass attenuation coefficient, scanning geometry information and real image are substituted into the data model to simulate the generation of noise-free dual-energy spectrum data, and then Gaussian white noise is added to generate noisy dual-energy spectrum data.

[0171] The noisy dual-energy spectrum data were reconstructed using the method described in the embodiment of the present invention, the iterative filtering back projection algorithm and the non-convex primal dual algorithm, the data domain interpolation method, and the nonlinear Kaczmarz method. After 20 iterations, the iterative filtering back projection algorithm still did not converge.

[0172] Figure 8 A relative error curve diagram of another base image provided by an embodiment of the present invention; Fig. 9 Another projection data relative error curve diagram provided by an embodiment of the present invention. Figure 8 and Fig. 9 It can be seen that in the image iterative reconstruction process of the method described in the embodiment of the present invention, the relative error of the base image and the relative error of the projection data both decrease with the increase of the number of iteration steps; in the image iteration process of the iterative filtered back projection algorithm, the relative error of the base image and the relative error of the projection data both decrease first, then increase sharply, and finally do not converge; in the image iteration process of the non-convex primal dual algorithm, the data domain interpolation method and the nonlinear Kaczmarz method, the relative error of the base image and the relative error of the projection data both decrease with the increase of the number of iteration steps, but the relative error reduction rate of these algorithms is much smaller than the corresponding relative error reduction rate of the method described in this embodiment.

[0173] Combination Figure 5 , Figure 6 , Figure 8 and Fig. 9 It can be seen that, in the method described in the embodiment of the present invention, the relative error of the base image and the relative error of the projection data decrease faster in the reconstruction of the spectral CT image with noise than in the reconstruction of the spectral CT image without noise.

[0174] The image combination results determined by these three methods are as follows Figure 7 As shown: the first row represents the real image, the second, third, fourth and fifth rows respectively represent the image combination results determined by the method described in the embodiment of the present invention and the non-convex primal dual algorithm, the data domain interpolation method, and the nonlinear Kaczmarz method; the image combination results are, from left to right, a water-based image, a bone-based image, a 60keV virtual monoenergetic image, and a 100keV virtual monoenergetic image. Since the iterative filtered back projection algorithm did not converge and the resulting image pixels are all nan, its results are not shown here. Figure 7 It can be seen that the image reconstructed by the method described in the embodiment of the present invention is almost the same as the real image; the non-convex primal-dual algorithm converges very slowly and the result is very blurry, the data domain interpolation method converges faster, but the image quality obtained is not as good as the method described in the embodiment of the present invention. The nonlinear Kaczmarz method also converges slowly, and its resulting image contains serious artifacts.

[0175] The index data of the image combination result determined based on the method of the embodiment of the present invention and the image combination result determined based on the non-convex primal dual algorithm and the data domain interpolation method are shown in Table 2:

[0176] Table 2 Evaluation index data display table

[0177]

[0178] It can be seen from the above table that under the same experimental conditions, for the same scale of noisy spectral CT image reconstruction problem, the accuracy of each index data of the image determined by the method described in the embodiment of the present invention is much higher than the accuracy of each index data of the image determined based on the non-convex primal dual algorithm. In addition, it was found through iterative timing that when the method described in the embodiment of the present invention is implemented, each iteration only takes 0.82 seconds, while the iterative filtered back projection algorithm and the non-convex primal dual algorithm require 108.29 seconds and 2.31 seconds per iteration, respectively, and the data domain interpolation method and the nonlinear Kaczmarz method require 307.49 seconds and 9031.75 seconds per iteration, respectively. In summary, it can be seen that the method described in the embodiment of the present invention can well balance image quality and reconstruction efficiency, and therefore can quickly and accurately reconstruct a high-quality target base image combination.

[0179] Fig.10 The schematic diagram of the structure of the energy spectrum CT image reconstruction device provided by the embodiment of the present invention is as follows. Fig.10 As shown, the device comprises:

[0180] An acquisition module 61 is used to acquire a combination of energy spectrum data corresponding to Q energy spectra and a combination of actual projection data, and a combination of mass attenuation coefficients corresponding to D basis materials and a combination of a current basis image, where Q is equal to D, and D is greater than or equal to 2;

[0181] A matrix determination module 62, configured to determine a target information extraction matrix according to the mass attenuation coefficient combination and the energy spectrum data combination;

[0182] A residual data determination module 63 is used to determine the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension during the image iterative reconstruction process;

[0183] an information extraction module 64, configured to determine an intermediate image combination based on the residual data combination if the residual data combination does not meet the iteration termination condition, and extract corresponding image information from each intermediate image of the intermediate image combination using the target information extraction matrix for each base material in the D base materials, so as to generate a target image corresponding to the current base material;

[0184] The iterative return module 65 is used to update the current base image combination according to all target images, and return to the step of determining the projection data combination corresponding to the current base image combination in the energy spectrum dimension, until the residual data combination meets the iteration termination condition, and the current base image combination is used as the target base image combination.

[0185] In one embodiment, the matrix determination module 62 is specifically used for:

[0186] Taking the product of the mass attenuation coefficient combination and the energy spectrum data combination as an intermediate information extraction matrix;

[0187] The inverse matrix of the intermediate information extraction matrix is ​​used as the target information extraction matrix.

[0188] In one embodiment, the acquisition module 61 is further used for:

[0189] Obtain a constant matrix, wherein the constant matrix includes a matrix block corresponding to each energy spectrum in the Q energy spectra, each matrix block is a non-zero matrix, and each element is greater than or equal to 0;

[0190] The matrix determination module 62 is specifically used for:

[0191] For each energy spectrum data in the energy spectrum data combination, determining a scale coefficient of the current energy spectrum data according to the mass attenuation coefficient corresponding to each base material in the mass attenuation coefficient combination and a matrix block corresponding to the current energy spectrum;

[0192] Reorder each energy spectrum data according to each energy spectrum data in the energy spectrum data combination and a scale coefficient corresponding to each energy spectrum data;

[0193] Normalizing each reorganized energy spectrum data to obtain a normalized reorganized energy spectrum data combination;

[0194] An intermediate information extraction matrix is ​​determined according to the mass attenuation coefficient combination and the normalized and reshaped energy spectrum data combination, and an inverse matrix of the intermediate information extraction matrix is ​​used as a target information extraction matrix.

[0195] In one embodiment, the information extraction module 64 is used to:

[0196] Performing image reconstruction on each residual data in the residual data combination to obtain an intermediate image combination;

[0197] The iterative return module 65 is used to:

[0198] Determine the sum of corresponding pixel values ​​of each current base image and the corresponding target image in the current base image combination to obtain an image merging result;

[0199] The combined results of all images are used as the updated current base image combination.

[0200] In one embodiment, the information extraction module 64 is used to:

[0201] Taking the sum of the residual data combination and the current iteration projection data combination as the updated current iteration projection data combination;

[0202] Performing image reconstruction on each current iteration projection data in the updated current iteration projection data combination to obtain an intermediate image combination;

[0203] The iterative return module 65 is used to:

[0204] The combination results of all target images are used as the updated current base image combination;

[0205] The iterative projection data combination in the first iteration round includes initial iterative projection data corresponding to each energy spectrum, and the initial iterative projection data is determined based on a nonlinear Taylor expansion of the data model at the initial basis image combination corresponding to the energy spectrum.

[0206] In one embodiment, the nonlinear Taylor expansion includes a projection data item, an actual projection data item, and an intermediate projection data item;

[0207] The result of the projection data item is the projection data corresponding to the initial basis image of the corresponding energy spectrum;

[0208] The result of the actual projection data item is the actual projection data corresponding to the corresponding energy spectrum in the actual projection data combination;

[0209] The result of the intermediate projection data item is determined based on target data, and the target data includes scanning geometry information under the corresponding energy spectrum, initial basis images of each basis material, and a sub-matrix corresponding to each basis material in the intermediate information extraction matrix.

[0210] In one embodiment, the residual data determination module 63 is further configured to:

[0211] If the current iteration round is the first iteration round, the initial base image combination is updated based on the target information extraction matrix and the initial iteration projection data, and the updated initial base image combination is used as the current base image combination of the current iteration round;

[0212] If the current iteration round is not the first iteration round, the current base image combination updated in the previous iteration round is used as the current base image combination of the current iteration round.

[0213] The technical solution of the energy spectrum CT image reconstruction device provided by the embodiment of the present invention determines the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension, so as to realize the determination of the residual data of each energy spectrum based on the actual projection data of each energy spectrum, which is helpful to improve the speed and accuracy of determining the residual data combination; by using the target information extraction matrix to extract the corresponding image information from the intermediate image combination, all the image information is integrated to generate the target image corresponding to the current base material, so as to improve the speed and accuracy of extracting the image information from the intermediate image, thereby improving the speed and accuracy of generating the target image, and the speed and accuracy of updating the current base image combination; moreover, in the image iteration process, the residual data combination is determined based on the actual projection data, and the corresponding image information is extracted from each intermediate information using the target information extraction matrix to generate the target image corresponding to the current base material, and all the target images are used to update the current base image combination. The coordinated use of technical features such as the above makes it possible to quickly and accurately reconstruct the energy spectrum data and the actual projection data with inconsistent scanning geometry without the need for image reconstruction based on the energy spectrum data and the actual projection data with consistent scanning geometry to obtain a high-quality target base image combination.

[0214] The energy spectrum CT image reconstruction device provided in the embodiment of the present invention can execute the energy spectrum CT image reconstruction method provided in any embodiment of the present invention, and has the corresponding functional modules and beneficial effects of the execution method.

[0215] Fig.11 A schematic diagram of the structure of an electronic device 10 that can be used to implement an embodiment of the present invention is shown. The electronic device is intended to represent various forms of digital computers, such as laptop computers, desktop computers, workstations, personal digital assistants, servers, blade servers, mainframe computers, and other suitable computers. The electronic device can also represent various forms of mobile devices, such as personal digital processing, cellular phones, smart phones, wearable devices (such as helmets, glasses, watches, etc.) and other similar computing devices. The components shown herein, their connections and relationships, and their functions are merely examples and are not intended to limit the implementation of the present invention described and / or required herein.

[0216] like Fig.11 As shown, the electronic device 10 includes at least one processor 11, and a memory connected to the at least one processor 11, such as a read-only memory (ROM) 12, a random access memory (RAM) 13, etc., wherein the memory stores a computer program that can be executed by at least one processor, and the processor 11 can perform various appropriate actions and processes according to the computer program stored in the read-only memory (ROM) 12 or the computer program loaded from the storage unit 18 to the random access memory (RAM) 13. In the RAM 13, various programs and data required for the operation of the electronic device 10 can also be stored. The processor 11, the ROM 12, and the RAM 13 are connected to each other through a bus 14. An input / output (I / O) interface 15 is also connected to the bus 14.

[0217] A number of components in the electronic device 10 are connected to the I / O interface 15, including: an input unit 16, such as a keyboard, a mouse, etc.; an output unit 17, such as various types of displays, speakers, etc.; a storage unit 18, such as a disk, an optical disk, etc.; and a communication unit 19, such as a network card, a modem, a wireless communication transceiver, etc. The communication unit 19 allows the electronic device 10 to exchange information / data with other devices through a computer network such as the Internet and / or various telecommunication networks.

[0218] The processor 11 may be a variety of general and / or special processing components with processing and computing capabilities. Some examples of the processor 11 include, but are not limited to, a central processing unit (CPU), a graphics processing unit (GPU), various special artificial intelligence (AI) computing chips, various processors running machine learning model algorithms, a digital signal processor (DSP), and any appropriate processor, controller, microcontroller, etc. The processor 11 executes the various methods and processes described above, such as the spectral CT image reconstruction method.

[0219] In some embodiments, the spectral CT image reconstruction method may be implemented as a computer program, which is tangibly contained in a computer-readable storage medium, such as a storage unit 18. In some embodiments, part or all of the computer program may be loaded and / or installed on the electronic device 10 via the ROM 12 and / or the communication unit 19. When the computer program is loaded into the RAM 13 and executed by the processor 11, one or more steps of the spectral CT image reconstruction method described above may be performed. Alternatively, in other embodiments, the processor 11 may be configured to perform the spectral CT image reconstruction method in any other appropriate manner (e.g., by means of firmware).

[0220] Various implementations of the systems and techniques described above herein can be implemented in digital electronic circuit systems, integrated circuit systems, field programmable gate arrays (FPGAs), application specific integrated circuits (ASICs), application specific standard products (ASSPs), systems on chips (SOCs), load programmable logic devices (CPLDs), computer hardware, firmware, software, and / or combinations thereof. These implementations may include: being implemented in one or more computer programs that can be executed and / or interpreted on a programmable system including at least one programmable processor, which can be a special purpose or general purpose programmable processor that can receive data and instructions from a storage system, at least one input device, and at least one output device, and transmit data and instructions to the storage system, the at least one input device, and the at least one output device.

[0221] Computer programs for implementing the methods of the present invention may be written in any combination of one or more programming languages. These computer programs may be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing device, so that when the computer program is executed by the processor, the functions / operations specified in the flow chart and / or block diagram are implemented. The computer program may be executed entirely on the machine, partially on the machine, partially on the machine and partially on a remote machine as a stand-alone software package, or entirely on a remote machine or server.

[0222] In the context of the present invention, a computer-readable storage medium may be a tangible medium that may contain or store a computer program for use by or in combination with an instruction execution system, device or equipment. A computer-readable storage medium may include, but is not limited to, electronic, magnetic, optical, electromagnetic, infrared, or semiconductor systems, devices or equipment, or any suitable combination of the foregoing. Alternatively, a computer-readable storage medium may be a machine-readable signal medium. A more specific example of a machine-readable storage medium may include an electrical connection based on one or more lines, a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing.

[0223] To provide interaction with a user, the systems and techniques described herein may be implemented on an electronic device having: a display device (e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor) for displaying information to the user; and a keyboard and a pointing device (e.g., a mouse or trackball) through which the user can provide input to the electronic device. Other types of devices may also be used to provide interaction with the user; for example, the feedback provided to the user may be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback); and input from the user may be received in any form (including acoustic input, voice input, or tactile input).

[0224] The systems and techniques described herein may be implemented in a computing system that includes backend components (e.g., as a data server), or a computing system that includes middleware components (e.g., an application server), or a computing system that includes frontend components (e.g., a user computer with a graphical user interface or a web browser through which a user can interact with implementations of the systems and techniques described herein), or a computing system that includes any combination of such backend components, middleware components, or frontend components. The components of the system may be interconnected by any form or medium of digital data communication (e.g., a communication network). Examples of communication networks include: a local area network (LAN), a wide area network (WAN), a blockchain network, and the Internet.

[0225] A computing system may include a client and a server. The client and the server are generally remote from each other and usually interact through a communication network. The client and server relationship is generated by computer programs running on the corresponding computers and having a client-server relationship with each other. The server may be a cloud server, also known as a cloud computing server or cloud host, which is a host product in the cloud computing service system to solve the defects of difficult management and weak business scalability in traditional physical hosts and VPS services.

[0226] It should be understood that the various forms of processes shown above can be used to reorder, add or delete steps. For example, the steps described in the present invention can be executed in parallel, sequentially or in different orders, as long as the desired results of the technical solution of the present invention can be achieved, and this document does not limit this.

[0227] The above specific implementations do not constitute a limitation on the protection scope of the present invention. It should be understood by those skilled in the art that various modifications, combinations, sub-combinations and substitutions can be made according to design requirements and other factors. Any modification, equivalent substitution and improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A spectral CT image reconstruction method, characterized in that: include: Obtaining a combination of energy spectrum data corresponding to Q energy spectra and a combination of actual projection data, and a combination of mass attenuation coefficients corresponding to D basis materials and a combination of a current basis image, where Q is equal to D, and D is greater than or equal to 2; Determine a target information extraction matrix according to the mass attenuation coefficient combination and the energy spectrum data combination; During the iterative image reconstruction process, determining the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension; If the residual data combination does not meet the iteration termination condition, determining an intermediate image combination based on the residual data combination, and for each base material in the D base materials, extracting corresponding image information from each intermediate image of the intermediate image combination using the target information extraction matrix to generate a target image corresponding to the current base material; The current base image combination is updated according to all target images, and the step of determining the projection data combination corresponding to the current base image combination in the energy spectrum dimension is returned until the residual data combination meets the iteration termination condition, and the current base image combination is used as the target base image combination.

2. The method according to claim 1, characterized in that: The step of determining a target information extraction matrix according to the mass attenuation coefficient combination and the energy spectrum data combination includes: Taking the product of the mass attenuation coefficient combination and the energy spectrum data combination as an intermediate information extraction matrix; The inverse matrix of the intermediate information extraction matrix is ​​used as the target information extraction matrix.

3. The method according to claim 1, characterized in that: The step of acquiring the energy spectrum data combination corresponding to the Q energy spectra and the actual projection data combination, and the mass attenuation coefficient combination corresponding to the D basis materials and the current basis image combination, further includes: Obtain a constant matrix, wherein the constant matrix includes a matrix block corresponding to each energy spectrum in the Q energy spectra, each matrix block is a non-zero matrix, and each element is greater than or equal to 0; The step of determining a target information extraction matrix according to the mass attenuation coefficient combination and the energy spectrum data combination includes: For each energy spectrum data in the energy spectrum data combination, determining a scale coefficient of the current energy spectrum data according to the mass attenuation coefficient corresponding to each base material in the mass attenuation coefficient combination and a matrix block corresponding to the current energy spectrum; Reorder each energy spectrum data according to each energy spectrum data in the energy spectrum data combination and a scale coefficient corresponding to each energy spectrum data; Normalizing each reorganized energy spectrum data to obtain a normalized reorganized energy spectrum data combination; An intermediate information extraction matrix is ​​determined according to the mass attenuation coefficient combination and the normalized and reshaped energy spectrum data combination, and an inverse matrix of the intermediate information extraction matrix is ​​used as a target information extraction matrix.

4. The method according to claim 1, characterized in that: The determining of the intermediate image combination based on the residual data combination comprises: Performing image reconstruction on each residual data in the residual data combination to obtain an intermediate image combination; The updating of the current base image combination according to all target images comprises: Determine the sum of the pixel values ​​of each current base image and the corresponding target image to obtain an image merging result; The combined results of all images are used as the updated current base image combination.

5. The method according to claim 2 or 3, characterized in that: The determining of the intermediate image combination based on the residual data combination comprises: Taking the sum of the residual data combination and the current iterative projection data combination as the updated current iterative projection data combination; Performing image reconstruction on each current iteration projection data in the updated current iteration projection data combination to obtain an intermediate image combination; The updating of the current base image combination according to all target images comprises: The combination results of all target images are used as the updated current base image combination; The iterative projection data combination in the first iteration round includes initial iterative projection data corresponding to each energy spectrum, and the initial iterative projection data is determined based on a nonlinear Taylor expansion of the data model at the initial basis image combination corresponding to the energy spectrum.

6. The method according to claim 5, characterized in that The nonlinear Taylor expansion includes a projection data item, an actual projection data item and an intermediate projection data item; The result of the projection data item is the projection data corresponding to the initial basis image of the corresponding energy spectrum; The result of the actual projection data item is the actual projection data corresponding to the corresponding energy spectrum in the actual projection data combination; The result of the intermediate projection data item is determined based on target data, and the target data includes scanning geometry information under the corresponding energy spectrum, initial basis images of each basis material, and a sub-matrix corresponding to each basis material in the intermediate information extraction matrix.

7. The method according to claim 5, characterized in that Before determining the projection data combination corresponding to the current base image combination, the method further includes: If the current iteration round is the first iteration round, the initial base image combination is updated based on the target information extraction matrix and the initial iteration projection data, and the updated initial base image combination is used as the current base image combination of the current iteration round; If the current iteration round is not the first iteration round, the current base image combination updated in the previous iteration round is used as the current base image combination of the current iteration round.

8. A spectral CT image reconstruction device, characterized in that: include: An acquisition module, used to acquire a combination of energy spectrum data corresponding to Q energy spectra and a combination of actual projection data, and a combination of mass attenuation coefficients corresponding to D basis materials and a combination of a current basis image, wherein Q is equal to D, and D is greater than or equal to 2; A matrix determination module, used to determine a target information extraction matrix according to the mass attenuation coefficient combination and the energy spectrum data combination; A residual data determination module, used to determine the projection data combination corresponding to the current base image combination and the residual data combination between the projection data combination and the actual projection data combination in the energy spectrum dimension during the image iterative reconstruction process; an information extraction module, configured to determine an intermediate image combination based on the residual data combination if the residual data combination does not meet the iteration termination condition, and extract corresponding image information from each intermediate image of the intermediate image combination using the target information extraction matrix for each base material in the D base materials, so as to generate a target image corresponding to the current base material; An iterative return module is used to update the current base image combination according to all target images, and return to the step of determining the projection data combination corresponding to the current base image combination in the energy spectrum dimension, until the residual data combination meets the iteration termination condition, and the current base image combination is used as the target base image combination.

9. An electronic device, characterized in that: The electronic device comprises: at least one processor; and a memory communicatively connected to the at least one processor; wherein, The memory stores a computer program executable by the at least one processor, and the computer program is executed by the at least one processor so that the at least one processor can perform the spectral CT image reconstruction method according to any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores computer instructions, and the computer instructions are used to enable a processor to implement the energy spectrum CT image reconstruction method according to any one of claims 1 to 7 when executed.

Images