CN103390285B - The incomplete angle reconstruction method of Cone-Beam CT based on margin guide - Google Patents

Abstract

The present invention relates to a kind of incomplete angle reconstruction method of Cone-Beam CT based on margin guide, specifically containing following steps: step 1: estimate initial pictures: utilize the data for projection scanned 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 reconstruction of sparse optimization Cone-Beam CT; Step 6: judge that reconstruction quality reaches requirement? in this way, then step 7 is performed; If not, then perform step 2; Step 7: terminate; The invention provides a kind of incomplete angle reconstruction of the Cone-Beam CT based on the margin guide method that efficiency is high, reconstructed image quality is good.

Description

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

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

(2), background technology: Computed tomography (ComputedTomography, CT) has been widely used in the field such as medical science, industry as a kind of modern imaging techniques.But, in a lot of practical application, owing to retraining by the geometric position of data acquisition time or imaging system scans, can only obtain data at incomplete angular range or at less projection angle, these all belong to incomplete angle problem (IncompleteDataProblem).Be lifted at the quality of CT image reconstruction under incomplete angle scanning and have important theoretical research and practical engineering value, the method for CT image reconstruction under high-precision incomplete angle that how to design also is focus and the difficulties of research.

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

Although can obtain good reconstructed results under comparatively sparse sampling based on the reconstruction algorithm of CS theoretical model, but not excavating to a deeper level can advantageously in the image information improving reconstructed results.The marginal element of image is the important elements of image, has important application in the multiple image applications field such as image repair, segmentation.In image reconstruction, effectively utilize marginal information can improve reconstruction quality largely.

For incomplete angle CT image reconstruction problem, classical disposal route is algebraically iteration and statistics iterative reconstruction algorithm, such as algebraic reconstruction algorithm (AlgebraicReconstructionTechnique, ART) and expect maximum (expectationmaximization, EM) algorithm, iterative algorithm has good reconstruction quality than analytical algorithm usually when processing missing data and being less, has good reconstruction quality as 1984 by the FDK algorithm of people's propositions such as Feldkamp.In recent years, based under the guidance of compressive sensing theory, having grown up, some can be applicable to the reconstruction algorithm of sparse sampling, comparatively typically have ASD-POCS (adapt-steepest-descentprojectionontoconvexsets) algorithm that the people such as Sidky and Pan in 2008 propose.This algorithm, with TV norm minimum design optimization model, adopts POCS (ART) and TV steepest to decline the strategy alternately performed, can reconstruct the image of better quality under sparse sampling.

After ASD-POCS algorithm is suggested, a lot of scholar also been proposed multiple optimized algorithm on model solution, the Split-BregmanTV algorithm proposed by people such as Vandeghinste for such as 2011, the ADTVM algorithm etc. proposed by people such as Zhang for 2012.Other algorithms introduce other priori etc., such as PICCS (priorimageconstrainedcompressedsensing, PICCS), within 2010, Wang Lin unit waits the RRD (reconstruction-referencedifference of people's proposition, RRD) scheduling algorithm is all demand object to be reconstructed sparse expressions under certain priori, utilize sparse expression constitution optimization model and design corresponding derivation algorithm, thus reaching the object of image reconstruction.The limited angle of frequency domain aspect rebuilds (main MRI imaging) comparatively strong instrument is in addition RecPF (reconstructionfrompartialFourierspacesmpling) algorithm, this algorithm does not relate to storage and the calculating of extensive matrix, utilizes difference operator D _igood character and FFT technology realize rebuilding efficient and accurately, be the theoretical example in MRI successful Application of CS.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 method that the people such as Yin propose to utilize iteration support to detect (iterativesupportdetection, ISD) achieved preferably reconstructed results in sparse signal recovers.In the method, constantly to iterative process M signal do support set detection process, the objective function propping up set pair Optimized model according to detected part upgrades, take this alternative manner can on the basis of extreme lack sampling effective restoring signal.The pure method based on support detection can not directly use in the image reconstruction of Cone-Beam CT, and the support set under utilizing image sparse to represent to integrate can bring the lifting of reconstruction quality as CT image reconstruction.

(3), summary of the invention:

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

Technical scheme of the present invention:

The incomplete angle reconstruction method of Cone-Beam CT based on margin guide, containing following steps:

Step 1: estimate initial pictures: utilize the data for projection scanned 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 reconstruction of sparse optimization Cone-Beam CT;

Step 6: judge whether reconstruction quality reaches requirement, in this way, then performs step 7; If not, then perform step 2;

Step 7: terminate.

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

Step 1.1: 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 reconstructed, A is system matrix, and vectorial b is the data for projection scanned;

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

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\Σ}{\mathrm{\ρ}}_j{D}_j^T{D}_j)}^{+}(A^{T}b+\underset{j}{\Σ}{\mathrm{\ρ}}_j{D}_j^T\left({z_j}^k-{u_j}^k/{\mathrm{\ρ}}_j\right))\\ {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 in j direction, j represents three dimension directions of 3-D view, u _j ^kfor the renewal factor in the j direction after kth time iteration, u _j ^k+1for the renewal factor in the j direction after kth+1 iteration, ρ _jfor the L2 norm penalty factor in j direction, λ is L1 norm penalty factor, D _jfor the gradient operator matrix on j direction, D _j ^tfor the gradient operator transpose of a matrix on j direction, z _jfor j direction gradient image, z _j ^kfor the j direction gradient image after kth time iteration, z _j ^k+1for the j direction gradient image after kth+1 iteration, max{}, sgn{} are respectively and ask maximal function and ask sign function, and vectorial b is the data for projection scanned, x ^k+1for the reconstruction image after kth+1 iteration, the reconstruction image after N iteration is denoted as x ^(N), x ^(N)for the initial reconstructed image finally estimated.

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: set the intermediate reconstructed images of certain iteration generation in step 1.2 as x1, utilize noise reduction factor F to carry out convolution algorithm to intermediate reconstructed images x1: x1 _f=x1*F, utilizes classical edge operator E to ask for the preliminary edge x1 of intermediate reconstructed images x1 _e-mid: x1 _e-midthe edge x1 of=E [x1], intermediate reconstructed images x1 _efor: x1 _e=x1 _e-mid∩ x1 _lMI, wherein, x1 _lMIfor the local mutual information image of 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 _ethe weighting factor w (i1, j1) of intermediate reconstructed images x1 at coordinate (i1, j1) place is determined, this weighting factor w (i1, j1)=λ x1 with L1 norm penalty factor λ _e(i1, j1), wherein, x1 _e(i1, j1) is for intermediate reconstructed images x1 is at the edge at coordinate (i1, j1) place;

The concrete grammar of step 4 is: the weighting diagonal matrix M calculating j direction according to the weighting factor w (i1, j1) in step 3 _j, $M_j=diag\left(w{(}_{1,1}\right),w{(}_{1,2}),...w{(}_{i1,j1}),...,w{(}_{N_x,N_y}\left)\right),$ So Optimized model is updated to following expression:

min||z ₁|| ₁+||z ₂|| ₁+||z ₃|| ₁

$s.t.\left\{\begin{array}{c}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 x coordinate, N _yfor the scale of intermediate reconstructed images x1 on 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: adopt augmented vector approach to transfer the Optimized model of the belt restraining after the renewal in step 4 to unconfined Optimized model, concrete formula is as follows:

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

Wherein, for the transposition of the renewal factor in j direction;

Adopt variables separation to above-mentioned unconfined Optimized model, 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{\ρ}}_j{D_j}^T{M}_j{D}_j)}^{+}(A^{T}b+\underset{j}{\Σ}{\mathrm{\ρ}}_j{D_j}^T{M_j}^T\left(z_j^k-{u_j}^k/{\mathrm{\ρ}}_j\right))\\ z_j^{k+1}=max\{\left|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}\right|-\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 epicycle and rebuild image;

Reconstruction quality in step 6 reaches requirement and refers to: epicycle rebuilds image compared with last round of reconstruction image without marked change; During reconstruction in carry out step 5 first, last round of reconstruction 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 one 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.

This Cone-Beam CT incomplete angle reconstruction method is based on the minimized basis of total variation, the method marginal information of image combined with it is proposed, design effective edge detection operator, in the iterative process of rebuilding, constantly utilize the objective function of marginal information to the Optimized model rebuild detected to upgrade, make reconstruction algorithm can adapt to less image data and promote reconstruction quality.

Beneficial effect of the present invention:

The high precision image that the present invention utilizes the sparse Optimized Iterative reconstruction technique based on margin guide to achieve the not complete angle of Cone-Beam CT is rebuild; For under incomplete angle or sparse sampling, conventional analytic algorithm is rebuild with strong artifact, has a strong impact on the problem distinguished of useful information; The technology of sampling iterative approximation utilizes the effective marginal information of intermediate reconstructed images in the process of iteration, proposes the method for reconstructing of the incomplete angle of Cone-Beam CT based on margin guide; From existing to minimize sparse reconstruction method based on total variation different, present invention incorporates edge and the gradient sparse prior information of image, choose suitable marginal information often taking turns in iteration, and design suitable weighting factor, utilize weighting factor to upgrade the objective function optimized, carry out obtaining more high-quality reconstruction image based on the iterative strategy of alternating direction method to the objective function after upgrading.The present invention utilizes edge-detection algorithm accurately, in conjunction with efficient optimisation strategy, improves the efficiency of convergence and rebuilds the quality of image.

(4), accompanying drawing illustrates:

Fig. 1 is cone-beam CT scan model;

Fig. 2 is rim detection example;

Fig. 3 is actual three-dimensional reconstruction result;

Fig. 4 is emulation reconstruction convergence curve.

(5), embodiment:

The incomplete angle reconstruction method of Cone-Beam CT based on margin guide contains following steps:

Step 1: estimate initial pictures: utilize the data for projection scanned 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 reconstruction of sparse optimization Cone-Beam CT;

Step 6: judge whether reconstruction quality reaches requirement, in this way, then performs step 7; If not, then perform step 2;

Step 7: terminate.

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

Step 1.1: 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 reconstructed, A is system matrix, and vectorial b is the data for projection scanned;

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

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\Σ}{\mathrm{\ρ}}_j{D}_j^T{D}_j)}^{+}(A^{T}b+\underset{j}{\Σ}{\mathrm{\ρ}}_j{D}_j^T\left({z_j}^k-{u_j}^k/{\mathrm{\ρ}}_j\right))\\ {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 in j direction, j represents three dimension directions of 3-D view, u _j ^kfor the renewal factor in the j direction after kth time iteration, u _j ^k+1for the renewal factor in the j direction after kth+1 iteration, ρ _jfor the L2 norm penalty factor in j direction, λ is L1 norm penalty factor, D _jfor the gradient operator matrix on j direction, D _j ^tfor the gradient operator transpose of a matrix on j direction, z _jfor j direction gradient image, z _j ^kfor the j direction gradient image after kth time iteration, z _j ^k+1for the j direction gradient image after kth+1 iteration, max{}, sgn{} are respectively and ask maximal function and ask sign function, and vectorial b is the data for projection scanned, x ^k+1for the reconstruction image after kth+1 iteration, the reconstruction image after N iteration is denoted as x ^(N), x ^(N)for the initial reconstructed image finally estimated.

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: set the intermediate reconstructed images of certain iteration generation in step 1.2 as x1, utilize noise reduction factor F to carry out convolution algorithm to intermediate reconstructed images x1: x1 _f=x1*F, utilizes classical edge operator E to ask for the preliminary edge x1 of intermediate reconstructed images x1 _e-mid: x1 _e-midthe edge x1 of=E [x1], intermediate reconstructed images x1 _efor: x1 _e=x1 _e-mid∩ x1 _lMI, wherein, x1 _lMIfor the local mutual information image of 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 _ethe weighting factor w (i1, j1) of intermediate reconstructed images x1 at coordinate (i1, j1) place is determined, this weighting factor w (i1, j1)=λ x1 with L1 norm penalty factor λ _e(i1, j1), wherein, x1 _e(i1, j1) is for intermediate reconstructed images x1 is at the edge at coordinate (i1, j1) place;

The concrete grammar of step 4 is: the weighting diagonal matrix M calculating j direction according to the weighting factor w (i1, j1) in step 3 _j, $M_j=diag\left(w{(}_{1,1}\right),w{(}_{1,2}),...w{(}_{i1,j1}),...,w{(}_{N_x,N_y}\left)\right),$ So Optimized model is updated to following expression:

min||z ₁|| ₁+||z ₂|| ₁+||z ₃|| ₁

$s.t.\left\{\begin{array}{c}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 x coordinate, N _yfor the scale of intermediate reconstructed images x1 on 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: adopt augmented vector approach to transfer the Optimized model of the belt restraining after the renewal in step 4 to unconfined Optimized model, concrete formula is as follows:

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

Wherein, for the transposition of the renewal factor in j direction;

Adopt variables separation to above-mentioned unconfined Optimized model, 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{\ρ}}_j{D_j}^T{M}_j{D}_j)}^{+}(A^{T}b+\underset{j}{\Σ}{\mathrm{\ρ}}_j{D_j}^T{M_j}^T\left(z_j^k-{u_j}^k/{\mathrm{\ρ}}_j\right))\\ z_j^{k+1}=max\{\left|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}\right|-\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 epicycle and rebuild image;

Reconstruction quality in step 6 reaches requirement and refers to: epicycle rebuilds image compared with last round of reconstruction image without marked change; During reconstruction in carry out step 5 first, last round of reconstruction 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 one 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.

This Cone-Beam CT incomplete angle reconstruction method is based on the minimized basis of total variation, the method marginal information of image combined with it is proposed, design effective edge detection operator, in the iterative process of rebuilding, constantly utilize the objective function of marginal information to the Optimized model rebuild detected to upgrade, make reconstruction algorithm can adapt to less image data and promote reconstruction quality.

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

The incomplete angle reconstruction method of this Cone-Beam CT is further illustrated below by example:

Step 1: estimate initial pictures;

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

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

min||x|| _TV

s.t.Ax＝b

Transfer unconstrained problem to:

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

ADM method variables separation is used to obtain:

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

s.t.D _ix＝z _i

2. we are 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 of its correspondence:

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

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

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

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

Alternating direction method is used to solve L below _aminimum point, use represent the near-minimizer of kth wheel iteration.Part about x is:

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

I.e. " x-subproblem ", this is the problem that a quadratic form is minimized, and carries out differentiate and make d to it _kx ()=0 obtains f _kthe strict minimum point of (x)

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

Wherein M ⁺represent the Moore-Penrose pseudoinverse of M.Be to obtain strict minimum point by asking pseudoinverse to solve x-subproblem theoretically, but to calculate numerical evaluation expense that is inverse or pseudoinverse be very huge often taking turns in iteration, therefore usually uses alternative manner to solve.

Steepest descending method can iterative x-subproblem, but for fairly large problem, neither an efficient algorithm.In fact augmented lagrangian function is by alternately solving z-subproblem and x-subproblem reaches minimum.Therefore, be unnecessary each step Exact Solution x-subproblem, only need the steepest decline result of use one step instead.Select to use BB (BarzilaiandBorwein) type method in the step-length of descent direction.

Obtain x ^k+1after solve solve by following subproblem:

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

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

Z-subproblem can to each z _iask separately minimum, its closed solutions:

${z_j}^{*}=shrinkage(D_{j}x+u_j/{\mathrm{\ρ}}_j,\mathrm{\λ}/{\mathrm{\ρ}}_j);$

The iterative process of ADM is:

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{i}{\Σ}{\mathrm{\ρ}}_i{D}_i^T{D}_i)}^{+}(A^{T}b+\underset{i}{\Σ}{\mathrm{\ρ}}_i{D}_i^T\left({z_i}^k-{u_i}^k/{\mathrm{\ρ}}_i\right))\\ {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 reconstructed, () ⁺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 in j direction, λ is L1 norm penalty factor, D _jfor the gradient operator matrix on j direction, z _jfor j direction gradient image, max{}, sgn{} are respectively and ask maximal function and ask sign function.

Step 2: Edge extraction;

Noise reduction factor F is utilized 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 image preliminary edge x _e-mid=E [x], E is generally the boundary 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 for test body mould (a), more high frequency noise has been mixed in image before non-noise reduction, image (b) after utilizing median filter to carry out two-dimensional filtering, noise obviously reduces the validity that can promote rim detection, utilize classical edge detection operator to carry out detecting the edge image (c) obtaining pilot process, the edge image of pilot process is by obtaining the marginal information (d) comparing and meet body mould used after the mutual information enhancing of local.

Step 3: design weighting factor;

By step 2 extract edge x _e, selection percentage coefficient lambda, so weighting factor is decided to be w (i, j)=λ x _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), calculates diagonal matrix $M=diag(d_{1,1},d_{1,2},...d_{i,j},...,d_{N_x,N_y}),$ D _i,j=λ x _e(i, j), so Optimized model is updated to following expression:

min||z ₁|| ₁+||z ₂|| ₁

$s.t.\left\{\begin{array}{c}Ax=b\\ M{D}_{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;

Augmented vector approach is adopted 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|Ax-b|{|}_2^2+\underset{i}{\Σ}\left(u_i^T\right(M_i{D}_{i}x-z_i)+\frac{{\mathrm{\ρ}}_i}{2}||M_i{D}_{i}x-z_i|{|}_2^2+\mathrm{\λ}\left|\right|z_i\left|{|}_1\right)$

Separating variables is carried out to the objective function of unconfined Optimized model, resolves into two subproblems:

(1) X subproblem

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

Also make as null solution goes out x to above formula to x differentiate:

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

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

(2) Z subproblem

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

Explicit solution is obtained by shrinkage operator:

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

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{\ρ}}_i{D_i}^T{M}_i{D}_i)}^{+}(A^{T}b+\underset{i}{\Σ}{\mathrm{\ρ}}_i{D_i}^T{M_i}^T\left(z_i^k-{u_i}^k/{\mathrm{\ρ}}_i\right))\\ z_j^{k+1}=max\{\left|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}\right|-\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: show reconstructed results, checks 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, scatter artefacts, whether hardening 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 (2)

1., based on the incomplete angle reconstruction method of Cone-Beam CT of margin guide, it is characterized in that: containing following steps:

Step 1: estimate initial pictures: utilize the data for projection scanned 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 reconstruction of sparse optimization Cone-Beam CT;

Step 6: judge whether reconstruction quality reaches requirement, in this way, then performs step 7; If not, then perform step 2;

Step 7: terminate;

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

Step 1.1: 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 reconstructed, A is system matrix, and vectorial b is the data for projection scanned;

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

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\Σ}{\mathrm{\ρ}}_j{D}_j^T{D}_j)}^{+}(A^{T}b+\underset{j}{\Σ}{\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\mathrm{+1}}={u_j}^k+{\mathrm{\ρ}}_j(D_j{x}^{k+1}x-{z_j}^{k\mathrm{+1}})\end{array}\right.$

Wherein, () ⁺for the Moore-Penrose pseudoinverse of matrix, A ^tfor the transposition of system matrix, u _jfor the renewal factor in j direction, j represents three dimension directions of 3-D view, u _j ^kfor the renewal factor in the j direction after kth time iteration, u _j ^k+1for the renewal factor in the j direction after kth+1 iteration, ρ _jfor the L2 norm penalty factor in j direction, λ is L1 norm penalty factor, D _jfor the gradient operator matrix on j direction, D _j ^tfor the gradient operator transpose of a matrix on j direction, z _jfor j direction gradient image, z _j ^kfor the j direction gradient image after kth time iteration, z _j ^k+1for the j direction gradient image after kth+1 iteration, max{}, sgn{} are respectively and ask maximal function and ask sign function, and vectorial b is the data for projection scanned, x ^k+1for the reconstruction image after kth+1 iteration, the reconstruction image after N iteration is denoted as x ^(N), x ^(N)for the initial reconstructed image finally estimated;

The concrete grammar of step 2 is: set the intermediate reconstructed images of certain iteration generation in step 1.2 as x1, utilize noise reduction factor F to carry out convolution algorithm to intermediate reconstructed images x1: x1 _f=x1*F, utilizes classical edge operator E to ask for the preliminary edge x1 of intermediate reconstructed images x1 _e-mid: x1 _e-midthe edge x1 of=E [x1], intermediate reconstructed images x1 _efor: x1 _e=x1 _e-mid∩ x1 _lMI, wherein, x1 _lMIfor the local mutual information image of 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 _ethe weighting factor w (i1, j1) of intermediate reconstructed images x1 at coordinate (i1, j1) place is determined, this weighting factor w (i1, j1)=λ x1 with L1 norm penalty factor λ _e(i1, j1), wherein, x1 _e(i1, j1) is for intermediate reconstructed images x1 is at the edge at coordinate (i1, j1) place;

The concrete grammar of step 4 is: the weighting diagonal matrix M calculating j direction according to the weighting factor w (i1, j1) in step 3 _j, so Optimized model is updated to following expression:

min||z ₁|| ₁+||z ₂|| ₁+||z ₃|| ₁

$s.t.\left\{\begin{array}{c}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 x coordinate, N _yfor the scale of intermediate reconstructed images x1 on 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: adopt augmented vector approach to transfer the Optimized model of the belt restraining after the renewal in step 4 to unconfined Optimized model, concrete formula is as follows:

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

Wherein, for the transposition of the renewal factor in j direction;

Adopt variables separation to above-mentioned unconfined Optimized model, 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{\ρ}}_j{D_j}^T{D}_j)}^{+}(A^{T}b+\underset{j}{\Σ}{\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\mathrm{+1}}=u_j^k+{\mathrm{\ρ}}_j(M_j{D}_j{x}^{k+1}x-z_j^{k+1})\end{array}\right.$

X ^k+1be epicycle and rebuild image;

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

2. the incomplete angle reconstruction method of the Cone-Beam CT based on margin guide according to claim 1, it is characterized in that: described noise reduction factor F is gaussian kernel function or median filter, classical edge operator E is any one 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.