CN103150744A - X-ray multi-energy spectrum computed tomography (CT) projection data processing and image reconstruction method - Google Patents

Abstract

The invention discloses an X-ray multi-energy spectrum computed tomography (CT) projection data processing and image reconstruction method, which mainly comprises an X-ray energy spectrum CT projection sinogram processing method and a compressed sensing-based accelerated iterative convergent reconstruction algorithm. The X-ray energy spectrum CT projection sinogram processing method mainly comprises the following steps of: (1) restraining vertical streaking artifacts in a projection sinogram; (2) removing high-brightness noisy points in the projection sinogram. The compressed sensing-based accelerated iterative convergent reconstruction algorithm refers to that image total variation (TV) minimization-based optimal constraint conditions and the ordered-subset simultaneous algebraic reconstruction techniques (OS-SART) are combined. As a number of defects still exist in a traditional X-ray energy spectrum CT detection system (X-ray energy resolution photon counting detector), more noise and artifacts exist in acquired projection data. According to the X-ray multi-energy spectrum CT projection data processing and image reconstruction method, X-ray multi-energy spectrum CT projection data are effectively preprocessed by utilizing a preprocessing means, and meanwhile, a TV-based OS-SART algorithm is introduced into X-ray multi-energy spectrum CT image reconstruction, and therefore, image iterative convergence is accelerated, and the noise and the artifacts in a reconstructed image is well restrained.

Description

A kind of X ray multi-power spectrum CT data for projection is processed 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 and process and image rebuilding method.

Background technology

X ray computer tomoscan imaging (X-ray Computed Tomography, X-CT) technology, that a tangent plane (tomography) to three-dimensional body carries out imaging, the Density Distribution situation that can lossless detection goes out the measured object section intuitively represents the internal structure of body situation and material forms with image format.Because the X-CT technology is having huge advantage aspect the material detection, the context of detection such as biomedicine, Aero-Space apparatus, geology archaeological samples, military project weapon, bridge dyke building and radioactive contamination have been 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 to receive the x-ray photon of different-energy is whole.Yet from the physics angle analysis as can be known, the X ray of different-energy is that different attenuation characteristics is arranged, and these attenuation characteristics, but can reflect the different physical properties of checking matter material, for example, after X ray and checking matter effect, the K-edge characteristic of material etc.So after the incident X-rays and object to be detected effect of different-energy, the signal that the Transmission X ray is entrained can present the more abundant cognitive information of checking matter.Due to traditional X-CT system, its detector is not differentiated the ability of X ray energy, receives (or claiming the integral measurement technology) but the x-ray photon of different-energy is integrated, the average attenuation characteristic of its reflection X ray.Thus, the loss that this has not only caused the X ray dampening information is unfavorable for the judgement to the checking matter material physical properties, and especially for medical science X-CT technology, the CT image after it is rebuild is difficult to distinguish the contrast difference of different soft-tissue imagings.A kind of new CT technology based on X ray energy-resolved photon counting Detection Techniques has appearred in recent years---X ray multi-power spectrum CT technology, it is to utilize the incident X-rays of different-energy and the entrained information of the Transmission X ray after the object to be detected effect and the technology of carrying out Computed tomography, can represent the more X ray attenuation characteristic of horn of plenty, the information that are conducive to differentiate substance characteristics are provided more.X ray multi-power spectrum CT technology has milestone significance to the development of X-CT technology, and it is a kind of CT new technology, is also a study hotspot in present CT field.

X ray multi-power spectrum CT technology abroad is in the initial stage of research and development at present, domestic the research work is not yet really started to walk.Existing X ray multi-power spectrum CT detection system (X ray energy-resolved photon counting detector) also comes with some shortcomings, shared and the impact of the factor such as pile-up effect by Compton scattering, electric charge, cause the X ray multi-power spectrum CT data for projection that obtains to have more noise and pseudo-shadow.In addition, X ray multi-power spectrum CT often carries out imaging in specific X ray energy range, and the photon number that detects in specific X ray energy range is limited, causes and has certain quantum noise in projected image.X ray multi-power spectrum CT image rebuilding method has been continued to use the traditional CT image reconstruction algorithm at present, in most of experimental studies, mainly utilize (the Filtered Backprojection of filtered back projection in the analytic reconstruction method, FBP) algorithm is rebuild X ray multi-power spectrum CT image, and the deficiency of FBP algorithm maximum is the anti-noise jamming ability.Based on this, in order to improve X ray multi-power spectrum CT image reconstruction quality, when utilizing some preprocessing means to process data for projection, also need to introduce the strong image rebuilding method of anti-noise jamming ability.

Therefore, design a kind of noise and pseudo-shadow that can effectively suppress in X ray multi-power spectrum CT data for projection, reconstruct the CT image of better quality, just become institute of the present invention problems of concern.

Summary of the invention

What the present invention need to solve is the X ray multi-power spectrum CT image quality issues that how to suppress noise and the pseudo-shadow in X ray multi-power spectrum CT data for projection and how to improve reconstruction.Consider the large characteristics of X ray multi-power spectrum CT data for projection noise, the purpose of this invention is to provide a kind of effective X ray multi-power spectrum CT data for projection processes and image rebuilding method, can suppress well noise and pseudo-shadow in X ray multi-power spectrum CT data for projection, and reconstruct the CT image of better quality.

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

This method comprises that X ray multi-power spectrum CT projection sinogram is processed and rebuild for two steps greatly based on the acceleration of iterative convergence of compressed sensing; Described projection sinogram is processed and is comprised that in sinogram, the pseudo-shadow of vertical wire suppresses and high brightness noise two steps of removal, and the pseudo-shadow of vertical wire suppresses to adopt the mode of frequency domain filtering, and the mode of spatial domain identification, filtering is adopted in the removal of high brightness noise; Rebuilding based on the acceleration of iterative convergence of compressed sensing is with based on image total variation (Total Variation, TV) minimized optimization constraint condition and order subset while algebraic reconstruction technique (Ordered Subset Simultaneous Algebraic Reconstruction Technique, OS-SART) combine, be divided into inside and outside two-layer iterative loop step, outer iteration loop is carried out the OS-SART iterative reconstruction algorithm, and the internal layer iterative loop is carried out the minimization process of the total variation (TV) of reconstructed image f.

The concrete steps of this method are as follows:

Step 1:X ray multi-power spectrum CT projection sinogram is processed

1. suppress in projection sinogram the pseudo-shadow of vertical wire: use in frequency domain filtering method offset of sinusoidal figure the pseudo-shadow of vertical wire to carry out filtering and process;

Suppose that the every row of detector has X probe unit, the sinogram under Y scanning angle is s (k, l), k=1 wherein ..., 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, the pseudo-shadow of vertical wire is in v=0 place's generation intensity values, defective pixel is at higher horizontal frequency u place's generation peak value, v and u represent respectively picture frequency conversion S (v in following formula, u) two independents variable: vertical frequency and horizontal frequency, utilize wave filter with the intensity values filtering at these frequencies place, suppress the pseudo-shadow of vertical wire in sinogram; Described wave filter is selected the fertile hereby low-pass filter of Bart, 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: at first utilize between noise and neighbor the gray difference amount to identify these noises, establish that in sinogram, noise is f (i, j), the poor quadratic sum of noise and adjacent 8 pixel grey scales 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 define that between noise and neighbor, the gray difference amount is

$c=\sqrt{D/8}$

When if measures of dispersion c is worth σ greater than certain, can judge that this pixel is noise, and give this pixel with the adjacent 8 pixel grey scale mean values of this pixel, 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 ^NRepresent reconstructed image, A=(a _ij) represent the 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 row element sum,

J=1 ..., N represents matrix A j column element sum, n is iterations;

The set of supposing all projection sequence numbers is

B={1,···,M}

The set of all projection sequence numbers can be divided into T son set, and in every subset, the set expression of projection sequence number 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$

The SART reconstruction algorithm 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

The total variation of reconstructed image f (TV) minimizes

$\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 reconstructed image f, f _i,jBe the grey scale pixel value of reconstructed image f, d _i,jIt is a discrete gradient;

Utilize gradient descent method, formula is as follows:

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

Wherein, λ ' is the Gradient Descent relaxation factor,

Be a gradient direction,

Be the Gradient Descent scale parameter, m is the interior loop iterations, and in n and step 2,1. outer OS-SART iterative loop n is corresponding.

The invention has the beneficial effects as follows:

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

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

3, will be based on the OS-SART algorithm application of TV in X ray multi-power spectrum CT image reconstruction, this algorithm is to combine based on the minimized optimization constraint condition of TV and OS-SART algorithm, in acceleration of iterative convergence, suppressed well noise and the pseudo-shadow in the reconstructed image.

Description of drawings

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

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

Fig. 3 is the pseudo-movie queen's of the vertical wire of embodiment of the present invention correction sinogram;

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

Fig. 5 is the X ray multi-power spectrum CT image that embodiment of the present invention utilization is 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 shown in Figure 1, the inventive method comprises the following steps: correct the pseudo-shadow of vertical wire in projection sinogram (1), utilizes the frequency domain filtering method to carry out filtering to X ray multi-power spectrum CT projection sinogram and process, and suppresses the pseudo-shadow of vertical wire in sinogram; (2) high brightness noise in projection sinogram is corrected, utilize between noise pixel and neighbor the gray difference amount to identify the high brightness noise, and the neighbor average gray value is replaced the original gray-scale value of noise, thereby remove high bright spot noise in sinogram; (3) utilization is rebuild X ray multi-power spectrum CT image based on the OS-SRAT algorithm of TV, further suppresses noise and pseudo-shadow in reconstructed image.

Concrete steps are as follows:

(1) suppress the pseudo-shadow of vertical wire in projection sinogram

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

Suppose that the every row of detector has X probe unit, the sinogram under Y scanning angle is s (k, l), k=1 wherein ..., 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, the pseudo-shadow of vertical wire produces intensity values at v=0 place, and defective pixel produces peak value at higher horizontal frequency u place, utilizes wave filter with the intensity values filtering at these frequencies place, can effectively suppress in sinogram vertically wire puppet shadow.Here wave filter is selected the fertile hereby low-pass filter of Bart, 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 experiment, can according to the pseudo-shadow characteristics of vertical wire in sinogram, adjust filter parameter.

(2) remove high brightness noise in projection sinogram

Be subjected to the impact of random noise, tend to occur the noise of high brightness in sinogram, these high brightness noises cause and occur the pseudo-shadow of some bar shapeds in reconstructed image.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.At first the present invention utilizes between noise and neighbor the gray difference amount to identify these noises, supposes that in sinogram, noise is f (i, j), and the poor quadratic sum of noise and adjacent 8 pixel grey scales 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 define that between noise and neighbor, the gray difference amount is

$c=\sqrt{D/8}$

If when measures of dispersion c is worth σ greater than certain, can judges that this pixel is noise, and give this pixel with the adjacent 8 pixel grey scale mean values of this pixel, wherein σ chooses according to actual sinogram gamma characteristic.

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

Under CT Image Iterative reconstruction framework, the 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 ^NRepresent reconstructed image, A=(a _ij) represent the reconstruction iteration matrix.Algebraic reconstruction technique (ART) is to be used for the earliest 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, i=nmod (M)+the 1st, and the index of equation,

Be the capable euclideam norm of Iterative Matrix A i.

The ART algorithm easily is subject to the salt-pepper noise impact in the reconstructed image process.Reduce relaxation factor λ and can suppress noise effect, but iterative convergence speed slows down.The SART algorithm is by the ART algorithm development, and it has kept the ART algorithm the convergence speed, has simultaneously noise inhibiting ability preferably.The 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 row element sum, J=1 ..., N represents matrix A j column element sum.

The set of supposing all projection sequence numbers is

B={1,…,M}

The set of all projection sequence numbers can be divided into T son set, and in every subset, the set expression of projection sequence number 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$

The SART reconstruction algorithm 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 constituents of inspected object have identical or approximate attenuation characteristic, just reflect the identical or close of gray scale in the CT image, so image can rarefactions.Wherein, the gradient conversion is a kind of sparse conversion commonly used, and the TV of image is through the l of its gradient image commonly used ₁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 reconstructed image f, f _i,jBe the grey scale pixel value of reconstructed image f, d _i,jIt is a discrete gradient.

Execution can be divided into inside and outside two iterative loop steps based on the Image Iterative reconstruction algorithm of TV.Outer iteration loop is carried out the OS-SART iterative reconstruction algorithm, and the internal layer iterative loop is carried out the minimization process of the total variation (TV) of reconstructed image f.When carrying out the internal layer iterative loop, can utilize gradient descent method, can be expressed as

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

Wherein, λ ' is the Gradient Descent relaxation factor,

Be a gradient direction,

Be the Gradient Descent scale parameter, m is the interior loop iterations.

Embodiment:

The present embodiment utilizes a X ray multi-power spectrum CT at a specific X ray energy range (chicken wings that is marked with contrast medium (iodine solution) of 25～40keV) scannings.

Fig. 2 is the projection sinogram of a chicken wings section of embodiment, occurs the pseudo-shadow of some vertical wire in this sinogram because detector defective (bad pixel or dead pixel) causes.Utilize the described method of step (1) to carry out filtering to the sinogram in Fig. 2 and process, suppress the pseudo-shadow of vertical wire in sinogram.

Fig. 3 is the pseudo-movie queen's of the vertical wire of embodiment correction sinogram, but still has some fragmentary high brightness noises to intersperse among in sinogram, utilizes the high brightness noise in the described method offset of sinusoidal of step (2) figure 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 the described method of step (3) respectively the sinogram in Fig. 2, sinogram and the sinogram in Fig. 4 in Fig. 3 to be carried out the CT image reconstruction, wherein carry out based on the OS-SART iterative reconstruction algorithm of TV and can divide following step: 1. input measurement data for projection b and initial pictures f=0;

2. utilize the OS-SART algorithm to upgrade reconstructed image;

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

4. repeating step 2. and 3., until satisfy rebuild convergence constraint condition till.

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

Fig. 5 is the X ray multi-power spectrum CT image that the embodiment utilization is 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, contain more ring artifact in this CT image, this is because the pseudo-shadow of some vertical wire in initial sinusoids figure causes.Fig. 5 b contains some strip artifacts for utilizing sinogram (Fig. 3) the reconstruction CT image out of revising the pseudo-movie queen of vertical wire in this CT image, this is because high brightness noise in sinogram causes.Fig. 5 c is sinogram (Fig. 4) the reconstruction CT image out after the high brightness noise is revised in utilization, through the processing to the pseudo-shadow of vertical wire and high brightness noise in initial sinusoids figure, utilization can further suppress noise and pseudo-shadow in CT image reconstruction based on the OS-SART iterative reconstruction algorithm of TV, obtains reconstructed results preferably.

Claims (2)

1. an X ray multi-power spectrum CT data for projection is processed and image rebuilding method, comprises the processing of X ray multi-power spectrum CT projection sinogram and rebuild two based on the acceleration of iterative convergence of compressed sensing going on foot greatly; Described projection sinogram is processed and is comprised that in sinogram, the pseudo-shadow of vertical wire suppresses and high brightness noise two steps of removal, and the pseudo-shadow of vertical wire suppresses to adopt the mode of frequency domain filtering, and the mode of spatial domain identification, filtering is adopted in the removal of high brightness noise; Rebuilding based on the acceleration of iterative convergence of compressed sensing is to combine with order subset while algebraic reconstruction technique (OS-SART) based on the minimized optimization constraint condition of image total variation (TV), be divided into inside and outside two-layer iterative loop step, outer iteration loop is carried out the OS-SART iterative reconstruction algorithm, and the internal layer iterative loop is carried out the minimization process of the total variation (TV) of reconstructed image f.

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

Step 1:X ray multi-power spectrum CT projection sinogram is processed

1. suppress in projection sinogram the pseudo-shadow of vertical wire: use in frequency domain filtering method offset of sinusoidal figure the pseudo-shadow of vertical wire to carry out filtering and process;

Suppose that the every row of detector has X probe unit, the sinogram under Y scanning angle is s (k, l), k=1 wherein ..., 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, the pseudo-shadow of vertical wire is in v=0 place's generation intensity values, defective pixel is at higher horizontal frequency u place's generation peak value, v and u represent respectively picture frequency conversion S (v in following formula, u) two independents variable: vertical frequency and horizontal frequency, utilize wave filter with the intensity values filtering at these frequencies place, suppress the pseudo-shadow of vertical wire in sinogram; Described wave filter is selected the fertile hereby low-pass filter of Bart, 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/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: at first utilize between noise and neighbor the gray difference amount to identify these noises, establish that in sinogram, noise is f (i, j), the poor quadratic sum of noise and adjacent 8 pixel grey scales 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 define that between noise and neighbor, the gray difference amount is

$c=\sqrt{D/8}$

When if measures of dispersion c is worth σ greater than certain, can judge that this pixel is noise, and give this pixel with the adjacent 8 pixel grey scale mean values of this pixel, 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, the CT imaging system is expressed as:

Af=b

Wherein, b=(b ¹..., b ^M) ∈ R ^MRepresent data for projection, f=(f ₁..., f _M) ∈ R ^NRepresent reconstructed image, A=(a _ij) represent the 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 row element sum, J=1 ..., N represents matrix A j column element sum, n is iterations;

The set of supposing all projection sequence numbers is

B={1,…,M}

The set of all projection sequence numbers can be divided into T son set, and in every subset, the set expression of projection sequence number 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$

The SART reconstruction algorithm 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

The total variation of reconstructed image f (TV) minimizes

$\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 reconstructed image f, f _i,jBe the grey scale pixel value of reconstructed image f, d _i,jIt is a discrete gradient;

Utilize gradient descent method, formula is as follows:

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

Wherein, λ ' is the Gradient Descent relaxation factor, Be a gradient direction,

Be the Gradient Descent scale parameter, m is the interior loop iterations.

Images