FLASH_simulator.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. Rewindow button seems to make image very bright\n",
    "# 2. move slider with mouse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import scipy.io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.widgets import Slider, Button"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flash simulator\n",
    "\n",
    "This simulator loads T1, T2 and PD parameter maps, visualizes them, and then runs the simulation to construct an estimated MR image based on the users input of TR, TE and Flip Angle (FA) for a chosen slice.\n",
    "\n",
    "The simulator implements the following equation:\n",
    "\n",
    "$$\n",
    "     Mxy = \\text{PD}*\\frac{(1-e^{-\\text{TR}/\\text{T1}})}{(1-\\cos{(FA)}*e^{-TR/T1})}* \\sin{(FA)}*e^{-\\text{TE}/\\text{T2}}\n",
    "$$\n",
    "<br>\n",
    "*Code History:*<br>\n",
    "- *CP:  October 2021* <br>\n",
    "&emsp; First Python version <br><br>\n",
    "- *PJV: v.1.4: November 2017* <br>\n",
    "   &emsp; Changes: - It loads the parameter maps directly, rather than computing them from the original images. It makes it                less confusing. <br>\n",
    "   &emsp;&emsp;&emsp;&emsp;&emsp; $\\;\\;\\;$ - If the script fails to open any of the image files it returns an error asking the user to cd to the                  mriTools/mriTutorial folder.<br><br>\n",
    "- *PJV: v.1.3: December 2015* <br>\n",
    "   &emsp; Minor change: - limit the minimum FA to 10 deg<br><br>\n",
    "- *PJV: v.1.2: November 2013*<br>\n",
    "    &emsp;Minor changes: - limit the minimum TR to TE+1 ms.<br>\n",
    "    &emsp;&emsp;&emsp;&emsp;&emsp; $\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;$- increase the minimum TR and TE to 1 ms.<br>\n",
    "    &emsp;&emsp;&emsp;&emsp;&emsp; $\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;$- improve functionality of the \"window\" button.<br><br>\n",
    "- PJV: v.1.1: August 2011 <br>\n",
    "    &emsp;Minor changes: - use \"single\" variables, for better memory performance.<br>\n",
    "    &emsp;&emsp;&emsp;&emsp;&emsp; $\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;$- only update the image when there have been changes in the parameter values.<br> <br>\n",
    "- *PJV: v.1.0: March 2009* <br>\n",
    "     &emsp;Created from MR_simulator.m (Written by G.M. Boynton for Psych 555C, Spring 2007 at U. Washington) Tissue parameters are based on the McGill &emsp;brain database (http://www.bic.mni.mcgill.ca/brainweb/tissue_mr_parameters.txt)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Load parameter maps and show T1, T2 and PD images in a figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Path for data\n",
    "mri_data_path = os.getcwd() + '/data/mri_data/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load T1, T2, and PD\n",
    "\n",
    "T1 = scipy.io.loadmat(mri_data_path+'T1.mat')['T1']\n",
    "T2 = scipy.io.loadmat(mri_data_path+'T2.mat')['T2']\n",
    "PD = scipy.io.loadmat(mri_data_path+'PD.mat')['PD']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(181, 217, 181)\n",
      "(181, 217, 181)\n",
      "(181, 217, 181)\n"
     ]
    }
   ],
   "source": [
    "# Check their sizes\n",
    "print(T1.shape)\n",
    "print(T2.shape)\n",
    "print(PD.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's pick an axial slice and show it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = 99\n",
    "\n",
    "#Plot it\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=[20,10])\n",
    "fig.subplots_adjust(wspace=0.04)\n",
    "\n",
    "#T1\n",
    "ax0 = axes[0].imshow(np.flipud(np.transpose(PD[:,:,z])/np.transpose(T1[:,:,z])), cmap='Greys_r')\n",
    "ax0.set_clim(0, 0.0015) \n",
    "axes[0].set_title('T1', fontsize = 20, fontweight='bold', fontname='Arial')\n",
    "axes[0].axis('off')\n",
    "\n",
    "#T2\n",
    "ax1 = axes[1].imshow(np.flipud(np.transpose(T2[:,:,z])), cmap='Greys_r')\n",
    "ax1.set_clim(50, 120) \n",
    "axes[1].set_title('T2', fontsize = 20, fontweight='bold', fontname='Arial')\n",
    "axes[1].axis('off')\n",
    "\n",
    "#PD\n",
    "ax2 = axes[2].imshow(np.flipud(np.transpose(PD[:,:,z])), cmap='Greys_r')\n",
    "axes[2].set_title('PD', fontsize = 20, fontweight='bold', fontname='Arial')\n",
    "axes[2].axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Set up the simulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The parametrized function to be plotted\n",
    "def f(PD, TR, T1, FA, TE, T2, z):\n",
    "    Mxy = PD[:,:,z]*(1-np.exp(-TR/T1[:,:,z]))/(1-np.cos(np.deg2rad(FA))*np.exp(-TR/T1[:,:,z]))* np.sin(np.deg2rad(FA))*np.exp(-TE/T2[:,:,z])\n",
    "    Mxy = np.flipud(np.transpose(Mxy))\n",
    "    return 64*Mxy/np.max(Mxy[:])\n",
    "\n",
    "#Initial parameters:\n",
    "TR = 500  #in msec\n",
    "TE = 15   #in msec\n",
    "FA = 90 #in degrees\n",
    "z = 99\n",
    "\n",
    "# Create the figure \n",
    "fig, ax = plt.subplots(figsize = (7,7))\n",
    "simulator = ax.imshow(f(PD, TR, T1, FA, TE, T2, z), cmap='Greys_r')\n",
    "ax.axis('off')\n",
    "axcolor = 'silver'\n",
    "ax.margins(x=0)\n",
    "# adjust the main plot to make room for the sliders\n",
    "plt.subplots_adjust(left=0.1, bottom=0.2)\n",
    "\n",
    "# Make a horizontal slider to control the FA.\n",
    "ax_FA = plt.axes([0.3, 0.03, 0.4, 0.03], facecolor=axcolor)\n",
    "FA_slider = Slider(\n",
    "    ax=ax_FA,\n",
    "    label='FA    ',\n",
    "    valmin=10,\n",
    "    valmax=180,\n",
    "    valstep = 1,\n",
    "    valinit=FA,\n",
    "    valfmt='    %i'\n",
    ")\n",
    "\n",
    "\n",
    "# Make a horizontal slider to control the TR.\n",
    "ax_TR = plt.axes([0.3, 0.07, 0.4, 0.03], facecolor=axcolor)\n",
    "TR_slider = Slider(\n",
    "    ax=ax_TR,\n",
    "    label='TR    ',\n",
    "    valmin=1,\n",
    "    valmax=2000,\n",
    "    valstep = 1,\n",
    "    valinit=TR,\n",
    "    valfmt='    %i'\n",
    ")\n",
    "\n",
    "# Make a horizontal slider to control the TE.\n",
    "ax_TE = plt.axes([0.3, 0.11, 0.4, 0.03], facecolor=axcolor)\n",
    "TE_slider = Slider(\n",
    "    ax=ax_TE,\n",
    "    label='TE    ',\n",
    "    valmin=1,\n",
    "    valmax=100,\n",
    "    valstep = 1,\n",
    "    valinit=TE,\n",
    "    valfmt='    %i'\n",
    ")\n",
    "\n",
    "# Make a horizontal slider to control the slice z.\n",
    "ax_slice = plt.axes([0.3, 0.15, 0.4, 0.03], facecolor=axcolor)\n",
    "slice_slider = Slider(\n",
    "    ax=ax_slice,\n",
    "    label='Slice    ',\n",
    "    valmin=0,\n",
    "    valmax=180,\n",
    "    valstep = 1,\n",
    "    valinit=z,\n",
    "    valfmt='    %i'\n",
    ")\n",
    "\n",
    "\n",
    "# The function to be called anytime a slider's value changes\n",
    "def update(val):\n",
    "    if int(TR_slider.val) < int(TE_slider.val):\n",
    "        TR_slider.set_val(int(TE_slider.val)+1)\n",
    "    simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    fig.canvas.draw_idle()\n",
    "    \n",
    "    \n",
    "# Register the update function with each slider\n",
    "TR_slider.on_changed(update)\n",
    "TE_slider.on_changed(update)\n",
    "slice_slider.on_changed(update)\n",
    "FA_slider.on_changed(update)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to increase the FA slider value (by 2 each time)\n",
    "FA_right_button_ax = plt.axes([0.7, 0.03, 0.014, 0.03])\n",
    "FA_right_button = Button(FA_right_button_ax, '>', color=axcolor, hovercolor='0.975')\n",
    "def increase_FA(event):\n",
    "    if FA_slider.val in range(10,179):\n",
    "        FA_slider.set_val(FA_slider.val+2)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    elif FA_slider.val ==179:\n",
    "        FA_slider.set_val(180)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, 180, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    fig.canvas.draw_idle()\n",
    "FA_right_button.on_clicked(increase_FA)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to decrease the FA slider value (by 2)\n",
    "FA_left_button_ax = plt.axes([0.286, 0.03, 0.014, 0.03])\n",
    "FA_left_button = Button(FA_left_button_ax, '<', color=axcolor, hovercolor='0.975')\n",
    "def decrease_FA(event):\n",
    "    if FA_slider.val in range(12,181):\n",
    "        FA_slider.set_val(FA_slider.val-2)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    elif FA_slider.val==11:\n",
    "        FA_slider.set_val(10)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, 10, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    fig.canvas.draw_idle()\n",
    "FA_left_button.on_clicked(decrease_FA)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to increase the TR slider value (by 20 each time)\n",
    "TR_right_button_ax = plt.axes([0.7, 0.07, 0.014, 0.03])\n",
    "TR_right_button = Button(TR_right_button_ax, '>', color=axcolor, hovercolor='0.975')\n",
    "def increase_TR(event):\n",
    "    if TR_slider.val in range(1,1981):\n",
    "        TR_slider.set_val(TR_slider.val+20)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    elif TR_slider.val in range(1981,2000):\n",
    "        TR_slider.set_val(2000)\n",
    "        simulator.set_data(f(PD, 2000, T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    fig.canvas.draw_idle()\n",
    "TR_right_button.on_clicked(increase_TR)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to decrease the TR slider value (by 20)\n",
    "TR_left_button_ax = plt.axes([0.286, 0.07, 0.014, 0.03])\n",
    "TR_left_button = Button(TR_left_button_ax, '<', color=axcolor, hovercolor='0.975')\n",
    "def decrease_TR(event):\n",
    "    if int(TR_slider.val) > int(TE_slider.val)+1:\n",
    "        if TR_slider.val in range(21,2001):\n",
    "            TR_slider.set_val(TR_slider.val-20)\n",
    "            simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "        elif TR_slider.val in range(1,21):\n",
    "            TR_slider.set_val(1)\n",
    "            simulator.set_data(f(PD, 1, T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    fig.canvas.draw_idle()\n",
    "TR_left_button.on_clicked(decrease_TR)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to increase the TE slider value (by 1)\n",
    "TE_right_button_ax = plt.axes([0.7, 0.11, 0.014, 0.03])\n",
    "TE_right_button = Button(TE_right_button_ax, '>', color=axcolor, hovercolor='0.975')\n",
    "def increase_TE(event):\n",
    "    if TE_slider.val in range(1,100):\n",
    "        TE_slider.set_val(TE_slider.val+1)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    fig.canvas.draw_idle()\n",
    "TE_right_button.on_clicked(increase_TE)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to decrease the TE slider value (by 1)\n",
    "TE_left_button_ax = plt.axes([0.286, 0.11, 0.014, 0.03])\n",
    "TE_left_button = Button(TE_left_button_ax, '<', color=axcolor, hovercolor='0.975')\n",
    "def decrease_TE(event):\n",
    "    if TE_slider.val in range(2,101):\n",
    "        TE_slider.set_val(TE_slider.val-1)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    fig.canvas.draw_idle()\n",
    "TE_left_button.on_clicked(decrease_TE)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to increase the slice slider value (by 2)\n",
    "slice_right_button_ax = plt.axes([0.7, 0.15, 0.014, 0.03])\n",
    "slice_right_button = Button(slice_right_button_ax, '>', color=axcolor, hovercolor='0.975')\n",
    "def increase_slice(event):\n",
    "    if slice_slider.val in range(1,179):\n",
    "        slice_slider.set_val(slice_slider.val+2)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)))\n",
    "    elif slice_slider.val in range(179,180):\n",
    "        slice_slider.set_val(180)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, 180))\n",
    "    fig.canvas.draw_idle()\n",
    "slice_right_button.on_clicked(increase_slice)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to decrease the slice slider value (by 2)\n",
    "slice_left_button_ax = plt.axes([0.286, 0.15, 0.014, 0.03])\n",
    "slice_left_button = Button(slice_left_button_ax, '<', color=axcolor, hovercolor='0.975')\n",
    "def decrease_slice(event):\n",
    "    if slice_slider.val in range(3,181):\n",
    "        slice_slider.set_val(slice_slider.val-2)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, int(slice_slider.val)-2))\n",
    "    elif slice_slider.val in range(1,3):\n",
    "        slice_slider.set_val(1)\n",
    "        simulator.set_data(f(PD, int(TR_slider.val), T1, FA_slider.val, int(TE_slider.val), T2, 1))\n",
    "    fig.canvas.draw_idle()\n",
    "    fig.canvas.draw_idle()\n",
    "slice_left_button.on_clicked(decrease_slice)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to reset the sliders to initial values.\n",
    "resetax = plt.axes([0.8, 0.025, 0.1, 0.04])\n",
    "button = Button(resetax, 'Reset', color=axcolor, hovercolor='0.975')\n",
    "def reset(event):\n",
    "    TR_slider.reset()\n",
    "    TE_slider.reset()\n",
    "    slice_slider.reset()    \n",
    "button.on_clicked(reset)\n",
    "\n",
    "# Create a `matplotlib.widgets.Button` to re-window the image.\n",
    "rewindowax = plt.axes([0.8, 0.1, 0.1, 0.04])\n",
    "rewindow_button = Button(rewindowax, 'Re-window', color=axcolor, hovercolor='0.975')\n",
    "def rewindow(event):\n",
    "    Mxy = f(PD, int(TR_slider.val), T1, int(FA_slider.val), int(TE_slider.val), T2, int(slice_slider.val))\n",
    "    Mxy_sorted = np.sort(Mxy[Mxy[:]>0]) # don't take zeros into account\n",
    "    low_5_pct = Mxy_sorted[int(np.round(0.05*len(Mxy_sorted)))] # low 5% level\n",
    "    bkgr = np.mean(Mxy_sorted[Mxy_sorted<=low_5_pct], axis=0) # mean background level\n",
    "    Mxy_no_bkgr = Mxy_sorted[Mxy_sorted>=10*bkgr]      #cut-off at 10x bkgr\n",
    "    window_low  = Mxy_no_bkgr[int(round(0.10*len(Mxy_no_bkgr)))]  # low 5% level\n",
    "    window_high = Mxy_no_bkgr[int(round(0.95*len(Mxy_no_bkgr)))]  # high 5% level\n",
    "    Mxy = 64*(Mxy-window_low)/(window_high-window_low)\n",
    "    simulator.set_data(Mxy)\n",
    "    fig.canvas.draw_idle()\n",
    "rewindow_button.on_clicked(rewindow)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}