CN103559699B - A kind of Multi-energy-spectruCT CT image reconstruction method estimated based on projection - Google Patents

Abstract

The present invention discloses a kind of Multi-energy-spectruCT CT image reconstruction method estimated based on projection, in order to rebuild the multiple sill density image of testee.The method includes: utilize tradition monoenergetic CT algorithm for reconstructing to be directly reconstructed the energy diagram picture of testee respectively by the polychrome data for projection gathered；Calculate each energy diagram picture line integral along the directions of rays of other power spectrums all, estimate that the polychrome of current power spectrum in those directions projects, try to achieve the polychrome projection that geometric parameter is consistent；Demarcate many sills analytic function, polychrome consistent for the geometric parameter estimated is projected, is decomposed into the line integral of multiple sill；Corresponding many sills density image is reconstructed respectively by many sills line integral.The inventive method is simple and practical, it is adaptable to the multi-power spectrum CT image reconstruction in the case of polychrome projection geometry is inconsistent.Compared with prior art, estimate, merely with surveyed data for projection, the polychrome projection that geometric parameter is consistent, just can reconstruct high-quality image.

Description

A kind of Multi-energy-spectruCT CT image reconstruction method estimated based on projection

Technical field

The present invention relates to X ray CT technical field of imaging, estimate based on projection in particular to one The Multi-energy-spectruCT CT image reconstruction method of meter.

Background technology

X ray computer chromatography imaging technique (abbreviation X ray CT) is that one utilizes X-ray and thing The interaction principle of matter, carries out a kind of technology of imaging to interior of articles information.X-ray ct technology It is not only the medical consultations detection means of a kind of routine, is also an important industrial support technology, extensively It is applied to the Non-Destructive Testing of the critical component of multiple fields.The most conventional X-ray tube can not sent out preferably Monoergic X-ray beam, send is the heterogeneous x ray bundle obeying certain Spectral structure, the number of scanning collection According to for polychrome data for projection, therefore traditional CT is also referred to as monoenergetic spectrum CT.Monoenergetic spectrum CT only uses one Testee is scanned by X-ray energy spectrum, and the CT image of the different material rebuild is likely to be of phase Same or close CT value, it is difficult to distinguish different materials.In addition traditional CT algorithm have ignored X to penetrate The pleochroism of line, causes producing hardening artifact in rebuilding image.During multi-power spectrum CT scan testee, Use multiple X-ray energy spectrum, it is possible to record testee information more more than traditional monoenergetic spectrum CT.Profit Testee equivalent atom ordinal sum electron density image can be reconstructed by these information, or reconstruct The density image of specific base material material, makes a distinction material based on these information.Owing to multi-power spectrum CT has There is more preferable material separating capacity, be therefore with a wide range of applications.

In the case of high-resolution imaging, multi-power spectrum CT image problem be a high dimension non-linear instead Problem.Preferably the data acquisition request of multi-power spectrum CT is " homology, simultaneously, in the same direction ".So-called same Source, i.e. ensures that the quality when X-ray of every kind of power spectrum gathers is completely the same；The while of so-called, i.e. ensure quilt Survey the object displacement when the X-ray of every kind of power spectrum gathers completely the same；What is called in the same direction, i.e. ensures every kind It is completely the same that the X-ray of power spectrum gathers path.In the case of meeting ideal conditions, can be at projector space It is the monochromatic projection of each sill by polychrome data for projection Direct Resolution, then rebuilds by traditional CT and calculate Method such as filter back-projection algorithm etc. can reconstruct the sill image of high-quality testee.But mesh Front conventional multi-power spectrum CT hardware system is all difficult to meet ideal data acquisition condition, is especially difficult to meet " in the same direction " condition, the geometric parameter causing the polychrome of the various power spectrums gathered to project is inconsistent.At this In the case of Zhong, multi-power spectrum imaging equation has the highest pathosis.Existing method for solving substantially can divide It is three classes: iterative method, projection domain decomposition method and image area decomposition method.Iterative method utilize numerical method or Optimization method structure iteration structure, obtains many sills density map by progressively revising image reconstruction result The information such as picture.Although the method reconstructed image quality is higher, but Accurate Calibration in advance is needed to go out x-ray source Multiple power spectrums and the attenuation quotient etc. of sill under different-energy, and computationally intensive, convergence speed Degree is slow, therefore, it is difficult to apply in actual multi-power spectrum CT imaging system.Projection domain decomposes formation method first Decomposite the projection value of sill density image in projection domain, then utilize traditional CT method to rebuild Go out the density image of sill.The method step is simple, image quality is high, but is suitable only for various energy The consistent situation of polychrome projection geometry of spectrum, is not therefore suitable for current multi-power spectrum CT system. Image area decomposition method first calculates high energy and the mental retardation attenuation quotient of object respectively with traditional CT method for reconstructing Image, then carries out linear combination to reconstruction image at image area and reconstructs the density image of sill.Figure Image field decomposition method is a kind of approximation formation method, resolve into as time have ignored the product term of multi-power spectrum projection, The sill density image rebuild has serious artifact.In view of the most conventional dual intensity spectrum CT imaging system And the polychrome projection that the multiple CT scan type collection such as spiral cone-beam tomography is consistent less than geometric parameter, Image area decomposes the multi-power spectrum CT formation method that formation method remains the most frequently used, it is therefore desirable to exploitation one Kind new be applicable to polychrome projection geometry inconsistent in the case of high-quality multi-power spectrum CT image reconstruction Method.

Summary of the invention

The present invention provides a kind of Multi-energy-spectruCT CT image reconstruction method estimated based on projection, in order to improve Multi-power spectrum polychrome projection geometric parameter inconsistent in the case of reconstructed image quality.

For reaching above-mentioned purpose, the invention provides a kind of multi-power spectrum CT image weight estimated based on projection Construction method, comprises the following steps:

S1: utilize monoenergetic CT algorithm for reconstructing the most straight by the polychrome data for projection of each power spectrum gathered Connect the energy diagram picture rebuilding testee；

S2: calculate each energy diagram picture line integral along the directions of rays of other power spectrums all, estimates The polychrome projection of current power spectrum in those directions, obtains the polychrome projection that geometric parameter is consistent；

S3: demarcate many sills analytic function, is multiple base by polychrome Projective decomposition consistent for geometric parameter The line integral of material；

S4: utilize monoenergetic CT image rebuilding method respectively by the line integral of multiple sills reconstruct right The many sills density image answered.

Further, in step s3, the step demarcating many sills analytic function specifically includes:

S31: making one with multiple sill and demarcate die body, in demarcation die body, various sills are each side Various sills line integral combination in all directions in testee is contained in line integral combination upwards, Utilize multi-power spectrum CT to measure with the scan pattern identical with testee and demarcate die body, it is thus achieved that multi-power spectrum is many Color data for projection；

S32: utilizing polychrome projection method of estimation to calculate, to demarcate the consistent multi-power spectrum of geometric parameter of die body many Color projects；

S33: by minimizing between reconstruction image and the die body image that polychrome projection is combined in image area Distance, calculates the coefficient of many sills analytic function, obtains the many sills analytic function demarcated.

Further, in step S31, calibration mold body is combined by multiple sill, wherein sill kind Number is less than or equal to power spectrum number, when sill kind number is less than power spectrum number, by combining different power spectrums Data for projection rebuilds sill density image.

Further, before polychrome projection is estimated, the main component in testee is carried out in advance line Propertyization corrects, to improve projection estimated accuracy.

Further, monoenergetic CT image rebuilding method is analytic method or iterative reconstruction approach.

Further, utilize monoenergetic CT image rebuilding method to carry out in image reconstruction, add physics about Bundle condition, makes the density image currently rebuild meet physical constraint condition, and physical constraint condition includes pixel Value nonnegativity restrictions and/or piecewise smooth constraint.

Further, step S1 utilize monoenergetic CT algorithm for reconstructing by the polychrome data for projection of each power spectrum The formula directly reconstructing testee image respectively is as follows:

$f_{q_m}\left(x\right)=R^{-1}\left(q_m\left(L\right)\right)$

Here q_m(L) representing the polychrome data for projection of m kind power spectrum, m is natural number；R^-1It is Radon Inverse transformation operator, represents monoenergetic CT image rebuilding method；Represent under this kind of power spectrum rebuild The energy diagram picture of testee.

Further, by calculating each energy diagram picture directions of rays along other power spectrum in step S2 Line integral, estimate that the formula of the polychrome projection of current power spectrum in those directions is as follows:

${q^{\′}}_m\left(L\right)={\∫}_L{f}_{q_m}\left(x\right)\mathrm{dl},L\∉{\mathrm{\Ω}}_m$

Here parameter m represents the sequence number of power spectrum, m=1,2 ..., M, M are natural number, i.e. have M energy Spectrum scanning testee；Directions of rays L is directions of rays during other energy-spectrum scanning testee, does not belongs to Ray path set omega in current m kind power spectrum_m；q′_m(L) it is along directions of rays L through current the M energy diagram pictureLine integral, be to m kind power spectrum in L direction polychrome projection estimate.

Further, step S3 utilizes the polychrome projection that estimated geometric parameter is consistent, is decomposed into many The formula of individual sill line integral is:

$F_m\left(L\right)=C_m({q^{\′}}_1\left(L\right),{q^{\′}}_2\left(L\right),...,{q^{\′}}_m\left(L\right))$

$=\underset{k_1,k_2,...,k_M}{\mathrm{\Σ}}{c}_{m,k_1,k_2,...,k_M}{q^{\′}}_1{\left(L\right)}^{k_1}{q^{\′}}_2{\left(L\right)}^{k_2}...{q^{\′}}_M{\left(L\right)}^{k_M}$

Here F_m(L) it is the line integral of wherein m-th sill, m=1,2 ..., M, M are natural number；C_m It is many sills analytic function to be calibrated, represents with many order polynomial function, parameterIt is repeatedly The coefficient of polynomial function；k₁,k₂,…,k_MIt is the index of corresponding polychrome projection respectively；Ray path L∈Ω₁∪Ω₂∪…∪Ω_M, i.e. all set of all x-ray path of M power spectrum；If currently penetrated Thread path L ∈ Ω_m, then q '_m(L)=q_m(L), i.e. use measured obtain polychrome projection q_m(L) conduct q′_m(L)；Otherwise use the polychrome projection q ' estimated_m(L)。

Further, step S4 utilize monoenergetic CT image rebuilding method to be amassed by multiple base material stocklines respectively The formula point reconstructing corresponding many sills density image is:

f_m(x)=R^-1(F_m(L))

Here f_mX () is the many sills density image reconstructed；F_m(L) it is wherein m-th sill Line integral, m=1,2 ..., M, M are natural number；R^-1It is Radon inverse transformation operator, represents monoenergetic CT Image rebuilding method.

The present invention is applicable to the situation that multi-power spectrum polychrome data for projection geometric parameter is the most inconsistent, merely with Surveyed data for projection estimate geometric parameter consistent polychrome projection, just can make image reconstruction mass ratio based on The reconstruction quality of the traditional images domain decomposition method of the inconsistent data for projection of geometric parameter is greatly improved, close to The reconstruction effect of conventional projection domain decomposition method based on the consistent data for projection of geometric parameter.Side the most of the present invention Method need not measured in advance spectral information, and noiseproof feature is strong, and for same class testee, material decomposes Function coefficients only needs to demarcate once.

Accompanying drawing explanation

In order to be illustrated more clearly that the embodiment of the present invention or technical scheme of the prior art, below will be to reality Execute the required accompanying drawing used in example or description of the prior art to be briefly described, it should be apparent that below, Accompanying drawing in description is only some embodiments of the present invention, for those of ordinary skill in the art, On the premise of not paying creative work, it is also possible to obtain other accompanying drawing according to these accompanying drawings.

Fig. 1 is the Multi-energy-spectruCT CT image reconstruction method flow chart of one embodiment of the invention；

Fig. 2 is the faultage image of experiment test die body；

Fig. 3 is the image that die body is demarcated in experiment；

Fig. 4 is the experiment X-ray energy spectrum schematic diagram comprising high and low two kinds of power spectrums；

Fig. 5 a is the projected data image gathered under 80kV tube voltage；

Fig. 5 b is the inconsistent with low power spectrum polychrome projection geometry of collection under 140kV tube voltage High power spectrum polychrome projected data image；

Fig. 5 c is to utilize the high energy consistent with low power spectrum polychrome projection geometry estimated by the inventive method Spectrum polychrome projected data image；

Fig. 6 a by surveyed polychrome data for projection muting in the case of use traditional images domain decomposition method rebuild Bone sill density image；

Fig. 6 b by surveyed polychrome data for projection muting in the case of use traditional images domain decomposition method weight The water material density image built；

Fig. 6 c by surveyed polychrome data for projection muting in the case of use traditional images domain decomposition method rebuild The linear attenuation coefficient image of 70keV photon；

Fig. 6 d by survey polychrome data for projection muting in the case of use the inventive method rebuild bone Sill density image；

Fig. 6 e by survey polychrome data for projection muting in the case of use the inventive method rebuild water base Density of material image；

Fig. 6 f by survey polychrome data for projection muting in the case of use the inventive method rebuild The linear attenuation coefficient image of 70keV photon；

Fig. 6 g by surveyed polychrome data for projection muting in the case of use conventional projection domain decomposition method weight The bone sill density image built；

Fig. 6 h by surveyed polychrome data for projection muting in the case of use conventional projection domain decomposition method weight The water material density image built；

Fig. 6 i by surveyed polychrome data for projection muting in the case of use conventional projection domain decomposition method rebuild The linear attenuation coefficient image of 70keV photon；

Fig. 7 a by surveyed polychrome data for projection noisy in the case of use traditional images domain decomposition method rebuild Bone sill density image；

Fig. 7 b by surveyed polychrome data for projection noisy in the case of use traditional images domain decomposition method weight The water material density image built；

Fig. 7 c by surveyed polychrome data for projection noisy in the case of use traditional images domain decomposition method rebuild The linear attenuation coefficient image of 70keV photon；

Fig. 7 d by survey polychrome data for projection noisy in the case of use the inventive method rebuild bone Sill density image；

Fig. 7 e by survey polychrome data for projection noisy in the case of use the inventive method rebuild water base Density of material image；

Fig. 7 f by survey polychrome data for projection noisy in the case of use the inventive method rebuild The linear attenuation coefficient image of 70keV photon；

Fig. 7 g by surveyed polychrome data for projection noisy in the case of use conventional projection domain decomposition method weight The bone sill density image built；

Fig. 7 h by surveyed polychrome data for projection noisy in the case of use conventional projection domain decomposition method weight The water material density image built；

Fig. 7 i by surveyed polychrome data for projection noisy in the case of use conventional projection domain decomposition method rebuild The linear attenuation coefficient image of 70keV photon.

Detailed description of the invention

Below in conjunction with the accompanying drawing in the embodiment of the present invention, the technical scheme in the embodiment of the present invention is carried out Clearly and completely describe, it is clear that described embodiment is only a part of embodiment of the present invention, and It is not all, of embodiment.Based on the embodiment in the present invention, those of ordinary skill in the art are not paying Go out all other embodiments obtained under creative work premise, broadly fall into the scope of protection of the invention.

Multi-power spectrum CT imaging problem can be expressed as the mathematical model of formula (1):

$q_m\left(L\right)=-\ln \∫S_m\left(E\right)e^{{-\∫}_L\mathrm{\μ}(x,E)\mathrm{dl}}\mathrm{dE},L\∈{\mathrm{\Ω}}_m.---\left(1\right)$

Formula (1) represents with M energy-spectrum scanning testee, obtains M group polychrome data for projection.Wherein parameter M represents the sequence number of power spectrum, m=1,2 ..., M, M are natural number；q_m(L) it is to penetrate with the X of m-th power spectrum During line scanning testee, the X-ray polychromatic projection detected along ray path L；μ (x, E) is About the linear attenuation coefficient of ENERGY E at some x；S_m(E) it is the normalized X-ray energy spectrum of m-th；Ω_m Set for all x-ray path of m-th power spectrum.In the case of preferably " in the same direction ", along appointing Anticipate an x-ray path, the polychrome projection of all power spectrums, the i.e. geometric parameters of polychrome projection can be collected Number is consistent.In a general case, each power spectrum can not be collected along all of ray many Color projects, and even under some multi-power spectrum CT scan patterns, the polychrome of the different power spectrums recorded is projected in Geometric parameter is the most inconsistent.

When with M kind X-ray energy spectrum scanning testee, multi-power spectrum CT obtains M set testee Polychrome data for projection.Utilize these data for projection, can weigh using testee as the combination of M kind material Build out the exploded view picture of this M kind material in testee.To this end, the μ (x, E) in equation (1) is decomposed into M kind as shown in formula (2) and the linear combination of substance characteristics correlative:

Wherein f_mX () represents the density spatial distribution of m kind sill, the most relevant with locus, i.e. base Density of material image, referred to as basic image；Represent the mass attentuation coefficient of m kind sill, only Relevant with energy, can be with measured in advance.Formula (2) is substituted into formula (1) obtain:

Here F_m(L) represent along ray L and pass f_mThe line integral of (x):

F_m(L)=∫_Lf_m(x)dl (4)

Then, multi-power spectrum imaging problem becomes one by M set projection number q_i(L) M sill density is solved Image f_mX the problem of (), i.e. needs the inverse transformation formula of first solution formula (3), then by F_m(L) utilize Tradition monoenergetic CT algorithm for reconstructing reconstructs f_m(x).In the case of high-resolution imaging, multi-power spectrum CT schemes As Problems of Reconstruction is the Nonlinear inverse problem of a high dimension, formula (3) being difficult to analytical Calculation, to obtain it inverse Conversion.Available its inverse transformation of multinomial Function Fitting, then by measurement standard die body to polynary many Item formula function is demarcated, and directly obtains the inverse transformation utilizing polychrome backprojection reconstruction sill image.

When the geometric parameter of the polychrome projection of different power spectrums is consistent, available following multinomial function is intended Conjunction inverse transformation:

$\begin{array}{c}{F}_m\left(L\right)=C_m(q_1\left(L\right),q_2\left(L\right),...,q_M\left(L\right))\\ =\underset{k_1,k_2,...,k_M}{\mathrm{\Σ}}{c}_{m,k_1,k_2,...,k_M}{q}_1{\left(L\right)}^{k_1}{q}_2{\left(L\right)}^{k_2}...q_M{\left(L\right)}^{k_M}\end{array}---\left(5\right)$

Here function C_mIt is that multi-power spectrum polychrome data for projection is decomposed into sill line integral F_m(L) decomposition Function, parameterRepresent the coefficient of polynary many order polynomial function.Formula (5) takes full advantage of many The character of item formula function, experiment shows have the highest fitting precision.Have not yet with in this formula The product term projected with the polychrome of power spectrum, is only applicable to the polychrome projection that geometric parameter is consistent.When different energy Spectrum polychrome be projected in geometric parameter inconsistent time, can not collect not along same x-ray path L Project with the polychrome of power spectrum, therefore can not comprise the product term of the polychrome projection of different power spectrum, at this moment formula (5) become:

$\begin{array}{c}{F}_m\left(L\right)=C_m(q_1\left(L\right),q_2\left(L\right),...,q_M\left(L\right))\\ =\underset{k_1}{\mathrm{\Σ}}{c}_{m,k_1,0,...,0}{q}_1{\left(L\right)}^{k_1}+\underset{k_2}{\mathrm{\Σ}}{c}_{m,0,k_2,...,0}{q}_2{\left(L\right)}^{k_2}+...+\underset{k_M}{\mathrm{\Σ}}{c}_{m,\mathrm{0,0},...,k_M}{q}_m{\left(L\right)}^{k_M}\end{array}---\left(6\right)$

From formula (6), this formula have ignored the product term of all of polychrome projection, it is adaptable to geometry The polychrome projection scanning pattern that parameter is inconsistent, but the fitting precision of formula (3) inverse transformation is substantially reduced, Therefore the method solves error greatly, rebuilds in image and has serious artifact.

To this end, the present invention proposes a kind of Multi-energy-spectruCT CT image reconstruction method estimated based on projection.Should Method is for the high-quality multi-power spectrum CT image reconstruction in the case of polychrome projection geometry is inconsistent. The method it is crucial that first estimate geometry with the inconsistent polychrome data for projection collection of geometric parameter measured The polychrome projection that parameter is consistent, then rebuilds image, specifically comprises the following steps that

1) utilize tradition monoenergetic CT algorithm for reconstructing by the polychrome data for projection of each power spectrum the most directly weight Build the image of corresponding testee:

$f_{q_m}\left(x\right)=R^{-1}\left(q_m\left(L\right)\right)---\left(7\right)$

Here q_m(L) the polychrome data for projection of m kind power spectrum is represented；R^-1It is Radon inverse transformation operator, Represent tradition monoenergetic CT image rebuilding method；Represent the testee under this kind of power spectrum rebuild Image, this be referred to as energy diagram picture by the polychrome image that directly reconstructs of projection.

2) calculate each energy diagram picture line integral along the directions of rays of other power spectrums all, estimate The polychrome projection of current power spectrum in those directions, tries to achieve the polychrome projection that geometric parameter is consistent.

${q^{\′}}_m\left(L\right)={\∫}_L{f}_{q_m}\left(x\right)\mathrm{dl},L\∉{\mathrm{\Ω}}_m---\left(8\right)$

Here directions of rays L is directions of rays during other energy-spectrum scanning testee；q′_m(L) be along Directions of rays L passes current m-th energy diagram pictureLine integral, be in L side to m kind power spectrum To polychrome projection estimate.When the height power spectrum projection gathered is abundant, formula (7) is utilized to rebuild Go out high-quality energy diagram pictureAccording to the character of Radon conversion, calculated by formula (8) Line integral is one to the projection of high-energy power spectrum polychrome and estimates the most accurately, thus can obtain all penetrating The polychrome projection that on line direction, various power spectrums geometric parameter each other is consistent, therefore the inventive method is in theory It is rational.When utilizing sill scaling method to carry out multi-power spectrum image reconstruction, need to gather big angulation Degree polychrome projection, this be estimate high-precision geometric parameter consistent polychrome projection create condition, Therefore the inventive method has the strongest practicality simultaneously.While it is true, due to q_m(L) it is polychrome projection, It not the monochromatic projection obtained under preferable monoergic X-ray scanning, hence with tradition monoenergetic CT weight Build the image that algorithm is rebuildIn there is certain hardening artifact.It is tighter when low energy image has During the hardening artifact weighed, the present invention can only select the polychrome projection gathered during higher-energy energy-spectrum scanning Data estimate the projection of low power spectrum directions of rays, to improve image reconstruction quality.This is because utilize height During energy-spectrum scanning, hardening phenomenon is inconspicuous, the image rebuildCloser to utilizing monochromatic projection The ideal image rebuild, estimated polychrome projection is more accurate.

3) utilize estimated by geometric parameter consistent polychrome projection, be decomposed into M sill line integral.

$\begin{array}{c}{F}_m\left(L\right)=C_m({q^{\′}}_1\left(L\right),{q^{\′}}_2\left(L\right),...,{q^{\′}}_M\left(L\right))\\ =\underset{k_1,k_2,...,k_M}{\mathrm{\Σ}}{c}_{m,k_1,k_2,...,k_M}{q^{\′}}_1{\left(L\right)}^{k_1}{q^{\′}}_2{\left(L\right)}^{k_2}...{q^{\′}}_M{\left(L\right)}^{k_M}\end{array}---\left(9\right)$

In formula (9), F_m(L) it is the line integral of wherein m-th sill；Ray path L∈Ω₁∪Ω₂∪…∪Ω_M, i.e. all set of all x-ray path of M power spectrum.If currently penetrated Thread path L ∈ Ω_m, then q '_m(L)=q_m(L), i.e. use measured obtain polychrome projection q_m(L) conduct q′_m(L)；Otherwise use the polychrome projection q ' that (8) formula estimates_m(L)。

4) utilize tradition monoenergetic CT image rebuilding method respectively by sill line integral F_m(L) reconstruct right The many sills density image f answered_m(x)。

f_m(x)=R^-1(F_m(L)) (10)

In the inventive method, the material analytic function coefficient of formula (9) needs to demarcate, concrete demarcation side Method is:

1) first making one with multiple sill and demarcate die body, in this die body, various sills are at each Line integral combination on direction should be able to be fully contemplated by various sills in testee in all directions Line integral combines.Utilize multi-power spectrum CT to measure with the scan pattern identical with testee and demarcate die body, Obtain multi-power spectrum polychrome data for projection.

2) above-mentioned polychrome projection method of estimation is utilized to calculate the consistent multipotency of geometric parameter demarcating die body Spectrum polychrome projection；

3) figure of various polychromes projection combination in tradition monoenergetic CT algorithm for reconstructing reconstruction formula (9) is utilized Picture.Utilize the linear behavio(u)r that Radon converts, by image area by minimize polychrome projection combination The distance rebuild between image and die body image, calculate the coefficient of analytic function in formula (9)Obtain all analytic function C demarcated_m。

Fig. 1 is the Multi-energy-spectruCT CT image reconstruction method stream estimated based on projection of one embodiment of the invention Cheng Tu；As it can be seen, wherein the flow chart in left side is the mistake demarcating multi-power spectrum CT sill analytic function Journey, the flow chart on right side, for utilizing sill analytic function, carries out the mistake that many sills density image is rebuild Journey.This Multi-energy-spectruCT CT image reconstruction method specifically includes following steps:

S1: utilize monoenergetic CT algorithm for reconstructing the most straight by the polychrome data for projection of each power spectrum gathered Connect the energy diagram picture rebuilding testee；

The polychrome data for projection of each power spectrum passes through multi-power spectrum CT system acquisition, and multi-power spectrum CT system is permissible It is to include the CT system with two kinds of X-ray energy spectrum scannings, the most double spectral CT system, it is also possible to be bag Include the CT system with the scanning of two or more X-ray energy spectrums.Double spectral CT system are multi-power spectrum CT imagings Modal imaging system in field, in following part, the present invention is with dual intensity spectrum CT image reconstruction As a example by enforcement, further illustrate the inventive method implements process.

X-ray energy spectrum projection had both included the broad x-ray power spectrum projection that Traditional x-ray detector detects, Also the narrow X-ray energy spectrum projection that photon counting detector detects is included.The scan pattern bag that can be suitable for Include the consistent scanning of multi-power spectrum geometric parameter and the inconsistent scanning of geometric parameter.

Wherein, monoenergetic CT image rebuilding method can be analytic method or iterative reconstruction approach.

Step S1 utilize monoenergetic CT algorithm for reconstructing by the polychrome data for projection of each power spectrum the most directly weight The formula building testee image is as follows:

$f_{q_m}\left(x\right)=R^{-1}\left(q_m\left(L\right)\right)$

Here q_m(L) representing the polychrome data for projection of m kind power spectrum, m is natural number；R^-1It is Radon Inverse transformation operator, represents monoenergetic CT image rebuilding method；Represent under this kind of power spectrum rebuild The energy diagram picture of testee.

S2: calculate each energy diagram picture line integral along the directions of rays of other power spectrums all, estimates The polychrome projection of current power spectrum in those directions, obtains the polychrome projection that geometric parameter is consistent；

By calculating each energy diagram picture line integral along the directions of rays of other power spectrum in step S2, Estimate that the formula that the polychrome of current power spectrum in those directions projects is as follows:

${q^{\′}}_m\left(L\right)={\∫}_L{f}_{q_m}\left(x\right)\mathrm{dl},L\∉{\mathrm{\Ω}}_m$

Here parameter m represents the sequence number of power spectrum, m=1,2 ..., M, M are natural number, i.e. have M energy Spectrum scanning testee；Directions of rays L is directions of rays during other energy-spectrum scanning testee, does not belongs to Ray path set omega in current m kind power spectrum_m；q′_m(L) it is along directions of rays L through current the M energy diagram pictureLine integral, be to m kind power spectrum in L direction polychrome projection estimate.

Further, step S3 utilizes the polychrome projection that estimated geometric parameter is consistent, is decomposed into many The formula of individual sill line integral is:

$\begin{array}{c}{F}_m\left(L\right)=C_m({q^{\′}}_1\left(L\right),{q^{\′}}_2\left(L\right),...,{q^{\′}}_M\left(L\right))\\ =\underset{k_1,k_2,...,k_M}{\mathrm{\Σ}}{c}_{m,k_1,k_2,...,k_M}{q^{\′}}_1{\left(L\right)}^{k_1}{q^{\′}}_2{\left(L\right)}^{k_2}...{q^{\′}}_M{\left(L\right)}^{k_M}\end{array}$

Here F_m(L) it is the line integral of wherein m-th sill, m=1,2 ..., M, M are natural number；C_mIt is Many sills analytic function to be calibrated, represents with many order polynomial function, parameterIt is the most The coefficient of item formula function；k₁,k₂,…,k_MIt is the index of corresponding polychrome projection respectively；Ray path L∈Ω₁∪Ω₂∪…∪Ω_M, i.e. all set of all x-ray path of M power spectrum.If currently penetrated Thread path L ∈ Ω_m, then q '_m(L)=q_m(L), i.e. use measured obtain polychrome projection q_m(L) conduct q′_m(L)；The polychrome projection q ' otherwise estimated_m(L)。

Wherein, before polychrome projection is estimated, it is also possible to the main component in testee is carried out in advance line Propertyization corrects, to improve projection estimated accuracy.As in medicine CT, carry out " water precorrection ", i.e. Obtain testee and demarcate the equivalent water thickness of die body.

When projection is estimated, it is also possible to adjust projecting direction, generate the projection being more easy to calculate such as parallel beam etc. and sweep Retouch the data for projection under pattern, and make image reconstruction according to this.

S3: demarcate many sills analytic function, is multiple base by polychrome Projective decomposition consistent for geometric parameter The line integral of material；

Wherein, in step s3, the step demarcating many sills analytic function specifically includes:

S31: making one with multiple sill and demarcate die body, in demarcation die body, various sills are each side Various sills line integral combination in all directions in testee is contained in line integral combination upwards, Utilize multi-power spectrum CT to measure with the scan pattern identical with testee and demarcate die body, it is thus achieved that multi-power spectrum is many Color data for projection；

Demarcation die body in step S31 is combined by multiple sill, and wherein sill kind number is less than In power spectrum number, when sill kind number is less than power spectrum number, by combining different power spectrum data for projection Rebuild sill density image.Demarcate the thickness in different CT scan directions of the multiple sill in die body Combination can contain these sills all possible thickness combination in testee.

S32: utilizing polychrome projection method of estimation to calculate, to demarcate the consistent multi-power spectrum of geometric parameter of die body many Color projects；

S33: by minimizing between reconstruction image and the die body image that polychrome projection is combined in image area Distance, calculates the coefficient of many sills analytic function, obtains the many sills analytic function demarcated.

Wherein for same class testee, many sills analytic function is demarcated only to be needed to demarcate once.

S4: utilize monoenergetic CT image rebuilding method respectively by the line integral of multiple sills reconstruct right The many sills density image answered.

Wherein, monoenergetic CT image rebuilding method is being utilized to carry out in image reconstruction, for improving image reconstruction Quality can also add physical constraint condition, makes the density image currently rebuild meet physical constraint condition, Physical constraint condition includes pixel value nonnegativity restrictions and/or piecewise smooth constraint.

Step S4 utilize monoenergetic CT image rebuilding method reconstructed institute by multiple sill line integrals respectively The formula of corresponding many sills density image is:

f_m(x)=R^-1(F_m(L))

Here f_mX () is the many sills density image reconstructed；F_m(L) it is wherein m-th sill Line integral, m=1,2 ..., M, M are natural number；R^-1It is Radon inverse transformation operator, represents monoenergetic CT Image rebuilding method.

It is the most inconsistent that the above embodiment of the present invention is applicable to multi-power spectrum polychrome data for projection geometric parameter Situation, estimates, merely with surveyed data for projection, the polychrome projection that geometric parameter is consistent, just can make image weight Build the reconstruction quality of mass ratio traditional images domain decomposition method based on the inconsistent data for projection of geometric parameter significantly Improve, close to the reconstruction effect of conventional projection domain decomposition method based on the consistent data for projection of geometric parameter. Additionally the inventive method need not measured in advance spectral information, and noiseproof feature is strong, for same class measured object Body, material analytic function coefficient only needs to demarcate once.

Fig. 2 is the faultage image of test die body used in one embodiment of the invention, and this test die body is The Forbild thoracic cavity model of one standard.Selected water and osseous tissue are double-basis material, and the density of water is 1.0g/cm³, the density of osseous tissue is 1.92g/cm³.Fig. 3 is demarcation die body used in the present embodiment, For demarcating the coefficient of analytic function.This die body comprises a bigger disk and a less disk, Wherein a diameter of 120mm of larger disk, less a diameter of 70mm.Big disk is filled with water, And roundel is osseous tissue.Double spectral CT system send high and low X-ray energy spectrum as shown in Figure 4.These are two years old Individual power spectrum is that the tube voltage of X-ray tube is respectively 80kV and 140kV(1mm copper filter plate) under send. When Numerical Implementation, the sampling interval of power spectrum is set to 1keV.In dual intensity spectrum CT scan, fan-beam is used to sweep Retouch.The sweep parameter of system is: the distance of radiographic source to detector is 1200mm, and radiographic source is to turntable The distance at center is 1000mm, and detector is made up of the probe unit of 1024 0.3mm, the visual field straight Footpath is 256mm.Radiographic source rotates a circle, and the staggered high energy uniformly gathering 1440 angles and mental retardation are thrown Each one group of shadow data, the most adjacent high energy projection and mental retardation projection have the angle of 0.125 degree, so obtain Inconsistent high and low of geometric parameter can polychrome data for projection.Rebuilding image size is 512 × 512.? Collect respectively under 80kV and 140kV tube voltage test die body polychrome data for projection as it is shown in figure 5, The low power spectrum polychrome data for projection gathered under wherein Fig. 5 a is 80kV tube voltage；Fig. 5 b is 140kV pipe electricity The high power spectrum polychrome data for projection inconsistent with low power spectrum polychrome projection geometry that pressure gathers；Fig. 5 c For utilizing the high power spectrum polychrome consistent with low power spectrum polychrome projection geometry estimated by the inventive method to throw Shadow data.With normalization mean square distance (NMSB) and normalization average absolute distance (NMAD) and the party Comparing as quantization to real high power spectrum polychrome data for projection, the value of NMSB and NMAD is respectively 0.012356 and 0.009745, show that the high power spectrum polychrome data for projection precision that the inventive method is estimated is the highest.

Fig. 6 by surveyed polychrome data for projection muting in the case of, traditional images domain decomposition method, this Bright method and conventional projection domain decomposition method these three method rebuild the comparison of effect, wherein Fig. 6 a, 6b, 6c The image rebuild for traditional images domain decomposition method, Fig. 6 d, 6e, 6f are the image that this method is rebuild, Fig. 6 g, 6h, 6i are the image that conventional projection domain decomposition method is rebuild, and Fig. 6 a, 6d, 6g are water material density image, Fig. 6 b, 6e, 6h are bone sill density images, Fig. 6 c, 6f, 6i be energy be the photon of 70keV Linear attenuation coefficient image.Here linear attenuation coefficient image is by the water base image that will rebuild and water pair The product of the x-ray photon mass attentuation coefficient of prescribed energy and the bone basic image of reconstruction and bone are to appointment The product summation of the x-ray photon mass attentuation coefficient of energy obtains.In order to highlight rebuild image Feature, has used narrow display gray scale window in Fig. 6, wherein the display gray scale window of Fig. 6 a, 6d, 6g is [0.951.15], The display gray scale window of Fig. 6 b, 6e, 6h is [0.71.2], and the display gray scale window of Fig. 6 c, 6f, 6i is [-150HU 150HU].In Fig. 6, traditional images domain decomposition method and the data for projection used by the inventive method are all geometric parameters The multi-power spectrum polychrome projection that number is inconsistent, and the data for projection used by the decomposition method of conventional projection territory is geometry The multi-power spectrum polychrome projection that parameter is consistent, as the reference of a Perfect Reconstruction result.In comparison diagram 6 three The image that the method for kind is rebuild respectively, it can be seen that the method for the embodiment of the present invention can decomposite double well Sill density image, only has some slight artifacts in sill density image.It is far superior to tradition figure The reconstructed results of image field decomposition method, particularly water material density image.Throw with based on geometric parameter is consistent The conventional projection domain decomposition method of shadow is compared, although the reconstruction effect of the inventive method is the most closely, special The quality not being synthesized linear attenuation coefficient image is the most suitable.

Fig. 7 by surveyed polychrome data for projection noisy in the case of, traditional images domain decomposition method, this Bright method and conventional projection domain decomposition method these three method rebuild the comparison of effect.Equally, Fig. 7 a, 7b, 7c is the image that traditional images domain decomposition method is rebuild, and Fig. 7 d, 7e, 7f are the image that this method is rebuild, figure 7g, 7h, 7i are the image that conventional projection domain decomposition method is rebuild, and Fig. 7 a, 7d, 7g are bone sill density Image, Fig. 7 b, 7e, 7h are water material density images, Fig. 7 c, 7f, 7i be energy be 70keV The linear attenuation coefficient image of photon.The noise added in polychrome data for projection is corresponding every ray The poisson noise of 1,000,000 photons.In order to highlight the feature of rebuild image, Fig. 7 still uses The narrow display gray scale window identical with in Fig. 6.As can see from Figure 7, at surveyed data for projection by noise Under pollution condition, the inventive method remains able to reconstruct high-quality sill image and linear attenuation system Number image, reconstructed results is still better than the reconstructed results of traditional images domain decomposition method, close to conventional projection The reconstructed results of domain decomposition method.But three kinds of methods are owing to being disturbed by noise, in tested die body The feature of the most tiny low contrast all can not be rebuild out.

One of ordinary skill in the art will appreciate that: accompanying drawing is the schematic diagram of an embodiment, in accompanying drawing Module or flow process not necessarily implement necessary to the present invention.

One of ordinary skill in the art will appreciate that: the module in device in embodiment can be according to enforcement Example describes in the device being distributed in embodiment, it is also possible to carries out respective change and is disposed other than the present embodiment In one or more devices.The module of above-described embodiment can merge into a module, it is also possible to further Split into multiple submodule.

Last it is noted that above example is only in order to illustrate technical scheme, rather than to it Limit；Although the present invention being described in detail with reference to previous embodiment, the ordinary skill of this area Personnel it is understood that the technical scheme described in previous embodiment still can be modified by it, or Wherein portion of techniques feature is carried out equivalent；And these amendments or replacement, do not make relevant art The essence of scheme departs from the spirit and scope of embodiment of the present invention technical scheme.

Claims (9)

1. one kind based on projection estimate Multi-energy-spectruCT CT image reconstruction method, it is characterised in that include with Lower step:

S1: utilize monoenergetic CT algorithm for reconstructing the most straight by the polychrome data for projection of each power spectrum gathered Connect the energy diagram picture rebuilding testee；

S2: calculate each energy diagram picture line integral along the directions of rays of other power spectrums all, estimates The polychrome projection of current power spectrum in those directions, obtains the polychrome projection that geometric parameter is consistent；

S3: demarcate many sills analytic function, is many by polychrome Projective decomposition consistent for described geometric parameter The line integral of individual sill；

S4: utilize monoenergetic CT image rebuilding method respectively by the line integral of multiple sills reconstruct right The many sills density image answered；

In step s3, the step demarcating many sills analytic function specifically includes:

S31: making one with multiple sill and demarcate die body, in demarcation die body, various sills are each side Various sills line integral combination in all directions in testee is contained in line integral combination upwards, Utilize multi-power spectrum CT to measure with the scan pattern identical with testee and demarcate die body, it is thus achieved that multi-power spectrum is many Color data for projection；

S32: utilizing polychrome projection method of estimation to calculate, to demarcate the consistent multi-power spectrum of geometric parameter of die body many Color projects；

S33: by minimizing between reconstruction image and the die body image that polychrome projection is combined in image area Distance, calculates the coefficient of many sills analytic function, obtains the many sills analytic function demarcated.

Multi-energy-spectruCT CT image reconstruction method the most according to claim 1, it is characterised in that step Demarcating die body described in S31 to be combined by multiple sill, wherein sill kind number is less than or equal to power spectrum Number, when sill kind number is less than power spectrum number, the power spectrum data for projection different by combination rebuilds base Density of material image.

Multi-energy-spectruCT CT image reconstruction method the most according to claim 1, it is characterised in that many Before color projection is estimated, the main component in testee is carried out in advance Linearized correction, to improve throwing Shadow estimated accuracy.

Multi-energy-spectruCT CT image reconstruction method the most according to claim 1, it is characterised in that described Monoenergetic CT image rebuilding method is analytic method or iterative reconstruction approach.

Multi-energy-spectruCT CT image reconstruction method the most according to claim 1, it is characterised in that in profit Carry out in image reconstruction with monoenergetic CT image rebuilding method, add physical constraint condition, make currently to rebuild Density image meet described physical constraint condition, described physical constraint condition includes pixel value nonnegativity restrictions And/or piecewise smooth constraint.

Multi-energy-spectruCT CT image reconstruction method the most according to claim 1, it is characterised in that step S1 utilize monoenergetic CT algorithm for reconstructing directly reconstructed measured object respectively by the polychrome data for projection of each power spectrum The formula of body image is as follows:

$f_{q_m}\left(x\right)=R^{-1}\left(q_m\left(L\right)\right)$

Here q_m(L) representing the polychrome data for projection of m kind power spectrum, m is natural number；R^-1It is Radon Inverse transformation operator, represents monoenergetic CT image rebuilding method；Represent under this kind of power spectrum rebuild The energy diagram picture of testee.

Multi-energy-spectruCT CT image reconstruction method the most according to claim 1, it is characterised in that step By calculating each energy diagram picture line integral along the directions of rays of other power spectrums all in S2, estimate The formula of the polychrome projection of current power spectrum in those directions is as follows:

${q^{\′}}_m\left(L\right)={\∫}_L{f}_{q_m}\left(x\right)dl,L\∉{\mathrm{\Ω}}_m$

Here parameter m represents the sequence number of power spectrum, m=1,2 ..., M, M are natural number, i.e. have M energy Spectrum scanning testee；Directions of rays L is directions of rays during other energy-spectrum scanning testee, does not belongs to Ray path set omega in current m kind power spectrum_m；q'_m(L) it is along directions of rays L through current the M energy diagram pictureLine integral, be to m kind power spectrum in L direction polychrome projection estimate.

Multi-energy-spectruCT CT image reconstruction method the most according to claim 1, it is characterised in that step S3 utilizes the polychrome projection that estimated geometric parameter is consistent, is decomposed into the public affairs of multiple sill line integral Formula is:

$\begin{array}{c}{F}_m\left(L\right)=C_m\left({q^{\′}}_1\right(L),{q^{\′}}_2(L),\mathrm{...},{q^{\′}}_M(L\left)\right)\\ =\underset{k_1,k_2,\mathrm{...},k_M}{\Σ}{c}_{m,k_1,k_2,\mathrm{...},k_M}{q^{\′}}_1{\left(L\right)}^{k_1}{q^{\′}}_2{\left(L\right)}^{k_2}\mathrm{...}{q^{\′}}_M{\left(L\right)}^{k_M}\end{array}$

Here F_m(L) it is the line integral of wherein m-th sill, m=1,2 ..., M, M are natural number；C_m It is many sills analytic function to be calibrated, represents with many order polynomial function, parameterIt is repeatedly The coefficient of polynomial function；k₁,k₂,...,k_MIt is the index of corresponding polychrome projection respectively；Ray path L∈Ω₁∪Ω₂∪...∪Ω_M, i.e. all set of all x-ray path of M power spectrum；If currently penetrated Thread path L ∈ Ω_m, then q'_m(L)=q_m(L), i.e. use measured obtain polychrome projection q_m(L) conduct q'_m(L)；Otherwise use the polychrome projection q' estimated_m(L)。

Multi-energy-spectruCT CT image reconstruction method the most according to claim 1, it is characterised in that step S4 utilize monoenergetic CT image rebuilding method reconstructed corresponding many by multiple sill line integrals respectively The formula of sill density image is:

f_m(x)=R^-1(F_m(L))

Here f_mX () is the many sills density image reconstructed；F_m(L) it is wherein m-th sill Line integral, m=1,2 ..., M, M are natural number；R^-1It is Radon inverse transformation operator, represents monoenergetic CT Image rebuilding method.