# Patent

## Title
Fast generalized autocalibrating partially parallel acquisition image reconstruction algorithm for magnetic resonance imaging

## 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.

## Background
BACKGROUND OF THE INVENTION
1. Field of the Invention
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).
2. Description of the Prior Art
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 duty cycle 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.
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 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.
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.
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).
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 combination 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.
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
An main object of the present invention is to provide a method for MRI that greatly increases the calculation speed of GRAPPA image reconstruction.
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.
These objects are achieved in accordance with the present invention by a fast GRAPPA image reconstruction algorithm for MRI, including the following steps:
(i) expressing the fitting calculation of a channel in the frequency domain as convolutions of channels, which are transformed into the image domain;
(ii) converting the channel combination calculation in the image domain into a linear calculation of images of individual channels and the sensitivities of the coils; and
(iii) merging the linear calculations in steps (i) and (ii) as a weighting matrix.
In step (i), the formula for expressing fitting calculation of a channel as convolution calculations of channels is




  
    
      rawdata
      kCH
      full
    
    =
    
      
        ∑
        
          iCh
          =
          
            1
            ⁢
            ...
            ⁢
            
              
            
            ⁢
            nCh
          
        
      
      ⁢
      
        
          rawdata
          iCh
          ppa
        
        ⊗
        
          
            F
            
              iCh
              ,
              kCh
            
          
          .
        
      
    
  
  ,




where rawdatakChfull represents the fully-sampled k-space data of the kChth kChth channel fitted after the convolution, rawdataiChppa represents the under-sampled k-space data of the iChth channel used to fit, and FiCh,kCh are the convolution kernels derived based on the fitting parameter.

By performing a Fourier transform on FiCh,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 imagekChfull is the fully-sampled image of the kChth channel, imageiChppa is the under-sampled image of the iChth channel, and fiCh,kCh is the Fourier transform of FiCh,kCh.

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 profile of the iChth channel is SiCh, the combination calculation by a sum of squares (SOS) can also be expressed as





  
    
      image
      SOS
    
    =
    
      
        
          
          image
          
        
        ⁢
        
          
            
              ∑
              
                iCh
                =
                
                  1
                  ⁢
                  ...
                  ⁢
                  
                    
                  
                  ⁢
                  nCh
                
              
            
            ⁢
            
              
                s
                iCh
              
              ·
              
                s
                iCh
                *
              
            
          
        
      
      =
      
        
          
          image
          
        
        ·
        
          I
          SOS
        
      
    
  
  ,




where ISOS is the SOS profile obtained after performing the combination calculation by a sum of squares on the sensitivities of the coils and the imagesos is achieved by multiplying the original image by ISOS of the coils.

In step (ii), the merging calculation by a sum of squares can also employ a combination method based on multiplication by the sum of complex conjugate (SCS) sensitivity weighted coil image, and the formula thereof is




  
    
      image
      SOS
    
    =
    
      
        
          ∑
          
            iCh
            =
            
              1
              ⁢
              ...
              ⁢
              
                
              
              ⁢
              nCh
            
          
        
        ⁢
        
          
            image
            iCh
          
          ·
          
            s
            iCH
            SOS
          
        
      
      =
      
        
          
          image
          
        
        ·
        
          I
          SOS
        
      
    
  
  ,
  
    where
    ⁢
    
      
    
    ⁢
    
      image
      SCS
    
  




where imageSCS is the image obtained by the combination method employing SCS, siChSOS is the sensitivities of the coils educed based on ISOS, and ISOS is as described above, and the formula for calculating





  
    s
    iCh
    SOS
  
  -
  
    
      image
      iCh
    
    /
    
      
        
          
            ∑
            iCh
          
          ⁢
          
            
              image
              iCh
            
            ·
            
              image
              iCh
              *
            
          
        
      
      .
    
  




The calculations in steps (i) and (ii), for being linear, can be expressed as a consecutive operation of two sparse matrixes. Therefore, the combination 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
            *
          
        
        .
      
    
  




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 combined 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 siChSOS and siChSOS are a Fourier transform pair.

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
                
              
            
          
        
      
    
    .
  




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
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
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 combination 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.
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 convolutions of channels and then transforming them into the image domain.
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 convolutions of the four channels:




  
    
      
        
          rawdata
          kCh
          full
        
        =
        
          
            ∑
            
              iCh
              =
              
                1
                ⁢
                ...
                ⁢
                
                  
                
                ⁢
                nCh
              
            
          
          ⁢
          
            
              rawdata
              iCh
              ppa
            
            ⊗
            
              F
              
                iCh
                ,
                kCh
              
            
          
        
      
    
    
      
        (
        1
        )
      
    
  




where rawdatakChfull represents the fully-sampled k-space data of the kChth channel fitted by the convolution. In the FIGURE kCh is 4. rawdataiChppa represents the under-sampled k-space data of the ichth channel used to fit. FiCh,kCh is the convolution kernel educed based on the fitting parameter. By performing Fourier transform on FiCh,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 imageiChfull is the fully-sampled image of the kChth channel, imageiChppa is the under-sampled image of the iChth channel, and fiCh,kCh is the Fourier transform of FiCh,kCh.

The second step of the method is converting the SOS channel combination calculation in the image domain into a linear calculation of images of individual channels and sensitivities of the coils.
The formula for the combination calculation by a sum of squares is:




  
    
      
        
          image
          SOS
        
        =
        
          
            
            image
            
          
          ⁢
          
            
              
                ∑
                
                  iCh
                  =
                  
                    1
                    ⁢
                    ...
                    ⁢
                    
                      
                    
                    ⁢
                    nCh
                  
                
              
              ⁢
              
                
                  image
                  iCh
                
                ·
                
                  image
                  iCh
                  *
                
              
            
          
        
      
    
    
      
        (
        3
        )
      
    
  




If the sensitivity of the iChth channel is siCh, the combination calculation via sum of squares can also be expressed as:




  
    
      
        
          image
          SOS
        
        =
        
          
            
              
              image
              
            
            ⁢
            
              
                
                  ∑
                  
                    iCh
                    =
                    
                      1
                      ⁢
                      ...
                      ⁢
                      
                        
                      
                      ⁢
                      nCh
                    
                  
                
                ⁢
                
                  
                    S
                    iCh
                  
                  ·
                  
                    S
                    iCh
                    *
                  
                
              
            
          
          =
          
            
              
              image
              
            
            ·
            
              I
              SOS
            
          
        
      
    
    
      
        (
        4
        )
      
    
  




where ISOS 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 combination calculation by a sum of squares is achieved by multiplying the original image by ISOS of coils.

The formula of a combination method based on the sum of complex conjugate sensitivity (SCS) weighted coil image is.




  
    
      
        
          image
          SCS
        
        =
        
          
            
              ∑
              
                iCh
                =
                
                  1
                  ⁢
                  ...
                  ⁢
                  
                    
                  
                  ⁢
                  nCh
                
              
            
            ⁢
            
              
                image
                iCh
              
              ·
              
                S
                iCh
                
                  SOS
                  *
                
              
            
          
          =
          
            
              
              image
              
            
            ·
            
              I
              SOS
            
          
        
      
    
    
      
        (
        5
        )
      
    
  




where imageSCS is the image obtained by the combination method employing SCS, and SiChSOS is the sensitivity of the coils educed based on ISOS. The formula is as follows:





  
    
      
        
          S
          iCh
          SOS
        
        =
        
          
            image
            iCh
          
          
            
              
                ∑
                iCh
              
              ⁢
              
                
                  image
                  iCh
                
                ·
                
                  image
                  iCh
                  *
                
              
            
          
        
      
    
    
      
        (
        6
        )
      
    
  




The equivalence of the two combination methods is easy to prove.
The third step in the inventive method is to merge the linear calculations in steps (i) and (ii) as a weighting matrix.
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.
The weighting matrix is expressed as follows:




  
    
      
        
          w
          kCh
        
        =
        
          
            ∑
            
              iCh
              =
              
                1
                ⁢
                ...
                ⁢
                
                  
                
                ⁢
                nCh
              
            
          
          ⁢
          
            
              f
              
                iCh
                ,
                kCh
              
            
            ·
            
              S
              iCh
              
                SOS
                *
              
            
          
        
      
    
    
      
        (
        7
        )
      
    
  




The accelerated calculation can be performed in the frequency domain.
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 SiChSOS and SiChSOS 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 SiChSOS in the frequency domain.

The weighting matrix can be used to calculate SNR losses of images.
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
        )
      
    
  




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.
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.

## Claims

1. A fast GeneRalized Autocalibrating Partially Parallel Acquisition (GRAPPA) image reconstruction method for magnetic resonance imaging, comprising the steps of:
(i) expressing s fitting calculation of a channel in the frequency domain as convolutions of channels and transforming the convolution into the image domain;
(ii) converting a SOS channel combination calculation in the image domain into a linear calculation of images of individual channels and sensitivities of the coils; and
(iii) merging the linear calculations in steps (i) and (ii) as a weighting matrix.

2. A method as claimed in claim 1, comprising in step (i), expressing the fitting calculation of a channel as convolutions of channels




  
    
      rawdata
      kCh
      full
    
    =
    
      
        ∑
        
          iCh
          =
          
            1
            ⁢
            ...
            ⁢
            
              
            
            ⁢
            nCh
          
        
      
      ⁢
      
        
          rawdata
          iCh
          ppa
        
        ⊗
        
          F
          
            iCh
            ,
            kCh
          
        
      
    
  
  ,




where rawdatakChfull is the fully-sampled k-space data of the kChth channel fitted after the convolution calculation, rawdataiChppa is the under-sampled k-space data of the iChth channel used to fit, and FiCh,kCh is the convolution kernel derived based on fitting parameters.

3. A method as claimed in claim 2 comprising performing a Fourier transform on FiCh,kCh, to transform the convolution into a point multiplication calculation in the image domain as




  
    
      image
      kCh
      image
    
    =
    
      
        ∑
        
          iCh
          =
          
            1
            ⁢
            ...
            ⁢
            
              
            
            ⁢
            nCh
          
        
      
      ⁢
      
        
          image
          iCh
          ppa
        
        ·
        
          f
          
            iCh
            ,
            kCh
          
        
      
    
  
  ,




where imagekChfull is the fully-sampled image of the kChth channel, imageiChppa is the under-sampled image of the iChth channel, and fiCh,kCh is the Fourier transform of FiCh,kCh.

4. A method as claimed in claim 1, comprising in step (ii), performing the combination calculation as a sum of squares as




  
    image
    SOS
  
  =
  
    
      
        
          ∑
          
            iCh
            =
            
              1
              ⁢
              
                
              
              ⁢
              ...
              ⁢
              
                
              
              ⁢
              nCh
            
          
        
        ⁢
        
          
            image
            iCh
          
          ·
          
            image
            iCh
            *
          
        
      
    
    .

5. A method as claimed in claim 4 comprising performing steps (i) and (ii) linearly as a consecutive action of two sparse matrices and, in step (iii), merging the linear calculations in steps (i) and (ii) by multiplying said two sparse matrices, as an accelerated calculation.

6. A method as claimed in claim 4, comprising multiplying the two sparse matrixes to obtain the weighting matrix as




  
    w
    kCh
  
  =
  
    
      ∑
      
        iCh
        =
        
          1
          ⁢
          
            
          
          ⁢
          ...
          ⁢
          
            
          
          ⁢
          nCh
        
      
    
    ⁢
    
      
        f
        
          iCh
          ,
          kCh
        
      
      ·
      
        
          s
          iCh
          
            SOS
            *
          
        
        .

7. A method as claimed in claim 4, wherein sensitivity of the iChth channel is siCh, and performing the combination calculation as a sum of squares as




  
    
      image
      SOS
    
    =
    
      
        
          
          image
          
        
        ⁢
        
          
            
              ∑
              
                iCh
                =
                
                  1
                  ⁢
                  
                    
                  
                  ⁢
                  ...
                  ⁢
                  
                    
                  
                  ⁢
                  nCh
                
              
            
            ⁢
            
              
                s
                iCh
              
              ·
              
                s
                iCh
                *
              
            
          
        
      
      =
      
        
          
          image
          
        
        ·
        
          I
          SOS
        
      
    
  
  ,




where ISOS is the SOS profile obtained after performing the combination calculation via sum of squares on the sensitivities of the coils, and imageSOS is achieved by multiplying the original image by ISOS of coils.

8. A method as claimed in claim 1 comprising, in step (ii), performing the merging calculation as a sum of squares as a combination method based on the sum of complex conjugate sensitivity (SCS) weighted coil images, as




  
    
      image
      SCS
    
    =
    
      
        
          ∑
          
            iCh
            =
            
              1
              ⁢
              
                
              
              ⁢
              ...
              ⁢
              
                
              
              ⁢
              nCh
            
          
        
        ⁢
        
          
            image
            iCh
          
          ·
          
            s
            iCh
            
              SOS
              *
            
          
        
      
      =
      
        
          
          image
          
        
        ·
        
          I
          SOS
        
      
    
  
  ,




where imageSCS is the SOS profile.

9. A method as claimed in claim 8, comprising calculating siChSOS as




  
    s
    iCh
    SOS
  
  =
  
    
      image
      iCh
    
    /
    
      
        
          
            ∑
            iCh
          
          ⁢
          
            
              image
              iCh
            
            ·
            
              image
              iCh
              *
            
          
        
      
      .

10. A method as claimed in claim 8 comprising performing steps (i) and (ii) linearly as a consecutive action of two sparse matrices and, in step (iii), merging the linear calculations in steps (i) and (ii) by multiplying said two sparse matrices, as an accelerated calculation.

11. A method as claimed in claim 10, comprising performing the accelerated calculation in the frequency domain.

12. A method as claimed in claim 11, comprising performing the accelerated calculation in the frequency domain by after step (i), performing a convolution of data of individual channels and the Fourier transform of the sensitivity functions of coils again, and further performing a Fourier transform to obtain the merged image.

13. A method as claimed in claim 12, comprising performing the accelerated calculation in the frequency domain as




  
    
      rawdata
      SCSfull
    
    =
    
      
        ∑
        
          iCh
          =
          
            1
            ⁢
            
              
            
            ⁢
            ...
            ⁢
            
              
            
            ⁢
            nCh
          
        
      
      ⁢
      
        
          rawdata
          iCh
          full
        
        ⊗
        
          S
          iCh
          SOS
        
      
    
  
  ,




where SiChSOS and SiChsoso are a Fourier transform pair.

14. A method as claimed in claim 11, comprising using the weighting matrix to compare respective signal-to-noise ratio losses of images obtained by different types of reconstruction methods in the image domain and the frequency domain to obtain a relative signal-to-noise ratio:




  
    SNR
    rel
  
  =
  
    
      
        
          ∑
          
            k
            =
            
              1
              ⁢
              
                
              
              ⁢
              ...
              ⁢
              
                
              
              ⁢
              nCh
            
          
        
        ⁢
        
          
            w
            kCh
          
          ·
          
            s
            kCh
            SOS
          
        
      
      
        
          
            ∑
            
              k
              =
              
                1
                ⁢
                
                  
                
                ⁢
                ...
                ⁢
                
                  
                
                ⁢
                nCh
              
            
          
          ⁢
          
            
              
                
                
                  s
                  kCh
                  SOS
                
                
              
              2
            
            ·
            
              
                
                  ∑
                  
                    k
                    =
                    
                      1
                      ⁢
                      
                        
                      
                      ⁢
                      ...
                      ⁢
                      
                        
                      
                      ⁢
                      nCh
                    
                  
                
                ⁢
                
                  
                    
                    
                      w
                      kCh
                    
                    
                  
                  2
                
              
            
          
        
      
    
    .