CN109875593A - CT imaging method, storage medium and device - Google Patents

CT imaging method, storage medium and device Download PDF

Info

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
measured object
decomposition model
matrix
projection
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
Application number
CN201910107265.4A
Other languages
Chinese (zh)
Other versions
CN109875593B (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

Landscapes

  • Analysing Materials By The Use Of Radiation (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

本发明提供一种CT成像方法、存储介质和装置，该方法包括：步骤11：设定被测物的成像电压范围[V1,V2]；步骤12：在[V1,V2]的范围内取N个不同的扫描电压，N大于被测物所包含的材质总数K，且任意两个扫描电压之间的差值大于第一预设值，在每个扫描电压下对被测物进行全角度扫描成像，获得符合成像灰度要求的N个投影图像序列；步骤13：基于N个投影图像序列，求解满足泊松极大似然原理的的变能量多谱衰减表达式中的被测物的所有材质对应的投影厚度矩阵D；步骤14：基于D，计算X射线中每个窄能谱段的分离投影pr，r＝1,2…R；步骤15：重建pr得到每个窄能谱段的重建图像ar。本发明的方法可以实现低成本、能谱分辨率高的多能谱CT成像。

The present invention provides a CT imaging method, storage medium and device. The method includes: step 11: setting the imaging voltage range [V1, V2] of the measured object; step 12: taking N within the range [V1, V2] different scanning voltages, N is greater than the total number K of materials contained in the measured object, and the difference between any two scanning voltages is greater than the first preset value, and the measured object is scanned at a full angle under each scanning voltage Image, obtain N projection image sequences that meet the imaging grayscale requirements; Step 13: Based on the N projection image sequences, solve all the measured objects in the variable energy multispectral attenuation expression that satisfies the Poisson maximum likelihood principle The projection thickness matrix D corresponding to the material; Step 14: Based on D, calculate the separation projection pr of each narrow energy spectrum in the X-ray, r =1,2... R ; Step 15: Reconstruct pr to obtain each narrow energy spectrum The reconstructed image a r of the segment. The method of the invention can realize multi-energy spectral CT imaging with low cost and high spectral resolution.

Description

CT imaging method, storage medium and device
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.一种CT成像方法，其特征在于，所述方法包括：1. A CT imaging method, wherein the method comprises: 步骤11：设定被测物的成像电压范围[V1,V2]，当V∈[V1,V2]时，所述被测物的穿透率符合第一预设条件、且所述被测物的任意两种材质在所述V时的线性衰减系数的相对偏差满足预设要求对应的能量范围符合第二预设条件；Step 11: Set the imaging voltage range [V1, V2] of the measured object, when V∈[V1, V2], the penetration rate of the measured object meets the first preset condition, and the measured object The relative deviation of the linear attenuation coefficients of any two materials at V meets the preset requirement, and the corresponding energy range meets the second preset condition; 步骤12：在[V1,V2]的范围内取N个不同的扫描电压，所述N大于所述被测物所包含的材质总数K，且任意两个扫描电压之间的差值大于第一预设值，在每个所述扫描电压下对所述被测物进行全角度扫描成像，获得符合成像灰度要求的N个投影图像序列，所述全角度扫描包括T个不同的成像角度；Step 12: Take N different scanning voltages within the range of [V1, V2], where N is greater than the total number K of materials contained in the measured object, and the difference between any two scanning voltages is greater than the first The preset value is to perform full-angle scanning imaging on the measured object under each of the scanning voltages to obtain N projection image sequences that meet the imaging grayscale requirements, and the full-angle scanning includes T different imaging angles; 步骤13：基于所述N个投影图像序列，求解满足泊松极大似然原理的变能量多谱衰减表达式中的所述被测物的所有材质对应的投影厚度矩阵D；Step 13: Based on the N projection image sequences, solve the projection thickness matrix D corresponding to all materials of the measured object in the variable energy multispectral attenuation expression satisfying the Poisson maximum likelihood principle; 步骤14：基于所述D，计算所述X射线中每个窄能谱段的分离投影pr，r＝1,2…R，R为所述X射线的窄能谱段总数；Step 14: Based on the D, calculate the separation projection pr of each narrow energy spectral segment in the X-ray, r =1,2...R, where R is the total number of narrow energy spectral segments of the X-ray; 步骤15：重建所述pr得到所述每个窄能谱段的重建图像ar。Step 15: Reconstruct the pr to obtain the reconstructed image ar of each narrow energy spectrum segment. 2.根据权利要求1所述的方法，其特征在于，所述步骤13包括：2. The method according to claim 1, wherein the step 13 comprises: 将所述N个投影图像序列转化为F＝(fn,m)N×M输入第一分解模型，M为每个扫描电压的投影图像序列的子单元总数，计算满足所述第一分解模型的所述被测物所有材质对应的投影厚度矩阵D；Convert the N projection image sequences into F=(f n,m ) N×M and input it into the first decomposition model, where M is the total number of subunits of the projection image sequence of each scanning voltage, and the calculation satisfies the first decomposition model The projected thickness matrix D corresponding to all materials of the measured object; 所述第一分解模型为The first decomposition model is 其中：in: S＝(sn,r)N×R为光谱分解系数矩阵；S=(s n,r ) N×R is the spectral decomposition coefficient matrix; D＝(dk,m)K×M为投影厚度矩阵；D=(d k,m ) K×M is the projected thickness matrix; μk(Er)为第k种材质在Er时的线性衰减系数，Er为第r个窄能谱段的光子能量。μ k (E r ) is the linear attenuation coefficient of the kth material at Er , and Er is the photon energy of the rth narrow energy spectrum. 3.根据权利要求2所述的方法，其特征在于，所述计算满足所述第一分解模型的所述被测物所有材质对应的投影厚度矩阵D包括：3. The method according to claim 2, wherein the calculating the projection thickness matrix D corresponding to all materials of the measured object satisfying the first decomposition model comprises: 步骤21：设定所述S和D的初始值： Step 21: Set the initial values of S and D: 步骤22：固定所述D，使用非负矩阵分解方法对满足所述第一分解模型的S进行更新；Step 22: fixing the D, and using a non-negative matrix decomposition method to update the S that satisfies the first decomposition model; 步骤23：固定所述S，使用加速近端梯度算法对满足所述第一分解模型的D进行更新；Step 23: fix the S, and use the accelerated proximal gradient algorithm to update the D that satisfies the first decomposition model; 步骤24：计算D值更新前后的差值，判断所述差值是否满足第三预设条件；如果是，输出更新后的D，如果否，则返回步骤22。Step 24 : Calculate the difference before and after the D value is updated, and determine whether the difference satisfies the third preset condition; if so, output the updated D; if not, return to step 22 . 4.根据权利要求1所述的方法，其特征在于，所述步骤13包括：4. The method according to claim 1, wherein the step 13 comprises: 将所述N个投影图像序列转化为F＝(fn,m)N×M输入第二分解模型，M为每个扫描电压的投影图像序列的子单元总数，计算满足所述第二分解模型的所述被测物所有材质对应的投影厚度矩阵D；Convert the N projection image sequences into F=(f n,m ) N×M and input the second decomposition model, where M is the total number of subunits of the projection image sequence of each scanning voltage, and the calculation satisfies the second decomposition model The projected thickness matrix D corresponding to all materials of the measured object; 所述第二分解模型为 The second decomposition model is 其中：in: S＝(sn,r)N×R为光谱分解系数矩阵；S=(s n,r ) N×R is the spectral decomposition coefficient matrix; D＝(dk,m)K×M为投影厚度矩阵；D=(d k,m ) K×M is the projected thickness matrix; μk(Er)为第k种材质在Er时的线性衰减系数，Er为第r个窄能谱段的光子能量；μ k (E r ) is the linear attenuation coefficient of the kth material at Er , and Er is the photon energy of the rth narrow energy spectrum; α为正则化参数；α is the regularization parameter; Dk＝(dk,m)M×1，Q为二维或一维离散梯度算子； D k =(d k,m ) M×1 , Q is a two-dimensional or one-dimensional discrete gradient operator; 为Tikhonov(L2)正则化函数，βk为正则化参数。 is the Tikhonov (L2) regularization function, and β k is the regularization parameter. 5.根据权利要求4所述的方法，其特征在于，所述计算满足所述第二分解模型的所述被测物所有材质对应的投影厚度矩阵D包括：5. The method according to claim 4, wherein the calculating the projection thickness matrix D corresponding to all materials of the measured object satisfying the second decomposition model comprises: 步骤31：设定所述S和D的初始值： Step 31: Set the initial values of S and D: 步骤32：固定所述D，使用非负矩阵分解方法对满足所述第二分解模型的S进行更新；Step 32: fixing the D, and using the non-negative matrix decomposition method to update the S that satisfies the second decomposition model; 步骤33：固定所述S，使用加速近端梯度算法对满足所述第二分解模型的D进行更新；Step 33: fix the S, and use the accelerated proximal gradient algorithm to update the D that satisfies the second decomposition model; 步骤34：计算D值更新前后的差值，判断所述差值是否满足第三预设条件；如果是，输出更新后的D，如果否，则返回步骤32。Step 34: Calculate the difference before and after the D value is updated, and determine whether the difference satisfies the third preset condition; if so, output the updated D; if not, return to step 32 . 6.根据权利要求3或5所述的方法，其特征在于，所述使用非负矩阵分解方法对满足所述第一分解模型或所述第二分解模型的S进行更新包括：6. The method according to claim 3 or 5, wherein the updating S satisfying the first decomposition model or the second decomposition model using a non-negative matrix decomposition method comprises: 其中：in: i表示迭代次数，初始值为1。 i represents the number of iterations, and the initial value is 1. 7.根据权利要求3所述的方法，其特征在于，所述使用加速近端梯度算法对满足所述第一分解模型的D进行更新包括：7. The method according to claim 3, wherein the updating D satisfying the first decomposition model using an accelerated proximal gradient algorithm comprises: 步骤41：Ci＝0；Step 41: C i =0; 步骤42：更新 Step 42: Update 其中：in: i表示迭代次数，初始值为1；i represents the number of iterations, and the initial value is 1; ti为搜索步长，ξ∈(0,1)， t i is the search step size, ξ∈(0,1), 所述第四预设条件为第(i-η)～第(i-1)次的搜索步长不变；The fourth preset condition is that the search step size of the (i-n)th to (i-1)th times is unchanged; 步骤43：判断Step 43: Judgment 是否成立，如果是，执行步骤24，如果否，则令Ci＝Ci+1，返回步骤42。 If yes, go to step 24; if no, set C i =C i +1, and return to step 42 . 8.根据权利要求5所述的方法，其特征在于，所述使用加速近端梯度算法对满足所述第二分解模型的D进行更新包括：8. The method according to claim 5, wherein the updating D satisfying the second decomposition model using an accelerated proximal gradient algorithm comprises: 步骤51：Ci＝0；Step 51: C i =0; 步骤52：更新 Step 52: Update 其中：in: i表示迭代次数，初始值为1； i represents the number of iterations, and the initial value is 1; ti为搜索步长，ξ∈(0,1)， t i is the search step size, ξ∈(0,1), 所述第四预设条件为第(i-η)～第(i-1)次的搜索步长不变；The fourth preset condition is that the search step size of the (i-n)th to (i-1)th times is unchanged; 为梯度算子Q的转置矩阵的第m行向量。 is the mth row vector of the transpose matrix of the gradient operator Q. 步骤53：判断Step 53: Judgment 是否成立，如果是，执行步骤34，如果否，则令Ci＝Ci+1，返回步骤52。 If yes, go to step 34; if no, set C i =C i +1, and return to step 52 . 9.根据权利要求1所述的方法，其特征在于，所述T大于等于180。9 . The method of claim 1 , wherein the T is greater than or equal to 180. 10 . 10.根据权利要求9所述的方法，其特征在于，所述T个不同的成像角度包括：初始成像角度为θ0、角度间隔为△θ的T个成像角度。10 . The method according to claim 9 , wherein the T different imaging angles comprise: T imaging angles with an initial imaging angle of θ 0 and an angular interval of Δθ. 11 . 11.一种非瞬时计算机可读存储介质，所述非瞬时计算机可读存储介质存储指令，其特征在于，所述指令在由处理器执行时使得所述处理器执行如权利要求1至10中任一所述的CT成像方法中的步骤。11. A non-transitory computer-readable storage medium storing instructions, wherein the instructions, when executed by a processor, cause the processor to perform as in claims 1 to 10 Any of the steps in the CT imaging method. 12.一种CT成像装置，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至10中任一所述的CT成像方法中的步骤。12. A CT imaging apparatus, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, when the processor executes the computer program, the computer program as claimed in the claim is implemented Steps in the CT imaging method of any one of claims 1 to 10.
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 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)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111134709A (en) * 2020-01-17 2020-05-12 清华大学 A kind of 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)

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

Patent Citations (11)

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

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

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111134709A (en) * 2020-01-17 2020-05-12 清华大学 A kind of 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
CN103501702B (en) Medical image-processing apparatus, medical image processing method
CN108010099B (en) X-ray multi-energy spectrum CT finite angle scanning and image iterative reconstruction method
Duerinckx et al. Polychromatic streak artifacts in computed tomography images
EP3050028B1 (en) Joint reconstruction of electron density images.
US9128194B2 (en) Pileup correction method for a photon-counting detector
JP6640498B2 (en) X-ray computed tomography apparatus, image reconstruction method, and image reconstruction program
KR102424145B1 (en) Method and apparatus for reconstructing 3D image from spatiotemporal overlapping X-rays
EP2843623A2 (en) X-ray dual-energy CT reconstruction method
US10314556B2 (en) Optimal energy weighting of dark field signal in differential phase contrast X-ray imaging
KR20120055468A (en) Image processing apparatus, image processing method, and non-transitory storage medium
CN109875593A (en) CT imaging method, storage medium and device
JP2016049455A (en) X-ray computer tomographic apparatus and image reconstitution apparatus
WO2009093305A1 (en) Positron ct device
US20220028127A1 (en) Energy weighting of photon counts for conventional imaging
Bubba et al. Tomographic X-ray data of carved cheese
CN111278362A (en) System and method for low dose multi-spectral X-ray tomography
WO2017006764A1 (en) Image computing device, image computing method, and tomograph
CN109448071A (en) A kind of power spectrum image rebuilding method and system
Rigaud et al. Reconstruction algorithm for 3D Compton scattering imaging with incomplete data
Hou et al. Analysis of compressed sensing based CT reconstruction with low radiation
Spronk et al. Feasibility of a stationary head CT scanner using a CNT x-ray source array
US9508164B2 (en) Fast iterative image reconstruction method for 3D computed tomography
KR20140000823A (en) Method for reconstructing of image for polychromatic x-ray tomography
Ramyar et al. Establishing a method to measure bone structure using spectral CT
Ji et al. An experimental method to directly measure DQE $(k) $ at k= 0 for 2D x-ray imaging systems

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