[1]:
# -*- coding: utf-8 -*-
# Copyright 2023 Physikalisch-Technische Bundesanstalt (PTB)
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Authored by: Christoph Kolbitsch (PTB),
# Edited by: Margaret Duff (STFC- UKRI), Letizia Protopapa (STFC- UKRI)
[2]:
# Make sure figures appears inline
%matplotlib widget
Dynamic MR Reconstruction#
MR can be used to capture dynamic processes. These dynamic processes can be used to obtain information about quantitative (bio-)physical parameters. A common parameter important for clinical diagnosis is the longitudinal relaxation time \(T_1\). In order to quantify \(T_1\) a MR-preparation pulse is used (e.g. inversion pulse) to disturb the equilibrium state of the MR spin system. After the preparation pulse the spin system relaxes back to the equilibrium state and this dynamic process can be described by the relaxation time \(T_1\). The dynamic acquisition after the preparation pulse therefore needs to be fast enough to capture the dynamic changes accurately.
Here we are going to reconstruct images covering such a dynamic process. Data is acquired continuously using a Golden radial acquisition scheme and we are going to split this data into different time frames and then reconstruct this data as 2D+t data block. The motivation for reconstructing one block is that we can utilise regularisation between different time points to reduce undersampling artefacts.
Because we are going to do MR image reconstruction we also need SIRF for this.
The imaged object is a phantom which consists of 9 tubes with different types of gel in them. Each tube has a different \(T_1\) time and hence the signal from each tube behaves differently over time.
CIL Version 22.2.0#
SIRF Version 3.4.0#
[3]:
import cil
print(cil.__version__)
import sirf
print(sirf.__version__)
22.2.0
3.4.0
The Dataset#
This requires the dataset 2D_Dyn_GRad.h5 from https://doi.org/10.5281/zenodo.7903233.
Please download the data and update the filename variable below to point to where you have the data saved:
[4]:
filename = './2D_Dyn_GRad.h5'
[5]:
# import engine module
import sirf.Gadgetron as pMR
# import CIL functionality for visualisation and iterative reconstruction
from cil.optimisation.algorithms import FISTA
from cil.optimisation.functions import LeastSquares, ZeroFunction
from cil.framework import BlockDataContainer
from cil.optimisation.operators import BlockOperator, ZeroOperator
# import further modules
import numpy as np
import matplotlib.pyplot as plt
Load and inspect the data#
Let’s load the data and have a look at the data acquisition over time.
[6]:
# Load in the data
acq_data = pMR.AcquisitionData(filename)
Started reading acquisitions from ./2D_Dyn_GRad.h5
0%..10%..20%..31%..40%..50%..61%..70%..80%..90%..
Finished reading acquisitions from ./2D_Dyn_GRad.h5
Data is aquired by sampling radial angles in k-space in an order determined by the golden ratio. Here we plot the first 1,2,4, and 20 angles.
[18]:
# Get k-space dimensions
kdim = acq_data.shape
# Get trajectory
ktraj = pMR.get_data_trajectory(acq_data)
ktraj = np.reshape(ktraj, (kdim[0], kdim[2], 2))
# Plot first 1,2,4 and 20 radial lines
plot_nrad = [1, 2, 4, 20]
titles=['1st angles', "First 2 angles", "First 4 angles", "First 20 angles"]
fig, ax = plt.subplots(1,4,figsize=(8,2))
for cax in ax.flatten():
cax.axis([-0.6, 0.6, -0.6, 0.6])
cax.set_yticks([])
for jnd in range(4):
for ind in range(plot_nrad[jnd]):
ax[jnd].plot(ktraj[ind,:,0], ktraj[ind,:,1], 'k-')
ax[jnd].set_title(titles[jnd])
Before we continue with the dynamic reconstruction we need to calculate the coil sensitivity maps:
[8]:
# Calculate coil sensitivity maps
csm = pMR.CoilSensitivityData()
csm.smoothness = 100
csm.calculate(acq_data)
Prepare dynamic reconstruction#
We are going to split the data into Ndyn dynamic frames each with Nrad_per_dyn radial lines. We chose N_dyn small here to speed things up and won’t be using all the data (1440 radial spokes) but feel free to adapt this.
[9]:
# Parameters for dynamics
Ndyn = 12
Nrad_per_dyn = 16
Now we are going to select the data acq_data and define the acquisition model E_dyn for each dynamic frame. A preliminary (i.e. pseudo-inverse) reconstruction is calculated in im_dyn_inv.
[10]:
# Go through each dynamic, create corresponding k-space and acquisition model
acq_dyn = [0] * Ndyn
E_dyn = [0] * Ndyn
im_dyn_inv = [0] * Ndyn
for ind in range(Ndyn):
acq_idx_dyn = np.linspace(ind*Nrad_per_dyn, (ind+1)*Nrad_per_dyn-1, Nrad_per_dyn)
acq_dyn[ind] = acq_data.get_subset(acq_idx_dyn)
acq_dyn[ind].sort_by_time()
# Create acquisition model
E_tmp = pMR.AcquisitionModel(acqs=acq_dyn[ind], imgs=csm)
E_tmp.set_coil_sensitivity_maps(csm)
im_dyn_inv[ind] = E_tmp.inverse(acq_dyn[ind])
E_dyn[ind] = pMR.AcquisitionModel(acqs=acq_dyn[ind], imgs=im_dyn_inv[ind])
E_dyn[ind].set_coil_sensitivity_maps(csm)
In order to be able to reconstruct all the dynamics at the same time, we have to create a diagonal matrix with all the acquisition models along the diagonal and zero-operators everywhere else. The acquisition data is also stacked in a block data container. The problem we want to solve in the end is:
\[\begin{split}\min_{[x_0\cdots x_{NDyn}]^T} \left\Vert \begin{bmatrix}
E_0 & 0 & 0\\
0 & \ddots & 0 \\
0 & 0 & E_{NDyn}
\end{bmatrix}
\begin{bmatrix}
x_0 \\
\vdots \\
x_{NDyn}
\end{bmatrix} -
\begin{bmatrix}
a_0 \\
\vdots \\
a_{NDyn}
\end{bmatrix}
\right\Vert_2^2\end{split}\]
where \(a\) is acq_data, \(E\) is E_dyn and \(x\) are the reconstructed dynamic images.
[11]:
# AcquisitionModel for each dynamic
# Create zero operator
ZOp = ZeroOperator(E_dyn[0].domain_geometry(), E_dyn[0].range_geometry())
E_dyn_diag = [ZOp,]*(Ndyn*Ndyn)
for ind in range(Ndyn):
E_dyn_diag[ind*Ndyn+ind] = E_dyn[ind]
A = BlockOperator(*E_dyn_diag, shape=(Ndyn, Ndyn))
# Put together all the raw k-space data for each motion state in a BlockDataContainer
acq_dyn_block = BlockDataContainer(*acq_dyn)
Initialise the reconstruction with the psuedo-inverse:
[12]:
# Starting image
x_init = A.adjoint(acq_dyn_block)
Dynamic reconstruction without regularisation#
We use CIL and the FISTA algorithm to solve the least squares problem
[13]:
# Set up optimisation
x_init_zero = x_init.copy()
x_init_zero.fill(0.0)
# Objective function
f = LeastSquares(A, acq_dyn_block, c=1)
f.calculate_Lipschitz()
G = ZeroFunction()
# Set up FISTA for least squares
fista = FISTA(initial=x_init_zero, f=f, g=G)
fista.max_iteration = 60
fista.update_objective_interval = 10
# Run FISTA
fista.run(verbose=True)
# Get result
im_fista = fista.get_output()
Iter Max Iter Time/Iter Objective
[s]
0 60 0.000 8.05277e-03
10 60 7.840 1.42176e-04
20 60 7.683 2.79802e-05
30 60 7.713 1.93923e-05
40 60 7.783 1.72206e-05
50 60 7.757 1.67181e-05
60 60 7.674 1.65991e-05
-------------------------------------------------------
60 60 7.674 1.65991e-05
Stop criterion has been reached.
With early stopping the results look look better:
[14]:
# Set up optimisation
x_init_zero = x_init.copy()
x_init_zero.fill(0.0)
# Objective function
f = LeastSquares(A, acq_dyn_block, c=1)
f.calculate_Lipschitz()
G = ZeroFunction()
# Set up FISTA for least squares
fista_early_stopping = FISTA(initial=x_init_zero, f=f, g=G)
fista_early_stopping.max_iteration = 10
fista_early_stopping.update_objective_interval = 10
# Run FISTA
fista_early_stopping.run(verbose=True)
# Get result
im_fista_early_stopping = fista_early_stopping.get_output()
Iter Max Iter Time/Iter Objective
[s]
0 10 0.000 8.05277e-03
10 10 7.387 1.42176e-04
-------------------------------------------------------
10 10 7.387 1.42176e-04
Stop criterion has been reached.
Dynamic reconstruction with TV regularisation#
In order to add TV-regularisation we need a bit of code which handles the diagonal structure of our image reconstruction problem. The main idea is to take the block data container and transform it into a 3d image, i.e. translate the “block-dimension” into the “channel-dimension” of a single image. Then we can apply the TV operator and afterwards we split it up again into different blocks. The code is found in tv_for_bdc.py so we simply import it.
[15]:
import sys
sys.path.append('CIL-User-Showcase/dynamic_mr_recon/')
from tv_for_bdc import FGP_TV_BDC
[16]:
# Set up optimisation
x_init_zero = x_init.copy()
# Objective function
f = LeastSquares(A, acq_dyn_block, c=1)
f.calculate_Lipschitz()
alpha = 2e-5
G = FGP_TV_BDC(x_init_zero, alpha)
# Set up FISTA for least squares
fista = FISTA(initial=x_init_zero, f=f, g=G)
fista.max_iteration = 50
fista.update_objective_interval = 10
# Run FISTA
fista.run(verbose=True)
# Get result
im_fista_tv = fista.get_output()
Iter Max Iter Time/Iter Objective
[s]
/opt/SIRF-SuperBuild/INSTALL/python/cil/plugins/ccpi_regularisation/functions/regularisers.py:125: ComplexWarning: Casting complex values to real discards the imaginary part
in_arr = np.asarray(x.as_array(), dtype=np.float32, order='C')
0 50 0.000 8.13893e-03
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
10 50 7.901 1.57843e-04
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
20 50 8.441 4.00597e-05
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
30 50 8.611 3.31214e-05
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
40 50 8.736 3.09825e-05
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
doing proximal_numpy, type of in_arr is float32
doing _fista_on_dual_rof, type of in_arr is float32
50 50 8.860 3.05365e-05
-------------------------------------------------------
50 50 8.860 3.05365e-05
Stop criterion has been reached.
Visualise results#
Here we plot in the rows the reconstruction from the psuedo inverse, least squares reconstruction with fista, and the TV regularised solution. The columns show reconstructions from 6 time frames, spread across the total range of time frames. The final column shows one column of the reconstruction at each step in time.
[17]:
# List of images to 2D+t matrix
im_list = [im_dyn_inv, im_fista, im_fista_early_stopping, im_fista_tv]
im_2dt = []
for ind in range(len(im_list)):
cim = []
for jnd in range(len(im_list[ind])):
cim.append(np.abs(np.squeeze(im_list[ind][jnd].as_array()))[:,:,np.newaxis])
im_2dt.append(np.concatenate(cim, axis=2))
# Normalise images
for ind in range(len(im_2dt)):
im_2dt[ind] = (im_2dt[ind] - np.min(im_2dt[ind]))/(np.max(im_2dt[ind]) - np.min(im_2dt[ind]))
# Plot images
nplots = min(Ndyn, 6)
fig, ax = plt.subplots(len(im_2dt), nplots+1, figsize=(16,12), constrained_layout = True, sharey=True)
for ind in range(nplots):
cidx = int(ind / nplots * Ndyn)
for jnd in range(len(im_2dt)):
ax[jnd,ind].imshow(im_2dt[jnd][:,:,cidx], vmin=0, vmax=1.0, cmap='inferno')
ax[jnd,ind].set_xlabel('x')
ax[jnd,ind].set_ylabel('y')
if ind == 0:
ax[0,ind].set_ylabel('Pseudo Inv \n y')
ax[1,ind].set_ylabel('FISTA \n y')
ax[2,ind].set_ylabel('FISTA early stopping \n y')
ax[3,ind].set_ylabel('FISTA + TV \n y')
# Last column is a plot of space vs time
for jnd in range(len(im_2dt)):
ax[jnd,nplots].imshow(im_2dt[jnd][64,:,:], vmin=0, vmax=1.0, cmap='inferno', aspect=0.2)
ax[jnd,nplots].set_xlabel('Time')
ax[jnd,nplots].set_ylabel('y')
plt.tight_layout=True
[ ]: