CN103390285A - Cone beam computed tomography (CT) incomplete angle rebuilding method based on edge guide - Google Patents

Abstract

The invention relates to a cone beam computed tomography (CT) incomplete angle rebuilding method based on edge guide. The method specifically includes the following steps: 1 estimating an initial image and utilizing scanned projection data to estimate an initial rebuilt image; 2 extracting the image edge; 3 designing a weighting factor; 4 updating and optimizing a model; 5 rebuilding a cone beam CT incomplete angle based on sparse optimization; 6 judging whether rebuilding quality meets the requirement, executing step 7 on yes judgment and executing the step 2 on no judgment; 7 finishing. The cone beam CT incomplete angle rebuilding method based on the edge guide is high in efficiency and good in rebuilt image quality.

Description

The incomplete angle method for reconstructing of Cone-Beam CT based on margin guide

(1), technical field: the present invention relates to the incomplete angle method for reconstructing of a kind of Cone-Beam CT, particularly relate to the incomplete angle method for reconstructing of a kind of Cone-Beam CT based on margin guide.

(2), background technology: Computed tomography (Computed Tomography, CT) has been widely used in the fields such as medical science, industry as a kind of modern imaging technique.Yet, in a lot of practical applications, because the geometric position that is subjected to data acquisition time or imaging system scanning retrains, can only or at less projection angle, obtain data at incomplete angular range, these all belong to incomplete angle problem (Incomplete Data Problem).The quality that is lifted at CT image reconstruction under incomplete angle scanning has important theoretical research and engineering practice meaning, and the method that how to design CT image reconstruction under high-precision incomplete angle is also focus and the difficulties of research.

For the incomplete angle problem of Cone-Beam CT, its data acquisition does not meet Exact Reconstruction data completeness condition, resolves the class reconstruction algorithm and can't obtain the reconstructed image of better quality.The Class of Iterative reconstruction algorithm does not have strict requirement to the data completeness, can obtain the more excellent reconstruction quality of relative parsing class algorithm.Classic algorithm is the algebraically iterative technique, and algebraic reconstruction algorithm has certain anti-shortage of data, and on same quantity of data, result is better than analytical algorithm usually, but takies the calculating storage resources, needs stronger support, rebuilds slower in reality.Based on the CT reconstruction algorithm of compressive sensing theory by the priori of excavating object to be rebuild and the sparse characteristic of portraying object, reconstruction quality that can be more excellent than classical iterative algorithm under the situation of compression sampling.

Although the reconstruction algorithm based on the CS theoretical model can, than obtaining reconstructed results preferably under sparse sampling, but not excavated the image information that can more be conducive to improve reconstructed results to a deeper level.The marginal element of image is the important elements of image, in image repair, the multiple image applications field such as cut apart important application arranged.In image reconstruction, effectively utilize marginal information can improve largely reconstruction quality.

For incomplete angle CT image reconstruction problem, classical disposal route is algebraically iteration and statistics iterative reconstruction algorithm, algebraic reconstruction algorithm (Algebraic Reconstruction Technique for example, ART) and maximum (the expectation maximization of expectation, EM) algorithm, iterative algorithm has reconstruction quality preferably than analytical algorithm usually when the processing missing data is less, as the FDK algorithm by the people such as Feldkamp proposition in 1984, reconstruction quality is preferably arranged.In recent years, under guidance based on compressive sensing theory, having grown up, some can be applicable to the reconstruction algorithm of sparse sampling, ASD-POCS (the adapt-steepest-descent projection onto convex sets) algorithm that comparatively typically has the people such as Sidky in 2008 and Pan to propose.This algorithm, with TV Norm minimum design optimization model, adopts POCS(ART) descend and replace the strategy of carrying out with the TV steepest, can reconstruct the image of better quality under sparse sampling.

After the ASD-POCS algorithm was suggested, a lot of scholars had proposed again multiple optimized algorithm on model solution, such as the Split-BregmanTV algorithm by the people such as Vandeghinste proposition in 2011, and the ADTVM algorithms by the people such as Zhang proposition in 2012 etc.Other algorithms are introduced other priori etc., PICCS(prior image constrained compressed sensing for example, PICCS), Wang Lin units in 2010 wait the RRD(reconstruction-reference difference of people's proposition, RRD) scheduling algorithm is all demand object to be rebuild sparse expressions under certain priori, utilize sparse expression constitution optimization model and design corresponding derivation algorithm, thereby reaching the purpose of image reconstruction.The comparatively strong instrument of the limited angle of frequency domain aspect reconstruction in addition (main MRI imaging) is RecPF(reconstruction from partial Fourier space smpling) algorithm, this algorithm does not relate to storage and the calculating of extensive matrix, utilizes difference operator D _iGood character and realize to rebuild efficient of FFT technology and accurately, be the example of CS theory in the MRI successful Application.RecPF also has stronger application prospect in CT image reconstruction field, and the people such as Zhang Hanming successfully realize the CT image reconstruction of degree of precision based on over-sampling and sparse interpolation technique.

2010, the people such as Yin proposed the method for utilizing iteration support to survey (iterative support detection, ISD) and obtained more excellent reconstructed results in sparse signal recover.In the method, constantly the iterative process M signal is done the processing that support set is surveyed, the objective function that props up the set pair Optimized model according to the part that detects upgrades, take this alternative manner can be on the basis of extremely owing to sample effective restoring signal.The pure method of surveying based on support can not directly use in the image reconstruction of Cone-Beam CT, utilizes the support set of image sparse under representing to integrate can bring as the CT image reconstruction lifting of reconstruction quality.

(3), summary of the invention:

The technical problem to be solved in the present invention is: overcome the defect of prior art, provide the Cone-Beam CT based on margin guide that a kind of efficiency is high, reconstructed image quality is good incomplete angle method for reconstructing.

Technical scheme of the present invention:

The incomplete angle method for reconstructing of a kind of Cone-Beam CT based on margin guide contains and has the following steps:

Step 1: estimate initial pictures: utilize the data for projection that scans to estimate initial reconstructed image;

Step 2: Edge extraction;

Step 3: design weighting factor;

Step 4: upgrade Optimized model;

Step 5: based on the incomplete angle of sparse optimization Cone-Beam CT, rebuild;

Step 6: judge that reconstruction quality reaches requirement? in this way, perform step 7; , if not being, perform step 2;

Step 7: finish.

The concrete method of estimation of step 1 is as follows:

Step 1.1: the image reconstruction problem is portrayed as following sparse model:

min||x|| _TV

s.t.Ax=b

Wherein, || x|| _TVFor total variation (TV) norm of object x to be rebuild, A is system matrix, and vectorial b is the data for projection that scans;

Step 1.2: with alternating direction method, the sparse model in step 1.1 is solved, iterations is N, obtains solution formula as follows:

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T({z_j}^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ {z_j}^{k+1}=\max \left\{\right|D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j})\\ {u_j}^{k+1}={u_j}^k+{\mathrm{\ρ}}_j(D_j{x}^{k+1}-{z_j}^{k+1})\end{array}\right.$

Wherein, () ⁺For the Moore-Penrose pseudoinverse of matrix, A ^TFor the transposition of system matrix, u _jFor the renewal factor of j direction, u _j ^kBe the renewal factor of the k time j direction after iteration, u _j ^k+1Be the renewal factor of the k+1 time j direction after iteration, ρ _jFor the L2 norm penalty factor of j direction, λ is L1 norm penalty factor, D _jFor the gradient operator matrix on the j direction, D _j ^TFor the gradient operator transpose of a matrix on the j direction, z _jFor j direction gradient image, z _j ^kBe the k time j direction gradient image after iteration, z _j ^k+1Be the k+1 time j direction gradient image after iteration, max{}, sgn{} are respectively and ask maximal function and ask sign function, vectorial b is the data for projection that scans, x ^k+1Be the k+1 time reconstructed image after iteration, through the reconstructed image after N iteration, be denoted as x ^(N), x ^(N)For the last initial reconstructed image of estimating.

Iteration in step 1.2 time N is 0.5～0.2 times of total convergence wheel number.

The concrete grammar of step 2 is: establishing the intermediate reconstructed images that in step 1.2, certain iteration produces is x1, utilizes noise reduction factor F to carry out convolution algorithm: x1 to middle reconstructed image x1 _F=x1***F, utilize classical edge operator E to ask for the preliminary edge x1 of intermediate reconstructed images x1 _E-mid: x1 _E-mid=E[x1], the edge x1 of intermediate reconstructed images x1 _eFor: x1 _e=x1 _E-mid∩ x1 _LMI, wherein, x1 _LMILocal mutual information image for intermediate reconstructed images x1;

The concrete grammar of step 3 is: according to the edge x1 of the intermediate reconstructed images x1 in step 2 _eλ determines the weighting factor w (i1, j1) that intermediate reconstructed images x1 locates at coordinate (i1, j1) with L1 norm penalty factor, this weighting factor w (i1, j1)=λ x1 _e(i1, j1), wherein, x1 _e(i1, j1) is intermediate reconstructed images x1 at the edge that coordinate (i1, j1) is located;

The concrete grammar of step 4 is: the weighting diagonal matrix M that calculates the j direction according to the weighting factor w (i1, j1) in step 3 _j, M _j=diag (w (1,1), w (1,2) ... w (i1, j1) ..., w (N _x, N _y)), Optimized model is updated to following expression so:

$\min {\left|\right|z_1\left|\right|}_1+{\left|\right|z_2\left|\right|}_1+{\left|\right|z_3\left|\right|}_1$

$s.t.\left\{\begin{array}{c}\mathrm{Ax}=b\\ M_j{D}_{1}x=z_1\\ M_j{D}_{2}x=z_2\\ M_j{D}_{3}x=z_3\end{array}\right.$

Wherein, N _xFor the scale of intermediate reconstructed images x1 on the x coordinate, N _yFor the scale of intermediate reconstructed images x1 on the y coordinate, D ₁Be the gradient operator matrix on 1 direction, D ₂Be the gradient operator matrix on 2 directions, D ₃It is the gradient operator matrix on 3 directions;

The concrete grammar of step 5 is: the Optimized model of the belt restraining after adopting the augmentation method of Lagrange multipliers with the renewal in step 4 transfers unconfined Optimized model to, and concrete formula is as follows:

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{j}{\mathrm{\Σ}}(u_j^T(M_j{D}_{j}x-z_j)+\frac{{\mathrm{\ρ}}_j}{2}{{\left|\right|M_j{D}_{j}x-z_j\left|\right|}_2^2+\mathrm{\λ}\left|\right|z_j\left|\right|}_1)$

Wherein,

Transposition for the renewal factor of j direction;

Above-mentioned unconfined Optimized model is adopted variables separation, utilize alternating direction method to ask minimum, iterative formula is as follows:

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M}_j{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M_j}^T(z_j^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ z_j^{k+1}=\max \left\{\right|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}[M_j{D}_j{x}^{k+1}+\frac{u_j^k}{{\mathrm{\ρ}}_j}]\\ u_j^{k+1}=u_j^k+{\mathrm{\ρ}}_j(M_j{D}_j{x}^{k+1}-z_j^{k+1})\end{array}\right.$

x ^k+1Be the epicycle reconstructed image;

Reconstruction quality in step 6 reaches requirement and refers to: the epicycle reconstructed image is compared without marked change with last round of reconstructed image; During reconstruction in carry out step 5 first, last round of reconstructed image refers to the initial reconstructed image in step 1.

Noise reduction factor F is gaussian kernel function or median filter, classical edge operator E is any in canny boundary operator, soble boundary operator, prewitt boundary operator, roberts boundary operator, log boundary operator, zeroscross boundary operator, and the value of L1 norm penalty factor λ is between 0～1.

The incomplete angle method for reconstructing of this Cone-Beam CT is on the minimized basis based on total variation, proposition is with the marginal information of the image method of combination with it, design effective edge detection operator, the marginal information that continuous utilization detects in the iterative process of rebuilding is upgraded the objective function of the Optimized model of reconstruction, makes reconstruction algorithm can adapt to image data still less and promote reconstruction quality.

Beneficial effect of the present invention:

Utilization of the present invention has realized the not high precision image reconstruction of complete angle of Cone-Beam CT based on the sparse Optimized Iterative reconstruction technique of margin guide; Under incomplete angle or sparse sampling, traditional analytical algorithm is rebuild with strong artifact, has a strong impact on the problem of distinguishing of useful information; The technology of sampling iterative approximation is utilized the effective marginal information of intermediate reconstructed images in the process of iteration, proposed the method for reconstructing based on the incomplete angle of Cone-Beam CT of margin guide; From existing to minimize sparse reconstruction method based on total variation different, the present invention combines edge and the gradient sparse prior information of image, choose suitable marginal information every in taking turns iteration, and design suitable weighting factor, utilize weighting factor to upgrade the objective function of optimizing, the objective function after upgrading is carried out obtaining more high-quality reconstructed image based on the iterative strategy of alternating direction method.The present invention utilizes edge-detection algorithm accurately,, in conjunction with efficient optimisation strategy, has promoted the efficiency of convergence and the quality of reconstructed image.

(4), description of drawings:

Fig. 1 is the cone-beam CT scan model;

Fig. 2 is the rim detection example;

Fig. 3 is actual three-dimensional reconstruction result;

Fig. 4 is that convergence curve is rebuild in emulation.

(5), embodiment:

Contain and have the following steps based on the incomplete angle method for reconstructing of Cone-Beam CT of margin guide:

Step 1: estimate initial pictures: utilize the data for projection that scans to estimate initial reconstructed image;

Step 2: Edge extraction;

Step 3: design weighting factor;

Step 4: upgrade Optimized model;

Step 5: based on the incomplete angle of sparse optimization Cone-Beam CT, rebuild;

Step 6: judge that reconstruction quality reaches requirement? in this way, perform step 7; , if not being, perform step 2;

Step 7: finish.

The concrete method of estimation of step 1 is as follows:

Step 1.1: the image reconstruction problem is portrayed as following sparse model:

min||x|| _TV

s.t.Ax=b

Wherein, || x|| _TVFor total variation (TV) norm of object x to be rebuild, A is system matrix, and vectorial b is the data for projection that scans;

Step 1.2: with alternating direction method, the sparse model in step 1.1 is solved, iterations is N, obtains solution formula as follows:

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T({z_j}^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ {z_j}^{k+1}=\max \left\{\right|D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j})\\ {u_j}^{k+1}={u_j}^k+{\mathrm{\ρ}}_j(D_j{x}^{k+1}-{z_j}^{k+1})\end{array}\right.$

Wherein, () ⁺For the Moore-Penrose pseudoinverse of matrix, A ^TFor the transposition of system matrix, u _jFor the renewal factor of j direction, u _j ^kBe the renewal factor of the k time j direction after iteration, u _j ^k+1Be the renewal factor of the k+1 time j direction after iteration, ρ _jFor the L2 norm penalty factor of j direction, λ is L1 norm penalty factor, D _jFor the gradient operator matrix on the j direction, D _j ^TFor the gradient operator transpose of a matrix on the j direction, z _jFor j direction gradient image, z _j ^kBe the k time j direction gradient image after iteration, z _j ^k+1Be the k+1 time j direction gradient image after iteration, max{}, sgn{ } be respectively that to ask maximal function and ask sign function, vectorial b be the data for projection that scans, x ^k+1Be the k+1 time reconstructed image after iteration, through the reconstructed image after N iteration, be denoted as x ^(N), x ^(N)For the last initial reconstructed image of estimating.

Iteration in step 1.2 time N is 0.5～0.2 times of total convergence wheel number.

The concrete grammar of step 2 is: establishing the intermediate reconstructed images that in step 1.2, certain iteration produces is x1, utilizes noise reduction factor F to carry out convolution algorithm: x1 to middle reconstructed image x1 _F=x1***F, utilize classical edge operator E to ask for the preliminary edge x1 of intermediate reconstructed images x1 _E-mid: x1 _E-mid=E[x1], the edge x1 of intermediate reconstructed images x1 _eFor: x1 _e=x1 _E-mid∩ x1 _LMI, wherein, x1 _LMILocal mutual information image for intermediate reconstructed images x1;

The concrete grammar of step 3 is: according to the edge x1 of the intermediate reconstructed images x1 in step 2 _eλ determines the weighting factor w (i1, j1) that intermediate reconstructed images x1 locates at coordinate (i1, j1) with L1 norm penalty factor, this weighting factor w (i1, j1)=λ x1 _e(i1, j1), wherein, x1 _e(i1, j1) is intermediate reconstructed images x1 at the edge that coordinate (i1, j1) is located;

The concrete grammar of step 4 is: the weighting diagonal matrix M that calculates the j direction according to the weighting factor w (i1, j1) in step 3 _j, M _j=diag (w (1,1), w (1,2) ... w (i1, j1) ..., w (N _x, N _y)), Optimized model is updated to following expression so:

$\min {\left|\right|z_1\left|\right|}_1+{\left|\right|z_2\left|\right|}_1+{\left|\right|z_3\left|\right|}_1$

$s.t.\left\{\begin{array}{c}\mathrm{Ax}=b\\ M_j{D}_{1}x=z_1\\ M_j{D}_{2}x=z_2\\ M_j{D}_{3}x=z_3\end{array}\right.$

Wherein, N _xFor the scale of intermediate reconstructed images x1 on the x coordinate, N _yFor the scale of intermediate reconstructed images x1 on the y coordinate, D ₁Be the gradient operator matrix on 1 direction, D ₂Be the gradient operator matrix on 2 directions, D ₃It is the gradient operator matrix on 3 directions;

The concrete grammar of step 5 is: the Optimized model of the belt restraining after adopting the augmentation method of Lagrange multipliers with the renewal in step 4 transfers unconfined Optimized model to, and concrete formula is as follows:

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{j}{\mathrm{\Σ}}(u_j^T(M_j{D}_{j}x-z_j)+\frac{{\mathrm{\ρ}}_j}{2}{{\left|\right|M_j{D}_{j}x-z_j\left|\right|}_2^2+\mathrm{\λ}\left|\right|z_j\left|\right|}_1)$

Wherein,

Transposition for the renewal factor of j direction;

Above-mentioned unconfined Optimized model is adopted variables separation, utilize alternating direction method to ask minimum, iterative formula is as follows:

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M}_j{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M_j}^T(z_j^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ z_j^{k+1}=\max \left\{\right|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}[M_j{D}_j{x}^{k+1}+\frac{u_j^k}{{\mathrm{\ρ}}_j}]\\ u_j^{k+1}=u_j^k+{\mathrm{\ρ}}_j(M_j{D}_j{x}^{k+1}-z_j^{k+1})\end{array}\right.$

x ^k+1Be the epicycle reconstructed image;

Reconstruction quality in step 6 reaches requirement and refers to: the epicycle reconstructed image is compared without marked change with last round of reconstructed image; During reconstruction in carry out step 5 first, last round of reconstructed image refers to the initial reconstructed image in step 1.

Noise reduction factor F is gaussian kernel function or median filter, classical edge operator E is any in canny boundary operator, soble boundary operator, prewitt boundary operator, roberts boundary operator, log boundary operator, zeroscross boundary operator, and the value of L1 norm penalty factor λ is between 0～1.

The incomplete angle method for reconstructing of this Cone-Beam CT is on the minimized basis based on total variation, proposition is with the marginal information of the image method of combination with it, design effective edge detection operator, the marginal information that continuous utilization detects in the iterative process of rebuilding is upgraded the objective function of the Optimized model of reconstruction, makes reconstruction algorithm can adapt to image data still less and promote reconstruction quality.

In order to make in actual applications the Cone-Beam CT under incomplete angle scanning can have high-precision imaging results, adapt to better the demand of practical application, need the Cone-Beam CT data for projection high precision under incomplete angle scanning is rebuild and studied.The achievement in research that the incomplete angle method for reconstructing of this Cone-Beam CT is surveyed based on existing compressive sensing theory and iteration support, adopt advanced sparse optimized algorithm to propose the incomplete angle scanning cone-beam CT reconstruction of high precision algorithm based on margin guide.The incomplete angle method for reconstructing of this Cone-Beam CT is from the initial estimation image, by add edge detection in process of reconstruction, and utilize the institute edge of surveying to improve the objective function of optimization, and by continuous iteration and circulative metabolism, raising reconstruction precision and speed of convergence.

Further illustrate the incomplete angle method for reconstructing of this Cone-Beam CT below by example:

Step 1: estimate initial pictures;

The cone-beam CT scan model as shown in Figure 1, utilizes the Cone-Beam CT data for projection under the incomplete angle scanning that scans, and estimates initial reconstructed image, and method of estimation is:

1. the incomplete pyramidal CT image Problems of Reconstruction under angle scanning:

min||x|| _TV

s.t.Ax=b

Transfer unconstrained problem to:

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}^2+{\mathrm{\λ}\left|\right|x\left|\right|}_{\mathrm{TV}}$

Use ADM method variables separation to obtain:

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}^2+\mathrm{\λ}\left(\underset{i}{\mathrm{\Σ}}{\left|\right|z_i\left|\right|}_1\right)$

s.t.D _ix=z _i

2. we on the basis of alternating direction method framework, provide a kind of concrete reconstruction algorithm for following TV minimum model, namely based on the total variation minimization algorithm of the alternating direction method of augmentation Lagrange.

The argument Lagrange function that it is corresponding:

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}^2+\underset{i}{\mathrm{\Σ}}(u_i^T(D_{i}x-z_i)+\frac{{\mathrm{\ρ}}_i}{2}{\left|\right|D_{i}x-z_i\left|\right|}^2{+\mathrm{\λ}\left|\right|z_i\left|\right|}_1)$

Make (x ^*, z ^*) be L _AMinimum point, multiplier more new formula be:

${\tilde{v}}_i=v_i-{\mathrm{\β}}_i(D_i{x}^{*}-z_i^{*})$

$\widetilde{\mathrm{\λ}}=\mathrm{\λ}-\mathrm{\μ}({\mathrm{Ax}}^{*}-b).$

Below with alternating direction method, solve L _AMinimum point, use

The near-minimizer that represents k wheel iteration.Part about x is:

$\min {\left|\right|\mathrm{Ax}-b\left|\right|}^2+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{\left|\right|D_{i}x-z_i+u_i/{\mathrm{\ρ}}_i\left|\right|}^2\right)$

I.e. " x-subproblem ", this is a problem that quadratic form is minimized, and it is carried out differentiate and makes d _k(x)=0 obtains f _k(x) strict minimum point

$x^{*}={(A^{T}A+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{D}_i^T{D}_i\right))}^{+}(A^{T}b+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{D}_i^T(z_i-u_i/{\mathrm{\ρ}}_i)\right)),$

M wherein ⁺The Moore-Penrose pseudoinverse of expression M.To obtain strict minimum point by asking pseudoinverse to solve the x-subproblem theoretically, yet, be very huge in every numerical evaluation expense of calculating contrary or pseudoinverse in taking turns iteration, therefore usually use alternative manner to solve.

The steepest descending method can iterative x-subproblem, yet for fairly large problem, neither an efficient algorithm.In fact augmented lagrangian function reaches minimum by alternately solving z-subproblem and x-subproblem.Therefore, be unnecessary each step Exact Solution x-subproblem, only need to use a steepest decline result that goes on foot instead.Select to use BB(Barzilai andBorwein in the step-length of descent direction) the type method.

Obtain x ^k+1After solve

Can solve by following subproblem:

$\underset{z_i}{\min }{L}_A(x^k,z_i)=\underset{i}{\mathrm{\Σ}}\left(\right|\left|z_i\right||-v_i^T(D_i{x}^k-z_i)+\frac{{\mathrm{\β}}_i}{2}\·{\left|\right|D_i{x}^k-z_i\left|\right|}_2^2)-{\mathrm{\λ}}^T({\mathrm{Ax}}^k-b)+\frac{\mathrm{\μ}}{2}\·{\left|\right|A{x}^k-b\left|\right|}_2^2$ I.e. " z-subproblem "

$\underset{w_i}{\min }\underset{i}{\mathrm{\Σ}}\left(\right|\left|w_i\right||-v_i^T(D_i{\stackrel{\→}{f}}^k-w_i)+\frac{{\mathrm{\β}}_i}{2}\·{\left|\right|D_i{\stackrel{\→}{f}}^k-w_i\left|\right|}_2^2)$

The z-subproblem can be to each z _iAsk separately minimum, its closed solutions:

z _j ^*=shrinkage(D _jx+u _jρ _j,λ/ρ _j)；

The iterative process of ADM is:

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_i{D}_i^T{D}_i)}^{+}(A^{T}b+\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_i{D}_i^T({z_i}^k-{u_i}^k/{\mathrm{\ρ}}_i))\\ {z_j}^{k+1}=\max \left\{\right|D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j})\\ {u_j}^{k+1}={u_j}^k+{\mathrm{\ρ}}_j(D_j{x}^{k+1}-{z_j}^{k+1})\end{array}\right.$

Wherein, || x|| _TVFor total variation (TV) norm of object x to be rebuild, () ⁺For the Moore-Penrose pseudoinverse of matrix, A is system matrix, and b is observation data, A ^TFor the transposition of system matrix, u _jFor upgrading the factor, ρ _jFor the L2 norm penalty factor of j direction, λ is L1 norm penalty factor, D _jFor the gradient operator matrix on the j direction, z _jFor j direction gradient image, max{}, sgn{} are respectively and ask maximal function and ask sign function.

Step 2: Edge extraction;

Utilize noise reduction factor F to carry out two-dimensional convolution x to image x _F=X**F, F are generally gaussian kernel function or median filter etc.Utilize classical edge operator E, ask for the preliminary edge x of image _E-mid=E[x], E is generally the edge operators such as canny, soble, prewitt, roberts, log, zeroscross.The edge of image is: x _e=x _E-mid∩ x _LMI, wherein, x _LMIFor image local mutual information image.

Illustration is analyzed: as Fig. 2, shown in for the test phantom (a), do not sneaked into more high frequency noise in image before noise reduction, utilize the image (b) after median filter carries out two-dimensional filtering, the obvious reduction of noise can promote the validity of rim detection, after strengthening by local mutual information, the edge image (c) that utilizes classical edge detection operator to detect to obtain pilot process, the edge image of pilot process obtained relatively meeting the marginal information (d) of phantom used.

Step 3: design weighting factor;

Extracted edge x by step 2 _e, the selection percentage coefficient lambda, weighting factor is decided to be w (i, j)=λ x so _e(i, j), wherein the value of scale-up factor λ is generally 0～1.

Step 4: upgrade Optimized model;

By step 3 gained weighting factor w (i, j)=λ x _e(i, j), calculate diagonal matrix

$M=\mathrm{diag}(d_{\mathrm{1,1}},d_{\mathrm{1,2}},...d_{i,j},...,d_{N_x,N_y}),$ d _i,j=λ x _e(i, j), Optimized model is updated to following expression so:

$\min {\left|\right|z_1\left|\right|}_1+{\left|\right|z_2\left|\right|}_1$

$s.t.\left\{\begin{array}{c}\mathrm{Ax}=b\\ {\mathrm{MD}}_{i}x=z_i\end{array}\right.$

Step 5: based on the Cone-Beam CT limited angle image reconstruction under the incomplete angle scanning of sparse optimization;

Adopt the augmentation method of Lagrange multipliers to transfer unconfined Optimized model to the Optimized model of the belt restraining after the renewal in step 4:

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{i}{\mathrm{\Σ}}(u_i^T(M_i{D}_{i}x-z_i)+\frac{{\mathrm{\ρ}}_i}{2}{{\left|\right|M_i{D}_{i}x-z_i\left|\right|}_2^2+\mathrm{\λ}\left|\right|z_i\left|\right|}_1)$

Objective function to unconfined Optimized model carries out the variable separation, resolves into two sub-problems:

(1) X subproblem

$\min {\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{\left|\right|M_i{D}_{i}x-z_i+u_i/{\mathrm{\ρ}}_i\left|\right|}_2^2\right)$

To following formula to x differentiate open order and go out x for null solution: $x^{*}={(A^{T}A+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{M_{i}D}_i^T{D}_i\right))}^{+}(A^{T}b+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{M}_i{D}_i^T(z_i-\frac{u_i}{{\mathrm{\ρ}}_i})\right))$

This separates in form is analytic solution, and in fact computation complexity is still very high, adopts a step steepest to descend when actual the realization.

(2) Z subproblem

$\min u_j^T(D_{j}x-z_j)+\frac{{\mathrm{\ρ}}_j}{2}{\left|\right|D_{j}x-z_j\left|\right|}^2+\mathrm{\λ}{\left|\right|z_j\left|\right|}_1$

Obtain explicit solution by the shrinkage operator:

$z_j^{*}=\mathrm{shrinkage}\{M_j{D}_{j}x+\frac{u_j}{{\mathrm{\ρ}}_j},\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j}\}$

The analytic solution form:

${z_j}^{k+1}=\max \left\{\right|M_j{D}_j{x}^{*}+\frac{u_j}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(M_j{D}_{j}x+\frac{u_j}{{\mathrm{\ρ}}_j})$

In sum, iterative step is as follows:

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_i{D_i}^T{M}_i{D}_i)}^{+}(A^{T}b+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{D_i}^T{M_i}^T(z_i^k-{u_i}^k/{\mathrm{\ρ}}_i)\right))\\ z_j^{k+1}=\max \left\{\right|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}[M_j{D}_j{x}^{k+1}+\frac{u_j^k}{{\mathrm{\ρ}}_j}]\\ u_j^{k+1}=u_j^k+{\mathrm{\ρ}}_j(M_j{D}_j{x}^{k+1}-z_j^{k+1})\end{array}\right.$

Step 6: reconstructed results is shown, check whether rebuild effect meets the demands;

The investigation of picture quality is mainly carried out from the following aspect: 1, the geometry artifact of image, and scatter artefacts, whether the sclerosis artifact is effectively suppressed.2, whether the diplopia of image affects or the distinguishing of interfering picture details.Whether the level of resolution of 3, rebuilding meets application demand.

Utilize this algorithm to rebuild real data, reconstructed results as shown in Figure 3.Utilize the constringency performance of emulated data testing algorithm, as shown in Figure 4.

Claims (5)

1. incomplete angle method for reconstructing of the Cone-Beam CT based on margin guide is characterized in that: contain and have the following steps:

Step 1: estimate initial pictures: utilize the data for projection that scans to estimate initial reconstructed image;

Step 2: Edge extraction;

Step 3: design weighting factor;

Step 4: upgrade Optimized model;

Step 5: based on the incomplete angle of sparse optimization Cone-Beam CT, rebuild;

Step 6: judge that reconstruction quality reaches requirement? in this way, perform step 7; , if not being, perform step 2;

Step 7: finish.

2. the incomplete angle method for reconstructing of the Cone-Beam CT based on margin guide according to claim 1, it is characterized in that: the concrete method of estimation of described step 1 is as follows:

Step 1.1: the image reconstruction problem is portrayed as following sparse model:

min||x|| _TV

s.t.Ax=b

Wherein, || x|| _TVFor the total variation norm of object x to be rebuild, A is system matrix, and vectorial b is the data for projection that scans;

Step 1.2: with alternating direction method, the sparse model in step 1.1 is solved, iterations is N, obtains solution formula as follows:

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T({z_j}^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ {z_j}^{k+1}=\max \left\{\right|D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j})\\ {u_j}^{k+1}={u_j}^k+{\mathrm{\ρ}}_j(D_j{x}^{k+1}-{z_j}^{k+1})\end{array}\right.$

Wherein, () ⁺For the Moore-Penrose pseudoinverse of matrix, A ^TFor the transposition of system matrix, u _jFor the renewal factor of j direction, u _j ^kBe the renewal factor of the k time j direction after iteration, u _j ^k+1Be the renewal factor of the k+1 time j direction after iteration, ρ _jFor the L2 norm penalty factor of j direction, λ is L1 norm penalty factor, D _jFor the gradient operator matrix on the j direction, D _j ^TFor the gradient operator transpose of a matrix on the j direction, z _jFor j direction gradient image, z _j ^kBe the k time j direction gradient image after iteration, z _j ^k+1Be the k+1 time j direction gradient image after iteration, max{}, sgn{} are respectively and ask maximal function and ask sign function, vectorial b is the data for projection that scans, x ^k+1Be the k+1 time reconstructed image after iteration, through the reconstructed image after N iteration, be denoted as x ^(N), x ^(N)For the last initial reconstructed image of estimating.

3. the incomplete angle method for reconstructing of the Cone-Beam CT based on margin guide according to claim 2 is characterized in that: time N of the iteration in described step 1.2 is 0.5～0.2 times of total convergence wheel number.

4. the incomplete angle method for reconstructing of the Cone-Beam CT based on margin guide according to claim 2, it is characterized in that: the concrete grammar of described step 2 is: establishing the intermediate reconstructed images that in step 1.2, certain iteration produces is x1, utilizes noise reduction factor F to carry out convolution algorithm: x1 to middle reconstructed image x1 _F=x1***F, utilize classical edge operator E to ask for the preliminary edge x1 of intermediate reconstructed images x1 _E-mid: x1 _E-mid=E[x1], the edge x1 of intermediate reconstructed images x1 _eFor: x1 _e=x1 _E-mid∩ x1 _LMI, wherein, x1 _LMILocal mutual information image for intermediate reconstructed images x1;

The concrete grammar of step 3 is: according to the edge x1 of the intermediate reconstructed images x1 in step 2 _eλ determines the weighting factor w (i1, j1) that intermediate reconstructed images x1 locates at coordinate (i1, j1) with L1 norm penalty factor, this weighting factor w (i1, j1)=λ x1 _e(i1, j1), wherein, x1 _e(i1, j1) is intermediate reconstructed images x1 at the edge that coordinate (i1, j1) is located;

The concrete grammar of step 4 is: the weighting diagonal matrix M that calculates the j direction according to the weighting factor w (i1, j1) in step 3 _j, M _j=diag (w (1,1), w (1,2) ... w (i1, j1) ..., w (N _x, N _y)), Optimized model is updated to following expression so:

$\min {\left|\right|z_1\left|\right|}_1+{\left|\right|z_2\left|\right|}_1+{\left|\right|z_3\left|\right|}_1$

$s.t.\left\{\begin{array}{c}\mathrm{Ax}=b\\ M_j{D}_{1}x=z_1\\ M_j{D}_{2}x=z_2\\ M_j{D}_{3}x=z_3\end{array}\right.$

Wherein, N _xFor the scale of intermediate reconstructed images x1 on the x coordinate, N _yFor the scale of intermediate reconstructed images x1 on the y coordinate, D ₁Be the gradient operator matrix on 1 direction, D ₂Be the gradient operator matrix on 2 directions, D ₃It is the gradient operator matrix on 3 directions;

The concrete grammar of step 5 is: the Optimized model of the belt restraining after adopting the augmentation method of Lagrange multipliers with the renewal in step 4 transfers unconfined Optimized model to, and concrete formula is as follows:

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{j}{\mathrm{\Σ}}(u_j^T(M_j{D}_{j}x-z_j)+\frac{{\mathrm{\ρ}}_j}{2}{{\left|\right|M_j{D}_{j}x-z_j\left|\right|}_2^2+\mathrm{\λ}\left|\right|z_j\left|\right|}_1)$

Wherein,

Transposition for the renewal factor of j direction;

Above-mentioned unconfined Optimized model is adopted variables separation, utilize alternating direction method to ask minimum, iterative formula is as follows:

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M}_j{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M_j}^T(z_j^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ z_j^{k+1}=\max \left\{\right|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}[M_j{D}_j{x}^{k+1}+\frac{u_j^k}{{\mathrm{\ρ}}_j}]\\ u_j^{k+1}=u_j^k+{\mathrm{\ρ}}_j(M_j{D}_j{x}^{k+1}-z_j^{k+1})\end{array}\right.$

x ^k+1Be the epicycle reconstructed image;

Reconstruction quality in step 6 reaches requirement and refers to: the epicycle reconstructed image is compared without marked change with last round of reconstructed image; During reconstruction in carry out step 5 first, last round of reconstructed image refers to the initial reconstructed image in step 1.

5. the incomplete angle method for reconstructing of the Cone-Beam CT based on margin guide according to claim 4, it is characterized in that: described noise reduction factor F is gaussian kernel function or median filter, classical edge operator E is any in canny boundary operator, soble boundary operator, prewitt boundary operator, roberts boundary operator, log boundary operator, zeroscross boundary operator, and the value of L1 norm penalty factor λ is between 0～1.

Images