CN109875593B - X-ray multi-spectrum separation imaging method, storage medium and device - Google Patents
X-ray multi-spectrum separation imaging method, storage medium and device Download PDFInfo
- Publication number
- CN109875593B CN109875593B CN201910107265.4A CN201910107265A CN109875593B CN 109875593 B CN109875593 B CN 109875593B CN 201910107265 A CN201910107265 A CN 201910107265A CN 109875593 B CN109875593 B CN 109875593B
- Authority
- CN
- China
- Prior art keywords
- imaging
- ray
- decomposition model
- updating
- measured object
- 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.)
- Active
Links
Images
Landscapes
- Analysing Materials By The Use Of Radiation (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
The invention discloses an X-ray multi-spectrum separation imaging method, a storage medium and a device, wherein the method comprises the following steps: step 11: setting imaging voltage range [ V1, V2] of measured object](ii) a Step 12: in [ V1, V2]]Taking N different scanning voltages within the range of (1), wherein N is larger than the total number K of materials contained in the object to be detected, and the difference value between any two scanning voltages is larger than a first preset value, and performing full-angle scanning imaging on the object to be detected under each scanning voltage to obtain N projection image sequences meeting the imaging gray scale requirement; step 13: based on the N projection image sequences, solving projection thickness matrixes D corresponding to all materials of the measured object in a variable energy multispectral attenuation expression which meets the Poisson maximum likelihood principle; step 14: calculating a separate projection p for each narrow energy band in the X-ray based on D r R =1,2 … R; step 15: reconstruction of p r Obtaining a reconstructed image a of each narrow spectral band r . The method of the invention can realize the multi-energy spectrum CT imaging with low cost and high energy spectrum resolution.
Description
Technical Field
The invention relates to the field of instruments, in particular to an X-ray multi-spectrum separation imaging method, a storage medium and a device.
Background
In the process of X-ray imaging, because the multispectral characteristics of X-rays and the attenuation coefficient of substances generally decrease with the increase of photon energy, the spectral center of the X-ray gradually moves to a high-energy direction along with the increase of the thickness of the X-ray passing through an object, namely, a hardening effect is generated. The hardening effect destroys the linear relation between the attenuation coefficient and the projection data, so that the CT image obtained by the conventional reconstruction method generates hardening artifacts, the phenomena of 'same-image foreign matters' and 'same-object different images' occur, and the quantitative analysis is not facilitated, such as the quantitative characterization of the microstructure of the novel material, such as the component content, the ore porosity and the like.
In response to the above problems, studies related to multispectral CT imaging techniques have emerged. The typical dual-energy CT can realize the nonlinear reconstruction of projection data by acquiring the projection data of two different energy spectrums, and mainly comprises two steps of projection dual-energy decomposition and image reconstruction. However, dual-energy CT has only two energy spectra and has low discrimination, so that a multi-spectrum CT imaging technology is proposed. The photon counting detector can realize imaging of a plurality of narrow energy spectrums, and has the defects of low time/space resolution and high technical cost; the synchrotron radiation source has good monochromaticity, but the device is huge, the manufacturing cost is expensive, and the synchrotron radiation source is difficult to be applied in practical engineering.
How to realize a multi-energy spectrum (multi-spectrum) CT imaging technology with low cost and high energy spectrum resolution is a technical problem to be solved urgently at present.
Disclosure of Invention
In view of the above, the present invention provides an X-ray multi-spectrum separation imaging method, a storage medium and an apparatus, so as to solve the problem of how to implement a multi-spectrum CT imaging technique with low cost and high energy spectrum resolution.
The invention provides an X-ray multi-spectrum separation imaging method, which comprises the following steps:
step 11: setting an imaging voltage range [ V1, V2] of a measured object, wherein when V belongs to [ V1, V2], the penetration rate of the measured object meets a first preset condition, and the relative deviation of linear attenuation coefficients of any two materials of the measured object in V meets an energy range corresponding to a preset requirement and meets a second preset condition;
step 12: taking N different scanning voltages within the range of [ V1, V2], wherein N is larger than the total number K of materials contained in the object to be measured, and the difference between any two scanning voltages is larger than a first preset value, carrying out full-angle scanning imaging on the object to be measured under each scanning voltage to obtain N projection image sequences meeting the imaging gray scale requirement, wherein the full-angle scanning comprises T different imaging angles;
step 13: based on the N projection image sequences, solving projection thickness matrixes D corresponding to all materials of the measured object in a variable energy multispectral attenuation expression which meets the Poisson maximum likelihood principle;
step 14: calculating a separate projection p for each narrow energy band in the X-ray based on D r R =1,2 … R, R being the total number of narrow spectral segments of the X-ray;
step 15: reconstruction of p r Obtaining a reconstructed image a of each narrow energy spectrum band r 。
The present invention also provides a non-transitory computer readable storage medium storing instructions that, when executed by a processor, cause the processor to perform the steps in the above-described X-ray multi-spectral separation imaging method.
The invention also provides an X-ray multispectral separation imaging device, which comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, wherein the processor realizes the steps in the X-ray multispectral separation imaging method when executing the computer program.
The X-ray multispectral separation imaging method can realize the blind separation of the X-ray projection images, and selects reasonable imaging parameters to obtain an X-ray projection image sequence according to the material and the size of an object; then, according to the projection image sequence, solving a projection thickness matrix D corresponding to all materials of the measured object in a variable energy multispectral attenuation expression which meets the Poisson maximum likelihood principle; and then realizing the projection separation of the multispectral projection image according to the matrix D. Compared with the X-ray CT imaging technology based on hardware and other multi-energy spectrum reconstruction, the method realizes the energy spectrum separation of multi-spectrum X-ray images from a projection domain, does not need to improve the existing high-dynamic CT imaging equipment, has the advantages of low cost, high energy spectrum resolution and the like, and has important significance for popularizing the application of the quantitative CT imaging technology.
Drawings
FIG. 1 is a flow chart of an X-ray multispectral separation imaging method of the present invention;
FIG. 2 is a schematic view of a sequence of projection images obtained by scanning voltages;
FIG. 3 is a diagram of a sub-element subdivision of FIG. 2;
FIG. 4 is a first embodiment of solution D of step 13 of FIG. 1;
FIG. 5 is a second embodiment of solution D for step 13 of FIG. 1;
fig. 6 is a structural diagram of an X-ray multi-spectrum separation imaging device of the invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in detail with reference to the accompanying drawings and specific embodiments.
As shown in fig. 1, the present invention provides an X-ray multi-spectral separation imaging method, comprising:
step 11 (S11): setting an imaging voltage range [ V1, V2] of a measured object, wherein when V belongs to [ V1, V2], the penetration rate of the measured object meets a first preset condition, and the relative deviation of linear attenuation coefficients of any two materials of the measured object in V meets an energy range corresponding to a preset requirement and meets a second preset condition;
for example, the first preset condition: selecting proper voltage interval V according to the material and size of the measured object, the dynamic range of the detector imaging and avoiding the information loss caused by underexposure and overexposure 1 ,V 2 ]The penetration rate of the object can reach 30-70% or other ranges.
The second preset condition is as follows: from the angle of component analysis of the measured object, the linear attenuation coefficient curve of the material is used to satisfy any voltage V epsilon [ V ∈ [ ] 1 ,V 2 ]Time of the arbitrary k 1 、k 2 The energy range omega and the energy range omega generated by V of the material satisfying the formula (1) 1 The intersection of (a) occupies [ omega ] 1 More than 60% to ensure the components are distinguishable over the selected energy range.
In the formula mu k1 (E) And mu k2 (E) Respectively represent the k-th 1 、k 2 The linear attenuation coefficient of the material at the energy E. The value of 0.05 in the formula (1) can be changed to other values according to the strictness of the preset requirement.
Step 12 (S12): taking N different scanning voltages within the range of [ V1, V2], wherein N is larger than the total number K of materials contained in the object to be measured, and the difference between any two scanning voltages is larger than a first preset value, carrying out full-angle scanning imaging on the object to be measured under each scanning voltage to obtain N projection image sequences meeting the imaging gray scale requirement, wherein the full-angle scanning comprises T different imaging angles;
the N different scanning voltages may be equally spaced voltage sequences or non-equally spaced voltage sequences, and whether equally spaced or not, the N is greater than the total number K of materials included in the object to be measured, and a difference between any two scanning voltages is greater than a first preset value, where the first preset value is set according to a minimum distinguishing requirement of a ratio of intensity of a ray to incident intensity on the same ray path passing through the object to be measured, for example, the minimum distinguishing requirement is that a relative deviation of the ratio is greater than 2%.
Likewise, the T imaging angles may also be equally or unequally spaced, for example, the T different imaging angles include: initial imaging angle of theta 0 And T imaging angles with the angle interval delta theta. The value of T is determined by step 13, and usually T is more than or equal to 180, so that the requirement of step 13 for solving D can be met.
When each angle of each scanning voltage is imaged, the current of the X-ray tube and the integration time of the detector are adjusted to enable the gray value of the acquired X-ray image to be within the effective imaging range of the detector, namely the imaging gray value meets the requirement.
Step 13 (S13): based on the N projection image sequences, solving a projection thickness matrix D corresponding to all materials of the measured object in a variable energy multi-spectrum attenuation expression meeting the Poisson maximum likelihood principle;
step 14 (S14): calculating a separate projection p for each narrow energy band in the X-ray based on D r R =1,2 … R, R being the total number of narrow spectral segments of the X-ray;
matrix D = (D) obtained in step 13 k,m ) K×M Where M is the total number of subunits of the projected image sequence per scan voltage, the split projection p is obtained by equation (2) r :
Step 15 (S15): reconstruction of p r Obtaining a reconstructed image a of each narrow energy spectrum band r 。
For example by pairing p r General linear reconstruction (Yan \38228, li Lei. CT image reconstruction algorithm [ M]Scientific press 2014) to obtain a reconstructed image a for each narrow spectral band r 。
A schematic description of step 13 is given below.
Firstly, the N projection image sequences obtained in step 12 need to be converted into F = (F) n,m ) N×M (ii) a M is the total number of sub-elements of the projection image sequence per scan voltage, and M can also be understood as the total number of ray paths of X-rays per scan voltage, since each sub-element corresponds to one ray path.
Fig. 2 is a sequence of projection images obtained at a certain scanning voltage, where each image has image pixels as sub-units as shown in fig. 3, and an index is established for each sub-unit, then M = T × n1 × M1.
Taking the image pixel as a subunit of the projection image sequence obtained at each scanning voltage, and establishing an index for each subunit, an F matrix can be obtained as follows.
F is an N M matrix. Wherein f is n,m The ratio of the ray intensity (corresponding to the gray value of the cell in the image) on the mth path in the image sequence to the incident intensity is projected for the nth scan voltage. Incident intensity, generally refers to the gray value on the ray path that does not pass through the object region.
In the X-ray imaging process, the multispectral characteristic of the X-ray enables projection image data obtained by the detector to contain the attenuation of multispectral segments, and therefore the basis is provided for projection decomposition by utilizing the multispectral projection image.
Let the intensity in the X-ray path L be I and the initial intensity be I 0 Normalized spectrum is S (E) and maximum energy is E max The attenuation coefficient of an object at a spatial position x is mu (x, E), and the expression is shown in the specification according to multispectral attenuation
Let E max =max{E 1,max ,…,E n,max ,…,E N,max },E n,max Representing the maximum photon energy for the nth scan voltage. According to the actual situation, the energy interval [0,E max ]Is divided into R narrow energy spectrum segment subintervals [ E' r-1 ,E' r ]R =1,2 … R, of which 0=E' 0 <E' 1 <…<E' R =E max Usually at equal intervalsEnergy E in interval r =t·E' r +(1-t)·E' r-1 Typically, t =0.5 is taken.
If the image forming substance is composed of K materials, an energy interval [ E' r-1 ,E' r ]Sufficiently narrow that the attenuation coefficient of the kth material in this interval is constant at mu k (E r ),E r After the linear attenuation coefficient mu of each material is determined k (E r ) Available with reference to the national NIST website.
Let the thickness of the kth material of the object to be measured on the path L be d k Then, the formula (3) can be converted into
WhereinCalled spectral decomposition coefficient, satisfiesIt is usually unknown, andi.e. the decomposed projection to be solved.
And d k In relation to the ray paths, that is, the thickness of each material on different ray paths is different, the requirement of obtaining the multi-spectrum projection of the thickness of multiple materials on a single path only under a single scanning voltage is a serious underdetermined problem, so that the projection images with different energies under multiple scanning voltages need to be acquired. At the same time, the coefficient of spectral decomposition is notTherefore, the total number N of the collected scanning voltages is generally larger than the number K of the materials contained in the material.
Each projection image sequence of scanning voltage contains M ray paths in total, f n,m Ratio of ray intensity to incident intensity, s, on the mth path in the sequence of projection images for the nth scan voltage n,r Generating for the nth scan an r-th spectral decomposition coefficient, d, in the X-ray spectrum k,m The thickness of the kth material on the mth path. Assuming that the ratio of ray intensity to incident intensity on each path is an independently distributed Poisson random variable, i.e.
Let S = (S) n,r ) N×R Is a spectral decomposition coefficient matrix, and R is the total number of narrow energy spectrum segments of the X-ray; d = (D) k,m ) K×M Is a projected thickness matrix.
Step 13 is according to the known F = (F) n,m ) N×M Solving variable energy multispectral attenuation expression satisfying Poisson maximum likelihood principleAnd (4) a projection thickness matrix D corresponding to all materials of the measured object.
A first embodiment of step 13 is given below.
Measured value f for N × M Poisson multispectral projection n,m The actual data f is measured as a negative log-likelihood function in the form of a generalized KL divergence n,m And ideal valueThe likelihood function of the difference is in the form:
s, D is an independent variable, and each line of F is arranged in sequence according to the pixel distribution of the projection image sequence and contains projection data of different projection angles. The following first decomposition model of the Poisson prior based variable energy X-ray projection image is then obtained:
wherein:
S=(s n,r ) N×R for the spectral decomposition coefficient matrix, D = (D) k,m ) K×M Is a projected thickness matrix.
μ k (E r ) The k material is in E r Linear attenuation coefficient of time, E r Is the photon energy in the r-th narrow band.
Namely solving variable energy multispectral attenuation expression satisfying Poisson maximum likelihood principleAnd D and S can be converted into D and S which satisfy the first decomposition model.
Obtaining solving variables S and D meeting non-negative constraints according to the fact that the spectral decomposition coefficient and the thickness of the material on the ray path are non-negative values, and marking that S is more than or equal to 0,D is more than or equal to 0 to represent each element S in the matrix n,r ≥0,d k,m Not less than 0, in S.1 R =1 N Represents the spectral attenuation coefficient sum at each voltage of 1,1 R Representing an R-dimensional all-1-column vector.
As shown in fig. 4, in the first embodiment of step 13, calculating the projection thickness matrix D corresponding to all the materials of the measured object satisfying the first decomposition model includes:
E n,max the maximum photon energy generated for the nth scan voltage.
Step 22: fixing D, updating S meeting the first decomposition model by using a non-negative matrix decomposition method;
step 23: fixing S, and updating D meeting the first decomposition model by using an accelerated near-end gradient algorithm;
and step 24: calculating a difference value before and after the D value is updated, and judging whether the difference value meets a third preset condition; if yes, the updated D is output, and if no, the step 22 is returned.
For example, the third preset condition isOrOr | | | D i -D i-1 || F < ε; the value of epsilon is selected according to the actual situation, such as epsilon =10 -3 。
Wherein, the updating S satisfying the first decomposition model by using the non-negative matrix factorization method in step 22 includes:
wherein:
i denotes the number of iterations, and the initial value is 1.
Wherein the updating D satisfying the first decomposition model using the accelerated near-end gradient algorithm in step 23 comprises:
step 41: c i =0;
Wherein:
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
step 43: judgment of
A second embodiment of step 13 is given below.
Note f n,m And f n,m+1 The ratio of the ray intensity to the incident intensity on the adjacent ray paths is shown, and the projected thickness value d of each material on the adjacent ray paths in the formula (7) is obtained because the components of the object are distributed continuously k,m And d k,m+1 The difference is small, so smooth priors can be added to the projection thickness matrix, and the following penalty functions are established
In the formula D k =(d k,m ) M×1 Q denotes a two-or one-dimensional discrete gradient operator,is a Tikhonov (L2) regularization function, beta k For regularization parameters, a different degree of smoothness prior may be added for each material using this parameter.
The addition of equation (9) to equation (7) results in a second decomposition model as follows, where α is the regularization parameter.
Based on this, step 13 may further include:
converting the N projection image sequences to F = (F) n,m ) N×M Inputting a second decomposition model, M being the throw of each scanning voltageCalculating the total number of the subunits of the shadow image sequence, and calculating a projection thickness matrix D corresponding to all materials of the measured object meeting the second decomposition model;
S=(s n,r ) N×R is a spectral decomposition coefficient matrix; d = (D) k,m ) K×M Is a projected thickness matrix;
μ k (E r ) The k material is in E r Linear attenuation coefficient of time, E r Is the photon energy of the r narrow band;
α is a regularization parameter, which can be set to 10 -7 ≤α≤10 -3 ;
As shown in fig. 5, in the second embodiment of step 13, calculating the projection thickness matrix D corresponding to all materials of the measured object satisfying the second decomposition model includes:
step 32: fixing D, updating S meeting the second decomposition model by using a non-negative matrix decomposition method;
step 33: fixing S, and updating D meeting the second decomposition model by using an accelerated near-end gradient algorithm;
step 34: calculating a difference value before and after the D value is updated, and judging whether the difference value meets a third preset condition or not; if yes, the updated D is output, and if no, the step 32 is returned.
For example, the third preset condition isOrOr | | | D i -D i-1 || F < epsilon; the value of epsilon is selected according to the actual situation, such as epsilon =10 -3 。
Wherein, the updating S satisfying the second decomposition model by using the non-negative matrix factorization method in step 32 includes:
i denotes the number of iterations, and the initial value is 1.
Wherein the updating D satisfying the first decomposition model using the accelerated near-end gradient algorithm in step 33 comprises:
step 51: c i =0;
Wherein:
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
Step 53: judgment of
In addition, the solution of step 13 in fig. 1 may refer to the derivation of a basic solution method of a non-negative matrix model, which has strong dependency on an initial value and converges slowly around a minimum value.
The above is a description of the X-ray multispectral separation imaging method of the present invention.
The X-ray multispectral separation imaging method can realize the blind separation of the X-ray projection images, and selects reasonable imaging parameters to obtain an X-ray projection image sequence according to the material and the size of an object; then, according to the projection image sequence, solving a projection thickness matrix D corresponding to all materials of the measured object in a variable energy multispectral attenuation expression which meets the Poisson maximum likelihood principle; and then realizing the projection separation of the multispectral projection image according to the matrix D. Compared with the X-ray CT imaging technology based on hardware and other multi-energy spectrum reconstruction, the method realizes the energy spectrum separation of multi-spectrum X-ray images from a projection domain, does not need to improve the existing high-dynamic CT imaging equipment, has the advantages of low cost, high energy spectrum resolution and the like, and has important significance for popularizing the application of the quantitative CT imaging technology.
The present invention also provides a non-transitory computer readable storage medium storing instructions that, when executed by a processor, cause the processor to perform the steps in the above-described X-ray multi-spectral separation imaging method.
The invention also provides an X-ray multispectral separation imaging device, which comprises a memory, a processor and a computer program stored in the memory and capable of running on the processor, wherein the processor realizes the steps in the X-ray multispectral separation imaging method when executing the computer program.
As shown in fig. 6, the apparatus includes:
scanning voltage range setting module: setting an imaging voltage range [ V1, V2] of a measured object, wherein when V belongs to [ V1, V2], the penetration rate of the measured object meets a first preset condition, and the relative deviation of linear attenuation coefficients of any two materials of the measured object in V meets an energy range corresponding to a preset requirement and meets a second preset condition;
the variable energy full-scanning module: taking N different scanning voltages within the range of [ V1, V2], wherein N is larger than the total number K of materials contained in the object to be measured, and the difference between any two scanning voltages is larger than a first preset value, carrying out full-angle scanning imaging on the object to be measured under each scanning voltage to obtain N projection image sequences meeting the imaging gray scale requirement, wherein the full-angle scanning comprises T different imaging angles;
and D, a resolving module: based on the N projection image sequences, solving projection thickness matrixes D corresponding to all materials of the measured object in a variable energy multispectral attenuation expression which meets the Poisson maximum likelihood principle;
a separate projection calculation module: calculating a separate projection p for each narrow spectral band in the X-ray based on D r R =1,2 … R, R being the total number of narrow spectral segments of the X-ray;
a separate projection image reconstruction module: reconstruction of p r Obtaining a reconstructed image a of each narrow energy spectrum band r 。
Wherein, the D calculating module may further include:
converting the N projection image sequences to F = (F) n,m ) N×M Inputting a first decomposition model, wherein M is the total number of subunits of a projection image sequence of each scanning voltage, and calculating a projection thickness matrix D corresponding to all materials of a measured object meeting the first decomposition model;
Wherein:
S=(s n,r ) N×R is a spectral decomposition coefficient matrix;
D=(d k,m ) K×M is a projected thickness matrix;
μ k (E r ) The k material is in E r Linear attenuation coefficient of time, E r Is the photon energy of the r-th narrow band.
Then, calculating the projection thickness matrix D corresponding to all materials of the measured object meeting the first decomposition model comprises:
s updating module 1: fixing D, updating S meeting the first decomposition model by using a non-negative matrix decomposition method;
d, updating the module 1: fixing S, and updating D meeting the first decomposition model by using an accelerated near-end gradient algorithm;
a judging module 1: calculating a difference value before and after the D value is updated, and judging whether the difference value meets a third preset condition; if yes, outputting the updated D, and if not, returning to the S updating module 1.
In the S update module 1: updating S that satisfies the first decomposition model using a non-negative matrix factorization method includes:
i represents the number of iterations, and the initial value is 1.
D, updating module 1: updating D satisfying the first decomposition model using an accelerated near-end gradient algorithm comprises:
initial value setting module 2: c i =0;
Wherein:
i represents the iteration number, and the initial value is 1;
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
and a judging module 2: judgment of
If yes, executing the judging module 1, if not, making C i =C i +1, return to the D updating module 2.
Wherein, the D calculating module may further include:
converting the N projection image sequences into F = (F) n,m ) N×M Inputting a second decomposition model, wherein M is the total number of subunits of a projection image sequence of each scanning voltage, and calculating a projection thickness matrix D corresponding to all materials of the measured object meeting the second decomposition model;
S=(s n,r ) N×R is a spectral decomposition coefficient matrix;
D=(d k,m ) K×M is a projected thickness matrix;
μ k (E r ) The k material is in E r Linear attenuation coefficient of time, E r Is the photon energy of the r narrow band;
alpha is a regularization parameter;
Then, calculating the projection thickness matrix D corresponding to all the materials of the measured object meeting the second decomposition model comprises:
and an S updating module 3: fixing D, updating S meeting the second decomposition model by using a non-negative matrix decomposition method;
d, updating the module 3: fixing S, and updating D meeting the second decomposition model by using an accelerated near-end gradient algorithm;
a judging module 3: calculating a difference value before and after the D value is updated, and judging whether the difference value meets a third preset condition or not; if yes, the updated D is output, and if not, the S is returned to the updating module 3.
In the S update module 3: updating S satisfying the second decomposition model using a non-negative matrix factorization method comprises:
i denotes the number of iterations, and the initial value is 1.
D, updating module 3: updating D satisfying the second decomposition model using the accelerated near-end gradient algorithm comprises:
initial value setting module 4: c i =0;
Wherein:
i represents the iteration number, and the initial value is 1;
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
And a judging module 4: judgment of
If yes, executing a judgment module 3, if not, making C i =C i +1, and returning to the D updating module 4.
Optionally, T is greater than or equal to 180.
Optionally, the T different imaging angles include: initial imaging angle of theta 0 And T imaging angles with an angle interval of delta theta.
It should be noted that the embodiments of the X-ray multi-spectrum separation imaging apparatus according to the present invention are the same in principle as the embodiments of the X-ray multi-spectrum separation imaging method, and the relevant points can be referred to each other.
The above description is only exemplary of the present invention and should not be taken as limiting the scope of the present invention, and any modifications, equivalents, improvements and the like made within the spirit and principle of the present invention should be included in the scope of the present invention.
Claims (13)
1. An X-ray multispectral separation imaging method, the method comprising:
step 11: setting imaging voltage range [ V1, V2] of measured object]When V is equal to [ V1, V2]]When the energy range meets the first preset condition, the penetration rate of the measured object meets the first preset condition, and the relative deviation of the linear attenuation coefficients of any two materials of the measured object at the V time meets the second preset condition; wherein the first preset condition is that: selecting a voltage interval [ V ] according to the material and size of the object to be measured, the imaging dynamic range of the detector and avoiding information loss caused by underexposure and overexposure 1 ,V 2 ]Make the objectThe penetration rate of the composite material reaches 30 to 70 percent; the second preset condition is as follows: from the angle of component analysis of the measured object, the linear attenuation coefficient curve of the material is used to satisfy any voltage V epsilon [ V ∈ [ ] 1 ,V 2 ]Time arbitrary k 1 、k 2 The seed material satisfies the energy range omega generated by the energy range omega and the energy range V 1 The intersection of (a) occupies [ omega ] 1 More than 60 percent of the total amount of the components to ensure that each component has differentiability in the selected energy range;
step 12: taking N different scanning voltages within the range of [ V1, V2], wherein N is larger than the total number K of materials contained in the measured object, and the difference between any two scanning voltages is larger than a first preset value, carrying out full-angle scanning imaging on the measured object under each scanning voltage to obtain N projection image sequences meeting the imaging gray scale requirement, wherein the full-angle scanning comprises T different imaging angles;
step 13: based on the N projection image sequences, solving projection thickness matrixes D corresponding to all materials of the measured object in a variable energy multispectral attenuation expression meeting the Poisson maximum likelihood principle;
step 14: calculating a separate projection p for each narrow spectral band in the X-ray based on the D r R =1,2 … R, R being the total number of narrow spectral segments of the X-ray;
step 15: reconstructing said p r Obtaining a reconstructed image a of each narrow energy spectrum section r 。
2. The method of claim 1, wherein the step 13 comprises:
converting the N projection image sequences into F = (F) n,m ) N×M Inputting a first decomposition model, wherein M is the total number of subunits of a projection image sequence of each scanning voltage, and calculating a projection thickness matrix D corresponding to all materials of the measured object meeting the first decomposition model;
the first decomposition model is
Wherein:
S=(s n,r ) N×R is a spectral decomposition coefficient matrix;
D=(d k,m ) N×M is a projected thickness matrix;
μ k (E r ) The k material is in E r Linear attenuation coefficient of time, E r Is the photon energy of the r-th narrow band.
3. The method of claim 2, wherein the calculating the projected thickness matrix D corresponding to all materials of the measured object that satisfy the first decomposition model comprises:
step 22: fixing the D, and updating S meeting the first decomposition model by using a non-negative matrix decomposition method;
step 23: fixing the S, and updating the D meeting the first decomposition model by using an accelerated near-end gradient algorithm;
step 24: calculating a difference value before and after the D value is updated, and judging whether the difference value meets a third preset condition or not; if yes, the updated D is output, and if no, the step 22 is returned.
5. The method according to claim 1, wherein the step 13 comprises:
converting the N projection image sequences into F = (F) n,m ) N×M Inputting a second decomposition model, wherein M is the total number of subunits of a projection image sequence of each scanning voltage, and calculating a projection thickness matrix D corresponding to all materials of the measured object meeting the second decomposition model;
Wherein:
S=(s n,r ) N×R is a spectral decomposition coefficient matrix;
D=(d k,m ) K×M is a projected thickness matrix;
μ k (E r ) The k material is in E r Linear attenuation coefficient of time, E r Is the photon energy of the r narrow band;
alpha is a regularization parameter;
6. The method of claim 5, wherein the calculating the projected thickness matrix D corresponding to all the materials of the measured object satisfying the second decomposition model comprises:
step 32: fixing the D, and updating the S meeting the second decomposition model by using a non-negative matrix decomposition method;
step 33: fixing the S, and updating the D meeting the second decomposition model by using an accelerated near-end gradient algorithm;
step 34: calculating a difference value before and after the D value is updated, and judging whether the difference value meets a third preset condition or not; if yes, the updated D is output, and if no, the step 32 is returned.
8. The method of claim 3, wherein the updating D satisfying the first decomposition model using an accelerated near-end gradient algorithm comprises:
step 41: c i =0;
Wherein:
i represents the iteration number, and the initial value is 1;
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
step 43: judgment of
9. The method of claim 6, wherein the updating D satisfying the second decomposition model using an accelerated near-end gradient algorithm comprises:
step 51: c i =0;
Wherein:
the fourth preset condition is that the search step length of the (i-eta) th time to the (i-1) th time is not changed;
step 53: judgment of
10. The method of claim 1, wherein T is equal to or greater than 180.
11. The method of claim 10, wherein the T different imaging angles comprise: initial imaging angle of theta 0 And T imaging angles with an angle interval of delta theta.
12. A non-transitory computer readable storage medium storing instructions which, when executed by a processor, cause the processor to perform the steps in the X-ray multispectral separation imaging method of any one of claims 1 to 11.
13. An X-ray multispectral separation imaging apparatus comprising a memory, a processor and a computer program stored in the memory and executable on the processor, characterized in that the processor implements the steps in the X-ray multispectral separation imaging method as claimed in any one of claims 1 to 11 when executing the computer program.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910107265.4A CN109875593B (en) | 2019-02-02 | 2019-02-02 | X-ray multi-spectrum separation imaging method, storage medium and device |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910107265.4A CN109875593B (en) | 2019-02-02 | 2019-02-02 | X-ray multi-spectrum separation imaging method, storage medium and device |
Publications (2)
Publication Number | Publication Date |
---|---|
CN109875593A CN109875593A (en) | 2019-06-14 |
CN109875593B true CN109875593B (en) | 2023-03-21 |
Family
ID=66927867
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910107265.4A Active CN109875593B (en) | 2019-02-02 | 2019-02-02 | X-ray multi-spectrum separation imaging method, storage medium and device |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN109875593B (en) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111134709B (en) * | 2020-01-17 | 2021-09-14 | 清华大学 | Multi-energy CT-based material decomposition method |
CN112697821B (en) * | 2020-12-02 | 2022-12-02 | 赛诺威盛科技(北京)股份有限公司 | Multi-energy spectrum CT scanning method and device, electronic equipment and CT equipment |
CN113392509B (en) * | 2021-05-27 | 2022-08-23 | 南方医科大学 | Accurate estimation method for actual energy spectrum of X-ray |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2013063577A1 (en) * | 2011-10-28 | 2013-05-02 | The University Of Chicago | Color x-ray histology for multi-stained biologic sample |
CN103876772A (en) * | 2014-03-20 | 2014-06-25 | 中北大学 | Multi-spectrum imaging method and device |
CN106483152A (en) * | 2016-11-21 | 2017-03-08 | 中北大学 | A kind of X-ray energy spectrum imaging method |
CN108010099A (en) * | 2017-12-04 | 2018-05-08 | 首都师范大学 | A kind of limited angle sweep of X-ray multi-power spectrum CT and image iterative reconstruction method |
Family Cites Families (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6950492B2 (en) * | 2003-06-25 | 2005-09-27 | Besson Guy M | Dynamic multi-spectral X-ray projection imaging |
US7760848B2 (en) * | 2006-09-08 | 2010-07-20 | General Electric Company | Method and system for generating a multi-spectral image of an object |
US8031831B2 (en) * | 2009-05-28 | 2011-10-04 | Kabushiki Kaisha Toshiba | Voltage and or current modulation in dual energy computed tomography |
US8862206B2 (en) * | 2009-11-12 | 2014-10-14 | Virginia Tech Intellectual Properties, Inc. | Extended interior methods and systems for spectral, optical, and photoacoustic imaging |
KR20150043630A (en) * | 2013-10-14 | 2015-04-23 | 삼성전자주식회사 | X-ray image apparatus and control method for the same |
US9672638B2 (en) * | 2014-06-16 | 2017-06-06 | The University Of Chicago | Spectral X-ray computed tomography reconstruction using a vectorial total variation |
US9870628B2 (en) * | 2015-03-18 | 2018-01-16 | Prismatic Sensors Ab | Image reconstruction based on energy-resolved image data from a photon-counting multi-bin detector |
-
2019
- 2019-02-02 CN CN201910107265.4A patent/CN109875593B/en active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2013063577A1 (en) * | 2011-10-28 | 2013-05-02 | The University Of Chicago | Color x-ray histology for multi-stained biologic sample |
CN103876772A (en) * | 2014-03-20 | 2014-06-25 | 中北大学 | Multi-spectrum imaging method and device |
CN106483152A (en) * | 2016-11-21 | 2017-03-08 | 中北大学 | A kind of X-ray energy spectrum imaging method |
CN108010099A (en) * | 2017-12-04 | 2018-05-08 | 首都师范大学 | A kind of limited angle sweep of X-ray multi-power spectrum CT and image iterative reconstruction method |
Non-Patent Citations (7)
Title |
---|
Material decomposition with the multi-energy attenuation coefficient ratio by using a multiple discriminant analysis;Woo-Jin Lee;《Journal of the Korean Physical Society》;20160727(第7期);231-240 * |
Narrow-Energy-Width CT Based on Multivoltage X-Ray Image Decomposition;Jiaotong Wei等;《Int J Biomed Imaging》;20171107(第11期);1-9 * |
基于MARS系统的X射线能谱CT研究;何鹏;《中国优秀硕士论文全文数据库信息科技辑》;20140215(第2期);I138-22 * |
基于变电压图像序列盲分离的X射线多谱CT成像;魏交统;《中国优秀博士论文全文数据库信息科技辑》;20170315(第3期);I138-68 * |
基于投影替代与矩阵低秩稀疏分解的多光谱图像融合;戎凯旋;《中国优秀硕士论文全文数据库信息科技辑》;20170215(第2期);I140-159 * |
基于能谱滤波分离的彩色CT成像方法研究;牛素鋆;《中国优秀硕士论文全文数据库信息科技辑》;20150715(第7期);I138-1268 * |
多电压X射线图像分解的高对比度CT成像;魏交统等;《光电工程》;20160831;第43卷(第8期);59-63 * |
Also Published As
Publication number | Publication date |
---|---|
CN109875593A (en) | 2019-06-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109875593B (en) | X-ray multi-spectrum separation imaging method, storage medium and device | |
CN108010098B (en) | Double-energy-spectrum CT (computed tomography) base material image iterative reconstruction method | |
EP2843623B1 (en) | X-ray dual-energy CT reconstruction method | |
US9406154B2 (en) | Iterative reconstruction in image formation | |
Ding et al. | Image‐domain multimaterial decomposition for dual‐energy CT based on prior information of material images | |
JP6630841B2 (en) | X-ray detection method and X-ray detector | |
WO2014171539A1 (en) | X-ray computed tomography device and correction method | |
JP6793469B2 (en) | Data processing equipment, X-ray CT equipment and data processing method | |
Zhang et al. | Constrained total generalized p-variation minimization for few-view X-ray computed tomography image reconstruction | |
Bubba et al. | Tomographic X-ray data of carved cheese | |
Cuadros et al. | Compressive spectral X-ray tomography based on spatial and spectral coded illumination | |
CN110415307A (en) | A kind of multipotency CT imaging method based on tensor completion, device and its storage equipment | |
Zhao et al. | An oblique projection modification technique (OPMT) for fast multispectral CT reconstruction | |
Kingston et al. | Inherent dose-reduction potential of classical ghost imaging | |
Kadu et al. | A convex formulation for binary tomography | |
CN109916933B (en) | X-ray computed tomography energy spectrum estimation method based on convolutional neural network | |
Titarenko | 1-D filter for ring artifact suppression | |
US9330456B2 (en) | Systems and methods for regularized Fourier analysis in x-ray phase contrast imaging | |
Bobin et al. | Sparse BSS From Poisson Measurements | |
Odinaka et al. | Spectrally grouped total variation reconstruction for scatter imaging using ADMM | |
CN109448071A (en) | A kind of power spectrum image rebuilding method and system | |
Wu et al. | Iterative CT reconstruction using curvelet-based regularization | |
Ramyar et al. | Establishing a method to measure bone structure using spectral CT | |
Brito et al. | Application of neural networks to model the signal-dependent noise of a digital breast tomosynthesis unit | |
Chen et al. | Ring artifacts reduction in cone-beam CT images based on independent component analysis |
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 | ||
GR01 | Patent grant | ||
GR01 | Patent grant |