notebooks/MR/d_undersampled_reconstructions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Iterative reconstruction of undersampled MR data\n",
    "\n",
    "This demonstration shows how to hande undersampled data\n",
    "and how to write a simple iterative reconstruction algorithm with\n",
    "the acquisition model.\n",
    "\n",
    "This demo is a 'script', i.e. intended to be run step by step in a\n",
    "Python notebook such as Jupyter. It is organised in 'cells'. Jupyter displays these\n",
    "cells nicely and allows you to run each cell on its own."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First version: 27th of March 2019  \n",
    "Updated: 26nd of June 2021  \n",
    "Author: Johannes Mayer, Christoph Kolbitsch\n",
    "\n",
    "CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF).  \n",
    "Copyright 2015 - 2017 Rutherford Appleton Laboratory STFC.  \n",
    "Copyright 2015 - 2017 University College London.  \n",
    "Copyright 2015 - 2017, 2019, 2021 Physikalisch-Technische Bundesanstalt.\n",
    "\n",
    "This is software developed for the Collaborative Computational\n",
    "Project in Positron Emission Tomography and Magnetic Resonance imaging\n",
    "(http://www.ccppetmr.ac.uk/).\n",
    "\n",
    "SPDX-License-Identifier: Apache-2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% make sure figures appears inline and animations works\n",
    "%matplotlib widget\n",
    "\n",
    "# Setup the working directory for the notebook\n",
    "import notebook_setup\n",
    "from sirf_exercises import cd_to_working_dir\n",
    "cd_to_working_dir('MR', 'd_undersampled_reconstructions')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "__version__ = '0.1.0'\n",
    "\n",
    "# import engine module\n",
    "import sirf.Gadgetron as pMR\n",
    "from sirf.Utilities import examples_data_path\n",
    "from sirf_exercises import exercises_data_path\n",
    "\n",
    "# import further modules\n",
    "import os, numpy\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "\n",
    "from tqdm.auto import trange"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Undersampled Reconstruction\n",
    "#### Goals of this notebook:\n",
    "\n",
    "- Write a fully sampled reconstruction on your own.\n",
    "- Obtain knowledge of how to deal with undersampled data.\n",
    "- User SIRF and Gadgetron to perform a GRAPPA reconstruction.\n",
    "- Implement an iterative parallel imaging SENSE reconstruction algorithm from scratch.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is just an auxiliary function\n",
    "def norm_array( arr ):\n",
    "    min_a = abs(arr).min()\n",
    "    max_a = abs(arr).max()\n",
    "\n",
    "    return (arr - min_a)/(max_a - min_a)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time to get warmed up again:\n",
    "Since we deal with undersampled data in this last section, we need to compare it to a reference.\n",
    "So we need to reconstruct the fully sampled dataset we encountered before.\n",
    "\n",
    "This is an ideal opportunity to test what we learned and employ the `pMR.FullSampledReconstructor` class from before."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Programming Task: Fully sampled reconstruction\n",
    "\n",
    "__Please write code that does the following:__\n",
    "- create a variable called `full_acq_data` of type `pMR.AcquisitionData` from the file `ptb_resolutionphantom_fully_ismrmrd.h5`\n",
    "- create a variable called `prep_full_data` and assign it the preprocessed data by calling the function `pMR.preprocess_acquisition_data` on our variable `full_acq_data`\n",
    "- create a variable called `recon` of type `pMR.FullySampledReconstructor()`\n",
    "- call the `set_input` method of `recon` on `prep_full_data` to assign our fully sampled dataset to our reconstructor\n",
    "- call the `process()` method of `recon` without arguments. \n",
    "- create a variable called `fs_image` and assign it the output of the `get_output` method of `recon`\n",
    "\n",
    "__Hint:__ if you call a function without arguments, don't forget the empty parenthesis.\n",
    "\n",
    "#### Don't look at the solution before you tried! \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# YOUR CODE GOES HERE\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# VALIDATION CELL\n",
    "fs_image_array = fs_image.as_array()\n",
    "fs_image_array = norm_array(fs_image_array)\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.set_cmap('gray')\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "ax.imshow( abs(fs_image_array[0,:,:]), vmin=0, vmax=1)\n",
    "ax.set_title('Fully sampled reconstruction')\n",
    "ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution Cell. Don't look if you didn't try!\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "data_path = exercises_data_path('MR', 'PTB_ACRPhantom_GRAPPA')\n",
    "filename_full_file = os.path.join(data_path, 'ptb_resolutionphantom_fully_ismrmrd.h5')\n",
    "\n",
    "full_acq_data = pMR.AcquisitionData(filename_full_file)\n",
    "prep_full_data = pMR.preprocess_acquisition_data(full_acq_data)\n",
    "\n",
    "recon = pMR.FullySampledReconstructor()\n",
    "recon.set_input(prep_full_data)\n",
    "recon.process()\n",
    "fs_image = recon.get_output()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# VALIDATION CELL\n",
    "fs_image_array = fs_image.as_array()\n",
    "fs_image_array = norm_array(fs_image_array)\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.set_cmap('gray')\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "ax.imshow( abs(fs_image_array[0,:,:]), vmin=0, vmax=1)\n",
    "ax.set_title('Fully sampled reconstruction')\n",
    "ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# LOADING AND PREPROCESSING DATA FOR THIS SET\n",
    "filename_grappa_file = os.path.join(data_path, 'ptb_resolutionphantom_GRAPPA4_ismrmrd.h5')\n",
    "acq_data = pMR.AcquisitionData(filename_grappa_file)\n",
    "preprocessed_data = pMR.preprocess_acquisition_data(acq_data)\n",
    "preprocessed_data.sort()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Is the data we loaded undersampled? %s' % preprocessed_data.is_undersampled())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% RETRIEVE K-SPACE DATA\n",
    "k_array = preprocessed_data.as_array()\n",
    "print('Size of k-space %dx%dx%d' % k_array.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% SELECT VIEW DATA FROM DIFFERENT COILS\n",
    "print('Size of k-space %dx%dx%d' % k_array.shape)\n",
    "\n",
    "num_channels = k_array.shape[1]\n",
    "print(num_channels)\n",
    "k_array = norm_array(k_array)\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.set_cmap('gray')\n",
    "for c in range(num_channels):\n",
    "    ax = fig.add_subplot(2,num_channels//2,c+1)\n",
    "    ax.imshow(abs(k_array[:,c,:]), vmin=0, vmax=0.05)\n",
    "    ax.set_title('Coil '+str(c+1))\n",
    "    ax.axis('off')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks pretty similar to what we had before, so we should be good."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we define a quick fft with sos coil combination to do a standard reconstruction.\n",
    "def our_fft( k_array ):\n",
    "    image_array = numpy.zeros(k_array.shape, numpy.complex128)\n",
    "    for c in range(num_channels):\n",
    "        image_array[:,c,:] = numpy.fft.fftshift( numpy.fft.ifft2( numpy.fft.ifftshift(k_array[:,c,:])))\n",
    "#     image_array = image_array/image_array.max()    \n",
    "    image_array = numpy.sqrt(numpy.sum(numpy.square(numpy.abs(image_array)),1))\n",
    "\n",
    "    return image_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we make a FFT of the data we looked at and compare it to our fully sampled image\n",
    "image_array_sos = our_fft(k_array)\n",
    "image_array_sos = norm_array(image_array_sos)\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,2,1)\n",
    "ax.imshow(abs(image_array_sos), vmin=0, vmax=1)\n",
    "ax.set_title('Undersampled')\n",
    "ax.axis('off')\n",
    "\n",
    "ax = fig.add_subplot(1,2,2)\n",
    "ax.imshow(abs(fs_image_array[0,:,:]), vmin=0, vmax=1)\n",
    "ax.set_title('Fully sampled')\n",
    "ax.axis('off')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question: \n",
    "Please answer the following question for yourself:\n",
    " - Why is the undersampled reconstruction squeezed, but covers the whole FOV?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOW LET'S LOOK WHICH PARTS ARE SAMPLED AND WHICH ARE LEFT OUT\n",
    "# kspace_encode_step_1 gives the phase encoding steps which were performed\n",
    "which_pe = preprocessed_data.get_ISMRMRD_info('kspace_encode_step_1')\n",
    "\n",
    "print('The following Phase Encoding (PE) points were sampled: \\n')\n",
    "print(which_pe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Observation: approximately 3 out of 4 phase encoding steps are missing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fill an array with 1 only if a datapoint was acquired\n",
    "sampling_mask = numpy.zeros([256,256])\n",
    "\n",
    "for pe in which_pe:\n",
    "    sampling_mask[pe,:] = 1\n",
    "\n",
    "# PLOT THE SAMPLING MASK    \n",
    "fig = plt.figure()\n",
    "plt.set_cmap('gray')\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "ax.imshow( sampling_mask, vmin=0, vmax=1)\n",
    "ax.set_title('Sampling pattern of a GRAPPA acquisition')\n",
    "plt.xlabel('Frequency encoding')\n",
    "plt.ylabel('Phase encoding')\n",
    "#ax.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We notice that $112 = \\frac{256}{4} + 48$.\n",
    "This means the around the center of k-space is densely sampled containing 48 readout lines.\n",
    "The outside of k-space is undersampled by a factor of 4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workaround: 'zero-fill' the k-space data\n",
    "Since the artifacts seem to be related to the shape of the data, let's just unsqueeze the shape into the correct one.  \n",
    "Somehow we need to supply more datapoints to our FFT. But can we just add datapoints? This will corrupt the data!\n",
    "\n",
    "Actually, a Fourier transform is just a sum weighted by a phase:\n",
    "\n",
    "$$\n",
    "f(x) = \\sum_k e^{j k x} \\cdot F(k)\n",
    "$$\n",
    "\n",
    "So it does not mind if we add a data point $F(k) = 0$ at a position we didn't sample!\n",
    "\n",
    "This means we make a larger array, __add our data in the correct spots__, and the gaps we fill with zeros. The correct spots are given by our sampling mask!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_shape = [sampling_mask.shape[0], num_channels, sampling_mask.shape[1]]\n",
    "zf_k_array = numpy.zeros(k_shape, numpy.complex128)\n",
    "\n",
    "for i in range(k_array.shape[0]):\n",
    "    zf_k_array[which_pe[i],:,:] = k_array[i,:,:]\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.set_cmap('gray')\n",
    "for c in range(num_channels):\n",
    "    ax = fig.add_subplot(2,num_channels//2,c+1)\n",
    "    ax.imshow(abs(zf_k_array[:,c,:]), vmin=0, vmax=0.05)\n",
    "    ax.set_title('Coil '+str(c+1))\n",
    "    ax.axis('off')\n",
    "\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reconstruct the zero-filled data and take a look\n",
    "zf_recon = our_fft( zf_k_array)\n",
    "zf_recon = norm_array(zf_recon)\n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,2,1)\n",
    "ax.imshow(abs(zf_recon), vmin=0, vmax=1)\n",
    "ax.set_title('Zero-filled Undersampled ')\n",
    "ax.axis('off')\n",
    "\n",
    "\n",
    "ax = fig.add_subplot(1,2,2)\n",
    "ax.imshow(abs(fs_image_array[0,:,:]), vmin=0, vmax=1)\n",
    "ax.set_title('Fully Sampled ')\n",
    "ax.axis('off')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Observation: \n",
    "Bummer. Now the shape is correct, however the artifacts are still present. \n",
    "\n",
    " - What artifacts appear in the zero-filled reconstruction?\n",
    " - Why are they artifacts all fine-lined?\n",
    " - How come they only appear in one direction?\n",
    "\n",
    "To get rid of these we will need some parallel imaging techniques."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coil Sensitivity Map computation\n",
    "Parallel imaging, this has something to do with exploiting the spatially varying coil sensitivities.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WHICH COILMAPS DO WE GET FROM THIS DATASET?\n",
    "csm = pMR.CoilSensitivityData()\n",
    "csm.smoothness = 50\n",
    "csm.calculate(preprocessed_data)\n",
    "csm_array = numpy.squeeze(csm.as_array())\n",
    "\n",
    "csm_array = csm_array.transpose([1,0,2])\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.set_cmap('jet')\n",
    "for c in range(csm_array.shape[1]):\n",
    "    ax = fig.add_subplot(2,num_channels//2,c+1)\n",
    "    ax.imshow(abs(csm_array[:,c,:]))\n",
    "    ax.set_title('Coil '+str(c+1))\n",
    "    ax.axis('off')\n",
    "plt.set_cmap('gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question: \n",
    "In practice we would want to use a weighted sum (WS) coil combination technique so any artifacts in the coilmap would directly translate into the combined image. But we didn't see any of the high frequency artifacts!\n",
    "\n",
    "Please answer the following question:\n",
    "\n",
    "\n",
    "- Why are there are artifacts in the reconstruction but not in the coilmaps? \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We learned before, that parallel imaging is easily able to get rid of undersampling factor R=4.\n",
    "But to have enough information to estimate coilmaps the center must be fully sampled.  \n",
    "\n",
    "Ergo a perfect acceleration by R=4 is not possible, one needs to spends some time to sample the center densely, to obtain coil sensitivities. Still, $\\frac{112}{256} = 0.44$, we acquired 56% faster.\n",
    "\n",
    "### GRAPPA Reconstruction\n",
    "GRAPPA is a parallel imaging technique which promises to get rid of undersampling artifacts.  \n",
    "So let's use one of our SIRF classes and see what it can do!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WE DO A GRAPPA RECONSTRUCTION USING SIRF\n",
    "\n",
    "recon = pMR.CartesianGRAPPAReconstructor()\n",
    "recon.set_input(preprocessed_data)\n",
    "recon.compute_gfactors(False)\n",
    "print('---\\n reconstructing...')\n",
    "\n",
    "recon.process()\n",
    "# for undersampled acquisition data GRAPPA computes Gfactor images\n",
    "# in addition to reconstructed ones\n",
    "grappa_images = recon.get_output()\n",
    "grappa_images_array = grappa_images.as_array()\n",
    "grappa_images_array = norm_array(grappa_images_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLOT THE RESULTS\n",
    "fig = plt.figure(figsize=(9, 4))\n",
    "\n",
    "ax = fig.add_subplot(1,3,1)\n",
    "ax.imshow(abs(zf_recon), vmin=0, vmax=1)\n",
    "ax.set_title('Zero-filled Undersampled ')\n",
    "ax.axis('off')\n",
    "\n",
    "ax = fig.add_subplot(1,3,2)\n",
    "ax.imshow(abs(grappa_images_array[0,:,:]), vmin=0, vmax=1)\n",
    "ax.set_title('GRAPPA')\n",
    "ax.axis('off')\n",
    "\n",
    "ax = fig.add_subplot(1,3,3)\n",
    "ax.imshow(abs(fs_image_array[0,:,:]), vmin=0, vmax=1)\n",
    "ax.set_title('Fully Sampled')\n",
    "ax.axis('off')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Well, that was very little code to perform a difficult task!__ That is because we sent our data off to The Gadgetron and they did all the work.\n",
    "### Question:\n",
    "In what respect did a GRAPPA reconstruction:\n",
    "\n",
    " * improve the resulting image?\n",
    " * deterioate the resulting image?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GREAT! Now we want to develop our own algorithm and be better than the GRAPPA reconstruction!\n",
    "## Urgh, let's rather not because we are annoyed by how much code we have to write all the time! Zero filling, coil combining, inverse FFTs. Frankly: terrible!\n",
    "\n",
    "We want to capture our entire imaging and reconstruction process in one single object and don't care about data structure. Also we don't want to have to sum over coil channels all the time and take care of zero filling, this is just too much work!  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we want to develop some reconstruction algorithm, we want to be able to go from an image $x$ __forward__ to k-space data $y$:\n",
    "$$ E: x \\rightarrow y,$$\n",
    "implicitly performing multiplication of the image with the coil sensitivities $C_c$ for each channel $c$ and performing an FFT:\n",
    "\n",
    "$$\n",
    "E x = y_c = \\mathcal{F}( C_c \\cdot x).\n",
    "$$\n",
    "\n",
    "In iterative image reconstruction we often to apply the so-called __backward__  or __adjoint__ to transform the k-space data into image space. This bundles doing the zero-filling, inverse FFT, and coil combination into one operation:\n",
    "$$ E^H: y \\rightarrow x,$$\n",
    "implicitly performing everything:  \n",
    "$$\n",
    "E^H y = x = \\sum_c C_c^*\\mathcal{F}^{-1}(y) \n",
    "$$\n",
    "\n",
    "$E^H$ in this case is the hermitian conjugate of the complex valued operator $E$. It is the combination of transposing and complex conjugation of a matrix: $ E^H = (E^T)^* $. Note, that this is not generally the inverse: $ E^H \\neq E^{-1}$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Enter: AcquisitionModel\n",
    "\n",
    "In SIRF there exists something called `AcquisitionModel`, in the literature also referenced as and Encoding operator $E$, *E* for encoding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOW WE GENERATE THE ACQUISITION MODEL\n",
    "E = pMR.AcquisitionModel(preprocessed_data, grappa_images)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We need help again to see what this thing here can do\n",
    "help(E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to supply coil info to the acquisition model we use the dedicated method\n",
    "E.set_coil_sensitivity_maps(csm)\n",
    "\n",
    "# Now we can hop back from k-space into image space in just one line:\n",
    "aq_model_image = E.backward( preprocessed_data )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well this is not much code any more. Suddenly implementing our own iterative algorithm seems feasible!\n",
    "### QUESTION\n",
    "BEFORE YOU RUN THE NEXT CELL AND LOOK AT THE PLOT:\n",
    "In the next plot the image stored in `aq_model_image_array` will be shown, i.e. $x = E^H y$.  \n",
    "Based on the discussion what the AcquisitionModel E does, what do you expect the reconstruction to look like?\n",
    "- Is it squeezed or is it the correct size?\n",
    "- Does it contain artifacts? If so, which ones?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aq_model_image_array = norm_array(aq_model_image.as_array())\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.set_cmap('gray')\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "ax.imshow(abs(aq_model_image_array[0,:,:]))\n",
    "ax.set_title('Result Backward Method of E ')\n",
    "ax.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Well, bummer again, the artifacts are still there!__ Of course, the acquisition model is just a compact version of our above code. We need something smarter to kill them off. But it got a bit more homogeneous, due to the weighted coil combination. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Reconstruction as an Inverse Problem\n",
    "## Iterative Parallel Imaging Reconstruction\n",
    "\n",
    "In order to employ parallel imagaing, we should look at image reconstruction as an inverse problem.\n",
    "By \"image reconstruction\" we actually mean to achieve the following equality:\n",
    "$$ E x = y,$$ or equivalently\n",
    "$$ E^H E \\, x = E^H y,$$ \n",
    "where $E$ is the encoding operator $x$ is the true image object and $y$ is the MR raw data we acquired. \n",
    "The task of image reconstruction boils down to optimizing the following function:\n",
    "$$ \\mathcal{C}(x) = \\frac{1}{2} \\bigl{|} \\bigl{|}  E \\, x - y \\bigr{|} \\bigr{|}_2^2  \\\\\n",
    "\\tilde{x} = \\min_x \\mathcal{C}(x)\n",
    "$$\n",
    "To iteratively find the minimum of this cost function, $\\tilde{x}$, we need to go through steps of the kind:\n",
    "1. have a first guess for our image (usually an empty image)\n",
    "2. generate k-space data from this guess and compute the discrepancy to our acquired data (i.e. evaluate the cost fuction).\n",
    "3. update our image guess based on the computed discrepancy s.t. the cost will be lowered.  \n",
    "\n",
    "__By iterating steps 2 and 3 many times we will end up at $\\tilde{x}$.__\n",
    "\n",
    "Is that going to be better than a GRAPPA reconstruction? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing Conjugate Gradient Descent SENSE\n",
    "\n",
    "Conjugate Gradient (CG) optimization is such an iterative procedure which will quickly result in finding $\\tilde{x}$.\n",
    "We can study the corresponding [Wikipedia Article](https://en.wikipedia.org/wiki/Conjugate_gradient_method#Description_of_the_problem_addressed_by_conjugate_gradients).\n",
    "\n",
    "This looks like our thang! \n",
    "\n",
    "For that we need to write a bit of code:\n",
    "- We already have this encoding operator `E` defined where we can go from image to k-space and back.\n",
    "- Now we need to implement our first guess\n",
    "- and somehow we loop through steps 2 and 3 updating our image s.t. the costs are lowered.\n",
    "\n",
    "They want to compute x in $Ax = b$, we want to compute x in $E^H E x = E^H y$.  \n",
    "This means we need to translate what it says on Wikipedia:\n",
    "- $x$ = reconstructed image\n",
    "- $A$ = $E^HE$\n",
    "- $b$ = $E^H y$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Programming task\n",
    "Please write code executing the following task:\n",
    "- define a fuction named `A_operator`\n",
    "- it should have one single argument `image`\n",
    "- it should return $E^H( E (image))$  \n",
    "\n",
    "__Hint 1:__ We defined `E` already. Use it methods `forward` and `backward`. `forward` goes from image space to k-space and `backward` the other way round.  \n",
    "__Hint 2:__ Short reminder on the syntax. The function should look like: \n",
    "```\n",
    "def function_name( arugment_name):\n",
    "    variable = code_that_does_something_with ( argument_name )\n",
    "    return variable\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here (this is as much space as you need!)\n",
    "# make sure the name of your function is A_operator\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Don't look at the solution before you tried! \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# With this guy we will write our optimization\n",
    "def A_operator( image ):\n",
    "    return E.backward( E.forward(image) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Back to Wikipedia!\n",
    "Now we have all the tools we need. Now let's write the code to optimize our cost function iteratively.  \n",
    "We don't care too much about maths, but we want the [algorithm](https://en.wikipedia.org/wiki/Conjugate_gradient_method#The_resulting_algorithm).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# our images should be the same shape as the GRAPPA output\n",
    "recon_img =  grappa_images \n",
    "\n",
    "# since we have no knowledge at all of what the image is we start from zero\n",
    "zero_array = recon_img.as_array()\n",
    "zero_array.fill(0)\n",
    "recon_img.fill(zero_array)\n",
    "\n",
    "# now name the variables the same as in the Wikipedia article:\n",
    "x = recon_img\n",
    "y = preprocessed_data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Programming task: Initialize Iterative Reconstruction\n",
    "Please write code executing the following task:\n",
    "- Initialize a variable `r` with `r` = $b - Ax$ (`r` stands for residual).   \n",
    "__Hint 0:__ Remember: $b=E^H y$. Don't forget we just defined the A operator!\n",
    "- Print the type of `r` using Pythons built-in `type` and `print` function. What type of r do you expect? Is it an image, or is it acquisition data?\n",
    "- After you wrote these two lines run your cell pressing `Ctrl+Enter`, to get the output of the print statement. This will tell you the class of `r`.\n",
    "- Afterwards, initialize a variable named `rr` with `rr` = $r^\\dagger r$. (`rr` stands for r times r).  `rr` is the value of the cost function by the way.  \n",
    "__Hint 1:__ No need to access any numpy arrays! Objects of type `sirf.Gadgetron.ImageData` have the method called ` norm()` giving you the square root of the quantity we are looking for.  \n",
    "__Hint 2:__ Python-Power: $c=a^b$ $\\equiv$ `c = a**b`.  \n",
    "- Initialize a variable `rr0` with the value of `rr` to store the starting norm of the residuals.\n",
    "- Initialize a variable `p` with the value of `r`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## WRITE YOUR CODE IN THIS CELL\n",
    "## Please make sure to name the variables correctly\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Don't look at the solution before you tried! #############################\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "############################################################################\n",
    "# this is our first residual\n",
    "r = E.backward( y ) - A_operator(x)\n",
    "\n",
    "# print the type\n",
    "print('The type of r is: ' + str( type(r) ) ) \n",
    "\n",
    "# this is our cost function at the start\n",
    "rr = r.norm() ** 2\n",
    "rr0 = rr\n",
    "\n",
    "# initialize p\n",
    "p = r\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we write down the algorithm\n",
    "\n",
    "# how many iterative steps do we want\n",
    "# how low should the cost be\n",
    "num_iter = 15\n",
    "sufficiently_small = 1e-7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#prep a container to store the updated image after each iteration, this is just for plotting reasons!\n",
    "data_shape = numpy.array( x.as_array().shape )\n",
    "data_shape[0] = num_iter\n",
    "array_with_iterations = numpy.zeros(data_shape, numpy.complex128)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# HERE WE RUN THE LOOP.\n",
    "\n",
    "print('Cost for k = 0: '  + str( rr/ rr0) )\n",
    "with trange(num_iter) as iters:\n",
    "  for k in iters:\n",
    "\n",
    "    Ap = A_operator( p )\n",
    "\n",
    "    alpha = rr / Ap.dot(p)\n",
    "\n",
    "    x = x + alpha * p\n",
    "\n",
    "    r = r - alpha * Ap\n",
    "\n",
    "    beta  = r.norm()**2 / rr\n",
    "    rr = r.norm()**2\n",
    "\n",
    "    p = r + beta * p\n",
    "\n",
    "    relative_residual = numpy.sqrt(rr/rr0)\n",
    "\n",
    "    array_with_iterations[k,:,:] = x.as_array()\n",
    "    iters.write('Cost for k = ' +str(k+1) + ': ' + str(relative_residual) )\n",
    "    iters.set_postfix(cost=relative_residual)\n",
    "\n",
    "    if( relative_residual  < sufficiently_small ):\n",
    "        iters.write('We achieved our desired accuracy. Stopping iterative reconstruction')\n",
    "        break\n",
    "\n",
    "\n",
    "if k is num_iter-1:    \n",
    "    print('Reached maximum number of iterations. Stopping reconstruction.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# See how the reconstructed image evolves\n",
    "fig = plt.figure()\n",
    "ims = []\n",
    "for i in range(k):\n",
    "    im = plt.imshow(abs( array_with_iterations[i,:,:]), animated=True)\n",
    "    ims.append([im])\n",
    "\n",
    "ani = animation.ArtistAnimation(fig, ims, interval=500, blit=True, repeat_delay=0)\n",
    "plt.show() \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## now check out the final result as a still image\n",
    "recon_arr = norm_array( x.as_array())\n",
    "\n",
    "plt.set_cmap('gray')\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "ax.imshow(abs(recon_arr[0,:,:]))\n",
    "ax.set_title('SENSE RECONSTRUCTION')\n",
    "ax.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's Plot\n",
    "recon_arr = x.as_array()\n",
    "\n",
    "fig = plt.figure(figsize=(9, 9))\n",
    "\n",
    "ax = fig.add_subplot(2,2,1)\n",
    "ax.imshow(abs(aq_model_image_array[0,:,:]))\n",
    "ax.set_title('UNDERSAMPLED RECONSTRUCTION')\n",
    "ax.axis('off')\n",
    "\n",
    "ax = fig.add_subplot(2,2,2)\n",
    "ax.imshow(abs(grappa_images_array[0,:,:]))\n",
    "ax.set_title('GRAPPA RECONSTRUCTION')\n",
    "ax.axis('off')\n",
    "\n",
    "ax = fig.add_subplot(2,2,3)\n",
    "ax.imshow(abs(recon_arr[0,:,:]))\n",
    "ax.set_title('SENSE RECONSTRUCTION')\n",
    "ax.axis('off')\n",
    "\n",
    "ax = fig.add_subplot(2,2,4)\n",
    "ax.imshow(abs(fs_image_array[0,:,:]))\n",
    "ax.set_title('FULLY SAMPLED RECONSTRUCTION')\n",
    "ax.axis('off')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question: Evaluation SENSE Reconstruction\n",
    "[So what is better, GRAPPA or SENSE?](https://www.youtube.com/watch?v=XVCtkzIXYzQ)  \n",
    "\n",
    "Please answer the following questions:\n",
    "- Where is the noise coming from? \n",
    "- Why has not every high-frequency artifact vanished?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And if you had typed the above code into your computer in 2001 and written a [paper](https://scholar.google.de/scholar?hl=de&as_sdt=0%2C5&q=Advances+in+sensitivity+encoding+with+arbitrary+k%E2%80%90space+trajectories&btnG=) on it, then 18 years later you had a good 1000 citations (plus 6k from the [previous one](https://scholar.google.de/scholar?hl=de&as_sdt=0%2C5&q=SENSE%3A+sensitivity+encoding+for+fast+MRI&btnG=)).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Undersampled Reconstruction\n",
    "#### Recap:\n",
    "\n",
    "In this notebook you\n",
    "- wrote your own fully sampled recon using SIRF.\n",
    "- saw the undersampling structure of GRAPPA files\n",
    "- discovered high-frequency undersampling artifacts.\n",
    "- implemented our own version of SENSE and beat (?) GRAPPA.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fin\n",
    "\n",
    "This was the last exercise. We hoped you learned some new things about MRI and had a pleasant experience with SIRF and Python.  \n",
    "See you later!\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython"
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}