Skip to content

MRI: generalized iterative reconstruction

Jump to bottom
Fa-Hsuan Lin edited this page Aug 20, 2023 · 4 revisions

This is the first example of using a collection of routines to reconstruct MR images using an iterative solver. These routines are expected to reconstruct images based on data encoded by either spatially linear or nonlinear gradient fields, including PatLoc or O-space imaging.

Introduction

Formulating the MRI data acquisition (encoding) and image reconstruction (decoding) as a linear equation, an iterative numerical solver can be used to reconstruct images that involve a versatile combination of receiver coil arrays, k-space trajectories, spatial encoding magnetic fields, and the presence of nuisance fields. Specifically, the iterative solver is based on the conjugated-gradient algorithm, which has been reported to be a generalized approach in sensitivity encoded MRI.

Data:

Definitions

To satisfy the Nyquist sampling, the range of the (linear) gradient field and the k-space coordinates must be within the following ranges. Note that the maximum and the minimum of the (linear) gradient field depends on the image matrix size.

For an image matrix of n-by-n, the ranges for (linear) gradient fields and k-space coordinates are:

  • Linear gradient fields: [-n/2 n/2-1]
  • K-space coordinate: [-1 1)

For example, for an image of 128-by-128 pixels, the minimum and the maximum of the (linear) gradient fields are -64 and 63, respectively. The minimum and the maximum of the k-space coordinates are -1 and 0.9844 (1 - 1/64), respectively.

Applications

In addition to the simple example below, these routines can be used to various applications.

MRI reconstruction: fully sampled data in a Cartesian trajectory with a volume coil

This example reconstructs an image using a one-channel coil with a homogeneous sensitivity map. It samples data in a rectilinear Cartesian k-space with full sampling.

Codes explanation

Here we simulated that the image x was spatially modulated by a volume RF coil with homogeneous sensitivity (all ones).

        for i=1:size(b1,3)
            ORIG(:,:,i)=imresize(x.*b1(:,:,i),matrix{m_idx});
        end;

The spatial encoding magnetic fields were simulated as linearly varying magnetic field along horizontal and vertical directions.

        n_freq=matrix{m_idx}(2);
        n_phase1=matrix{m_idx}(1);
        
        [grid_freq,grid_phase]=meshgrid([-floor(n_freq/2):ceil(n_freq/2)-1],[-floor(n_phase1/2) :1:ceil(n_phase1/2)-1]);
        grid_freq=fmri_scale(grid_freq,ceil(n_freq/2)-1,-floor(n_freq/2));
        grid_phase=fmri_scale(grid_phase,ceil(n_phase1/2)-1,-floor(n_phase1/2));
        
        G_general(:,:,1)=grid_phase;
        G_general(:,:,2)=grid_freq;

K-space sampling coordinates were specified to cover the [-1 1) range in kx and ky directions.

                K{g_idx}=zeros(n_phase1,n_freq);
                K{g_idx}(1:R{r_idx}(1):end,1:R{r_idx}(2):end)=1;
                
                [k_phase1,k_freq, k_phase2]=ind2sub(matrix{1},find(K{g_idx}));
                
                k_phase1=((k_phase1-1)-n_phase1./2)./n_phase1.*2;
                k_freq=((k_freq-1)-n_freq./2)./n_freq.*2;
                K_general=[k_phase1(:) k_freq(:)];

The actual MRI signal encoding and decoding were then simulated and reconstructed, respectively.

            tic;
            [ACC_enc]=mri_patloc_tdr_forward_op('n_freq',matrix{m_idx}(2),'n_phase',matrix{m_idx}(1),'S',b1,'flag_display',1,'K_general',K_general,'G_general',G_general,'X',x0);
            time_encoding_matrix=toc;
            
            % time-domain recon::decoding
            tic;
            [recon_opt,recon_history_opt,error_opt]=itdr4_core_ktraj_cg('n_freq',matrix{m_idx}(2),'n_phase',matrix{m_idx}(1),'Y',ACC_enc,'S',b1,'flag_display',0,'K_general',K_general,'G_general',G_general,'iteration_max',50,'epsilon',1e-80);
            time_recon=toc;

The last part of the code was about combining images across MRI channels. As there was only one channel in this simulation, it was trivial.

Results

The following shows the linear gradient fields along the x (left-right) and y (up-down) directions. Note that the maximum and the minimum values matches the requirement of Nyquist sampling ([-64 +63] for an image matrix of 128-by-128 pixels).

The k-space sampling pattern can be shown by the following.

plot(K_general(:,1),K_general(:,2),'b.')

The zoom-in figure better illustrates the discrete sampling of the k-space in kx (left-right) and ky (up-down) directions.

The image to be simulated (ground truth) is an axial slice of the brain. The reconstructed image shows a visually identical image. These were shown in the following figure left and right.

figure;
imagesc(x); axis off image;
figure;
imagesc(Recon_opt); axis off image;

The history of the iterative reconstruction was shown below. The first reconstruction is already the final reconstruction: no much difference in the four iterative calculations.

The reconstruction shows

Clone this wiki locally