CN103218795B - Based on the part K spatial sequence image reconstructing method of self-adaptation doubledictionary study - Google Patents

Abstract

The invention discloses a kind of part K spatial sequence image reconstructing method based on the study of self-adaptation doubledictionary, mainly solve the existing method problem that reconstructed image quality suppression ratio is more serious when 10 times of down-samplings.Its key step is: collecting part K space data, utilizes the correlativity between these part K space data, synthesizes complete K space data, obtains training image by complete K space data; With KSVD algorithm, training is carried out to training image again and obtain high-resolution and low-resolution dictionary; Utilize the part K space data of the relation between high-resolution and low-resolution dictionary to input to be reconstructed, and residual compensation is carried out to reconstructed image obtain reconstruction result more accurately.The present invention effectively can improve reconstructed image quality under the condition of 10 times of down-samplings, can be used for the MRI sequence image reconstruct at multiple position.

Description

Based on the part K spatial sequence image reconstructing method of self-adaptation doubledictionary study

Technical field

The invention belongs to technical field of image processing, relate to the method for Medical Image Processing, can be used for the MRI Image Reconstruction at multiple position.

Background technology

Partial K space image reconstruct is to accelerate a kind of data acquisition amount that reduces that magnetic resonance image (MRI) image taking speed proposes to reconstruct the problem of high-definition picture.In order to address this problem, a lot of classical method is had to be suggested:

The first is the most frequently used zero padding method, and the K space data namely do not gathered is filled up with zero, and then do the formation method that Fourier inversion obtains image space, this formation method can improve image taking speed, and defect has artifact in image.

The second is Phase Correction Method, and the phase place in these class methods hypothesis magnetic resonance image (MRI) space is slow variable condition.It estimates phase place with the image of part low frequency K spatial data reconstruct, and for phase correction, thus reach the object utilizing symmetry to supplement the K space data do not gathered.The POCS method that Stark, H. propose is exactly modal method in these class methods.But the condition slowly changed due to magnetic resonance image (MRI) phase place is usually difficult to meet in whole image space, causes phase estimation error large, cause larger reconstructed error, so that up to the present cannot apply in clinical medicine.

The third is Signal estimation method, and this approach application Signal estimation is theoretical, utilizes the interpolation of part K space data, spreads to the methods such as multiple-objection optimization to obtain the K space data do not gathered outward.M.Funderer proposed the method for image maximum likelihood in 1989, the constrained procedure that E.M.Haacke proposed in the same year, and be all the typical method of these class methods, this class methods imaging effect is far inferior to method for correcting phase.

4th kind is singular spectrum modelling, the thought of these class methods is the facts that can represent by the weighted sum of singular function based on any signal, set up new image expression model, by part K space data extraction model parameter, then reconstruct complete K space by model and parameter.This method is in most cases better than Phase Correction Method and Signal estimation method.

5th kind is the method for compressed sensing, and the method utilizes small echo, finite difference, and dictionary learning etc. carry out rarefaction representation to reconstructed image, and the effect of these class methods is in most cases all better than additive method.

Above-mentioned part K Space Reconstruction method all needs sampling rate to more than 30% usually, could obtain good quality reconstruction.But sampling rate is higher, required acquisition time is longer, and the person of being imaged will be trapped in Image-forming instrument for a long time, and the motion due to the person of being imaged will make image produce motion blur.But, reduce further sampling rate and reconstructed image can be caused to produce artifact and the low problem of image resolution ratio.

Summary of the invention

The object of the invention is to the deficiency for above-mentioned prior art, propose a kind of part K spatial sequence image reconstructing method based on the study of self-adaptation doubledictionary, to reduce data sampling rate, improve the quality of reconstructed image.

For achieving the above object, the present invention includes following steps:

(1) N width part K space data is gathered, with the K space data Q that this N width part K space data synthesis n is complete _i, i=1,2 ..., n; To Q _icarry out 10 times of down-samplings, obtain corresponding part K space data P _i; To Q _imake Fourier inversion, obtain high-definition picture H _i, to P _imake Fourier inversion, obtain low-resolution image L _i, this n to high-definition picture H _iwith low-resolution image L _ias training image;

(2) high resolving power training image H is inputted respectively _iwith low resolution training image L _i, and adopt nonoverlapping mode to get the fritter of 4 × 4 to every width training image, obtain initial high resolution dictionary H and initial low resolution dictionary L;

(3) utilize KSVD algorithm to train initial high resolution dictionary H and initial low resolution dictionary L, obtain new high resolving power dictionary D _hwith new low-resolution dictionary D _l, and high-definition picture H _isparse coefficient α _hiwith low-resolution image L _isparse coefficient α _li;

(4) part K space data P to be reconstructed is inputted _t, to this part K space data P _tadopt the process of zero padding method, obtain initial reconstructed image L _t,

(5) low-resolution dictionary D is utilized _lwith initial reconstructed image L _t, solve initial reconstructed image L _tsparse coefficient α _l;

(6) initial reconstructed image L is asked respectively _twith n width low resolution training image L _ierror: obtain initial reconstructed image L _twith the jth width training image L in n width low resolution training image _jleast error: $e{r}_j=\underset{i=1}{\overset{n}{\min }}\left\{e{r}_i\right\};$

(7) least error er is judged _jwhether be less than threshold value σ=0.1 of setting, if error e r _jbe less than threshold value σ, obtain high-definition picture H to be reconstructed _t' sparse coefficient α _h; If error e r _jbe greater than threshold value, return step (1), Resurvey N width part K space data, upgrade dictionary;

(8) high resolving power dictionary D is utilized _hwith high-definition picture H to be reconstructed _t' sparse coefficient α _h, try to achieve high-definition picture: H _t'=D _h* α _h; Again to changing high-definition picture H _t' carry out residual compensation, obtain final reconstructed image H _t.

The present invention has the following advantages compared with prior art:

1. the present invention utilizes the correlativity between part K spatial sequence data to synthesize several complete K space data, training image is obtained by these complete K space data, and utilize these training images to train dictionary, thus the quantity of information that comprises of dictionary is abundanter, can reconstruct the detailed information of image preferably;

2. the present invention is owing to upgrading dictionary during changing greatly between sequence image, makes the adaptivity of dictionary stronger, improves the robustness of reconstruct, namely all can obtain good quality reconstruction to 1800 width sequence datas;

Simulation result shows, the present invention just can carry out high-quality reconstruct to MRI image under the condition of 10 times of down-samplings, reduces data sampling rate, shortens data acquisition time.

Accompanying drawing explanation

Fig. 1 is general flow chart of the present invention;

Fig. 2 is with the quality reconstruction figure of the present invention to the 100th width chest test pattern;

Fig. 3 is with the quality reconstruction figure of the present invention to the 70th width belly test pattern.

Specific implementation method

With reference to accompanying drawing 1, concrete steps of the present invention comprise as follows:

Step 1. synthesizes complete K space data, obtains training image

1a) gather N width part K space data, with the K space data Q that this N width part K space data synthesis n is complete _i, i=1,2 ..., n, the method for generated data has based on the synthetic method of Pixel-level, the synthetic method of feature based level and the synthetic method etc. based on decision level, and this example adopts but the synthetic method be not limited to based on Pixel-level, and its building-up process is as follows:

A width K space data is synthesized, with the first width part K space data of this N/n width part K space data for standard by N/n width part K space data;

What collect in the second width part K space, and the data that the first width part K space does not collect are added to the first width part K spatially;

What collect in the 3rd width part K space, and the data that front two width part K spaces all do not collect are added to the first width part K spatially, by that analogy, synthesize n complete K space data Q _i.

1b) to above-mentioned spatial data Q _icarry out 10 times of down-samplings, obtain corresponding part K space data P _i; To Q _imake Fourier inversion, obtain high-definition picture H _i, to P _imake Fourier inversion, obtain low-resolution image L _i, this n to high-definition picture H _iwith low-resolution image L _ias training image.

Step 2. pair training image carries out pre-service

Input high resolving power training image H respectively _iwith low resolution training image L _i, and adopt nonoverlapping mode to get the fritter of 4 × 4 to every width training image, obtain initial high resolution dictionary H and initial low resolution dictionary L.

Step 3. trains high-resolution and low-resolution dictionary

Initial high resolution dictionary H and initial low resolution dictionary L is trained, obtains new high resolving power dictionary D _hwith new low-resolution dictionary D _l, and high-definition picture H _isparse coefficient α _hiwith low-resolution image L _isparse coefficient α _li, the method for training dictionary mainly contains two kinds, and be principal component analysis (PCA) PCA and K singular value decomposition method KSVD respectively, this example uses KSVD algorithm to train, and its training process is as follows:

3a) to total optimization formula of KSVD algorithm: $\min \left\{{\left|\right|Y-\mathrm{DX}\left|\right|}_F^2\right\}\mathrm{Subject\; to}\∀l,{\left|\right|X_l\left|\right|}_0\≤T_0,$ Be out of shape, be about to optimization formula wherein be deformed into:

${\left|\right|Y-\mathrm{DX}\left|\right|}_F^2={\left|\right|Y-\underset{m=1}{\overset{k}{\mathrm{\Σ}}}{d}_m{x}_T^m\left|\right|}_F^2={\left|\right|(Y-\underset{m\≠k}{\mathrm{\Σ}}{d}_m{x}_T^m)-d_k{x}_T^k\left|\right|}_F^2={\left|\right|E_k-d_k{x}_T^k\left|\right|}_F^2,$

Wherein, Y is the initial dictionary of input, and D is target training dictionary, and X is Its Sparse Decomposition matrix, for any l row, ‖ X _l‖ ₀for X _l0 norm, for solving 2 norms of Y-DX, T ₀for degree of rarefication control coefrficient; d _mfor the m row atom of D, for the m of X is capable, K is total columns of D, d _kfor the kth row atom of target training dictionary D, for the row k of X, E _kfor not using the kth row atom d of D _kcarry out the error matrix that signal Its Sparse Decomposition produces;

3b) give the optimization formula after distortion be multiplied by matrix Ω _k=P*| ω _k|, obtain goal decomposition formula ${\left|\right|E_k{\mathrm{\Ω}}_k-d_k{x}_T^k{\mathrm{\Ω}}_k\left|\right|}_F^2={\left|\right|E_k^R-d_k{x}_R^k\left|\right|}_F^2,$

Wherein Ω _ksize be P*| ω _k|, P is the columns of the initial dictionary Y of input, | ω _k| be ω _kmodulus value, and Ω _kat (ω _k(m), m) place is 1, other place is 0, wherein 1≤m entirely≤| ω _k|, ω _km () is ω _km number;

3c) to goal decomposition formula in error matrix carry out decomposition of singular matrix to obtain wherein U is left singular matrix, V ^tfor right singular matrix, Φ is singular value matrix;

3d) get k=1 successively, 2 ..., K, with the kth row atom of the first row of left singular matrix U more fresh target train word allusion quotation D, tries to achieve the dictionary D ' after renewal, obtains new high resolving power dictionary D _hwith new low-resolution dictionary D _l;

3e) utilize the initial dictionary Y of input and the dictionary D ' after upgrading, try to achieve Its Sparse Decomposition matrix X ', obtain high-definition picture H _isparse coefficient α _hiwith low-resolution image L _isparse coefficient α _li.

The part K space data that step 4. pre-service is to be reconstructed

4a) input part K space data P to be reconstructed _t, to this part K space data P _tcarry out pre-service, obtain initial reconstructed image L _t, pretreated method can use zero padding method, Phase Correction Method, Signal estimation method etc., this example by zero padding method to this part K space data P _tin after the data padding that do not gather, then carry out Fourier inversion, obtain initial reconstructed image L _t;

4b) utilize low-resolution dictionary D _lwith initial reconstructed image L _t, solve initial reconstructed image L _tsparse coefficient α _l, solve sparse coefficient α _lcan use matching pursuit algorithm, base tracing algorithm, orthogonal matching pursuit algorithm etc., this example uses orthogonal matching pursuit algorithm, tries to achieve initial reconstructed image L _tsparse coefficient α _l, its solution formula is: L _t=D _lα _l.

Step 5. solves the sparse coefficient of high-definition picture

5a) ask initial reconstructed image L respectively _twith n width low resolution training image L _ierror:

$e{r}_i=\frac{{\left|\right|L_t-L_i\left|\right|}_2}{{\left|\right|L_t\left|\right|}_2}$

Obtain initial reconstructed image L _twith the jth width training image L in n width low resolution training image _jleast error: $e{r}_j=\underset{i=1}{\overset{n}{\min }}\left\{e{r}_i\right\};$

5b) judge least error er _jwhether be less than threshold value σ=0.1 of setting,

If error e r _jbe greater than threshold value σ, then return step 1, Resurvey N width part K space data, upgrade dictionary;

If error e r _jbe less than threshold value σ, obtain high-definition picture H to be reconstructed _t' sparse coefficient α _h, solution procedure is:

5b1) ask low resolution difference matrix: Δ l=α _l-α _lj, wherein, α _linitial reconstructed image L _tsparse coefficient, α _ljfor low resolution training image L _jsparse coefficient;

5b2) obtain high resolving power difference matrix Δ h by low resolution difference matrix Δ l:

Make Δ h be an element being full the matrix of zero, matrix size is equal with Δ l, obtains the average a of all non-vanishing elements in Δ l;

5b3) find out high resolving power training image H _jsparse coefficient α _hjthe position that middle all elements is non-vanishing, makes the element in Δ h in same position equal a;

5b4) utilize high resolving power difference matrix Δ h and high resolving power training image H _jsparse coefficient α _hj, try to achieve and treat reconstruct high-definition picture H _t' sparse coefficient: α _h=α _hj+ Δ h.

Step 6. reconstructs high-definition picture

6a) utilize high resolving power dictionary D _hwith high-definition picture H to be reconstructed _t' sparse coefficient α _h, try to achieve high-definition picture: H _t'=D _h* α _h;

6b) to high-definition picture H _t' carry out residual compensation, the method for residual compensation has inverse iteration residual compensation method, and based on the residual compensation method etc. of localized mass, this example is the residual compensation method based on entire image, and detailed process is as follows:

6b1) carry out Fourier transform, 10 times of down-samplings, inversefouriertransforms successively, obtain low-resolution image L _t';

6b2) try to achieve low-resolution image L _t' with initial reconstructed image L _tresidual error: Δ=L _t-L _t';

6b3) according to high-definition picture H _t' and residual delta, try to achieve final reconstructed image H _t=H _t'+Δ.

Effect of the present invention can be further illustrated by following experiment:

1) experiment condition

This experiment adopts 1800 width chest MRI test patterns respectively, and size is 192 × 160, and 1800 width belly MRI test patterns, and size is 192 × 176 as experimental data.

2) experiment content

Utilize zero padding method, PBDW algorithm, TVCMRI algorithm and the present invention respectively, the test pattern of input be reconstructed:

First, 1800 width chest images are reconstructed, wherein, to the reconstruction result of the 100th width chest image as shown in Figure 2, reconstruction result, Fig. 2 (f) reconstruction result of the present invention that the reconstruction result that the reconstruction result that wherein Fig. 2 (a) is the 100th width chest test pattern, Fig. 2 (b) is the 100th width part K space data image of input, Fig. 2 (c) is zero padding method, Fig. 2 (d) are PBDW, Fig. 2 (e) are TVCMRI;

Secondly, 1800 width abdomen images are reconstructed, wherein, to the reconstruction result of the 70th width abdomen images as shown in Figure 3, reconstruction result, Fig. 3 (f) reconstruction result of the present invention that the reconstruction result that the reconstruction result that wherein Fig. 3 (a) is the 70th width belly test pattern, Fig. 3 (b) is the 70th width part K space data image of input, Fig. 3 (c) is zero padding method, Fig. 3 (d) are PBDW, Fig. 3 (e) are TVCMRI;

3) interpretation

As can be seen from Figures 2 and 3, the present invention is better than other method in the visual effect of reconstructed image, it is relatively good that the detailed information of image all keeps, and for the input picture at each position as chest image, abdomen images can obtain good quality reconstruction.

Claims (7)

1., based on a part K spatial sequence image reconstructing method for self-adaptation doubledictionary study, comprise the steps:

(1) N width part K space data is gathered, with the K space data Q that this N width part K space data synthesis n is complete _i, i=1,2 ..., n; To Q _icarry out 10 times of down-samplings, obtain corresponding part K space data P _i; To Q _imake Fourier inversion, obtain high-definition picture H _i, to P _imake Fourier inversion, obtain low-resolution image L _i, this n to high-definition picture H _iwith low-resolution image L _ias training image;

(2) high resolving power training image H is inputted respectively _iwith low resolution training image L _i, and adopt nonoverlapping mode to get the fritter of 4 × 4 to every width training image, obtain initial high resolution dictionary H and initial low resolution dictionary L;

(3) utilize KSVD algorithm to train initial high resolution dictionary H and initial low resolution dictionary L, obtain new high resolving power dictionary D _hwith new low-resolution dictionary D _l, and high-definition picture H _isparse coefficient α _hiwith low-resolution image L _isparse coefficient α _li;

(4) part K space data P to be reconstructed is inputted _t, to this part K space data P _tadopt the process of zero padding method, obtain initial reconstructed image L _t,

(5) low-resolution dictionary D is utilized _lwith initial reconstructed image L _t, solve initial reconstructed image L _tsparse coefficient α _l;

(6) initial reconstructed image L is asked respectively _twith n width low resolution training image L _ierror: obtain initial reconstructed image L _twith the jth width training image L in n width low resolution training image _jleast error:

(7) least error er is judged _jwhether be less than threshold value σ=0.1 of setting, if error e r _jbe less than threshold value σ, obtain high-definition picture H to be reconstructed _t' sparse coefficient α _h; If error e r _jbe greater than threshold value, return step (1), Resurvey N width part K space data, upgrade dictionary;

(8) high resolving power dictionary D is utilized _hwith high-definition picture H to be reconstructed _t' sparse coefficient α _h, try to achieve high-definition picture: H _t'=D _h* α _h; Again to changing high-definition picture H _t' carry out residual compensation, obtain final reconstructed image H _t.

2. the part K spatial sequence image reconstructing method based on the study of self-adaptation doubledictionary according to claim 1, the K space data Q complete with N width part K space data synthesis n gathered wherein described in step (1) _i, i=1,2 ..., n, carries out as follows:

2a) synthesize a width K space data by N/n width part K space data, with the first width part K space data of this N/n width part K space data for standard;

2b) collecting in the second width part K space, and the data that the first width part K space does not collect are added to the first width part K spatially;

2c) collecting in the 3rd width part K space, and the data that front two width part K spaces all do not collect are added to the first width part K spatially, by that analogy, synthesize n complete K space data Q _i.

3. the part K spatial sequence image reconstructing method based on the study of self-adaptation doubledictionary according to claim 1, initial high resolution dictionary H and initial low resolution dictionary L is trained wherein described in step (3), carry out as follows:

3a) to the optimization formula of KSVD algorithm: restrictive condition: be out of shape, be about to wherein be expressed as:

Wherein, Y is the initial dictionary of input, and D is target training dictionary, and X is Its Sparse Decomposition matrix, for any l row, || X _l|| ₀for X _l0 norm, for solving 2 norms of Y-DX, T ₀for degree of rarefication control coefrficient; d _mfor the m row atom of D, for the m of X is capable, K is total columns of D, d _kfor the kth row atom of target training dictionary D, for the row k of X, E _kfor not using the kth row atom d of D _kcarry out the error matrix that signal Its Sparse Decomposition produces;

3b) to the formula after distortion be multiplied by matrix Ω _k, obtain goal decomposition formula

Wherein Ω _ksize be P*| ω _k|, P is the columns of the initial dictionary Y of input, | ω _k| be ω _kmodulus value, and Ω _kat (ω _k(m), m) place is 1, other place is 0, wherein 1≤m entirely≤| ω _k|, ω _km () is ω _km number;

3c) to goal decomposition formula in error matrix carry out decomposition of singular matrix to obtain wherein U is left singular matrix, V ^tfor right singular matrix, Φ is singular value matrix;

3d) get k=1 successively, 2 ..., K, with the kth row atom of the first row of left singular matrix U more fresh target train word allusion quotation D, tries to achieve the dictionary D ' after renewal, obtains new high resolving power dictionary D _hwith new low-resolution dictionary D _l;

3e) utilize the initial dictionary Y of input and the dictionary D ' after upgrading, try to achieve Its Sparse Decomposition matrix X ', obtain high-definition picture H _isparse coefficient α _hiwith low-resolution image L _isparse coefficient α _li.

4. according to claim 1 based on self-adaptation doubledictionary study part K spatial sequence image reconstructing method, wherein described in step (4) to part K space data P _tadopt the process of zero padding method, be to after the data padding do not gathered, then carry out Fourier inversion, obtain reconstructed image.

5. the part K spatial sequence image reconstructing method based on the study of self-adaptation doubledictionary according to claim 1, utilizes low-resolution dictionary D wherein described in step (5) _lwith initial reconstructed image L _t, solve initial reconstructed image L _tsparse coefficient α _l, its solution formula is: L _t=D _lα _l.

6. the part K spatial sequence image reconstructing method based on the study of self-adaptation doubledictionary according to claim 1, obtains high-definition picture H to be reconstructed wherein described in step (7) _t' sparse coefficient α _h, carry out as follows:

6a) ask low resolution difference matrix: Δ l=α _l-α _lj, wherein, α _linitial reconstructed image L _tsparse coefficient, α _ljfor low resolution training image L _jsparse coefficient;

6b) obtain high resolving power difference matrix Δ h by low resolution difference matrix Δ l:

Make Δ h be an element being full the matrix of zero, matrix size is equal with Δ l, obtains the average a of all non-vanishing elements in Δ l;

Find out high resolving power training image H _jsparse coefficient α _hjthe position that middle all elements is non-vanishing, makes the element in Δ h in same position equal a;

6c) utilize high resolving power difference matrix Δ h and high resolving power training image H _jsparse coefficient α _hj, try to achieve and treat reconstruct high-definition picture H _t' sparse coefficient: α _h=α _hj+ Δ h.

7. according to claim 1 based on self-adaptation doubledictionary study part K spatial sequence image reconstructing method, wherein described in step (8) to high-definition picture H _t' carry out residual compensation, carry out as follows:

7a) to high-definition picture H _t' carry out Fourier transform, 10 times of down-samplings, inversefouriertransforms successively, obtain low-resolution image L _t';

7b) try to achieve low-resolution image L _t' with initial reconstructed image L _tresidual error: Δ=L _t-L _t';

7c) try to achieve final reconstructed image H _t=H _t'+Δ.