Method of reconstruction from multi-energy CT scan and device thereof

1. A method of reconstruction from a multi-energy CT scan, comprising a first data collecting step, a second data collecting step, a multi-energy correction step and a multi-energy image reconstruction step, wherein

the first data collecting step comprises:

setting a first-class object comprising a first material and a second material as an object to be scanned; performing data collection at a first-class voltage with a method for collecting multi-energy CT data, to collect n sets of first-class scan data {y i } in n projection angles, wherein i is an integer, 1 ≦i ≦n , n≧2 the first-class voltage is varied periodically during the data collection, m different sampling points are distributed between a highest tube voltage and a lowest tube voltage of the first-class voltage, performing a data collection at a first-class voltage comprises performing the data collection periodically at the voltages of the m sampling points during the data collection;

the second data collecting step comprises:

scanning a double-cylinder correction phantom in n projection angles respectively, at the respective voltages of the m sampling points with the CT scanner, to obtain m groups of data, each group comprising n sets of second-class scan data, wherein the double-cylinder correction phantom comprises two materials to be discriminated in the first-class object, one cylinder comprises a first material, and the other cylinder comprises a second material;

the multi-energy correction step comprises:

reconstructing and correcting m second-class images from the m groups of data each group comprising n sets of second-class scan data, wherein each group including n sets of second-class scan data, corresponds to one of the second-class images;

obtaining third-class images of the first material and third-class images of the second material according to the position and size of the first material and the second material in the double-cylinder correction phantom; and

calculating respective combination coefficients M H first and M H second corresponding to the voltage values of the m sampling points, according to m equations: M H first third-class image of first material +M H second third-class image of second material =the H th second-class scan image, wherein H is H th sampling point of the m tube voltage, H=0,1, . . . , m−1;

the multi-energy image reconstruction step comprises:

determining the image vectors X first and X second for which a difference between the first-class scan data {y i } and combination projection data C i first P i* X first +C i second P i* X second reaches its minimum value,

wherein C i first and C i second are combination coefficients corresponding to the case that the first-class scan data is collected at the first-class voltage in the i th projection angle, and C i first and C i second are obtained according to the respective combination coefficients M H first and M H second corresponding to the voltage values of the m sampling points; and

P i* is the i th row of a projection matrix; and

determining the image vectors X first and X second for which the difference between the first-class scan data {y i } and the combination projection data C i first P i* X first +C i second P i* X second has the minimum value, as the transpose of the first-class image vector of the first material and the transpose of the first-class image vector of the second material respectively;

displaying the reconstructed image on a display unit.

2. The method according to claim 1 , wherein the number m of different sampling points of the first-class voltage between the highest tube voltage and the lowest tube voltage is equal to or greater than 3.

3. The method according to claim 1 , wherein the determining the image vectors X first and X second for which a difference between the first-class scan data {y i } and combination projection data C i first P i* X first +C i second P i* X second reaches its minimum value comprises:

calculating a minimum value of an objective function

ψ

⁡

(

X

first

,

X

second

)

=

1

2

⁢

∑

i

⁢

(

y

i

-

(

C

first

i

⁢

P

i

*

X

first

+

C

second

i

⁢

P

i

*

X

second

)

)

2

+

β

⁡

(



X

first



TV

+



X

second



TV

)

,

wherein β is a compromise parameter.

4. The method according to claim 3 , wherein calculating the minimum value of the objective function

ψ

⁡

(

X

first

,

X

second

)

=

1

2

⁢

∑

i

⁢

(

y

i

-

(

C

first

i

⁢

P

i

*

X

first

+

C

second

i

⁢

P

i

*

X

second

)

)

2

+

β

⁡

(



X

first



TV

+



X

second



TV

)

comprises:

indicating X first by an image vector X first r ;

indicating X second by an image vector X second r ;

assuming r =0, X first 0 and X second 0 being initial image vectors;

performing an orthographic projection on X first 0 by the projection matrix P to obtain initial orthographic projection data q first 0 =P·X first 0 ; performing an orthographic projection on X second 0 by the projection matrix P to obtain initial orthographic projection data q second 0 =P·X second 0 ;

if an output condition is not satisfied, repeating iterative steps, wherein the current iteration is the (r+1) th iteration, X first r+1 indicates the image vector X first obtained in the (r+1) th iteration, and X first r indicates the image vector X first before the (r+1) th iteration, X second r+1 indicates the image vector X second obtained in the (r+1) th iteration, and X second r indicates the image vector X second before the (r+1) th iteration; the value of r is increased by 1 each time the iterative steps are performed, the iterative steps comprise: updating the image vector X first r and updating the image vector X second r , wherein

the step of updating the X first r comprises:

repeating the following steps until all pixels of X first r are updated;

selecting one pixel X first r (j) which is not updated in the (r+1) th iteration step, from the image vector X first r ;

substituting X first r+1 (j) into the equation

ψ

⁡

(

X

first

,

X

second

)

=

1

2

⁢

∑

i

⁢

(

y

i

-

(

C

first

i

⁢

P

i

*

X

first

+

C

second

i

⁢

P

i

*

X

second

)

)

2

+

β

⁡

(



X

first



TV

+



X

second



TV

)

to obtain an equation in which a partial derivative of ψ with respect to X first r+1 (j) is equal to 0:

∑

i

⁢

C

first

i

⁢

P

ij

⁡

(

y

i

-

(

C

first

i

⁢

P

i

*

X

first

r

+

C

second

i

⁢

P

i

*

X

second

r

)

)

+

(

X

first

r

⁡

(

j

)

-

X

first

r

+

1

⁡

(

j

)

)

⁢

∑

i

⁢

(

C

first

i

⁢

P

ij

)

2

-

β

⁢

∂



X

first

r



TV

∂

X

first

r

+

1

⁡

(

j

)

=

0

;

assuming

θ

1

first

=

∑

i

⁢

C

first

i

⁢

P

ij

⁡

(

y

i

-

(

C

first

i

⁢

P

i

*

X

first

r

+

C

second

i

⁢

P

i

*

X

second

r

)

)

,

⁢

θ

2

first

=

∑

i

⁢

(

C

first

i

⁢

P

ij

)

2

,

where P ij is a projection coefficient in the i th row and the j th column of the projection matrix P;

calculating θ 1 first and θ 2 first , under the conditions that P i* X first r equals the i th row of the projection data a q first r and P i* X second r equals the i th row of the projection data q second r ;

calculating X first r+1 (j) from

θ

1

first

+

θ

2

first

⁡

(

x

first

r

⁡

(

j

)

-

X

first

r

+

1

⁡

(

j

)

)

-

β

⁢

∂



X

first



TV

∂

X

first

=

0

in a dichotomy method to update the pixel X first r+1 (j); and

updating the orthographic projection data by q first r+1 =q first r +P·(X first r+1 −X first r );

the step of updating the X second r comprises:

repeating the following steps until all pixels of X second r are updated:

selecting one pixel X second r (j) which is not updated in the (r+1) th iteration step, from the image vector x second r ;

substituting X second r+1 (j) into the equation

ψ

⁡

(

X

first

,

X

second

)

=

1

2

⁢

∑

i

⁢

(

y

i

-

(

C

first

i

⁢

P

i

*

X

first

+

C

second

i

⁢

P

i

*

X

second

)

)

2

+

β

⁡

(



X

first



TV

+



X

second



TV

)

to obtain an equation in which a partial derivative of ψwith respect to X second r+1 (j) equals 0:

∑

i

⁢

C

second

i

⁢

P

ij

⁡

(

y

i

-

(

C

first

i

⁢

P

i

*

X

first

r

+

C

second

i

⁢

P

i

*

X

second

r

)

)

+

(

X

second

r

⁡

(

j

)

-

X

second

r

+

1

⁡

(

j

)

)

⁢

∑

i

⁢

(

C

second

i

⁢

P

ij

)

2

-

β

⁢

∂



X

second

r



TV

∂

X

second

r

+

1

⁡

(

j

)

=

0

assuming

θ

1

second

=

∑

i

⁢

C

second

i

⁢

P

ij

⁡

(

y

i

-

(

C

first

i

⁢

P

i

*

X

first

r

+

C

second

i

⁢

P

i

*

X

second

r

)

)

,

⁢

θ

2

second

=

∑

i

⁢

(

C

second

i

⁢

P

ij

)

2

,

where P ij is a projection coefficient in the i th row and the j th column of the projection matrix P;

calculating θ 1 second and θ 2 second under the conditions that equals the i th row of the projection data q first r and P i* X second r equals the i th row of the projection data q second r ;

calculating X X second r+1 (j) from

θ

1

second

+

θ

2

second

⁡

(

x

second

r

⁡

(

j

)

-

X

second

r

+

1

⁡

(

j

)

)

-

β

⁢

∂



X

second



TV

∂

X

second

=

0

in a dichotomy method to update the pixel X second r (j), and

updating the orthographic projection data by q second r+1 =q second r +P·(X second r+1 −X second r ).

5. The method according to claim 4 , wherein the output condition is ∥X first r+1 −X first r ∥<εand ∥X second r+1 −X second r ∥<ε, wherein εis a preset small threshold.

6. The method according to claim 5 , wherein the preset small threshold ε is at least less than 10 −3 .

7. The method according to claim 1 , wherein the step of calculating respective combination coefficients M H first and M H second corresponding to the voltage values of the m sampling points according to m equations M H first third-class image of first material +M second H third-class image of second material =the H th second-class scan image is performed in a least square method.

8. A device of a reconstruction from a multi-energy CT scan, comprising:

a first data collecting unit adapted to set a first-class object comprising a first material and a second material as object to be scanned; perform data collection at a first-class voltage with a method for collecting multi-energy CT data, to collect n sets of first-class scan data {y i } in n projection angles, wherein i is an integer, 1≦i≦n, n≧2, the first-class voltage is varied periodically during the data collection, m different sampling points are distributed between a highest tube voltage and a lowest tube voltage of the first-class voltage, and performing the data collection at the first-class voltage comprises performing the data collection periodically at the voltages of the m sampling points during the data collection; and send the n sets of the first-class scan data {y i } to a multi-energy image reconstruction unit,

a second data collecting unit adapted to scan a double-cylinder correction phantom in n projection angles respectively, at the respective voltages of the m sampling points with the CT scanner, to obtain m groups of data, each group including n sets of second-class scan data, wherein the double-cylinder correction phantom comprises two materials to be discriminated in the first-class object, one cylinder comprises a first material, and the other cylinder comprises a second material; and send the m groups each including n sets of second-class scan data to a multi-energy correction unit;

the multi-energy correction unit adapted to reconstruct and correct m second-class scan images from the m groups of n sets of second-class scan data, wherein each group, including n sets of the second-class scan data, corresponds to one of reconstructed and corrected second-class scan images; obtain a third-class image vector of the first material and a third-class image vector of the second material according to the position and size of the first material and the second material in the double-cylinder correction phantom; calculate respective combination coefficients M H first and M H second corresponding to the voltage values of the m sampling points, according to m equations: M H first third-class image vector of first material +M H second third-class image vector of second material =the H th second-class scan image; and send the respective combination coefficients M H first and M H second corresponding to the voltage values of the m sampling points to the multi-energy reconstruction unit, where H is H th sampling point of the m tube voltage, H=0,1, . . . , m− 1 ; and

a multi-energy image reconstruction unit adapted to obtain image vectors X first and X second for which a difference between the first-class scan data {y i } and combination projection data C i first P i* X first +C i second P i* X second has a minimum value, wherein C first i and C i second are combination coefficients corresponding to the case that the first-class scan data is collected at the first-class voltage in the i th projection angle, and C i first and C i second are obtained according to the respective combination coefficients M H first and M H second corresponding to the voltage values of the m sampling points, and P i* is the i th row of a projection matrix; determine the image vectors X first and X second for which a difference between the first-class scan data {y i } and the combination projection data C i first P i* X first +C i second P i* X second has the minimum value, as the transpose of the first-class image vector of the first material and the transpose of the first-class image vector of the second material respectively;

a display unit adapted to display the reconstructed image.