CN116168102A - Multi-energy spectrum CT iterative reconstruction method with scatter correction and self-adaptive direction and step length - Google Patents
Multi-energy spectrum CT iterative reconstruction method with scatter correction and self-adaptive direction and step length Download PDFInfo
- Publication number
- CN116168102A CN116168102A CN202310021609.6A CN202310021609A CN116168102A CN 116168102 A CN116168102 A CN 116168102A CN 202310021609 A CN202310021609 A CN 202310021609A CN 116168102 A CN116168102 A CN 116168102A
- Authority
- CN
- China
- Prior art keywords
- base material
- energy spectrum
- formula
- projection
- energy
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000001228 spectrum Methods 0.000 title claims abstract description 55
- 238000000034 method Methods 0.000 title claims abstract description 43
- 238000012937 correction Methods 0.000 title claims abstract description 34
- 239000000463 material Substances 0.000 claims abstract description 57
- 238000009826 distribution Methods 0.000 claims abstract description 24
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 19
- 230000003044 adaptive effect Effects 0.000 claims description 10
- 239000011159 matrix material Substances 0.000 claims description 10
- 238000005259 measurement Methods 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- OIGNJSKKLXVSLS-VWUMJDOOSA-N prednisolone Chemical compound O=C1C=C[C@]2(C)[C@H]3[C@@H](O)C[C@](C)([C@@](CC4)(O)C(=O)CO)[C@@H]4[C@@H]3CCC2=C1 OIGNJSKKLXVSLS-VWUMJDOOSA-N 0.000 claims description 3
- 239000000758 substrate Substances 0.000 abstract description 11
- 238000000354 decomposition reaction Methods 0.000 abstract description 10
- 238000004445 quantitative analysis Methods 0.000 abstract description 2
- 238000005457 optimization Methods 0.000 description 7
- 239000000126 substance Substances 0.000 description 6
- 230000000694 effects Effects 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 4
- 210000000988 bone and bone Anatomy 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 3
- 238000002591 computed tomography Methods 0.000 description 3
- 238000013170 computed tomography imaging Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000012360 testing method Methods 0.000 description 3
- 238000001514 detection method Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000012805 post-processing Methods 0.000 description 2
- 230000003595 spectral effect Effects 0.000 description 2
- 230000017105 transposition Effects 0.000 description 2
- 229920002430 Fibre-reinforced plastic Polymers 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 125000003275 alpha amino acid group Chemical group 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 239000013078 crystal Substances 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 238000004836 empirical method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000011151 fibre-reinforced plastic Substances 0.000 description 1
- 235000003642 hunger Nutrition 0.000 description 1
- 230000036039 immunity Effects 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000002203 pretreatment Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000005316 response function Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000037351 starvation Effects 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
The invention discloses a multi-energy spectrum CT iterative reconstruction method with scatter correction and self-adaptive direction and step length, which comprises the following steps: step 1, setting a multi-energy spectrum CT projection reconstruction model with scattering correction; step 2, solving the increment of the distribution projection value of the base material by angle by using an iterative algorithm, and updating the current base material image into a new base material image; wherein, the convergent substrate projection image is obtained rapidly by the self-adaptive direction and step length. By using the scheme provided by the invention, after scattering is removed, the base material decomposition result is more accurate, the reconstruction result has fewer artifacts, and a foundation is laid for quantitative analysis.
Description
Technical Field
The invention relates to the field of computed tomography (Computed Tomography, CT), in particular to a multi-energy-spectrum CT iterative reconstruction method with self-adaptive direction and step length and scattering correction.
Background
The computer tomography technology is used for reconstructing an image of the internal structure of the measured object through a specific algorithm and a computer by collecting X-ray attenuation data of the measured object under a series of different angles without damaging the measured object. CT imaging has the advantages of non-contact, no image overlapping, high resolution, easy processing, storage, transmission and the like, and has been widely applied to a plurality of fields of medicine, industry, biology, petroleum, materials, security inspection and the like.
The traditional CT image reconstruction model and algorithm are based on the assumption of single-energy rays, and the reconstructed image is the spatial distribution of linear attenuation coefficients of the measured object to X-rays under certain energy. However, the X-rays generated by the X-ray machines used in CT currently widely used in clinic, industry, etc. are all multi-functional, and the corresponding CT image is actually reconstructed from a set of scan data of rays under a certain energy spectrum by a single energy CT image reconstruction algorithm. The reconstructed CT image can be approximately regarded as a spatial distribution of linear attenuation coefficients of the object under test corresponding to X-rays at a certain energy. Because the image reconstruction model is inconsistent with the physical process of data acquisition, and a plurality of factors such as photon scattering, detector efficiency, signal to noise ratio of data and the like are added, the CT value of the reconstructed image can be distorted to generate a homogeneous heterology (the same substance corresponds to different CT values) and a heterology homology (the same CT value corresponding to different substances). Typical image CT value distortions are "hardening artifacts", "ringing artifacts", "photon starvation artifacts", "scattering artifacts", "half-volume effect artifacts", etc. Therefore, conventional clinical CT and industrial CT are relatively limited in their material discrimination capability.
The distinguishing capability of the dual-energy spectrum CT on substances is obviously higher than that of the traditional CT, and a series of successful applications are achieved in the fields of medical science and industrial CT imaging. The linear attenuation coefficient of the object is approximately expressed as the linear combination of two basis functions in the dual-energy spectrum CT imaging, and the coefficient distribution (namely the basis image) of the two basis functions is reconstructed by acquiring two groups of X-ray scanning data corresponding to the two energy spectrums. The equivalent atomic number and equivalent electron density of the measured object can be obtained through the two base images, and the substances can be distinguished and quantitatively analyzed.
The problem of dual-energy spectrum CT image reconstruction can be reduced to a solving problem of a nonlinear integral equation system in mathematics, namely, two basic images are reconstructed by two groups of nonlinear projections corresponding to two energy spectrums. The approximation of the linear attenuation coefficient of a typical measured object is represented by two classes, one is to select two known materials, the linear attenuation coefficient depending only on the X-ray energy as a basis function (often called substrate decomposition), and the other is to select the attenuation coefficient of the substance for photoelectric absorption of X-rays and the attenuation coefficient of compton scattering as a basis function (often called effect decomposition). Energy spectrum CT image reconstruction methods can be divided into three main categories: a pre-processing method, a post-processing method and a direct iterative reconstruction method. Three process differences are illustrated by the substrate decomposition, with similar effects. The pretreatment method comprises the steps of decomposing substances in a projection domain to obtain base image projection data of corresponding base materials, and reconstructing a base image of the base materials according to a conventional CT reconstruction algorithm. Such methods typically use known spectra to orthographically project substrates of different thicknesses, establish a look-up table of the relationship between thickness and projection, or use higher order polynomials to establish the relationship between line integrals from attenuated projections to the base image, and then use different thickness combinations of the base material to perform a calibration fit. The method is only suitable for the condition that projection data of different energy spectrums are compatible, such as energy spectrum projection data acquired by a double-layer detector or a photon counting detector. The post-processing method is to reconstruct images of nonlinear projection data under each energy spectrum respectively, and then realize material decomposition in an image domain. Classical representatives of this class of methods are the EDEC empirical method, the MDIR method, and various algorithms developed based on photon counting detectors. The direct iterative reconstruction method directly solves the base image in a nonlinear integral equation in an iterative mode. Representative of such methods are ASD-POCS, ASD-NC-POCS and iterative solutions based on linear approximation of nonlinear integral equations. Typical examples of approximate iterative solutions that linearize nonlinear integral equations are the E-ART algorithm, the E-SART algorithm, the IFBP algorithm, and the tilt correction algorithm OPMT that accelerates the E-ART extension.
However, the existing approximate iteration solving method for linearizing the nonlinear integral equation has the problems of slow convergence rate, weak noise immunity, inadaptation to the condition of inconsistent ray paths, neglecting influence of scattering signals on the multi-energy-spectrum CT image and the like.
Disclosure of Invention
It is an object of the present invention to provide a multi-energy spectral CT iterative reconstruction method with scatter correction with adaptive direction and step size, which overcomes or at least alleviates at least one of the above-mentioned drawbacks of the prior art.
In order to achieve the above object, the present invention provides a multi-energy spectrum CT iterative reconstruction method with scatter correction and adaptive direction and step size, comprising:
step 1, a multi-energy spectrum CT projection reconstruction model with scattering correction is represented by the following formula:
wherein ：
for polychromatic projection data of X-rays with scatter correction along path l in the kth energy spectrum, F (l) = (F) 1 (l),F 2 (l),…,F M (l)) τ ＝(∫ l f 1 (x)dl,∫ l f 2 (x)dl,…,∫ l f M (x)dl) τ Representing the projection values of the distribution of M base materials along the ray path l, F m (l)＝∫ l f m (x) dl, m=1, 2, …, M is the number of base materials, f m (x) Representing the density distribution of the mth substrate in the measured object, and τ represents the transposition of the vector; e (E) min and Emax The lowest energy and the highest energy of the photons, respectively; s is S k (E) Is the normalized equivalent energy spectrum of the kth energy spectrum; mu (mu) m (E) M=1, 2, …, M being the linear attenuation coefficient of the mth base material; />
For the kth energy spectrum S k (E) A kind of electronic deviceThe collection of ray paths, K is the number of energy spectrums;
recording the current distribution projection value of the base material as
The value of the distribution projection of the base material estimated after n iterations is represented, n represents the number of iterations, and n=0 during initialization;
for all rays at the current angle
Calculated to obtain
And the scattering intensity Sc l And by making->
Calculating DeltaF (n) (l)；
Let F (n+1) (l)＝F (n) (l)+ΔF (n) (l) Updating to obtain new base material distribution projection value F (n+1) (l) And back projecting to obtain the material image after the n+1st iteration update.
Preferably, ΔF is calculated in step 2 (n) (l) Comprising the following steps:
for function P k (F (n+1) (l) At point F) (n) (l) At this point, a first-order taylor expansion is performed, resulting in the following expression:
wherein ,
as a function P k Regarding F (l) at F (n) (l) Partial derivatives at the points, representing the respective base materialsA weighted average of the attenuation coefficients of (a) under the current path;
p pair P k,F (F (n) (l) Normalized to obtain
wherein ,
wherein ,
for the current point F (n) (l) A directed distance to a hyperplane represented by formula (6);
solving according to formulas (6) and (7) to obtain delta F (n) (l)。
Preferably, the ΔF is obtained by solving according to formulas (6) and (7) (n) (l) Comprising the following steps:
let the ray under the current angle have r compatible measurement data, r is less than or equal to K, and is obtained by the formula (7):
wherein the formula (9) is a matrix form of the formula (8),
taking the direction of the minimum solution of equation (9) as the main direction, solving the minimum solution of equation (9) is equivalent to solving the following equation (10):
wherein ,ΔF1 (n) (l) The expression is that the formula (10) is established
Minimum ΔF (n) (l) A value;
the main direction is expressed as:
the following formulas (15) and (16) are set in the auxiliary direction:
wherein the formula (16) is a matrix form of the formula (15),
and (3) making:
According toPrincipal direction
And auxiliary direction->
To obtain self-adaptive descending direction
Wherein a is a constant of 0 or more;
setting DeltaF (n) (l) Is constrained by the furthest distance of the current point from each hyperplane, i.e
k=1, 2, …, K, c is a constant;
step length s (n) (l) Is F (n) (l) The weighted average of the absolute distances to r hyperplanes is expressed by:
wherein ,
represents F (n) (l) Directed distance to kth hyperplane, ω k Representation->
Weights of (2);
with adaptive descent direction and step size s (n) (l) Obtaining
ΔF (n) (l)＝s (n) (l)D (n) (l)……(20)。
Due to the adoption of the technical scheme, the invention has the following advantages:
by using the scheme provided by the invention, after scattering is removed, the base material decomposition result is more accurate, the reconstruction result has fewer artifacts, and a foundation is laid for quantitative analysis.
Drawings
Fig. 1 is a schematic flow chart of a multi-energy spectrum CT iterative reconstruction method with scatter correction and adaptive direction and step size according to an embodiment of the present invention.
Fig. 2 is a schematic diagram of the result of applying a multi-energy spectrum CT iterative reconstruction method with scatter correction with adaptive direction and step size to a water-based image according to an embodiment of the present invention.
Fig. 3 is a schematic diagram of a result of applying a multi-energy spectrum CT iterative reconstruction method with scatter correction in a self-adaptive direction and step size to a bone-based image according to an embodiment of the present invention.
Fig. 4 is a schematic diagram of a result of applying a multi-energy spectrum CT iterative reconstruction method with scatter correction with adaptive direction and step size to a pseudo-single-energy image according to an embodiment of the present invention.
Detailed Description
In the drawings, the same or similar reference numerals are used to denote the same or similar elements or elements having the same or similar functions. Embodiments of the present invention will be described in detail below with reference to the accompanying drawings.
In the description of the present invention, the terms "center", "longitudinal", "lateral", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", and the like indicate an orientation or a positional relationship based on that shown in the drawings, only for convenience of description and simplification of the description, and do not indicate or imply that the apparatus or element to be referred to must have a specific orientation, be configured and operated in a specific orientation, and therefore should not be construed as limiting the scope of protection of the present invention.
In the case of no conflict, the technical features in the embodiments and the implementation modes of the present invention may be combined with each other, and are not limited to the embodiments or implementation modes where the technical features are located.
The invention will be further described with reference to the drawings and the specific embodiments, it being noted that the technical solution and the design principle of the invention will be described in detail with only one optimized technical solution, but the scope of the invention is not limited thereto.
The following terms are referred to herein, and for ease of understanding, the meaning thereof is described below. It will be understood by those skilled in the art that other names are possible for the following terms, but any other name should be construed to be consistent with the terms set forth herein without departing from their meaning.
The embodiment of the invention provides a multi-energy spectrum CT iterative reconstruction method with self-adaptive direction and step length and scattering correction, which is shown in figure 1 and comprises the following steps:
step 1, a multi-energy spectrum CT projection reconstruction model with scattering correction is represented by the following formula:
wherein ：
for polychromatic projection data of X-rays with scatter correction along path l in the kth energy spectrum, F (l) = (F) 1 (l),F 2 (l),…,F M (l)) τ ＝(∫ l f 1 (x)dl,∫ l f 2 (x)dl,…,∫ l f M (x)dl) τ Representing the projection values of the distribution of M base materials along the ray path l, F m (l)＝∫ l f m (x) dl, m=1, 2, …, M is the number of base materials, f m (x) Representing the density distribution of the mth substrate in the measured object, and τ represents the transposition of the vector; e (E) min and Emax The lowest energy and the highest energy of the photons, respectively; s is S k (E) Is the normalized equivalent energy spectrum of the kth energy spectrum; mu (mu) m (E) M=1, 2, …, M being the linear attenuation coefficient of the mth base material; />
For the kth energy spectrum S k (E) K is the number of energy spectrums.
Wherein, the multi-energy spectrum projection model expression (1) with scattering correction is obtained by the following way:
when the multi-energy X-ray passes through the object, the linear attenuation coefficient of the multi-energy X-ray is related to the measured object and the energy of the X-ray. Thus, the intensity of X-rays passing through the object under test can be expressed as
wherein I0 (l) Representing an initial flow intensity of the X-rays along path l; i (l) represents the intensity of flow along the X-ray path l after passing through the object under test; μ (X, E) is the linear attenuation coefficient of an object with a spatial position X for X-rays with energy E; s is S k (E) Normalized equivalent energy spectrum (related to radiation source, filter, detector crystal, etc.); e (E) min and Emax The lowest energy and the highest energy of the photons, respectively; l is the ray path; sc (Sc) l Is the intensity of the scattered signal received by the detector pixel corresponding to path l;
for the kth energy spectrum S k (E) K is the number of energy spectrums.
Considering the influence of scattering signals on the quality of multi-energy-spectrum CT images, the measured attenuation data are firstly measured
And performing scattering correction. After each update of the projection image of the base material, the scatter signal intensity Sc received by the detector pixels corresponding to the ray path l is simulated by using a method gQMC fd based on the combination of quasi-Monte Carlo and forced detection l . Recording the attenuation data after scattering correction as
The above-mentioned steps are logarithmically converted, so that the following multi-energy spectrum projection model with scattering correction can be obtained
wherein ,
is polychromatic projection data along path l with scatter correction. Mu (x, E) is decomposed into the following forms based on the method of decomposing the base material
wherein μm (E) M=1, 2, …, M is the linear attenuation coefficient of the M-th base material, f m (x) M=1, 2, …, M is the density distribution of the mth base material in the measured object, and M is the number of base materials. Note F (l) = (F) 1 (l),F 2 (l),…,F M (l)) τ ＝(∫ l f 1 (x)dl,∫ l f 2 (x)dl,…,∫ l f M (x)dl) τ The projection values of the distribution of the M base materials along the ray path l. Therefore, the multi-energy spectrum projection model with scattering correction can be converted into an M-element nonlinear equation set with the following form
And 2, solving the increment of the distribution projection value of the base material by angle by using an iterative algorithm, and updating the current base material image into a new base material image.
The invention solves the material image f by using the multi-energy spectrum CT iterative reconstruction method ADS-SC iteration with the self-adaptive direction and step length and the scattering correction m (x) The linear attenuation coefficient mu of the base material is predetermined before m=1, 2, …, M m (E) And at each iterationUpdating the base material image f m (x) Then, the intensity of scattered signals received by the detector pixels corresponding to the ray path l is simulated, and then the energy spectrum S is estimated by using an EM algorithm with scattering correction k (E) A. The invention relates to a method for producing a fibre-reinforced plastic composite For example, the linear attenuation coefficient μ of the base material is given in advance by table look-up m (E) M=1, 2, …, M, and updates the material image f at each iteration m (x) After m=1, 2, …, M, the scatter signal intensities Sc received by the detector pixels corresponding to the ray path l are rapidly modeled by a method gQMCFFD based on a combination of quasi-monte carlo and forced detection l Energy spectrum S is estimated by using EM algorithm with scattering correction k (E) K=1, 2, …, K, the detector response function is known.
Recording the current distribution projection value of the base material as
The value of the distribution projection of the base material estimated after n iterations is represented, n represents the number of iterations, and n=0 at the time of initialization.
For all rays at the current angle
Calculated to obtain
And the scattering intensity Sc l And by making->
Calculating DeltaF (n) (l)；
Further, let F (n+1) (l)＝F (n) (l)+ΔF (n) (l) Updating to obtain new base material distribution projection value F (n+1) (l) And back-projecting by using SART algorithm to obtain the n+1th iteration updated base material image
/>
Wherein DeltaF is calculated (n) (l) The process of (1) specifically comprises:
for function P k (F (n+1) (l) At point F) (n) (l) At this point, a first-order taylor expansion is performed, resulting in the following expression:
wherein ,
as a function P k Regarding F (l) at F (n) (l) The partial derivative of the position, also called the normal direction of the hyperplane (6), represents the weighted average of the attenuation coefficients of the respective base materials under the current path;
p pair P k,F (F (n) (l) Normalized to obtain
wherein ,
from equations (6) and (7), ΔF can be obtained by solving (n) (l)。
In the formula (7), the amino acid sequence of the compound,
is the current point F (n) (l) Directed distance to the hyperplane (the hyperplane represented by equation (6)).
Solving according to formulas (6) and (7) to obtain delta F (n) (l) Comprising the following steps:
by adapting ΔF (n) (l) Direction D of (2) (n) (l) Sum step s (n) (l) To quickly obtain a converged substrate projection image F (l) = (F) 1 (l),F 2 (l),…,F M (l)) τ 。
In the present inventionThrough the main direction
And auxiliary direction->
Is combined to achieve DeltaF (n) (l) Direction D (n) (l) Is provided. The main direction is given below +.>
Auxiliary direction->
Step s (n) (l) Is a calculation method of (a).
Let the ray under the current angle have r compatible measurement data, r is less than or equal to K, and is obtained by the formula (7):
wherein the formula (9) is a matrix form of the formula (4),
taking the direction of the minimum solution of equation (9) as the main direction, solving the minimum solution of equation (9) is equivalent to solving the optimization problem shown in equation (10) below:
wherein ,ΔF1 (n) (l) The expression is such that the expression (10) is established
Minimum ΔF (n) (l) Values.
The main direction is expressed as:
wherein the process of deriving the expression of the main direction according to formula (10) includes:
when matrix
When the rank of (2) is r, the optimization problem is quadratic convex optimization, and a unique solution exists. A lagrangian function solution can be constructed:
wherein λ＝(λ1 ,λ 2 ,…λ r ) τ The partial derivative of each variable is calculated for the function, and the derivative is zero, thus obtaining
Delta F is given by formula (12) (n) (l) In the form of solution, substituting formula (12) into formula (13) to obtain the equation set of parameters
wherein 
Covariance matrix of>
The matrix is positive determined for symmetry. Solving the equation set (14), substituting the equation set (12) to obtain the solution of the optimization problem (10) and obtain the main direction
Is an expression of (2).
The auxiliary direction is described below
Is a solution to (c). Wherein at the current point F (n) (l) The direction of the individual equations of (1) is known, i.e. function P k (F (n) (l) Partial derivative->
Is computable. Projection data with scatter correction, which differ in the first r functions +.>
Known then K-r measurement +.>
Unknown, a reasonable idea is to keep the unknown function value stationary (moving only along the tangent plane) while the increment is calculated.
Based on the above analysis, in the present invention, the following formulas (15) and (16) are set in the auxiliary direction:
wherein the formula (16) is a matrix form of the formula (15),
and (3) making: />
Wherein equation (16) may be directly solved, or may be indefinite, or overdetermined. The present invention contemplates the form of a generic solution that introduces unconstrained optimization with parameter lambda
If the optimization model is not directly solved, the method can be used for solving the condition that the formula (16) is overdetermined (m is more than K) by a least square method or manually discarding m-K equations, and then obtaining the optimization model according to a main direction solving method
Then, according to the main direction
And auxiliary direction->
To obtain self-adaptive descending direction
Wherein a is a constant equal to or greater than 0.
In the present invention, step s (n) (l) Is F (n) (l) The weighted average of the absolute distances to r hyperplanes is expressed by:
wherein ,
represents F (n) (l) Directed distance to kth hyperplane, ω k Representation->
Weights of (2);
with adaptive descent direction and step size s (n) (l) Obtaining
ΔF (n) (l)＝s (n) (l)D (n) (l)……(20)。
In the present invention, the adaptive step length means the correction amount Δf (n) (l) Can be constrained by the furthest distance of the current point from each hyperplane, i.e
k=1, 2, …, K, c is a constant. Step s is selected herein (n) (l) The above equation (19) is obtained as a weighted average of the absolute distances to the respective hyperplanes.
In the invention, in the initial stage of iterative calculation, in order to accelerate the convergence speed of the reconstructed image, the step length s can be given (n) (l) Multiplying a constant greater than 1, but requiring an empirical calculation
Where ε is an arbitrarily small number greater than 0. If there is energy spectrum data
If the drop condition (21) is not satisfied, the step s is reduced (n) (l)。
In summary, the steps of the multi-energy spectrum CT iterative reconstruction method with scatter correction in the adaptive direction and step size provided by the present invention can be summarized as follows:
experimental results
In order to illustrate the technical effects of the method provided by the invention, experimental results are now provided as follows.
The model of this experiment was a Shepp-Logan model, consisting of water and bone base (ideal images on the left in fig. 2 and 3). Fig. 2 is a water-based image result (gray window [0,1.55 ]), fig. 3 is a bone-based image result (gray window [0,1.05 ]), and fig. 4 is a pseudo-mono-energetic (50 keV) image result (gray window [0.0113,0.0352 ]).
And simulating and generating high-low energy projection images with scattering and noise by using an MCGPU program, and then carrying out a substrate decomposition experiment. The simulation parameters are as follows: sod=500 mm, sdd=1000 mm, detector matrix 512x512, pixel size 0.8mm, projection angle number 360, and 2-38 photons simulated with two spectral curves of 80kv and 125kv, respectively, to obtain projection images at each angle. Where SOD represents the distance from the X-ray source S to the object to be measured, and SDD represents the distance from the X-ray source S to the object to be measured. The high and low energy data are incompatible and are spaced 0.5 degrees apart.
The projected data is directly reconstructed by using a substrate decomposition reconstruction algorithm, one is directly reconstructed, and the other is directly reconstructed by adding scattering removal in an iterative process, namely, the iterative algorithm ADS-SC with the self-adaptive direction and step length and the scattering correction is used for reconstruction. The second and third columns of fig. 2 and 3 show the reconstruction results using the conventional substrate decomposition reconstruction method and the reconstruction results using the ADS-SC method of the present invention, respectively. From the results of fig. 2 and 3, it can be seen that the image reconstructed with ADS-SC is closer to the ideal image.
The algorithm of substrate iteration carries out 10 decomposition iterations, and the time is 21s. The method provided by the invention is added with a scattering removal module on the basis, 120s is used for carrying out one round of scattering removal, two rounds of iteration is generally carried out, and 282s is used for total consumption. From the decomposed water-based result and the synthesized pseudo-single-energy image (figure 4), the reconstructed result after scattering removal is closer to the original image, which shows that the multi-energy spectrum CT iterative reconstruction method ADS-SC with the self-adaptive direction and step length and the scattering correction is effective.
Finally, it should be pointed out that: the above embodiments are only for illustrating the technical solution of the present invention, and are not limiting. Those of ordinary skill in the art will appreciate that: the technical schemes described in the foregoing embodiments may be modified or some of the technical features may be replaced equivalently; such modifications and substitutions do not depart from the spirit and scope of the technical solutions of the embodiments of the present invention.
Claims (3)
1. The multi-energy spectrum CT iterative reconstruction method with the scatter correction and the self-adaption direction and step length is characterized by comprising the following steps of:
step 1, a multi-energy spectrum CT projection reconstruction model with scattering correction is represented by the following formula:
wherein ：
for polychromatic projection data of X-rays with scatter correction along path l in the kth energy spectrum, F (l) = (F) 1 (l),F 2 (l),…,F M (l)) τ ＝(∫ l f 1 (x)dl,∫ l f 2 (x)dl,…,∫ l f M (x)dl) τ Representing the projection values of the distribution of M base materials along the ray path l, F m (l)＝∫ l f m (x) dl, m=1, 2, …, M is the number of base materials, f m (x) Represents the density distribution of the mth base material in the measured object, τ represents the vector quantityIs a transpose of (2); e (E) min and Emax The lowest energy and the highest energy of the photons, respectively; s is S k (E) Is the normalized equivalent energy spectrum of the kth energy spectrum; mu (mu) m (E) M=1, 2, …, M being the linear attenuation coefficient of the mth base material; />
For the kth energy spectrum S k (E) K is the number of energy spectrums;
step 2, solving the increment of the distribution projection value of the base material by angle by using an iterative algorithm, and updating the current base material image into a new base material image; comprising the following steps:
recording the current distribution projection value of the base material as
The value of the distribution projection of the base material estimated after n iterations is represented, n represents the number of iterations, and n=0 during initialization;
for all rays at the current angle
Calculated to obtain
And the scattering intensity Sc l And by making
Calculating DeltaF (n) (l)；
Let F (n+1) (l)＝F (n) (l)+ΔF (n) (l) Updating to obtain new base material distribution projection value F (n+1) (l) And back projecting to obtain the material image after the n+1st iteration update.
2. The method according to claim 1, wherein Δf is calculated in step 2 (n) (l) Comprising the following steps:
for function P k (F (n+1) (l) At point F) (n) (l) At this point, a first-order taylor expansion is performed, resulting in the following expression:
wherein ,
as a function P k Regarding F (l) at F (n) (l) Partial derivative of each base material represents a weighted average of attenuation coefficients of each base material under the current path;
p pair P k,F (F (n) (l) Normalized to obtain
wherein ,
wherein ,
for the current point F (n) (l) A directed distance to a hyperplane represented by formula (6);
solving according to formulas (6) and (7) to obtain delta F (n) (l)。
3. The method according to claim 2, wherein Δf is obtained by solving according to formulas (6) and (7) (n) (l) Comprising the following steps:
let the ray under the current angle have r compatible measurement data, r is less than or equal to K, and is obtained by the formula (7):
wherein the formula (9) is a matrix form of the formula (8),
taking the direction of the minimum solution of equation (9) as the main direction, solving the minimum solution of equation (9) is equivalent to solving the following equation (10):
wherein ,ΔF1 (n) (l) The expression is that the formula (10) is established
Minimum ΔF (n) (l) A value;
the main direction is expressed as:
the following formulas (15) and (16) are set in the auxiliary direction:
wherein the formula (16) is a matrix form of the formula (15),
and (3) making:
According to the main direction
And auxiliary direction->
To obtain self-adaptive descending direction
Wherein a is a constant of 0 or more;
setting DeltaF (n) (l) Is constrained by the furthest distance of the current point from each hyperplane, i.e
c is a constant;
step length s (n) (l) Is F (n) (l) The weighted average of the absolute distances to r hyperplanes is expressed by:
wherein ,
represents F (n) (l) Directed distance to kth hyperplane, ω k Representation->
Weights of (2);
with adaptive descent direction and step size s (n) (l) Obtaining
ΔF (n) (l)＝s (n) (l)D (n) (l)……(20)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310021609.6A CN116168102A (en) | 2023-01-07 | 2023-01-07 | Multi-energy spectrum CT iterative reconstruction method with scatter correction and self-adaptive direction and step length |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310021609.6A CN116168102A (en) | 2023-01-07 | 2023-01-07 | Multi-energy spectrum CT iterative reconstruction method with scatter correction and self-adaptive direction and step length |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116168102A true CN116168102A (en) | 2023-05-26 |
Family
ID=86412644
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310021609.6A Pending CN116168102A (en) | 2023-01-07 | 2023-01-07 | Multi-energy spectrum CT iterative reconstruction method with scatter correction and self-adaptive direction and step length |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116168102A (en) |
-
2023
- 2023-01-07 CN CN202310021609.6A patent/CN116168102A/en active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107356615B (en) | Method and system for dual-energy X-ray CT | |
| US8306180B2 (en) | Image reconstruction method for high-energy, dual-energy CT system | |
| CN108010099B (en) | X-ray multi-energy spectrum CT finite angle scanning and image iterative reconstruction method | |
| JP7042806B2 (en) | Spectral Computed Tomography (CT) Spectral Calibration | |
| JP2018140165A (en) | Medical image generation device | |
| WO2019067524A1 (en) | Monochromatic ct image reconstruction from current-integrating data via machine learning | |
| CN110175957B (en) | Multi-energy CT-based material decomposition method | |
| EP3612095B1 (en) | Beam hardening correction in x-ray dark-field imaging | |
| US12073491B2 (en) | Energy weighting of photon counts for conventional imaging | |
| CN110415307B (en) | Tensor completion-based multi-energy CT imaging method and device and storage equipment thereof | |
| CN104411245A (en) | Dynamic modeling of imperfections for photon counting detectors | |
| Zhao et al. | An oblique projection modification technique (OPMT) for fast multispectral CT reconstruction | |
| Yu et al. | Medipix-based spectral Micro-CT | |
| Ding et al. | A high-resolution photon-counting breast CT system with tensor-framelet based iterative image reconstruction for radiation dose reduction | |
| Jumanazarov et al. | Significance of the spectral correction of photon counting detector response in material classification from spectral x-ray CT | |
| Higuchi et al. | X-ray energy spectrum estimation based on a virtual computed tomography system | |
| CN116168102A (en) | Multi-energy spectrum CT iterative reconstruction method with scatter correction and self-adaptive direction and step length | |
| Wei et al. | Blind separation model of multi-voltage projections for the hardening artifact correction in computed tomography | |
| Chukalina et al. | A hardware and software system for tomographic research: reconstruction via regularization | |
| Thierry et al. | Hybrid simulation of scatter intensity in industrial cone-beam computed tomography | |
| Viganò et al. | Physically corrected forward operators for induced emission tomography: a simulation study | |
| Haase et al. | Estimation of statistical weights for model-based iterative CT reconstruction | |
| Wei et al. | A Multienergy Computed Tomography Method Based on a Blind Decomposition Model for Multivoltage X-Ray Transmission Images | |
| US20240193827A1 (en) | Determining a confidence indication for deep learning image reconstruction in computed tomography | |
| Alsaffar et al. | Multi-Material Blind Beam Hardening Correction Based on Non-Linearity Adjustment of Projections |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |