CN101051388A - Magnetic resonant part K data image reestablishing method based on compound two dimension singular sprectrum analysis - Google Patents

Abstract

The present invention relates to a magnetic resonance partial K data image reconstruction method based on complex two-dimensional singular spectral analysis. Said method includes the following steps: collecting partial K data from preset image space range in magnetic resonance imaging scanner, making zeroization imaging treatment of said partial K data, according to the approximate image obtained by zeroization imaging treatment and partial K data making model parameter estimation, and utilizing mathematical model of magnetic resonance image of complex coefficient weighted two-dimensional singular function and complex two-dimensional singular spectral analysis model making reconstruction of magnetic resonance complex image.

Description

Magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis

Technical field

The present invention relates to medical imaging detection technique field, particularly quick mr imaging technique field specifically is meant a kind of magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis.

Background technology

Along with the continuous development of modern medicine technology, Magnetic resonance imaging (MRI) technology has become means indispensable in the medical imaging detection range.Magnetic resonance signal space (original data space) is called the K space, is the Fourier transform space, and K spatial sampling signal promptly obtains nuclear magnetic resonance (MR) image through delivery again after the Fourier inversion.The method that realizes fast imaging mainly contains raising hardware device class (main field is strong such as improving), adopt tactful fast (as FLASH, EPI or the like), Partial K spacescan (as half spectrum scanning, non-centrosymmetry scanning or the like) and on-right angle network scanning methods such as (as spiral scans).Wherein the imaging of one dimension Partial K spatial data can be under the constant prerequisite of hardware and scan mode, improve sweep velocity exponentially and (see also document: P.Margosian, F.Schmitt, D.Purdy, " Faster MR Imaging:Imaging with Halfthe Data; " Health Care Instrum., vol.1, pp.195-197,1986., with J.van Cuppen and A.van Est, " Reducing MR imaging time by one-sided reconstruction, " Mag.Reso.Imag., vol.5, pp.526-527,1987.).Popular strategy is based on the Partial K spatial data image reconstruction of phase correction at present, typical method has half spectrum POCS, phase correction conjugation balanced method and HM method or the like (see also document: E.M.Haacke, E.D.Lindskog, W.Lin, " A fast; iterative partial Fourier technique capable of local phase recovery; " J.Magn.Reson., vol.92, pp.126-145,1991., V.A.Stenger, D.C.Noll, F.E.Boada, " Partial k-space reconstruction for 3Dgradient echo functional MRI:A comparison of phase correction methods; " Magn.Reson.Med., vol.40, pp.481-490,1998., D.C.Noll, D.G.Nishimura, A.Macovski, " Homodyne detection inmagnetic resonance imaging; " IEEE Trans.Med.Imag., vol.10, pp.154-163,1991., and G.McGibney, M.R.Smith, S.T.Nicholas and A.Crawley, " Quantitative evaluation of several partial-Fourierreconstruction algorithms used in MRI; " Magn.Reson.Med., vol.30, pp.51-59,1993.).

Meanwhile, two-dimentional Partial K spatial data scanning can be saved half sweep time than the scanning of one dimension Partial K spatial data, but its formation method has only zero padding method (FZI method), promptly after the K spatial data zero padding of not gathering, uses the imaging of discrete fourier inverse transformation again.Because the imaging of zero padding method always has pseudo-shadow, this is owing to this method is to utilize sacrifice picture quality to exchange image taking speed for.The distortion phenomenon that these shortcomings caused is enough to make the clinical diagnosis doctor to produce mistaken diagnosis, so that hamper them all the time and enter the gate that clinical medicine is used, bring very big inconvenience so just for people's work and life, and limited further developing of medical imaging detection technique to a certain extent.

Therefore, the key issue that improves Partial K spatial data image taking speed is to seek a kind of new graphical representation method, and promptly a kind of new iconic model under this model, can come presentation video with less variable.The present inventor (specifically sees also document: Luo JH once in the literature, Zhu YM, MR image reconstruction from truncated k-space usinga layer singular point extraction technique, IEEE TRANSACTIONS ON NUCLEAR SCIENCE 51 (1): 157-169Part 1 FEB 2004) point out that any one real digital signal can be represented with the weighted sum of singular function.Therefore, such model K spatial data of being fit to the one dimension Partial K spatial data of real signal and having only a direction in two-dimentional K space is by the image reconstruction problem under the truncation condition.

Summary of the invention

The objective of the invention is to have overcome above-mentioned shortcoming of the prior art, provide a kind of can improve the speed of rebuilding the magnetic resonance complex pattern, effectively reduce image error, improve magnetic resonance image (MRI) quality, highly effective, stable and reliable working performance, the scope of application comparatively widely based on the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis.

In order to realize above-mentioned purpose, the magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis of the present invention is as follows:

Should be based on the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis, its principal feature is that described method may further comprise the steps:

(1) collecting part K data G (k in the image space scope of from magnetic resonance imagine scanner, presetting _x, k _y);

(2) these Partial K data are carried out the zero padding imaging processing;

(3) approximate image and the Partial K data that obtain according to the zero padding imaging are carried out model parameter estimation;

(4), utilize the mathematical model and the multiple two-dimentional singular spectrum analysis model of the magnetic resonance image (MRI) of complex coefficient weighting two dimension singular function to carry out the reconstruct of magnetic resonance complex pattern according to the result of model parameter estimation.

This magnetic resonance image (MRI) mathematical model based on the complex coefficient weighting two dimension singular function of the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis is:

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

Wherein, (x, y), 0≤x, y＜N are that pixel is the complex pattern signal of the two-dimentional magnetic resonance image (MRI) of N * N to g, (x _i, y _i) be g (x, i singular point y), a _iBe the multiple singular value on this i singular point, Q is g (x, the number of singular point y), u _{δ x}(x-x _i, y-y _i) be with (x _i, y _i) be the two-dimentional singular function of singular point.

This multiple two-dimentional singular spectrum analysis model based on the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis is:

$G{(k}_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}{(k}_x,k_y),0\&le;k_x,k_y<N;$

Wherein, G (k _x, k _y)=DFT[g (x, y)], $U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ DFT[] be two-dimensional discrete fourier transform operator.

This zero padding imaging processing of carrying out based on the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis may further comprise the steps:

(1) K space Ω has been divided into data subspace Ω _kWith no datat subspace Ω _k

(2) with Ω _kData in the subspace replace Ω with zero _kThe space remains unchanged, and obtains G (k respectively according to following formula _x, k _y) and U _{δ x}, (k _x, k _y) frequency spectrum data after the zero padding imaging

With

$\tilde{G}(k_x,k_y)=\left\{\begin{array}{cc}G(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(3) calculate G according to following formula _Δ(k _x, k _y) and U _{Δ δ x, i}(k _x, k _y):

$U_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[{\mathrm{\&Delta;u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

$G_{\mathrm{\&Delta;}}(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),0\&le;k_x,k_y<N;$

Wherein, Δ u _{δ x}(x-x _i, y-y _i) be g (x, two-dimentional singular function u y) _{δ x}(x-x _i, y-y _i) in the difference of y direction;

(4) obtain G respectively according to following formula _Δ(k _x, k _y) and U _{Δ δ x, i}(k _x, k _y) frequency spectrum data after the zero padding imaging

With

${\tilde{G}}_{\mathrm{\&Delta;}}(k_x,k_y)=\left\{\begin{array}{cc}{G}_{\mathrm{\&Delta;}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,{\tilde{U}}_{\mathrm{\&Delta;}{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(5) obtain respectively according to following formula

With

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

$\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{\mathrm{\&Delta;}{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

Wherein, IDFT[] be two-dimensional discrete Fuli leaf inverse transformation operator;

(6) according to following formula obtain respectively g (x, y) and g (x is y) at the difference delta g of y direction (x, the complex pattern signal after zero padding imaging y)

With

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

$\mathrm{\&Delta;}\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N.$

This model parameter estimation of carrying out based on the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis may further comprise the steps:

(1) basis

With

Use two-dimentional chromatography to carry out g (x, the estimation processing of unusual point set y);

(2) (x, unusual point set y) carry out the estimation of corresponding singular value to be handled according to g;

(3) (x, unusual point set y) and corresponding singular value are returned as the result of model parameter estimation with resulting g.

Should based on use two dimension chromatography of the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis carry out g (x, the estimation of unusual point set y) is handled and be may further comprise the steps:

(1 ') exists

In, find

The absolute value maximum, promptly

Position coordinates (x _i, y _i), and with (x _i, y _i) add among the singular point formation Q;

Calculate according to following formula (2 ')

$\mathrm{\&Delta;}\tilde{g}(x,y)=\mathrm{\&Delta;}\tilde{g}(x,y)-\mathrm{\&alpha;\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

Wherein $\mathrm{\&alpha;}=\frac{\left|\mathrm{\&Omega;}\right|\mathrm{\&Delta;}\tilde{g}(x_i,y_i)}{\left|{\mathrm{\&Omega;}}_k\right|},$ | Ω | be the number of element among the Ω of space, | Ω _k| be space Ω _kThe number of middle element;

(3 ') are judged $\underset{0\&le;x,y<N}{\max }\left(\right|\mathrm{\&Delta;}\tilde{g}(x,y)\left|\right)>T$ Whether set up, if, then return above-mentioned steps (1 '), wherein T is the threshold values relevant with noise default in the system;

(4 ') on the contrary then with this singular point formation Q as unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) output.

Should based on the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis carry out g (x, the estimation of the singular value of unusual point set y) is handled and can be may further comprise the steps:

(1) calculates according to following formula

$\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right];$

(2) according to unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) choose singular function u _{δ x}(x-x _i, y-y _i), i=1～Q wherein, and calculate according to following formula:

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

(3) according to following formula connection row singular function equation:

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

(4) solve described singular function equation with the pseudo inverse matrix method, obtain a least error and separate, obtain Q plural singular value { a ₁, a ₂..., a _QAnd output.

Should based on the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis carry out g (x, the estimation of the singular value of unusual point set y) is handled and also can be may further comprise the steps:

(1) according to unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) choose singular function u _{δ x}(x-x _i, y-y _i), i=1～Q wherein;

(2) calculate U according to following formula _{δ x, i}(k _x, k _y):

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ I=1～Q wherein;

(3) according to following formula connection row singular spectrum equation:

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),k_x,k_y\&Element;{\mathrm{\&Omega;}}_k;$

(4) solve described singular spectrum equation with the pseudo inverse matrix method, obtain a least error and separate, obtain Q plural singular value { a ₁, a ₂..., a _QAnd output.

Should based on the reconstruct of carrying out the magnetic resonance complex pattern of the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis can for:

The { (x of singular point as a result based on model parameter estimation ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) and singular value { a ₁, a ₂..., a _Q, according to the described complex pattern signal of following formula reconstruct g (x, y):

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

Perhaps also can may further comprise the steps:

(1) based on the { (x of singular point as a result of model parameter estimation ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) and singular value { a ₁, a ₂..., a _Q, according to the described Fourier spectrum data of following formula reconstruct G (k _x, k _y):

$G\left(k\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{W}_{b_q}\left(k\right),k=\mathrm{0,1},...,N-1;$

(2) according to following formula obtain described complex pattern signal g (x, y):

g(x，y)＝IDFT[G(k _x，k _y)]，0≤x，y＜N。

Adopted the magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis of this invention, owing to the image space scope collecting part K spatial data that at first basis is preset from actual magnetic resonance equipment, then these Partial K data are carried out the zero padding imaging processing, then approximate image and the Partial K data message that obtains according to this process zero padding imaging carries out model parameter estimation, last result according to model parameter estimation, utilize the mathematical model and the multiple two-dimentional singular spectrum analysis model of the magnetic resonance image (MRI) of complex coefficient weighting two dimension singular function to carry out the reconstruct of magnetic resonance complex pattern, thereby saved sweep time greatly with respect to one dimension magnetic resonant part K spatial data image reconstruction process, realize fast imaging, guaranteed the high s/n ratio of image simultaneously, high resolving power and pinpoint accuracy; And the two-dimentional Partial K spatial data image reconstructing method of the prior art of comparing, can overcome existing pseudo-shadow in the imaging of zero padding method, effectively reduce image error, accurately show former magnetic resonance image (MRI), provide high-quality reliable image information for the medical nmr imaging detects; Simultaneously, method highly effective of the present invention, stable and reliable working performance, the scope of application are comparatively extensive, bring great convenience for people's work and life, and have also established solid theories and practical basis for further developing with popularization and application on a large scale of medical imaging detection technique.

Description of drawings

Fig. 1 is the singular function image of singular point of the present invention for (58,36).

Fig. 2 a, 2b, 2c be respectively two-dimentional complex pattern signal g in the object simulation magnetic resonance imaging test (x, y), K data image G (k _x, k _y) and g (x is y) at the difference delta g of y direction (x, image y).

Fig. 2 d, 2e, 2f are respectively and use multiple two-dimentional singular spectrum analysis (2DSSA) method of the present invention that Fig. 2 a, 2b, 2c are carried out after the zero padding imaging

With Image.

Fig. 3 a, 3b are respectively and use multiple two-dimentional singular spectrum analysis (2DSSA) method of the present invention to g (x, two-dimentional singular function u y) _{δ x}(x-x _i, y-y _i) carry out singular function after the zero padding imaging of truncated spectrum

Real part functions and imaginary part functional image.

Fig. 3 c, 3d are respectively and use multiple two-dimentional singular spectrum analysis (2DSSA) method of the present invention to g (x, two-dimentional singular function u y) _{δ x}(x-x _i, y-y _i) at the difference delta u of y direction _{δ x}(x-x _i, y-y _i) carry out after the zero padding imaging of truncated spectrum

Real part functions and imaginary part functional image.

Fig. 4 is the course of work principle schematic based on the magnetic resonant part K image reconstruction data method of answering two-dimentional singular spectrum analysis of the present invention.

Fig. 5 is for being modulated the phase diagram of the slow complex image that changes of phase place that forms by emulating image in the emulation experiment.

Fig. 6 a, 6b, 6c be respectively original image in the emulation experiment, use the image after FZI method of the prior art is reconstructed the image of Fig. 6 a and use image after multiple two-dimentional singular spectrum analysis (2DSSA) method of the present invention is rebuild the image of Fig. 6 a.

Fig. 7 is the central point changing condition synoptic diagram of the standard deviation STD between the image of the image that uses multiple two-dimentional singular spectrum analysis (2DSSA) method of the present invention and ZFI method of the prior art in the actual human body magnetic resonance imaging test and respectively the Partial K spatial data is reconstructed and complete K spatial spectrum data reconstruction with the Partial K data.

Fig. 8 is the changing condition synoptic diagram of the standard deviation STD between the image of complete K spatial spectrum data reconstruction in the image that uses multiple two-dimentional singular spectrum analysis (2DSSA) method of the present invention and ZFI method of the prior art in the actual human body magnetic resonance imaging test and respectively section middle part, three-dimensional K space sub image data is reconstructed and this corresponding section with different sections.

Fig. 9 a is for using the image of complete K spatial spectrum data reconstruction in the 88th width of cloth transverse section of three-dimensional K space in the actual human body magnetic resonance imaging test.

Fig. 9 b is a Partial K spatial data synoptic diagram in the 88th width of cloth transverse section of three-dimensional K space in this actual human body magnetic resonance imaging test.

Fig. 9 c is for using the image after FZI method of the prior art is reconstructed according to the data among Fig. 9 b.

Fig. 9 d is for using the image after multiple two-dimentional singular spectrum analysis (2DSSA) method of the present invention is reconstructed according to the data among Fig. 9 b.

Embodiment

In order more to be expressly understood technology contents of the present invention, describe in detail especially exemplified by following examples.

The present invention is collecting part K data in the image space scope of presetting from actual magnetic resonance equipment, and adopt the Partial K spatial data formation method of two-dimentional right angle grid configuration, thereby claim the multiple direct reconstruct magnetic resonance image (MRI) of the two-dimentional unusual spectral analysis method method of utilization to be multiple two-dimentional unusual spectral analysis method (2DSSA, 2 Dimension Singular Spectrum Analysis).

Before setting forth overall work process of the present invention and principle of work,, at first need to provide to give a definition for clearer and more definite its art-recognized meanings:

Definition 1: a given real or multiple digital signal, the non-vanishing point of its difference is a singular point, and the difference value on the singular point is a singular value, and singular value can be that real number also can be a plural number.

Definition 2: real digital signal w (x), x=0,1 ..., N-1 has a unique singular point, and singular value is 1, claims that then w (x) is a singular function.

In order to adapt to the image reconstruction of two-dimentional Partial K spatial data, we define two-dimentional singular function is u _{δ x}(x-x _i, y-y _i):

$u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{\&delta;}(x-x_i)u(y-y_i)---\left(1\right)$

δ (x-x wherein _i) and u (y-y _i) be respectively unit pulse and unit-step function:

$\mathrm{\&delta;}(x-x_i)=\left\{\begin{array}{cc}1,& x={\mathrm{<figure-callout\; id="0"\; label="x\; i"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">x</figure-callout>}}_i\\ 0,& x\&NotEqual;x_i\end{array}\right.,u(y-y_i)=\left\{\begin{array}{cc}1,& y\&GreaterEqual;y_i\\ 0,& y<y_i\end{array}\right..$

As singular point is that the singular function of (58,36) sees also shown in Figure 1.

For the random two-dimensional pixel be N * N magnetic resonance image (MRI) g (x, y), 0≤x, y＜N, its multiple two-dimentional singular function magnetic resonance image (MRI) model is:

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N---\left(2\right)$

That is, (x y) can use u to g _{δ x}(x-x _i, y-y _i) complex weighted and the expression.Wherein, (x _i, y _i) and a _iBe called the singular value of i singular point (to adapt to the situation that magnetic resonance image (MRI) mostly is plural number) and i singular point, Q is unusual counting.Formula (2) both sides are got the y difference, then image g (x, y direction difference delta g y) (x y) is expressed as:

$\mathrm{\&Delta;g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i\mathrm{\&Delta;}{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)$

$=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i\mathrm{\&delta;}(x-x_i)\mathrm{\&delta;}(y-y_i),0\&le;x,y<N$ ......(3)

By following formula as seen, multiple two-dimentional singular function magnetic resonance image (MRI) pattern die shape parameter singular point is the non-vanishing position of y direction difference, and singular value is exactly a difference value.

Discrete fourier transform is got together in formula (2) both sides, must magnetic resonance K space G (k _x, k _y) multiple two-dimentional singular spectrum U _{δ x, i}(k _x, k _y) the analysis image model:

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),0\&le;k_x,k_y<N---\left(4\right)$

G (k wherein _x, k _y)=DFT[g (x, y)], $U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right].$ Be that magnetic resonance K spatial data can be represented by the complex coefficient weighted sum of the spectral function of two-dimentional singular function: DFT[wherein] expression two-dimensional discrete fourier transform operator.In like manner, formula (3) both sides are had with leaf transformation in the fetching:

$G_{\mathrm{\&Delta;}}(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),0\&le;k_x,k_y<N---\left(5\right)$

Wherein, G _Δ(k _x, k _y)=DFT[Δ g (x, y)], $U_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[{\mathrm{\&Delta;u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right].$

Also can define the two-dimentional singular function u of x direction similarly _{δ y}(x-x _i, y-y _i), all theoretical methods of being discussed among the present invention can roughly the same be applied to u _{δ y}(x-x _i, y-y _i) mathematical model of weighted sum presentation video, no longer explanation later on.

If an only known Partial K spatial data, the singular point and the singular value parameter conventional method that obtain model are difficult.Strategy of the present invention is earlier the Partial K spatial data to be carried out the imaging of zero padding method, finds out singular point with chromatography from this image then, obtains singular value with the singular spectrum analysis method at last.

For this reason, the present invention needs earlier K space Ω to be divided into data subspace Ω _kWith no datat subspace Ω _kIn the imaging of zero padding method, with Ω _kSpatial data replaces Ω with zero _kThe space remains unchanged.G (k like this _x, k _y), U _{δ x, i}(k _x, k _y), G _Δ(k _x, k _y) and U _{Δ δ x, i}(k _x, k _y) adopt the frequency spectrum data after the imaging of zero padding method to be expressed as:

$\tilde{G}(k_x,k_y)=\left\{\begin{array}{cc}G(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.---\left(6\right)$

${\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.---\left(7\right)$

${\tilde{G}}_{\mathrm{\&Delta;}}(k_x,k_y)=\left\{\begin{array}{cc}{G}_{\mathrm{\&Delta;}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.---\left(8\right)$

${\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.---\left(9\right)$

Correspondingly, can be with the signal indication of reconstruct after the imaging of Partial K spatial data zero padding method:

$\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right]---\left(10\right)$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right]---\left(11\right)$

$\mathrm{\&Delta;}\tilde{g}(x,y)=\mathrm{IDFT}\left[{\tilde{G}}_{\mathrm{\&Delta;}}(k_x,k_y)\right]---\left(12\right)$

$\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)\right]---\left(13\right)$

IDFT[wherein] expression two-dimensional discrete Fuli leaf inverse transformation operator.To formula (4) both sides truncated spectrum simultaneously, and with the imaging of zero padding method, the image function of then zero padding method

Also can use the zero padding method singular function of truncated spectrum Weighted sum represent:

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N---\left(14\right)$

Equally for formula (5), also can become for:

$\mathrm{\&Delta;}\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N---\left(15\right)$

See also shown in Fig. 2 a～2f, wherein clearly provided G (k _x, k _y) with

G (x, y) with

And Δ g (x, y) with Comparison.

See also again shown in Fig. 3 a～3d, wherein clearly provided Ω _k: 0≤k _x, k _yUnder＜72 the situation

Real part functions and the imaginary part function and Real part functions and the image of imaginary part function.

Next need examination

See also shown in Fig. 3 c, the functional value that maximum is arranged on former singular point position, and on two directions Gibbs' effect is arranged in length and breadth.And

On Gibbs' effect be by all

The coefficient result of Gibbs' effect.So,

It can be singular point that maximum place has most.In view of the above, it is as follows that two-dimentional chromatography is asked for the algorithm of singular point:

The first step,

In find

The absolute value maximum, promptly

Position coordinates (x _i, y _i) (being singular point), and add among the singular point formation Q;

In second step, calculate according to following formula:

$\mathrm{\&Delta;}\tilde{g}(x,y)=\mathrm{\&Delta;}\tilde{g}(x,y)-\mathrm{\&alpha;\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N---\left(16\right)$

Wherein $\mathrm{\&alpha;}=\frac{\left|\mathrm{\&Omega;}\right|\mathrm{\&Delta;}\tilde{g}(x_i,y_i)}{\left|{\mathrm{\&Omega;}}_k\right|},$ | Ω | and | Ω _k| difference representation space Ω and Ω _kThe number of middle element;

In the 3rd step, judge $\underset{}{\max }$ $\underset{0\&le;x,y<N}{}\left(\right|\mathrm{\&Delta;}\tilde{g}(x,y)\left|\right)>T$ Whether set up, wherein T is the given threshold values relevant with noise default in the system, if set up, then returns the first step and continues to look for singular point; Otherwise with this singular point formation Q as unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) output.

Through after the above-mentioned step of asking for singular point, the method for asking for the pairing singular value of this unusual point set has following two kinds:

(1) in image space, (x, unusual point set y) have just determined (x, the singular function after zero padding imaging y) about g according to above-mentioned g

Collection, the singular value of unusual point set can be determined according to formula (14) or (15).Its algorithm is as follows:

The first step is asked for $\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right];$

Second step is according to unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) choose singular function, i=1～Q, and calculate according to following formula:

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

The 3rd step, connection row singular function equation:

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N$

Solve described singular function equation with the pseudo inverse matrix method, obtain a least error and separate, determine singular value { a ₁, a ₂..., a _Q.

(2) in the K space, (x, unusual point set y) have just determined (x, singular spectrum function U y) about g according to above-mentioned g _{δ x, i}(k _x, k _y) collection, utilize formula (15) can calculate the singular value that the singular point set pair is answered, its algorithm is as follows:

The first step is by unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) choose singular function u _{δ x}(x-x _i, y-y _i);

In second step, i=1～Q calculates:

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

The 3rd step, connection row singular spectrum equation:

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),k_x,k_y\&Element;{\mathrm{\&Omega;}}_k.---\left(17\right)$

Solve described singular spectrum equation with the pseudo inverse matrix method, obtain a least error and separate, determine singular value { a ₁, a ₂..., a _Q.

Be to carry out image reconstruction at last, following two kinds of methods arranged according to above-mentioned singular point and singular value information:

(1) according to singular point { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) and singular value { a ₁, a ₂..., a _Q, by formula (2) reconstructed image g (x, y):

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N$

(2) by formula the K spatial data G (k that (4) reconstruct is not gathered _x, k _y):

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),0\&le;k_x,k_y<N;$

And then with the imaging of discrete fourier inverse transformation:

g(x，y)＝IDFT[G(k _x，k _y)]，0≤x，y＜N。

See also shown in Figure 4ly, should may further comprise the steps based on the magnetic resonant part K image reconstruction data method of multiple two-dimentional singular spectrum analysis:

(1) collecting part K data G (k in the image space scope of from magnetic resonance imagine scanner, presetting _x, k _y);

(2) these Partial K data are carried out the zero padding imaging processing, may further comprise the steps:

(a) K space Ω has been divided into data subspace Ω _kWith no datat subspace Ω _k

(b) with Ω _kData in the subspace replace Ω with zero _kThe space remains unchanged, and obtains G (k respectively according to following formula _x, k _y) and U _{δ x, i}(k _x, k _y) frequency spectrum data after the zero padding imaging With

$\tilde{G}(k_x,k_y)=\left\{\begin{array}{cc}G(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(c) calculate G according to following formula _Δ(k _x, k _y) and U _{Δ δ x, i}(k _x, k _y):

$U_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[{\mathrm{\&Delta;u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

$G_{\mathrm{\&Delta;}}(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),0\&le;k_x,k_y<N;$

Wherein, Δ u _{δ x}(x-x _i, y-y _i) be g (x, two-dimentional singular function u y) _{δ x}(x-x _i, y-y _i) in the difference of y direction;

(d) obtain G respectively according to following formula _Δ(k _x, k _y) and U _{Δ δ x, i}(k _x, k _y) frequency spectrum data after the zero padding imaging

With

${\tilde{G}}_{\mathrm{\&Delta;}}(k_x,k_y)=\left\{\begin{array}{cc}{G}_{\mathrm{\&Delta;}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,{\tilde{U}}_{\mathrm{\&Delta;}{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{\mathrm{\&Delta;}{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_{\mathrm{<figure-callout\; id="0"\; label="k"\; filenames="A20071004065500231.png,A20071004065500241.png"\; state="\{\{state\}\}">k</figure-callout>}}\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(e) obtain respectively according to following formula

With

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

$\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{\mathrm{\&Delta;}{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

Wherein, IDFT[] be two-dimensional discrete Fuli leaf inverse transformation operator;

(f) according to following formula obtain respectively g (x, y) and g (x is y) at the difference delta g of y direction (x, the complex pattern signal after zero padding imaging y)

With

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

$\mathrm{\&Delta;}\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

(3) the Partial K data message according to this process zero padding imaging processing carries out model parameter estimation, may further comprise the steps:

(a) basis

With Use two-dimentional chromatography carry out g (x, the estimation of unusual point set y) is handled, and may further comprise the steps:

(i) exist

In, find

The absolute value maximum, promptly

Position coordinates (x _i, y _i), and with (x _i, y _i) add among the singular point formation Q;

(ii) calculate according to following formula

$\mathrm{\&Delta;}\tilde{g}(x,y)=\mathrm{\&Delta;}\tilde{g}(x,y)-\mathrm{\&alpha;\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

Wherein $\mathrm{\&alpha;}=\frac{\left|\mathrm{\&Omega;}\right|\mathrm{\&Delta;}\tilde{g}(x_i,y_i)}{\left|{\mathrm{\&Omega;}}_k\right|},$ | Ω | be the number of element among the Ω of space, | Ω _k| be space Ω _kThe number of middle element;

(iii) judge $\underset{}{\max }$ $\underset{0\&le;x,y<N}{}\left(\right|\mathrm{\&Delta;}\tilde{g}(x,y)\left|\right)>T$ Whether set up, if, then return above-mentioned steps (i), wherein T is the threshold values relevant with noise default in the system;

Otherwise (iv) then with this singular point formation Q as unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) output.

(b) according to g (x, unusual point set y) carry out the estimation of corresponding singular value to be handled, and can may further comprise the steps:

(i) calculate according to following formula

$\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right];$

(ii) according to unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) choose singular function u _{δ x}(x-x _i, y-y _i), i=1～Q wherein, and calculate according to following formula:

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

(iii) according to following formula connection row singular function equation:

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

(iv) solve described singular function equation, obtain a least error and separate, obtain Q plural singular value { a with the pseudo inverse matrix method ₁, a ₂..., a _QAnd output;

Also can may further comprise the steps:

(i) according to unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) choose singular function u _{δ x}(x-x _i, y-y _i), i=1～Q wherein;

(ii) calculate U according to following formula _{δ x, j}(k _x, k _y):

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ I=1～Q wherein;

(iii) according to following formula connection row singular spectrum equation:

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),k_x,k_y\&Element;{\mathrm{\&Omega;}}_k;$

(iv) solve described singular spectrum equation, obtain a least error and separate, obtain Q plural singular value { a with the pseudo inverse matrix method ₁, a ₂..., a _QAnd output;

(c) with resulting g (x, unusual point set y) and corresponding singular value are returned as the result of model parameter estimation;

(4), utilize the mathematical model and the multiple two-dimentional singular spectrum analysis model of the magnetic resonance image (MRI) of complex coefficient weighting two dimension singular function to carry out the reconstruct of magnetic resonance complex pattern according to the result of model parameter estimation; Wherein, the magnetic resonance image (MRI) mathematical model of this complex coefficient weighting two dimension singular function is:

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

Wherein, (x, y), 0≤x, y＜N are that pixel is the complex pattern signal of the two-dimentional magnetic resonance image (MRI) of N * N to g, (x _i, y _i) be g (x, i singular point y), a _iBe the multiple singular value on this i singular point, Q is g (x, the number of singular point y), u _{δ x}(x-x _i, y-y _i) be with (x _i, y _i) be the two-dimentional singular function of singular point;

Should multiple two-dimentional singular spectrum analysis model be:

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),0\&le;k_x,k_y<N;$

Wherein, G (k _x, k _y)=DFT[g (x, y)], $U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ DFT[] be two-dimensional discrete fourier transform operator;

The reconstruct of this magnetic resonance complex pattern can for:

The { (x of singular point as a result based on model parameter estimation ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) and singular value { a ₁, a ₂..., a _Q, according to the described complex pattern signal of following formula reconstruct g (x, y):

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

Perhaps also can may further comprise the steps:

(a) based on the { (x of singular point as a result of model parameter estimation ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) and singular value { a ₁, a ₂..., a _Q, according to the described Fourier spectrum data of following formula reconstruct G (k _x, k _y):

$G\left(k\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{W}_{b_q}\left(k\right),k=\mathrm{0,1},...,N-1;$

(b) according to following formula obtain described complex pattern signal g (x, y):

g(x，y)＝IDFT[G(k _x，k _y)]，0≤x，y＜N。

The emulation experiment that below is to use method of the present invention to carry out:

At first the appliance computer emulated data is carried out test of heuristics, and the image that emulation experiment is used is that a width of cloth tonal range is 0～255, and picture size is 128 * 128.Considering magnetic resonance image (MRI) in most cases, is slow variation in the phase place of image-region, and variation range generally (if the phase change frequency surpasses this scope, can be proofreaied and correct by K space center's point shift method) within [0 °, 360 °]; So in the experiment, getting the phase change scope is 0 °～360 °.

Emulation experiment one: the susceptibility experiment of noise

With the complex image that becomes a width of cloth phase place slowly to change in the emulating image modulation, its phase diagram sees also shown in Figure 5, and the standard deviation that adds 0 average again is respectively 1,2,3,4,5,6,7,8,9 white Gaussian noise, as original image.Then this original image is carried out Fourier transform, get K spatial data Ω again _k={ 24＜k _x, k _y＜40} carries out image reconstruction with method of the present invention, calculates the standard deviation of reconstructed image and original image at last.Experimental result is as shown in table 1.

STD	1	2	3	4	5	6	7	8	9
										2DSSA	1.3582	2.3576	3.6088	4.2717	4.3534	5.5481	6.0089	7.2448	7.6032
FZI	9.9308	9.9842	9.9928	10.051	10.142	10.148	10.187	10.315	10.332

Table 1. compares with the reconstruction accuracy that noise changes

Data result from table 1 as can be seen, the standard deviation (STD) of 2DSSA method reconstructed image of the present invention is approached the STD that adds noise, noise is very little to 2DSSA algorithm affects of the present invention, and for strong noise image the effect that improves signal to noise ratio (S/N ratio) is arranged.This mainly be since 2DSSA by chromatography avoiding introducing false singular point because of noise, by finding the solution the false singular point of sneaking in the further elimination chromatography of pseudoinverse equation, thereby suppressed the influence of noise greatly to 2DSSA.2DSSA contains the high high frequency K space of noise component because of clipping, and then high frequency component signal has been reduced in reconstruct again, thereby has improved signal to noise ratio (S/N ratio).And in the FZI of prior art method, in the low noise section, its STD mainly comes from and blocks artefact influence, and in the double influence that the strong noise section is subjected to blocking artefact and noise, so higher STD is all arranged, reconstruction quality is relatively poor relatively.

Emulation experiment two: complex image phase change susceptibility experiment

The adding standard deviation is that 5 average is 0 Gaussian noise on above-mentioned emulating image, emulating image behind the plus noise is carried out the original image that phase modulation (PM) constitutes test usefulness, the phase change scope of modulation is respectively 40 °, 80 °, 120 °, 160 °, 200 °, 240 °, 280 °, 320 °, 360 °, then this original image is carried out Fourier transform, get K spatial data Ω again _k={ 24＜k _x, k _y＜40} carries out image reconstruction with method of the present invention, calculates the standard deviation of reconstructed image and original image at last.Experimental result is as shown in table 2.

PHASE	40	80	120	160	200	240	280	320	360
										2DSSA	4.6015	4.5875	4.5502	4.7380	5.1505	4.7065	4.115	4.2908	4.3534
FZI	10.157	10.025	10.167	10.604	10.125	10.213	10.213	10.186	10.142

Table 2. compares with the reconstruction accuracy that image phase changes

By the data of table 2 as can be known, without any correlativity, the ZFI of 2DSSA method just of the present invention and prior art is not subjected to the influence of the slow phase change of complex image to phase change to the 2DSSA reconstruction accuracy slowly.Because slowly the phase change that changes produces the image difference size, with respect to the differential variation that noise produces, want little much slight.Find the solution in chromatography, pseudoinverse under the effect of the false singular point of compacting of singular value method, phase change is difficult for producing the influence to the 2DSSA reconstruction accuracy.But actual phase change is not that the SPA sudden phase anomalies situation happens occasionally as so simple in the above-mentioned experiment.

See also shown in Figure 6ly, wherein provide the example of an emulating image reconstruct.Ω _k={ 24＜k _x, k _y＜40}, the STD of its additivity Gauss zero mean noise is 5.In the figure of FZI method, to block artefact and be seen everywhere, the 2DSSA method is not almost blocked artefact, more only has the details of some areas to become fuzzy slightly with original image.

Below be to use method of the present invention to carry out the test experiments of actual MRI data.

Experiment three: the space parallax opposite sex of the Partial K data acquisition of analytical algorithm

The actual magnetic resonance image of experiment usefulness is that a width of cloth tonal range is 0～255, and picture size is 256 * 176.The Partial K space size that is used to test is 128 * 88.The experimentation design is as follows:

Allow the central point of Partial K data space, from initial point, with horizontal direction step pitch 4, vertically step pitch 6 moves to the lower right corner, and whenever moves the ZFI method reconstructed image of once using 2DSSA method of the present invention and prior art, the image of the image of reconstruct and complete K spatial data reconstruct, compare, and provided the central point situation of change of standard deviation STD, see also shown in Figure 7 with the Partial K data.From Fig. 7, can find:

The first, the ZFI of 2DSSA method comparison prior art of the present invention has much better reconstruction accuracy, the partial frequency spectrum data that the better reconstruct of this explanation 2DSSA energy is not got;

The second, under same Partial K spatial data condition, the center position in collection space is deviation from origin slightly, unusual useful to the extraction of 2DSSA method of the present invention with singular value, but should not be too asymmetric, otherwise spectrum energy is lost too much, and introduce than mistake.

Experiment four: 2DSSA is to the susceptibility of image anatomical structure

Actual magnetic resonance data volume is 256 * 256 * 276 three-dimensional K spatial data, and experimental design is as follows:

According to the result of above-mentioned experiment three, getting Partial K data acquisition scope is Ω _k={ 28＜k _x＜60 ,-42＜k _y＜86}, so use the ZFI method of 2DSSA method of the present invention and prior art to be reconstructed to each picture, at last with reconstructed image with compose the data reconstruction image fully and compare, and provided STD situation of change with different transverse sections, the result as shown in Figure 8, wherein show the same influence that is subjected to the image anatomical structure of ZFI method of 2DSSA method of the present invention and prior art, but the 2DSSA method always has than ZFI higher reconstructed image precision is arranged.

In order to observe this method effect intuitively from vision, see also shown in Fig. 9 a～9d, wherein provided the 88th width of cloth cross-sectional view picture in the image volume, wherein Fig. 9 a is the image that complete K spatial data is done, Fig. 9 b is a Partial K space synoptic diagram, Fig. 9 c is according to the data in the Partial K space among Fig. 9 b image with the reconstruct of ZFI method of the prior art institute, and Fig. 9 d is according to the data in the Partial K space among Fig. 9 b image with the reconstruct of 2DSSA method of the present invention institute.As can be seen directly perceived from Fig. 9 c and Fig. 9 d in the ZFI of prior art method, blocked artefact and is seen everywhere, and in the image of 2DSSA method of the present invention almost difficulty seek this artefact that blocks, we can say that 2DSSA method of the present invention is an efficient ways.

From as can be seen above, the basic thought of method of the present invention is: at first provide two-dimentional singular spectrum analysis image reconstruction model, for the phase problem that solves magnetic resonance image (MRI) is introduced complex weighted coefficient, utilization chromatography, pseudo inverse matrix are determined the singular value method, thereby can suppress noise preferably, and eliminated phase change to restructuring graph picture element amount caused adverse effect.In the test experiments of emulation experiment and actual magnetic resonance image data, 2DSSA method of the present invention all shows the result more much better than prior art ZFI.2DSSA method of the present invention is generalized to the one dimension Partial K spatial data reconstruct problem of studying at present in the two dimension and goes, and this will provide a kind of new thinking methods to the image reconstruction problem that solves Partial K space MR data.

Adopted above-mentioned magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis, owing to the image space scope collecting part K spatial data that at first basis is preset from actual magnetic resonance equipment, then these Partial K data are carried out the zero padding imaging processing, then approximate image and the Partial K data that obtain according to the zero padding imaging are carried out model parameter estimation, last result according to model parameter estimation, utilize the mathematical model and the multiple two-dimentional singular spectrum analysis model of the magnetic resonance image (MRI) of complex coefficient weighting two dimension singular function to carry out the reconstruct of magnetic resonance complex pattern, thereby saved sweep time greatly with respect to one dimension magnetic resonant part K spatial data image reconstruction process, realize fast imaging, guaranteed the high s/n ratio of image simultaneously, high resolving power and pinpoint accuracy; And the two-dimentional Partial K spatial data image reconstructing method of the prior art of comparing, can overcome existing pseudo-shadow in the imaging of zero padding method, effectively reduce image error, accurately show former magnetic resonance image (MRI), provide high-quality reliable image information for the medical nmr imaging detects; Simultaneously, method highly effective of the present invention, stable and reliable working performance, the scope of application are comparatively extensive, bring great convenience for people's work and life, and have also established solid theories and practical basis for further developing with popularization and application on a large scale of medical imaging detection technique.

In this instructions, the present invention is described with reference to its certain embodiments.But, still can make various modifications and conversion obviously and not deviate from the spirit and scope of the present invention.Therefore, instructions and accompanying drawing are regarded in an illustrative, rather than a restrictive.

Claims (9)

1, a kind of magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis is characterized in that described method may further comprise the steps:

(1) collecting part K data G (k in the image space scope of from magnetic resonance imagine scanner, presetting _x, k _y);

(2) these Partial K data are carried out the zero padding imaging processing;

(3) approximate image and the Partial K data that obtain according to the zero padding imaging are carried out model parameter estimation;

(4), utilize the mathematical model and the multiple two-dimentional singular spectrum analysis model of the magnetic resonance image (MRI) of complex coefficient weighting two dimension singular function to carry out the reconstruct of magnetic resonance complex pattern according to the result of model parameter estimation.

2, the magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis according to claim 1 is characterized in that, the magnetic resonance image (MRI) mathematical model of described complex coefficient weighting two dimension singular function is:

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

Wherein, (x, y), 0≤x, y＜N are that pixel is the complex pattern signal of the two-dimentional magnetic resonance image (MRI) of N * N to g, (x _i, y _i) be g (x, i singular point y), a _iBe the multiple singular value on this i singular point, Q is g (x, the number of singular point y), u _{δ x}(x-x _i, y-y _i) be with (x _i, y _i) be the two-dimentional singular function of singular point.

3, the magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis according to claim 2 is characterized in that, described multiple two-dimentional singular spectrum analysis model is:

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),0{\&le;k}_x,k_y<N;$

Wherein, G (k _x, k _y)=DFT[g (x, y)], $U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ DFT[] be two-dimensional discrete fourier transform operator.

4, the magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis according to claim 3 is characterized in that the described zero padding imaging processing of carrying out may further comprise the steps:

(1) K space Ω has been divided into data subspace Ω _kWith no datat subspace Ω _k

(2) with Ω _kData in the subspace replace Ω with zero _kThe space remains unchanged, and obtains G (k respectively according to following formula _x, k _y) and U _{δ xi}(k _x, k _y) frequency spectrum data after the zero padding imaging

With

$\tilde{G}(k_x,k_y)=\left\{\begin{array}{cc}G(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,$ ${\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(3) calculate G according to following formula _Δ(k _x, k _y) and U _{Δ δ x, j}(k _x, k _y):

$U_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[{\mathrm{\&Delta;u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

$G_{\mathrm{\&Delta;}}(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),0\&le;k_x,k_y<N;$

Wherein, Δ u _{δ x}(x-x _i, y-y _i) be g (x, two-dimentional singular function u y) _{δ x}(x-x _i, y-y _i) in the difference of y direction;

(4) obtain G respectively according to following formula _Δ(k _x, k _y) and U _{Δ δ x, j}(k _x, k _y) frequency spectrum data after the zero padding imaging

With

${\tilde{G}}_{\mathrm{\&Delta;}}(k_x,k_y)=\left\{\begin{array}{cc}{G}_{\mathrm{\&Delta;}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,$ ${\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(5) obtain respectively according to following formula

With

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

$\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)\right];$

Wherein, IDFT[] be two-dimensional discrete Fuli leaf inverse transformation operator;

(6) according to following formula obtain respectively g (x, y) and g (x is y) at the difference delta g of y direction (x, the complex pattern signal after zero padding imaging y)

With

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

$\mathrm{\&Delta;}\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N.$

5, the magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis according to claim 4 is characterized in that the described model parameter estimation of carrying out may further comprise the steps:

(1) basis With

Use two-dimentional chromatography to carry out g (x, the estimation processing of unusual point set y);

(2) (x, unusual point set y) are constructed two-dimentional singular spectrum equation, carry out the estimation of corresponding singular value and handle according to g;

(3) (x, unusual point set y) and corresponding singular value are returned as the result of model parameter estimation with resulting g.

6, the magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis according to claim 5 is characterized in that, the two-dimentional chromatography of described use carry out g (x, the estimation of unusual point set y) is handled and be may further comprise the steps:

(1 ') exists In, find

The absolute value maximum, promptly

Position coordinates (x _i, y _i), and with (x _i, y _i) add among the singular point formation Q;

Calculate according to following formula (2 ')

$\mathrm{\&Delta;}\tilde{g}(x,y)=\mathrm{\&Delta;}\tilde{g}(x,y)-\mathrm{\&alpha;\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

Wherein $\mathrm{\&alpha;}=\frac{\left|\mathrm{\&Omega;}\right|\mathrm{\&Delta;}\tilde{g}(x_i,y_i)}{\left|{\mathrm{\&Omega;}}_k\right|},$ | Ω | be the number of element among the Ω of space, | Ω _k| be space Ω _kThe number of middle element;

(3 ') are judged $\underset{0\&le;x,y<N}{\max }\left(\right|\mathrm{\&Delta;}\tilde{g}(x,y)\left|\right)>T$ Whether set up, if, then return above-mentioned steps (1 '), wherein T is the threshold values relevant with noise default in the system;

(4 ') on the contrary then with this singular point formation Q as unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) output.

7, the magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis according to claim 6 is characterized in that, described carry out g (x, the estimation of the singular value of unusual point set y) is handled and be may further comprise the steps:

(1) calculates according to following formula

$\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right];$

(2) according to unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) choose singular function u _{δ x}(x-x _i, y-y _i), i=1～Q wherein, and calculate according to following formula:

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

(3) according to following formula connection row singular function equation:

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

(4) solve described singular function equation with the pseudo inverse matrix method, obtain a least error and separate, obtain Q plural singular value { a ₁, a ₂..., a _QAnd output.

8, the magnetic resonant part K image reconstruction data method based on multiple two-dimentional singular spectrum analysis according to claim 6 is characterized in that, described carry out g (x, the estimation of the singular value of unusual point set y) is handled and be may further comprise the steps:

(1) according to unusual point set { (x ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) choose singular function u _{δ x}(x-x _i, y-y _i), i=1～Q wherein;

(2) calculate U according to following formula _{δ x, j}(k _x, k _y):

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ I=1～Q wherein;

(3) according to following formula connection row singular spectrum equation:

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),k_x,k_y\&Element;{\mathrm{\&Omega;}}_k;$

(4) solve described singular spectrum equation with the pseudo inverse matrix method, obtain a least error and separate, obtain Q plural singular value { a ₁, a ₂..., a _QAnd output.

9, according to claim 7 or 8 described magnetic resonant part K image reconstruction data methods, it is characterized in that, describedly carry out being reconstructed into of magnetic resonance complex pattern based on multiple two-dimentional singular spectrum analysis:

The { (x of singular point as a result based on model parameter estimation ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) and singular value { a ₁, a ₂..., a _Q, according to the described complex pattern signal of following formula reconstruct g (x, y):

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),0\&le;x,y<N;$

Perhaps may further comprise the steps:

(1) based on the { (x of singular point as a result of model parameter estimation ₁, y ₁), (x ₂, y ₂) ..., (x _Q, y _Q) and singular value { a ₁, a ₂..., a _Q, according to the described Fourier spectrum data of following formula reconstruct G (k _x, k _y):

$G\left(k\right)=\underset{q=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_q{W}_{b_q}\left(k\right),k=\mathrm{0,1},...,N-1;$

(2) according to following formula obtain described complex pattern signal g (x, y):

g(x，y)＝IDFT[G(k _x，k _y)]，0≤x，y＜N。

Images