CN104778667A - Level-set-based correction method for cupping artifact in cone-beam CT - Google Patents

Abstract

The invention discloses a level-set-based correction method for a cupping artifact in cone-beam CT. The level-set-based correction method is applied to the correction of a cone-beam CT slice image; the correction of the cupping artifact in the cone-beam CT can be carried out in a self-adaptive manner, and automatically completed without manual intervention; repeated scanning of an object to be detected is not required; the complexity of a cone-beam CT system is not increased; the level-set-based correction method is aimed at the reconstructed slice image and can be directly oriented to users, and the correction work can be completed without any change in original equipment of the original cone-beam CT; through the method, the correction of the cupping artifact in the cone-beam CT can be efficiently carried out, and the image contrast can also be improved.

Description

A kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set

Technical field

The present invention relates to a kind of bearing calibration using cupping artifact in the Cone-Beam CT sectioning image of level set algorithm, belong to technical field of image processing.

Background technology

Cone-Beam CT has good sweep velocity and radiation utilization factor, and the load that effectively can reduce X-ray tube exports, and reduces scanning cost, also can obtain high-resolution three-dimension faultage image data fast.

At present, cone of influence beam CT reconstruction picture quality a lot of because have, as: x-ray scattering, noise, geometric error, power spectrum, probe unit non_uniform response etc.But because the dull and stereotyped CT of cone-beam uses large-scale X-ray flat panel detector, this makes image quality be more vulnerable to the impact of X ray scattering and beam hardening compared with traditional CT.The Analysis and judgments had a strong impact on rebuilding image such as inaccurate of the artifact formed because of scattering and beam hardening, CT number.For in the cone-beam CT reconstruction image of human body, these artifact main manifestations are CT value unevenness artifact, the artifact mostly showing as cup-shaped of this artifact.These artifacts affect very serious for the visual display aspect based on threshold value and the pyramidal CT image segmentation aspect based on threshold value.Therefore, the correction for Cone-Beam CT cupping artifact seems particularly necessary.

In order to reduce the impact of cupping artifact (that is: CT value unevenness artifact), at present, the cupping artifact because scattering causes mainly is considered in most bearing calibration, and the research of prior art or documents and materials mainly concentrates on the scatter correction on projected image.These methods can be divided into two classes: a class is based on software, another kind of is based on hardware.The method of great majority based on software is based on Monte-Carlo Simulation, and it is that one corrects very effective method to cone beam computed tomography (CT) scattering.But it is consuming time especially.In recent years, people have proposed some Monte Carlo simulation algorithms improved, as: based on the method etc. of GPU.But even if adopt fast algorithm, heavy calculating still hinders the application of its reality.Many hardware based bearing calibrations are also good, such as: part ray blocks, elementary modulation etc.

The shortcoming of prior art mainly comprises:

(1) prior art is mainly for projection image correction, does not have the direct correction for rebuilding rear sectioning image.

(2) prior art great majority concentrate in the method for the artifact correction caused because of scattering, and mostly need to add hardware device, as: patent 200710019084 and 201310039298, these two patents all need to add hardware device on the Cone-Beam CT equipment of costliness, add the complicacy of operation and cause potential security risk to equipment.Particularly patent 200710019084 needs twice sweep testee, the radiant quantity of the measured object increased so undoubtedly.

In sum, in the method for prior art or documents and materials, Monte-carlo Simulation Method expends time in very much, the structure that result is limited to modulation panel self is corrected in primary ray modulator approach, based on partial dispersion radionetric survey method, some needs to increase exposure dose, and some method is not high to the accuracy of estimation of scatter distributions.And the present invention can solve problem above well.

Summary of the invention

The object of the invention is to propose a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set, and the method is applied to Cone-Beam CT sectioning image and corrects.The cupping artifact that the method adaptively can carry out Cone-Beam CT corrects, and just automatically can complete correction without the need to manual intervention.The method does not need multiple scanning testee; Do not increase the complexity of cone-beam CT system; For the sectioning image after reconstruction, can direct user oriented, do not need to make any change to the existing equipment of original Cone-Beam CT, just can complete correction work, the cupping artifact that the method can carry out Cone-Beam CT efficiently corrects, and can also improve the contrast of image simultaneously.

The present invention solves the technical scheme that its technical matters takes: a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set, the method comprises the steps:

Step 1: obtain the Cone-Beam CT slice of data with artifact;

Step 2: calculate overall flow function ε (φ, f according to formula (7) and formula (8) _s, p) with level set equation of constraint F (φ, p, f _s);

Step 3: to the Cone-Beam CT slice of data with artifact, according to formula (10), fixing p and f _s, use finite difference method iteration to develop

Step 4: according to formula (11), fixing φ and f _s, calculate the estimated value of p variate-value

Step 5: according to formula (12), fixing φ and p, calculates f _sthe estimated value of variate-value

Step 6: if do not restrain or reach and do not arrive iterations, then make and get back to step 3;

Step 7: according to formula (13), calculates the Cone-Beam CT sectioning image after correcting.

Effective effect:

1, the present invention can directly correct for the cupping artifact of the sectioning image after reconstruction.

2, calculated amount of the present invention is relatively little, while can carrying out the correction of Cone-Beam CT sectioning image cupping artifact efficiently, also can improve the contrast of image well.

3, the present invention can directly towards CT section demand user, and not needing makes any change to original Cone-Beam CT existing equipment just can complete correction work.

Accompanying drawing explanation

Fig. 1 is method flow diagram of the present invention.

Fig. 2 is the schematic diagram of image slice.

Identifier declaration: (a) represents CTP486 original slice figure; B () represents the CTP486 slice map after using the inventive method to correct; C () represents CTP Top original slice figure; D () represents the CTP Top slice map after using the inventive method to correct.

Fig. 3 is the horizontal sectional drawing of image slice.

Identifier declaration: the horizontal sectional drawing of row in the middle of (a), (b) represent.

Fig. 4 is CTP486 cupping artifact index τ _cupthe schematic diagram of zoning.

Fig. 5 is CTP Top cupping artifact index τ _cupzoning.

Fig. 6 is the schematic diagram of image (right side) after image (left side) and gray scale inhomogeneity correction before gray scale inhomogeneity correction.

Fig. 7 is 1 dimension horizontal sectional drawing of Fig. 6 head Cone-Beam CT section.

Fig. 8 is the experiment of breast Cone-Beam CT sectioning image gray scale inhomogeneity correction, and left column is image (a, c) before correction, the schematic diagram of image (b, d) after the right side is classified as and corrects.

Embodiment

Below in conjunction with Figure of description, the invention is described in further detail.

As shown in Figure 1, a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set, the method comprises the steps:

Step 1: obtain the Cone-Beam CT slice of data with artifact;

Step 2: calculate overall flow function ε (φ, f according to formula (7) and formula (8) _s, p) with level set equation of constraint F (φ, p, f _s);

Step 3: to the Cone-Beam CT slice of data with artifact, according to formula (10), fixing p and f _s, use finite difference method iteration to develop

Step 4: according to formula (11), fixing φ and f _s, calculate the estimated value of p variate-value

Step 5: according to formula (12), fixing φ and p, calculates f _sthe estimated value of variate-value

Step 6: if do not restrain or reach and do not arrive iterations, then make and get back to step 3;

Step 7: according to formula (13), calculates the Cone-Beam CT sectioning image after correcting.

Bearing calibration of the present invention specifically comprises:

Rebuild sectioning image decomposition to comprise:

In Cone-Beam CT, reconstruction algorithm is based on FDK algorithm, and rebuilding image set f can be write as

$\begin{array}{c}f=\frac{1}{{4\mathrm{\π}}^2}{\∫}_0^{2\mathrm{\π}}\frac{{d_{\mathrm{so}}}^2}{{(d_{\mathrm{so}}+r\·s)}^2}{\∫}_{-\∞}^{\∞}\frac{d_{\mathrm{so}}}{{({d_{\mathrm{so}}}^2+t^2+z^2)}^{1/2}}\\ \·I_{3D}(t,z\left(r\right),\mathrm{\φ})\·{\∫}_{-\∞}^{\∞}\left|\mathrm{\ω}\right|e^{j2\mathrm{\π\ω}(t\left(r\right)-t)}\mathrm{d\ωdtd\φ}\end{array}---\left(1\right)$

Wherein d _soexpression source to the distance of turning axle, I _3D(t, z (r), φ) represents the sequence of projected image.Here, projected image I _3Dbe decomposed into as follows:

I _3D＝P _3D+S _3D+n (2)

Wherein P _3Dthe true picture of never artifact, S _3Dbe the artifact sections caused by scattering and beam hardening, n is noise.P _3Dbe the intrinsic physical attribute of object slice object, if there is material in N, so can suppose that it is divided into the constant region of N number of gray scale (that is: CT value).According to the observation, the present invention can find Cone-Beam CT section artifact S _3Dit is the cupping artifact of slowly change.Rewrite reconstruction image formula as follows,

$\begin{array}{c}f=\frac{1}{{4\mathrm{\π}}^2}{\∫}_2^{2\mathrm{\π}}\frac{{d_{\mathrm{so}}}^2}{{(d_{\mathrm{so}}+r\·s)}^2}{\∫}_{-\∞}^{\∞}\frac{d_{\mathrm{so}}}{{({d_{\mathrm{so}}}^2+t^2+z^2)}^{1/2}}\·\\ (P_{3D}(t,z\left(r\right),\mathrm{\φ})+S_{3D}(t,z\left(r\right),\mathrm{\φ})+n)\·{\∫}_{-\∞}^{\∞}\left|\mathrm{\ω}\right|e^{j2\mathrm{\π\ω}(t\left(t\right)-t)}\mathrm{d\ωdtd\φ}\\ =f_p+f_s+f_n\end{array}---\left(3\right)$

Wherein

$\begin{array}{c}{f}_s=\frac{1}{{4\mathrm{\π}}^2}{\∫}_0^{2\mathrm{\π}}\frac{{d_{\mathrm{so}}}^2}{{(d_{\mathrm{so}}+r\·s)}^2}{\∫}_{-\∞}^{\∞}\frac{d_{\mathrm{so}}}{{({d_{\mathrm{so}}}^2+t^2+z^2)}^{1/2}}\\ \·S_{3D}(t,z\left(r\right),\mathrm{\φ})\·{\∫}_{-\∞}^{\∞}\left|\mathrm{\ω}\right|e^{j2\mathrm{\π\ω}(t\left(r\right)-t)}\mathrm{d\ωdtd\φ}\end{array},$

That is rebuild image and can be expressed as three independent components additions.

Definition neighborhood self-energy function comprises:

The present invention considers the region Θ of the circle of a radius ρ _y∈ Ω, and define Θ _y={ x:|x-y|≤ρ }, y ∈ Ω.Whole territory sectioning image is designated as Ω, and Ω can be divided into N number of subregion, is respectively due to f _san image slowly changed equally, so at circular domain Θ _yin, the value f of all X _sx () is close to value f _s(y).Therefore, at each subregion Θ _y∩ Ω _i, intensity f _p(x)+f _sx () is close to coefficient p _i+ f _s(y).

Therefore, for region Θ _y∩ Ω _ithe present invention has:

f(x)≈f _s(y)+p _i， (4)

To neighborhood Θ _yin, the present invention defines energy function and is

${\mathrm{\ϵ}}_{{\mathrm{\Θ}}_y}=\underset{{\mathrm{\Ω}}_i}{\∫}K(y-x){|\mathrm{ff}\left(x\right)-(f_s\left(y\right)+p_i)|}^2\mathrm{dx},---\left(5\right)$

Wherein α is a normaliztion constant, makes ∫ K (u)=1, and δ is the standard variance (or scale parameter) of Gaussian function, and ρ is Θ _ythe radius of circular domain.

Definition global energy function and level set function comprise:

For whole image-region, energy function can be expressed as:

$\begin{array}{c}\mathrm{\ϵ}=\∫{\mathrm{\ϵ}}_{{\mathrm{\Θ}}_y}\mathrm{dy}\\ =\∫\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{{\mathrm{\Ω}}_i}{\∫}K(y-x){|f\left(x\right)-(f_s\left(y\right)+p_i)|}^2\mathrm{dx}\right)\mathrm{dy}\end{array}---\left(6\right)$

The present invention uses level set to represent energy function, and for the level set of N=2, other multilevel collection situations can similarly be derived.For simplicity, the vectorial p=(p of the present invention ₁... p _n) represent constant p _1,p _n.Therefore, level set function φ, vectorial p, and cup-shaped f _sbe the variable of energy ε, therefore ε can be write as ε (φ, f _s, following form p):

$\begin{array}{c}\mathrm{\ϵ}(\mathrm{\φ},f_s,p)=\∫(\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\∫K(y-x){|f\left(x\right)-(f_s\left(y\right)+p_i)|}^2{M}_i\left(\mathrm{\φ}\left(x\right)\right)\mathrm{dx})\mathrm{dy}\\ =\∫(\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\∫K(y-x){|f\left(x\right)-(f_s\left(y\right)+p_i)|}^2\mathrm{dy})M_i\left(\mathrm{\φ}\left(x\right)\right)\mathrm{dx}\end{array},---\left(7\right)$

Wherein, M ₁(φ)=H (φ) and M ₂(φ)=1-H (φ), H is Heaviside function, and φ is a level set function, and its zero level integrates as C ₀={ x: φ (x)=0}, its component image field is two disjoint Ω ₁={ x: φ (x) > 0} and Ω ₂={ x: φ (x) < 0} region.For 2 level sets, C ₀={ x: φ (x)=0} and C _k={ x: φ _kx ()=0} is then divided into three regions, i.e. N=3 image.

In order to retrain zero level set function, the present invention adds bound term to whole level set function.New height collection equation is defined as:

F(φ,p,f _s)＝ε(φ,f _s,p)+vL(φ)+μR _q(φ) (8)

Wherein be used for retraining the slickness of zero level outline line, H is Heaviside function,

$R_p\left(\mathrm{\&phi;}\right)=\&Integral;q(\&dtri;\mathrm{\&phi;})\mathrm{dx}$ Be used for maintaining symbolic distance attribute, wherein $q=\frac{{(s-1)}^2}{2}.$

7.2.4 predictor φ, p, f _s

Fixing p and f _s, minimize F (φ, p, f _s) use normal gradients descending method, separate gradient current equation obtain:

$\begin{array}{c}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{{\&PartialD;M}_i}{\&PartialD;\mathrm{\&phi;}}{e}_i+\mathrm{v\&delta;}\left(\mathrm{\&phi;}\right)\mathrm{div}\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)\\ +\mathrm{\&mu;}\mathrm{div}(\mathrm{dq}(\&dtri;\mathrm{\&phi;})\&dtri;\mathrm{\&phi;})\end{array},---\left(9\right)$

Here e _i=∫ K (y-x) | f (x)-(f _s(y)+p _i) | ²dy, be terraced operator, div () is divergence operator, and δ is Dirac delta function.

We use finite difference method iteration to develop φ, and result is used represent.Method is described below: establish for the value of φ discretize, for the right formula approximate value of formula (9).Then have and then obtain iterative formula:

${\mathrm{\&phi;}}_{i,j}^{k+1}={\mathrm{\&phi;}}_{i,j}^k+\mathrm{\&Delta;tL}\left({\mathrm{\&phi;}}_{i,j}^k\right),---\left(10\right)$

Wherein make c ₀for being greater than the constant of 0.

Fixing φ and f _s, use be worth after representing the optimization of p.Minimize ε (φ, f _s, p), then $\frac{\&PartialD;\mathrm{\&epsiv;}(\mathrm{\&phi;},f_s,p)}{\&PartialD;p}=0,$ Obtain

${\hat{p}}_i=\frac{\&Integral;(\&Integral;\left(K(y-x)(f\left(x\right)-f_s\left(y\right))\right)\mathrm{dy})M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}}{\&Integral;(\&Integral;K(y-x)\mathrm{dy})M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}},---\left(11\right)$

Fixing φ and p.Minimize ε (φ, f _s, p), then obtain

${\hat{f}}_s=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{\&Integral;(\&Integral;\left(K(y-x)(f\left(x\right)-p_i)\right)\mathrm{dy})M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}}{\&Integral;(\&Integral;K(y-x)\mathrm{dy})M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}}---\left(12\right)$

By process of iteration, all minimize F (φ, p, f each time _s).Use last estimated result to carry out iteration in current iteration, carry out until Data Convergence or till reaching iterations.Finally, the present invention obtains correcting rear image and is:

$f_p=f-{\hat{f}}_s---\left(13\right)$

Experimentation of the present invention and result specifically comprise:

Quantitative test index definition comprises:

The present invention defines cupping artifact τ _cup=100 (u _{m, edge}-u _{m, center})/u _{m, edge}, wherein u _{m, center}and u _{m, edge}it is the CT value (that is: gray-scale value) of die body center and peripheral.

Root mean square contrast is expressed as wherein I _ijit is two dimensional image (i, j) position pixel value.

Contrast signal to noise ratio (S/N ratio) degree CNR computing formula is CNR=|u _{m, 1}-u _{m, 2}|/σ _m, wherein u _{m, 1}, u _{m, 2}be the average of two contrast districts, σ _{m, 1}, σ _{m, 2}for the standard deviation of contrast district, pixel noise σ _mbe defined as σ _m=(σ _{m, 1}+ σ _{m, 2})/2.

Catphan 500 die body Cone-Beam CT section cupping artifact corrects experiment.

Experiment of the present invention uses Catphan500 die body to test, and uses CTP486 and die body top (that is: CTP Top) in die body.Can find out from Fig. 2, the present invention can eliminate cupping artifact and be difficult to perceive level to human eye; What represent in figure 3 is the 1 dimension horizontal cross-section in corresponding diagram 2, therefrom can find out that the present invention significantly can eliminate the cupping artifact of sectioning image.To scatter correction effect quantitatively analysis in table 1 and table 2.Can obtain from table, the method that the present invention proposes makes cupping artifact τ _cupaverage decline about 91.8%.

Quantitative test before and after table 1:CTP486 sectioning image corrects, the figure that cuts before correction is expressed as PI_NONE, and the slice map after correction is expressed as PI_SC

Quantitative test before and after table 2:CTP Top sectioning image corrects, the figure that cuts before correction is expressed as PI_NONE, and the slice map after correction is expressed as PI_SC.

The Cone-Beam CT section cupping artifact of human skull corrects experiment

Cone-Beam CT cupping artifact for human skull corrects front and back slice map as shown in Figure 6.Fig. 7 is the middle row sectional view of Fig. 6 image.There is gray scale inequality (on the left of Fig. 7) in the cerebral tissue before not correcting, shows as the rising of fringe region brightness.Uniform gray level after correcting, reaches ideal effect.

Breast Cone-Beam CT cupping artifact corrects experiment

As shown in Figure 8, breast Cone-Beam CT experiment display method of the present invention can be reduced to the imperceptible state of human eye gray scale unevenness artifact, and the present invention improves the contrast of image.RMSC after correction is about 1.3 times before correcting.

Claims (8)

1. based on a bearing calibration for cupping artifact in the Cone-Beam CT of level set, it is characterized in that, described method comprises the steps:

Step 1: obtain the Cone-Beam CT slice of data with artifact;

Step 2: calculate overall flow function ε (φ, f according to formula (7) and formula (8) _s, p) with level set equation of constraint F (φ, p, f _s);

Step 3: to the Cone-Beam CT slice of data with artifact, according to formula (10), fixing p and f _s, use finite difference method iteration to develop

Step 4: according to formula (11), fixing φ and f _s, calculate the estimated value of p variate-value

Step 5: according to formula (12), fixing φ and p, calculates f _sthe estimated value of variate-value

Step 6: if do not restrain or reach and do not arrive iterations, then make and get back to above-mentioned steps 3;

Step 7: according to formula (13), calculates the Cone-Beam CT sectioning image after correcting.

2. a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set according to claim 1, is characterized in that, the formula (7) of described step 2 is:

$\begin{array}{c}\mathrm{\&epsiv;}(\mathrm{\&phi;},f_s,p)=\&Integral;(\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}\&Integral;K(y-x)|f\left(x\right)-(f_s\left(y\right)+p_i){|}^2{M}_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx})\mathrm{dy}\\ =\&Integral;(\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}\&Integral;K(y-x)|f\left(x\right)-(f_s\left(y\right)+p_i){|}^2\mathrm{dy})M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}\end{array};$

The formula (8) of described step 2 is:

F(φ,p,f _s)＝ε(φ,f _s,p)+vL(φ)+μR _q(φ)。

3. a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set according to claim 1, is characterized in that, the formula (10) of described step 3 is:

${\mathrm{\&phi;}}_{i,j}^{k+1}={\mathrm{\&phi;}}_{i,j}^k+\mathrm{\&Delta;tL}\left({\mathrm{\&phi;}}_{i,j}^k\right).$

4. a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set according to claim 1, is characterized in that, the formula (11) of described step 4 is:

${\hat{p}}_i=\frac{\&Integral;(\&Integral;\left(K(y-x)(f\left(x\right)-f_s\left(y\right))\right)\mathrm{dy})M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}}{\&Integral;(\&Integral;K(y-x)\mathrm{dy})M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}}.$

5. a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set according to claim 1, is characterized in that, the formula (12) of described step 5 is:

${\hat{f}}_s=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\frac{\&Integral;(\&Integral;\left(K(y-x)(f\left(x\right)-p_i)\right)\mathrm{dy})M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}}{\&Integral;(\&Integral;K(y-x)\mathrm{dy})M_i\left(\mathrm{\&phi;}\left(x\right)\right)\mathrm{dx}}.$

6. a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set according to claim 1, is characterized in that, the formula (13) of described step 7 is:

$f_p=f-{\hat{f}}_s.$

7. a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set according to claim 1, is characterized in that, the adaptive cupping artifact carrying out Cone-Beam CT of described method corrects, and just automatically can complete correction without the need to manual intervention; Described method does not need multiple scanning testee; Do not increase the complexity of cone-beam CT system; For the sectioning image after rebuilding, can direct user oriented, do not need to make any change to the existing equipment of original Cone-Beam CT, just can complete correction work.

8. a kind of bearing calibration based on cupping artifact in the Cone-Beam CT of level set according to claim 1, is characterized in that, described method is applied to Cone-Beam CT sectioning image and corrects.