A method and device for estimating energy spectrum using CT scanning data

Abstract

The present invention discloses a method and apparatus for estimating energy spectrum using CT scan data, comprising: step 1, constructing a nonlinear equation for the X-ray energy spectrum by combining a nonlinear model of the intensity I of multi-energy X-rays after passing through a phantom; step 2, converting the nonlinear equation for the X-ray energy spectrum into a spectrum estimation optimization problem by minimizing the error between the estimated attenuation data AS and the original attenuation data I. Incorporating information about the absence of hardening artifacts in the reconstructed image, the total variation minimum of the reconstructed discrete phantom image F is set as a constraint s.tmin||F|| for the spectrum estimation optimization problem. TV Step 3: Solve the spectrum estimation optimization problem to obtain the energy spectrum to be estimated S. The present invention utilizes a large amount of CT attenuation data, which is conducive to improving the robustness of energy spectrum estimation. In addition, the information that the reconstructed image has no hardening artifacts is taken into account in the energy spectrum estimation, which can improve the accuracy of energy spectrum estimation.

Description

Technical Field

[0001] The present invention relates to the technical field of X-ray computed tomography, and in particular to a method and device for estimating energy spectrum using CT scanning data. Background Art

[0002] X-ray energy spectrum estimation is widely used in some X-ray computed tomography (CT) imaging technologies. For example, X-ray energy spectrum distribution plays a vital role in CT image hardening correction, dual-spectrum CT (Computer Tomography) imaging, and radiation dose calculation. In clinical CT systems, the X-ray flux generated by the X-ray machine (X-ray intensity) is high, while the detector count rate is limited, making the X-ray energy spectrum distribution difficult to measure directly using an X-ray spectrometer. Therefore, the ability to quickly and accurately estimate the X-ray energy spectrum distribution is of great research significance and application value. Summary of the Invention

[0003] An object of the present invention is to provide a method and apparatus for estimating energy spectrum using CT scan data to overcome or at least alleviate at least one of the above-mentioned drawbacks of the prior art.

[0004] To achieve the above object, the present invention provides a method for estimating energy spectrum using CT scan data, comprising:

[0005] Step 1: Combine the nonlinear model of the intensity I of the multi-energy X-ray after passing through the phantom and construct the nonlinear equation (1) about the X-ray energy spectrum:

[0006]

[0007] Among them, I u,β represents the CT data collected by the detector unit with a rotation angle of β and coordinates u during fan-beam CT scanning; u represents the coordinates of the detector, u=1,...,U, and U is the total number of detectors in the CT system; β represents the rotation angle, β=1,...,V, and V represents the total number of CT scanning angles; S(E i ) represents the energy spectrum value after the energy spectrum S to be estimated is discretized; μ(E i ) represents the discretized attenuation coefficient value of the linear attenuation coefficient μ; F represents the discrete phantom image; represents the Radon transform of F; △E represents the interval of the discrete energy spectrum, which has a fixed value and is 1 in this example. Step 2, by minimizing the error between the estimated attenuation data AS and the original attenuation data I, the nonlinear equation (1) of the X-ray energy spectrum is transformed into the spectrum estimation optimization problem described by equation (3). Combined with the information that the reconstructed image has no hardening artifacts, the total variation minimum of the reconstructed discrete phantom image F is set as the constraint condition st min||F|| of the spectrum estimation optimization problem. TV ;

[0008] (F,S)=argmin[||AS-I|| 2 ] (3)

[0009] Where A is the CT system matrix, K=U×V,|||| 2 represents the square of the two-norm, ||F|| TV represents the total variation of the reconstructed discrete phantom image F;

[0010] Step 3: Solve the spectrum estimation optimization problem to obtain the energy spectrum S to be estimated.

[0011] Furthermore, the discrete image F of the single-material phantom satisfies the property shown in the following formula (2):

[0012]

[0013] Among them, f j is the sampling value of the image F at the jth pixel; Ω represents the area where the motif is located in the image.

[0014] Furthermore, the spectrum estimation optimization problem in step 3 is decomposed into a first subproblem and a second subproblem, and S is solved by alternating iterations;

[0015] Among them, the transition energy spectrum value S~ is obtained through the first sub-problem described by formula (4); the reconstructed image is solved by solving the second sub-problem described by formula (6);

[0016] (F0,S)=argmin[||I-AS|| 2 ] (4)

[0017]

[0018] Among them, F0 represents the initial value of the reconstructed image F, It means that the solution of the optimization problem of formula (4) is obtained.

[0019] Furthermore, in the first sub-problem, the SART algorithm provided by the following formula (5) is used to solve the equation group A -1 AS=A -1 I, obtain

[0020]

[0021] Among them, S 0 represents the initial value of the energy spectrum iteration, S n It's S 0 The result after n iterations; λ is the relaxation factor; I β It is the collection of CT data collected by the CT system under the scanning angle of β;.

[0022] Furthermore, in the second sub-problem, the extended algebraic reconstruction algorithm of the single color shown in the following equation (7) is combined with the total variation minimization to solve it:

[0023]

[0024] Among them, F k represents the result of reconstructing the image after iterating k-1 times with the initial value F0; p u,β Indicates I u,β Projection data, p u,β =-lnI u,β ; represents the projection vector of the ray path determined by u, β, where R represents the projection contribution of the j-th pixel in F to the ray path determined by u, β. u,β and The value of R is determined by the CT system. When the scanning parameters are determined, u,β and The value is determined.

[0025] The present invention also provides a device for estimating energy spectrum using CT scan data, comprising:

[0026] The nonlinear model construction unit is used to combine the nonlinear model of the intensity I of the multi-energy X-ray after passing through the phantom to construct the nonlinear equation (1) about the X-ray energy spectrum:

[0027]

[0028] Among them, I u,β represents the CT data collected by the detector unit with a rotation angle of β and coordinates u during fan-beam CT scanning; u represents the coordinates of the detector, u=1,...,U, and U is the total number of detectors in the CT system; β represents the rotation angle, β=1,...,V, and V represents the total number of CT scanning angles; S(E i ) represents the energy spectrum value after the energy spectrum S to be estimated is discretized; μ(E i) represents the discretized attenuation coefficient value of the linear attenuation coefficient μ; F represents the discrete phantom image, Indicates the Radon transform of F; △E represents the interval of the discrete energy spectrum, which has a fixed value and is 1 in this example.

[0029] The spectrum estimation optimization problem description unit is used to transform the nonlinear equation (1) of the X-ray energy spectrum into the spectrum estimation optimization problem described by equation (3) by minimizing the error between the estimated attenuation data AS and the original attenuation data I. Combined with the information that the reconstructed image has no hardening artifacts, the total variation minimum of the reconstructed discrete phantom image F is set as the constraint condition st min||F|| of the spectrum estimation optimization problem. TV ;

[0030] (F,S)=argmin[||AS-I|| 2 ]st min||F|| TV (3)

[0031] Where A is the CT system matrix, K=U×V,|||| 2 represents the square of the two-norm, ||F|| TV Represents the total variation of the reconstructed image F;

[0032] The computing unit is used to solve the spectrum estimation optimization problem and obtain the energy spectrum S to be estimated.

[0033] Furthermore, the discrete image F of the single-material phantom satisfies the property shown in the following formula (2):

[0034]

[0035] Among them, f j is the sampling value of the image F at the jth pixel; Ω represents the area where the motif is located in the image.

[0036] Furthermore, the spectrum estimation optimization problem in the computing unit is decomposed into a first subproblem and a second subproblem, and S is solved alternately and iteratively;

[0037] Among them, the transition energy spectrum value is obtained by the first sub-problem described by formula (4): The reconstructed image is solved by solving the second sub-problem described by equation (6);

[0038] (F0,S)=argmin[||I-AS|| 2 ] (4)

[0039]

[0040] Among them, F0 represents the initial value of the reconstructed image F, It means that the solution of the optimization problem of formula (4) is obtained.

[0041] Furthermore, in the first sub-problem, the SART algorithm provided by the following formula (5) is used to solve the equation group A -1 AS=A -1 I, obtain

[0042]

[0043] Among them, S 0 represents the initial value of the energy spectrum iteration, S n It's S 0 The result after n iterations; λ is the relaxation factor; I β It is the collection of CT data collected by the CT system under the scanning angle β.

[0044] Furthermore, in the second sub-problem, the extended algebraic reconstruction algorithm of the single color shown in the following equation (7) is combined with the total variation minimization to solve it:

[0045]

[0046] Among them, F k represents the result of reconstructing the image after iterating k-1 times with the initial value F0; p u,β Indicates I u,β Projection data, p u,β =-lnI u,β ; represents the projection vector of the ray path determined by u, β, where R represents the projection contribution of the j-th pixel in F to the ray path determined by u, β. u,β and The value of R is determined by the CT system. When the scanning parameters are determined, u,β and The value is determined.

[0047]

[0048] The present invention has the following advantages due to the adoption of the above technical solution:

[0049] 1. The present invention utilizes a large amount of CT attenuation data, thereby improving the robustness of energy spectrum estimation.

[0050] 2. The energy spectrum estimation of the present invention takes into account the information that the reconstructed image has no hardening artifacts, which greatly improves the accuracy of the energy spectrum estimation. DETAILED DESCRIPTION

[0051] The present invention is described in detail below with reference to the embodiments.

[0052] The method for estimating energy spectrum using CT scan data according to an embodiment of the present invention includes:

[0053] Step 1: Combine the nonlinear model of the intensity I of the multi-energy X-ray after passing through the phantom and construct the nonlinear equation (1) about the X-ray energy spectrum:

[0054]

[0055] Among them, I u,β represents the CT data collected by the detector unit with a rotation angle of β and coordinates u during fan-beam CT scanning; u represents the coordinates of the detector, u=1,...,U, and U is the total number of detectors in the CT system; β represents the rotation angle, β=1,...,V, and V is the total number of CT scanning angles; I represents the intensity of the multi-energy X-ray after passing through the phantom, which is I u,β K-dimensional column vector composed of; S represents the energy spectrum to be estimated, which consists of N discretized energy spectrum values S(E k ), S=(S(E1),...,S(E N )); μ(E i ) represents the discretized attenuation coefficient value of the linear attenuation coefficient μ, which is composed of N discretized attenuation coefficient values μ(E i ), μ=(μ(E1),...,μ(E N )); △E represents the interval of discrete energy spectrum, that is: E i -E i-1 For example, the energy spectrum is 80keV and is discrete at 1keV energy intervals. In this case, △E=1; F represents the discrete model image. represents the Radon transform of F; the discrete phantom image F of the single-material phantom satisfies the properties shown in the following formula (2):

[0056]

[0057] Among them, f j is the sampling value of the image F at the jth pixel; Ω represents the area where the motif is located in the image.

[0058] Step 2: By minimizing the error between the estimated attenuation data AS and the original attenuation data I, the nonlinear equation (1) of the X-ray energy spectrum is transformed into the spectrum estimation optimization problem described by equation (3). Combined with the information provided by equation (2) that the reconstructed image has no hardening artifacts, the total variation minimum of the reconstructed discrete phantom image F is set as the constraint condition st min||F|| of the spectrum estimation optimization problem. TV .

[0059] (F,S)=argmin[||AS-I|| 2 ]

[0060] st min||F|| TV (3)

[0061] Where A is the CT system matrix, K=U×V, U is the total number of detectors in the CT system, V is the total number of CT scanning angles, |||| 2 represents the square of the two-norm, ||F|| TV represents the total variation of the reconstructed discrete phantom image F.

[0062] Step 3: Solve the spectrum estimation optimization problem to obtain the X-ray energy spectrum distribution S to be estimated.

[0063] This embodiment uses a method similar to ADMM to decompose the spectrum estimation optimization problem into a first sub-problem and a second sub-problem, and solves them alternately and iteratively.

[0064] The first sub-problem is to solve the transition energy spectrum value based on the known F, using the measured intensity I of the multi-energy X-ray after passing through the phantom. The specific description is as follows:

[0065] (F0,S)=argmin[||I-AS|| 2 ] (4)

[0066] Where F0 is the initial value of the reconstructed image F. F0 can be obtained using the projection data p after the X-ray passes through the phantom using existing image reconstruction methods and then performing threshold segmentation. It can also be manually pre-set, for example, by setting each pixel value to 0 or another value.

[0067] There are many ways to solve the sub-problems. For example, in the first sub-problem, the SART algorithm provided by the following formula (5) can be used to solve the equation group A. -1 AS=A -1 I, obtain the solution of the optimization problem of formula (4), that is, the transition energy spectrum value The first sub-problem can also be solved using existing methods such as the EM method and the SART method.

[0068]

[0069] in, It represents the result of SART iteration n times for the i-th discrete energy spectrum value; I β It represents the set of CT data collected by the CT system under the scanning angle β. The iterative initial value S of the energy spectrum 0It can be obtained by the open source software Spectrum GUI, S n It's S 0 The result after n iterations; λ is the relaxation factor, and its specific value can be between (0,1).

[0070] The second sub-problem is to solve the reconstructed image F based on the transition energy spectrum value S~ and use the measured intensity I of the multi-energy X-ray after passing through the phantom to minimize the total variation of the reconstructed discrete phantom image F, which is specifically described as the following formula (6):

[0071]

[0072] In the second sub-problem, this embodiment can be solved by using, but not limited to, a monochrome extended algebraic reconstruction algorithm (E-ART) combined with total variation minimization (TVM). The algorithm can be abbreviated as E-ART-TVM, as shown in the following equation (7):

[0073]

[0074] Among them, F k represents the result of reconstructing the image after iterating k-1 times with the initial value F0; p u,β Indicates I u,β Projection data, p u,β =-lnI u,β ; represents the projection vector of the ray path determined by u, β, where R represents the projection contribution of the j-th pixel in F to the ray path determined by u, β. u,β and The value of R is determined by the CT system. When the scanning parameters are determined, u,β and The value is determined. and It represents an intermediate quantity used to simplify the formula and has no actual physical meaning.

[0075] The above-mentioned total variation minimization (TVM) can be solved by, but is not limited to, an optimization algorithm such as the steepest descent algorithm. The second word problem can also be solved by existing methods such as the POCS-TVM algorithm.

[0076] In summary, the specific steps for solving the spectrum estimation optimization problem are as follows:

[0077] Step 1, select initial values S0, F0;

[0078] Step 2: Set F n-1 is the result after the n-1th iteration. Solving the first subproblem yields the energy spectrum Sn ;

[0079] Step 3, S n Substitute into the second subproblem and use the E-ART-TVM algorithm to solve the second subproblem, and get F n ;

[0080] Step 4, return to step 2 until ||I-AS|| 2 <ε or n> N MAX , ε is the threshold set in advance, N MAX The maximum number of iterations of the optimization algorithm.

[0081] The method for estimating X-ray energy spectra proposed in this invention significantly increases the number of equations required for spectrum estimation, thereby improving the well-posedness of the system. Furthermore, the proposed spectrum estimation model considers the absence of hardening artifacts in the reconstructed image, significantly improving the accuracy of spectrum estimation. Numerical and practical experiments demonstrate that the proposed method can accurately estimate the X-ray energy spectrum, and that the images reconstructed using this method are free of hardening artifacts.

[0082] The present invention also provides a device for estimating energy spectrum using CT scan data, which includes a nonlinear model construction unit, a spectrum estimation optimization problem description unit and a calculation unit, wherein:

[0083] The nonlinear model construction unit is used to combine the nonlinear model of the intensity I of the multi-energy X-ray after passing through the phantom to construct a nonlinear equation (1) about the X-ray energy spectrum.

[0084] The spectrum estimation optimization problem description unit is used to transform the nonlinear equation (1) of the X-ray energy spectrum into the spectrum estimation optimization problem described by equation (3) by minimizing the error between the estimated attenuation data AS and the original attenuation data I. Combined with the information that the reconstructed image has no hardening artifacts, the total variation minimum of the reconstructed discrete image F is set as the constraint condition of the spectrum estimation optimization problem.

[0085] The computing unit is used to solve the spectrum estimation optimization problem and obtain the energy spectrum S to be estimated.

[0086] In one embodiment, the spectrum estimation optimization problem in the calculation unit is decomposed into a first sub-problem and a second sub-problem, and S is solved alternately and iteratively. Among them, the transition energy spectrum value is obtained by the first sub-problem described by formula (4). The reconstructed discrete image F is solved by the second sub-problem described by equation (6).

[0087] Finally, it should be noted that the above embodiments are intended only to illustrate the technical solutions of the present invention and are not intended to limit them. Those skilled in the art will appreciate that the technical solutions described in the aforementioned embodiments may be modified, or some of the technical features thereof may be replaced with equivalents; such modifications or replacements do not deviate from the spirit and scope of the technical solutions of the various embodiments of the present invention.

Claims

1. A method for estimating energy spectrum using CT scan data, characterized in that: include: Step 1: Combine the nonlinear model of the intensity I of the multi-energy X-ray after passing through the phantom and construct the nonlinear equation (1) about the X-ray energy spectrum: Among them, I u,β represents the CT data collected by the detector unit with a rotation angle of β and coordinates u during fan-beam CT scanning; u represents the coordinates of the detector, u=1,...,U, and U is the total number of detectors in the CT system; β represents the rotation angle, β=1,...,V, and V represents the total number of CT scanning angles; S(E i ) represents the energy spectrum value after the energy spectrum S to be estimated is discretized; μ(E i ) represents the discretized attenuation coefficient value of the linear attenuation coefficient μ; F=(f1,…,f j ,…,f J ) represents the discrete phantom image, f j is the sampling value of image F at the jth pixel; Indicates the Radon transform of F; ΔE represents the interval of the discrete energy spectrum, which has a fixed value and is 1 in this example; Step 2: By minimizing the error between the estimated attenuation data AS and the original attenuation data I, the nonlinear equation (1) of the X-ray energy spectrum is transformed into the spectrum estimation optimization problem described by equation (3). Combined with the information that the reconstructed image has no hardening artifacts, the total variation minimum of the reconstructed discrete phantom image F is set as the constraint condition of the spectrum estimation optimization problem. (F,S)=argmin[||AS-I|| 2 ] s.t min||F|| TV (3) Where A is the CT system matrix, K=U×V,|||| 2 represents the square of the two-norm, ||F|| TV represents the total variation of the reconstructed discrete phantom image F; Step 3: Solve the spectrum estimation optimization problem to obtain the energy spectrum S to be estimated.

2. The method for estimating energy spectrum using CT scan data according to claim 1, wherein: The discrete image F of a single-material phantom satisfies the following property (2): Where Ω represents the region in the image where the motif is located.

3. The method for estimating energy spectrum using CT scan data according to claim 1 or 2, wherein: Decompose the spectrum estimation optimization problem in step 3 into the first subproblem and the second subproblem, and solve S alternately and iteratively; Among them, the transition energy spectrum value is obtained by the first sub-problem described by formula (4): The reconstructed image is solved by solving the second sub-problem described by equation (6); (F0,S)=argmin[||I-AS|| 2 ](4) Among them, F0 represents the initial value of the reconstructed image F, It means that the solution of the optimization problem of formula (4) is obtained.

4. The method for estimating energy spectrum using CT scan data according to claim 3, wherein: In the first sub-problem, the SART algorithm provided by equation (5) is used to solve the equation group A -1 AS=A -1 I, obtain Among them, S 0 represents the initial value of the energy spectrum iteration, S n It's S 0 The result after n iterations; λ is the relaxation factor; I β It is the collection of CT data collected by the CT system under the scanning angle β.

5. The method for estimating energy spectrum using CT scan data according to claim 3, wherein: In the second sub-problem, the extended algebraic reconstruction algorithm of the single color shown in Equation (7) is combined with the total variation minimization to solve it: Among them, F k 、F k-1 They represent the results of reconstructing the image after the initial value F0 is iterated k times and k-1 times respectively; p u,β Indicates I u,β Projection data, p u,β =-lnI u,β ; represents the projection vector of the ray path determined by u, β, where represents the projection contribution of the j-th pixel in F to the ray path determined by u, β; 6. A device for estimating energy spectrum using CT scan data, characterized in that: include: The nonlinear model construction unit is used to combine the nonlinear model of the intensity I of the multi-energy X-ray after passing through the phantom to construct the nonlinear equation (1) about the X-ray energy spectrum: Among them, I u,β represents the CT data collected by the detector unit with a rotation angle of β and coordinates u during fan-beam CT scanning; u represents the coordinates of the detector, u=1,...,U, and U is the total number of detectors in the CT system; β represents the rotation angle, β=1,...,V, and V represents the total number of CT scanning angles; S(E i ) represents the energy spectrum value after the energy spectrum S to be estimated is discretized; μ(E i ) represents the discretized attenuation coefficient value of the linear attenuation coefficient μ; F represents the discrete phantom image, Indicates the Radon transform of F; ΔE represents the interval of discrete energy spectrum; The spectrum estimation optimization problem description unit is used to transform the nonlinear equation (1) of the X-ray energy spectrum into the spectrum estimation optimization problem described by equation (3) by minimizing the error between the estimated attenuation data AS and the original attenuation data I. Combined with the information that the reconstructed image has no hardening artifacts, the total variation minimum of the reconstructed discrete phantom image F is set as the constraint condition st min||F|| of the spectrum estimation optimization problem. TV ; (F,S)=argmin[||AS-I|| 2 ] st min||F|| TV (3) Where A is the CT system matrix, K=U×V,|||| 2 represents the square of the two-norm, ||F|| TV Represents the total variation of the reconstructed image F; The computing unit is used to solve the spectrum estimation optimization problem and obtain the energy spectrum S to be estimated.

7. In the apparatus for estimating energy spectrum using CT scan data according to claim 6, the discrete image F of the single-material phantom satisfies the property shown in the following equation (2): in, Ω represents the region in the image where the motif is located.

8. The device for estimating energy spectrum using CT scan data according to claim 6 or 7, characterized in that: Decompose the spectrum estimation optimization problem in the computing unit into the first subproblem and the second subproblem, and solve S alternately and iteratively; Among them, the transition energy spectrum value is obtained by the first sub-problem described by formula (4): The reconstructed image is solved by solving the second sub-problem described by equation (6); (F0,S)=argmin[||I-AS|| 2 ](4) Among them, F0 represents the initial value of the reconstructed image F, It means that the solution of the optimization problem of formula (4) is obtained.

9. The device for estimating energy spectrum using CT scan data according to claim 8, wherein: In the first sub-problem, the SART algorithm provided by equation (5) is used to solve the equation group A -1 AS=A -1 I, obtain Among them, S 0 represents the initial value of the energy spectrum iteration, S n It's S 0 The result after n iterations; λ is the relaxation factor; I β It is the collection of CT data collected by the CT system under the scanning angle β.

10. The device for estimating energy spectrum using CT scan data according to claim 8, wherein: In the second sub-problem, the extended algebraic reconstruction algorithm of the single color shown in Equation (7) is combined with the total variation minimization to solve it: Among them, F k represents the result of reconstructing the image after iterating k-1 times with the initial value F0; p u,β Indicates I u,β Projection data, p u,β =-lnI u,β ; represents the projection vector of the ray path determined by u, β;

Images