notebooks/MR/f_create_undersampled_kspace.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "86ea4a11",
   "metadata": {},
   "source": [
    "# Create undersampled k-space\n",
    "This demonstration shows how to create different undersampled k-space data which can be used either directly for image reconstruction or used to simulate MR data acquisition of a new object.\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",
   "id": "1df941f7",
   "metadata": {},
   "source": [
    "First version: 18th of June 2021\n",
    "Author: Christoph Kolbitsch\n",
    "\n",
    "CCP SyneRBI Synergistic Image Reconstruction Framework (SIRF).  \n",
    "Copyright 2015 - 2021 Rutherford Appleton Laboratory STFC.  \n",
    "Copyright 2015 - 2021 University College London.  \n",
    "Copyright 2015 - 2021 Physikalisch-Technische Bundesanstalt.\n",
    "\n",
    "This is software developed for the Collaborative Computational Project in Synergistic Reconstruction for Biomedical Imaging \n",
    "(http://www.ccpsynerbi.ac.uk/).\n",
    "\n",
    "SPDX-License-Identifier: Apache-2.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dcd9fd7",
   "metadata": {},
   "source": [
    "In the previous MR notebooks such as `e_advanced_recon` we used an existing MR raw data file with regular Cartesian sampling and an undersampling factor (R) of 4. If we would like to repeate the notebook for a dataset with R =  6, we would need to go to the MR scanner and acquire new raw data. This of course is one of the strengths of __SIRF__ because it means all our developed code can be applied to real data in a straight forward way, but sometimes it would be nice to switch from R = 4 to R = 6 by adapting a parameter in our script. Or go from regular Cartesian undersampling to random Cartesian undersampling without having to implement a new trajectory on the MR scanner.\n",
    "\n",
    "This notebook will show how we can achieve this at least to a certain degree. The idea is to start with a fully sampled Cartesian data set and then select only a subset of acquired k-space lines for the image reconstruction. We will use a 2D Cartesian data set and hence we can select a different subset of $k_y$ points. \n",
    "\n",
    "Of course this approach has several limitations. We cannot got from a Cartesian to a non-Cartesian (e.g. radial) sampling pattern and we cannot adapt the overall FOV of the scan, but it is a start."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3b5ae3a",
   "metadata": {},
   "source": [
    "So let's think about what we need to do.\n",
    "\n",
    " * (A) We need a fully sampled Cartesian data set. Let's take `ptb_resolutionphantom_fully_ismrmrd.h5` which is in `exercises_data_path('MR', 'PTB_ACRPhantom_GRAPPA')`. We need to load the data and we will already call `preprocess_acquisition_data()` for this data.\n",
    " * (B) Then we need to find out, which $k_y$ points have been acquired and where the centre of k-space is (i.e. $k_y$ = 0), because for any undersampled MR acquisition it is a good idea to have a small fully sampled region around the k-space centre and carry out the undersampling in the higher k-space frequencies. \n",
    " * (C) Define a subset of the orignal $k_y$ points.\n",
    " * (D) Create a new `AcquisitionData` object with only the subset defined in (C).\n",
    " * (E) Do a simple reconstruction to check we did the right thing.\n",
    " \n",
    "After completing all the previous notebooks you are already MR reconstruction experts. Therefore, you should be able to do these steps by yourself. Give it a try and try to create a new k-space with an undersampling factor of 4 and a fully sampled central region of 10 $k_y$ points. Only have a look at the example solution below if you are stuck. A few hints to get you started:\n",
    "\n",
    " * Details on how to get information about the acquired k-space (e.g. which $k_y$ points have been acquired) can be find in the notebook `d_undersampled_reconstructions`.\n",
    " * We can define an empty `AcquisitionData` object using `acq_new = preprocessed_data.new_acquisition_data(empty=True)`.\n",
    " * We can select an `Acquisition` object from an existing `AcquisitionData` object using `cacq = existing_acq_data.acquisition(acq_index)`.\n",
    " * Finally we can add this `Acquisition` object to our so far empty `AcquisitionData` object using `acq_new.append_acquisition(cacq)`.\n",
    " * Once we have added all your `Acquisition` objects, we have to sort the data again: `acq_new.sort()`.\n",
    " * In order to find out, how to do a simple reconstruction have a look at the MR part of the introductory notebook `acquisition_model_mr_pet_ct`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f9ad622",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% make sure figures appears inline and animations works\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "add8c747",
   "metadata": {},
   "outputs": [],
   "source": [
    "__version__ = '0.1.1'\n",
    "\n",
    "import notebook_setup\n",
    "import numpy\n",
    "\n",
    "# import engine module\n",
    "import sirf.Gadgetron as mr\n",
    "from sirf_exercises import exercises_data_path\n",
    "\n",
    "# import further modules\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "657fc281",
   "metadata": {},
   "source": [
    "You should already have called the following script. If not, uncomment BOTH lines and run it now!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13f32125",
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%bash \n",
    "#bash ../../scripts/download_data.sh -m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dbf9057",
   "metadata": {},
   "source": [
    "# Utilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b0505d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First define some handy function definitions\n",
    "# To make subsequent code cleaner, we have a few functions here. You can ignore\n",
    "# ignore them when you first see this demo.\n",
    "\n",
    "def plot_2d_image(idx,vol,title,clims=None,cmap=\"viridis\"):\n",
    "    \"\"\"Customized version of subplot to plot 2D image\"\"\"\n",
    "    plt.subplot(*idx)\n",
    "    plt.imshow(vol,cmap=cmap)\n",
    "    if not clims is None:\n",
    "        plt.clim(clims)\n",
    "    plt.colorbar()\n",
    "    plt.title(title)\n",
    "    plt.axis(\"off\")\n",
    "\n",
    "def crop_and_fill(templ_im, vol):\n",
    "    \"\"\"Crop volumetric image data and replace image content in template image object\"\"\"\n",
    "    # Get size of template image and crop\n",
    "    idim_orig = templ_im.as_array().shape\n",
    "    idim = (1,)*(3-len(idim_orig)) + idim_orig\n",
    "    offset = (numpy.array(vol.shape) - numpy.array(idim)) // 2\n",
    "    vol = vol[offset[0]:offset[0]+idim[0], offset[1]:offset[1]+idim[1], offset[2]:offset[2]+idim[2]]\n",
    "    \n",
    "    # Make a copy of the template to ensure we do not overwrite it\n",
    "    templ_im_out = templ_im.copy()\n",
    "    \n",
    "    # Fill image content \n",
    "    templ_im_out.fill(numpy.reshape(vol, idim_orig))\n",
    "    return(templ_im_out)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "859cb547",
   "metadata": {},
   "source": [
    "## (A) Fully sampled k-space data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "061d9bb1",
   "metadata": {},
   "source": [
    "Load in fully sampled k-space data and preprocess it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d49b62f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load MR AcquisitionData\n",
    "mr_acq = mr.AcquisitionData(exercises_data_path('MR', 'PTB_ACRPhantom_GRAPPA') \n",
    "                            + '/ptb_resolutionphantom_fully_ismrmrd.h5' )\n",
    "preprocessed_data = mr.preprocess_acquisition_data(mr_acq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb6f4183",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate image\n",
    "recon = mr.FullySampledReconstructor()\n",
    "recon.set_input(preprocessed_data)\n",
    "recon.process()\n",
    "im_mr = recon.get_output()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55729b96",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display it\n",
    "plt.figure();\n",
    "plot_2d_image([1,1,1], numpy.abs(im_mr.as_array())[0,:,:], 'Original image', cmap=\"Greys_r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14ea559f",
   "metadata": {},
   "source": [
    "## (B) Find out which k-space points have been acquired"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7cbe868",
   "metadata": {},
   "source": [
    "We will get the information about the $k_y$ position for each `Acquisition`. Because we have to go through all the acquired data, this can take a bit of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66ecb338",
   "metadata": {},
   "outputs": [],
   "source": [
    "ky_index = preprocessed_data.get_ISMRMRD_info('kspace_encode_step_1')\n",
    "print(ky_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c12aa8e7",
   "metadata": {},
   "source": [
    "So we have got 256 phase encoding points $k_y$. Because this is a fully sampled Cartesian acquisition we can savely assume that the k-space centre is located bang in the middle, i.e. $k_y$ = 0 for `ky_index` = 128."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bef042c8",
   "metadata": {},
   "source": [
    "## (C) Define a subset of k-space data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e8c5660",
   "metadata": {},
   "source": [
    "Let's start with something easy. Define a subset for a regular undersampling factor R = 2 but with a fully sampled central k-space region of 20 $k_y$ points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1e3d979",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define an undersampling factor\n",
    "R = 4\n",
    "\n",
    "# Define the number of fully sampled k-space points in the k-space centre\n",
    "N_ctr = 10\n",
    "\n",
    "# and your k-space centre was in the middle, so \n",
    "ky0_index = len(ky_index)//2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5df19f0a",
   "metadata": {},
   "source": [
    "Let's first select the fully sampled k-space centre"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "332ca3be",
   "metadata": {},
   "outputs": [],
   "source": [
    "ky_index_subset =  numpy.arange(ky0_index-N_ctr//2, ky0_index+N_ctr//2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c6b86f2",
   "metadata": {},
   "source": [
    "Now we can add the rest of the data with an undersampling factor of R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "604dc36a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ky_index_subset = numpy.concatenate((ky_index_subset, numpy.arange(start=0, stop=len(ky_index), step=R)), axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dad8757",
   "metadata": {},
   "source": [
    "Of course, now we might have added points from the centre of k-space again. To make sure that no $k_y$ index occurs twice, we simply call `numpy.unique`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ce6b919",
   "metadata": {},
   "outputs": [],
   "source": [
    "ky_index_subset = numpy.unique(ky_index_subset)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e53226b",
   "metadata": {},
   "source": [
    "Now we can plot the original fully sampled and the new undersampled indices. Hint: zoom into the figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5375a95d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(ky_index, numpy.ones(ky_index.shape), 'bo')\n",
    "plt.plot(ky_index[ky_index_subset], numpy.ones(ky_index[ky_index_subset].shape), 'r.');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4eed69a",
   "metadata": {},
   "source": [
    "## (D) Create new acquisition data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c91e4550",
   "metadata": {},
   "source": [
    "Now we know which k-space points to select, we need to select them and create a new `AcquisitionData` object. We will follow the steps detailed at the beginning of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62074fe0",
   "metadata": {},
   "outputs": [],
   "source": [
    "acq_new = preprocessed_data.new_acquisition_data(empty=True)\n",
    "\n",
    "# Create raw data\n",
    "for jnd in range(len(ky_index_subset)):\n",
    "    cacq = preprocessed_data.acquisition(ky_index_subset[jnd])\n",
    "    acq_new.append_acquisition(cacq)\n",
    "\n",
    "acq_new.sort()     "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd2f1b4b",
   "metadata": {},
   "source": [
    "## (E) Simple reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dadab562",
   "metadata": {},
   "source": [
    "Now we will do a simple reconstruction by defining and `AcquisitionModel` based on the `AcquisitionData` and then call `backward()` (i.e. Fourier transform)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b80fdd47",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original data\n",
    "csm_orig = mr.CoilSensitivityData()\n",
    "csm_orig.smoothness = 200\n",
    "csm_orig.calculate(preprocessed_data)\n",
    "\n",
    "A_orig = mr.AcquisitionModel(preprocessed_data, im_mr)\n",
    "A_orig.set_coil_sensitivity_maps(csm_orig)\n",
    "\n",
    "im_orig = A_orig.backward(preprocessed_data)\n",
    "\n",
    "\n",
    "# Undersampled data\n",
    "csm_new = mr.CoilSensitivityData()\n",
    "csm_new.smoothness = 200\n",
    "csm_new.calculate(acq_new)\n",
    "\n",
    "A_new = mr.AcquisitionModel(acq_new, im_mr)\n",
    "A_new.set_coil_sensitivity_maps(csm_orig)\n",
    "\n",
    "im_new = A_orig.backward(acq_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0e63856",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display it\n",
    "plt.figure();\n",
    "plot_2d_image([1,2,1], numpy.abs(im_orig.as_array())[0,:,:], 'Original image', cmap=\"Greys_r\")\n",
    "plot_2d_image([1,2,2], numpy.abs(im_new.as_array())[0,:,:], 'Undersampled image', cmap=\"Greys_r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d97c4b54",
   "metadata": {},
   "source": [
    "Here we used the coil sensitivity maps from the original data. This is of course not correct, but we should calculate the coil sensitivity maps again from our new `AcquisitionData`. Nevertheless, for a Cartesian sampling scheme, the coil sensitivity maps need to be calculated from the fully sampled k-space centre. For the algorithm to detect which $k_y$ points belong to this fully sampled k-space centre, these points need to be labelled. The MR scanner does this automatically but here we cannot do it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4304239",
   "metadata": {},
   "source": [
    "Now you have got a framework to simulate different undersampling patterns. Now you can try\n",
    "\n",
    "   * different undersmapling factors\n",
    "   * different number of points for the fully sampled central k-space region\n",
    "   * calculate the coil maps from the newly created data and compare them to the original data\n",
    "   * random undersampling and compare it to the regular undersampling\n",
    "   * ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eab06850",
   "metadata": {},
   "source": [
    "In the notebook `cil_joint_tv_mr` we will use what we learned here again, but the code above will be packaged into a function to make it easier to use."
   ]
  }
 ],
 "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": 5
}