Grappa algorithm matlab code

Introduce MRI parallel Imaging

Parallel imaging is the use of phased array coils (sometimes called multi coil arrays) for the purpose of faster scanning. For example, the purpose of using parallel Imaging is to reduce the scan time, such as magnetic resonance imaging scan in abdominal if use the parallel imaging method can reduce the patient’s time that holding breath.

Parallel imaging application

Parallel imaging is most useful for applications that have abundant SNR. One of the most successful applications has been to contrast-enhanced MR angiography.
Reducing the scan time is beneficial for catching the peak of the bolus and shorter breath-hold times.
Another general way to use parallel imaging is for artifact reduction for each train pulse sequence such as EPI or RARE. EPI suffers from geometric distortion due to off-resonant spins, whereas RARE suffers from blurring due to T2 decay. Parallel imaging reduces the total number of k-space lines, which in turn can reduce the echo train length.

Parallel imaging reconstruction algorithm catgorized

Parallel imaging methods can generally be categorized as either k-space methods or image space methods

• k-space lines method:
  Such as SMASH, GRAPPA, in which missing k-space lines are restored prior to the Fourier transform.
• image space method:
  Such as SENSE, in which aliasing is removed in image space after the Fourier transform.

What is Grappa

The Grappa algorithm means additional gradient data beyond the original set are measured to obtain the weights. This is called a ’self-calibrating’ technique; it utilizes a limited set of ‘auto-calibration signal’ (ACS) points in k-space to get the weights.

Picture Source: http://mriquestions.com/grappaarc.html

We can explain grappa algorithm that how many coils we have will calculate how many three-dimensional kernels. Each coil’s kernel will be used to fit the k-space data lines that are not sampled per channel, respectively. For example, we have eight coils, then we have to calculate eight three-dimensional kernels, and these eight three-dimensional will be used to fit the eight channel’s k-space data lines that not sampled.
The coil arrangement and choice of kernel size, all of them will influence the effect of reconstruction. From a physical point of view, The image domain changes gently, so the image domain can be represented by a large kernel, and transformed into the frequency domain, we can use a small kernel to fit the unsampled k-space data lines.

How to implement Grappa algorithm

Coil n …

ACQUIRE: Partial k-space data with centeral oversampling.
ESTIMATE: Missing lines of k-space and regional coil sensitivities using harmonics.
GENERTE: Individual coil images.
COMBINE: Into final magnitude image using sum of squares.

Grappa reconstruction

• Sample central k-space lines fully (smooth, low frequency)
• Calculate weighting factors for each coil
• Estimate missing k-space lines
• Combine corrected images coming from different coil

Matlab code

Load raw-data

%####################################################################
%#######################  Copyright (c) 2019   ######################
%#####  Created on 2019-07-24                         ###############
%#####  Author:Yanglei_Wu                             ###############
%#####  Email:yanglei.wuu@gmail.com                   ###############
%#####  GRAPPA                                        ###############
%#####  Version 1.0                                   ###############
%#####  Copyright (c) 2019 Yanglei.Wu. All rights reserved  #########
%####################################################################
%%
clear all;
close all;
%% read rawdata
data = readBRF(  '..\(20190711125332033)-2019030201-C289-flash2dme_mecho-[T1_FLASH_AXIAL_mecho_256_256]\BatchRecon.brf');

About readBRF is a function for load row-data, you will use your function to load raw-data.

Initialization

%% initialization
Input.acs = 32;        %auto-calibration scan
Input.R = 3;          %accelerated factor
kernel.size_p = 4;    %kernel size in phase encoding direction
kernel.size_f = 5;    %kernel size in frequence encoding direction
dataG.Npe = data.m_pStatic.Views;     %number of phase encoding lines
dataG.Nfe = data.m_pStatic.al / 2.0;  %number of frequency encoding  lines
dataG.Nco = data.m_pStatic.Partition;     %number of coils

Undersampling the k-space data of each coils

%% undersampling the k-space data of each coil
Input.acs = fix(Input.acs / 2) * 2;
dataG.first = floor(dataG.Npe / 2) - fix(Input.acs / 2);      %the fist  acs line
dataG.last = dataG.first + Input.acs -1;       %the last acs line
dataG.kspace = zeros(dataG.Nco,dataG.Npe,dataG.Nfe);
dataG.kspace(:,1:Input.R:dataG.first,:) =  squeeze(data.kspace(:,:,:,1,1,1:Input.R:dataG.first,1:2:end));
dataG.kspace(:,dataG.first:dataG.last,:) =  squeeze(data.kspace(:,:,:,1,1,dataG.first:dataG.last,1:2:end));
dataG.kspace(:,dataG.last:Input.R:end,:) =  squeeze(data.kspace(:,:,:,1,1,dataG.last:Input.R:end,1:2:end));

Calculate kernel in ACS

%% calculateing kernel in ACS
kernel.size_p_center = 2 * ceil(kernel.size_p / 6);
kernel.size_f_center = 1 + floor(kernel.size_f / 2.5);
kernel.res = [];
for i = dataG.first:(dataG.last - kernel.size_p + 1)
    for j = 1:(dataG.Nfe - kernel.size_f + 1)
        kernel.a = dataG.kspace(:,i:Input.R:i + kernel.size_p - 1,j:j +  kernel.size_f - 1);
        kernel.b = squeeze(reshape(kernel.a,1,dataG.Nco *  ceil(kernel.size_p / Input.R),kernel.size_f));
        kernel.c = permute(reshape(kernel.b,dataG.Nco * ceil(kernel.size_p  / Input.R) * kernel.size_f,1),[2 1]);
        kernel.res = [kernel.res;kernel.c];
    end
end
for i = 1:(Input.R - 1)
    kernel.calres1 = dataG.kspace(:,(dataG.first + floor(kernel.size_p/3)  + i - 1):(dataG.last + floor(kernel.size_p/3) - kernel.size_p +  i),ceil(kernel.size_f/2):dataG.Nfe - floor(kernel.size_f/2));
    kernel.calres(:,:,i) =  permute(squeeze(reshape(permute(kernel.calres1,[1 3  2]),dataG.Nco,1,(dataG.Nfe - kernel.size_f + 1)*(Input.acs -  kernel.size_p + 1))),[2 1]);
    for j=1:dataG.Nco
        kernel.result(:,j,i) = inv(kernel.res' * kernel.res) * kernel.res'  * kernel.calres(:,j,i);
    end
end

Estimate missing k-space lines

%% Estimate missing k-space lines
out.kspace = dataG.kspace;
for i = 1:Input.R:(dataG.first - kernel.size_p + 1)
    for j = 1:(dataG.Nfe - kernel.size_f + 1)
        estimate.a = dataG.kspace(:,i:Input.R:i + kernel.size_p - 1,j:j +  kernel.size_f - 1);
        estimate.b = squeeze(reshape(estimate.a,1,dataG.Nco *  ceil(kernel.size_p / Input.R),kernel.size_f));
        estimate.c = permute(reshape(estimate.b,dataG.Nco *  ceil(kernel.size_p / Input.R) * kernel.size_f,1),[2 1]);
        for k = 1:dataG.Nco
            for l = 1:Input.R - 1
                out.kspace(k,i + kernel.size_p_center - 2 + l,j +  kernel.size_f_center - 1) = estimate.c * kernel.result(:,k,l);
            end
        end
    end
end
for i = (dataG.last):Input.R:(dataG.Npe - kernel.size_p + 1)
    for j = 1:(dataG.Nfe - kernel.size_f + 1)
        estimate.a = dataG.kspace(:,i:Input.R:i + kernel.size_p - 1,j:j +  kernel.size_f - 1);
        estimate.b = squeeze(reshape(estimate.a,1,dataG.Nco *  ceil(kernel.size_p / Input.R),kernel.size_f));
        estimate.c = permute(reshape(estimate.b,dataG.Nco *  ceil(kernel.size_p / Input.R) * kernel.size_f,1),[2 1]);
        for k = 1:dataG.Nco
            for l = 1:Input.R - 1
                out.kspace(k,i + kernel.size_p_center - 2 + l,j +  kernel.size_f_center - 1) = estimate.c * kernel.result(:,k,l);
            end
        end
    end
end

Reconstruction

image = permute(out.kspace,[2 3 1]);
im = fftshift(fft2(fftshift(image)));
figure,imshow(abs(im(:,:,1)),[]);
im_comb = sqrt(sum(abs(im).^2,3));
figure,imshow(im_comb,[]);

Result

Different ACS, different acceleration speeds will reconstruct images of different quality with different signal to noise ratios.

Reference

[1] Handbook of MRI Pulse sequence 1st Edition. Matt Bernstein etc. 2004
[2] Magnetic Resonance Imaging Physical Principles and Sequence Design 2e Robert W. Brown 2014