CN109875593A - CT imaging method, storage medium and device - Google Patents
CT imaging method, storage medium and device Download PDFInfo
- Publication number
- CN109875593A CN109875593A CN201910107265.4A CN201910107265A CN109875593A CN 109875593 A CN109875593 A CN 109875593A CN 201910107265 A CN201910107265 A CN 201910107265A CN 109875593 A CN109875593 A CN 109875593A
- Authority
- CN
- China
- Prior art keywords
- imaging
- meeting
- decomposition model
- measured object
- matrix
- 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.)
- Granted
Links
Landscapes
- Analysing Materials By The Use Of Radiation (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
The present invention provides a kind of CT imaging method, storage medium and device, this method comprises: step 11: setting the imaging voltage range [V1, V2] of measured object；Step 12: in [V1, V2] in the range of take N number of different scanning voltage, N is greater than the material sum K that measured object is included, and the difference between any two scanning voltage is greater than the first preset value, full angle scanning imagery is carried out to measured object under each scanning voltage, obtains and meets N number of projection image sequence that imaging gray scale requires；Step 13: being based on N number of projection image sequence, solve the corresponding projection thickness matrix D of all materials for becoming the measured object in the multispectral decaying expression formula of energy for meeting Poisson maximum likelihood principle；Step 14: being based on D, the separation for calculating each narrow energy spectral coverage in X-ray projects pr, r=1,2 ... R；Step 15: rebuilding prObtain the reconstruction image a of each narrow energy spectral coverager.The multi-power spectrum CT imaging of low cost, power spectrum high resolution may be implemented in method of the invention.
Description
Technical field
The present invention relates to instrument field, in particular to a kind of CT imaging method, storage medium and device.
Background technique
During x-ray imaging, since the multispectral characteristic of X-ray and the attenuation coefficient of substance are usually with photon energy
Increase and reduce, ray when passing through object its spectral centroid with the increase for passing through thickness it is gradually mobile to high energy direction,
Generate hardening effect.Hardening effect destroys the linear relationship between attenuation coefficient and data for projection, so that conventional reconstruction method
Gained CT image generates hardening artifact, " with shadow foreign matter " and " the different shadow of jljl " phenomenon occurs, is unfavorable for quantitative analysis, such as novel
The quantification of the microstructures such as constituent content, the ore porosity of material characterizes.
In view of the above-mentioned problems, there is the correlative study of multispectral CT imaging technique.Typical dual intensity CT is by obtaining two
The data for projection of different power spectrums can be realized the non-linear reconstruction to data for projection, it mainly include projection dual intensity decompose and
Two steps of image reconstruction, the more commonly used algorithm for reconstructing of double-energy CT system common at present is based on substratess matter decomposition model
Double substance decomposition algorithms, due to introducing the information of power spectrum, dual intensity CT has better substance separating capacity compared with traditional CT.But dual intensity
Only there are two power spectrums by CT, and discrimination is not high, and then there has been proposed multispectral CT imaging techniques.Wherein, application and research are more
It is the multi-power spectrum CT imaging technique based on photon counting-type detector and based on synchrotron radiation source, photon counting-type detector can be with
Multiple narrow spectral imagings are realized, the disadvantage is that its time/spatial resolution is lower, technical costs is higher；Synchrotron radiation source has good
Good monochromaticjty, but its equipment is huge, involves great expense, it is difficult to it applies in practical projects.
How low cost is realized, multi-power spectrum (multispectral) CT imaging technique of power spectrum high resolution is skill urgently to be resolved at present
Art problem.
Summary of the invention
In view of this, the present invention provides a kind of CT imaging method, storage medium and device, with solve how to realize it is low at
Originally, the problem of the multi-power spectrum CT imaging technique of power spectrum high resolution.
The present invention provides a kind of CT imaging method, this method comprises:
Step 11: setting the imaging voltage range [V1, V2] of measured object, as V ∈ [V1, V2], the penetrance of measured object
Meet linear attenuation coefficient of any two kinds of materials of the first preset condition and measured object in V relative deviation meet it is default
It is required that corresponding energy range meets the second preset condition；
Step 12: N number of different scanning voltage is taken in the range of [V1, V2], it is total that N is greater than the material that measured object is included
Number K, and the difference between any two scanning voltage is greater than the first preset value, carries out under each scanning voltage to measured object complete
Angle scanning imaging obtains and meets N number of projection image sequence that imaging gray scale requires, full angle scanning include T it is different at
Image angle degree；
Step 13: being based on N number of projection image sequence, solve the multispectral attenuation meter of change energy for meeting Poisson maximum likelihood principle
Up to the corresponding projection thickness matrix D of all materials of the measured object in formula；
Step 14: being based on D, the separation for calculating each narrow energy spectral coverage in X-ray projects pr, r=1,2 ... R, R are X-ray
Narrow energy spectral coverage sum；
Step 15: rebuilding prObtain the reconstruction image a of each narrow energy spectral coverager。
The present invention also provides a kind of non-transitory computer-readable storage medium, non-transitory computer-readable storage medium storages
Instruction, instruction make processor execute the step in above-mentioned CT imaging method when executed by the processor.
The present invention also provides a kind of CT imaging device, including memory, processor and storage are in memory and can be
The computer program run on processor, processor realize the step in above-mentioned CT imaging method when executing computer program.
X-ray projection image analogy may be implemented in CT imaging method of the invention, according to the material of object, size, choosing
Reasonable imaging parameters are taken to obtain x-ray projection image sequence；Then according to projection image sequence, solution meets Poisson greatly seemingly
The corresponding projection thickness matrix D of all materials for becoming the measured object in the multispectral decaying expression formula of energy of right principle；Further according to square
Battle array D realizes the projection separation of multispectral projected image.Skill is imaged compared to the X ray CT rebuild based on hardware and other multi-power spectrums
Art, this method from projection domain realize multispectral radioscopic image power spectrum separation, without to existing high dynamic CT imaging device into
Row improves, and has many advantages, such as at low cost, power spectrum high resolution, and the application for promoting quantification CT imaging technique has important meaning
Justice.
Detailed description of the invention
Fig. 1 is the flow chart of CT imaging method of the present invention；
Fig. 2 is the projection image sequence schematic diagram that a scanning voltage obtains；
The subelement that Fig. 3 is Fig. 2 segments schematic diagram；
Fig. 4 is the first embodiment of the solution D of step 13 in Fig. 1；
Fig. 5 is the second embodiment of the solution D of step 13 in Fig. 1；
Fig. 6 is the structure chart of CT imaging device of the present invention.
Specific embodiment
To make the objectives, technical solutions, and advantages of the present invention clearer, right in the following with reference to the drawings and specific embodiments
The present invention is described in detail.
As shown in Figure 1, the present invention provides a kind of CT imaging method, comprising:
Step 11 (S11): the imaging voltage range [V1, V2] of measured object is set, as V ∈ [V1, V2], measured object is worn
The relative deviation that saturating rate meets linear attenuation coefficient of any two kinds of materials of the first preset condition and measured object in V meets
The corresponding energy range of preset requirement meets the second preset condition；
For example, the first preset condition: according to material, size, the detector image-forming dynamic range of measured object, and avoiding owing exposure
Light and overexposure phenomenon cause the missing of information, choose suitable voltage range [V1,V2] make the penetrance of object reach 30%~
70% or other ranges.
Second preset condition are as follows: from the angle of measured object component analysis, according to the linear attenuation coefficient curve of material,
Meet free voltage V ∈ [V1,V2] Shi Renyi kth1、k2Kind material meets energy caused by the energy range Ω and V of formula (1)
Range Ω1Intersection account for Ω160% or more, with guarantee each component choose energy range on have ga s safety degree.
μ in formulak1(E) and μk2(E) kth is respectively indicated1、k2Linear attenuation coefficient of the kind material in ENERGY E.According to default
It is required that it is stringent whether, 0.05 in formula (1) can also change other values into.
Step 12 (S12): taking N number of different scanning voltage in the range of [V1, V2], and N is greater than measured object and is included
Material sum K, and the difference between any two scanning voltage is greater than the first preset value, to measured object under each scanning voltage
Full angle scanning imagery is carried out, acquisition meets N number of projection image sequence that imaging gray scale requires, and full angle scanning includes T a not
Same imaging angle；
N number of different scanning voltage can be contact potential series at equal intervals, or unequal interval contact potential series, either
No all to need to meet at equal intervals, N is greater than the material sum K that measured object is included, and the difference between any two scanning voltage is big
In the first preset value, the first preset value is according to transmitted intensity in the same ray path for passing through measured object and incident intensity ratio
Minimum distinguishes requirement setting, such as minimum distinguish requires to be greater than 2% for the relative deviation of ratio.
Similarly, T imaging angle may be that at equal intervals or unequal interval, such as T different imaging angles include:
Initial imaging angle is θ0, the T imaging angle of △ θ is divided between angle.The value of T is determined that usual T >=180 are equal by step 13
The requirement that step 13 solves D can be met.
When each angle imaging of each scanning voltage, by adjusting the electric current of X-ray tube and the time of integration of detector
Make collected radioscopic image gray value in effective areas imaging of detector, i.e. imaging gray scale meets the requirements.
Step 13 (S13): being based on N number of projection image sequence, and the change energy that solution meets Poisson maximum likelihood principle is multispectral
The corresponding projection thickness matrix D of all materials of measured object in the expression formula that decays；
Step 14 (S14): being based on D, and the separation for calculating each narrow energy spectral coverage in X-ray projects pr, r=1,2 ... R, R X
The narrow energy spectral coverage sum of ray；
The matrix D that step 13 obtains=(dk,m)K×M, wherein M is the subelement of the projection image sequence of each scanning voltage
Sum can must separate projection p by formula (2)r:
Step 15 (S15): p is rebuildrObtain the reconstruction image a of each narrow energy spectral coverager。
Such as by prProgress general linear reconstruction (Yan Bin, Li Lei .CT image reconstruction algorithm [M] Science Press,
2014) the reconstruction image a of each narrow energy spectral coverage is obtainedr。
The principle explanation of step 13 is given below.
The N number of projection image sequence for obtaining step 12 is needed to be converted into F=(f firstn,m)N×M；M is each scanning voltage
The subelement sum of projection image sequence, due to the corresponding ray path of each subelement, M is it also will be understood that each scanning is electric
Depress the ray path sum of X-ray.
Fig. 2 is the projection image sequence that obtains under some scanning voltage, will wherein each image all as shown in Figure 3 with image
Pixel establishes index as subelement, and for each subelement, then M=T × n1 × m1.
It by the projection image sequence obtained under each scanning voltage using image pixel as subelement, and is every height list
Member establishes index, so that it may it is as follows to obtain F matrix.
F is the matrix of a N × M.Wherein fn,mFor penetrating on m paths in n-th of scanning voltage projection image sequence
Line intensity (gray value of the unit in correspondence image) and incident intensity ratio.Incident intensity refers generally to not pass through object area
Ray path on gray value.
During x-ray imaging, the multispectral characteristic of X-ray make include in detector projecting image data obtained
Multispectral section of decaying, this for using multispectral projected image carry out Projective decomposition provide foundation.
If the intensity on x-ray path L is I, initial strength I0, normalization spectrum is S (E), ceiling capacity Emax,
Object attenuation coefficient at the x of spatial position is μ (x, E), is had according to multispectral decaying expression formula
Enable Emax=max { E1,max,…,En,max,…,EN,max, En,maxIndicate the maximum photon energy of n-th of scanning voltage
Amount.According to actual conditions by energy section [0, Emax] it is divided into R narrow energy spectral coverage subinterval [E'r-1,E'r], r=1,2 ... R,
Wherein 0=E'0<E'1<…<E'R=Emax, usually take at equal intervalsSection self-energy Er=t
E'r+(1-t)·E'r-1, generally take t=0.5.
If imaging substance is made of K kind material, it is assumed that energy section [E'r-1,E'r] sufficiently narrow so that kth kind material is at this
Attenuation coefficient on section is constant for μk(Er), ErIt takes after determining, the linear attenuation coefficient μ of each materialk(Er) referring to U.S. NIST net
Standing can obtain.
Enable on the L of path the kth kind material of measured object with a thickness of dk, then formula (3) can be converted into
WhereinReferred to as spectral resolution coefficient meetsIt is usually unknown, andDecomposition projection as to be asked.
And dkRelated with ray path, i.e., the thickness of every kind of material is different in different ray paths, it is desirable that obtains single-pathway
The multispectral projection that the thickness of upper various material only leans under single sweep voltage is serious underdetermined problem, so needing to acquire multiple scannings
The projected image of different-energy under voltage.Simultaneously because spectral resolution coefficient is unknown, therefore generally require acquisition scans voltage amounts N
The number K of material included greater than substance.
It altogether include M ray path, f in the projection image sequence of each scanning voltagen,mFor the throwing of n-th of scanning voltage
Transmitted intensity and incident intensity ratio in shadow image sequence on m paths, sn,rIt is generated the on X-ray spectrum for n-th of scanning
R spectral resolution coefficient, dk,mFor the thickness of kth kind material on m paths.Assuming that transmitted intensity in each path with enter
Penetrating intensity rate is the Poisson stochastic variable being independently distributed, i.e.,
Enable S=(sn,r)N×RFor spectral resolution coefficient matrix, R is the narrow energy spectral coverage sum of X-ray；D=(dk,m)K×MTo throw
Shadow thickness matrix.
NoteFor ideal transmitted intensity and incident intensity
Ratio.
Then step 13 is according to known F=(fn,m)N×M, solve and meet the change energy of Poisson maximum likelihood principle and multispectral decline
Subtract expression formulaIn measured object the corresponding projection thickness matrix D of all materials.
The first embodiment of step 13 is given below.
Cephalometry f multispectral for N × M Poissonn,m, weighed with the negative log-likelihood function of broad sense KL Divergence Form
Measure real data fn,mWith ideal valueBetween difference, likelihood function form is as follows:
Wherein, S, D are independent variable, and every a line of F is to be arranged successively with projection image sequence pixel distribution and includes difference
The data for projection of projection angle.Then change Energy X-ray the first decomposition model of projected image based on Poisson priori as follows is obtained:
Wherein:
S=(sn,r)N×RFor spectral resolution coefficient matrix, D=(dk,m)K×MFor projection thickness matrix.
μk(Er) it is kth kind material in ErWhen linear attenuation coefficient, ErFor the photon energy of r-th narrow energy spectral coverage.
Solve the multispectral decaying expression formula of change energy for meeting Poisson maximum likelihood principleMiddle D
And S, the D and S for solving and meeting the first decomposition model can also be converted to.
Thickness according to material in spectral resolution coefficient and ray path is nonnegative value, obtains solving variable S and D satisfaction
S >=0, each element s in the representing matrix of D >=0 are remembered in nonnegativity restrictionsn,r>=0, dk,m>=0, with S1R=1NIndicate each electricity
The Spectrum attenuation coefficient of pressure and be 1,1RIndicate that R ties up complete 1 column vector.
As shown in figure 4, calculating all materials of measured object for meeting the first decomposition model in the first embodiment of step 13
Corresponding projection thickness matrix D includes:
Step 21: the initial value of setting S and D:
En,maxFor maximum photon energy caused by n-th of scanning voltage.
Step 22: fixed D is updated the S for meeting the first decomposition model using non-negative matrix factorization method；
Step 23: fixed S is updated the D for meeting the first decomposition model using acceleration proximal end gradient algorithm；
Step 24: calculating the difference that D value updates front and back, judge whether difference meets third preset condition；If so, output
Updated D, if it is not, then return step 22.
For example, third preset condition isOrOr | | Di-Di-1||F< ε；ε value according to
Actual conditions are selected, such as ε=10-3。
Wherein, the S for meeting the first decomposition model is updated using non-negative matrix factorization method in step 22 and includes:
Wherein:
I indicates the number of iterations, initial value 1.
Wherein, the D for meeting the first decomposition model is updated using acceleration proximal end gradient algorithm in step 23 and includes:
Step 41:Ci=0；
Step 42: updating
Wherein:
tiFor step-size in search,ξ ∈ (0,1),
4th preset condition is that (i- η)~the (i-1) secondary step-size in search is constant；
Step 43: judgement
It is
No establishment, if so, step 24 is executed, if it is not, then enabling Ci=Ci+ 1, return step 42.
The second embodiment of step 13 is given below.
Remember fn,mAnd fn,m+1The transmitted intensity and incident intensity ratio on adjacent ray path are indicated, since compositions of matter is in
Continuity distribution, therefore in formula (7) on adjacent ray path every kind of material projection thickness value dk,mAnd dk,m+1Otherness is smaller, institute
Smooth priori can be added to projection thickness matrix, following penalty is established
D in formulak=(dk,m)M×1, Q indicate two dimension or one-dimensional discrete gradient operator,For
Tikhonov (L2) Regularization function, βkFor regularization parameter, slickness degree can be added for every kind of material using the parameter
Different priori.
It is as follows can to obtain the second decomposition model for increase formula (9) in formula (7), and α is regularization parameter in formula.
Based on this, step 13 can also include:
F=(f is converted by N number of projection image sequencen,m)N×MThe second decomposition model is inputted, M is the throwing of each scanning voltage
The subelement sum of shadow image sequence, calculates the corresponding projection thickness matrix of all materials of measured object for meeting the second decomposition model
D；
Second decomposition model is
Wherein:
S=(sn,r)N×RFor spectral resolution coefficient matrix；D=(dk,m)K×MFor projection thickness matrix；
μk(Er) it is kth kind material in ErWhen linear attenuation coefficient, ErFor the photon energy of r-th narrow energy spectral coverage；
α is regularization parameter, can be set as 10-7≤α≤10-3；
Q is two dimension or one-dimensional discrete gradient operator；
For Tikhonov (L2) Regularization function, βkFor regularization parameter.
As shown in figure 5, calculating all materials of measured object for meeting the second decomposition model in the second embodiment of step 13
Corresponding projection thickness matrix D includes:
Step 31: the initial value of setting S and D:
Step 32: fixed D is updated the S for meeting the second decomposition model using non-negative matrix factorization method；
Step 33: fixed S is updated the D for meeting the second decomposition model using acceleration proximal end gradient algorithm；
Step 34: calculating the difference that D value updates front and back, judge whether difference meets third preset condition；If so, output
Updated D, if it is not, then return step 32.
For example, third preset condition isOrOr | | Di-Di-1||=< ε；ε value according to
Actual conditions are selected, such as ε=10-3。
Wherein, the S for meeting the second decomposition model is updated using non-negative matrix factorization method in step 32 and includes:
Wherein:
I indicates the number of iterations, initial value 1.
Wherein, the D for meeting the first decomposition model is updated using acceleration proximal end gradient algorithm in step 33 and includes:
Step 51:Ci=0；
Step 52: updating
Wherein:
tiFor step-size in search,ξ ∈ (0,1),
4th preset condition is that (i- η)~the (i-1) secondary step-size in search is constant；
For the m row vector of the transposed matrix of gradient operator Q.
Step 53: judgement
Whether
It sets up, if so, step 34 is executed, if it is not, then enabling Ci=Ci+ 1, return step 52.
In addition, the solution of step 13 is also referred to the derivation of the basic method for solving of nonnegative matrix model, this method in Fig. 1
It is stronger to the dependence of initial value, and restrained slowly near minimum.
It is the explanation to CT imaging method of the present invention above.
X-ray projection image analogy may be implemented in CT imaging method of the invention, according to the material of object, size, choosing
Reasonable imaging parameters are taken to obtain x-ray projection image sequence；Then according to projection image sequence, solution meets Poisson greatly seemingly
The corresponding projection thickness matrix D of all materials for becoming the measured object in the multispectral decaying expression formula of energy of right principle；Further according to square
Battle array D realizes the projection separation of multispectral projected image.Skill is imaged compared to the X ray CT rebuild based on hardware and other multi-power spectrums
Art, this method from projection domain realize multispectral radioscopic image power spectrum separation, without to existing high dynamic CT imaging device into
Row improves, and has many advantages, such as at low cost, power spectrum high resolution, and the application for promoting quantification CT imaging technique has important meaning
Justice.
The present invention also provides a kind of non-transitory computer-readable storage medium, non-transitory computer-readable storage medium storages
Instruction, instruction make processor execute the step in above-mentioned CT imaging method when executed by the processor.
The present invention also provides a kind of CT imaging device, including memory, processor and storage are in memory and can be
The computer program run on processor, processor realize the step in above-mentioned CT imaging method when executing computer program.
As shown in fig. 6, the device includes:
Scanning voltage range setting module: setting the imaging voltage range [V1, V2] of measured object, as V ∈ [V1, V2],
The penetrance of measured object meets the phase of linear attenuation coefficient of any two kinds of materials of the first preset condition and measured object in V
The corresponding energy range of preset requirement is met to deviation and meets the second preset condition；
Become energy full scan module: taking N number of different scanning voltage in the range of [V1, V2], N is greater than measured object and is wrapped
The material sum K contained, and the difference between any two scanning voltage is greater than the first preset value, to quilt under each scanning voltage
It surveys object and carries out full angle scanning imagery, obtain and meet N number of projection image sequence that imaging gray scale requires, full angle scanning includes T
A different imaging angle；
D resolves module: being based on N number of projection image sequence, solves and meet the change energy of Poisson maximum likelihood principle and multispectral decline
Subtract the corresponding projection thickness matrix D of all materials of the measured object in expression formula；
Separation projection computing module: being based on D, and the separation for calculating each narrow energy spectral coverage in X-ray projects pr, r=1,2 ... R,
R is the narrow energy spectral coverage sum of X-ray；
It separates backprojection image reconstruction module: rebuilding prObtain the reconstruction image a of each narrow energy spectral coverager。
Wherein, D, which resolves module, to include:
F=(f is converted by N number of projection image sequencen,m)N×MThe first decomposition model is inputted, M is the throwing of each scanning voltage
The subelement sum of shadow image sequence, calculates the corresponding projection thickness matrix of all materials of measured object for meeting the first decomposition model
D；
First decomposition model is
Wherein:
S=(sn,r)N×RFor spectral resolution coefficient matrix；
D=(dk,m)K×MFor projection thickness matrix；
μk(Er) it is kth kind material in ErWhen linear attenuation coefficient, ErFor the photon energy of r-th narrow energy spectral coverage.
In turn, the corresponding projection thickness matrix D of all materials of measured object of the first decomposition model of calculating satisfaction includes:
Initial value design module 1: the initial value of setting S and D:
S update module 1: fixed D is updated the S for meeting the first decomposition model using non-negative matrix factorization method；
D update module 1: fixed S is updated the D for meeting the first decomposition model using acceleration proximal end gradient algorithm；
Judgment module 1: the difference that D value updates front and back is calculated, judges whether difference meets third preset condition；If so,
Updated D is exported, if it is not, then returning to S update module 1.
In S update module 1: being updated to the S for meeting the first decomposition model using non-negative matrix factorization method and include:
Wherein:
I indicates the number of iterations, initial value 1.
In D update module 1: including: using accelerating proximal end gradient algorithm to be updated the D for meeting the first decomposition model
Initial value design module 2:Ci=0；
D update module 2: it updates
Wherein:
I indicates the number of iterations, initial value 1；
tiFor step-size in search,ξ ∈ (0,1),
4th preset condition is that (i- η)~the (i-1) secondary step-size in search is constant；
Judgment module 2: judgement
It is
No establishment, if so, judgment module 1 is executed, if it is not, then enabling Ci=Ci+ 1, return to D update module 2.
Wherein, D, which resolves module, to include:
F=(f is converted by N number of projection image sequencen,m)N×MThe second decomposition model is inputted, M is the throwing of each scanning voltage
The subelement sum of shadow image sequence, calculates the corresponding projection thickness matrix of all materials of measured object for meeting the second decomposition model
D；
Second decomposition model is
Wherein:
S=(sn,r)N×RFor spectral resolution coefficient matrix；
D=(dk,m)K×MFor projection thickness matrix；
μk(Er) it is kth kind material in ErWhen linear attenuation coefficient, ErFor the photon energy of r-th narrow energy spectral coverage；
α is regularization parameter；
Dk=(dk,m)M×1, Q is two dimension or one-dimensional discrete gradient operator；
For Tikhonov (L2) Regularization function, βkFor regularization parameter.
In turn, the corresponding projection thickness matrix D of all materials of measured object of the second decomposition model of calculating satisfaction includes:
Initial value design module 3: the initial value of setting S and D:
S update module 3: fixed D is updated the S for meeting the second decomposition model using non-negative matrix factorization method；
D update module 3: fixed S is updated the D for meeting the second decomposition model using acceleration proximal end gradient algorithm；
Judgment module 3: the difference that D value updates front and back is calculated, judges whether difference meets third preset condition；If so,
Updated D is exported, if it is not, then returning to S update module 3.
In S update module 3: being updated to the S for meeting the second decomposition model using non-negative matrix factorization method and include:
Wherein
I indicates the number of iterations, initial value 1.
In D update module 3: including: using accelerating proximal end gradient algorithm to be updated the D for meeting the second decomposition model
Initial value design module 4:Ci=0；
D update module 4: it updates
Wherein:
I indicates the number of iterations, initial value 1；
tiFor step-size in search,ξ ∈ (0,1),
4th preset condition is that (i- η)~the (i-1) secondary step-size in search is constant；
；
For the m row vector of the transposed matrix of gradient operator Q.
Judgment module 4: judgement
Whether
It sets up, if so, judgment module 3 is executed, if it is not, then enabling Ci=Ci+ 1, return to D update module 4.
Optionally, T is more than or equal to 180.
Optionally, it is θ that T different imaging angles, which include: initial imaging angle,0, the T imaging of △ θ is divided between angle
Angle.
It should be noted that the embodiment of CT imaging device of the invention, identical as the embodiment principle of CT imaging method,
Related place can mutual reference.
The foregoing is merely illustrative of the preferred embodiments of the present invention, not to limit scope of the invention, it is all
Within the spirit and principle of technical solution of the present invention, any modification, equivalent substitution, improvement and etc. done should be included in this hair
Within bright protection scope.
Claims (12)
1. a kind of CT imaging method, which is characterized in that the described method includes:
Step 11: setting the imaging voltage range [V1, V2] of measured object, as V ∈ [V1, V2], the penetrance of the measured object
Meet the relative deviation of linear attenuation coefficient of any two kinds of materials of the first preset condition and the measured object in the V
Meet the corresponding energy range of preset requirement and meets the second preset condition；
Step 12: N number of different scanning voltage is taken in the range of [V1, V2], the N is greater than the material that the measured object is included
Matter sum K, and the difference between any two scanning voltage is greater than the first preset value, to described under each scanning voltage
Measured object carries out full angle scanning imagery, obtains and meets N number of projection image sequence that imaging gray scale requires, the full angle scanning
Including T different imaging angles；
Step 13: being based on N number of projection image sequence, solve the multispectral attenuation meter of change energy for meeting Poisson maximum likelihood principle
Up to the corresponding projection thickness matrix D of all materials of the measured object in formula；
Step 14: being based on the D, calculate the separation projection p of each narrow energy spectral coverage in the X-rayr, r=1,2 ... R, R are described
The narrow energy spectral coverage sum of X-ray；
Step 15: rebuilding the prObtain the reconstruction image a of each narrow energy spectral coverager。
2. the method according to claim 1, wherein the step 13 includes:
F=(f is converted by N number of projection image sequencen,m)N×MThe first decomposition model is inputted, M is the throwing of each scanning voltage
The subelement sum of shadow image sequence, calculates the corresponding projection of all materials of the measured object for meeting first decomposition model
Thickness matrix D；
First decomposition model is
Wherein:
S=(sn,r)N×RFor spectral resolution coefficient matrix；
D=(dk,m)K×MFor projection thickness matrix；
μk(Er) it is kth kind material in ErWhen linear attenuation coefficient, ErFor the photon energy of r-th narrow energy spectral coverage.
3. according to the method described in claim 2, it is characterized in that, described calculate the quilt for meeting first decomposition model
Surveying the corresponding projection thickness matrix D of all materials of object includes:
Step 21: set the initial value of the S and D:
Step 22: the fixed D is updated the S for meeting first decomposition model using non-negative matrix factorization method；
Step 23: the fixed S is updated the D for meeting first decomposition model using acceleration proximal end gradient algorithm；
Step 24: calculating the difference that D value updates front and back, judge whether the difference meets third preset condition；If so, output
Updated D, if it is not, then return step 22.
4. the method according to claim 1, wherein the step 13 includes:
F=(f is converted by N number of projection image sequencen,m)N×MThe second decomposition model is inputted, M is the throwing of each scanning voltage
The subelement sum of shadow image sequence, calculates the corresponding projection of all materials of the measured object for meeting second decomposition model
Thickness matrix D；
Second decomposition model is
Wherein:
S=(sn,r)N×RFor spectral resolution coefficient matrix；
D=(dk,m)K×MFor projection thickness matrix；
μk(Er) it is kth kind material in ErWhen linear attenuation coefficient, ErFor the photon energy of r-th narrow energy spectral coverage；
α is regularization parameter；
Dk=(dk,m)M×1, Q is two dimension or one-dimensional discrete gradient operator；
For Tikhonov (L2) Regularization function, βkFor regularization parameter.
5. according to the method described in claim 4, it is characterized in that, described calculate the quilt for meeting second decomposition model
Surveying the corresponding projection thickness matrix D of all materials of object includes:
Step 31: set the initial value of the S and D:
Step 32: the fixed D is updated the S for meeting second decomposition model using non-negative matrix factorization method；
Step 33: the fixed S is updated the D for meeting second decomposition model using acceleration proximal end gradient algorithm；
Step 34: calculating the difference that D value updates front and back, judge whether the difference meets third preset condition；If so, output
Updated D, if it is not, then return step 32.
6. the method according to claim 3 or 5, which is characterized in that described to use non-negative matrix factorization method to meeting
The S for stating the first decomposition model or second decomposition model, which is updated, includes:
Wherein:
I indicates the number of iterations, initial value 1.
7. according to the method described in claim 3, it is characterized in that, it is described using accelerate proximal end gradient algorithm to meeting described the
The D of one decomposition model, which is updated, includes:
Step 41:Ci=0；
Step 42: updating
Wherein:
I indicates the number of iterations, initial value 1；
tiFor step-size in search,ξ ∈ (0,1),
4th preset condition is that (i- η)~the (i-1) secondary step-size in search is constant；
Step 43: judgement
Whether at
It is vertical, if so, step 24 is executed, if it is not, then enabling Ci=Ci+ 1, return step 42.
8. according to the method described in claim 5, it is characterized in that, it is described using accelerate proximal end gradient algorithm to meeting described the
The D of two decomposition models, which is updated, includes:
Step 51:Ci=0；
Step 52: updating
Wherein:
I indicates the number of iterations, initial value 1；
tiFor step-size in search,ξ ∈ (0,1),
4th preset condition is that (i- η)~the (i-1) secondary step-size in search is constant；
For the m row vector of the transposed matrix of gradient operator Q.
Step 53: judgement
It is whether true, if so, step 34 is executed, if it is not, then enabling Ci=Ci+ 1, return step 52.
9. the method according to claim 1, wherein the T is more than or equal to 180.
10. according to the method described in claim 9, it is characterized in that, the T different imaging angles include: initial imaging
Angle is θ0, the T imaging angle of △ θ is divided between angle.
11. a kind of non-transitory computer-readable storage medium, the non-transitory computer-readable storage medium store instruction is special
Sign is that described instruction executes the processor as described in any in claims 1 to 10
Step in CT imaging method.
12. a kind of CT imaging device, including memory, processor and storage are in the memory and can be in the processor
The computer program of upper operation, which is characterized in that the processor realized when executing the computer program as claim 1 to
Step in 10 in any CT imaging method.
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 true CN109875593A (en) | 2019-06-14 |
| CN109875593B 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) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111134709A (en) * | 2020-01-17 | 2020-05-12 | 清华大学 | Multi-energy CT-based material decomposition method |
| CN112697821A (en) * | 2020-12-02 | 2021-04-23 | 赛诺威盛科技（北京）有限公司 | Multi-energy spectrum CT scanning method and device, electronic equipment and CT equipment |
| CN113392509A (en) * | 2021-05-27 | 2021-09-14 | 南方医科大学 | Accurate estimation method for actual energy spectrum of X-ray |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040264628A1 (en) * | 2003-06-25 | 2004-12-30 | Besson Guy M. | Dynamic multi-spectral imaging with wideband seletable source |
| US20080063135A1 (en) * | 2006-09-08 | 2008-03-13 | General Electric Company | Method and system for generating a multi-spectral image of an object |
| US20100303196A1 (en) * | 2009-05-28 | 2010-12-02 | Kabushiki Kaisha Toshiba | Voltage modulation in dual energy computed tomography |
| US20110282181A1 (en) * | 2009-11-12 | 2011-11-17 | Ge Wang | Extended interior methods and systems for spectral, optical, and photoacoustic imaging |
| 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 |
| US20150103976A1 (en) * | 2013-10-14 | 2015-04-16 | Samsung Electronics Co., Ltd. | X-ray imaging apparatus and control method for the same |
| US20150363947A1 (en) * | 2014-06-16 | 2015-12-17 | The University Of Chicago | Spectral x-ray computed tomography reconstruction using a vectorial total variation |
| CN106483152A (en) * | 2016-11-21 | 2017-03-08 | 中北大学 | A kind of X-ray energy spectrum imaging method |
| US20170270692A1 (en) * | 2015-03-18 | 2017-09-21 | Prismatic Sensors Ab | Image reconstruction based on energy-resolved image data from a photon-counting multi-bin detector |
| 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 |
-
2019
- 2019-02-02 CN CN201910107265.4A patent/CN109875593B/en active Active
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040264628A1 (en) * | 2003-06-25 | 2004-12-30 | Besson Guy M. | Dynamic multi-spectral imaging with wideband seletable source |
| US20080063135A1 (en) * | 2006-09-08 | 2008-03-13 | General Electric Company | Method and system for generating a multi-spectral image of an object |
| US20100303196A1 (en) * | 2009-05-28 | 2010-12-02 | Kabushiki Kaisha Toshiba | Voltage modulation in dual energy computed tomography |
| US20110282181A1 (en) * | 2009-11-12 | 2011-11-17 | Ge Wang | Extended interior methods and systems for spectral, optical, and photoacoustic imaging |
| WO2013063577A1 (en) * | 2011-10-28 | 2013-05-02 | The University Of Chicago | Color x-ray histology for multi-stained biologic sample |
| US20150103976A1 (en) * | 2013-10-14 | 2015-04-16 | Samsung Electronics Co., Ltd. | X-ray imaging apparatus and control method for the same |
| CN103876772A (en) * | 2014-03-20 | 2014-06-25 | 中北大学 | Multi-spectrum imaging method and device |
| US20150363947A1 (en) * | 2014-06-16 | 2015-12-17 | The University Of Chicago | Spectral x-ray computed tomography reconstruction using a vectorial total variation |
| US20170270692A1 (en) * | 2015-03-18 | 2017-09-21 | Prismatic Sensors Ab | Image reconstruction based on energy-resolved image data from a photon-counting multi-bin detector |
| 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 |
|---|
| JIAOTONG WEI等: "Narrow-Energy-Width CT Based on Multivoltage X-Ray Image Decomposition", 《INT J BIOMED IMAGING》 * |
| WOO-JIN LEE: "Material decomposition with the multi-energy attenuation coefficient ratio by using a multiple discriminant analysis", 《JOURNAL OF THE KOREAN PHYSICAL SOCIETY》 * |
| 何鹏: "基于MARS系统的X射线能谱CT研究", 《中国优秀硕士论文全文数据库信息科技辑》 * |
| 戎凯旋: "基于投影替代与矩阵低秩稀疏分解的多光谱图像融合", 《中国优秀硕士论文全文数据库信息科技辑》 * |
| 牛素鋆: "基于能谱滤波分离的彩色CT成像方法研究", 《中国优秀硕士论文全文数据库信息科技辑》 * |
| 魏交统: "基于变电压图像序列盲分离的X射线多谱CT成像", 《中国优秀博士论文全文数据库信息科技辑》 * |
| 魏交统等: "多电压X射线图像分解的高对比度CT成像", 《光电工程》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111134709A (en) * | 2020-01-17 | 2020-05-12 | 清华大学 | Multi-energy CT-based material decomposition method |
| CN111134709B (en) * | 2020-01-17 | 2021-09-14 | 清华大学 | Multi-energy CT-based material decomposition method |
| CN112697821A (en) * | 2020-12-02 | 2021-04-23 | 赛诺威盛科技（北京）有限公司 | Multi-energy spectrum CT scanning method and device, electronic equipment and CT equipment |
| CN112697821B (en) * | 2020-12-02 | 2022-12-02 | 赛诺威盛科技(北京)股份有限公司 | Multi-energy spectrum CT scanning method and device, electronic equipment and CT equipment |
| CN113392509A (en) * | 2021-05-27 | 2021-09-14 | 南方医科大学 | Accurate estimation method for actual energy spectrum of X-ray |
| CN113392509B (en) * | 2021-05-27 | 2022-08-23 | 南方医科大学 | Accurate estimation method for actual energy spectrum of X-ray |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109875593B (en) | 2023-03-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US9128194B2 (en) | Pileup correction method for a photon-counting detector | |
| JP6640498B2 (en) | X-ray computed tomography apparatus, image reconstruction method, and image reconstruction program | |
| Duerinckx et al. | Polychromatic streak artifacts in computed tomography images | |
| CN109875593A (en) | CT imaging method, storage medium and device | |
| KR102424145B1 (en) | Method and apparatus for reconstructing 3D image from spatiotemporal overlapping X-rays | |
| CN108577876A (en) | A kind of static CT of polygon and its working method | |
| Faulkner et al. | Analysis of x-ray computed tomography images using the noise power spectrum and autocorrelation function | |
| CN107822652B (en) | Method for rebuilding spectral results image data | |
| WO2009093305A1 (en) | Positron ct device | |
| CN110415307B (en) | Tensor completion-based multi-energy CT imaging method and device and storage equipment thereof | |
| US12073491B2 (en) | Energy weighting of photon counts for conventional imaging | |
| CN111788499B (en) | Dead time correction method in quantitative Positron Emission Tomography (PET) reconstruction of various objects and radioactivity distributions | |
| Liu et al. | Spectral response model for a multibin photon-counting spectral computed tomography detector and its applications | |
| US20140218362A1 (en) | Monte carlo modeling of field angle-dependent spectra for radiographic imaging systems | |
| CN105581804B (en) | Throughput sub-count detector optimizes signal detection | |
| Zhang et al. | Reconstruction method for DECT with one half-scan plus a second limited-angle scan using prior knowledge of complementary support set (Pri-CSS) | |
| Spronk et al. | Feasibility of a stationary head CT scanner using a CNT x-ray source array | |
| JP6462397B2 (en) | X-ray computed tomography apparatus and image reconstruction method | |
| Hou et al. | Analysis of compressed sensing based CT reconstruction with low radiation | |
| CN111526796B (en) | System and method for image scatter correction | |
| US9508164B2 (en) | Fast iterative image reconstruction method for 3D computed tomography | |
| CN109685865A (en) | It is suitble to the cone-beam cross sectional reconstruction method of linear scanning track | |
| Boquet-Pujadas et al. | PET rebinning with regularized density splines | |
| Tang et al. | Joint Regularized-based Image Reconstruction by Combining Super-Resolution Sinogram for Computed Tomography Imaging | |
| Lee et al. | Design optimization of multi-pinhole micro-SPECT configurations by signal detection tasks and system performance evaluations for mouse cardiac imaging |
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 |