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 PDF

Info

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
Application number
CN201910107265.4A
Other languages
Chinese (zh)
Other versions
CN109875593A (en
Inventor
魏交统
赵晓杰
陈平
刘宾
潘晋孝
韩焱
苏新彦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
North University of China
Original Assignee
North University of China
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by North University of China filed Critical North University of China
Priority to CN201910107265.4A priority Critical patent/CN109875593B/en
Publication of CN109875593A publication Critical patent/CN109875593A/en
Application granted granted Critical
Publication of CN109875593B publication Critical patent/CN109875593B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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

X-ray multi-spectrum separation imaging method, storage medium and device
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.
Figure GDA0002115178310000041
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 ：
Figure GDA0002115178310000042
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.
Figure GDA0002115178310000051
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
Figure GDA0002115178310000052
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 intervals
Figure GDA0002115178310000053
Energy 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
Figure GDA0002115178310000061
Wherein
Figure GDA0002115178310000062
Called spectral decomposition coefficient, satisfies
Figure GDA0002115178310000063
It is usually unknown, and
Figure GDA0002115178310000064
i.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.
Figure GDA0002115178310000065
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.
Note the book
Figure GDA0002115178310000066
The ideal ratio of the ray intensity to the incident intensity is obtained.
Figure GDA0002115178310000067
Step 13 is according to the known F = (F) n,m ) N×M Solving variable energy multispectral attenuation expression satisfying Poisson maximum likelihood principle
Figure GDA0002115178310000068
And (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 value
Figure GDA0002115178310000071
The likelihood function of the difference is in the form:
Figure GDA0002115178310000072
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:
Figure GDA0002115178310000073
wherein:
Figure GDA0002115178310000074
Figure GDA0002115178310000075
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 principle
Figure GDA0002115178310000076
And 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:
step 21: setting initial values of S and D:
Figure GDA0002115178310000077
Figure GDA0002115178310000078
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 is
Figure GDA0002115178310000081
Or
Figure GDA0002115178310000082
Or | | | 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:
Figure GDA0002115178310000083
wherein:
Figure GDA0002115178310000084
Figure GDA0002115178310000085
Figure GDA0002115178310000086
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；
Step 42: updating
Figure GDA0002115178310000087
Wherein:
Figure GDA0002115178310000088
Figure GDA0002115178310000091
t i in order to search for the step size,
Figure GDA0002115178310000092
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
Figure GDA0002115178310000093
Figure GDA0002115178310000094
step 43: judgment of
Figure GDA0002115178310000095
If true, go to step 24, if true, let C go to i ＝C i +1, return to step 42.
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
Figure GDA0002115178310000096
In the formula D k ＝(d k,m ) M×1 Q denotes a two-or one-dimensional discrete gradient operator,
Figure GDA0002115178310000097
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.
Figure GDA0002115178310000098
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;
the second decomposition model is
Figure GDA0002115178310000101
Wherein:
Figure GDA0002115178310000102
Figure GDA0002115178310000103
Figure GDA0002115178310000104
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 ；
Figure GDA0002115178310000105
D k ＝(d k,m ) M×1 Q is a two-dimensional or one-dimensional discrete gradient operator;
Figure GDA0002115178310000106
is a Tikhonov (L2) regularization function, beta k Is a regularization parameter.
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 (ii) of31: setting initial values of S and D:
Figure GDA0002115178310000107
Figure GDA0002115178310000108
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 is
Figure GDA0002115178310000111
Or
Figure GDA0002115178310000112
Or | | | 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:
Figure GDA0002115178310000113
wherein:
Figure GDA0002115178310000114
Figure GDA0002115178310000115
Figure GDA0002115178310000116
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；
Step 52: updating
Figure GDA0002115178310000117
Wherein:
Figure GDA0002115178310000118
Figure GDA0002115178310000119
t i in order to search for the step size,
Figure GDA00021151783100001110
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
Figure GDA00021151783100001111
Figure GDA0002115178310000121
Figure GDA0002115178310000122
is the m-th row vector of the transposed matrix of the gradient operator Q.
Step 53: judgment of
Figure GDA0002115178310000123
If true, go to step 34, if false, let C i ＝C i +1, return to step 52.
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;
the first decomposition model is
Figure GDA0002115178310000131
Wherein:
Figure GDA0002115178310000132
Figure GDA0002115178310000133
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:
initial value setting module 1: setting initial values of S and D:
Figure GDA0002115178310000134
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:
Figure GDA0002115178310000141
wherein:
Figure GDA0002115178310000142
Figure GDA0002115178310000143
Figure GDA0002115178310000144
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；
D, updating the module 2: updating
Figure GDA0002115178310000145
Wherein:
Figure GDA0002115178310000146
Figure GDA0002115178310000147
i represents the iteration number, and the initial value is 1;
Figure GDA0002115178310000148
Figure GDA0002115178310000149
t i in order to search for the step size,
Figure GDA0002115178310000151
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
Figure GDA0002115178310000152
Figure GDA0002115178310000153
and a judging module 2: judgment of
Figure GDA0002115178310000154
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;
the second decomposition model is
Figure GDA0002115178310000155
Wherein:
Figure GDA0002115178310000156
Figure GDA0002115178310000157
Figure GDA0002115178310000158
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;
Figure GDA0002115178310000159
D k ＝(d k,m ) M×1 q is a two-dimensional or one-dimensional discrete gradient operator;
Figure GDA00021151783100001510
is a Tikhonov (L2) regularization function, beta k 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:
initial value setting module 3: setting initial values of S and D:
Figure GDA0002115178310000161
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:
Figure GDA0002115178310000162
wherein
Figure GDA0002115178310000163
Figure GDA0002115178310000164
Figure GDA0002115178310000165
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；
D, updating the module 4: updating
Figure GDA0002115178310000166
Wherein:
Figure GDA0002115178310000167
Figure GDA0002115178310000168
i represents the iteration number, and the initial value is 1;
Figure GDA0002115178310000171
Figure GDA0002115178310000172
t i in order to search for the step size,
Figure GDA0002115178310000173
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
Figure GDA0002115178310000174
Figure GDA0002115178310000175
Figure GDA0002115178310000176
is the m-th row vector of the transposed matrix of the gradient operator Q.
And a judging module 4: judgment of
Figure GDA0002115178310000177
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
Figure FDA0004047147300000021
Wherein:
Figure FDA0004047147300000022
Figure FDA0004047147300000023
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 21: setting initial values of the S and the D:
Figure FDA0004047147300000024
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.
4. The method of claim 3, wherein updating S satisfying the first decomposition model using a non-negative matrix factorization method comprises:
Figure FDA0004047147300000025
wherein:
Figure FDA0004047147300000026
Figure FDA0004047147300000031
Figure FDA0004047147300000032
i represents the number of iterations, and the initial value is 1.
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;
the second decomposition model is
Figure FDA0004047147300000033
Wherein:
Figure FDA0004047147300000034
Figure FDA0004047147300000035
Figure FDA0004047147300000036
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;
Figure FDA0004047147300000037
D k ＝(d k，m ) M×1 q is a two-dimensional or one-dimensional discrete gradient operator;
Figure FDA0004047147300000038
is a Tikhonov (L2) regularization function, beta k 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 31: setting initial values of the S and the D:
Figure FDA0004047147300000039
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.
7. The method of claim 6, wherein the updating S satisfying the second decomposition model using a non-negative matrix factorization method comprises:
Figure FDA0004047147300000041
wherein:
Figure FDA0004047147300000042
Figure FDA0004047147300000043
Figure FDA0004047147300000044
i represents the number of iterations, and the initial value is 1.
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；
Step 42: updating
Figure FDA0004047147300000045
Wherein:
Figure FDA0004047147300000046
Figure FDA0004047147300000047
i represents the iteration number, and the initial value is 1;
Figure FDA0004047147300000048
Figure FDA0004047147300000049
t i in order to search for the step size,
Figure FDA0004047147300000051
ξ∈(0，1)，
Figure FDA0004047147300000052
the fourth preset condition is that the search step length from (i-eta) to (i-1) is unchanged;
Figure FDA0004047147300000053
Figure FDA0004047147300000054
step 43: judgment of
Figure FDA0004047147300000055
If true, go to step 24, if true, let C go to i ＝C i +1, return to step 42.
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；
Step 52: updating
Figure FDA0004047147300000056
Wherein:
Figure FDA0004047147300000057
Figure FDA0004047147300000058
i represents the iteration number, and the initial value is 1;
Figure FDA0004047147300000059
Figure FDA00040471473000000510
t i in order to search for the step size,
Figure FDA00040471473000000511
ξ∈(0，1)，
Figure FDA00040471473000000512
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;
Figure FDA0004047147300000061
Figure FDA0004047147300000062
Figure FDA0004047147300000063
is the m-th row vector of the transposed matrix of the gradient operator Q;
step 53: judgment of
Figure FDA0004047147300000064
If true, go to step 34, if false, let C i ＝C i +1, return to step 52.
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.
CN201910107265.4A 2019-02-02 2019-02-02 X-ray multi-spectrum separation imaging method, storage medium and device Active CN109875593B (en)

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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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
JP6585182B2 (en) * 2015-03-18 2019-10-02 プリズマティック、センサーズ、アクチボラグＰｒｉｓｍａｔｉｃ Ｓｅｎｓｏｒｓ Ａｂ Image reconstruction based on energy-resolved image data from a photon counting multi-bin detector

Patent Citations (4)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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
Liu et al. Total variation-stokes strategy for sparse-view X-ray CT image reconstruction
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
Barutcu et al. Limited-angle computed tomography with deep image and physics priors
Bubba et al. Tomographic X-ray data of carved cheese
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
CN106157345B (en) Method for generating an image
Odinaka et al. Spectrally grouped total variation reconstruction for scatter imaging using ADMM
Yazdanpanah et al. Sparse-view CT reconstruction using curvelet and TV-based regularization
Wu et al. Iterative CT reconstruction using curvelet-based regularization
Bobin et al. Sparse BSS From Poisson Measurements
Zeegers et al. A multi-channel dart algorithm
Chen et al. Ring artifacts reduction in cone-beam CT images based on independent component analysis
Melli et al. A sparsity-based iterative algorithm for reconstruction of micro-CT images from highly undersampled projection datasets obtained with a synchrotron X-ray source
Pineda et al. What does DQE say about lesion detectability in digital radiography?

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