Energy spectrum CT image reconstruction method and storage medium

Abstract

The embodiment of the invention discloses an energy spectrum CT image reconstruction method and a storage medium. The method comprises the following steps: in the current iterative reconstruction process of an energy spectrum CT image, for each energy spectrum, determining transformation projection data corresponding to the energy spectrum, substituting the transformation projection data into a sub-problem pre-constructed for the energy spectrum, and solving the sub-problem after substitution based on a first constraint preset for the sub-problem, the transformation base image meets the first constraint, and the error between the transformation base image and the transformation projection data is minimum; and combining the transformation base images corresponding to the energy spectrums to obtain a combined base image, and substituting the combined base image into a pre-constructed optimization problem to reconstruct a current base image with a minimum error with the combined base image. According to the technical scheme provided by the embodiment of the invention, the energy spectrum CT image can be efficiently and accurately reconstructed, and particularly, the energy spectrum CT image can be efficiently and accurately reconstructed under a low-dose condition and / or a condition that scanning geometric parameters are inconsistent under different energy spectrums.

Description

Technical Field

[0001] The embodiments of the present invention relate to the technical field of medical image processing, and particularly to an energy spectrum CT image reconstruction method and a storage medium. Background Art

[0002] Compared with traditional Computed Tomography (CT), as an advanced imaging technology, energy spectrum CT has obvious advantages in suppressing artifacts and differentiating different substances, and thus has a wide range of applications in the fields of medical imaging and industrial imaging.

[0003] In energy spectrum CT, by using X-rays with different energy spectra to scan a target object, and then reconstructing an energy spectrum CT image based on the detected projection data, or reconstructing a basis material decomposition image (hereinafter referred to as a basis image for short). Moreover, in clinical applications, the X-ray dose is often reduced by reducing the tube voltage or tube current, or reducing the sampling angle, etc., so as to ensure the safety of the scanned target object.

[0004] However, reducing the tube voltage or tube current will result in more noise in the projection data, and reducing the sampling angle will result in sparser projection data. Through practical verification, projection data with more noise or sparser data will affect the reconstruction accuracy and reconstruction efficiency of energy spectrum CT images, which urgently need to be solved. Summary of the Invention

[0005] The embodiments of the present invention provide an energy spectrum CT image reconstruction method and a storage medium, which solve the problem of low reconstruction accuracy and reconstruction efficiency of energy spectrum CT images, and particularly solve the problem of low reconstruction accuracy and reconstruction efficiency of low-dose energy spectrum CT images.

[0006] According to one aspect of the present invention, an energy spectrum CT image reconstruction method is provided, including:

[0007] In the current iterative reconstruction process of an energy spectrum CT image, for each of at least two energy spectra, according to the projection matrix and the original projection data corresponding to the energy spectrum, and the previous basis image reconstructed in the previous iterative reconstruction process, determine the transformed projection data; substitute the transformed projection data and the projection matrix into the sub-problem pre-constructed for the energy spectrum, and solve the substituted sub-problem based on the first constraint preset for the sub-problem, so as to obtain a transformed basis image that satisfies the first constraint and has the smallest error from the transformed projection data; combine the transformed basis images corresponding to each energy spectrum to obtain a combined basis image, and substitute the combined basis image into the pre-constructed optimization problem to reconstruct the current basis image with the smallest error from the combined basis image, so as to realize the reconstruction of the energy spectrum CT image.

[0008] According to another aspect of the present invention, there is provided a computer-readable storage medium having stored thereon computer instructions for causing a processor to perform, when executed, the energy spectrum CT image reconstruction method provided in any embodiment of the present invention.

[0009] The technical solution of the embodiment of the present invention can efficiently and accurately reconstruct an energy spectrum CT image, especially in the case of low dose, an energy spectrum CT image can be efficiently and accurately reconstructed.

[0010] It should be understood that the content described in this part is not intended to identify the key or important features of the embodiments of the present invention, nor is it used to limit the scope of the present invention. Other features of the present invention will become easily understood through the following description. Description of the Drawings

[0011] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the following will briefly introduce the drawings required for the description of the embodiments. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can be obtained based on these drawings.

[0012] Figure 1 is a flowchart of an energy spectrum CT image reconstruction method provided according to an embodiment of the present invention;

[0013] Figure 2 is a flowchart of another energy spectrum CT image reconstruction method provided according to an embodiment of the present invention;

[0014] Figure 3 is a flowchart of yet another energy spectrum CT image reconstruction method provided according to an embodiment of the present invention;

[0015] Figure 4a is a graph showing the change of the relative error of the basis image with the number of iterations in yet another energy spectrum CT image reconstruction method provided according to an embodiment of the present invention under sparse noise-free data, corresponding to Experiment 1;

[0016] Figure 4b is a graph showing the change of the relative error of the projection data with the number of iterations in yet another energy spectrum CT image reconstruction method provided according to an embodiment of the present invention under sparse noise-free data, corresponding to Experiment 1;

[0017] Figure 5 is a schematic diagram of a head phantom image reconstructed by different algorithms in yet another energy spectrum CT image reconstruction method provided according to an embodiment of the present invention under sparse noise-free data, corresponding to Experiment 1;

[0018] Figure 6aIt is a graph showing the change of the relative error of the basis image with the number of iterations under sparse noise-free data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 2;

[0019] Figure 6b It is a graph showing the change of the relative error of the projection data with the number of iterations under sparse noise-free data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 2;

[0020] Figure 7 It is a schematic diagram of the XCAT phantom image reconstructed by different algorithms under sparse noise-free data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 2;

[0021] Figure 8a It is a graph showing the change of the relative error of the basis image with the number of iterations under sparse noisy data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 3;

[0022] Figure 8b It is a graph showing the change of the relative error of the projection data with the number of iterations under sparse noisy data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 3;

[0023] Figure 9a It is a graph showing the change of the relative error of the basis image with time under sparse noisy data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 3;

[0024] Figure 9b It is a graph showing the change of the relative error of the projection data with time under sparse noisy data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 3;

[0025] Figure 10a It is a graph showing the change of the relative error of the basis image between two adjacent steps with the number of iterations under sparse noisy data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 3;

[0026] Figure 10b It is a graph showing the change of the relative error of the projection data between two adjacent steps with the number of iterations under sparse noisy data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 3;

[0027] Figure 11 It is a schematic diagram of the lung CT image reconstructed by different algorithms under sparse noisy data in another spectral CT image reconstruction method provided by an embodiment of the present invention, corresponding to Experiment 3;

[0028] Figure 12 It is a structural block diagram of an energy spectrum CT image reconstruction device provided according to an embodiment of the present invention;

[0029] Figure 13 It is a structural block diagram of an electronic device for implementing the energy spectrum CT image reconstruction method according to an embodiment of the present invention. Specific embodiments

[0030] In order to enable those skilled in the art to better understand the solution of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0031] It should be noted that the terms "first", "second", etc. in the specification and claims of the present invention and the above drawings are used to distinguish similar objects, and do not necessarily need to describe a specific order or sequence. It should be understood that such data can be interchanged under appropriate circumstances so that the embodiments of the present invention described here can be implemented in an order other than those illustrated or described here. The situations of "target", "original", etc. are similar and will not be repeated here. In addition, the terms "comprising" and "having" and any variations thereof are intended to cover non-exclusive inclusion. For example, a process, method, system, product or device comprising a series of steps or units does not necessarily have to be limited to those steps or units clearly listed, but may include other steps or units not clearly listed or inherent to these processes, methods, products or devices.

[0032] Before introducing the embodiments of the present invention, the application scenario and the specific reasons for the problem that the related solutions have low reconstruction accuracy and reconstruction efficiency will be exemplarily described, so as to better understand how the energy spectrum CT image reconstruction solution proposed by the embodiments of the present invention improves the reconstruction accuracy and reconstruction efficiency.

[0033] Exemplarily, when X-rays are incident on a target object, the target object will absorb a part of the X-ray energy. By detecting the intensity attenuation change before and after the X-rays penetrate, the linear attenuation coefficient distribution of the substances inside the target object can be calculated, and thus a CT image can be reconstructed. Physically speaking, the linear attenuation coefficient of a substance is correlated with both its own composition and the incident photon energy.

[0034] For ease of processing, researchers usually need to decompose the linear attenuation coefficient based on physical effects (such as the photoelectric effect and Compton scattering) or material composition. Different from traditional CT which treats X-rays as a single-energy photon beam, spectral CT fully considers the spectral characteristics of X-rays. The rays emitted by the X-ray tube actually contain photons of various different energies, and the energy distribution is called the spectrum of X-rays. In particular, dual-energy CT uses two X-rays with different spectra (respectively called the low-energy spectrum and the high-energy spectrum) to scan the target object.

[0035] Available Represents the th spectrum, where Is the proportion of the number of photons in the th energy interval (or energy node) to the total number of photons. The data model of spectral CT can be expressed as:

[0036] ,

[0037] Where, Is the data noise generation operator, Represents all projection data or observed data, Represents all spectral CT imaging operators, defined as:

[0038] ,

[0039] Where, Represents the number of spectra, Represents the projection data of all detector units and all projection angles under the th spectrum, Is the imaging operator of the th spectrum, that is:

[0040] .

[0041] According to Beer's law, the imaging operator corresponding to the th X-ray of the th spectrum can be modeled as:

[0042] ,

[0043] Where, Corresponds to the index of all X-rays under the th spectrum. There are X-rays under the th spectrum. Is the index of the energy node, Is the number of energy nodes. Is the index of the decomposed base material, The quantity of the base materials. is the th base material's mass attenuation coefficient at the th energy node. It should be noted that this is a known physical quantity that can reflect the absorption ability of a substance to X-rays of a specific energy. is the th energy spectrum's th X-ray's path length passing through the th pixel block (or voxel) (e.g., length for 2D imaging and area for 3D imaging). is the th pixel value in the th base image. There are a total of pixel values. Then is the vector corresponding to the th base image, that is , where represents the transpose of the vector. Physically, the th base image represents the discrete spatial distribution of the th base material. represents the vector formed by stitching all the base images together. Denote as the discrete X-ray transform vector.

[0044] In the field of spectral CT image reconstruction, the core problem is based on the physical model , and through designing a reconstruction algorithm to accurately solve the base images from the projection data . In addition, in clinical applications, as elaborated above, the X-ray dose is often reduced by lowering the tube voltage or tube current, or reducing the sampling angle, etc., to ensure the safety of the scanned target object. However, the projection data obtained based on this has more noise or is sparser, which will not only lead to low reconstruction accuracy of the spectral CT image, but also significantly increase the complexity of modeling due to the sparse projection data, thereby resulting in low reconstruction efficiency of the spectral CT image.

[0045] On this basis, to solve the problems of low reconstruction accuracy and low reconstruction efficiency, especially to solve the problem of low reconstruction accuracy, relevant solutions introduce the prior information of the base images to be reconstructed, so that more accurate base images can be reconstructed under the guidance of the prior information. However, the original problem to be solved in the relevant reconstruction process is already a high-dimensional problem. On this basis, introducing prior information will further increase the dimension of this original problem, making it difficult to solve the original problem. That is, although the relevant solutions can improve the reconstruction accuracy to a certain extent, the improvement effect is not good, and it will further reduce the reconstruction efficiency.

[0046] In response to this, through careful research, it is found that based on certain characteristics of the original problem, the original problem can be decomposed into two or more independent sub-problems, especially sub-problems corresponding to each energy spectrum. Therefore, a reconstruction idea is proposed to decompose the original high-dimensional problem into low-dimensional sub-problems, introduce corresponding prior information for each sub-problem respectively, and perform reconstruction on this basis. Compared with the high-dimensional original problem, the low-dimensional sub-problems are easier to solve, thereby improving the solution efficiency and accuracy of the overall problem. Furthermore, each sub-problem is independent in structure and is naturally suitable for parallel computing, which can thus greatly shorten the reconstruction time and meet the clinical requirements for fast imaging. Next, this reconstruction idea will be elaborated in detail.

[0047] Before reconstruction, the known data can include projection data of energy spectra discrete X-ray transform matrices of energy spectra energy spectrum distributions and mass attenuation coefficients of the base materials used Among them,

[0048] can also be referred to as the projection matrix or the system matrix.

[0049] ,

[0050] and define the X-ray energy spectrum distribution matrix:

[0051] .

[0052] During the reconstruction process, what needs to be solved is the unknown number of basis images represents the number of energy spectra used, represents the number of base materials, and regardless of and what the situation is between them, the embodiments of the present invention can achieve the decomposition of the base materials.

[0053] Specifically, the embodiments of the present invention first utilize an important property of the spectral CT imaging operator : its Jacobian matrix has the following special form at certain special points :

[0054] ,

[0055] where , which represents a diagonal matrix composed of and is a constant matrix, related to and can be pre-computed. For example, it can be taken as . Specifically,

[0056] ,

[0057] where , can be defined as the points that satisfy the following conditions:

[0058] , for any ,

[0059] where is a set of real constants. In particular, when are all taken as 0, .

[0060] Utilizing this property, the embodiments of the present invention consider the following optimization problem (i.e., the original problem described above):

[0061] ,

[0062] where is the constraint set corresponding to the prior information. The physical meaning of this optimization problem is that in each iterative reconstruction process, by measuring the at in the approximate first-order Taylor expansion (i.e., ) and the error between the projection data , the basis images are gradually updated to approximate the true basis images. By minimizing the above formula, the reconstruction error can be effectively reduced and the image quality can be improved on the premise of satisfying .

[0063] However, since the dimension of the space where the basis images are located (i.e., ) is usually high, it is difficult to directly solve the above optimization problem. Utilizing the special form of the Jacobian matrix , through derivation, the objective function of the above optimization problem can be re-expressed as:

[0064] ,

[0065] Here, Since they are obtained by transformation based on the corresponding basis images, they can be called transformed basis images. Denote as the vector formed by splicing these transformed basis images, and they have the following relationship with the basis images :

[0066] ,

[0067] Among them, is the block matrix multiplication.

[0068] In addition, in Since it is obtained by transformation based on the projection data , it can be called the transformed projection data, which is defined as:

[0069] .

[0070] Furthermore, after defining and , the following low-dimensional sub-problems can be considered:

[0071] ,

[0072] Among them, is the corresponding constraint set. The above sub-problems are used to measure the and error between.

[0073] After updating , according to its relationship with the basis image , the following optimization problem can be considered to update the basis image :

[0074] ,

[0075] Among them, represents the constraint set corresponding to the prior information of the basis image , can represent a positive definite weight matrix. The physical meaning of this optimization problem is: by minimizing the and the transformed basis image error between, ensuring that the two are physically and mathematically consistent.

[0076] Based on the above derivation process, the spectral CT image reconstruction method described in the following embodiments is proposed.

[0077] Figure 1The flowchart of a spectral CT image reconstruction method provided by an embodiment of the present invention. This embodiment is applicable to the situation of spectral CT image reconstruction, especially applicable to the situation of spectral CT image reconstruction under low-dose conditions. This method can be executed by the spectral CT image reconstruction device provided by the embodiment of the present invention. The device can be implemented in software and / or hardware, and can be integrated on an electronic device, which can be various user terminals or servers.

[0078] See Figure 1 , the method of the embodiment of the present invention specifically includes the following steps:

[0079] S110. During the current iterative reconstruction process of the spectral CT image, for each of at least two spectra, determine the transformed projection data according to the projection matrix and the original projection data corresponding to the spectrum, and the previous basis image reconstructed during the previous iterative reconstruction process.

[0080] Among them, the spectral CT image is obtained through multiple iterative reconstructions. Taking the current iterative reconstruction as an example, during the current iterative reconstruction process, the transformed projection data and the transformed basis image corresponding to each of all spectra are determined in sequence. Taking the th spectrum as an example, the specific implementation is as follows:

[0081] Obtain the projection matrix corresponding to the th spectrum and the projection data . It should be noted that, in order to distinguish from the projection data obtained by subsequent transformation, the projection data [[ID= and the projection data are called the original projection data, and the projection data is called the transformed projection data. And obtain the basis image reconstructed during the previous iterative reconstruction process of the current iterative reconstruction process. Similarly, in order to distinguish from the basis image reconstructed during the current iterative reconstruction process, the basis image is called the previous basis image, and the basis image is called the current basis image. On this basis, according to the above description, the previous basis image can be considered as the basis image obtained by splicing the basis images corresponding to each basis material reconstructed during the previous iterative reconstruction process. Therefore, according to the previous basis image , the basis images corresponding to each basis material can be deduced inversely. Of course, if the current iterative reconstruction process is the first iterative reconstruction process, that is, there is no previous iterative reconstruction process, the preset initial basis image As the previous basis image Application

[0082] Furthermore, based on the obtained projection matrix and the original projection data , as well as the previous basis image , the transformed projection data corresponding to the th energy spectrum can be determined. That is, through the projection matrix and the previous basis image , the original projection data is transformed to obtain the transformed projection data . This transformed projection data can be understood as the projection data corresponding to the transformed basis image under the th energy spectrum to be solved in S120 .

[0083] Exemplarily, a determination scheme for the transformed projection data is given here:

[0084] Obtain the projection matrix, original projection data, and imaging operator corresponding to the energy spectrum; obtain the previous basis image reconstructed in the previous iterative reconstruction process; based on the imaging operator and the previous basis image, obtain the second projection data, and based on the projection matrix and the previous basis image, obtain the third projection data; based on the second projection data, original projection data, and third projection data, obtain the transformed projection data. Among them, after obtaining the projection matrix , original projection data , imaging operator , and previous basis image , the previous basis image can be processed based on the imaging operator to obtain the second projection data. For example, multiplying the two can obtain the second projection data ; and based on the projection matrix and the previous basis image , the third projection data can be obtained. For example, the basis images corresponding to each basis material can be obtained based on the previous basis image , and then based on the projection matrix and each basis image , the third projection data can be obtained. Furthermore, based on the second projection data, original projection data, and third projection data, the transformed projection data can be obtained. For example, the addition and subtraction result of these three projection data can be used as the transformed projection data .

[0085] On this basis, to better understand the above technical solution, the transformed projection data can be obtained through the following formula​ :

[0086] ,

[0087] wherein, represents the second projection data; represents the third projection data; is a constant in the constant matrix, and its specific meaning has been described above;

[0088] On this basis, the transformed projection data corresponding to each energy spectrum can be expressed by the following formula:

[0089] .

[0090] S120. Substitute the transformed projection data and the projection matrix into the sub-problem pre-constructed for the energy spectrum, and solve the substituted sub-problem based on the first constraint preset for the sub-problem, so as to obtain a transformed basis image that satisfies the first constraint and has the smallest error from the transformed projection data.

[0091] wherein, here continue to take the th energy spectrum as an example, and pre-construct a corresponding sub-problem for it. This sub-problem is used to measure the error between the transformed basis image after being projected by the projection matrix and the transformed projection data . By solving this sub-problem, a transformed basis image with the smallest error between the transformed projection data after being projected by the projection matrix and the transformed projection data can be obtained, that is, the obtained transformed basis image has the smallest error between the transformed projection data after being projected by the projection matrix and the transformed projection data . And for the sake of easy distinction, the transformed basis image with the smallest error here can be called the transformed basis image , that is, it is the optimal solution of the transformed basis image .

[0092] In an embodiment of the present invention, optionally, the above sub-problem can be expressed by minimizing a fourth formula, and the fourth formula can be obtained according to a fifth formula. The fifth formula can represent the error between a sixth formula and the transformed projection data , and the sixth formula can be expressed by the projection matrix and the transformed basis image to be solved. Further optionally, the sub-problems pre-constructed for each energy spectrum can be expressed by the following formula:

[0093] ,

[0094] Among them, represents the sixth formula, represents the fifth formula, represents the fourth formula. It can be known from this that each sub-problem is not related to each other and has a natural parallel solvable structure, thereby improving the reconstruction efficiency.

[0095] Furthermore, in order to improve the reconstruction accuracy, a first constraint can be preset for the sub-problems , and this first constraint can be used to constrain the transformed basis image in the sub-problem to ensure that the finally solved transformed basis image is consistent with the actual situation. In the embodiment of the present invention, optionally, the number of this first constraint can be one, two or more; further optionally, this first constraint can be at least one of non-negativity constraint and total variation (TV) constraint, etc. Among them, the non-negativity constraint can ensure the physical meaning of the transformed basis image , and the total variation constraint helps to suppress noise and enhance image edges, etc., further improving the reconstruction accuracy; the above content can be set according to actual needs and is not specifically limited here. On this basis, taking the first constraint including non-negativity constraint and total variation constraint as an example, the first constraint can be represented by the following formula:

[0096] ,

[0097] Among them, represents the discrete form of the gradient operator; is preset for the th energy spectrum; represents the set composed of all non-negative vectors.

[0098] After obtaining the sub-problems and the corresponding first constraints, the transformed projection data and the projection matrix can be substituted into the sub-problems, and the substituted sub-problems are solved based on the first constraint , so as to obtain a transformed basis image that satisfies the first constraint and has the smallest error with the transformed projection data after being projected by the projection matrix , that is, the transformed basis image . In an embodiment of the present invention, optionally, this solving process can be implemented based on any one of the original dual hybrid gradient algorithm, the alternating direction method of multipliers, and the semi-smooth Newton method, etc., which can be set according to actual needs and are not specifically limited herein.

[0099] On this basis, by way of example, combining the above various examples, here is a solving example of a sub-problem:

[0100] For the candidate transform basis images that satisfy the first constraint, calculate the first projection data according to the candidate transform basis images and the projection matrix, and determine the error between the first projection data and the transform projection data, so as to solve the candidate transform basis image corresponding to the smallest error among the errors corresponding to each candidate transform basis image, and use it as the transform basis image. Among them, the candidate transform basis images can be understood as the transform basis images that satisfy the first constraint of the transform basis images , and the number thereof can be one, two or more, which is related to the actual situation and is not specifically limited herein. For each candidate transform basis image, the first projection data can be calculated according to the candidate transform basis image and the projection matrix , for example, according to the above example, the product result of the candidate transform basis image and the projection matrix can be used as the first projection data; then, determine the error between the first projection data and the transform projection data , for example, the difference between the first projection data and the transform projection data can be used as the error, or the difference can be processed to obtain the error, etc. Thus, after obtaining the errors corresponding to each candidate transform basis image, the candidate transform basis image corresponding to the smallest error can be used as the optimal solution of the transform basis image , that is, used as the transform basis image , realizing the accurate solution of the transform basis image .

[0101] In addition, it should be noted that in spectral CT, the scanning geometric parameters of X-rays under different energy spectra often have inconsistent or mismatched situations, which will also affect the reconstruction accuracy of spectral CT images. In an embodiment of the present invention, corresponding sub-problems are constructed for each energy spectrum, that is, only one projection matrix is involved in one sub-problem , so as to separately process each projection matrix , thereby solving the problem of low reconstruction accuracy caused by inconsistent or mismatched scanning geometric parameters.

[0102] S130. Combine the transform basis images corresponding to each energy spectrum to obtain a combined basis image, and substitute the combined basis image into the pre-constructed optimization problem to reconstruct the current basis image with the smallest error from the combined basis image, realizing the reconstruction of spectral CT images.

[0103] Among them, after solving the transformation basis images corresponding to each energy spectrum these transformation basis images can be combined to obtain a combined basis image , that is, the combined basis image obtained through iterative processes. This combined basis image can be expressed as . .

[0104] The optimization problem can be understood as a problem of pre - constructing to measure the error between the combined basis image and the basis image to be solved , in order to solve the basis image to be solved with the smallest error between it and the combined basis image , that is, to solve the current basis image , and this current basis image is the optimal solution of the basis image to be solved . In the embodiments of the present invention, optionally, this optimization problem can be represented in various ways. Here, a representation example is given: this optimization problem can be represented by minimizing a first expression; this first expression is obtained by taking the weighted norm of a second expression through at least a weight matrix , such as directly obtaining it through the weighted norm or performing quadratic processing on the basis of the weighted norm. This weight matrix is obtained through presetting. For example, it can take the identity matrix or when Q = D, it can take etc., which can be set according to actual needs and are not specifically limited here. Among them, is a constant matrix; this second expression can represent the error between a third expression and the combined basis image , and this third expression is represented by a preset constant matrix and the basis image to be solved . Further optionally, this optimization problem can be represented by the following expression: where

[0105] ,

[0106] where represents the first expression, represents the second expression, represents the third expression, can be defined as , represents the transpose of the vector of the basis image of the th base material, represents block matrix multiplication, represents the basis image to be solved The second constraint to be satisfied, which will be described in detail later.

[0107] Further optionally, for the above constant matrix , assuming that the above spectral CT image reconstruction method involves spectral channels, energy nodes, and basis materials, where , and are all integers greater than 1, then the constant matrix can be obtained according to the X-ray spectral distribution matrix and the mass attenuation coefficient matrix , for example . Among them, the X-ray spectral distribution matrix may include ratios, which represent the ratio of the number of photons at the corresponding spectral channel and energy node to the total number of photons. Here, taking the ratio as an example, it represents the ratio of the number of photons at the -rd spectral channel and the -th energy node to the total number of photons. The mass attenuation coefficient matrix includes mass attenuation coefficients, which represent the mass attenuation coefficients of the corresponding basis materials at the corresponding energy nodes. Here, taking the mass attenuation coefficient as an example, it represents the mass attenuation coefficient of the -rd basis material at the -th energy node.

[0108] After obtaining the optimization problem, the combined basis image can be substituted into the optimization problem to reconstruct the current basis image with the smallest error between the combined basis image , especially the current basis image with the smallest error after being transformed by the constant matrix and the combined basis image , that is, the obtained current basis image has the smallest error with the combined basis image after being transformed by the constant matrix . From another perspective, the optimization problem can be solved to update the previous basis image to obtain the current basis image , thereby realizing the reconstruction of the spectral CT image. Exemplarily, after obtaining the current basis image , if the current iteration number has reached the maximum iteration number or the current basis image and the previous basis image the error between When the error is less than a preset error threshold, the current base image can be output , the spectral CT image reconstruction is completed; otherwise, iterative reconstruction can be performed again through the above steps to achieve spectral CT image reconstruction. On this basis, optionally, the above error can be represented by the difference between the current base image and the previous base image , or can also be represented by the formula , and so on, which can be set according to actual needs and are not specifically limited here.

[0109] In the technical solution of the embodiment of the present invention, during the current iterative reconstruction process of the spectral CT image, for each of at least two energy spectra, through the projection matrix corresponding to the energy spectrum and the original projection data , as well as the previous base image reconstructed in the previous iterative reconstruction process , the transformed projection data is determined; the transformed projection data and the projection matrix are substituted into the sub-problem pre-constructed for the energy spectrum, and based on the first constraint preset for the sub-problem, the substituted sub-problem is solved to obtain a transformed base image that satisfies the first constraint and has the smallest error between the projection after passing through the projection matrix and the transformed projection data ; then, the transformed base images corresponding to each energy spectrum are combined to obtain a combined base image , and the combined base image is substituted into the pre-constructed optimization problem, so as to reconstruct a current base image whose error between the transformation by the constant matrix and the combined base image is the smallest, realizing the reconstruction of the spectral CT image. The above technical solution, by defining the transformed base image , decomposes the high-dimensional original problem into low-dimensional sub-problems, and sets corresponding first constraints (i.e., prior information or regularization information) for each sub-problem, and thus, by using the first constraint On the basis of improving the reconstruction accuracy, the original problem can be decomposed to reduce the problem dimension, obtaining a low-dimensional, especially low-dimensional and parallel-solvable sub-problem, and the spectral CT image reconstruction is carried out by solving the sub-problem, thereby improving the reconstruction efficiency and further improving the reconstruction accuracy, achieving the effect of efficient and accurate reconstruction, especially achieving the effect of efficient and accurate reconstruction in the case of low dose.

[0110] An optional technical solution is to solve the substituted sub-problem based on the first constraint preset for the sub-problem, including:

[0111] Obtain the preset primal-dual hybrid gradient algorithm, and use the primal-dual hybrid gradient algorithm to solve the substituted sub-problem based on the first constraint preset for the sub-problem.

[0112] In this technical solution, the sub-problem can be solved based on the primal-dual hybrid gradient algorithm. Exemplarily, here the following th sub-problem (i.e., the th sub-problem) is taken as an example, and the solution steps are as follows:

[0113] ,

[0114] wherein, a step size parameter and an extrapolation parameter are given; a regularization parameter is given; let the inner iteration count index of the sub-problem , that is represents the number of iterations of the sub-problem, represents assignment; initialize , for example can be initialized according to its relationship with the base image , while can use the dual variable recorded after the previous iteration of the th sub-problem as the initial value, can use the dual variable recorded after the previous iteration of the th sub-problem as the initial value, which is used to solve the variables in the sub-problem and is continuously updated during the iteration.

[0115] It can be carried out in a cold start manner or a warm start manner. Among them, the core idea of the warm start manner is to use the results or relevant information of the previous inner iteration as the initial value for solving the current sub-problem, thereby accelerating convergence and improving computational efficiency; while the cold start manner is to completely not rely on any existing information when solving each sub-problem and re-initialize the variables.

[0116] Specifically, given the initial dual variable , the hot start initialization corresponds to:

[0117] .

[0118] The cold start initialization corresponds to:

[0119] .

[0120] Sub-step one: Update

[0121] ,

[0122] where, represents the orthogonal projection of the vector onto the set of all non-negative vectors . This is an update of the gradient descent type, is the gradient information from the data fidelity term, is the gradient information from the TV constraint term, passed through the dual variable. Finally, the updated vector is projected onto the non-negative set to ensure the non-negativity of the pixels in the image.

[0123] Sub-step two: Let

[0124] ,

[0125] where, this is an acceleration technique that uses the difference between the previous two iterations for "ahead" update, which helps to accelerate convergence.

[0126] Sub-step three: Update

[0127] ,

[0128] where, this is the update of the dual variable which is related to the data fidelity term " ". When differs greatly from , usually has a large adjustment. is the step size parameter.

[0129] Sub-step four: Update

[0130] ,

[0131] where, represents the projection of the vector onto the set . This is the update of the dual variable Update, which is related to the total variation constraint Related

[0132] Sub-step five: If Reaches the maximum number of inner iterations, then let:

[0133] ,

[0134] Otherwise, increment the iteration count by 1, i.e., , and continue the iteration, jumping to sub-step one.

[0135] The above technical solution realizes the efficient and accurate solution of the sub-problem through the primal-dual hybrid gradient algorithm.

[0136] Figure 2 is a flowchart of another energy spectrum CT image reconstruction method provided in an embodiment of the present invention. This embodiment is optimized based on the above technical solutions. In this embodiment, optionally, the combined basis image is substituted into the pre-constructed optimization problem to reconstruct the current basis image with the smallest error from the combined basis image, including: obtaining the pre-constructed optimization problem and the second constraint preset for the optimization problem; substituting the combined basis image into the optimization problem, and based on the second constraint, solving the substituted optimization problem to reconstruct the current basis image that satisfies the second constraint and has the smallest error from the combined basis image. Among them, the explanations of the same or corresponding terms as those in the above embodiments are not repeated here.

[0137] See Figure 2 , the method of this embodiment can specifically include the following steps:

[0138] S210. During the current iteration reconstruction process of the energy spectrum CT image, for each of at least two energy spectra, determine the transformed projection data according to the projection matrix and the original projection data corresponding to the energy spectrum, and the previous basis image reconstructed during the previous iteration reconstruction process.

[0139] S220. Substitute the transformed projection data and the projection matrix into the sub-problem pre-constructed for the energy spectrum, and solve the substituted sub-problem based on the first constraint preset for the sub-problem to obtain the transformed basis image that satisfies the first constraint and has the smallest error from the transformed projection data.

[0140] Among them, this step can obtain the transformed basis image that satisfies the first constraint and after being projected by the projection matrix has the smallest error from the transformed projection data .

[0141] S230. Combine the transformed basis images corresponding to each energy spectrum to obtain a combined basis image.​

[0142] S240. Substitute the combined basis image into the pre-constructed optimization problem, and solve the substituted optimization problem based on the second constraint preset for the optimization problem, so as to reconstruct the current basis image that satisfies the second constraint and has the minimum error from the combined basis image, thereby realizing the reconstruction of the spectral CT image.

[0143] Among them, in order to further improve the reconstruction accuracy, a second constraint can be preset for the optimization problem. This second constraint is used to constrain the basis image to be solved in the optimization problem so as to utilize this second constraint to solve the optimization problem, thereby obtaining the current basis image that satisfies this second constraint and has the minimum error from the combined basis image after being transformed by the constant matrix . .

[0144] In an embodiment of the present invention, optionally, the above-mentioned second constraint may be at least one of those in the real vector space, non-negativity constraint, total variation constraint, etc., which can be set according to actual needs and are not specifically limited herein.

[0145] The technical solution of the embodiment of the present invention, by presetting a second constraint for the optimization problem in advance , and then using this second constraint to solve the optimization problem, which helps to further improve the reconstruction accuracy of the spectral CT image.

[0146] An optional technical solution is that the optimization problem is represented by minimizing a first expression;

[0147] The first expression is obtained by at least taking the weighted norm of a second expression through a weight matrix;

[0148] The second expression represents the error between a third expression and the combined basis image;

[0149] The third expression is represented by a preset constant matrix and the basis image to be solved, where the current basis image is the basis image obtained by solving the basis image to be solved.

[0150] The above representation scheme of the optimization problem has been described above and will not be elaborated here.

[0151] On this basis, several optional selection methods for the second constraint and the weight matrix are given, as well as the solution method for the optimization problem, especially when the constant matrix Given in the reversible case, whereby adjustments can be made and , flexibly determine the specific update method of the base image, and through adjust the importance of different base materials or spatial regions, so as to better adapt to different application scenarios and data characteristics.

[0152] Optionally, when the second constraint is located in the real vector space , , and the weight matrix is a positive definite matrix, the current base image can be obtained according to the inverse matrix of the constant matrix and the combined base image .

[0153] Exemplarily (i.e., Example 1), when is taken as , , and the constant matrix has full column rank, regardless of what kind of positive definite matrix it is, the current base image is all:

[0154] ,

[0155] At this time, no prior constraint is imposed on the base image to be solved , then the current base image can be directly obtained by linearly combining the transformed base images corresponding to different energy spectra in the combined base image . This represents an idealized and direct material decomposition process. At this time, the optimization problem can be simplified to the unconstrained linear least squares problem shown by the following formula, and the above linear combination can be considered as the standard solution to solve this unconstrained linear least squares problem, and the weight matrix W does not change the solution:

[0156] .

[0157] Optionally again, when the second constraint is a non - negative constraint , , the constant matrix is invertible, and the weight matrix is the product of the above - mentioned inverse matrix and the transpose matrix of the inverse matrix , the current base image can be obtained by non - negative projection of the estimated base image, and the estimated base image is based on the inverse matrix and the combined base image Obtained.

[0158] Exemplary (i.e., Example 2), when is taken as , , and the constant matrix is invertible, i.e.,

[0159] ,

[0160] is taken as the following closed - form solution exists for the above - mentioned optimization problem:

[0161] ,

[0162] where represents the estimated basis image. The above Example 2 can be understood as first performing a direct material decomposition with Example 1 to obtain a preliminary estimated basis image ; then, applying non - negative projection to this estimated basis image to force all calculated negative pixel values to zero. This ensures that the reconstructed current basis image is physically realistic (i.e., non - negative). The above Example 2 makes the optimization problem equivalent to:

[0163]

[0164] ,

[0164] The optimal solution to this optimization problem is to orthogonally project the unconstrained solution onto the feasible set This can be considered as a computationally efficient way to introduce non - negative constraints.

[0165] On this basis, the non - negative constraint described above is further explained:

[0166] ,

[0167] This constraint reflects the non - negative property of image pixel values and has a clear physical meaning. That is, in a clinical scenario, the pixel values of the basis image to be solved physically correspond to non - negative quantities. Therefore, by introducing the non - negative constraint not only can it better conform to the physical characteristics of the actual problem, but also effectively narrow the solution space range, thereby improving the solution efficiency and reconstruction accuracy.

[0168] In addition, when in order to ensure the existence and uniqueness of the solution to the optimization problem, prior information needs to be additionally introduced for the basis image to be solved. For example, when When, it is possible to take as

[0169] ,

[0170] This is a physical assumption that conforms to reality.

[0171] Alternatively, when the second constraint includes a non - negative constraint and a total variation constraint, and the weight matrix is the identity matrix, the current basis image can be solved based on the preset primal - dual hybrid gradient algorithm.

[0172] Exemplarily (i.e., Example 3), when and can have any size relationship and is taken as

[0173] ,

[0174] is taken as the identity matrix, the above - mentioned optimization problem has no closed - form solution. Here, it can be solved based on the primal - dual hybrid gradient algorithm. Of course, it can also be solved based on the alternating direction method of multipliers or the semi - smooth Newton method, etc., which is not specifically limited here. Here, a set of non - negative basis images are sought. These basis images should not only correspond well to the transformed basis image in terms of the data fidelity term:

[0175] ,

[0176] but each basis image itself should have the piece - wise constant property (i.e., satisfy the total variation constraint) to suppress noise and protect edges.

[0177] On this basis, the above - mentioned total variation constraint is further explained:

[0178] ,

[0179] wherein, represents the discrete - form gradient operator, represents the 1 - norm of the vector, represents the regularization parameter, which is preset for each energy spectrum and controls the strength of this total variation constraint. The above - mentioned total variation constraint corresponds to the piece - wise constant property of the basis image to be solved, that is, in most regions of the basis image to be solved, it shows piece - wise constant changes, and large gradient changes are allowed at the boundaries or edges. In practical applications, the total variation constraint can effectively characterize the basis image The spatial distribution characteristics are particularly suitable for describing images with sparse gradients (such as tissue structures in medical images). By restricting the 1-norm of the gradient, the total variation constraint can suppress the influence of noise on the reconstruction result while retaining the important edge information in the image.

[0180] Furthermore, the implementation process of solving the optimization problem based on the primal-dual hybrid gradient algorithm described above is further explained:

[0181] Given the step size parameter , and the extrapolation parameter . Initialize the dual variable . When solving the optimization problem here, the two initialization methods of warm start or cold start described above can also be referred to, and will not be elaborated again. Let . Sub-step one: Update :

[0182] ,

[0183] where, under the condition of the current dual variable given, update the base image by gradient descent and immediately project it onto the non-negative set.

[0184] Sub-step two: Let

[0185] .

[0186] Sub-step three: Update

[0187] ,

[0188] where, according to the residual between the extrapolated value of the current base image and the transformed base image , adjust the dual variable .

[0189] Sub-step four: Update

[0190] ,

[0191] where, according to the extrapolated value of the gradient of the current base image, adjust the dual variable related to the total variation constraint.

[0192] Sub-step five: If the maximum number of inner iterations is reached, then output ; otherwise, increment the iteration count by 1, that is , continue the iteration, and jump to sub-step one.

[0193] The above example realizes the efficient and accurate solution of the optimization problem.

[0194] Figure 3 It is a flowchart of another energy spectrum CT image reconstruction method provided in an embodiment of the present invention. This embodiment is optimized based on the above technical solutions. Among them, the explanations of the same or corresponding terms as those in the above embodiments will not be repeated here.

[0195] See Figure 3 , the method of this embodiment may specifically include the following steps:

[0196] S310. Obtain the initial basis image , and the original projection data corresponding to each energy spectrum , the X-ray energy spectrum distribution matrix and the mass attenuation coefficient matrix , and set the constant matrix and the current iteration number .

[0197] S320. Calculate the transformed projection data corresponding to each energy spectrum , that is, obtain .

[0198] Among them, is to construct the projection data corresponding to for the sub-problem in the following S330.

[0199] S330. Solve sub-problems in parallel to solve the predefined transformed basis image , and after each sub-problem is solved, set to obtain the transformed basis image .

[0200] S340. Update the previous basis image by solving the optimization problem to obtain the current basis image .

[0201] S350. If the current iteration number has reached the maximum iteration number or the error between the current basis image and the previous basis image is less than the preset error threshold, then output the current basis image ; otherwise, increment the current iteration number by 1, that is , and jump to S320 to continue the iteration.

[0202] The technical solution of the embodiment of the present invention achieves the effect of accurately and efficiently reconstructing low-dose spectral CT images.

[0203] Verified by experiments, for the above technical solution, under the condition of sparse and noise-free data, the method described in the embodiment of the present invention can quickly converge to a result that is highly close to the real image (see details in Figure 4a , Figure 4b , Figure 6a and Figure 6b ), and there is no visual difference (see details in Figure 5 and Figure 7 ); in terms of numerical indicators such as relative error, it is significantly better than the results of related solutions (i.e., using the non-convex primal-dual algorithm for spectral CT image reconstruction) (see Table 1 and Table 2 for details). In particular, the relative error between the image reconstructed by the method described in the embodiment of the present invention and the real image has approached the machine precision.

[0204] In addition, under the condition of sparse and high-noise data, the method described in the embodiment of the present invention can quickly converge to a result that is relatively close to the real image (see details in Figure 8a , Figure 8b , Figure 9a , Figure 9b , Figure 10a and Figure 10b ); the obtained image is highly consistent with the real image in details (see details in Figure 11 ), and has strong robustness to noise and is not easily affected by noise interference; in terms of numerical indicators such as relative error, it is significantly better than the results of related solutions (see Table 3 for details).

[0205] In addition, the method described in the embodiment of the present invention can be considered as an algorithm framework, and its specific implementation modules can flexibly adopt a variety of design schemes, with a wide range of applications and good scalability. For example, the second constraint and the weight matrix can be selected in different forms according to requirements; the number of inner iterations of the sub-problem can also be flexibly adjusted. When the number of inner iterations is large, different sub-problems can be solved simultaneously in a parallel manner, which can improve the overall operation efficiency of the algorithm.

[0206] Here, through the dual-energy CT simulation experiment under sparse projection angles, the advantages of the method described in the embodiments of the present invention are verified. For a fair comparison with related solutions, considering that the method described in the embodiments of the present invention needs to solve sub-problems through inner iterations, while the related solutions do not set inner iterations, the number of inner iterations for the sub-problems in the method described in the embodiments of the present invention can be set to 1 time. In fact, the method described in the embodiments of the present invention allows flexible setting of different numbers of inner iterations, and can effectively reconstruct the basis images under different settings. This will not be further experimentally demonstrated. It should be noted that in the method described in the embodiments of the present invention, multiple independent sub-problems have a natural parallel solution structure. When the number of inner iterations is large, parallel computing can be used for acceleration to further improve the computing efficiency. In addition, the running time is related to the number of iterations, while the memory occupancy is not affected by the number of iterations.

[0207] In addition, there are some hyperparameters in the method described in the embodiments of the present invention. Taking the total variation regularization parameter as an example, its physical meaning is the total variation value of the basis image to be solved. In practice, according to the target imaging task (such as organ parts like the head or chest) and the type of imaging device, statistical analysis can be performed by referring to similar image datasets, and a reasonable parameter range can be empirically set based on this.

[0208] Here, the relative error of the basis image and the relative error of the projection data are defined:

[0209] ,

[0210] where, is the true image; is the basis image determined after k iterations, such as the basis images and the virtual monoenergetic images at the corresponding energies; represents the 2-norm of the vector. When reconstructing noisy data, to test the numerical convergence of the method described in the embodiments of the present invention, the relative error of the basis images in two adjacent steps (i.e., the residual described above) and the relative error of the projection data in two adjacent steps can be considered, and are respectively defined as:

[0211] .

[0212] In addition, the following quantitative index can also be used to measure the quality of the reconstructed image. Given two image vectors and , their normalized root mean-squared error (NRMSE) is defined as:

[0213] ,

[0214] And the structural similarity (SSIM) is defined as:

[0215] ,

[0216] where corresponds to the th window, corresponds to the th window, represents the number of windows, then the local similarity index is defined as:

[0217] ,

[0218] where represents the mean of , represents the variance of , represents the mean of , represents the variance of , represents the covariance between and . In addition . The peak signal-to-noise ratio (PSNR) is defined as:

[0219] .

[0220] On this basis, the following Experiments 1, 2, and 3 are carried out (the following experiments are carried out taking as an example):

[0221] Experiment 1: Reconstructing the head phantom using sparse noise-free data

[0222] Experiment 1 uses the head phantom as the real image. The head phantom is composed of water-based materials and bone-based materials, distributed in a region of [-15 cm, 15 cm] × [-15 cm, 15 cm], and the image size is 256 × 256, as Figure 5As shown in the first line. In the experiment, a low-energy spectrum of 80 kVp and a high-energy spectrum of 140 kVp were used. Among them, the X-rays of the high-energy spectrum were processed by a 1-mm-thick copper filter. Specifically, for the low-energy spectrum, 30 projection angles evenly distributed between 0 degrees and 180 degrees were used, and 301 equally spaced X-ray beams within the detection range of [-15 cm, 15 cm] were set at each projection angle. The high-energy spectrum maintained the same number of rays, but its 30 projection angles were evenly distributed between 3 degrees and 183 degrees. The energy spectrum distribution, mass attenuation coefficient data, and real image were substituted into the imaging model to generate noiseless dual-energy projection data.

[0223] In this experiment, the initial basis image was set to zero, and all dual variables were also initialized to zero, that is

[0224] ,

[0225] Take the identity matrix, and the regularization parameter was taken as the total variation value of the real transformation basis image, that is

[0226] ,

[0227] The step size parameter was taken as 0.08 for both, and the extrapolation parameter was taken as 1.0. In addition, in the non-convex primal-dual algorithm in the related scheme, the initial basis image was set to zero, all dual variables were also initialized to zero, and the parameter of the total variation regularization term was taken as the total variation value of the real image, that is , and the step size parameter was also taken as 0.08 for both, and the extrapolation parameter was also taken as 1.0.

[0228] In the experiment, both the method described in the embodiment of the present invention and the non-convex primal-dual algorithm in the related scheme were iterated 200,000 times, and the curves of the relative error of the basis image and the relative error of the projection data changing with the number of iterations were plotted. The two curves are as shown in Figure 4a and Figure 4b . It can be seen from the curves of the relative error changing with the number of iterations that in the iterative reconstruction process of the method described in the embodiment of the present invention, both the relative error of the basis image and the relative error of the projection data decreased significantly with the increase of the number of iterations; in contrast, although the relative error of the basis image and the relative error of the projection data of the related scheme also decreased with the increase of the number of iterations, the rate of error decrease was much lower than that of the method described in the embodiment of the present invention. Finally, the relative error level achieved by the method described in the embodiment of the present invention was significantly lower than that of the related scheme. This fully reflects the advantages of the method described in the embodiment of the present invention in terms of speed and accuracy.

[0229] Furthermore, we presented the reconstruction results of the method described in the embodiment of the present invention and the related scheme in Figure 5Among them, the first row is the real image, the second row is the reconstruction result of the method according to the embodiment of the present invention, and the third row is the reconstruction result of the related solution; from left to right are the water-based image, the bone-based image, the 60 keV virtual monoenergetic image, and the 100 keV virtual monoenergetic image. From Figure 5 It can be seen that the image reconstructed by the method according to the embodiment of the present invention is almost completely consistent with the real image. In the reconstruction result of the related solution, there are obvious strip artifacts caused by undersampling in the base image, and the virtual monoenergetic image is blurred, and the overall image quality is poor. This comparison further verifies the advantages of the method according to the embodiment of the present invention in terms of reconstruction accuracy and image quality.

[0230] Table 1 shows the display results of the evaluation indexes provided by the method according to the embodiment of the present invention. For each cell in Table 1, from top to bottom are the normalized root mean square error, the structural similarity, and the peak signal-to-noise ratio. It can be seen from Table 1 that under the same experimental conditions, for the dual-energy CT image reconstruction problem of the head phantom with the same scale of sparse noise-free, the method according to the embodiment of the present invention is significantly better than the non-convex primal-dual algorithm in the related solution in terms of various numerical indexes.

[0231] Table 1 Quantitative evaluation of image reconstruction quality (head phantom)

[0232]

[0233] Experiment 2: Reconstruct the XCAT phantom using sparse noise-free data

[0234] In Experiment 2, the XCAT phantom was used to construct a digital phantom as the real image. The phantom was composed of water-based materials and bone-based materials, distributed in a region of [-15 cm, 15 cm] × [-15 cm, 15 cm], and the image size was 256×256, as Figure 7 shown in the first row. The imaging system used a low energy spectrum of 80 kVp and a high energy spectrum of 140 kVp, and the X-rays of the high energy spectrum were processed by a 1-mm-thick copper filter. Specifically, the low energy spectrum used 60 projection angles evenly distributed from 0 degrees to 180 degrees, and each projection angle was set with 301 equally spaced X-ray beams distributed in the detection range of [-15 cm, 15 cm]. The high energy spectrum maintained the same number of rays, but the 60 projection angles were evenly distributed from 1.5 degrees to 181.5 degrees. Substituting the energy spectrum distribution, the mass attenuation coefficient data, and the real image into the imaging model can generate noise-free dual-energy projection data.

[0235] In this experiment, the initial base image was set to zero, and all dual variables were also initialized to zero, that is

[0236] ,

[0237] Take the identity matrix, and the regularization parameter All are taken as the total variation values of the true transformation basis images, i.e.,

[0238] ,

[0239] The step size parameter is all taken as 0.08, and the extrapolation parameter is all taken as 1.0. In addition, in the non-convex primal-dual algorithm in the related scheme, the initial basis image is set to zero, all dual variables are also initialized to zero, and the parameter of the total variation regularization term is taken as the total variation value of the true image, i.e., . The step size parameters in the related scheme are all taken as 0.1, and the extrapolation parameter is taken as 1.0.

[0240] In the experiment, both the method described in the embodiment of the present invention and the non-convex primal-dual algorithm (including the total variation regularization term) in the related scheme are iterated 150,000 times, and the curve graphs of the relative error of the basis image and the relative error of the projection data changing with the number of iterations are plotted. The two curve graphs are as Figure 6a and Figure 6b shown.

[0241] It can be seen from the curve graph of the relative error changing with the number of iterations that in the iterative reconstruction process of the method described in the embodiment of the present invention, both the relative error of the basis image and the relative error of the projection data decrease with the increase of the number of iterations; while in the iterative reconstruction process of the related scheme, although its relative error of the basis image and the relative error of the projection data also decrease with the increase of the number of iterations, the error reduction speed is much smaller than that of the method described in the embodiment of the present invention. The relative error level of the method described in the embodiment of the present invention is significantly lower than that of the related scheme, which shows the advantages of the method described in the embodiment of the present invention in terms of speed and accuracy.

[0242] We display the results reconstructed by the method described in the embodiment of the present invention and the related scheme in Figure 7 , where the first row represents the true image, the second represents the reconstruction result of the method described in the embodiment of the present invention, and the third row represents the reconstruction result of the related scheme; from left to right are the water basis image, the bone basis image, the 60 keV virtual monoenergetic image, and the 100 keV virtual monoenergetic image. From Figure 7 it can be seen that the image reconstructed by the method described in the embodiment of the present invention has almost no difference from the true image, while there are obvious artifacts in the basis image reconstructed by the related scheme.

[0243] Table 2 shows the evaluation index display results provided by the method according to the embodiments of the present invention. For each cell in Table 2, from top to bottom are the normalized root mean square error, structural similarity, and peak signal-to-noise ratio. As can be seen from Table 2, under the same experimental conditions, for the dual-energy CT image reconstruction problem of the XCAT phantom with the same scale and sparse noise-free, the method according to the embodiments of the present invention is significantly superior to the non-convex primal-dual algorithm in the related solutions in terms of various numerical indexes.

[0244] Table 2 Quantitative evaluation of image reconstruction quality (XCAT phantom)

[0245]

[0246] Experiment 3: Reconstructing real patient lung CT images using sparse noisy data

[0247] In Experiment 3, real patient lung CT images in the LIDC-IDRI (The Lung Image Database Consortium and Image Database Resource Initiative) dataset are used as the real images. This image is composed of water-based materials and bone-based materials, distributed in a region of [-15 cm, 15 cm] × [-15 cm, 15 cm], and the image size is 256×256, as Figure 11 shown in the first row. The imaging system uses a low energy spectrum of 80 kVp and a high energy spectrum of 140 kVp, where the X-rays of the high energy spectrum are processed by a 1-mm thick copper filter. Specifically, the low energy spectrum uses 120 projection angles evenly distributed from 0 degrees to 180 degrees, and each projection angle is set with 301 equally spaced X-ray beams distributed in the detection range of [-15 cm, 15 cm]. The high energy spectrum maintains the same number of rays, but the 120 projection angles are evenly distributed from 0.75 degrees to 180.75 degrees. Substituting the energy spectrum distribution, mass attenuation coefficient data, and real image into the imaging model can generate noise-free dual-energy projection data. Subsequently, noisy dual-energy projection data is generated by adding Gaussian white noise. In this experiment, the signal-to-noise ratio of the dual-energy projection data is 23.0 dB.

[0248] In this experiment, the initial basis image is set to zero, and all dual variables are also initialized to zero, that is

[0249] .

[0250] In step 4, the weight matrix is taken as the identity matrix, and the second constraint is taken as

[0251] ,

[0252] where The step size parameter is taken as 0.04 for all, and the extrapolation parameter is taken as 1.0, and the base image is updated using 50 iterations. In addition, in the non-convex primal-dual algorithm in the related scheme, the initial base image is set to zero, all dual variables are also initialized to zero, and the parameter of the total variation regularization term is taken as the total variation value of the real image, that is . The step size parameters in the related scheme are also taken as 0.04 for all, and the extrapolation parameter is also taken as 1.0.

[0253] In this experiment, the method described in the embodiments of the present invention and the non-convex primal-dual algorithm in the related scheme are both iterated 100,000 times, and the curves of the relative error of the base image and the relative error of the projection data varying with the number of iterations are plotted, as can be seen in Figure 8a and Figure 8b . It can be seen from the curves of the relative error varying with the number of iterations that in the iterative reconstruction process of the method described in the embodiments of the present invention, both the relative error of the base image and the relative error of the projection data gradually decrease with the increase of the number of iterations and converge rapidly; in contrast, although the related scheme shows a trend that the relative error of the base image and the relative error of the projection data gradually decrease with time in the initial stage of iteration, after a certain number of iterations, the relative error of its base image starts to rise instead. In addition, the relative error level of the method described in the embodiments of the present invention is significantly lower than that of the related scheme.

[0254] In addition, the curves of the relative error of the base image and the relative error of the projection data varying with time are also plotted, as can be seen in Figure 9a and Figure 9b . It can be seen from the figure that within the shown time range, both the relative error of the base image and the relative error of the projection data of the method described in the embodiments of the present invention gradually decrease with the running time and quickly converge to a lower level. Although the relative error of the related scheme also gradually decreases with the running time, its decreasing speed is slower. In addition, the relative error level achieved by the method described in the embodiments of the present invention is significantly lower than that of the related scheme. These results fully demonstrate that the method described in the embodiments of the present invention has obvious advantages in the convergence speed and reconstruction accuracy, can achieve high-quality image reconstruction in a shorter time, effectively suppress noise at the same time, and shows higher computational efficiency and robustness.

[0255] Figure 10a shows the change of the relative error of the base image between adjacent two steps, while Figure 10b shows the change of the relative error of the projection data between adjacent two steps. It can be seen from the figure that the error between adjacent two steps of the method described in the embodiments of the present invention has approached the machine precision. In contrast, the non-convex primal-dual algorithm still needs more iterations to converge.

[0256] Figure 11The figure shows a comparison chart of the reconstruction results obtained by the method according to the embodiments of the present invention and related solutions. Among them, the first row is the real image, the second row is the reconstruction result of the method according to the embodiments of the present invention, and the third row is the reconstruction result of the related solution; from left to right are the water-based image, bone-based image, 60 keV virtual monoenergetic image, and 100 keV virtual monoenergetic image. It can be clearly seen from the figure that there is almost no difference between the image reconstructed by the method according to the embodiments of the present invention and the real image, while there are obvious artifacts caused by noise in the base images reconstructed by the related solution.

[0257] In addition, Figure 11 an enlarged comparison chart of the local details of the real image and the reconstruction result is also provided. From these details, it can be further observed that the water-based image, bone-based image, 60 keV virtual monoenergetic image, and 100 keV virtual monoenergetic image reconstructed by the method according to the embodiments of the present invention are highly consistent with the real image in details, with extremely small differences; while there are significant noise interferences in the water-based image, 60 keV virtual monoenergetic image, and 100 keV virtual monoenergetic image reconstructed by the related solution, resulting in blurred details of some tissues, excessive noise in the image, and it is difficult to distinguish fine structures. These results further verify the advantages of the method according to the embodiments of the present invention in terms of image quality and anti-noise performance.

[0258] Table 3 shows the display results of the evaluation indexes provided by the method according to the embodiments of the present invention. For each cell in Table 3, from top to bottom are the normalized root mean square error, structural similarity, and peak signal-to-noise ratio. It can be seen from Table 3 that under the same experimental conditions, for the dual-energy CT image reconstruction problem with sparse and noisy images of the same scale, the method according to the embodiments of the present invention is significantly better than the non-convex primal-dual algorithm in the related solution in terms of various numerical indexes.

[0259] Table 3 Quantitative evaluation of image reconstruction quality (real patient lung CT image)

[0260]

[0261] The above 3 experiments all prove that the method according to the embodiments of the present invention can efficiently and accurately reconstruct the spectral CT image.

[0262] Figure 12 is a structural block diagram of the spectral CT image reconstruction device provided by the embodiments of the present invention. This device is used to execute the spectral CT image reconstruction method provided by any of the above embodiments. This device and the spectral CT image reconstruction methods of the above embodiments belong to the same inventive concept. For the details not described in detail in the embodiments of the spectral CT image reconstruction device, reference can be made to the embodiments of the above spectral CT image reconstruction method. See Figure 12 , specifically, this device may include: a transformed projection data determination module 410, a transformed base image obtaining module 420, and a spectral CT image reconstruction module 430.

[0263] Among them, the transformed projection data determination module 410 is configured to determine transformed projection data for each of at least two energy spectra during the current iteration reconstruction process of the spectral CT image, according to the projection matrix corresponding to the energy spectrum, the original projection data, and the previous basis image reconstructed during the previous iteration reconstruction process.

[0264] The transformed basis image obtaining module 420 is configured to substitute the transformed projection data and the projection matrix into a sub-problem pre-constructed for the energy spectrum, and solve the substituted sub-problem based on a first constraint preset for the sub-problem, so as to obtain a transformed basis image that satisfies the first constraint and has the minimum error from the transformed projection data.

[0265] The spectral CT image reconstruction module 430 is configured to combine the transformed basis images corresponding to each energy spectrum to obtain a combined basis image, and substitute the combined basis image into a pre-constructed optimization problem, so as to reconstruct a current basis image that has the minimum error from the combined basis image, thereby realizing the reconstruction of the spectral CT image.

[0266] Optionally, the first constraint is used to constrain the transformed basis image to be solved by the sub-problem, and the sub-problem is solved by the following unit:

[0267] The transformed basis image obtaining unit is configured to calculate first projection data according to a candidate transformed basis image that satisfies the first constraint and the projection matrix, and determine the error between the first projection data and the transformed projection data, so as to solve for the candidate transformed basis image corresponding to the minimum error among the errors corresponding to each candidate transformed basis image, and use it as the transformed basis image.

[0268] Optionally, the transformed basis image obtaining module 420 may include:

[0269] The sub-problem solving unit is configured to obtain a preset primal-dual hybrid gradient algorithm, and use the primal-dual hybrid gradient algorithm to solve the substituted sub-problem based on the first constraint preset for the sub-problem.

[0270] Optionally, the first constraint includes a non-negativity constraint and / or a total variation constraint.

[0271] Optionally, the spectral CT image reconstruction module 430 may include:

[0272] The second constraint obtaining unit is configured to obtain a pre-constructed optimization problem and a second constraint preset for the optimization problem;

[0273] The current basis image reconstruction unit is configured to substitute the combined basis image into the optimization problem, and solve the substituted optimization problem based on the second constraint, so as to reconstruct a current basis image that satisfies the second constraint and has the minimum error from the combined basis image.

[0274] Based on any of the above devices, optionally, the optimization problem is represented by minimizing a first expression;

[0275] The first expression is obtained by at least taking the weighted norm of a second expression through a weight matrix;

[0276] The second expression represents the error between a third expression and a combined basis image;

[0277] The third expression is represented by a preset constant matrix and a basis image to be solved, where the current basis image is the basis image obtained by solving the basis image to be solved.

[0278] On this basis, optionally, the above device involves an energy spectrum and a basis material, where and are both integers greater than 1. If the current basis image is reconstructed based on the optimization problem and a second constraint preset for the optimization problem, then:

[0279] When the second constraint is in the real vector space, , the constant matrix is invertible, and the weight matrix is a positive definite matrix, the current basis image is obtained according to the inverse matrix of the constant matrix and the combined basis image; and / or,

[0280] When the second constraint is a non-negativity constraint, , the constant matrix is invertible, and the weight matrix is the product of the inverse matrix and the transpose matrix of the inverse matrix, the current basis image is obtained by performing non-negative projection on an estimated basis image, and the estimated basis image is obtained according to the inverse matrix and the combined basis image; and / or,

[0281] When the second constraint includes a non-negativity constraint and a total variation constraint, and the weight matrix is an identity matrix, the current basis image is obtained by solving based on a preset primal-dual hybrid gradient algorithm.

[0282] Another option is that the above device involves an energy spectrum, energy nodes, and a basis material, where , and are all integers greater than 1; the constant matrix is obtained according to the X-ray energy spectrum distribution matrix and the mass attenuation coefficient matrix;

[0283] Among them, the X-ray energy spectrum distribution matrix includes proportions, and the proportion represents the proportion of the number of photons at the corresponding energy spectrum and energy node to the total number of photons;

[0284] The mass attenuation coefficient matrix includes a mass attenuation coefficient, where the mass attenuation coefficient represents the mass attenuation coefficient of the corresponding base material at the corresponding energy node.

[0285] Optionally, the transformed projection data determination module 410 may include:

[0286] an imaging operator acquisition unit, configured to acquire a projection matrix corresponding to the energy spectrum, original projection data, and an imaging operator;

[0287] a previous base image acquisition unit, configured to acquire the previous base image reconstructed in the previous iterative reconstruction process;

[0288] a third projection data obtaining unit, configured to obtain second projection data according to the imaging operator and the previous base image, and obtain third projection data according to the projection matrix and the previous base image;

[0289] a transformed projection data obtaining unit, configured to obtain transformed projection data according to the second projection data, the original projection data, and the third projection data.

[0290] The spectral CT image reconstruction device provided by the embodiments of the present invention can achieve efficient and accurate spectral CT image reconstruction through the mutual cooperation of each module, especially in the case of low dose, efficient and accurate spectral CT image reconstruction can be achieved.

[0291] The spectral CT image reconstruction device provided by the embodiments of the present invention can execute the spectral CT image reconstruction method provided by any embodiment of the present invention, and has the corresponding functional modules and beneficial effects for executing the method.

[0292] It should be noted that in the embodiments of the above spectral CT image reconstruction device, the included units and modules are only divided according to functional logic, but are not limited to the above division, as long as the corresponding functions can be achieved; in addition, the specific names of the functional units are only for the convenience of mutual distinction and do not limit the protection scope of the present invention.

[0293] Figure 13 FIG. shows a schematic structural diagram of an electronic device 10 that can be used to implement the embodiments of the present invention. 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 only examples and are not intended to limit the implementation of the present invention described and / or claimed herein.

[0294] As shown Figure 13 in FIG. 1, the electronic device 10 includes at least one processor 11 and a memory communicatively connected to the at least one processor 11, such as a read-only memory (ROM) 12, a random access memory (RAM) 13, etc. The memory stores a computer program executable by the at least one processor. 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 into 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 via a bus 14. An input / output (I / O) interface 15 is also connected to the bus 14.

[0295] Multiple 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 magnetic 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 via a computer network such as the Internet and / or various telecommunication networks.

[0296] The processor 11 can be various general-purpose and / or special-purpose 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 dedicated 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 performs the various methods and processes described above, such as the spectral CT image reconstruction method.

[0297] In some embodiments, the spectral CT image reconstruction method can be implemented as a computer program tangibly embodied in a computer-readable storage medium, such as the storage unit 18. In some embodiments, part or all of the computer program can be loaded and / or installed onto 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 can be performed. Alternatively, in other embodiments, the processor 11 can be configured to perform the spectral CT image reconstruction method by any other appropriate means (e.g., by means of firmware).

[0298] The various embodiments of the systems and techniques described above in this specification can be implemented in digital electronic circuitry, integrated circuit systems, field programmable gate arrays (FPGA), application specific integrated circuits (ASIC), application specific standard products (ASSP), systems on a chip or system-on-chip (SOC), complex programmable logic devices (CPLD), computer hardware, firmware, software, and / or combinations thereof. These various embodiments can include: being implemented in one or more computer programs that are executable and / or interpretable on a programmable system including at least one programmable processor, which can be a special-purpose or general-purpose programmable processor that receives data and instructions from, and transmits data and instructions to, a storage system, at least one input device, and at least one output device.

[0299] The computer programs for implementing the methods of the present invention can be written in any combination of one or more programming languages. These computer programs can be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus, such that the computer programs, when executed by the processor, cause the functions / operations specified in the flowchart and / or block diagram to be implemented. The computer programs can be executed entirely on the machine, partly on the machine, as a stand-alone software package partly on the machine and partly on a remote machine or entirely on the remote machine or server.

[0300] In the context of the present invention, a computer-readable storage medium can be a tangible medium that can contain or store a computer program for use by or in connection with an instruction execution system, apparatus, or device. The computer-readable storage medium can include, but is not limited to, electronic, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatus, or devices, or any suitable combination of the foregoing. Alternatively, the computer-readable storage medium can be a machine-readable signal medium. More specific examples of the machine-readable storage medium would include an electrical connection based on one or more wires, a portable computer diskette, 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 disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing.

[0301] To provide interaction with a user, the systems and techniques described herein can 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 a trackball) through which the user can provide input to the electronic device. Other kinds of devices can also be used to provide interaction with the user; for example, the feedback provided to the user can be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback); and input from the user can be received in any form (including acoustic input, voice input, or tactile input).

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

[0303] The computing system can include a client and a server. The client and the server are generally far from each other and usually interact through a communication network. The client-server relationship is created by computer programs running on the respective computers and having a client-server relationship with each other. The server can be a cloud server, also known as a cloud computing server or a cloud host, which is a host product in the cloud computing service system and solves the defects of difficult management and weak business scalability existing in traditional physical hosts and VPS services.

[0304] In particular, according to an embodiment of the present invention, the processes described above with reference to the flowcharts can be implemented as computer software programs. For example, an embodiment of the present invention includes a computer program product comprising a computer program carried on a non-transitory computer-readable medium, the computer program including program code for performing the methods shown in the flowcharts. In such an embodiment, the computer program can be downloaded and installed from a network through the communication unit 19, or installed from the storage unit 18, or installed from the ROM 12. When the computer program is executed by the processor 11, the above functions defined in the methods of the embodiments of the present invention are performed.

[0305] It should be understood that the various forms of processes shown above can be used, with steps reordered, added or deleted. 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 no limitation is imposed herein.

[0306] The above specific embodiments do not constitute a limitation on the protection scope of the present invention. Those skilled in the art should understand that various modifications, combinations, sub - combinations and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions and improvements made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A method for reconstructing an energy spectrum CT image, characterized in that, include: During the current iterative reconstruction process of the energy spectrum CT image, for each of the at least two energy spectra, determining transformed projection data based on the projection matrix and original projection data corresponding to the energy spectrum, and a previous basis image reconstructed in a previous iterative reconstruction process; Substituting the transformed projection data and the projection matrix into a pre-constructed sub-problem for the energy spectrum, and solving the sub-problem based on a first constraint preset for the sub-problem to obtain a transformed basis image that satisfies the first constraint and has a minimum error with the transformed projection data; The transformed basis images corresponding to the energy spectra are combined to obtain a combined basis image, and the combined basis image is substituted into a pre-constructed optimization problem to reconstruct a current basis image with the minimum error between the current basis image and the combined basis image, thereby realizing reconstruction of the energy spectrum CT image.

2. The method for reconstructing an energy spectrum CT image according to claim 1, wherein The first constraint is used to constrain the transformation basis image to be solved by the subproblem, and the subproblem is solved in the following manner: For the candidate transformation basis images that satisfy the first constraint, first projection data is calculated based on the candidate transformation basis images and the projection matrix, and the error between the first projection data and the transformation projection data is determined to solve the candidate transformation basis image corresponding to the smallest error among the errors corresponding to each of the candidate transformation basis images as the transformation basis image.

3. The energy spectrum CT image reconstruction method according to claim 1, characterized in that, Solving the substituted subproblem based on the first constraint preset for the subproblem includes: A preset primal-dual hybrid gradient algorithm is obtained, and the primal-dual hybrid gradient algorithm is used to solve the subproblem after substitution based on the first constraint preset for the subproblem.

4. The energy spectrum CT image reconstruction method according to claim 1, wherein, The first constraint includes a non-negativity constraint and / or a total variation constraint.

5. The spectral CT image reconstruction method according to claim 1, characterized in that: Substituting the combined base image into a pre-constructed optimization problem to reconstruct a current base image having a minimum error with the combined base image includes: Obtaining a pre-constructed optimization problem and a second constraint preset for the optimization problem; Substitute the combined base image into the optimization problem, and solve the optimization problem based on the second constraint to reconstruct a current base image that satisfies the second constraint and has the smallest error with the combined base image.

6. The spectral CT image reconstruction method according to any one of claims 1 to 5, characterized in that: The optimization problem is expressed by minimizing the first equation: The first formula is obtained by performing a weighted norm on the second formula at least through a weight matrix; The second equation represents the error between the third equation and the combined basis image; The third equation is represented by a preset constant matrix and a basis image to be solved, wherein the current basis image is a basis image obtained by solving the basis image to be solved.

7. The method for reconstructing an energy spectrum CT image according to claim 6, wherein The method involves energy spectra and base materials, where and are both integers greater than 1. If the current base image is reconstructed based on the optimization problem and a second constraint preset for the optimization problem, then: wherein the second constraint is located in a real vector space, when the constant matrix is invertible and the weight matrix is a positive definite matrix, the current basis image is obtained according to the inverse matrix of the constant matrix and the combined basis image; and / or, wherein the second constraint is a non - negative constraint, when the constant matrix is invertible and the weight matrix is the product of the inverse matrix and the transpose of the inverse matrix, the current basis image is obtained by performing non - negative projection on the estimated basis image, and the estimated basis image is obtained from the inverse matrix and the combined basis image; and / or, When the second constraint includes a non-negative constraint and a total variation constraint, and the weight matrix is an identity matrix, the current basis image is solved based on a preset primal-dual hybrid gradient algorithm.

8. The method for reconstructing an energy spectrum CT image according to claim 6, wherein, The method involves energy spectra,[[]] energy nodes, and base materials, where , and are all integers greater than 1. The constant matrix is obtained according to the X-ray energy spectrum distribution matrix and the mass attenuation coefficient matrix; Among them, the X-ray energy spectrum distribution matrix includes proportions, and the proportion represents the proportion of the number of photons at the corresponding energy spectrum and energy node to the total number of photons; The mass attenuation coefficient matrix includes mass attenuation coefficients, and the mass attenuation coefficient represents the mass attenuation coefficient of the corresponding base material at the corresponding energy node.

9. The spectral CT image reconstruction method according to claim 1, wherein, Determining transformed projection data based on a projection matrix corresponding to the energy spectrum, original projection data, and a previous basis image reconstructed in a previous iterative reconstruction process includes: Obtaining a projection matrix corresponding to the energy spectrum, original projection data, and an imaging operator; Obtaining a previous basis image reconstructed in a previous iterative reconstruction process; Obtaining second projection data according to the imaging operator and the previous basis image, and obtaining third projection data according to the projection matrix and the previous basis image; Obtaining transformed projection data according to the second projection data, the original projection data, and the third projection data.

10. A computer-readable storage medium, characterized in that, Computer instructions are stored in the computer-readable storage medium, and the computer instructions are used to cause a processor to implement the energy spectrum CT image reconstruction method according to any one of claims 1-9 when executed.

Images