3D-ISMRMRD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<h3>  1️⃣ Know your data (at least) <strong style=\"color:darkgreen\"><u>at the level you are dealing with it</u></strong> </h3></div>\n",
    "\n",
    "### NOt RApid Honestly JOyful NEvertheleSs (NORAH JONES) Readout\n",
    "\n",
    "<br>\n",
    "<center><code style=\"background:yellow;color:black\"> We are looking at the readout block of our NORAH JONES sequence <code>/Sequence/CartesianReadout3D.spv</code>. If you are an OSX user, you can inspect sequence components using <a href=\"https://www.heartvista.ai/spinbench\">SpinBench!</a></code></center>\n",
    "<br>\n",
    "\n",
    "![](assets/sequence_order.gif)\n",
    "\n",
    "\n",
    "We have **100 `acquisitions`** (or `TR`s or `Phase Encoding`s) that are ordered outside-in to sample the center of the k-space the latest (steady-state). At each **`acquisition`** we have **100 `samples`**. NORAH JONES sequence acquires all 4 images in this fashion, at the same location. You can see all the readout parameters below:\n",
    "\n",
    "\n",
    "![](assets/readout_params.png)\n",
    "\n",
    "### NORAH JONES Raw Data in ISMRM-RD Format\n",
    "\n",
    "![](assets/ismrmrd_structure.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<h3> 2️⃣ Mise en place </h3></div>\n",
    "\n",
    "Let's read and parse the data using `NumPy` relying purely on metadata provided by `python-ismrmrd`.\n",
    "\n",
    "### Import ISMRM-RD modules 🧲"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ismrmrd import Dataset as read_ismrmrd\n",
    "from ismrmrd.xsd import CreateFromDocument as parse_ismrmd_header\n",
    "import numpy as np\n",
    "import plotly.express as px\n",
    "import plotly.graph_objects as go"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a simple list storing file names 🗂"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = list()\n",
    "for chord in ['chord1','chord2','chord3','chord4']:\n",
    "    files.append('./KSpace/sub-ismrm_ses-sunrise_acq-' + chord + '.h5')\n",
    "files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Explore `acq-chord1` file 🕵️‍♀️\n",
    "\n",
    "#### Read data & parse header"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the first file to explore data\n",
    "dset = read_ismrmrd(files[0], 'dataset')\n",
    "\n",
    "# Parse header \n",
    "header = parse_ismrmd_header(dset.read_xml_header())\n",
    "\n",
    "# See what's in the header\n",
    "print('header has:',list(vars(header).keys()))\n",
    "print('----')\n",
    "print('header.encoding[0] has:',list(vars(header.encoding[0]).keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Grab some useful metadata values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nX = header.encoding[0].encodedSpace.matrixSize.x\n",
    "nY = header.encoding[0].encodedSpace.matrixSize.y\n",
    "nZ = header.encoding[0].encodedSpace.matrixSize.z\n",
    "nCoils = header.acquisitionSystemInformation.receiverChannels\n",
    "print('nX: ', nX,'nY: ', nY,'nZ: ', nZ, 'nCoils: ', nCoils)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Initialize a 3D `NumPy` matrix to store multi-channel raw data 🗃\n",
    "\n",
    "We already know the dimensions! Given that we expect complex data, we can set `dtype=complex64` 😎 For simplicity I will omit nZ (we have one slice). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw = np.zeros((nCoils, nX, nY), dtype=np.complex64)\n",
    "print('Raw data shape: ', raw.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loop over `acquisitions` to fill out `raw` \n",
    "\n",
    "The `dset` object we created from the ISMRM-RD `Dataset` (we used alias `read_ismrmrd`) class has an important method: **`read_acquisition`**. We will loop over `nY` to read k-space stripes (cartesian) one by one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for tr in range(nY):\n",
    "    raw[:,:,tr] = dset.read_acquisition(tr).data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Is not `Plotly` simply amazing? \n",
    "This feature has been introduced in `v4.14.0`: with a single line of code you can animate or grid display 3D data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.imshow(raw.real,animation_frame=0, color_continuous_scale='viridis', labels=dict(animation_frame=\"Channel\"),template='plotly_dark')\n",
    "fig.update_layout(title='Channel Raw')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "<h3> 3️⃣ Time to perform a modest reconstruction :) </h3></div>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.fft import fft2, fftshift\n",
    "from scipy import ndimage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loop over channels and `fft2`\n",
    "\n",
    "Tip: remove `fftshift` from `fourier_ellipsoid` function call, set filter size > `20`, then play the animation along with your remix from the first notebook 🎵🎶"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im = np.zeros(raw.shape)\n",
    "\n",
    "# Let's apply some ellipsoid filter. \n",
    "raw = ndimage.fourier_ellipsoid(fftshift(raw),size=2)\n",
    "#raw = ndimage.fourier_ellipsoid(raw,size=2)\n",
    "\n",
    "for ch in range(nCoils):\n",
    "    # Comment in and see what it gives \n",
    "    im[ch,:,:] = abs(fftshift(fft2(raw[ch,:,:])))\n",
    "    # Normalize \n",
    "    im[ch,:,:] /= im[ch,:,:].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.imshow(im,color_continuous_scale='viridis', animation_frame=0,template='plotly_dark')\n",
    "fig.update_layout(title='Channel Recon').show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im_sos = np.sqrt(np.sum(im*np.conj(im),axis=0))\n",
    "# Normalize\n",
    "im_sos /= im_sos.max()\n",
    "\n",
    "# Crop a bit \n",
    "fig = px.imshow(im_sos[10:90,0:80],template='plotly_dark',color_continuous_scale='viridis')\n",
    "fig.update_layout(title='Sum of Squares Combine')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "<h3>  Exercises </h3></div>\n",
    "\n",
    "* Rotate/flip images\n",
    "* Reconstruct `SoS` images for every TR by writing some functions and collect them into a 3D matrix.\n",
    "    * Display voxel intensity changes over TRs at the center of each (14) sphere "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}