CN1818973A - CT projection data beam hardening effect correcting method based on IIL consistency condition - Google Patents

Abstract

An emendation method of CT projective data beam vulcanized effect base on the relativity of the data, the method is construct a current emendation model of CT projective data-beam vulcanized effect base on the physics imaging module of the medical diagnosed X-ray; construct the relative matrix equation of the CT projective data according to the consistency condition of H-L; get the parameter derivative equation of the emendation module by combining the emendation module and the relative matrix; estimating the attenuation proportion of the high-density substance in the projective data. The invention is fit for achievement of any X-ray and CT equipment beam vulcanized emendation function. Compare with the existing module emendation method, it can automatic fit for the different objects and consider the difference of different density organization.

Description

CT data for projection beam hardening effect correction method based on the HL condition for consistence

Technical field

The invention belongs to X line CT reestablishment imaging technical field, relate to beam hardening antidote based on the X line sequence projecting image data of equivalent nonlinear model and H-L condition for consistence.

Background technology

The principle that CT rebuilds is based on RADON conversion and inverse transformation theory [1] thereof, and its basis is that data for projection comes from the linear projection model, such as, for two-dimentional fault imaging, following projection formula is arranged:

$p(\mathrm{\θ},l)={\∫}_{L(\mathrm{\θ},l)}\mathrm{\μ}(x,y)\mathrm{dxdy}$

Wherein, and L (θ, l) for ray penetrates the path, θ is a projection angle, l is the sampled point of data for projection;

Though the CT technology develops into the 5th current generation from the first generation, track while scan is from the earliest straight line, until present circle and helical trajectory and even any track, projection pattern has fan-shaped and pencil-beam, and the basis of its projection and reconstruction theory is always based on the linear projection model.The research of the reconstruction theory that carried out is always carried out round different track while scans and projection pattern.

In theory, if the X line source of CT is a monoenergetic, the requirement of so above-mentioned linear projection model is then satisfied, and this is because of the decay exponentially attenuated form of material to particular energy X linear light, can satisfy so-called linear projection model after taking the logarithm, as shown in the formula:

$I(\mathrm{\θ},l)=I_0\exp \{-{\∫}_{L(\mathrm{\θ},l)}u(x,y)\mathrm{dxdy}\}$

In the following formula, I ₀For bulb sends the sub-intensity of X linear light, (x y) multiply by material density for the mass attenuation coefficient of imaging object to u.

In fact, the X line spectrum sent of the bulb of modern CT is the wide range function with certain spectral distribution.Under the commonly used 120kVp voltage of Cranial Computed Tomography, the X line spectrum that bulb sends is an example, the sub-energy of its X linear light from 30keV to 120keV according to the certain rules continuous distribution, and also there are 4 discontinuous characteristic radiation spectrum in its distribution function at the 58-68keV place, shown in figure (1), the X line spectrum functional digraph (top) that the CT bulb sends among the figure, lower graph is that formed projection imaging data are as follows through the spectral function (having amplified 5 times) after the 35mm aluminium filtration (according to the common way of CT equipment):

$I(\mathrm{\&theta;},l)={\mathrm{<figure-callout\; id="0"\; label="I"\; filenames="A20061004186800191.png"\; state="\{\{state\}\}">I</figure-callout>}}_0{\&Integral;}_0^{V}P(V,E)\exp \{-{\&Integral;}_{L(\mathrm{\&theta;},l)}u(E,x,y)\mathrm{dxdy}\}\mathrm{dE}$

In the following formula, (V is when E) getting VkVp for the bulb power supply for P, the X line photon spectra of sending multiply by the normalized function of the absorption spectra of detector, u (E, x, y), but consider the relation that this parameter changes with the sub-ENERGY E of X linear light still for the mass attenuation coefficient of imaging object multiply by material density.After following formula is taken the logarithm, no longer satisfy above-mentioned linear projection model.With the common tissue-cortex of bone of human body, muscle and fat are example, the tissue of different-thickness is under the 120kVp voltage, funtcional relationship (supposed that the X line detector is the Na:CsI crystal after the attenuation data of wide range X linear light that bulb sends was taken the logarithm shown in figure (2), down together), bone tissue among the figure (2), muscle and fat are in the spectrum that 120kVp voltage bulb sends, the relation of damping capacity and thickness, wherein, thick line is represented bone tissue, below two fine dotted lines represent muscle and fat, the fine line of arch is represented the difference (having amplified 10 times) of bone tissue decay and linear relationship, and the two other fine dotted line has then reflected the situation (having amplified 10 times) of muscle and fat respectively.

It is form of straight lines that the linear projection model requires above-mentioned function, but in fact, this function does not satisfy its requirement.The result that so direct utilization reconstruction theory obtains is: for same tissue, the reconstruction numerical value in different spaces zone (so-called CT value) can occur different, because CT equipment provides main information to be reflected on the CT value, therefore this can bring vital error to diagnosis, formed CT value error that Here it is so-called " beam hardening effect (Beam Hardening) " is called as " the pseudo-shadow of beam hardening ", is one of pseudo-shadow of main reconstruction of CT equipment.

In the existing CT equipment, realize that the method for so-called beam hardening correction (Beam Hardening Correction) comes from engineering practice.Observe figure (2), under the wide range condition, though data for projection and thickness of tissue are not linear relationships, but be more or less the same with straight line, therefore, can consider to adopt a polynomial expression to describe this function, and realize the linearization correction according to this multinomial coefficient, consider correction accuracy, CT equipment adds certain thickness metal material (general equivalence is the aluminium of 35mm thickness) in the outlet of bulb, with the filtering as much as possible of lower energy photon in the spectrum, make the linearity of the X line projection data that penetrate human body quite a lot of as much as possible, carry out polynomial expression again and proofread and correct, therefore, obtain the beam hardening correction measure 1 of existing CT equipment.

Corrective action 1: after data for projection taken the logarithm, utilize a known polynomial expression to carry out linearization and proofread and correct, and this polynomial coefficient is to obtain by prior die body (water mould) experiment, owing to there are problems such as bulb wears out, needs regular calibration.

Notice shown in the figure (2), different tissue (bone skin tissue, soft tissue) has different correction multinomial coefficients, and at the data for projection timing, and do not know what of different tissues composition in the human body, therefore, for fear of the correction error that the human body different structure brings, existing C T equipment has adopted second kind of corrective action.

Corrective action 2: at bulb in the space of detector, except the shared space of tissue, a kind of filling material (resin) that is similar to the human body soft tissue decay characteristics is filled up in other space, like this, during the different angles projection, the all similar basically human body soft tissue (blood vessel of the material that ray penetrates, the attenuation properties and the soft tissue of tissue such as blood and muscle are suitable, the attenuation characteristic of adipose tissue is suitable with soft tissue basically, be the density difference), and the shared ratio of high density bone tissue is very little, therefore, can suppose that correction equation coefficient and projection angle and along slope coordinate (longitudinal) are irrelevant basically, so can adopt unified polynomial expression to carry out beam hardening correction, the error control of generation in allowed limits.

In fact, the beam hardening correction method of existing medicine CT is the combination of above-mentioned two kinds of measures.Obviously, there is following problem in such method.

The shortcoming that existing CT equipment beam hardening correction method exists:

1. each CT need carry out regular water mould experiment, to obtain the correction multinomial coefficient, increases and uses difficulty and cost.

2. in order to adapt to the difference that the partes corporis humani divides, need design die body respectively, experimentize, obtain different correction multinomial coefficients according to the position.

3. for same human body, there are individual difference and vertical (longitudinal) difference, therefore proofread and correct with the polynomial expression of unified coefficient, may have problems, most typical here is the pseudo-shadow of epidural in the Cranial Computed Tomography, still there is not desirable solution at present, in most CT equipment, need to judge differentiation by doctor's clinical experience, and in other number of C T equipment, the scheme that adopts post-processing approach to proofread and correct is arranged, but have the false risk of data reconstruction.

4. proofread and correct for the polynomial expression that adopts unified coefficient, and realize certain correction accuracy, CT equipment has to adopt filling material to make the data for projection of different angles even, and in fact this bring two shortcoming: A, the radiation dose that increases patient and the thermal capacity of bulb; B, uniform data for projection mean less dynamic range, and therefore, the corresponding increase of the noise of reconstruction has reduced the grey scale accuracy of rebuilding, and in order to reduce noise, has to increase X line dosage, and this has increased patient's radiation injury risk again.

Though the core of CT equipment is a reconstruction algorithm, in order to realize rebuilding accurately effect, numerous correcting algorithms is all kinds of CT equipment key core technology always.Except beam hardening correction, also have various N Reference Alignments, probe access gamma correction or the like.In these are proofreaied and correct, have only beam hardening correction relevant with the imaging object of scanning, other is then only relevant with equipment.Whether device-dependent correction no matter complexity all can accurately be carried out by a specific program.And the beam hardening correction relevant with imaging object if only adopt existing die body bearing calibration, can not be accomplished accurately so in theory.That is to say that this method is a kind of scheme of dealing with problems of through engineering approaches, lack strict theory support.

The applicant is in the further investigation process to CT equipment beam hardening correction, notice such fact: in the projection imaging of different angles, each coordinate points in the imaging object is passed by a plurality of rays, that is to say that the data for projection of each different angles has all comprised the information of this point, this means that the data for projection of different angles is correlated with.If there is a correction equation with limited parameter, so, may be able to solve the coefficient of correction equation by this correlativity, thereby realize the target of correction, and such correction comes from the data for projection of object itself, is adaptive, have better precision in theory, might avoid above-mentioned shortcoming.

The correlativity that is described as picture object space information and data for projection can be explained with the H-L condition for consistence.But what this condition was described is the situation of two-dimentional collimated beam projection, do not meet the projection situation of X line CT (XCT), the applicant has constructed the H-L condition for consistence formula of two-dimensional sector-shaped bundle, and on this basis, find the solution predefined nonlinear equation parameter, thereby reached the purpose of linearity correction.

Consider that existing result of study shows, the data for projection of XCT is proofreaied and correct and can be finished by polynomial expression, and therefore, this method can be used for the data for projection correction of wide range XCT in theory.More than work to give based on data for projection self correlativity undoubtedly and carry out the research work that the linearization of CT data is proofreaied and correct, reliable theoretical basis is provided.

Summary of the invention

The objective of the invention is to, a kind of CT data for projection beam hardening effect correction method based on the HL condition for consistence is provided, this method is applicable to the adjustment of data that all kinds of CT equipment are rebuild.

The technical solution that realizes the foregoing invention purpose is:

A kind of CT data for projection beam hardening effect correction method based on the HL condition for consistence, it is characterized in that, this method is according to equivalent beam hardening distortion correction model and H-L condition for consistence, projecting image data based on different angles, the correction that realizes data for projection of finding the solution by to equivalent beam hardening distortion correction model parameter specifically may further comprise the steps:

Step 1: the equivalent correction of non-linear distortions model that makes up data for projection

For medicine CT equipment, according to imaging physics model, the data for projection that observes is decomposed into the concrete model that high density tissue, low-density organize decay respectively to form, construct equivalent beam hardening distortion correction model:

$g(t,\mathrm{\&beta;})=P({\stackrel{\&RightArrow;}{\mathrm{\&alpha;}}}_k,r(t,\mathrm{\&beta;})f(t,\mathrm{\&beta;}))+P({\stackrel{\&RightArrow;}{\mathrm{\&eta;}}}_k,(1-r(t,\mathrm{\&beta;}))f(t,\mathrm{\&beta;}))\&CenterDot;\&CenterDot;\&CenterDot;\left(1\right)$

In the following formula, f (t, β) data for projection that obtains for detector, β is a projection angle, and t is the sampled point of data for projection, g (t, β) be the data for projection after proofreading and correct, (t β) represents that on the path of each bar projection ray, the high density bone tissue is to the contribution of projection values to r, span is (0,1), P () is the polynomial equation of linearity correction

Be the coefficient of polynomial equation, For the polynomial form of coefficient as follows:

$P({\stackrel{\&RightArrow;}{\mathrm{\&eta;}}}_k,x)={\mathrm{\&eta;}}_{1}x+{\mathrm{\&eta;}}_2{x}^2+{\mathrm{\&eta;}}_{1/2}{x}^{1/2}+{\mathrm{\&eta;}}_3{x}^3+{\mathrm{\&eta;}}_{1/3}{x}^{1/3}+\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;+{\mathrm{\&eta;}}_L{x}^L+{\mathrm{\&eta;}}_{1/L}{x}^{1/L}\&CenterDot;\&CenterDot;\&CenterDot;\left(2\right)$

Wherein,

For the polynomial form of coefficient with Identical;

Step 2: structure meets the expression formula of the H-L condition for consistence of CT projection mode

According to the HL consistance condition of known collimated beam projection, the expression formula that extends to circular scan track, fan-shaped projection pattern is:

$m_{i,k}=\&Integral;{\&Integral;}_C{x}^i{y}^{k}u(x,y)\mathrm{dxdy},i\&GreaterEqual;0,k\&GreaterEqual;0\&CenterDot;\&CenterDot;\&CenterDot;\left(3\right)$

$v_d\left(\mathrm{\&beta;}\right)={\&Integral;}_{-\&infin;}^{+\&infin;}{g}_f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\left(\frac{\mathrm{<figure-callout\; id="2"\; label="tD\; D"\; filenames="A20061004186800161.png"\; state="\{\{state\}\}">tD</figure-callout>}}{\sqrt{D^{}}}\right)}^d\left(\frac{D^3}{{({\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="A20061004186800161.png"\; state="\{\{state\}\}">D</figure-callout>}}^2+t^2)}^{3/2}}\right)\mathrm{dt}\&CenterDot;\&CenterDot;\&CenterDot;\left(4\right)$

$Q_d\left(\mathrm{\&beta;}\right)=\underset{r=0}{\overset{d}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}d\\ r\end{array}\right)m_{r,d-r}{\cos }^r\left(\mathrm{\&beta;}\right){\sin }^{d-r}\left(\mathrm{\&beta;}\right)\&CenterDot;\&CenterDot;\&CenterDot;\left(5\right)$

v _d(β)＝Q _d(β) (6)

Wherein, 0≤d≤N-1, and N is the number of projection angle;

In the above-listed formula, original faultage image be u (x, y), R is that it supports radius, C={ (x, y) | x ²+ y ²≤ R ², β is the projection angle of fan type bundle;

Step 3: according to equivalent correction of non-linear distortions model and H-L condition for consistence, the parameter that designs a model find the solution matrix, by calculating finding the solution of model parameter, thereby finish trimming process to data for projection.

Method of the present invention is suitable for the realization of the beam hardening correction function of all kinds of X line CT equipment.The beam hardening phenomenon is most important for application scenario meanings such as the diagnosis of many fine tissue pathologies and industry CT.Because beam hardening correction is relevant with the imaging object of scanning, among method of the present invention and the existing all kinds of CT, the method that only adopts die body to proofread and correct is compared, and can consider the difference of the beam hardening effect of different densities tissue respectively, automatically adapt to the difference of imaging object, have better correction accuracy.

Description of drawings

Fig. 1 is the X line spectrum functional digraph (top) that the CT bulb sends, and lower graph is through the spectral function (having amplified 5 times) after the 35mm aluminium filtration (according to the common way of CT equipment);

Fig. 2 is that bone tissue, muscle and fat are in the spectrum that 120kVp voltage bulb sends, the relation of damping capacity and thickness, wherein, thick line is represented bone tissue, below two fine dotted lines represent muscle and fat, the fine line of arch is represented the difference (having amplified 10 times) of bone tissue decay and linear relationship, and the two other fine dotted line has then reflected the situation (having amplified 10 times) of muscle and fat respectively.

Fig. 3 is the synoptic diagram of circular scan track, fan-shaped projected bundle, and β is a projection angle, and S is a radiographic source, g _f(t β) is data for projection, and O is a rotation center, and D is the distance of S to O;

Fig. 4 is the FORBILD head model of experiment usefulness, is used for the explanation of beam hardening correction, and wherein, bright part is a bone skin tissue, and grey color part is a soft tissue, and black partly is air.Projective parameter is: bulb voltage 120kVp, added filtration are the aluminium of 35mm thickness.Show among the figure and utilize the reconstructed image based on H-L condition for consistence bearing calibration acquisition that is proposed, method for reconstructing adopts FBP.In order to show the reconstruction effect before and after proofreading and correct, in follow-up figure, give the hatching line data of the indicated position of bright line in publishing picture.

Fig. 5 (a) is the locational contrast of proofreading and correct the front and back data reconstruction of horizontal hatching line.Article three, curve represent respectively through HL consistance correction equation proofread and correct result that the back rebuilds, without the reconstructed results of the adjustment of data, proofread and correct before and after the difference of data, the difference of data has been amplified 5 times before and after wherein proofreading and correct.

Fig. 5 (b) is the locational contrast of proofreading and correct the front and back data reconstruction of vertical hatching line.Article three, curve represent respectively through HL consistance correction equation proofread and correct result that the back rebuilds, without the reconstructed results of the adjustment of data, proofread and correct before and after the difference of data, the difference of data has been amplified 5 times before and after wherein proofreading and correct.

Fig. 6 utilizes the reconstructed image based on the acquisition of the H-L condition for consistence bearing calibration after cutting apart that is proposed, and method for reconstructing adopts FBP.

Fig. 7 (a) is the locational contrast of proofreading and correct the front and back data reconstruction of horizontal hatching line.Article three, curve represent respectively through formula (22) correction equation proofread and correct result that the back rebuilds, without the reconstructed results of the adjustment of data, proofread and correct before and after the difference of data.Wherein, r (β) cut apart back calculating by the precorrection reconstructed image and obtain, and the difference of data has been amplified 5 times before and after proofreading and correct by t.

The contrast of data reconstruction before and after locational correction of the vertical hatching line of Fig. 7 (b).Article three, curve represent respectively through formula (22) correction equation proofread and correct result that the back rebuilds, without the reconstructed results of the adjustment of data, proofread and correct before and after the difference of data.Wherein, r (β) cut apart back calculating by the precorrection reconstructed image and obtain, and the difference of data has been amplified 5 times before and after proofreading and correct by t.

The present invention is described in further detail below in conjunction with accompanying drawing.

Embodiment

1. make up the equivalent correction of non-linear distortions model of data for projection

For medicine CT equipment, according to imaging physics model, the data for projection that observes is decomposed into the concrete model that high density tissue, low-density organize decay respectively to form, construct equivalent beam hardening distortion correction model:

$g(t,\mathrm{\&beta;})=P({\stackrel{\&RightArrow;}{\mathrm{\&alpha;}}}_k,r(t,\mathrm{\&beta;})f(t,\mathrm{\&beta;}))+P({\stackrel{\&RightArrow;}{\mathrm{\&eta;}}}_k,(1-r(t,\mathrm{\&beta;}))f(t,\mathrm{\&beta;}))\&CenterDot;\&CenterDot;\&CenterDot;\left(1\right)$

In the following formula, f (t, β) data for projection that obtains for detector, β is a projection angle, and t is the sampled point of data for projection, g (t, β) be the data for projection after proofreading and correct, (t β) represents that on the path of each bar projection ray, the high density bone tissue is to the contribution of projection values to r, span is (0,1), P () is the polynomial equation of linearity correction

Coefficient for polynomial equation.

For the polynomial form of coefficient as follows:

$P({\stackrel{\&RightArrow;}{\mathrm{\&eta;}}}_k,x)={\mathrm{\&eta;}}_{1}x+{\mathrm{\&eta;}}_2{x}^2+{\mathrm{\&eta;}}_{1/2}{x}^{1/2}+{\mathrm{\&eta;}}_3{x}^3+{\mathrm{\&eta;}}_{1/3}{x}^{1/3}+\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;+{\mathrm{\&eta;}}_L{x}^L+{\mathrm{\&eta;}}_{1/L}{x}^{1/L}\&CenterDot;\&CenterDot;\&CenterDot;\left(2\right)$

And For the polynomial form of coefficient with

Similar:

$P({\stackrel{\&RightArrow;}{\mathrm{\&alpha;}}}_k,x)={\mathrm{\&alpha;}}_{1}x+{\mathrm{\&alpha;}}_2{x}^2+{\mathrm{\&alpha;}}_{1/2}{x}^{1/2}+{\mathrm{\&alpha;}}_3{x}^3+{\mathrm{\&alpha;}}_{1/3}{x}^{1/3}+\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;+{\mathrm{\&alpha;}}_L{x}^L+{\mathrm{\&alpha;}}_{1/L}{x}^{1/L}\&CenterDot;\&CenterDot;\&CenterDot;\left(3\right)$

2. the H-L condition for consistence expression formula under circular scan track, the fan-shaped projected bundle

Perspective view is asked for an interview figure (3), and the formula of derivation is as follows:

$m_{i,k}=\&Integral;{\&Integral;}_C{x}^i{y}^{k}u(x,y)\mathrm{dxdy},i\&GreaterEqual;0,k\&GreaterEqual;0\&CenterDot;\&CenterDot;\&CenterDot;\left(4\right)$

$v_d\left(\mathrm{\&beta;}\right)={\&Integral;}_{-\&infin;}^{+\&infin;}{g}_f(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\left(\frac{\mathrm{<figure-callout\; id="2"\; label="tD\; D"\; filenames="A20061004186800161.png"\; state="\{\{state\}\}">tD</figure-callout>}}{\sqrt{D^{}}}\right)}^d\left(\frac{D^3}{{({\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="A20061004186800161.png"\; state="\{\{state\}\}">D</figure-callout>}}^2+t^2)}^{3/2}}\right)\mathrm{dt}\&CenterDot;\&CenterDot;\&CenterDot;\left(5\right)$

$Q_d\left(\mathrm{\&beta;}\right)=\underset{r=0}{\overset{d}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}d\\ r\end{array}\right)m_{r,d-r}{\cos }^r\left(\mathrm{\&beta;}\right){\sin }^{d-r}\left(\mathrm{\&beta;}\right)\&CenterDot;\&CenterDot;\&CenterDot;\left(6\right)$

The H-L condition for consistence that popularization is then arranged: v _d(β)=Q _d(β) (7)

Wherein, 0≤d≤N-1, and N is the number of projection angle.

3. construct a kind of bearing calibration of general data for projection nonlinear distortion

Suppose that in the imaging process, X-ray machine rotates a circle around initial point, obtain N width of cloth cone-beam image altogether, image sequence can obtain fan-beam projection through pre-service

Wherein the β value is respectively β ₁, β ₂... β _NSuppose fan-beam projection under the β have M (the individual uniform sampling data of M＞N),

Discrete form can be expressed in matrix as:

Here Be the result after the nonlinear transformation, remember that real projection result is:

Do not consider The noise, then have:

${\hat{g}}_{i,j}=T(\mathrm{\&Theta;},g_{i,j})(1\&le;i\&le;M,1\&le;j\&le;N)\&CenterDot;\&CenterDot;\&CenterDot;\left(10\right)$

Wherein Θ is the K dimension parameter vector of nonlinear transformation T.

Investigate the H-L condition for consistence, have:

$\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_f(t_i,{\mathrm{\&beta;}}_j-\mathrm{arctg}\frac{t_i}{D}){\left(\frac{t_{\mathrm{<figure-callout\; id="2"\; label="i\; D\; D"\; filenames="A20061004186800161.png"\; state="\{\{state\}\}">i</figure-callout>}}D}{\sqrt{D^{}}}\right)}^d\left(\frac{D^3}{{({\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="A20061004186800161.png"\; state="\{\{state\}\}">D</figure-callout>}}^2+{t_i}^2)}^{3/2}}\right)$

$=\underset{r=0}{\overset{d}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}d\\ r\end{array}\right)m_{r,d-r}{\cos }^r\left({\mathrm{\&beta;}}_j\right){\sin }^{d-r}\left({\mathrm{\&beta;}}_j\mathrm{\&beta;}\right)(1\&le;j\&le;N,0\&le;d\&le;N-1)\&CenterDot;\&CenterDot;\&CenterDot;\left(11\right)$

If G, Θ, m _{R, d-r}(d≤N-1) as unknown number, total: MN+K+N (N+1)/2 variable, formula (10) then has MN+N with (11) simultaneous ²Individual equation, as long as K≤N (N-1)/2, then system of equations is separated, and that is to say to finish the gamma correction of image rotating sequence in the hope of Θ.

For two-dimensional slice image, suppose that the fan-beam projection that obtains is from the rotatable sequence image And real projection should be G.Because nonlinear transformation T is a monotonic increasing function, so T is reversible, has:

$g_{i,j}=T^{-1}\left({\hat{g}}_{i,j}\right)(1\&le;i\&le;M,1\&le;j\&le;N)\&CenterDot;\&CenterDot;\&CenterDot;\left(12\right)$

Here T ^-1The inverse transformation of expression T, it also is the nonlinear transformation of a monotone increasing.As long as we can determine T ^-1, just can finish rectification.Yet T ^-1Be the function of a form the unknown, by the approximation of function theory, it can approach with the combination of the basis function on certain function space, as long as exponent number is enough high, the effect of approaching is just enough good.

If S _kBe The function space of being opened, that is:

$S_k=\mathrm{span}\{1,x,{\mathrm{<figure-callout\; id="2"\; label="x"\; filenames="A20061004186800161.png"\; state="\{\{state\}\}">x</figure-callout>}}^2,\sqrt[2]{x},\&CenterDot;\&CenterDot;\&CenterDot;x^k,\sqrt[k]{x}\}\&CenterDot;\&CenterDot;\&CenterDot;\left(13\right)$

Then have: S ₂ S ₃ ... S _k ... S _∞(14)

So for any one nonlinear function T ^-1, can use S _kIn function approach, k is big more, it is good more to approach effect.That is to say T ^-1(x) can write:

$T^{-1}\left(x\right)=a_0+a_{1}x+a_2{x}^2+b_2\sqrt{x}+\&CenterDot;\&CenterDot;\&CenterDot;a_k{x}^k+b_k\sqrt[k]{x}+O\left(\min (x^k,\sqrt[k]{x})\right)\&CenterDot;\&CenterDot;\&CenterDot;\left(15\right)$

If can determine coefficient, just can reach the purpose of rectification.Because a is here known in the definition of T ₀=0.When 2k＜N-1, order:

$p(j,\mathrm{\&gamma;})=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}[{\hat{g}}_f(t_i,{\mathrm{\&beta;}}_j-\mathrm{arctg}\frac{t_i}{D}){]}^{\mathrm{\&gamma;}}\left(\frac{D^3}{{({\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="A20061004186800161.png"\; state="\{\{state\}\}">D</figure-callout>}}^2+{t_i}^2)}^{3/2}}\right)\&CenterDot;\&CenterDot;\&CenterDot;\left(16\right)$

X＝[a ₁?a ₂?b ₂?…?a _k?b _k] ^T (18)

E＝[1?1?1?…?1?1] ^T (19)

Wherein

Can obtain by interpolation, the H-L condition for consistence during according to d=0 can get:

AX＝m _0，0E (20)

Formula (20) is a 2k unit system of linear equations once, by least square method, can obtain its optimum solution:

X＝m _0，0(A ^TA) ^-1A ^TE (21)

As 2k＞N-1, can utilize d=0, the H-L constraint condition under 1 condition is similarly derived.In fact, general k very little (k≤4), the condition of 2k≤N-1 is easy to satisfy.As long as determined parameter, just rectification work can be finished by through type (15).Must be pointed out Coefficient m _0,0Though be uncertain, it can be absorbed by the grey scale mapping process.

4. construct a kind of beam hardening correction method based on the H-L condition for consistence

In order to obtain beam hardening effect correction method accurately, need to consider when same slice projection the difference of correction equation coefficient under the different projection angles, and the difference of the correction equation of vertical (longitudinal) different section.According to the physical model of imaging, in medical science X radiodiagnosis x spectral range, the human body material depends primarily on Compton scattering decay and photoelectric effect decay to the decay of ray, and degree of freedom wherein has only two dimension.This means that any material can replace with the linear combination of two kinds of different material decay the attenuation function of X linear light.Consider that the correction of the beam hardening effect of each material can proofread and correct with a polynomial function, therefore, we have at first constructed the correction equation that can accurately proofread and correct the beam hardening effect, and (this equation is at fan-beam projection, the correction equation of cone beam projection also can be constructed according to same principle), as formula (1).So-called trimming process is exactly the unknown quantity of finding the solution above-mentioned equation by the H-L condition for consistence:

R (t, β).

On the basis of die body shown in Figure 4, discuss in two kinds of situation.

1) .r (t, β)=0 o'clock bearing calibration

Here corresponding is the existing similar correcting scheme of CT equipment, promptly in the data for projection of same tomography, thinks that the data for projection of different angles can be proofreaied and correct with the polynomial expression of same coefficient.Fig. 4 is the image based on H-L condition for consistence bearing calibration acquisition that utilizes us to propose, and the method for reconstructing of employing is FBP, and adopts R-L type RAMP filter function.For calibration result is described, we select to show under two longitudinal and transverse hatching lines, the data contrast before and after proofreading and correct, and two bright lines wherein are used for the hatching line position of indicated number, and the result is as shown in Figure 5.

Fig. 5 has provided the contrast of the correction front and back data reconstruction of horizontal hatching line.Red solid line among the figure is represented the reconstructed results without beam hardening correction, and blue look solid line is represented the reconstructed results through beam hardening correction, the difference of data before and after blue look dotted line is represented to proofread and correct.The difference of data is as seen to the calibration result of the regional soft tissue in the centre soft tissue apparently higher than areas at both sides before and after proofread and correct, this has reflected on the one hand and has had typical beam hardening distortion phenomenon in the reconstructed image before proofreading and correct, and shows that on the other hand the bearing calibration of H-L condition for consistence plays calibration result to the beam hardening distortion of soft tissue area really.

Above-mentioned r (t, β)=0 o'clock correcting scheme is equivalent to the improvement of existing CT equipment calibration method, and its major advantage is as follows:

A. correction equation coefficient self-adaptation obtains, and accuracy is higher, and need not prior die body experiment, does not also just need regular CT parameter correction;

B. the correction equation coefficient is gone up adaptive change at along slope coordinate (longitudinal), and based on the feature of human body symmetry, this point can overcome the deficiency of classic method, will have using value preferably.

2). estimate that (t β), asks r

Bearing calibration

The key of this method is at first to obtain r, and (t β), on this basis, directly tries to achieve by the bearing calibration that we have constructed

(t β) can carry out: suppose that earlier (t β)=0, proofreaies and correct the back and rebuilds r like this to estimate r; Then in reconstructed image, the method that passing threshold is cut apart distinguishes highdensity bone sex organization zone with low-density soft tissue area; According to the physical features of result of cutting apart and different tissues estimate r (t, β); On this basis, proofread and correct and rebuild.

Fig. 6 is the reconstructed image based on image segmentation and H-L condition for consistence bearing calibration acquisition that utilizes the applicant to propose, and the method for reconstructing of employing is FBP, and adopts R-L type RAMP filter function.For calibration result is described, we select to show under two longitudinal and transverse hatching lines, the data contrast before and after proofreading and correct, and two bright lines wherein are used for the hatching line position of indicated number, and the result is as shown in Figure 7.

From Fig. 7 (a) as can be seen, inboard soft-tissue density value at left and right sides bone tissue has obtained obvious improvement, this point is better than the calibration result among Fig. 5 (a), and this illustrates that this method to the problem that pseudo-this CT equipment of shadow of epidural is difficult to overcome, will have solution approach preferably.

Claims (2)

1. CT data for projection beam hardening effect correction method based on the HL condition for consistence, it is characterized in that, this method is by making up the equivalent correction of non-linear distortions model of data for projection, structure meets the expression formula of the H-L condition for consistence of CT projection mode, according to equivalent beam hardening distortion correction model and H-L condition for consistence, projecting image data based on different angles, the correction that realizes data for projection of finding the solution by to equivalent beam hardening distortion correction model parameter specifically may further comprise the steps:

Step 1: the equivalent correction of non-linear distortions model that makes up data for projection

For medicine CT equipment, according to imaging physics model, the data for projection that observes is decomposed into the concrete model that high density tissue, low-density organize decay respectively to form, construct equivalent beam hardening distortion correction model:

$g(t,\mathrm{\&beta;})=P({\stackrel{\&RightArrow;}{\mathrm{\&alpha;}}}_k,r(t,\mathrm{\&beta;})f(t,\mathrm{\&beta;}))+P({\stackrel{\&RightArrow;}{\mathrm{\&eta;}}}_k,(1-r(t,\mathrm{\&beta;}))f(t,\mathrm{\&beta;}))---\left(1\right)$

In the following formula, f (t, β) data for projection that obtains for detector, β is a projection angle, and t is the sampled point of data for projection, g (t, β) be the data for projection after proofreading and correct, (t β) represents that on the path of each bar projection ray, the high density bone tissue is to the contribution of projection values to r, span is (0,1), P () is the polynomial equation of linearity correction Be the coefficient of polynomial equation,

For the polynomial form of coefficient as follows:

$P({\stackrel{\&RightArrow;}{\mathrm{\&eta;}}}_k,x)={\mathrm{\&eta;}}_{1}x+{\mathrm{\&eta;}}_2{x}^2+{\mathrm{\&eta;}}_{1/2}{x}^{1/2}+{\mathrm{\&eta;}}_3{x}^3+{\mathrm{\&eta;}}_{1/3}{x}^{1/3}+\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;\&CenterDot;+{\mathrm{\&eta;}}_L{x}^L+{\mathrm{\&eta;}}_{1/L}{x}^{1/L}---\left(2\right)$

Wherein,

For the polynomial form of coefficient with Identical;

Step 2: structure meets the expression formula of the H-L condition for consistence of CT projection mode

According to the HL consistance condition of known collimated beam projection, the expression formula that extends to circular scan track, fan-shaped projection pattern is:

m _i，k＝∫∫ _Cx ⁱy ^ku(x，y)dxdy?i≥0，k≥0 (3)

$v_d\left(\mathrm{\&beta;}\right)={\&Integral;}_{-\&infin;}^{+\&infin;}(t,\mathrm{\&beta;}-\mathrm{arctg}\frac{t}{D}){\left(\frac{\mathrm{tD}}{\sqrt{D^2+t^2}}\right)}^d\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}---\left(4\right)$

$Q_d\left(\mathrm{\&beta;}\right)=\underset{r=0}{\overset{d}{\mathrm{\&Sigma;}}}\left(\begin{array}{c}d\\ r\end{array}\right)m_{r,d-r}{\cos }^r\left(\mathrm{\&beta;}\right){\sin }^{d-r}\left(\mathrm{\&beta;}\right)---\left(5\right)$

v _d(β)＝Q _d(β) (6)

Wherein, 0≤d≤N-1, and N is the number of projection angle.

In the above-listed formula, original faultage image be u (x, y), R is that it supports radius, C={ (x, y) | x ²+ y ²≤ R ², β is the projection angle of fan type bundle;

Step 3: according to equivalent correction of non-linear distortions model and H-L condition for consistence, the parameter that designs a model find the solution matrix, by calculating finding the solution of model parameter, thereby finish trimming process to data for projection.

2. the method for claim 1 is characterized in that: the unknown parameter of finding the solution equivalent correction of non-linear distortions model

R (t β), can form following two kinds of bearing calibrations, and concrete grammar is as follows:

2.1 (t, β)=0, all the soft tissue with human body is suitable promptly to suppose to penetrate on the path exhausted big most of material by ray, finds the solution by descending column matrix and operational method at first to suppose r Obtain the beam hardening correction result of method I:

In the CT imaging process, X-ray machine rotates a circle around initial point can obtain fan-beam projection

Obtain through after the discretize

1≤i≤M wherein, 1≤j≤N;

$p(j,\mathrm{\&gamma;})=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{\left[{\hat{g}}_f(t_i,{\mathrm{\&beta;}}_j-\mathrm{arctg}\frac{t_i}{D})\right]}^{\mathrm{\&gamma;}}\left(\frac{D^3}{{(D^2+{t_i}^2)}^{3/2}}\right)---\left(7\right)$

X＝[a ₁?a ₂?b ₂?…?a _k?b _k] ^T (9)

E＝[1?1?1?…?1?1] ^T (10)

Wherein Can obtain by interpolation, the H-L condition for consistence during according to d=0 has:

AX＝m _0，0E (11)

X is a unknown quantity, can try to achieve the optimum solution of X according to many method of value solving;

X promptly

As long as determined this parameter,, must be pointed out Coefficient m just rectification work can be finished by through type (2) _0,0Though be uncertain, it can be absorbed by the grey scale mapping process;

2.2 on the basis as a result that bearing calibration I obtains, finish the reconstruction of CT image, use again following method obtain r (t, result of calculation β):

Suppose that the distortion of CT image beam hardening is mainly caused by the bone material, to the CT image that has obtained:

At first, cutting apart (perhaps other are through well-designed image partition method) through simple threshold binary image, is bone zone and non-bone zone with the CT image segmentation;

Then, along (t, β) direction with the numerical value in non-bone zone sue for peace Φ (t, β), this process also often is called re-projection, computing formula is:

$\mathrm{\&Phi;}(t,\mathrm{\&beta;})=\underset{i=1}{\overset{N_I}{\mathrm{\&Sigma;}}}{f}_I(x_i,y_i){\mathrm{\&Psi;}}_{(t,\mathrm{\&beta;})}(x_i,y_i)---\left(12\right)$

N wherein _IBe the pixel sum of whole C T reconstructed image, f _IBe the resulting CT reconstructed results of bearing calibration I, Ψ _{(t, β)}Be the space subordinate function, promptly

${\mathrm{\&Psi;}}_{(t,\mathrm{\&beta;})}(x,y)=\left\{\begin{array}{c}1,(x,y)\&Element;(M\&cap;l(t,\mathrm{\&beta;}))\\ 0,(x,y)\&NotElement;(M\&cap;l(t,\mathrm{\&beta;}))\end{array}\right.---\left(13\right)$

Soft tissue area after wherein M represents to cut apart, l (t, β) the specific X ray of expression path;

At last, calculate

$r(t,\mathrm{\&beta;})=\{{\hat{g}}_f(t,\mathrm{\&beta;})-\mathrm{\&Phi;}(t,\mathrm{\&beta;})\}/{\hat{g}}_f(t,\mathrm{\&beta;});---\left(14\right)$

2.3 (t on basis β), by following operational method, finds the solution at known r With Obtain the method II correction result that formula (1) is described:

At first, according to the method described in the step 2.1, with For projection is found the solution

(t β) finds the solution for projection with Φ

According to formula (1), incite somebody to action then

Correction result and Φ (t, correction result addition β) obtains the projection after overcorrection;

Projection after overcorrection rebuild obtain image through the CT of overcorrection;

At last, and the method renewal r in the employing step 2.2 (t, β);

Iteration is finished above step, converges to a determined value up to the result.

Images