The method and apparatus of parallel collection image reconstruction
Technical field
The present invention relates to magnetic resonance imaging (MRI, Magnetic Resonance Imaging) technology, particularly relate to the method and apparatus of a kind of parallel acquisition (PA, Parallel Acquisition) image reconstruction.
Background technology
In the MRI technology, the speed of imaging is unusual important parameters.Early stage inspection usually needs to spend several hours, and owing to the improvement of technology aspect field intensity, gradient hardware and pulse train, the speed of imaging has had a bigger lifting subsequently.But, continuous radio frequency (the RF of field gradient Fast transforms and high density, Radio Frequency) the pulse meeting brings the unaffordable specific absorption rate (SAR of the Human Physiology limit, Specific Absorption Rate) and the pyrogenicity amount of organ-tissue, therefore, the lifting of image taking speed has run into bottleneck.
Subsequently, the researchist finds that by using complicated computerized image reconstruction algorithm and the coil array that matches, the image taking speed of magnetic resonance can be promoted greatly, and this kind technology is commonly called parallel imaging (parallel imaging) technology.Parallel collection image reconstruction is a kind of image reconstruction (reconstruction) technology that is used for parallel acquisition, it utilizes spatial sensitivity (sensitivity) difference of phased-array coil to carry out space encoding, and gather simultaneously with phased-array coil, obtain doubly even higher image taking speed than the fast 2-6 of conventional magnetic resonance imaging.Adopt parallel imaging technique, magnetic resonance imaging system is proposed new requirement, as a plurality of receiving cables of needs, complex array coil and coil sensitivity calibration, with special data processing and image rebuilding method or the like.
Parallel imaging can improve image acquisition speed, and the raising of image taking speed realizes by reducing the K space filling factor.If being lower than the restriction of sampling thheorem, the K space filling factor can cause directly utilizing the image after Fourier is rebuild pseudo-shadow to occur.General magnetic resonance imaging is to obtain image by the frequency domain information of gathering object by Fourier transform.According to sampling thheorem, be inversely proportional to the sampling interval of frequency domain in the repetition period of image area object.If the space repetition period of image, the image after the reconstruction will overlap less than the size of image itself, this phenomenon is called aliasing in signal Processing.If equal the inverse of phase-encoding direction FOV in the sampling interval of phase-encoding direction, image just just can not overlap, claim this moment the K space be sampled as full sampling.The sampling of sampling interval during greater than 1/FOV is called owes sampling, that is to say that the K spatial information that collects is not enough to rebuild complete image, owes to sample to cause the overlapping of image; Otherwise be over-sampling, also can not cause the doubling of the image.The synoptic diagram that Fig. 1 shows full sampling and owes to sample, wherein left figure is full sampling, right figure owes sampling.In Fig. 1, vertical direction is a phase-encoding direction, and left side figure is full sampling, and the sampling interval on phase-encoding direction is Δ Ky; Right figure be 1/2 owe the sampling, the sampling interval on phase-encoding direction is 2 Δ Ky.
In magnetic resonance imaging, the magnetic resonance signal of tissue is by the receiving coil collection.Fig. 2 is the image that phase-encoding direction shown in Figure 1 is completely sampled and obtained when owing to sample.Among Fig. 2, vertical direction is a phase directional, and horizontal direction is a frequency direction.As shown in Figure 2, the image that left figure obtains when being full the sampling, sampling interval is Δ Ky; Right figure is 1/2 image of owing to sample and obtaining, and sampling interval is 2 Δ Ky.As shown at right, if (Coil1 and Coil2) gathers the sample among the figure respectively with two adjacent surface coils, and the sampling rate of phase-encoding direction is 2/FOV, will obtain having two width of cloth images of overlapping pseudo-shadow after directly rebuilding, the visual field of every width of cloth image is FOV/2, but the content that all comprises another width of cloth figure among every width of cloth figure, it is the contribution that the value of each pixel among the right figure comes from two volume elements among the former figure, for instance, the signal intensity of pixel 1 comes from two volume elements 1 among the left figure among the right figure, 2 contribution, equally, the signal intensity of pixel 2 also comes from two volume elements 2 among the left figure among the right figure, 1 contribution.From right figure, be not difficult to observe a little less than the signal intensity of signal intensity than original image of pseudo-shadow, this be since the sample spin signal that coil collects by the sensitivity function weighting of coil, and the amplitude of the sensitivity function of surface coils is far away regional less of distance coil.The aliasing that sampling causes depends on the filling mode in K space, the just track in K space.Pseudo-shadow above-mentioned is the pseudo-shadow that uniform grid is owed to sample and caused, for helical trajectory, the pseudo-shadow that causes of owing to sample will be very irregular.
For multi-thread circle collection, though the K spatial information that each coil collects is not enough, the difference that can utilize different coils to collect signal obtains a complete image through handling.The reconstruction algorithm of the pseudo-shadow of parallel imaging anti-aliasing can be divided into two classes haply: space harmonics parallel acquisition (SMASH, SiMultaneous Acquisition of Spatial Harmonics) and sensitivity encoding parallel acquisition technique (SENSE, SENSitivity Endoding parallel acquisitiontechnique).Wherein, the SMASH method is a kind of method of utilizing the sensitivity function composition space harmonics of each passage coil and carrying out auxiliaring coding.The sensitivity function of general coil is slowly to change in the space, and can be similar to and think gauss of distribution function, so can utilize the linear combination of the sensitivity function of each passage coil to constitute the space harmonics of certain frequency.And utilize the space harmonics function to gather out not have the phase encoding of actual acquisition line.
The scheme of handling at frequency domain with SMASH is different, and the SENSE method is eliminated the pseudo-shadow of owing to sample and causing by the method for finding the solution system of linear equations at image area.Because the eclipsing effects that causes of image space cycle, for instance, when sampling rate for full sampling 1/2 the time, directly utilize of the contribution of the value of each pixel in the image that Fourier transform obtains from two volume elements in the original image.Can explain with formula (1):
s
1=p
11m
1+p
12m
2(1)
Formula (1) can cooperate Fig. 2 to illustrate: the value of right figure mid point 1 by the magnetization m1 of volume elements 1 among the left figure multiply by coil 1 the magnetization m2 of the sensitivity P11 at point 1 place and volume elements 2 multiply by coil 1 point 2 places sensitivity P12's and decide.But P11 and P12 are measuring amount, want to obtain m1 and m2, need two equations at least, i.e. the data that need two absolute coils to measure at least.With matrix representation be exactly,
Wherein S2 is the value of right figure mid point 2, and the definition of P21 and P22 is the same.Formula (2) can further be expressed as:
S=PM (3)
Wherein S is the vector of Nc * 1, and Nc is a port number.M is the vector of Np * 1, and Np is for owing sampling rate.P inverted just obtain original image.For port number greater than the situation of owing sampling rate, the Nc>Np in the equation (3), this is an overdetermination system of linear equations, promptly the number of known conditions is greater than the number of unknown number.Generally ask the plus sige of P to reach optimum against making to separate.
Wherein, matrix W is called as weighting matrix, and superscript H represents conjugate transpose.
In radio frequency electromagnetic field, human body is regarded as load, and the sensitivity function of coil can change when measuring different human bodies to some extent, and these variations are enough to influence quality of reconstructed images, so the sensitivity of coil generally need be obtained in real time.Completely sample in the k central zone of space during SENSE method image data, owe sampling in outer peripheral areas.So original K spatial data is divided into two parts: evenly owe the data that the data of sampling and low frequency are completely sampled.Evenly owe the data of the sampling image that generates aliasing, and the data that low frequency is completely sampled are with generating fuzzy tissue image and then obtaining real-time coil sensitivity profiles and weighting matrix, to utilize data aliased image that generates and the weighting matrix that utilizes the full sampled data of low frequency to obtain of evenly owing to sample to synthesize at last, obtain high-resolution no aliased image.Here, the full sampled data of low frequency that is used to obtain coil sensitivity profiles and weighting matrix is called reference data, and the phase encoding line that K space medium and low frequency is completely sampled is called reference line.
The advantage of utilizing hyperchannel coil parallel acquisition to improve image taking speed is conspicuous.If under the identical sampling time, utilize parallel acquisition technique also can mention the resolution of image.Parallel acquisition also can attach some benefits, such as alleviating owing to the image artifacts that goes resonance (off-resonance) to cause.Simultaneously, because the raising of picking rate, parallel acquisition technique can also alleviate motion artifacts.
Yet the realization parallel acquisition need be paid certain cost.At first and since actual acquisition to data reduced, though utilize the parallel acquisition reconstruction technique can reconstruct the image of no aliasing, the signal to noise ratio (S/N ratio) of image is the full signal noise ratio (snr) of image of sampling and obtaining under the same hardware condition
N wherein
fBe speedup factor, therefore, speedup factor is big more, and picture quality reduces big more, and promptly signal to noise ratio (S/N ratio) (SNR, Signal-to-Noise Ratio) reduces big more.
Summary of the invention
One object of the present invention is to propose a kind of method of parallel collection image reconstruction, improves the SNR that parallel acquisition is rebuild the back image.
Another object of the present invention is to propose a kind of device of parallel collection image reconstruction, improve the SNR that parallel acquisition is rebuild the back image.
To achieve these goals, technical scheme of the present invention comprises:
A kind of method of parallel collection image reconstruction, it comprises: evenly owe the data reconstruction that sampled data and the full sampled data of low frequency combine according to mixing K space that sampling pattern generates magnetic resonance MR imaging; Calculate coil sensitivity profiles according to the full sampled data of described low frequency; According to described data reconstruction, described coil sensitivity profiles and described mixing sampling pattern reconstructed image.
Wherein, describedly evenly owe the data reconstruction that sampled data and the full sampled data of low frequency combine and comprise: generate the data reconstruction matrix according to mixing K space that sampling pattern generates magnetic resonance MR imaging, in described data reconstruction matrix, evenly owe sampled data and the full sampled data of low frequency for the K space of actual acquisition, fill with the actual acquisition value, for the data of not gathering, fill with 0.
Wherein, the line number of described data reconstruction matrix is acquisition channel number, phase encoding number and the product of reading the direction sampling number, and columns is 1.
Wherein, described data computation coil sensitivity profiles of completely sampling according to low frequency comprises: calculate the coil sensitivity profiles matrix according to the full sampled data of low frequency, and arrangement and the dimension of determining described sensitivity profile matrix according to the arrangement and the dimension of described data reconstruction matrix.
Wherein, described according to data reconstruction, coil sensitivity profiles and mix the sampling pattern reconstructed image and comprise: the mathematic(al) representation of determining described mixing sampling pattern; Mathematic(al) representation reconstructed image according to described data reconstruction, described coil sensitivity profiles matrix and described mixing sampling pattern.
Wherein, the mathematic(al) representation that described generation mixes sampling pattern comprises: generate and mix sampling matrix, on the diagonal line of described mixing sampling matrix, fill at interval with 0 and 1 according to owing sampling rate, and the relevant position of completely sampling in the centre fills with 1, and other positions except that diagonal line are all with 0 filling.
Wherein, describedly comprise: arrangement and the dimension of determining the inverse Fourier transform matrix according to the arrangement and the dimension of described data reconstruction matrix according to data reconstruction, coil sensitivity profiles matrix and the mathematic(al) representation reconstructed image that mixes sampling pattern; Utilize described inverse Fourier transform matrix, described mixing sampling matrix and coil sensitivity profiles matrix computations to mix the parallel imaging sampling matrix of sampling pattern: A=M
UndersamplingM
FourierM
Sensitivity, wherein, M
FourierExpression inverse Fourier transform matrix, M
SensitivityExpression coil sensitivity profiles matrix, M
UndersamplingExpression mixes sampling matrix, and A represents to mix the parallel imaging sampling matrix of sampling pattern: contrary as reconstruction matrix with the parallel imaging sampling matrix of described mixing sampling pattern; Utilize described reconstruction matrix and described data reconstruction matrix reconstructed image: ima=(A ' A)
-1A ' V
Rawdata, wherein, ima represents reconstructed image matrix, V
RawdataExpression data reconstruction matrix.
Wherein, parallel imaging sampling matrix and data reconstruction matrix reconstructed image that described utilization mixes sampling pattern comprise: with described data reconstruction matrix decomposition is N vector, each vector comprises M element, wherein M is the product of acquisition channel number and phase encoding line number, and N is for reading the direction sampling number; Respectively a described N vector is rebuild, obtained N image vector; With described N image vector combination, obtain the reconstructed image matrix.
A kind of device of parallel collection image reconstruction, it comprises: the data reconstruction generation unit is used for evenly owing the data reconstruction that sampled data and the full sampled data of low frequency combine according to mixing K space that sampling pattern generates magnetic resonance MR imaging; The sensitivity profile computing unit is used for calculating coil sensitivity profiles according to the full sampled data of described low frequency; Image reconstruction unit is used for according to described data reconstruction, described coil sensitivity profiles and described mixing sampling pattern reconstructed image.
Wherein, described data reconstruction generation unit generates the data reconstruction matrix, in described data reconstruction matrix, evenly owes sampled data and the full sampled data of low frequency for the K space of actual acquisition, fills with the actual acquisition value, for the data of not gathering, fills with 0.
Wherein, described sensitivity profile computing unit calculates according to the full sampled data of described low frequency and generates the coil sensitivity profiles matrix, and arrangement and the dimension of determining described coil sensitivity profiles matrix according to the arrangement and the dimension of described data reconstruction matrix.
Wherein, described image reconstruction unit comprises: mix the sampling matrix generation module, be used for generating the mixing sampling matrix; Inverse Fourier transform matrix generation module is used for arrangement and the dimension of determining the inverse Fourier transform matrix according to the arrangement and the dimension of described data reconstruction matrix, and generates described inverse Fourier transform matrix; The reconstruction matrix generation module, be used to utilize described inverse Fourier transform matrix, described sensitivity profile matrix and described mixing sampling matrix to generate the parallel imaging sampling matrix that mixes sampling pattern, and contrary as reconstruction matrix with the parallel imaging sampling matrix of described mixing sampling pattern; The image reconstruction computing module is used to utilize described reconstruction matrix and described data reconstruction matrix reconstructed image.
Wherein, described mixing sampling matrix generation module is when generating the mixing sampling matrix, fill at interval with 0 and 1 according to owing sampling rate on the diagonal line of described mixing sampling matrix, and the relevant position of completely sampling in the centre fills with 1, other positions except that diagonal line are all with 0 filling.
From above technical scheme as can be seen, in the present invention, combine as the parallel acquisition data that are used to rebuild with the full sampled data of low frequency evenly owing sampled data in the K space, because being positioned at K space center partial data comprises more useful information, therefore the data that K space low frequency is completely sampled can comprise more useful information.So, utilize the parallel acquisition data that combine the full sampled data of low frequency to carry out image reconstruction, more useful information can be provided, thereby effectively improve the SNR of reconstructed image.
Description of drawings
To make clearer above-mentioned and other feature and advantage of the present invention of those of ordinary skill in the art by describing the preferred embodiments of the present invention in detail with reference to accompanying drawing below, identical label is represented identical parts, in the accompanying drawing:
Fig. 1 is the synoptic diagram that K phase encode direction is completely sampled and owed to sample in the MRI scanning;
Fig. 2 is the image that phase-encoding direction shown in Figure 1 is completely sampled and obtained when owing to sample;
Fig. 3 is the method flow diagram of parallel collection image reconstruction according to an embodiment of the invention;
Fig. 4 (a) shows and evenly owes the phase encoding of sampling line;
Fig. 4 (b) shows the phase encoding line that low frequency is completely sampled;
Fig. 4 (c) shows the combination of the phase encoding line of uniform sampling and the phase encoding line that low frequency is completely sampled;
Fig. 5 is the synoptic diagram that mixes sampling matrix according to an embodiment of the invention;
Fig. 6 is the device process flow diagram of parallel collection image reconstruction according to an embodiment of the invention;
Fig. 7 be consider under the situation of reference line with the situation of not considering reference line under the SNR comparison diagram of reconstructed image.
Embodiment
In order to make purpose of the present invention, technical scheme and advantage clearer,, the present invention is further elaborated below in conjunction with drawings and Examples.Should be appreciated that specific embodiment described herein only in order to explanation the present invention, and be not used in qualification the present invention.
According to the present invention, expire sampled data and evenly owe the data reconstruction that sampled data combines according to the low frequency that mixes in the sampling pattern generation K space, and when image reconstruction, consider to mix sampling pattern.
Below by a specific embodiment method of the present invention is described in detail.
Fig. 3 is the method flow diagram of parallel collection image reconstruction according to an embodiment of the invention.As shown in Figure 3, in the present embodiment, the method for parallel collection image reconstruction comprises the steps: at least
Step 301: generate the full sampled data of K space low frequency and evenly owe the data reconstruction that sampled data combines according to mixing sampling pattern.
The K spatial data gathered comprise that collection evenly owes the data and the reference data of sampling.Reference data is the data that the low frequency part in K space is completely sampled and obtained.
Fig. 4 (a), Fig. 4 (b) and Fig. 4 (c) will evenly owe sampled data and the low frequency synoptic diagram that full sampled data combines.Wherein Fig. 4 (a) shows and evenly owes the phase encoding of sampling line, and Fig. 4 (b) shows the phase encoding line that low frequency is completely sampled, and Fig. 4 (c) shows the combination of the phase encoding line of uniform sampling and the phase encoding line that low frequency is completely sampled.Suppose that the speedup factor of owing to sample in the K space is 2, promptly sampling rate is 1/2, for instance, if full sampling needs to gather 17 phase encoding lines on the phase-encoding direction, is expressed as k respectively
8To k
-8, sampling rate is that 9 phase encoding lines need are gathered in 1/2 the sampling of evenly owing, and illustrates as solid line among Fig. 4 (a), these 9 phase encoding lines are expressed as [k respectively
8, k
6, k
4, k
2, k
0, k
-2, k
-4, k
-6, k
-8]; Reference data takies 5 phase encoding lines of the most close K space center, and promptly 5 reference lines illustrate as solid line among Fig. 4 (b), are expressed as [k respectively
2, k
1, k
0, k
-1, k
-2]; The data and the reference data combination of uniform sampling are obtained data reconstruction, consequently, union is got in the set of phase encoding line and the reference line set of evenly owing to sample, and the data that will be somebody's turn to do and concentrate are as data reconstruction, illustrate as solid line among Fig. 4 (c), the phase encoding line that should and concentrate is expressed as [k respectively
8, k
6, k
4, k
2, k
1, k
0, k
-1, k
-2, k
-4, k
-6, k
-8].
With the matrix representation of the data reconstruction of the data of evenly owing to sample and reference data be combined into is V
Rawdata, its length calculation formula is:
V
RawdataLength=nChnPEnRO (5)
Wherein, nCh is the acquisition channel number, and nPE is a phase encoding line number, and nRO reads the direction sampling number, at V
RawdataIn the matrix, the data on the phase encoding line of actual acquisition are constant, and the data on the phase encoding line of owing to adopt are 0.Preferably, the line number of data reconstruction matrix is acquisition channel number, phase encoding line number and the product of reading the direction sampling number, and columns is 1.
Step 302: calculate the coil sensitivity profiles matrix M according to the full sampled data of low frequency
Sensitivity
M
SensitivityArrangement and dimension according to V
RawdataArrangement and dimension determine.
Step 303: determine to mix the mathematic(al) representation of sampling pattern, preferably, determine to mix sampling matrix M
Undersampling
Mix sampling matrix M
UndersamplingExample as shown in Figure 5, as can be seen from Figure 5, at M
UndersamplingDiagonal line on, evenly owe the position of sampling and fill at interval with 1 and 0 according to sampling rate, fill with 1 the relevant position of middle full sampling, 0 filling is all used in the position beyond the diagonal line.
Step 304: according to the coil sensitivity profiles matrix M
Sensitivity, the data reconstruction matrix V
RawdataWith mixing sampling matrix M
UndersamplingReconstructed image.
At first need according to V
RawdataThe arrangement of matrix and dimension are determined the inverse Fourier transform matrix M
FourierArrangement and dimension;
Then, generate the parallel imaging sampling matrix that mixes sampling pattern according to formula (6):
A=M
Undersampling·M
Fourier·M
Sensitivity (6)
Wherein, M
FourierExpression inverse Fourier transform matrix, M
SensitivityExpression coil sensitivity profiles matrix, M
UndersamplingExpression mixes sampling matrix, and A represents to mix the parallel imaging sampling matrix of sampling pattern.
Then, will mix parallel imaging sampling matrix A contrary of sampling pattern as reconstruction matrix, and according to formula (7) reconstructed image:
ima=(A′·A)
-1·A′·V
rawdata(7)
Wherein, ima represents the reconstructed image matrix.
Preferably, in actual applications, can be with the data reconstruction matrix V
RawdataBe decomposed into N vector, being about to the data reconstruction matrix decomposition is N vector, and each vector comprises M element, and wherein M is the product of acquisition channel number and phase encoding line number, and N is for reading the direction sampling number; Respectively a described N vector is rebuild, obtained N image vector; Then, described N image vector combination obtained the reconstructed image matrix.
Fig. 6 shows the structural representation of the parallel collection image reconstruction device of one specific embodiment according to the present invention.As shown in Figure 6, in the present embodiment, the parallel collection image reconstruction device comprises data reconstruction generation unit 601, sensitivity profile computing unit 602 and image reconstruction unit 603 at least.
Data reconstruction generation unit 601 evenly owes sampled data according to the K space of mixing sampling pattern generation MR imaging and low frequency is expired the data reconstruction that sampled data combines.
Sensitivity profile computing unit 602 calculates coil sensitivity profiles according to the full sampled data of low frequency.
Image reconstruction unit 603 is according to data reconstruction, coil sensitivity profiles and mixing sampling pattern reconstructed image.
Preferably, data reconstruction generation unit 601 generates the data reconstruction matrix, in the data reconstruction matrix, evenly owes sampled data and the full sampled data of low frequency for the K space of actual acquisition, fills with the actual acquisition value, for the data of not gathering, fills with 0.
Preferably, sensitivity profile computing unit 602 calculates according to the full sampled data of low frequency and generates the coil sensitivity profiles matrix, and arrangement and the dimension of determining the coil sensitivity profiles matrix according to the arrangement and the dimension of data reconstruction matrix.
Further, image reconstruction unit 603 comprises mixing sampling matrix generation module 6031, inverse Fourier transform matrix generation module 6032, reconstruction matrix generation module 6033 and image reconstruction computing module 6043.
Mix sampling matrix generation module 6031 and generate the mixing sampling matrix;
Arrangement and dimension that inverse Fourier transform matrix generation module 6032 is determined the inverse Fourier transform matrix according to the arrangement and the dimension of data reconstruction matrix, and generate the inverse Fourier transform matrix;
Reconstruction matrix generation module 6033 utilizes inverse Fourier transform matrix, sensitivity profile matrix and mixes sampling matrix and generates the parallel imaging sampling matrix that mixes sampling pattern, and will mix sampling pattern the parallel imaging sampling matrix against as reconstruction matrix;
Image reconstruction computing module 6043 utilizes reconstruction matrix and data reconstruction matrix reconstructed image.
Preferably, mix sampling matrix generation module 6031 when generating the mixing sampling matrix, fill at interval with 0 and 1 according to owing sampling rate on the diagonal line of described mixing sampling matrix, and the relevant position of completely sampling in the centre fills with 1, other positions except that diagonal line are all with 0 filling.
As mentioned above, in the present invention, combine with the full sampled data of low frequency evenly owing sampled data in the K space, as the parallel acquisition data that are used to rebuild, because being positioned at K space center partial data comprises more useful information, therefore the data that K space low frequency is completely sampled can comprise more useful information.So the parallel acquisition data that combine the full sampled data of low frequency can provide more useful information, thereby effectively improve the SNR of reconstructed image behind Fourier transform.Be appreciated that the variation with the ratio of the quantity of total phase encoding line changes SNR improvement situation according to reference line, in general, the ratio that reference line accounts for total phase encoding line is big more, and SNR is high more.
Fig. 7 shows the SNR that ratio that reference line accounts for total phase encoding line is a reconstructed image under the SNR of reconstructed image under 10% the situation and the situation of not considering reference line.In Fig. 7, the ratio that the curve representation reference line of top accounts for total phase encoding line is the SNR of reconstructed image under 10% the situation, and the curve representation of below is not considered the SNR of reconstructed image under the situation of reference line.As can be seen from Figure 7, compare with the SNR of reconstructed image under the situation of not considering reference line, the ratio that reference line accounts for total phase encoding line is under 10% the situation, and SNR improves and is about 2%.
The above only is preferred embodiment of the present invention, not in order to restriction the present invention, all any modifications of being done within the spirit and principles in the present invention, is equal to and replaces and improvement etc., all should be included within protection scope of the present invention.