Abstract:
The present invention proposes a fast GeneRalized Autocalibrating Partially Parallel Acquisition (GRAPPA) image reconstruction algorithm for magnetic resonance imaging. The algorithm simplifies data fitting and channel merging in the process of reconstruction into a one-step linear calculation. Parameters needed to perform the linear calculation step can be pre-calculated and stored, thereby greatly increasing the image reconstruction speed and solving the problem of the relatively long image reconstruction time needed by prior art GRAPPA algorithms. Also, the algorithm can employ a weighting matrix to conveniently compare signal-to-noise ratio losses of images brought by different types of reconstruction methods in image domain and frequency domain.

Description:
BACKGROUND OF THE INVENTION  
       [0001]     1. Field of the Invention  
         [0002]     The present invention relates to a GeneRalized Autocalibrating Partially Parallel Acquisition (GRAPPA) image reconstruction algorithm, and more particularly to a fast GRAPPA image reconstruction algorithm for magnetic resonance imaging (MRI).  
         [0003]     2. Description of the Prior Art  
         [0004]     In MRI techniques, imaging speed is a very important parameter. Early examinations often took several hours, and later, as technologies in field intensity, gradient hardware and pulse sequence developed, a relatively large increase in the imaging speed was achieved. However, rapidly switching field gradient and high-density consecutive radio frequency (RF) pulses may result in a specific absorption rate (SAR) and heat amount in organs and tissues that exceed human physiological limits. Therefore, a further increase in the imaging speed encounters a limitation.  
         [0005]     It has been found that the imaging speed of MRI can be greatly increased with a technique that employs complicated computer image reconstruction algorithms and a matching coil array. Such techniques are known as parallel imaging techniques. Parallel imaging techniques include SiMultaneous Acquisition of Spatial Harmonics (SMASH), SENSitivity Encoding Parallel Acquisition Techniques (SENSE), GeneRalized Autocalibrating Partially Parallel Acquisitions (GRAPPA) etc. Parallel acquisition image reconstruction, an image reconstruction technique for fast acquisitions, utilizes differences in the spatial sensitivities of phased array coils to perform spatial encoding while using the phased array coils to perform signal acquisitions, achieving an imaging speed two to six times higher than that of conventional MRI, or even higher.  
         [0006]     Using parallel imaging techniques imposes new requirements on MRI systems. For example, multiple receiving channels and multi-element array coils are needed, coil sensitivity should be calibrated, and special data processing and image reconstruction methods must be adopted.  
         [0007]     SMASH is a parallel acquisition reconstruction method that employs coil sensitivity to fit a spatial harmonic function and fills under-sampled data. The algorithm thereof is characterized by data of all channels being added directly to serve as a basis and subject for multi-channel fitting. There is a relatively large error present in the fitting operation of said algorithm, causing severe artifacts in SMASH images as well as low signal-to-noise ratios (SNRs).  
         [0008]     GRAPPA is an enhanced SMASH sequence, the GRAPPA algorithm utilizes sampled data of all channels to fit and fill under-sampled data of each of the channels, and performs channel merging operation on fitted full-sampled image of each of the channels, i.e., calculating the sum of squares for each image, adding the sum, and obtaining the extraction of square root of the addition to obtain a resulting image. The GRAPPA algorithm reduces the fitting error and improves the image quality. GRAPPA, however, has a relatively long image reconstruction time since the time needed to fit full channel data is proportional to the number of channels. As the number of channels for parallel acquisitions in MRI devices is increased gradually, this disadvantage of GRAPPA in terms of the image reconstruction speed will become more and more apparent.  
         [0009]     Therefore, an urgent problem in the field of MRI is how to provide a fast GRAPPA image reconstruction algorithm that will greatly increase the calculation speed of the GRAPPA image reconstruction and calculate SNR losses of images obtained by parallel acquisition image reconstruction algorithms in the image domain and the frequency domain.  
       SUMMARY OF THE INVENTION  
       [0010]     An main object of the present invention is to provide a fast GRAPPA image reconstruction algorithm for MRI that greatly increases the calculation speed of GRAPPA image reconstruction.  
         [0011]     Another object of the present invention is to provide a fast GRAPPA image reconstruction algorithm for MRI that calculates SNR losses of images obtained by parallel acquisition image reconstruction algorithms in the image domain and the frequency domain.  
         [0012]     These objects are achieved in accordance with the present invention by a fast GRAPPA image reconstruction algorithm for MRI, including the following steps:  
         [0000]     (i) expressing the fitting calculation of a channel in the frequency domain as a convolution of channels, which is transformed into the image domain;  
         [0000]     (ii) converting the channel merging calculation in the image domain into a linear calculation of images of individual channels and the sensitivities of the coils; and  
         [0000]     (iii) merging the linear calculations in steps (i) and (ii) as a weighting matrix.  
         [0013]     In step (i), the formula for expressing fitting calculation of a channel as convolution calculation of channels is  
           rawdata   kCH   full     =       ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       rawdata   iCh   ppa     ⊗       F     iCh   ,   kCh       .           ,       
 
 where rawdata kCh   full  represents the fully-sampled k-space data of the kChth channel fitted after the convolution, rawdata iCh   ppa  represents the under-sampled k-space data of the iChth channel used to fit, and F iCh,kCh  is the convolution kernel derived based on the fitting parameter. 
 
         [0014]     By performing a Fourier transform on F iCh,kCh , the convolution is transformed into a point multiplication calculation in image domain as follows:  
           image   kCh   full     =       ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       image   iCh   ppa     ·     f     iCh   ,   kCh             ,       
 
 where image kCh   full  is the fully-sampled the image of the kChth channel, image iCh   ppa  is the under-sampled image of the iChth channel, and f iCh,kCh  is the Fourier transform of F iCh,kCh . 
 
         [0015]     In step (ii), the formula for the merging calculation by a sum of squares is  
           image   SOS     =         ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       image   iCh     ·     image   iCh   *             ;       
 
 and when the sensitivity function of the iChth channel is S iCh , the merging calculation by a sum of squares can also be expressed as  
           image   SOS     =            image        ⁢         ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       s   iCh     ·     s   iCh   *             =          image        ·     I   SOS           ,       
 
 where I SOS  is the image obtained after performing the merging calculation by a sum of squares on the sensitivities of the coils and this image is achieved by multiplying the original image by I SOS  of the coils. 
 
         [0016]     In step (ii), the merging calculation by a sum of squares can also employ a merging method based on multiplication by complex conjugate of the sensitivities of the coils (SCS), and the formula thereof is  
           image   SOS     =         ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       image   iCh     ·     s   iCH   SOS         =          image        ·     I   SOS           ,     where   ⁢           ⁢     image   SCS           
 
 is the image obtained by the merging method employing SCS, s iCh   SOS  is the sensitivities of the coils educed based on I SOS , and I SOS  is the image obtained after performing the merging calculation by a sum of squares on the sensitivities of the coils; and the formula for calculating s iCh   SOS  is  
         s   iCh   SOS     -       image   iCh     /           ∑   iCh     ⁢       image   iCh     ·     image   iCh   *           .           
 
         [0017]     The calculations in steps (i) and (ii), for being linear, can be expressed as a consecutive operation of two sparse matrixes. Therefore, the merging calculation of the two steps in step (iii) is equivalent to multiplying the two sparse matrixes, and thus middle steps thereof are omitted, achieving the effect of accelerating the calculation. The weighting matrix by multiplying the two sparse matrixes is  
         w   kCh     =       ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       f     iCh   ,   kCh       ·       s   iCh     SOS   *       .             
 
         [0018]     In the fast GRAPPA image reconstruction algorithm for MRI according to the present invention, the accelerated calculation can be performed in the frequency domain, for performing the accelerated calculation in the frequency domain, after the calculation in step (i), a convolution of data of individual channels is performed and the Fourier transformation of the sensitivities of the coils is done once again, and a further Fourier transformation is performed to obtain the merged image. The formula for performing the accelerated calculation in the frequency domain is  
           rawdata   SCSfull     =       ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       rawdata   iCh   full     ⊗     S   iCh     SOS   *             ,       
 
 where s iCh   SOS  and s iCh   SOS  are a Fourier transform pair. 
 
         [0019]     The fast GRAPPA image reconstruction algorithm for MRI according to the present invention can also use the weighting matrix to conveniently compare SNR losses of images brought by different types of reconstruction methods in image domain and frequency domain based on the formula as follows:  
         SNR   rel     =           ∑     k   =       1   ⁢   …   ⁢           ⁢   m     =   nCh         ⁢       w   kCh     ·     s   kCh   SOS               ∑     k   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢              s   kCh   SOS          2     ·         ∑     k   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢            w   kCh          2                 .         
 
         [0020]     In summary, the fast GRAPPA image reconstruction algorithm for MRI according to the present invention can simplify data fitting and channel merging in the process of reconstruction into one-step linear calculation, and parameters needed to perform the linear calculation step can be pre-calculated and stored, thereby greatly increasing the image reconstruction speed. Also, the inventive method can use the weighting matrix to conveniently compare SNR losses of images brought by different types of reconstruction methods in the image domain and the frequency domain, thereby solving the problem of the relatively long image reconstruction time needed by conventional GRAPPA algorithms. 
     
    
     DESCRIPTION OF THE DRAWING  
       [0021]     The single FIGURE is a schematic diagram for fitting fully-sampled data with under-sampled data by the fast GRAPPA image reconstruction algorithm for MRI according to the present invention. 
     
    
     DESCRIPTION OF THE PREFERRED EMBODIMENTS  
       [0022]     In functional imaging applications, there is a need for collecting multiple images of the same tissue in the same position at different points in time for several repetitions under the same imaging parameters. Under such conditions, the sensitivity of the coils is constant. Therefore, the fast GRAPPA image reconstruction algorithm for MRI according to the present invention simplifies data fitting and channel merging in the process of reconstruction into a one-step linear calculation. Parameters needed to perform the linear calculation step can be pre-calculated and stored, thereby greatly increasing the image reconstruction speed.  
         [0023]     The fast GRAPPA image reconstruction algorithm for MRI according to the present invention includes the initial step of expressing the fitting calculation of a channel in the frequency domain as a convolution of channels and then transforming it into the image domain.  
         [0024]     Referring to the FIGURE, the principle of expressing the fitting calculation of GRAPPA in the frequency domain as a convolution and transforming the convolution into the image domain is illustrated by an example, in which under-sampled data of four channels for four coils are used to fit fully-sampled data of the fourth channel, but the present invention is not limited to this example. In the FIGURE, white points represent under-sampled data that have not been filled, points with lines leading thereto represent data that are being filled by the fitting method, and black points represent actually sampled data. Since the data fitting is a linear calculation, the fitting calculation of the fourth channel can be converted into a convolution of the four channels:  
               rawdata   kCh   full     =       ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       rawdata   iCh   ppa     ⊗     F     iCh   ,   kCh                   (   1   )             
 
 where rawdata kCh   full  represents the fully-sampled k-space data of the kChth channel fitted by the convolution. In the FIGURE kCh is 4. rawdata iCh   ppa  represents the under-sampled k-space data of the iChth channel used to fit. F iCh,kCh  is the convolution kernel educed based on the fitting parameter. By performing Fourier transform on F iCh,kCh , the convolution calculation can also be transformed into point multiplication calculation in image domain as follows:  
               image   kCh   full     =       ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       image   iCh   ppa     ·     f     iCh   ,   kCh                   (   2   )             
 
 where image iCh   full  is the fully-sampled image of the kChth channel, image iCh   ppa  is the under-sampled image of the iChth channel, and f iCh,kCh  is the Fourier transform of F iCh,kCh . 
 
         [0025]     The second step of the method is converting the channel merging calculation in the image domain into a linear calculation of images of individual channels and sensitivities of the coils.  
         [0026]     The formula for the merging calculation by a sum of squares is:  
               image   SOS     =          image        ⁢         ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       image   iCh     ·     image   iCh   *                     (   3   )             
 
         [0027]     If the sensitivity of the iChth channel is s iCh , the merging calculation via sum of squares can also be expressed as:  
               image   SOS     =           ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       S   iCh     ·     S   iCh   *           =          image        ·     I   SOS                 (   4   )             
 
 where I SOS  is the image obtained after performing the merging calculation by a sum of squares on the sensitivities of the coils. It is not difficult to prove that the image obtained by the merging calculation by a sum of squares is achieved by multiplying the original image by I SOS  of coils. 
 
         [0028]     The formula of a merging method based on multiplication by the complex conjugate of the sensitivities of the coils (SCS) is.  
               image   SCS     =         ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       image   iCh     ·     S   iCh     SOS   *           =          image        ·     I   SOS                 (   5   )             
 
 where image SCS  is the image obtained by the merging method employing SCS, and s iCh   SOS  is the sensitivity of the coils educed based on I SOS . The formula is as follows:  
               S   iCh   SOS     =       image   iCh           ∑   iCh     ⁢       image   iCh     ·     image   iCh   *                     (   6   )             
 
         [0029]     The equivalence of the two merging methods is easy to prove.  
         [0030]     The third step in the inventive method is to merge the linear calculations in steps (i) and (ii) as a weighting matrix.  
         [0031]     The calculations in steps (i) and (ii), being linear, can be expressed as a consecutive operation of two sparse matrixes. Merging of the two steps is equivalent to multiplying the two sparse matrixes, and thus middle steps thereof are omitted, thereby accelerating the calculation.  
         [0032]     The weighting matrix is expressed as follows:  
               w   kCh     =       ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       f     iCh   ,   kCh       ·     S   iCh     SOS   *                   (   7   )             
 
         [0033]     The accelerated calculation can be performed in the frequency domain.  
         [0034]     The three calculation steps described above can be performed completely in the frequency domain. For this purpose, after the calculation of formula (1), the convolution is performed on data of individual channels with the Fourier transform of the sensitivities of the coils being done once again, and further performing a Fourier transformation to obtain the merged image:  
               rawdata   SCSfull     =       ∑     iCh   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       rawdata   iCh   full     ⊗     S   iCh     SOS   *                   (   8   )             
 
 where s iCH   SOS  and S iCh   SOS  are a Fourier transform pair. Performing the accelerated calculation in the image domain or doing so in the frequency domain are substantially equivalent and the main difference therebetween lies in the truncation error present when using the convolution S iCh   SOS  in the frequency domain. 
 
         [0035]     The weighting matrix can be used to calculate SNR losses of images.  
         [0036]     Using the weighting matrix conveniently allows a comparison SNR losses of images obtained by different types of reconstruction methods in the image domain and the frequency domain. Before the inventive method, the problem of SNR loss has not been addressed in parallel reconstruction methods in image domain. The formula is as follows:  
               SNR   rel     =         ∑     k   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢       w   kCh     ·     S   kCh   SOS                 ∑     k   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢            S   kCh   SOS          2         ·         ∑     k   =     1   ⁢   …   ⁢           ⁢   nCh         ⁢            w   kCh          2                     (   9   )             
 
         [0037]     In summary, in an embodiment of the fast GRAPPA image reconstruction algorithm for MRI according to the present invention, assuming the number of channels is eight and the image resolution for full sampling is 256×256, the amount of calculation needed to reconstruct an image using the fast GRAPPA image reconstruction algorithm according to the present invention, by performing said steps (i) to (v), is only that for eight times of two-dimensional Fourier transform, 256×256×8 times of multiplication of complex numbers and 256×256×7 times of addition of complex numbers. In addition, compared to conventional GRAPPA, the fast GRAPPA image reconstruction algorithm according to the present invention can shorten the reconstruction time by more than an order of magnitude since inverse calculation is not needed.  
         [0038]     Although modifications and changes may be suggested by those skilled in the art, it is the intention of the inventors to embody within the patent warranted hereon all changes and modifications as reasonably and properly come within the scope of their contribution to the art.