CN103150744B - A kind of X ray multi-power spectrum CT data for projection process and image rebuilding method - Google Patents

Abstract

The invention discloses a kind of X ray multi-power spectrum CT data for projection process and image rebuilding method.Mainly comprise X-ray energy spectrum CT projection sinogram disposal route and the acceleration of iterative convergence reconstruction algorithm based on compressed sensing.X-ray energy spectrum CT projection sinogram disposal route mainly comprises two aspects: 1. suppress vertical wire artifact in projection sinogram; 2. high brightness noise in projection sinogram is removed.Acceleration of iterative convergence reconstruction algorithm based on compressed sensing combines based on the minimized optimization constraint condition of image total variation (TV) with algebraic reconstruction technique (OS-SART) while of order subset.Because current X-ray energy spectrum CT detection system (X-ray energy resolving photon digital detector) also comes with some shortcomings, there is more noise and artifact in the data for projection of acquisition.This method is utilizing preprocessing means to carry out effectively pretreated to X ray multi-power spectrum CT data for projection while, OS-SART algorithm based on TV is introduced in X ray multi-power spectrum CT image reconstruction, accelerate Image Iterative convergence, inhibit the noise in reconstruction image and artifact well.

Description

A kind of X ray multi-power spectrum CT data for projection process and image rebuilding method

Technical field

The invention belongs to X ray multi-power spectrum CT (Computed Tomography) image reconstruction and Digital Image Processing research field, relate to a kind of X ray multi-power spectrum CT data for projection process and image rebuilding method.

Background technology

X ray computer Tomography (X-ray Computed Tomography, X-CT) technology, that imaging is carried out to a tangent plane (tomography) of three-dimensional body, the Density Distribution situation of measured object section can be gone out by lossless detection, intuitively represent internal structure of body situation and material composition with image format.Because X-CT technology has huge advantage in material detection, the context of detection such as biomedicine, Aero-Space apparatus, geology archaeological samples, military project weapon, bridge dyke building and radioactive contamination are widely used in.

For tradition or conventional X-CT system, what its x-ray source produced is polychrome (multi-energy) X ray (X ray continuum), and what its detector adopted is the integrated detection mode of X-ray energy, that is, it is received the x-ray photon of different-energy entirety.But from physics angle analysis, the X ray of different-energy has different attenuation characteristics, and these attenuation characteristics, the different physical properties of checking matter material but can be reflected, such as, after X ray and checking matter effect, the K-edge characteristic etc. of material.So when after the incident X-rays and object to be detected effect of different-energy, the signal entrained by transmission X-ray, can present the cognitive information that checking matter enriches more.Due to traditional X-CT system, the ability of X-ray energy do not differentiated by its detector, but the x-ray photon of different-energy integrated reception (or claiming integral measurement technology), the average attenuation characteristic of its reflection X ray.Thus, this not only causes the loss of X-ray attenuation information, is unfavorable for the judgement to checking matter material physical properties, and especially for medical science X-CT technology, the CT image after its reconstruction, is difficult to the contrast difference distinguishing different soft-tissue imaging.There is a kind of new CT technology based on X-ray energy resolving photon counting detection technology in recent years---X ray multi-power spectrum CT technology, it is the information entrained by transmission X-ray after utilizing the incident X-rays of different-energy and object to be detected effect and carries out the technology of Computed tomography, the X-ray attenuation characteristic of more horn of plenty can be represented, provide the information being conducive to differentiating substance characteristics more.The development of X ray multi-power spectrum CT technology to X-CT technology has milestone significance, and it is a kind of CT new technology, is also a study hotspot in current CT field.

X ray multi-power spectrum CT technology, is in the initial stage of research and development abroad at present, domesticly not yet really starts to walk to the research work.Existing X ray multi-power spectrum CT detection system (X-ray energy resolving photon digital detector) also comes with some shortcomings, share by Compton scattering, electric charge and the impact of the factor such as pile-up effect, cause the X ray multi-power spectrum CT data for projection of acquisition to there is more noise and artifact.In addition, X ray multi-power spectrum CT often carries out imaging in specific X-ray energy scope, and specific X-ray energy range detection to photon number be limited, cause in projected image and there is certain quantum noise.Current X ray Multi-energy-spectruCT CT image reconstruction method has continued to use conventional CT image reconstruction algorithm, in most of experimental study, mainly utilize filtered back projection (Filtered Backprojection in analytic reconstruction method, FBP) algorithm rebuilds X ray multi-power spectrum CT image, and the maximum deficiency of FBP algorithm is anti-noise jamming ability.Based on this, in order to improve X ray multi-power spectrum CT image reconstruction quality, while utilizing some preprocessing means process data for projection, also need to introduce the strong image rebuilding method of anti-noise jamming ability.

Therefore, design and a kind ofly effectively can suppress noise in X ray multi-power spectrum CT data for projection and artifact, reconstruct the CT image of better quality, just become the problem that the present invention pays close attention to.

Summary of the invention

What the present invention needed to solve is how to suppress the noise in X ray multi-power spectrum CT data for projection and artifact and how to improve the X ray multi-power spectrum CT image quality issues of reconstruction.Consider the feature that X ray multi-power spectrum CT data for projection noise is large, the object of this invention is to provide a kind of effective X ray multi-power spectrum CT data for projection process and image rebuilding method, the noise in X ray multi-power spectrum CT data for projection and artifact can be suppressed well, and reconstruct the CT image of better quality.

In order to reach above object, the present invention adopts following technical scheme:

This method comprises X ray multi-power spectrum CT projection sinogram process and walks greatly based on the acceleration of iterative convergence reconstruction two of compressed sensing; Described projection sinogram process comprises vertical wire artifact suppression and high brightness noise in sinogram and removes two steps, and vertical wire artifact suppresses the mode adopting frequency domain filtering, and high brightness noise removes the mode adopting spatial domain identification, filtering; It is by based on image total variation (Total Variation that acceleration of iterative convergence based on compressed sensing is rebuild, TV) minimized optimization constraint condition and order subset algebraic reconstruction technique (Ordered Subset Simultaneous Algebraic ReconstructionTechnique simultaneously, OS-SART) combine, be divided into inside and outside two-layer iterative loop step, outer iteration loop performs OS-SART iterative reconstruction algorithm, and internal layer iterative loop performs the minimization process of the total variation (TV) of rebuilding image f.

The concrete steps of this method are as follows:

The process of step 1:X ray multi-power spectrum CT projection sinogram

1. vertical wire artifact in projection sinogram is suppressed: use vertical wire artifact in frequency domain filtering method offset of sinusoidal figure to carry out filtering process;

Suppose detector often row have X probe unit, the sinogram under Y scanning angle is s (k, l), wherein k=1 ..., X, l=1 ..., Y, the Fourier transform of sinogram is expressed as

$S(u,v)=\frac{1}{\mathrm{XY}}\underset{k=0}{\overset{X-1}{\mathrm{\Σ}}}\underset{l=0}{\overset{Y-1}{\mathrm{\Σ}}}s(k,l)e^{-2\mathrm{\πj}(\mathrm{uk}/X+\mathrm{vl}/Y)}$

In sinogram, vertical wire artifact produces intensity values at v=0 place, defective pixel produces peak value at higher horizontal frequency u place, v and u represents image frequency conversion S (v in above formula respectively, u) two independents variable: vertical frequency and horizontal frequency, utilize wave filter by the intensity values filtering at these frequency places, suppress vertical wire artifact in sinogram; Butterworth lowpass filters selected by described wave filter, and its conversion expression formula is

$H(u,v)=\{\begin{array}{c}\frac{1}{1+{\left(\frac{u}{u_0}\right)}^{2n}},\mathrm{if}\left|v\right|\≤v_0\\ 1,\mathrm{otherwise}\end{array}$

Wherein, n=4, u ₀=8/N △ x, v ₀=1/M △ x, △ x is the width of each probe unit of detector here.

2. remove high brightness noise in projection sinogram: first utilize between noise and neighbor that gray difference amount is to identify these noises, if noise is f (i, j) in sinogram, then the quadratic sum of noise and adjacent 8 pixel grey scale differences is

D=(f(i,j)-f(i-1,j-1)) ²+(f(i,j)-f(i-1,j)) ²+(f(i,j)-f(i-1,j+1)) ²

+(f(i,j)-f(i,j-1)) ²+(f(i,j)-f(i,j+1)) ²+(f(i,j)-f(i+1,j-1)) ²

+(f(i,j)-f(i+1,j)) ²+(f(i,j)-f(i+1,j+1)) ²

Here the gray difference amount between noise and neighbor that defines is

$c=\sqrt{D/8}$

If when measures of dispersion c is greater than certain value σ, can judge that this pixel is as noise, and give this pixel by adjacent for this pixel 8 pixel grey scale mean values, wherein σ chooses according to actual sinogram gamma characteristic, namely chooses according to non-zero pixels average gray value in sinogram.

Step 2: the OS-SART based on TV rebuilds

Under CT Image Iterative reconstruction framework, CT imaging system expression formula is as follows

Af=b

Wherein, b=(b ¹..., b ^m) ∈ R ^mrepresent data for projection, f=(f ₁..., f _m) ∈ R ⁿimage is rebuild in representative, A=(a _ij) represent reconstruction iteration matrix;

1. outer OS-SART iterative loop

SART iterative algorithm following expression:

$f_j^{(n+1)}=f_j^{\left(n\right)}+\frac{1}{a_{+j}}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\frac{a_{\mathrm{ij}}}{a_{i+}}(b^i-A^i{f}^{\left(n\right)}),n=\mathrm{0,1},\·\·\·$

Wherein, i=1 ..., M, represents matrix A i-th row element sum, j=1 ..., N, represent matrix A jth column element sum, n is iterations;

Suppose that the set of all projection sequence numbers is

B={1,···,M}

The set of all projection sequence numbers can be divided into T subclass, and the set expression of the sequence number that projects in each subset is

$B_t=\{i_1^t,\·\·\·,i_{M\left(t\right)}^t\},t=1,\·\·\·,T$

Therefore, the set expression of all projection sequence numbers is:

$B=\{1,\·\·\·,M\}=\underset{1\≤t\≤T}{\mathrm{\∪}}{B}_t$

SART reconstruction algorithm then based on order subset has following expression:

$f_j^{(n+1)}=f_j^{\left(n\right)}+\underset{{i\∈B}_{\left[n\right]}}{\mathrm{\Σ}}\frac{a_{\mathrm{ij}}}{a_{+j}}\frac{b_i-A^i{f}^{\left(n\right)}}{a_{i+}},n=\mathrm{0,1},\·\·\·;$

2. internal layer minimizes iterative loop based on TV

Rebuild minimizing of total variation (TV) of image f

$\underset{f}{\min }{\left|\right|\▿f\left|\right|}_1$

Wherein,

${\left|\right|\▿f\left|\right|}_1=\underset{i,j}{\mathrm{\Σ}}{d}_{i,j},d_{i,j}=\sqrt{{(f_{i,j}-f_{i+1,j})}^2+(f_{i,j}-f_{i,j+1}{)}^2}$

Here be defined as the total variation (TV) of rebuilding image f, f _i,jfor rebuilding a grey scale pixel value of image f, d _i,jit is a discrete gradient;

Utilize gradient descent method, formula is as follows:

f ^(m+1)=f ^(m)-λ′ωυ

Wherein, λ ' is Gradient Descent relaxation factor, be a gradient direction, for Gradient Descent scale parameter, m is interior loop iterations, and 1. outer OS-SART iterative loop n is corresponding with step 2 for n.

The invention has the beneficial effects as follows:

1, utilize frequency domain filtering method to carry out filtering process to X ray multi-power spectrum CT projection sinogram, effectively can suppress vertical wire artifact in sinogram.

2, to utilize between noise pixel and neighbor gray difference amount to identify high brightness noise, and neighbor average gray value is replaced the original gray-scale value of noise, effectively can remove high bright spot noise in sinogram.

3, by the OS-SART algorithm application based on TV in X ray multi-power spectrum CT image reconstruction, this algorithm combines based on the minimized optimization constraint condition of TV and OS-SART algorithm, while acceleration of iterative convergence, inhibit the noise in reconstruction image and artifact well.

Accompanying drawing explanation

Fig. 1 is embodiment of the present invention process flow diagram;

Fig. 2 is the projection sinogram of an embodiment of the present invention chicken wings section;

Fig. 3 is the sinogram after the vertical wire artifact of embodiment of the present invention correction;

Fig. 4 is the sinogram after embodiment of the present invention correction high brightness noise;

Fig. 5 is that the embodiment of the present invention utilizes the X ray multi-power spectrum CT image rebuild based on the OS-SART algorithm of TV.

Embodiment

Mode below by embodiment further illustrates the present invention, does not therefore limit the present invention among described scope of embodiments.

Below in conjunction with accompanying drawing, illustrate technical conceive of the present invention, and the course of work under this design:

According to Fig. 1, the inventive method comprises the following steps: (1) corrects wire artifact vertical in projection sinogram, utilizes frequency domain filtering method to carry out filtering process to X ray multi-power spectrum CT projection sinogram, suppresses vertical wire artifact in sinogram; (2) high brightness noise in projection sinogram is corrected, to utilize between noise pixel and neighbor gray difference amount to identify high brightness noise, and neighbor average gray value is replaced the original gray-scale value of noise, thus remove high bright spot noise in sinogram; (3) utilize the OS-SRAT algorithm based on TV to rebuild X ray multi-power spectrum CT image, suppress further to rebuild the noise in image and artifact.

Concrete steps are as follows:

(1) vertical wire artifact in projection sinogram is suppressed

X ray multi-power spectrum CT explorer portion probe unit defect (bad pixel or dead pixel) often causes in sinogram exists some vertical wire artifacts, and these vertical wire artifacts cause to rebuild in image and occur a large amount of ring artifact.Therefore, before CT image reconstruction, need to carry out pretreatment operation, suppress vertical wire artifact in sinogram.The present invention utilizes vertical wire artifact in a kind of frequency domain filtering method offset of sinusoidal figure to carry out filtering process, and the method is described below:

Suppose detector often row have X probe unit, the sinogram under Y scanning angle is s (k, l), wherein k=1 ..., X, l=1 ..., Y, the Fourier transform of sinogram is expressed as

$S(u,v)=\frac{1}{\mathrm{XY}}\underset{k=0}{\overset{X-1}{\mathrm{\Σ}}}\underset{l=0}{\overset{Y-1}{\mathrm{\Σ}}}s(k,l)e^{-2\mathrm{\πj}(\mathrm{uk}/X+\mathrm{vl}/Y)}$

In sinogram, vertical wire artifact produces intensity values at v=0 place, and defective pixel produces peak value at higher horizontal frequency u place, utilizes wave filter by the intensity values filtering at these frequency places, effectively can suppress vertical wire artifact in sinogram.Here Butterworth lowpass filters selected by wave filter, and its conversion expression formula can be written as

$H(u,v)=\{\begin{array}{c}\frac{1}{1+{\left(\frac{u}{u_0}\right)}^{2n}},\mathrm{if}\left|v\right|\≤v_0\\ 1,\mathrm{otherwise}\end{array}$

Wherein, n=4, u ₀=8/N △ x, v ₀=1/M △ x, △ x is the width of each probe unit of detector here.In an experiment, according to wire artifact feature vertical in sinogram, filter parameter can be adjusted.

(2) high brightness noise in projection sinogram is removed

By the impact of random noise, often occur the noise of high brightness in sinogram, these high brightness noises cause to rebuild in image and occur some bar shaped artifacts.Therefore, before CT image reconstruction, need to remove high brightness noise in sinogram.In true sinogram, the isolated distribution of single pixel often of high brightness noise.First the present invention utilizes between noise and neighbor that gray difference amount is to identify these noises, and suppose that in sinogram, noise is f (i, j), then the quadratic sum of noise and adjacent 8 pixel grey scale differences is

D=(f(i,j)-f(i-1,j-1)) ²+(f(i,j)-f(i-1,j)) ²+(f(i,j)-f(i-1,j+1)) ²

+(f(i,j)-f(i,j-1)) ²+(f(i,j)-f(i,j+1)) ²+(f(i,j)-f(i+1,j-1)) ²

+(f(i,j)-f(i+1,j)) ²+(f(i,j)-f(i+1,j+1)) ²

Here the gray difference amount between noise and neighbor that defines is

$c=\sqrt{D/8}$

If when measures of dispersion c is greater than certain value σ, can judge that this pixel is as noise, and give this pixel by adjacent for this pixel 8 pixel grey scale mean values, wherein σ chooses according to actual sinogram gamma characteristic.

(3) based on the OS-SART reconstruction algorithm of TV

Under CT Image Iterative reconstruction framework, CT imaging system can be write as following expression

Af=b

Wherein, b=(b ¹..., b ^m) ∈ R ^mrepresent data for projection, f=(f ₁..., f _m) ∈ R ⁿimage is rebuild in representative, A=(a _ij) represent reconstruction iteration matrix.Algebraic reconstruction technique (ART) is the earliest for the iterative algorithm of CT image reconstruction, and the iteration form of this algorithm can be expressed as

$f_j^{(n+1)}=f_j^{\left(n\right)}+{\mathrm{\λ}}_n\frac{a_{\mathrm{ij}}}{{\left|\right|A^i\left|\right|}^2}(b^i-A^i{f}^{\left(n\right)}),n=\mathrm{0,1},\·\·\·$

Wherein, n represents iterations, and i=nmod (M)+1 is the index of equation, for the euclideam norm of Iterative Matrix A i-th row.

ART algorithm, in reconstruction image process, is easily subject to salt-pepper noise impact.Reduce relaxation factor λ can affect by restraint speckle, but iterative convergence speed slows down.SART algorithm, by ART algorithm development, it maintains ART algorithm the convergence speed, has good noise inhibiting ability simultaneously.SART iterative algorithm can be write as following expression

$f_j^{(n+1)}=f_j^{\left(n\right)}+\frac{1}{a_{+j}}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\frac{a_{\mathrm{ij}}}{a_{i+}}(b^i-A^i{f}^{\left(n\right)}),n=\mathrm{0,1},\·\·\·$

Wherein, , i=1 ..., M, represents matrix A i-th row element sum, j=1 ..., N, represents matrix A jth column element sum.

Suppose that the set of all projection sequence numbers is

B={1,…,M}

The set of all projection sequence numbers can be divided into T subclass, and the set expression of the sequence number that projects in each subset is

$B_t=\{i_1^t,\·\·\·,i_{M\left(t\right)}^t\},t=1,\·\·\·,T$

Therefore, the set of all projection sequence numbers also can be expressed as:

$B=\{1,\·\·\·,M\}=\underset{1\≤t\≤T}{\mathrm{\∪}}{B}_t$

SART reconstruction algorithm then based on order subset can have following expression:

$f_j^{(n+1)}=f_j^{\left(n\right)}+\underset{{i\∈B}_{\left[n\right]}}{\mathrm{\Σ}}\frac{a_{\mathrm{ij}}}{a_{+j}}\frac{b_i-A^i{f}^{\left(n\right)}}{a_{i+}},n=\mathrm{0,1},\·\·\·$

In actual CT image reconstruction Study on Problems, the inner a lot of constituent of inspected object has identical or approximate attenuation characteristic, and in CT image, just reflect the identical or close of gray scale, therefore image can rarefaction.Wherein, gradient conversion is a kind of conventional sparse transformation, and the TV of image is through commonly using the l of its gradient image ₁norm represents.Iterative approximation problem based on TV can be converted into following optimization problem

$\underset{f}{\min }{\left|\right|\▿f\left|\right|}_1,s.t.\mathrm{Af}=b$

Wherein,

${\left|\right|\▿f\left|\right|}_1=\underset{i,j}{\mathrm{\Σ}}{d}_{i,j},d_{i,j}=\sqrt{{(f_{i,j}-f_{i+1,j})}^2+(f_{i,j}-f_{i,j+1}{)}^2}$

Here be defined as the total variation (TV) of rebuilding image f, f _i,jfor rebuilding a grey scale pixel value of image f, d _i,jit is a discrete gradient.

The Image Iterative reconstruction algorithm performed based on TV can be divided into inside and outside two iterative loop steps.Outer iteration loop performs OS-SART iterative reconstruction algorithm, and internal layer iterative loop performs the minimization process of the total variation (TV) of rebuilding image f.When performing internal layer iterative loop, can gradient descent method be utilized, can be expressed as

f ^(m+1)=f ^(m)-λ′ωυ

Wherein, λ ' is Gradient Descent relaxation factor, be a gradient direction, for Gradient Descent scale parameter, m is interior loop iterations.

Embodiment:

The present embodiment utilizes a X ray multi-power spectrum CT to scan a specific X-ray energy scope (25 ~ 40keV) chicken wings that is marked with contrast medium (iodine solution).

Fig. 2 is the projection sinogram of an embodiment chicken wings section, occurs some vertical wire artifacts because detector defect (bad pixel or dead pixel) causes in this sinogram.Utilize step (1) described method to carry out filtering process to the sinogram in Fig. 2, suppress the vertical wire artifact in sinogram.

Fig. 3 is the sinogram after the vertical wire artifact of embodiment correction, but the high brightness noise still having some fragmentary intersperses among in sinogram, utilizes the high brightness noise in the described method offset of sinusoidal figure of step (2) to identify, and removes these high brightness noises.

Fig. 4 is the sinogram after embodiment correction high brightness noise, and the high brightness noise in sinogram is well suppressed.Utilize step (3) described method to carry out CT image reconstruction to the sinogram in the sinogram in Fig. 2, the sinogram in Fig. 3 and Fig. 4 respectively, the OS-SART iterative reconstruction algorithm execution wherein based on TV can divide following step: 1. input measurement data for projection b and initial pictures f=0;

2. utilize OS-SART algorithm to upgrade and rebuild image;

3. utilize gradient descent method that reconstruction image total variation (TV) is minimized;

2. and 3. 4. step is repeated, till meeting and rebuilding convergence constraint condition.

In execution based in the OS-SART iterative reconstruction algorithm process of TV, Gradient Descent relaxation factor λ ' is that to minimize iterations m be 30, OS-SART image reconstruction iterations n to 0.2, TV is 20.

Fig. 5 is that embodiment utilizes the X ray multi-power spectrum CT image rebuild based on the OS-SART algorithm of TV, wherein Fig. 5 a is for utilizing initial sinusoids figure (Fig. 2) reconstruction CT image out, containing more ring artifact in this CT image, this is because the vertical wire artifact of some in initial sinusoids figure causes.Fig. 5 b utilizes to revise the sinogram (Fig. 3) after vertical wire artifact and rebuild CT image out, and containing some strip artifacts in this CT image, this causes due to high brightness noise in sinogram.Fig. 5 c utilizes the reconstruction of the sinogram (Fig. 4) after revising high brightness noise CT image out, through the process to wire artifact vertical in initial sinusoids figure and high brightness noise, utilize and can suppress noise in CT image reconstruction and artifact further based on the OS-SART iterative reconstruction algorithm of TV, obtain good reconstructed results.

Claims (2)

1. X ray multi-power spectrum CT data for projection process and an image rebuilding method, comprises X ray multi-power spectrum CT projection sinogram process and walks greatly based on the acceleration of iterative convergence reconstruction two of compressed sensing; Described projection sinogram process comprises vertical wire artifact suppression and high brightness noise in sinogram and removes two steps, and vertical wire artifact suppresses the mode adopting frequency domain filtering, and high brightness noise removes the mode adopting spatial domain identification, filtering; It is combine based on the minimized optimization constraint condition of image total variation TV with algebraic reconstruction technique OS-SART while of order subset that acceleration of iterative convergence based on compressed sensing is rebuild, be divided into inside and outside two-layer iterative loop step, outer iteration loop performs OS-SART iterative reconstruction algorithm, and internal layer iterative loop performs the minimization process of the total variation TV rebuilding image f.

2. X ray multi-power spectrum CT according to claim 1 data for projection process and image rebuilding method, is characterized in that: the concrete steps of this method are as follows:

The process of step 1:X ray multi-power spectrum CT projection sinogram

1. vertical wire artifact in projection sinogram is suppressed: use vertical wire artifact in frequency domain filtering method offset of sinusoidal figure to carry out filtering process;

Suppose detector often row have X probe unit, the sinogram under Y scanning angle is s (k, l), wherein k=1 ..., X, l=1 ..., Y, the Fourier transform of sinogram is expressed as

$S(u,v)=\frac{1}{XY}\underset{k=0}{\overset{X-1}{\Σ}}\underset{l=0}{\overset{Y-1}{\Σ}}s(k,l)e^{-2\mathrm{\π}j(uk/X+vl/Y)}$

In sinogram, vertical wire artifact produces intensity values at v=0 place, defective pixel produces peak value at higher horizontal frequency u place, v and u represents S (u in above formula respectively, v) two independents variable: vertical frequency and horizontal frequency, utilize wave filter by the intensity values filtering at these frequency places, suppress vertical wire artifact in sinogram; Butterworth lowpass filters selected by described wave filter, and its conversion expression formula is

$H(u,v)=\left\{\begin{array}{c}\frac{1}{1+{\left(\frac{u}{u_0}\right)}^{2n}},if\left|v\right|\≤v_0\\ 1,otherwise\end{array}\right.$

Wherein, n=4, u ₀=8/X Δ x, v ₀=1/Y Δ x, Δ x is the width of each probe unit of detector here;

2. remove high brightness noise in projection sinogram: first utilize between noise and neighbor that gray difference amount is to identify these noises, if noise is f (i, j) in sinogram, then the quadratic sum of noise and adjacent 8 pixel grey scale differences is

D＝(f(i,j)-f(i-1,j-1)) ²+(f(i,j)-f(i-1,j)) ²+(f(i,j)-f(i-1,j+1)) ²

+(f(i,j)-f(i,j-1)) ²+(f(i,j)-f(i,j+1)) ²+(f(i,j)-f(i+1,j-1)) ²

+(f(i,j)-f(i+1,j)) ²+(f(i,j)-f(i+1,j+1)) ²

Here the gray difference amount between noise and neighbor that defines is

$c=\sqrt{D/8}$

If when measures of dispersion c is greater than certain value σ, can judge that this pixel is as noise, and give this pixel by adjacent for this pixel 8 pixel grey scale mean values, wherein σ chooses according to actual sinogram gamma characteristic, namely chooses according to non-zero pixels average gray value in sinogram;

Step 2: the OS-SART based on TV rebuilds

Under CT Image Iterative reconstruction framework, CT imaging system is expressed as:

Af＝b

Wherein, b=(b ¹..., b ^m) ∈ R ^mrepresent data for projection, f=(f ₁..., f _m) ∈ R ⁿimage is rebuild in representative, A=(a _ij) represent reconstruction iteration matrix;

1. outer OS-SART iterative loop

SART iterative algorithm following expression:

${f_j}^{(n+1)}={f_j}^{\left(n\right)}+\frac{1}{a_{+j}}\underset{i=1}{\overset{M}{\Σ}}\frac{a_{ij}}{a_{i+}}(b^i-A^i{f}^{\left(n\right)}),n=0,1,...$

Wherein, i=1 ..., M, represents matrix A i-th row element sum, represent matrix A jth column element sum, n is iterations;

Suppose that the set of all projection sequence numbers is

B＝{1,…,M}

The set of all projection sequence numbers can be divided into T subclass, and the set expression of the sequence number that projects in each subset is

$B_t=\{i_1^t,...,i_{M\left(t\right)}^t\},t=1,...,T$

Therefore, the set expression of all projection sequence numbers is:

$B=\{1,...,M\}=\underset{1\≤t\≤T}{\∪}{B}_t$

SART reconstruction algorithm then based on order subset has following expression:

${f_j}^{(n+1)}={f_j}^{\left(n\right)}+\underset{i\∈B_{\[n\]}}{\Σ}\frac{a_{ij}}{a_{+j}}\frac{b_i-A^i{f}^{\left(n\right)}}{a_{i+}},n=0,1,...$

2. internal layer minimizes iterative loop based on TV

Rebuild minimizing of the total variation TV of image f

$\underset{f}{\min }{||\▿f||}_1$

Wherein,

${||\▿f||}_1=\underset{i,j}{\Σ}{d}_{i,j},d_{i,j}=\sqrt{{(f_{i,j}-f_{i+1,j})}^2+{(f_{i,j}-f_{i,j+1})}^2}$

Here be defined as the total variation (TV) of rebuilding image f, f _i,jfor rebuilding a grey scale pixel value of image f, d _i,jit is a discrete gradient;

Utilize gradient descent method, formula is as follows:

f ^(m+1)＝f ^(m)-λ′ωυ

Wherein, λ ' is Gradient Descent relaxation factor, be a gradient direction, for Gradient Descent scale parameter, m is interior loop iterations.