Limited-angle CT reconstruction method based on anisotropic total variation

Abstract

The invention discloses a limited-angle CT reconstruction method based on Anisotropic Total Variation. According to the method, through an image reconstruction model using low dose and sparse-view-angle CT images, a fast iterative reconstruction algorithm is combined with an Anisotropic Total Variation method. The problems that in an existing limited-angle CT reconstruction method are effectively solved, such as partial boundary ambiguity, slow convergence speed and unable to accurately solve. In the process of solving the model, the slope filter is introduced in the Filtered Back-Projection to preprocess the iterative equation, and the Alternating Projection Proximal is used to solve the iterative equation, and the iteration is repeated until the termination condition is satisfied; the experimental comparison with the existing reconstruction methods shows that the invention can achieve better reconstruction effect.

Claims

1. A limited-angle computed tomography (CT) reconstruction method based on Anisotropic Total Variation, comprising the following steps: (1) using a detector to collect projection data of CT images of an object in different angle directions to form a projection data set {right arrow over (b)}, wherein the projection data set {right arrow over (b)} is composed of projection data vectors corresponding to each angle direction, the dimension of the projection data vectors is n, and the value of each element is the projection data measured by the corresponding detector, and n is the number of detectors; (2) establishing CT image reconstruction equations in a low dose mode and a sparse view mode, respectively, by a processing circuit; the expression of CT image reconstruction equation in low dose mode being: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }\left(\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{b}.Math.}^2\right),s.t.\stackrel{.fwdarw.}{x}\ge 0& \left(1.1\right)\end{array}$ the expression of CT image reconstruction equation in sparse view mode being: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }{.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}},s.t.A\stackrel{.fwdarw.}{x}=\stackrel{.fwdarw.}{b},\stackrel{.fwdarw.}{x}\ge 0& \left(2.1\right)\end{array}$ wherein, {right arrow over (x)} is the CT image data set and the dimension is k, k is the total number of pixels of the CT image, the value of each element {right arrow over (x)} is the X-ray absorption coefficient of the corresponding pixel in the CT image to be reconstructed, A is a system matrix, and β is a given weight ∥ ∥.sub.ATV indicates Anisotropic Total Variation; (3) pre-processing the CT image reconstruction equations in the low dose mode and the sparse view mode by the processing circuit to obtain the corresponding objective function, the pre-processing process comprising introducing a variable and a non-singular matrix P, and using Lagrange duality to convert the CT image reconstruction equations into a saddle point problem; (4) according to the low dose image reconstruction mode or the sparse view image mode, selecting and optimizing the objective function under the corresponding mode to obtain the CT image, and outputting the CT image to an output display; wherein, in the step (4), an Alternating Projection Proximal is used to optimize and solve the objective function to reconstruct the CT image, the parameters involved in the alternating projection approach algorithm include a total number of iterations N.sub.iter, a number of TV denoising iterations N.sub.ATV, and a iteration termination threshold δ and ∈, a TV term weighting coefficient β, a parameter α, a parameter τ, a parameter σ, an ATV weighting factors γ.sub.h and γ.sub.v; wherein, the total number of iterations N.sub.iter ranges from 1˜50000, the number of TV denoising iterations N.sub.ATV ranges from 1˜10000, the iteration termination threshold δ and ∈ ranges from 10.sup.−10˜1, the TV term weighting coefficient β ranges from 10.sup.−6˜1, the parameter α ranges from 0.01˜0.25, the parameter τ ranges from 10.sup.−3˜10.sup.3, the parameter σ ranges from 10.sup.−3˜10.sup.3, the ATV weighting factors γ.sub.h and γ.sub.v ranges from 10.sup.−6˜1.

2. The limited-angle CT reconstruction method according to claim 1, wherein, the pre-processing process of the CT image reconstruction equation in the low dose mode in the step (3) is as follows: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }\left(\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.\stackrel{.fwdarw.}{y}-\stackrel{.fwdarw.}{b}.Math.}^2\right),s.t.A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{y}=0,\stackrel{.fwdarw.}{x}\ge 0& \left(1.2\right)\end{array}$ first, introducing the variable {right arrow over (y)}, the CT image reconstruction equation in the low dose mode is rewritten to the following form: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }\left(\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.\stackrel{.fwdarw.}{y}-\stackrel{.fwdarw.}{b}.Math.}^2\right),s.t.A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{y}=0,\stackrel{.fwdarw.}{x}\ge 0& \left(1.2\right)\end{array}$ wherein, {right arrow over (y)} indicates the forward projection data calculated from the CT image data set {right arrow over (x)}; then, introducing the non-singular matrix P, equation (1.2) being further rewritten to the following form: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }\left(\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.\stackrel{.fwdarw.}{y}-\stackrel{.fwdarw.}{b}.Math.}^2\right),s.t.P^{1/2}\left(A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{y}\right)=0,\stackrel{.fwdarw.}{x}\ge 0& \left(1.3\right)\end{array}$ the Lagrange equation corresponding to equation (1.3) being as follows:
L({right arrow over (x)},{right arrow over (y)},{right arrow over (μ)})=β∥{right arrow over (x)}∥.sub.ATV+½∥{right arrow over (y)}−{right arrow over (b)}∥.sup.2+{right arrow over (μ)}.sup.TP.sup.1/2(A{right arrow over (x)}−{right arrow over (y)}) wherein, {right arrow over (μ)} is the Lagrange multiplier vector, T indicates transposition; $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{\mu }}{\max }\underset{\stackrel{.fwdarw.}{x}\ge 0,\stackrel{.fwdarw.}{y}}{\min }L\left(\stackrel{.fwdarw.}{x},\stackrel{.fwdarw.}{y},\stackrel{.fwdarw.}{\mu }\right)=\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.\stackrel{.fwdarw.}{y}-\stackrel{.fwdarw.}{b}.Math.}^2+{\stackrel{.fwdarw.}{\mu }}^T{P}^{1/2}\left(A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{y}\right)& \left(1.4\right)\end{array}$ the equation (1.3) is converted to the saddle point problem, the specific expression is as follows: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{\mu }}{\max }\underset{\stackrel{.fwdarw.}{x}\ge 0,\stackrel{.fwdarw.}{y}}{\min }L\left(\stackrel{.fwdarw.}{x},\stackrel{.fwdarw.}{y},\stackrel{.fwdarw.}{\mu }\right)=\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.\stackrel{.fwdarw.}{y}-\stackrel{.fwdarw.}{b}.Math.}^2+{\stackrel{.fwdarw.}{\mu }}^T{P}^{1/2}\left(A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{y}\right)& \left(1.4\right)\end{array}$ finally, minimizing the equation (1.4) and removing the variable {right arrow over (y)}, and the corresponding objective function being obtained as follows: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{\mu }}{\max }\underset{\stackrel{.fwdarw.}{x}\ge 0,\stackrel{.fwdarw.}{y}}{\min }L\left(\stackrel{.fwdarw.}{x},\stackrel{.fwdarw.}{\mu }\right)=\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}-\frac{1}{2}{.Math.{\left(P^{1/2}\right)}^T\stackrel{.fwdarw.}{\mu }.Math.}^2-{\stackrel{.fwdarw.}{\mu }}^T{P}^{1/2}\stackrel{.fwdarw.}{b}+{\left(P^{1/2}\right)}^T\stackrel{.fwdarw.}{\mu }.Math.A\stackrel{.fwdarw.}{x}& \left(1.5\right)\end{array}$ wherein, ● indicates inner product operation.

3. The limited-angle CT reconstruction method according to claim 1, wherein, the pre-processing process of the CT image reconstruction equation in sparse view mode in the step (3) is as follows: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }{.Math.\overline{x}.Math.}_{ATV}s.t.P^{1/2}\left(A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{b}\right)=0\stackrel{.fwdarw.}{x}\ge 0& \left(2.2\right)\end{array}$ first, introducing the variable {right arrow over (y)}, the CT image reconstruction equation in sparse view mode being rewritten to the following form: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }{.Math.\overline{x}.Math.}_{ATV}s.t.P^{1/2}\left(A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{b}\right)=0\stackrel{.fwdarw.}{x}\ge 0& \left(2.2\right)\end{array}$ the Lagrange equation corresponding to equation (2.2) being as follows:
L({right arrow over (x)},{right arrow over (μ)})=β∥{right arrow over (x)}∥.sub.ATV+{right arrow over (μ)}.sup.TP.sup.1/2(A{right arrow over (x)}−{right arrow over (b)}) the equation (1.3) being converted to the saddle point problem, and the corresponding objective function being obtained as follows: $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{\mu }}{\max }\underset{\stackrel{.fwdarw.}{x}\ge 0}{\min }L\left(\stackrel{.fwdarw.}{x},\stackrel{.fwdarw.}{\mu }\right)=\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+{\stackrel{.fwdarw.}{\mu }}^T{P}^{1/2}\left(A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{b}\right)& \left(2.3\right)\end{array}$ wherein, {right arrow over (μ)} is the Lagrange multiplier vector, T indicates transposition, β is the given weight.

4. The limited-angle CT reconstruction method according to claim 2, wherein, in the low dose mode, the calculation expression of the non-singular matrix P is as follows: $P=F^{-1}R\left(\omega \right)FR\left(\omega \right)=\frac{.Math.\omega .Math.}{2m\tau +.Math.\omega .Math.}$ wherein, F and F.sup.−1 are one-dimensional Fourier transform operator and inverse Fourier transform operator respectively, ω is the frequency domain variable after the Fourier transform of the projection data vector collected under the corresponding angle direction, m is the number of angle directions, and τ is a given parameter.

5. The limited-angle CT reconstruction method according to claim 3, wherein, in the sparse view mode, the calculation expression of the non-singular matrix P is as follows: $P=F^{-1}R\left(\omega \right)FR\left(\omega \right)=\frac{.Math.\omega .Math.}{2m\tau }$ wherein, F and F.sup.−1 are one-dimensional Fourier transform operator and inverse Fourier transform operator respectively, ω is the frequency domain variable after the Fourier transform of the projection data vector collected under the corresponding angle direction, m is the number of angle directions, and τ is a given parameter.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a schematic flowchart of a CT image reconstruction algorithm of the present invention.

(2) FIG. 2 shows the true value image of the Shepp & Logan phantom model, showing the gray range [1.00, 1.05].

(3) FIG. 3 shows the gray value image corresponding to the true value image of the Shepp & Logan phantom model.

(4) FIG. 4 shows the sonogram image (that is, a image composed of projection data at different angles, the angle range is 180°, and the angles is 180°) corresponding to the true value image of the Shepp & Logan phantom model.

(5) FIG. 5 shows the change of the root mean square error with the number of iterations of the reconstructed images obtained by the two methods (the method of the present invention and the Anisotropic Total Variation method) in the reconstruction process compared to the truth map, the angle range of the projection data is 120°, the angle is 120°.

(6) FIG. 6(a) shows the CT image reconstructed by the present invention by using the Shepp & Logan phantom model, the projection data angle range is 120% and the number of iterations is 3000 times.

(7) FIG. 6(b) shows the CT image reconstructed by the Anisotropic Total Variation method by using the Shepp & Logan phantom model, the projection data angle range is 120% and the number of iterations is 3000 times.

(8) FIG. 7 is a schematic diagram comparing the horizontal midline of the reconstructed image and the true value of the reconstructed image obtained by the two methods (the method of the present invention and the Anisotropic Total Variation method), and the number of iterations is 3000 times.

DETAILED DESCRIPTION OF THE INVENTION

(9) In order to describe the present invention more specifically, the technical solution of the present invention will be described in detail below with reference to the drawings and specific embodiments.

(10) As shown in FIG. 1, the limited-angle CT reconstruction method based on Anisotropic Total Variation of the present invention comprises the following steps:

(11) (1) collecting CT images in different angle direction by using a detector to form a projection data set {right arrow over (b)}, the projection data set {right arrow over (b)} is composed of projection data vectors corresponding to each angle directions, the dimension of the projection data vectors is n, and the value of each element is the projection data measured by the corresponding detector, and n is the number of detectors.

(12) (2) Establishing CT image reconstruction equations in a low dose mode and a sparse view mode respectively:

(13) The Low dose CT:

(14) $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }\left(\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{b}.Math.}^2\right),s.t.\stackrel{.fwdarw.}{x}\ge 0& \left(1\right)\end{array}$

(15) The Sparse view CT:

(16) $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }{.Math.\stackrel{.fwdarw.}{x}.Math.}_{ATV},s.t.A\stackrel{.fwdarw.}{x}=\stackrel{.fwdarw.}{b},\stackrel{.fwdarw.}{x}\ge 0& \left(2\right)\end{array}$

(17) Wherein, {right arrow over (x)} is the CT image data set and the dimension is k, k is the total number of pixels of the CT image, the value of each element {right arrow over (x)} is the X-ray absorption coefficient of the corresponding pixel in the CT image to be reconstructed, A is the system matrix, A is suitable for various forms such as parallel beam projection, fan beam projection and cone beam projection, and β is a given weight, ∥ ∥.sub.ATV indicates Anisotropic Total Variation.

(18) In the case of a low dose CT, assuming that nm>k and {right arrow over (b)} contains noise. In the case of a sparse view CT, assuming that nm<k and {right arrow over (b)} does contain noise, m is the number of angle directions.

(19) (3) Pre-processing the equation (1) and (2):

(20) For the low dose CT, introducing the variable {right arrow over (y)}, the equation (1) is rewritten to the following form:

(21) $\underset{\stackrel{.fwdarw.}{x}}{\min }\left(\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.\stackrel{.fwdarw.}{y}-\stackrel{.fwdarw.}{b}.Math.}^2\right),s.t.A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{y}=0,\stackrel{.fwdarw.}{x}\ge 0$

(22) Wherein, {right arrow over (y)} indicates the forward projection data calculated from the CT image data set {right arrow over (x)}.

(23) Introducing the non-singular matrix P.sup.1/2, the above equation is rewritten to the following form:

(24) $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }\left(\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.\stackrel{.fwdarw.}{y}-\stackrel{.fwdarw.}{b}.Math.}^2\right),s.t.P^{1/2}\left(A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{y}\right)=0,\stackrel{.fwdarw.}{x}\ge 0& \left(3\right)\end{array}$

(25) The Lagrange equation corresponding to equation (3) is as follows:
L({right arrow over (x)},{right arrow over (y)},{right arrow over (μ)})=β∥{right arrow over (x)}∥.sub.ATV+½∥{right arrow over (y)}−{right arrow over (b)}∥.sup.2+{right arrow over (μ)}.sup.TP.sup.1/2(A{right arrow over (x)}−{right arrow over (y)})

(26) Wherein, {right arrow over (μ)} is the Lagrange multiplier vector, equation (3) is converted to the saddle point problem:

(27) $\underset{\stackrel{.fwdarw.}{\mu }}{\max }\underset{\stackrel{.fwdarw.}{x}\ge 0,\stackrel{.fwdarw.}{y}}{\min }L\left(\stackrel{.fwdarw.}{x},\stackrel{.fwdarw.}{y},\stackrel{.fwdarw.}{\mu }\right)=\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.\stackrel{.fwdarw.}{y}-\stackrel{.fwdarw.}{b}.Math.}^2+{\stackrel{.fwdarw.}{\mu }}^T{P}^{1/2}\left(A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{y}\right).$

(28) The variable {right arrow over (y)} can be removed by minimizing the above equation, and the corresponding objective function is obtained as follows:

(29) $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{\mu }}{\max }\underset{\stackrel{.fwdarw.}{x}\ge 0}{\min }L\left(\stackrel{.fwdarw.}{x},\stackrel{.fwdarw.}{\mu }\right)=\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}-\frac{1}{2}{.Math.{\left(P^{1/2}\right)}^T\stackrel{.fwdarw.}{\mu }.Math.}^2-{\stackrel{.fwdarw.}{\mu }}^T{P}^{1/2}\stackrel{.fwdarw.}{b}+{\left(P^{1/2}\right)}^T\stackrel{.fwdarw.}{\mu }.Math.A\stackrel{.fwdarw.}{x}P=F^{-1}R\left(\omega \right)FR\left(\omega \right)=\frac{.Math.\omega .Math.}{2m\tau +.Math.\omega .Math.}& \left(4\right)\end{array}$

(30) Wherein, .circle-solid. indicates inner product operation, F and F.sup.−1 are one-dimensional Fourier transform operator and inverse Fourier transform operator respectively, ω is the frequency domain variable after the Fourier transform of the projection data vector collected under the corresponding angle direction.

(31) For sparse view CT, introducing the non-singular matrix P.sup.1/2, the equation (2) is rewritten to the following form:

(32) $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{x}}{\min }{.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}},s.t.P^{1/2}\left(A\stackrel{.fwdarw.}{x}-\stackrel{.fwdarw.}{\gamma }\right)=0,\stackrel{.fwdarw.}{x}\ge 0& \left(5\right)\end{array}$

(33) The Lagrange equation corresponding to equation (5) is as follows:
L({right arrow over (x)},{right arrow over (μ)})=β∥{right arrow over (x)}∥.sub.ATV+{right arrow over (μ)}.sup.TP.sup.1/2(A{right arrow over (x)}−{right arrow over (b)})

(34) The equation (5) is converted to the following form:

(35) $\begin{array}{cc}\underset{\stackrel{.fwdarw.}{\mu }}{\max }\underset{\stackrel{.fwdarw.}{x}\ge 0}{\min }L\left(\stackrel{.fwdarw.}{x},\stackrel{.fwdarw.}{\mu }\right)=\beta {.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+{\stackrel{.fwdarw.}{\mu }}^T{P}^{1/2}\left(A\stackrel{.fwdarw.}{x}-b\right)P=F^{-1}R\left(\omega \right)FR\left(\omega \right)=\frac{.Math.\omega .Math.}{2m\tau }& \left(6\right)\end{array}$

(36) (4) Using an alternating projection approach algorithm to solve the equation (4) and (6):

(37) 4-1 Initializing {right arrow over (x)}=0, {right arrow over (μ)}=0, and setting the value of each parameter: total number of iterations N.sub.iter, the number of TV denoising iterations N.sub.ATV, the iteration termination threshold δ and ∈, the TV term weighting coefficient β, the parameter α, the parameter τ, the parameter σ, and the ATV weighting factors γ.sub.h and γ.sub.v; wherein, the total number of iterations N.sub.iter ranges from 1˜50000, the number of TV denoising iterations N.sub.ATV ranges from 1˜10000, the iteration termination threshold δ and ∈ ranges from 10.sup.−10˜1, the TV term weighting coefficient β ranges from 10.sup.−6˜1, the parameter α ranges from 0.01˜0.25, the parameter τ ranges from 10.sup.−3˜10.sup.3, the parameter σ ranges from 10.sup.−3˜10.sup.3, the ATV weighting factors γ.sub.h and γ.sub.v ranges from 10.sup.−6˜1.

(38) 4-2 The iteration count k=1 is setted initially.

(39) 4-3 Updating {right arrow over (μ)}, when k is 1:
{right arrow over (μ)}.sup.(k+1)=−σP{right arrow over (b)}

(40) When k is not 1:

(41) The Low dose CT: {right arrow over (μ)}.sup.(k+1)=2{right arrow over (μ)}.sup.(k)−{right arrow over (μ)}.sup.(k−1)−σP({right arrow over (μ)}.sup.(k)−{right arrow over (μ)}.sup.(k−1))

(42) The Sparse view CT: {right arrow over (μ)}.sup.(k+1)=2{right arrow over (μ)}.sup.(k)−{right arrow over (μ)}.sup.(k−1)

(43) 4-4 Updating {right arrow over (x)}:

(44) $\begin{array}{cc}{\stackrel{.fwdarw.}{x}}^{\left(k+1\right)}=\arg \underset{\stackrel{.fwdarw.}{x}\ge 0}{\min }\left(\mathrm{\tau \beta }{.Math.\stackrel{.fwdarw.}{x}.Math.}_{\mathrm{ATV}}+\frac{1}{2}{.Math.\stackrel{.fwdarw.}{x}-\left({\stackrel{.fwdarw.}{x}}^{\left(k\right)}-\tau A^T{\stackrel{.fwdarw.}{\mu }}^{\left(k+1\right)}\right).Math.}^2\right)& \left(7\right)\end{array}$

(45) The equation (7) can be solved using the TV denoising algorithm of Chamball, {right arrow over (x)}.sup.(k)−τA.sup.T{right arrow over (μ)}.sup.(k+1) is constant in this step, {right arrow over (z)}={right arrow over (x)}.sup.(k)−τA.sup.T{right arrow over (μ)}.sup.(k+1) is introduced for convenience, then the solving equation of equation (7) is:

(46) 0$\begin{array}{cc}{p}_{s,t}^{\left(k+1\right)}=\frac{p_{s,t}^{\left(k\right)}+{\alpha \left(\nabla \left(div{\stackrel{.fwdarw.}{p}}^{\left(k\right)}-\stackrel{.fwdarw.}{z}/\left(\tau \beta \right)\right)\right)}_{s,t}}{1+\alpha .Math.{\left(\nabla \left(\mathrm{div}{\stackrel{.fwdarw.}{p}}^{\left(k\right)}-\stackrel{.fwdarw.}{z}/\left(\tau \beta \right)\right)\right)}_{s,t}.Math.}& \left(8\right)\\ {\stackrel{.fwdarw.}{x}}^{\left(k+1\right)}=\stackrel{.fwdarw.}{z}-\tau \beta div{\stackrel{.fwdarw.}{p}}^{\left(k+1\right)}& \left(9\right)\end{array}$

(47) Wherein, ∇ is a gradient operator, assuming that the size of the reconstructed image is N×N pixels. In this embodiment, it is defined as follows:

(48) ${\left(\nabla \stackrel{.fwdarw.}{x}\right)}_{s,t}=\left({\left(\nabla \stackrel{.fwdarw.}{x}\right)}_{s,t}^1,{\left(\nabla \stackrel{.fwdarw.}{x}\right)}_{s,t}^2\right)$${\left(\nabla \stackrel{.fwdarw.}{x}\right)}_{s,t}^1=\{\begin{array}{cc}\left(x_{s+1,t}-x_{s,t}\right){\gamma }_h,& \mathrm{if}s<N\\ 0,& \mathrm{if}s=N\end{array}{\left(\nabla \stackrel{.fwdarw.}{x}\right)}_{s,t}^2=\{\begin{array}{cc}\left(x_{s+1,t}-x_{s,t}\right){\gamma }_v,& \mathrm{if}t<N\\ 0,& \mathrm{if}t=N\end{array}$

(49) Then the Anisotropic Total Variation sub-item ∥{right arrow over (x)}∥.sub.ATV is calculated as follows:

(50) ${.Math.\stackrel{.fwdarw.}{x}.Math.}_{ATV}=\underset{1\le s,t\le N}{.Math.}.Math.{\left(\nabla \stackrel{.fwdarw.}{x}\right)}_{s,t}.Math.$

(51) For each p.sub.s,t=(p.sub.s,t.sup.1,p.sub.s,t.sup.2), the div operator is calculated as follows:

(52) ${\left(div\stackrel{.fwdarw.}{p}\right)}_{s,t}=\{\begin{array}{cc}{p}_{s,t}^1-p_{s-1,t}^1,& \mathrm{if}1<s<N\\ p_{s,t}^1,& \mathrm{if}s=1\\ -p_{s-1,t}^1,& \mathrm{if}s=N\end{array}+\{\begin{array}{cc}{p}_{s,t}^1-p_{s,t-1}^1,& \mathrm{if}1<t<N\\ p_{s,t}^1,& \mathrm{if}t=1\\ -p_{s,t-1}^1,& \mathrm{if}t=N\end{array}$

(53) Calculating ∈.sub.c=(∥{right arrow over (x)}.sup.(k+1)∥.sub.ATV−∥{right arrow over (x)}.sup.(k)∥.sub.ATV)/∥{right arrow over (x)}.sup.(k)∥.sub.ATV, if ∈.sub.c>∈ and the number of repetitions does not reach N.sub.ATV, the equation (8) and (9) is repeatedly calculated.

(54) 4-5 Updating {right arrow over (x)} again:

(55) The Low dose CT: {right arrow over (μ)}.sup.(k+1)={right arrow over (μ)}.sup.(k)σP(A{right arrow over (x)}.sup.(k+1)−{right arrow over (b)}−{right arrow over (μ)}.sup.(k))

(56) The Sparse view CT: {right arrow over (μ)}.sup.(k+1)={right arrow over (μ)}.sup.(k)+σP(A{right arrow over (x)}.sup.(k+1)−{right arrow over (b)})

(57) 4-6 Determine whether the sum of squares of all pixel values of the difference image of the reconstructed image {right arrow over (x)}.sup.(k+1) and he reconstructed image {right arrow over (x)}.sup.(k) is less than the iteration termination threshold δ or whether k is greater than N.sub.iter; if yes, execute the step (4-7); if not, let k=k+1, execute the step (4-3)˜(4-5).

(58) 4-7 Terminating the iteration and getting the final reconstructed image

(59) The practicability and reliability of the implementation method by reconstructing the sinogram of the Shepp & Logan phantom mode are verified. The true value of the Shepp & Logan phantom mode is shown if FIG. 2 (due to the gray value of most of its structures is very close, in order to clearly display its structure, select the gray range [1.00, 1.05] for display), the corresponding gray value graph is shown in FIG. 3, and the sonogram (the angle range is 180°,) is shown in FIG. 4. The proposed method and the Anisotropic Total Variation method (ART+ATV) are used for reconstruction. FIG. 5 shows the change of the root mean square error with the number of iterations of the reconstructed images obtained by the two methods in the reconstruction process compared to the truth map. The FIG. 5 shows both methods gradually converge with the increase of the number of iterations, the root mean square error of the reconstruction result of the method of the present invention is smaller, and the convergence speed is very fast. The reconstructed image at the maximum number of iterations (that is, 3000 times) is selected and displayed. FIG. 6 (a) is the CT image reconstructed by the Anisotropic Total Variation method, and FIG. 6 (b) is the reconstructed image by the method of the present invention. It can be seen that the image artifacts reconstructed by the method of the present invention are much smaller, and the detailed information of the image can be better recovered, and the reconstruction effect is significantly better than the anisotropic full variation method. In order to more intuitively prove the superiority of the method of the present invention, FIG. 7 shows the comparison between the horizontal midline profile of the reconstructed image shown in FIGS. 6 (a) and 6 (b) and the true value profile. The section line of the method of the present invention almost completely coincides with the true value, that is, it can obtain a reconstruction result closer to the true value.

(60) The above description of the embodiments is to facilitate those of ordinary skill in the art to understand and apply the present invention. Those skilled in the art can obviously make various modifications to the above-mentioned embodiments, and apply the general principles described here to other embodiments without creative work. Therefore, the present invention is not limited to the above-mentioned embodiments. According to the disclosure of the present invention, those skilled in the art should make improvements and modifications to the present invention within the protection scope of the present invention.