notebooks/MR/c_coil_combination.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reconstruction and combination of multi-coil MR data\n",
    "\n",
    "This demonstration shows how to reconstruct MR images on a coil-by-coil basis and how to combine the image information from the different receiver coils\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.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First version: 27th of May 2017  \n",
    "Updated 26nd of June 2021  \n",
    "Author: Christoph Kolbitsch, Johannes Mayer  \n",
    "\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, 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": [
    "%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', 'c_coil_combination')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coil Combination Methods\n",
    "### Goals of this notebook:\n",
    "- Explore ways to combine acquisition data from multiple receiver coils."
   ]
  },
  {
   "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",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coil Combination\n",
    "#### Goals of this notebook:\n",
    "- Explore ways to combine acquisition data from multiple receiver coils."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi Receiver Channel Imaging\n",
    "In principle one does not need multiple coils placed on a patient or phantom in an MR exam. Every scanner has a built-in __body coil__ which is able to receive the MRI signal and one can reconstruct an image from that.\n",
    "\n",
    "Images from the in-built body coil are rarely used diagnostically. Instead what is used are so-called __receiver coils__, aka \"coils\", or somethimes refered to as __phased arrays__. When talking about data recorded by different coils these are also referred to as __channels__ (in the literature it always says something of the kind: \"_... data were acquired using a 32-channel cardiac coil..._\").  \n",
    "\n",
    "This has several advantages, one being that the receiver can be placed very close to the signal source. Others we will discuss in the following! \n",
    "\n",
    "__Important:__ each of this coils also adds a global complex phase onto the recorded signal (hence the name __phased__ arrays!)\n",
    "\n",
    "\n",
    "#### Naming Convention\n",
    "The following expressions are usually used synonymously:\n",
    "- coil sensitivity profile\n",
    "- coil sensitivity maps (CSM)\n",
    "- coil maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% LOAD AND PREPROCESS RAW MR DATA\n",
    "data_path = exercises_data_path('MR', 'PTB_ACRPhantom_GRAPPA')\n",
    "filename = os.path.join(data_path, 'ptb_resolutionphantom_fully_ismrmrd.h5')\n",
    "acq_data = pMR.AcquisitionData(filename)\n",
    "preprocessed_data = pMR.preprocess_acquisition_data(acq_data)\n",
    "preprocessed_data.sort()\n",
    "\n"
   ]
  },
  {
   "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)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% PLOT K-SPACE DATA\n",
    "k_array = k_array / numpy.max(abs(k_array[:]))\n",
    "num_channels = k_array.shape[1]\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",
    "plt.tight_layout()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% APPLY INVERSE FFT TO EACH COIL AND VIEW IMAGES\n",
    "image_array = numpy.zeros(k_array.shape, numpy.complex128)\n",
    "for c in range(k_array.shape[1]):\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",
    "\n",
    "\n",
    "fig = plt.figure()   \n",
    "plt.set_cmap('gray')\n",
    "for c in range(image_array.shape[1]):\n",
    "    ax = fig.add_subplot(2,num_channels//2,c+1)\n",
    "    ax.imshow(abs(image_array[:,c,:]), vmin=0, vmax=0.4)\n",
    "    ax.set_title('Coil '+str(c+1))\n",
    "    ax.axis('off')\n",
    "plt.tight_layout()    \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question:\n",
    "- What differences appear in the individual channel reconstructions compared to the combined image we saw in the last notebook?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sum of Square (SOS) Coil Combination\n",
    "\n",
    "As you can see the individual receiver channels have a spatially varying intensity due to the coil sensitivity profiles. This information needs to be combined.\n",
    "\n",
    "When you have a set of independently reconstructed images for each channel $f_c$ where $c \\in \\{1, \\dots N_c \\}$ labels the individual coil channels.\n",
    "\n",
    "One way to combine the signal from all coil channels is to use a sum-of-squares approach:\n",
    "\n",
    "$$\n",
    "f_{sos} = \\sqrt{ \\sum_c \\bigl{|} \\, f_c \\bigr{|}^2 }\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% COMBINE COIL IMAGES USING SOS\n",
    "#image_array_sos = numpy.sqrt(abs(numpy.sum(image_array,1)))\n",
    "image_array_sos = numpy.sqrt(numpy.sum(numpy.square(numpy.abs(image_array)),1))\n",
    "\n",
    "image_array_sos = image_array_sos/image_array_sos.max()\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.set_cmap('gray')\n",
    "plt.imshow(image_array_sos, vmin=0, vmax=0.7)\n",
    "plt.title('Combined image using sum-of-squares') \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question:\n",
    "\n",
    "- Apart from fancy parallel imaging techniques, why is it even useful to have more than one receiver channel?  \n",
    "- What could be a possible disadvantage of this coil combination approach?\n",
    "- Why is SOS preferable to simply summing them without squaring: $\\;f_{combined} = \\sum_c f_c$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why Coil Sensitivity Maps?\n",
    "There are several reasons why one would need to compute the CSMs, e.g.:\n",
    "1. I want to do parallel imaging (i.e. use spatial encoding provided by CSMs)  \n",
    "2. I want to improve my SNR.\n",
    "\n",
    "\n",
    "##### Weighted Sum (WS) Coil Combination\n",
    "Instead of squaring each channel before adding them up there is another method of combining coils, by weighting each image with it's corresponding CSM.  \n",
    "If the spatial distribution of the coil sensitivity $C_c$ for each coil $c$ is known then combining images as:  \n",
    "\n",
    "$$\n",
    "f_{ws} =  \\frac{1}{{\\sum_{c'}{\\bigl{|} C_{c'} \\bigr{|}^2}}}{\\sum_c C_c^* \\cdot f_c}\n",
    "$$\n",
    "\n",
    "yields an optimal signal-to-noise ratio (SNR).    \n",
    "Note also, that this way of combining channels does not destroy the phase information.\n",
    "\n",
    "__However,__ for each coil one needs to either\n",
    "- estimate the coil map $C_c$ from the data itself.\n",
    "- measure them separately (e.g. in Philips systems).\n",
    "\n",
    "We will focus on computing them from the data in the following.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing Coil Sensitivities\n",
    "The blunt approach is to compute the CSMs using the __SRSS__ (Square Root of the Sum of Squares) approach:\n",
    "\n",
    "$$\n",
    "    C^{SRSS}_c = \\frac{f_c}{f_{sos}}\n",
    "$$\n",
    "\n",
    "and to apply some smoothing to the data.\n",
    "\n",
    "As you can imagine there will be no big SNR difference between $f_{ws}$ and $f_{sos}$ using these coil maps. __We didn't put in any additional effort!__ This works well if the SOS image is homogenous.\n",
    "\n",
    "__This seems a bit pointless!__ True, combining your images this way will __not give you a gigantic SNR gain__, but you __still get a CSM which you can use for parallel imaging__. And you can generate an coil-combined  image __without losing phase information__ (because we smooth the CSMs!)! \n",
    "\n",
    "This all works well when the SOS image is good to begin with. Otherwise, there are more \"sophisticated\" ways to estimate the coilmaps, e.g. methods named _Walsh_ or _Inati_ which lie beyond of the scope of this workshop.  \n",
    "__SIRF already provides this functionality__. For the task of estimating coil sensitivities from the acquisition data, in `pMR` there is a class called `CoilSensitivityData`. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% CALCULATE COIL SENSITIVITIES\n",
    "csm = pMR.CoilSensitivityData()\n",
    "#help(csm)\n",
    "\n",
    "csm.smoothness = 50\n",
    "csm.calculate(preprocessed_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell plotting the Coilmaps\n",
    "csm_array = numpy.squeeze(csm.as_array())\n",
    "\n",
    "# csm_array has orientation [coil, im_x, im_y]\n",
    "csm_array = csm_array.transpose([1,0,2])\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.set_cmap('jet')\n",
    "for c in range(num_channels):\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')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question\n",
    "Please answer the following questions:\n",
    "- Why is there noise in some regions of the coilmaps and outside the object? \n",
    "- Is this noise in the coilmap going to have a strong negative impact on the image quality in this region of the combined image?\n",
    "- In which organ in the human anatomy would you expect a coilmap to look similarly noisy?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% COMBINE COIL IMAGES USING WEIGHTED SUM\n",
    "image_array_ws = numpy.sum(numpy.multiply(image_array, numpy.conj(csm_array)),1)\n",
    "image_array_ws = abs(numpy.divide(image_array_ws, numpy.sum(numpy.multiply(csm_array, numpy.conj(csm_array)),1)))\n",
    "image_array_ws = image_array_ws/image_array_ws.max()\n",
    "\n",
    "\n",
    "diff_img_arr = abs(image_array_sos-image_array_ws)\n",
    "diff_img_arr = diff_img_arr/diff_img_arr.max()\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=[12, 4])\n",
    "plt.set_cmap('gray')\n",
    "\n",
    "ax = fig.add_subplot(1,3,1)\n",
    "ax.imshow(image_array_sos, vmin=0, vmax=0.7)\n",
    "ax.set_title('Sum-of-squares (SOS)')\n",
    "ax.axis('off')\n",
    "\n",
    "ax = fig.add_subplot(1,3,2)\n",
    "ax.imshow(image_array_ws, vmin=0, vmax=0.7)\n",
    "ax.set_title('Weighted sum (WS)')\n",
    "ax.axis('off')\n",
    "\n",
    "ax = fig.add_subplot(1,3,3)\n",
    "ax.imshow(diff_img_arr, vmin=-0, vmax=0.1)\n",
    "ax.set_title('SOS - WS')\n",
    "ax.axis('off')\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Image Quality Assessment\n",
    "In low-signal regions you can see from the difference image, a weighted coil combination will give you an improved SNR.  \n",
    "This dataset was acquired with a head-coil so which is very well-matched, and no flow artifacts appear so the difference is not huge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question: What's \"The mysterious next step?\":\n",
    "1. Reconstructing individual channels, add them up and get one image out.\n",
    "2. Reconstructing individual channels, weight them by their coil maps, add them up and get one image out.\n",
    "3. Do the mysterious next step, and get one image out. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recap Coil Combination\n",
    "\n",
    "In this exercise we learned:\n",
    "- how to combine multichannel image reconstructions\n",
    "- how to compute a simple coil sensitivity from our data.\n",
    "- how to employ SIRF to execute this computation for us.\n",
    "\n",
    "\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": 1
}