mri_acq_recon.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage import io\n",
    "from skimage.restoration import denoise_bilateral\n",
    "from scipy.io import loadmat\n",
    "\n",
    "\n",
    "def show_image(image, title='', is_fourier=False):\n",
    "    \"\"\"Show image with matplotlib\"\"\"\n",
    "    plt.figure()\n",
    "    plt.title(title)\n",
    "    if is_fourier:\n",
    "        plt.imshow(np.abs(image), cmap='gray', vmax=np.abs(image).max()*0.002)\n",
    "    else:\n",
    "        plt.imshow(np.abs(image), cmap='gray')\n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A bit of context\n",
    "MRI stands for Magnetic Resonance Imaging.\n",
    "It is a medical imaging modality that is used to probe soft tissues inside the human body non-invasively and without the use of ionizing radiations.\n",
    "\n",
    "The MRI data is acquired in the Fourier space of the image. Tradiailly, this acquisition domain is also called \"k\"-space, with \"k\" referring to spatial frequency.\n",
    "However, MRI is particularly slow because it relies on the Magnetic Resonance phenomenon, which is inherently slow.\n",
    "One way to speed up MRI is to acquire fewer signals (k-space shots ${\\bf y}$), i.e. undersample.\n",
    "If we do that, we find ourselves in a situation where we have not enough signals to reconstruct the image using simply the Inverse Fourier Transform ($F^{-1}$) as is typically done.\n",
    "In fact, we have an ill-posed inverse problem that one must solve using priors on the MR images ${\\bf x}$.\n",
    "\n",
    "$$\n",
    "F_{\\Omega}{\\bf x} = {\\bf y}\n",
    "$$\n",
    "\n",
    "where $F_{\\Omega}$ is a masked Fourier Transform with a mask $\\Omega$ that selects the sample locations in the k-space.\n",
    "\n",
    "This means that we need a function to tell apart all the possible solutions of the inverse problem, that tells which solution is more likely to be an MR image.\n",
    "\n",
    "In classical Compressed Sensing, these priors are handcrafted.\n",
    "Deep Learning allows us to learn these priors from data.\n",
    "\n",
    "In this exercise, we will load a basic synthetic MR image, understand the corresponding Fourier space (acquisition domain) and later undersample the signal. Finally, we reconstruct the image and denoise it."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load an Synthetic MRI image for our use\n",
    "Here we simply pull up an already saved synthetic brain MR image we have saved in the repository."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dirimg = os.path.join(os.getcwd(), \"data\")\n",
    "img_size = 512  \n",
    "FOV = 0.2 \n",
    "pixelSize = FOV/img_size\n",
    "\n",
    "#load data file corresponding to the target resolution\n",
    "filename = \"BrainPhantom\" + str(img_size) + \".png\"\n",
    "mri_filename = os.path.join(dirimg, filename)\n",
    "mri_img = io.imread(mri_filename, as_gray=True)\n",
    "show_image(mri_img, \"Synthetic MRI image\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MRI Acquisition\n",
    "Traditional MR acquisition takes place in the fourier domain of the image. \n",
    "In the most simplistic example, to acquire the above image, we sample its corresponding Fourier coefficients.\n",
    "Then the image is reconstructed by transforming this acquired coeffcients back to image domain.\n",
    "Note that the Fourier Coefficients is in complex domain.\n",
    "\n",
    "\n",
    "Below we write simple interfaces to transfer the image to fourier coefficients and back to image domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def image_to_coeff(image):\n",
    "    \"\"\"Convert image to Fourier coefficients\"\"\"\n",
    "    return np.fft.fftshift(np.fft.fft2(image))\n",
    "\n",
    "def coeff_to_image(coeff):\n",
    "    \"\"\"Convert Fourier coefficients to image\"\"\"\n",
    "    return np.fft.ifft2(np.fft.ifftshift(coeff))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualize the Fourier coefficients\n",
    "**TODO:** Can you please write the code to visualize the Fourier coefficients of the image? Try to use above defined functions and `show_image`.\n",
    "\n",
    "\n",
    "**PROTIP:** Maybe pass `is_fourier=True` to better visualize."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "<TODO: compute the Fourier coefficients of the image and display the result>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Undersampling the acquistion\n",
    "Now as the MRI scans are very slow, we traditionally undersample the image. \n",
    "\n",
    "**TODO:** Based on the above coefficients, what portions of the coefficients are most important? \n",
    "\n",
    "**Answer**: "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Undersampling mask\n",
    "Here we pull up an already generated mask. The mask is a binary image with 1s where we want to sample and 0s where we don't want to sample. We will use this mask to sample the Fourier coefficients of the image.\n",
    "\n",
    "**TODO**: Can you please get the fourier coefficients and sample is based on the mask?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "densities = loadmat(os.path.join(dirimg, \"2d_sampling_densities.mat\"))['distrib2d_N512_sym10']\n",
    "fourier_mask = densities > 15 * densities.min() \n",
    "\n",
    "show_image(fourier_mask, 'Sampling mask')\n",
    "\n",
    "\n",
    "# 1. Apply the mask to the Fourier coefficients of the image\n",
    "# 2. Visualize the sampled data\n",
    "sampled_coeff = <TODO: apply mask to Fourier coefficients>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modeling noise\n",
    "Generally, the acquired data at the scanner has noise in it. This noise is usually additive white gaussian noise. \n",
    "\n",
    "**TODO**: Can you please add this noise into the sampled data above. Use `np.random.random()` to generate random gaussian noise. Be careful, add noise only on sampled fourier coefficients, others must remain 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noisy_data = <TODO: Your code here>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstruct the image\n",
    "**TODO**: Can you please take the noisy samples back to image domain using `coeff_to_image` and visualize the result?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reconstructed_image = <TODO: Your code here>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Denoising the final reconstructed image.\n",
    "\n",
    "The reconstructed image above is noisy due to 1) sample noise we simulated, and also as we 2) undersampled the acquisition.\n",
    "\n",
    "**TODO**: Now we would like to clean this image, can you use Bilateral filter from `scikit-image` library to denoise the image?\n",
    "Have a look at the documentation of `denoise_bilatral`.\n",
    "\n",
    "**NOTE**: Bilateral filters work on real images, while reconstructed image is complex, maybe pull up only the magnitude image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cleanup_image = <TODO: Add code to renoise the image>\n",
    "show_image(cleanup_image, 'Cleanup image')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.12 64-bit ('p310')",
   "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.10.12"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "ecdd29a57db5b5201a11872982822c462aa68cfb0b0bb3165c3351781cd5f076"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}