mp2rage_blog/MP2RAGE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS",
    "tags": [
     "scratch"
    ]
   },
   "source": [
    "## Welcome to a qMRLab interactive blog post Jupyter Notebook!\n",
    "\n",
    "If this is your first time running a Juptyer Notebook, there's a lot of tutorials available online to help. [Here's one](https://www.dataquest.io/blog/jupyter-notebook-tutorial/) for your convenience.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This notebook contains everything needed to reproduce the MP2RAGE T<sub>1</sub> blog post on the [qMRLab website](). In fact, this notebook generated the HTML for the blog post too! This notebook is currently running on a MyBinder server that only you can access, but if you want to be kept up-to-date on any changes that the developpers make to this notebook, you should go to it's [GitHub repository](https://github.com/qMRLab/t1_notebooks) and follow it by clicking the \"Watch\" button in the top right (you may need to create a GitHub account, if you don't have one already).\n",
    "\n",
    "## Tips\n",
    "\n",
    "Here's a few things you can do in this notebook\n",
    "\n",
    "### Code\n",
    "* Run the entire processing by clicking above on the \"Kernel\" tab, then \"Restart & Run All\". It will be complete when none of the cells have an asterix \"\\*\" in the square brackets.\n",
    "* To change the code, you need to click once on code cells. To re-run that cell, click the \"Run\" button above when the cell is selected.\n",
    "  * **Note:** Cells can depend on previous cells, or even on previous runs of the cell itself, so it's best to run all the previous cells beforehand.\n",
    "* This binder runs on SoS, which allows the mixing of Octave (i.e. an open-source MATLAB) and Python cells. Take a look a the drop down menu on the top right of the cells to know which one you are running.\n",
    "* To transfer data from cells of one language to another, you need to create a new cell in the incoming language and run `%get (param name) --from (outgoing language)`. See cells below for several examples within this notebook.\n",
    "\n",
    "### HTML\n",
    "* To reproduce the HTML of the blog post, run the entire processing pipeline (see point one in the previous section), then save the notebook (save icon, top left). Now, click on the drop down menu on the left pannel, and select `%sossave --to html --force` . After a few seconds, it should output \"Workflow saved to VariableFlipAngle.html\" – click on the HTML name, and you're done!\n",
    "* Cells with tags called \"scratch\" are not displayed in the generated HTML.\n",
    "* Cells with the tag \"report_output\" display the output (e.g. figures) in the generated HTML.\n",
    "* Currently in an un-run notebook, the HTML is not formatted like the website. To do so, run the Python module import cell (`# Module imports`) and then very last cell (`display(HTML(...)`).\n",
    "\n",
    "**If you have any other questions or comments, please raise them in a [GitHub issue](https://github.com/qMRLab/t1_notebooks/issues).**\n",
    "\n",
    "# Note\n",
    "\n",
    "The following cell is meant to be displayed for instructional purposes in the blog post HTML when \"All cells\" gets displayed (i.e. the Octave code)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Octave"
   },
   "outputs": [],
   "source": [
    "% **Blog post code introduction**\n",
    "% \n",
    "% Congrats on activating the \"All cells\" option in this interactive blog post =D\n",
    "%\n",
    "% Below, several new HTML blocks have appears prior to the figures, displaying the Octave/MATLAB code that was used to generate the figures in this blog post.\n",
    "%\n",
    "% If you want to reproduce the data on your own local computer, you simply need to have qMRLab installed in your Octave/MATLAB path and run the \"startup.m\" file, as is shown below.\n",
    "%\n",
    "% If you want to get under the hood and modify the code right now, you can do so in the Jupyter Notebook of this blog post hosted on MyBinder. The link to it is in the introduction above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Python3",
    "tags": [
     "scratch"
    ]
   },
   "outputs": [],
   "source": [
    "# PYTHON CODE\n",
    "# Module imports\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import plotly.plotly as py\n",
    "import plotly.graph_objs as go\n",
    "import numpy as np\n",
    "from plotly import __version__\n",
    "from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\n",
    "config={'showLink': False, 'displayModeBar': False}\n",
    "\n",
    "init_notebook_mode(connected=True)\n",
    "\n",
    "from IPython.core.display import display, HTML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "<center><h1 style=\"font-family: timesnewroman;font-size: 40px;\">MP2RAGE T<sub>1</sub> Mapping</h1></center>\n",
    "<p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "<div class=blog_body>\n",
    "<p style=\"text-align:justify;\">\n",
    "Dictionary-based MRI techniques capable of generating T<sub>1</sub> maps are increasing in popularity, due to their growing availability on clinical scanners, rapid scan times, and fast post-processing computation time, thus making quantitative T<sub>1</sub> mapping accessible for clinical applications. Generally speaking, dictionary-based quantitative MRI techniques use numerical dictionaries—databases of pre-calculated signal values simulated for a wide range of tissue and protocol combinations—during the image reconstruction or post-processing stages. Popular examples of dictionary-based techniques that have been applied to T<sub>1</sub> mapping are MR Fingerprinting (MRF) (Ma et al. 2013), certain flavours of compressed sensing (CS) (Doneva et al. 2010; Li et al. 2012), and Magnetization Prepared 2 Rapid Acquisition Gradient Echoes (MP2RAGE) (Marques et al. 2010). Dictionary-based techniques can usually be classified into one of two categories: techniques that use information redundancy from parametric data to assist in accelerated imaging (e.g. CS, MRF), or those that use dictionaries to estimate quantitative maps using the MR images after reconstruction. Because MP2RAGE is a technique implemented primarily for T<sub>1</sub> mapping, and it is becoming increasingly available as a standard pulse sequence on many MRI systems, the remainder of this section will focus solely on this technique. However, many concepts discussed are shared by other dictionary-based techniques.\n",
    "\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "MP2RAGE is an extension of the conventional MPRAGE pulse sequence widely used in clinical studies (Haase et al. 1989; Mugler & Brookeman 1990). A simplified version of the MP2RAGE pulse sequence is shown in Figure 1. MP2RAGE can be seen as a hybrid between the inversion recovery and VFA pulse sequences: a 180° inversion pulse is used to prepare the magnetization for T<sub>1</sub> sensitivity at the beginning of each TRMP2RAGE, and then two images are acquired at different inversion times using gradient recalled echo (GRE) imaging blocks with low flip angles and short repetition times (TR). During a given GRE imaging block, each excitation pulse is followed by a constant in-plane (“y”) phase encode weighting (varied for each TRMP2RAGE), but with different 3D (“z”) phase encoding gradients (varied at each TR). The center of k-space for the 3D phase encoding direction is acquired at the TI for each GRE imaging block. The main motivation for developing the MP2RAGE pulse sequence was to provide a metric similar to MPRAGE, but with self-bias correction of the static (B<sub>0</sub>) and receive (B<sub>1</sub><sup>-</sup>) magnetic fields, and a first order correction of the transmit magnetic field (B<sub>1</sub><sup>+</sup>). However, because two images at different TIs are acquired (unlike MPRAGE, which only acquires data at a single TI), information about the T<sub>1</sub> values can also be inferred, thus making it possible to generate quantitative T<sub>1</sub> maps using this data.\n",
    "</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "<div class=figure_caption>\n",
    "<p style=\"text-align:justify;\">\n",
    "<b>\n",
    "Figure 1. Simplified diagram of an MP2RAGE pulse sequence. TR: repetition time between successive gradient echo readouts, TR<sub>MP2RAGE</sub>: repetition time between successive adiabatic 180° inversion pulses, TI<sub>1</sub> and TI<sub>2</sub>: inversion times, <i>θ<sub>1</sub></i> and <i>θ<sub>2</sub></i>: excitation flip angles. The imaging readout events occur within each TR using a constant in-plane phase encode (“y”) gradient set for each TR<sub>MP2RAGE</sub>, but varying 3D phase encode (“z”) gradients between each successive TR.\n",
    "</b>\n",
    "</p>\n",
    "</div>\n",
    "\n",
    "<p>\n",
    "<center><img src=\"mp2rage_pulsesequence.png\" style=\"width:800px;height:auto;\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "<center> <h2 style=\"font-family:timesnewroman;font-size:30px\">Signal Modelling</h2> </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "SoS"
   },
   "source": [
    "<div class=blog_body>\n",
    "<p style=\"text-align:justify;\">\n",
    "Prior to considering the full signal equations, we will first introduce the equation for the MP2RAGE parameter (<i>S</i><sub>MP2RAGE</sub>) that is calculated in addition to the T<sub>1</sub> map. For complex data (magnitude and phase, or real and imaginary), the MP2RAGE signal (<i>S</i><sub>MP2RAGE</sub>) is calculated from the images acquired at two TIs (<i>S</i><sub>GRE,TI1</sub> and <i>S</i><sub>GRE,TI2</sub>) using the following expression (Marques et al. 2010):\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "<center><img src=\"equation1.png\" style=\"width:800px;height:auto;margin-bottom: 50px;margin-top: 50px;\"></center>\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "This value is bounded between [-0.5, 0.5], and helps reduce some B<sub>0</sub> inhomogeneity effects using the phase data. For real data, or magnitude data with polarity restoration, this metric is instead calculated as:\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "<center><img src=\"equation2.png\" style=\"width:800px;height:auto;margin-bottom: 50px;margin-top: 50px;\"></center>\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Because MP2RAGE is a hybrid of pulse sequences used for inversion recovery and VFA, the resulting signal equations are more complex. Typically, a steady state is not achieved during the short train of GRE imaging blocks, so the signal at the center of k-space for each readout (which defines the contrast weighting) will depend on the number of phase-encoding steps. For simplicity, the equations presented here assume that the 3D phase-encoding dimension is fully sampled (no partial Fourier or parallel imaging acceleration). For this case (see appendix of (Marques et al. 2010) for derivation details), the signal equations are:\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "<center><img src=\"equation3.png\" style=\"width:800px;height:auto;margin-bottom: 50px;margin-top: 50px;\"></center>\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "<center><img src=\"equation4.png\" style=\"width:800px;height:auto;margin-bottom: 50px;margin-top: 50px;\"></center>\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "where B<sub>1</sub><sup>-</sup> is the receive field sensitivity, “eff” is the adiabatic inversion pulse efficiency, ER = exp(-TR/T<sub>1</sub>), EA = exp(-TA/T<sub>1</sub>), EB = exp(-TB/T<sub>1</sub>), EC = exp(-TC/T<sub>1</sub>). The variables TA, TB, and TC are the three different delay times (TA: time between inversion pulse and beginning of the GRE<sub>1</sub> block, TB: time between the end of GRE<sub>1</sub> and beginning of GRE<sub>2</sub>, TC: time between the end of GRE<sub>2</sub> and the end of the TR). If no k-space acceleration is used (e.g. no partial Fourier or parallel imaging acceleration), then these values are TA = TI<sub>1</sub> - (n/2)TR, TB = TI<sub>2</sub> - (TI<sub>1</sub> + nTR), and TC = TR<sub>MP2RAGE</sub> - (TI<sub>2</sub> + (n/2)TR), where n is the number of voxels acquired in the 3D phase encode direction varied within each GRE block. The value m<sub>z,ss</sub> is the steady-state longitudinal magnetization prior to the inversion pulse, and is given by:\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "<center><img src=\"equation5.png\" style=\"width:800px;height:auto;margin-bottom: 50px;margin-top: 50px;\"></center>\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "<center><img src=\"equation6.png\" style=\"width:800px;height:auto;margin-bottom: 50px;margin-top: 50px;\"></center>\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "    From Eqs. (3–6), it is evident that the MP2RAGE parameter <i>S</i><sub>MP2RAGE</sub> (Eqs. (1.11) and (1.12)) cancels out the effects of receive field sensitivity, T<sub>2</sub><sup>*</sup>, and M<sub>0</sub>. The signal sensitivity related to the transmit field (B<sub>1</sub><sup>+</sup>), hidden in Eqs. (3–6) within the flip angle values <i>θ<sub>1</sub></i> and <i>θ<sub>2</sub></i>, can also be reduced by careful pulse sequence protocol design (Marques et al. 2010), but not entirely eliminated (Marques & Gruetter 2013).\n",
    "</p>\n",
    "\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<center> <h2 style=\"font-family:timesnewroman;font-size:30px\">Data Fitting</h2> </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<div class=blog_body>\n",
    "<p style=\"text-align:justify;\">\n",
    "Dictionary-based techniques such as MP2RAGE do not typically use conventional minimization algorithms (e.g. Levenberg-Marquardt) to fit signal equations to observed data. Instead, the MP2RAGE technique uses pre-calculated signal values for a wide range of parameter values (e.g. T<sub>1</sub>), and then interpolation is done within this dictionary of values to estimate the T<sub>1</sub> value that matches the observed signal. This approach results in rapid post-processing times because the dictionaries can be simulated/generated prior to scanning and interpolating between these values is much faster than most fitting algorithms. This means that the quantitative image can be produced and displayed directly on the MRI scanner console rather than needing to be fitted offline.\n",
    "</p>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<div class=figure_caption>\n",
    "<p style=\"text-align:justify;\">\n",
    "<b>\n",
    "    Figure 2. T<sub>1</sub> lookup table as a function of B<sub>1</sub> and <i>S</i><sub>MP2RAGE</sub> value. Inversion times used to acquire this magnitude image dataset were 800 ms and 2700 ms, the flip angles were 4° and 5° (respectively),  TR<sub>MP2RAGE</sub> = 6000 ms, and TR = 6.7 ms. The code that was used were shared open sourced by the authors of the original MP2RAGE paper (<a href=\"url\">https://github.com/JosePMarques/MP2RAGE-related-scripts</a>).\n",
    "</b>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Octave"
   },
   "outputs": [],
   "source": [
    "%% MATLAB/OCTAVE CODE\n",
    "% Code used to generate the data required for Figure 5 of the blog post\n",
    "\n",
    "clear all\n",
    "cd ../qMRLab\n",
    "startup\n",
    "\n",
    "\n",
    "MP2RAGE.B0=7;           % in Tesla\n",
    "MP2RAGE.TR=6;           % MP2RAGE TR in seconds \n",
    "MP2RAGE.TRFLASH=6.7e-3; % TR of the GRE readout\n",
    "MP2RAGE.TIs=[800e-3 2700e-3];% inversion times - time between middle of refocusing pulse and excitatoin of the k-space center encoding\n",
    "MP2RAGE.NZslices=[35 72];% Slices Per Slab * [PartialFourierInSlice-0.5  0.5]\n",
    "MP2RAGE.FlipDegrees=[4 5];% Flip angle of the two readouts in degrees\n",
    "invEFF=0.99;\n",
    "\n",
    "B1_vector=0.005:0.05:1.9;\n",
    "T1_vector=0.5:0.05:5.2;\n",
    "\n",
    "npoints=40;\n",
    "\n",
    "\n",
    "%% creates a lookup table of MP2RAGE intensities as a function of B1 and T1\n",
    "\n",
    "k=0;\n",
    "for b1val=B1_vector\n",
    "    k=k+1;\n",
    "    [Intensity T1vector ]=MP2RAGE_lookuptable(2,MP2RAGE.TR,MP2RAGE.TIs,b1val*MP2RAGE.FlipDegrees,MP2RAGE.NZslices,MP2RAGE.TRFLASH,'normal',invEFF);\n",
    "    MP2RAGEmatrix(k,:)=interp1(T1vector,Intensity,T1_vector);\n",
    "end;\n",
    "\n",
    "%% make the matrix  MP2RAGEMatrix into T1_matrix(B1,ratio)\n",
    "MP2RAGE_vector=linspace(-0.5,0.5,npoints);\n",
    "k=0;\n",
    "for b1val=B1_vector\n",
    "    k=k+1;\n",
    "    try\n",
    "\n",
    "        T1matrix(k,:)=interp1(MP2RAGEmatrix(k,:),T1_vector,MP2RAGE_vector,'pchirp');\n",
    "    catch\n",
    "        temp=MP2RAGEmatrix(k,:);\n",
    "        temp(isnan(temp))=linspace(-0.5,-1,length(find(isnan(temp)==1)));\n",
    "        temp=interp1(temp,T1_vector,MP2RAGE_vector);\n",
    "\n",
    "        T1matrix(k,:)=temp;\n",
    "        \n",
    "    end\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Python3",
    "tags": [
     "scratch"
    ]
   },
   "outputs": [],
   "source": [
    "%get T1matrix --from Octave\n",
    "%get T1matrix --from Octave\n",
    "%get T1matrix --from Octave\n",
    "%get T1matrix --from Octave\n",
    "%get B1_vector --from Octave\n",
    "%get MP2RAGE_vector --from Octave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Python3",
    "tags": [
     "report_output",
     "scratch"
    ]
   },
   "outputs": [],
   "source": [
    "# PYTHON CODE\n",
    "\n",
    "init_notebook_mode(connected=True)\n",
    "# The polling here is to ensure that plotly.js has already been loaded before\n",
    "# setting display alignment in order to avoid a race condition.\n",
    "\n",
    "trace5 = go.Heatmap(x = MP2RAGE_vector,\n",
    "                   y = B1_vector,\n",
    "                   z=T1matrix,\n",
    "                   zmin=0,\n",
    "                   zmax=5,\n",
    "                   colorscale='Portland',\n",
    "                   xaxis='x2',\n",
    "                   yaxis='y2',\n",
    "                   visible=True,\n",
    "                   name = 'T1 values (ms)')\n",
    "\n",
    "data=[trace5]\n",
    "\n",
    "layout = dict(\n",
    "    width=560,\n",
    "    height=345,\n",
    "    margin = dict(\n",
    "                t=40,\n",
    "                r=50,\n",
    "                b=50,\n",
    "                l=80),\n",
    "    annotations=[\n",
    "         dict(\n",
    "            x=-0.14,\n",
    "            y=0.5,\n",
    "            showarrow=False,\n",
    "            text='B<sub>1</sub> values',\n",
    "            font=dict(\n",
    "                family='Times New Roman',\n",
    "                size=22\n",
    "            ),\n",
    "            textangle=-90,\n",
    "            xref='paper',\n",
    "            yref='paper'\n",
    "        ),\n",
    "          dict(\n",
    "            x=0.5,\n",
    "            y=-0.2,\n",
    "            showarrow=False,\n",
    "            text='S<sub>MP2RAGE</sub> values',\n",
    "            font=dict(\n",
    "                family='Times New Roman',\n",
    "                size=22\n",
    "            ),\n",
    "            xref='paper',\n",
    "            yref='paper'\n",
    "        ),\n",
    "        dict(\n",
    "            x=0.5,\n",
    "            y=1.15,\n",
    "            showarrow=False,\n",
    "            text='Lookup Table',\n",
    "            font=dict(\n",
    "                family='Times New Roman',\n",
    "                size=26\n",
    "            ),\n",
    "            xref='paper',\n",
    "            yref='paper'\n",
    "        ),\n",
    "        dict(\n",
    "            x=1.17,\n",
    "            y=1.15,\n",
    "            showarrow=False,\n",
    "            text='T<sub>1</sub> (ms)',\n",
    "            font=dict(\n",
    "                family='Times New Roman',\n",
    "                size=20\n",
    "            ),\n",
    "            xref='paper',\n",
    "            yref='paper'\n",
    "        ),\n",
    "    ],\n",
    "    xaxis = dict(range = [0,1], autorange = False,\n",
    "             showgrid = False, zeroline = False, showticklabels = False,\n",
    "             ticks = '', domain=[0, 0.5]),\n",
    "    yaxis = dict(range = [0,1], autorange = False,\n",
    "             showgrid = False, zeroline = False, showticklabels = False,\n",
    "             ticks = '', domain=[0, 1]),\n",
    "    showlegend = False,\n",
    "    autosize = False,\n",
    ")\n",
    "\n",
    "\n",
    "fig = dict(data=data, layout=layout)\n",
    "\n",
    "iplot(fig, filename = 'basic-heatmap', config = config)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<div class=blog_body>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "To produce T<sub>1</sub> maps with good accuracy and precision using dictionary-based interpolation methods, it is important that the signal curves are unique for each parameter value. MP2RAGE can produce good T<sub>1</sub> maps by using a dictionary with only dimensions (T<sub>1</sub>, <i>S</i><sub>MP2RAGE</sub>), since <i>S</i><sub>MP2RAGE</sub> is unique for each T<sub>1</sub> value for a given protocol  (Marques et al. 2010). However, as was noted above, <i>S</i><sub>MP2RAGE</sub> is also sensitive to B<sub>1</sub> because of <i>θ<sub>1</sub></i> and <i>θ<sub>2</sub></i> in Eqs. (1.13–1.16). The  B<sub>1</sub>-sensitivity can be reduced substantially with careful MP2RAGE protocol optimization (Marques et al. 2010), and further improved by including B<sub>1</sub> as one of the dictionary dimensions [T<sub>1</sub>, B<sub>1</sub>, <i>S</i><sub>MP2RAGE</sub>] (Figure 1.15).  This requires an additional acquisition of a B<sub>1</sub> map (Marques & Gruetter 2013), which lengthens the scan time. \n",
    "</p>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<div class=figure_caption>\n",
    "<p style=\"text-align:justify;\">\n",
    "<b>\n",
    "Figure 3. Example MP2RAGE dataset of a healthy adult brain at 7T and T<sub>1</sub> map. Inversion times used to acquire this magnitude image dataset were 800 ms and 2700 ms, the flip angles were 4° and 5° (respectively),  TR<sub>MP2RAGE</sub> = 6000 ms, and TR = 6.7 ms. The dataset and code that was used were shared open sourced by the authors of the original MP2RAGE paper (<a href=\"url\">https://github.com/JosePMarques/MP2RAGE-related-scripts</a>).\n",
    "\n",
    "</b>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Octave"
   },
   "outputs": [],
   "source": [
    "%% MATLAB/OCTAVE CODE\n",
    "% Download variable flip angle brain MRI data for Figure 7 of the blog post\n",
    "\n",
    "cmd = ['curl -L -o mp2rage.zip https://osf.io/8x2c9/download?version=1'];\n",
    "[STATUS,MESSAGE] = unix(cmd);\n",
    "unzip('mp2rage.zip');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Octave"
   },
   "outputs": [],
   "source": [
    "%% MATLAB/OCTAVE CODE\n",
    "% Code used to generate the data required for Figure 5 of the blog post\n",
    "\n",
    "clear all\n",
    "cd ../qMRLab\n",
    "startup\n",
    "\n",
    "MP2RAGE.B0=7;           % in Tesla\n",
    "MP2RAGE.TR=6;           % MP2RAGE TR in seconds \n",
    "MP2RAGE.TRFLASH=6.7e-3; % TR of the GRE readout\n",
    "MP2RAGE.TIs=[800e-3 2700e-3];% inversion times - time between middle of refocusing pulse and excitatoin of the k-space center encoding\n",
    "MP2RAGE.NZslices=[35 72];% Slices Per Slab * [PartialFourierInSlice-0.5  0.5]\n",
    "MP2RAGE.FlipDegrees=[4 5];% Flip angle of the two readouts in degrees\n",
    "MP2RAGE.filenameUNI='MP2RAGE_UNI.nii'; % file with UNI\n",
    "MP2RAGE.filenameINV1='MP2RAGE_INV1.nii'; % file with UNI \n",
    "MP2RAGE.filenameINV2='MP2RAGE_INV2.nii'; % file with INV2 \n",
    "\n",
    "% load the MP2RAGE data - it can be either the SIEMENS one scaled from\n",
    "% 0 4095 or the standard -0.5 to 0.5\n",
    "MP2RAGEimg=load_untouch_nii(MP2RAGE.filenameUNI, [], [], [], [], [], []);\n",
    "MP2RAGEINV1img=load_untouch_nii(MP2RAGE.filenameINV1, [], [], [], [], [], []);\n",
    "MP2RAGEINV2img=load_untouch_nii(MP2RAGE.filenameINV2, [], [], [], [], [], []);\n",
    "\n",
    "[T1map , M0map , R1map]=T1M0estimateMP2RAGE(MP2RAGEimg,MP2RAGEINV2img,MP2RAGE,0.96);\n",
    "\n",
    "% Code used to re-orient the images to make pretty figures, and to assign variables with the axis lengths.\n",
    "\n",
    "\n",
    "T1_map = imrotate(squeeze(T1map.img(:,130,:)),180);\n",
    "T1_map(T1map.img>5)=0;\n",
    "T1_map = T1_map*1000; % Convert to ms\n",
    "\n",
    "xAxis = [0:size(T1_map,2)-1];\n",
    "yAxis = [0:size(T1_map,1)-1];\n",
    "\n",
    "% Raw MRI data at different TI values\n",
    "S_INV1 = imrotate(squeeze(MP2RAGEINV1img.img(:,130,:)),180);\n",
    "S_INV2 = imrotate(squeeze(MP2RAGEINV2img.img(:,130,:)),180);\n",
    "B1map = imrotate(-0.5+1/4095*double(squeeze(MP2RAGEimg.img(:,130,:))),180);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Python3",
    "tags": [
     "scratch"
    ]
   },
   "outputs": [],
   "source": [
    "%get T1_map --from Octave\n",
    "%get S_INV1 --from Octave\n",
    "%get S_INV2 --from Octave\n",
    "%get B1map --from Octave\n",
    "%get xAxis --from Octave\n",
    "%get yAxis --from Octave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Python3",
    "tags": [
     "report_output",
     "scratch"
    ]
   },
   "outputs": [],
   "source": [
    "init_notebook_mode(connected=True)\n",
    "# The polling here is to ensure that plotly.js has already been loaded before\n",
    "# setting display alignment in order to avoid a race condition.\n",
    "\n",
    "trace1 = go.Heatmap(x = xAxis,\n",
    "                   y = yAxis,\n",
    "                   z=S_INV1,\n",
    "                   colorscale='Greys',\n",
    "                   showscale = False,\n",
    "                   visible=False,\n",
    "                   name = 'Signal')\n",
    "trace2 = go.Heatmap(x = xAxis,\n",
    "                   y = yAxis,\n",
    "                   z=S_INV2,\n",
    "                   colorscale='Greys',\n",
    "                   showscale = False,\n",
    "                   visible=False,\n",
    "                   name = 'Signal')\n",
    "trace3 = go.Heatmap(x = xAxis,\n",
    "                   y = yAxis,\n",
    "                   z=B1map,\n",
    "                   zmin=-0.5,\n",
    "                   zmax=0.5,\n",
    "                   colorscale='RdBu',\n",
    "                   showscale = False,\n",
    "                   visible=True,\n",
    "                   name = 'S_MP2RAGE')\n",
    "trace5 = go.Heatmap(x = xAxis,\n",
    "                   y = yAxis,\n",
    "                   z=T1_map,\n",
    "                   zmin=0.0,\n",
    "                   zmax=5000,\n",
    "                   colorscale='Portland',\n",
    "                   xaxis='x2',\n",
    "                   yaxis='y2',\n",
    "                   visible=True,\n",
    "                   name = 'T1 values (ms)')\n",
    "\n",
    "data=[trace1, trace2, trace3, trace5]\n",
    "\n",
    "\n",
    "updatemenus = list([\n",
    "    dict(active=2,\n",
    "         x = 0.09,\n",
    "         xanchor = 'left',\n",
    "         y = -0.15,\n",
    "         yanchor = 'bottom',\n",
    "         direction = 'up',\n",
    "         font=dict(\n",
    "                family='Times New Roman',\n",
    "                size=16\n",
    "            ),\n",
    "         buttons=list([   \n",
    "            dict(label = 'S<sub>INV1</sub>',\n",
    "                 method = 'update',\n",
    "                 args = [{'visible': [True, False, False, True]},\n",
    "                         ]),\n",
    "            dict(label = 'S<sub>INV2</sub>',\n",
    "                 method = 'update',\n",
    "                 args = [{'visible': [False, True, False, True]},\n",
    "                           ]),\n",
    "            dict(label = 'S<sub>MP2RAGE</sub>',\n",
    "                 method = 'update',\n",
    "                 args = [{'visible': [False, False, True, True]},\n",
    "                           ])\n",
    "        ])\n",
    "    )\n",
    "])\n",
    "\n",
    "layout = dict(\n",
    "    width=560,\n",
    "    height=345,\n",
    "    margin = dict(\n",
    "                t=40,\n",
    "                r=50,\n",
    "                b=10,\n",
    "                l=50),\n",
    "    annotations=[\n",
    "        dict(\n",
    "            x=0.055,\n",
    "            y=1.15,\n",
    "            showarrow=False,\n",
    "            text='Input Data',\n",
    "            font=dict(\n",
    "                family='Times New Roman',\n",
    "                size=26\n",
    "            ),\n",
    "            xref='paper',\n",
    "            yref='paper'\n",
    "        ),\n",
    "        dict(\n",
    "            x=0.6,\n",
    "            y=1.15,\n",
    "            showarrow=False,\n",
    "            text='T<sub>1</sub> map',\n",
    "            font=dict(\n",
    "                family='Times New Roman',\n",
    "                size=26\n",
    "            ),\n",
    "            xref='paper',\n",
    "            yref='paper'\n",
    "        ),\n",
    "        dict(\n",
    "            x=1.22,\n",
    "            y=1.15,\n",
    "            showarrow=False,\n",
    "            text='T<sub>1</sub> (ms)',\n",
    "            font=dict(\n",
    "                family='Times New Roman',\n",
    "                size=26\n",
    "            ),\n",
    "            xref='paper',\n",
    "            yref='paper'\n",
    "        ),\n",
    "    ],\n",
    "    xaxis = dict(range = [0,206], autorange = False,\n",
    "             showgrid = False, zeroline = False, showticklabels = False,\n",
    "             ticks = '', domain=[0, 0.58]),\n",
    "    yaxis = dict(range = [0,215], autorange = False,\n",
    "             showgrid = False, zeroline = False, showticklabels = False,\n",
    "             ticks = '', domain=[0, 1]),\n",
    "    xaxis2 = dict(range = [0,206], autorange = False,\n",
    "             showgrid = False, zeroline = False, showticklabels = False,\n",
    "             ticks = '', domain=[0.40, 0.98]),\n",
    "    yaxis2 = dict(range = [0,215], autorange = False,\n",
    "             showgrid = False, zeroline = False, showticklabels = False,\n",
    "             ticks = '', domain=[0, 1], anchor='x2'),\n",
    "    showlegend = False,\n",
    "    autosize = False,\n",
    "    updatemenus=updatemenus\n",
    ")\n",
    "\n",
    "\n",
    "fig = dict(data=data, layout=layout)\n",
    "\n",
    "iplot(fig, filename = 'basic-heatmap', config = config)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<div class=blog_body>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "The MP2RAGE pulse sequence is increasingly being distributed by MRI vendors, thus typically a data fitting package is also available to reconstruct the T<sub>1</sub> maps online. Alternatively, several open source packages to create T<sub>1</sub> maps from MP2RAGE data are available online (Marques 2017; de Hollander 2017), and for new users these are recommended—as opposed to programming one from scratch—as there are many potential pitfalls (e.g. adjusting the equations to handle partial Fourier or parallel imaging acceleration).\n",
    "</p>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<center> <h2 style=\"font-family:timesnewroman;font-size:30px\">Benefits and Pitfalls</h2> </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<div class=blog_body>\n",
    "<p style=\"text-align:justify;\">\n",
    "This widespread availability and its turnkey acquisition/fitting procedures are a main contributing factor to the growing interest for including quantitative T<sub>1</sub> maps in clinical and neuroscience studies. T<sub>1</sub> values measured using MP2RAGE showed  high levels of reproducibility for the brain in an inter- and intra-site study at eight sites (same MRI hardware/software and at 7 T) of two subjects (Voelker et al. 2016). Not only does MP2RAGE have one of the fastest acquisition and post-processing times among quantitative T<sub>1</sub> mapping techniques, but it can accomplish this while acquiring very high resolution T<sub>1</sub> maps (1 mm isotropic at 3T and submillimeter at 7T, both in under 10 min (Fujimoto et al. 2014)), opening the doors to cortical studies which greatly benefit from the smaller voxel size (Waehnert et al. 2014; Beck et al. 2018; Haast et al. 2018).\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Despite these benefits, MP2RAGE and similar dictionary-based techniques have certain limitations that are important to consider before deciding to incorporate them in a study. Good reproducibility of the quantitative T<sub>1</sub> map is dependent on using one pre-calculated dictionary. If two different dictionaries are used (e.g. cross-site with different MRI vendors), the differences in the dictionary interpolations will likely result in minor differences in T<sub>1</sub> estimates for the same data. Also, although the B1-sensitivity of the MP2RAGE T<sub>1</sub> maps can be reduced with proper protocol optimization, it can be substantial enough that further correction using a measured B<sub>1</sub> map should be done (Marques & Gruetter 2013; Haast et al. 2018). However B<sub>1</sub> mapping brings an additional potential source of error, so carefully selecting a B<sub>1</sub> mapping technique and accompanying post-processing method (e.g. filtering) should be done before integrating it in a T<sub>1</sub> mapping protocol (Boudreau et al. 2017). Lastly, the MP2RAGE equations (and thus, dictionaries) assume monoexponential longitudinal relaxation, and this has been shown to result in suboptimal estimates of the long T<sub>1</sub> component for a biexponential relaxation model (Rioux et al. 2016), an effect that becomes more important at higher fields.\n",
    "</p>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<center> <h2 style=\"font-family:timesnewroman;font-size:30px\">Works Cited</h2> </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<div class=biblio_body>\n",
    "<p style=\"text-align:justify;\">\n",
    "Beck, E.S. et al., 2018. Improved Visualization of Cortical Lesions in Multiple Sclerosis Using 7T MP2RAGE. <i>Am. J. Neuroradiol.</i> Available at: <a href=\"http://dx.doi.org/10.3174/ajnr.A5534\">http://dx.doi.org/10.3174/ajnr.A5534</a>.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Boudreau, M. et al., 2017. B<sub>1</sub> mapping for bias-correction in quantitative T<sub>1</sub> imaging of the brain at 3T using standard pulse sequences. <i>J. Magn. Reson. Imaging</i>, 46(6), pp.1673–1682.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Doneva, M. et al., 2010. Compressed sensing reconstruction for magnetic resonance parameter mapping. <i>Magn. Reson. Med.</i>, 64(4), pp.1114–1120.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Fujimoto, K. et al., 2014. Quantitative comparison of cortical surface reconstructions from MP2RAGE and multi-echo MPRAGE data at 3 and 7 T. <i>NeuroImage</i>, 90, pp.60–73.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Haase, A. et al., 1989. Inversion recovery snapshot FLASH MR imaging. <i>J. Comput. Assist. Tomogr.</i>, 13(6), pp.1036–1040.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Haast, R.A.M., Ivanov, D. & Uludağ, K., 2018. The impact of B<sub>1</sub><sup>+</sup> correction on MP2RAGE cortical T<sub>1</sub> and apparent cortical thickness at 7T. <i>Hum. Brain Mapp.</i>, 39(6), pp.2412–2425.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "de Hollander, G., 2017. PyMP2RAGE. Available at: <a href=\"https://github.com/Gilles86/pymp2rage\">https://github.com/Gilles86/pymp2rage</a> [Accessed January 2, 2019].\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Li, W., Griswold, M. & Yu, X., 2012. Fast cardiac T<sub>1</sub> mapping in mice using a model-based compressed sensing method. <i>Magn. Reson. Med.</i>, 68(4), pp.1127–1134.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Ma, D. et al., 2013. Magnetic resonance fingerprinting. <i>Nature</i>, 495(7440), pp.187–192.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Marques, J., 2017. MP2RAGE related scripts. Available at: <a href=\"https://github.com/JosePMarques/MP2RAGE-related-scripts\">https://github.com/JosePMarques/MP2RAGE-related-scripts</a> [Accessed January 2, 2019].\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Marques, J.P. et al., 2010. MP2RAGE, a self bias-field corrected sequence for improved segmentation and T<sub>1</sub>-mapping at high field. <i>NeuroImage</i>, 49(2), pp.1271–1281.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Marques, J.P. & Gruetter, R., 2013. New Developments and Applications of the MP2RAGE Sequence - Focusing the Contrast and High Spatial Resolution R<sub>1</sub> Mapping. <i>PloS one</i>, 8(7), p.e69294.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Mugler, J.P., 3rd & Brookeman, J.R., 1990. Three-dimensional magnetization-prepared rapid gradient-echo imaging (3D MP RAGE). <i>Magn. Reson. Med.</i>, 15(1), pp.152–157.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Rioux, J.A., Levesque, I.R. & Rutt, B.K., 2016. Biexponential longitudinal relaxation in white matter: Characterization and impact on T<sub>1</sub> mapping with IR-FSE and MP2RAGE. <i>Magn. Reson. Med.</i>, 75(6), pp.2265–2277.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Voelker, M.N. et al., 2016. The traveling heads: multicenter brain imaging at 7 Tesla. <i>Magma</i>, 29(3), pp.399–415.\n",
    "</p>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "Waehnert, M.D. et al., 2014. Anatomically motivated modeling of cortical laminae. <i>NeuroImage</i>, 93 Pt 2, pp.210–220.\n",
    "</p>\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<center> <h2 style=\"font-family:timesnewroman;font-size:30px\">Errata</h2> </center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "kernel": "Python3"
   },
   "source": [
    "<div class=biblio_body>\n",
    "<ol> \n",
    "<li><p style=\"text-align:justify;\">\n",
    "A previous version of this text incorrectly stated the equations for TB and TC. They were written as TB = TI<sub>2</sub> - TI<sub>1</sub> + (n/2)TR and TC = TR<sub>MP2RAGE</sub> - TI<sub>2</sub> + (n/2)TR, whereas they should have been TB = TI<sub>2</sub> - (TI<sub>1</sub> + nTR) and TC = TR<sub>MP2RAGE</sub> - (TI<sub>2</sub> + (n/2)TR). This has been corrected in the text. Thank you to Marc-Antoine Fortin (<a href=\"https://www.linkedin.com/in/marc-antoine-fortin-7b3777163/?originalSubdomain=ca\" target=\"_blank\">LinkedIn</a>, <a href=\"https://twitter.com/marcus_pulus\" target=\"_blank\">Twitter</a>) from the <a href=\"https://www.mcgill.ca/medphys/mri-methods-research-group\" target=\"_blank\">MRI Methods Research Group</a> at McGill University for notifying us of this error. \n",
    "<p style=\"text-align:justify;\">\n",
    "</li>\n",
    "</ol>\n",
    "\n",
    "<p style=\"text-align:justify;\">\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "kernel": "Python3"
   },
   "outputs": [],
   "source": [
    "# PYTHON CODE\n",
    "\n",
    "display(HTML(\n",
    "    '<style type=\"text/css\">'\n",
    "    '.output_subarea {'\n",
    "        'display: block;'\n",
    "        'margin-left: auto;'\n",
    "        'margin-right: auto;'\n",
    "    '}'\n",
    "    '.blog_body {'\n",
    "        'line-height: 2;'\n",
    "        'font-family: timesnewroman;'\n",
    "        'font-size: 18px;'\n",
    "        'margin-left: 0px;'\n",
    "        'margin-right: 0px;'\n",
    "    '}'\n",
    "    '.biblio_body {'\n",
    "        'line-height: 1.5;'\n",
    "        'font-family: timesnewroman;'\n",
    "        'font-size: 18px;'\n",
    "        'margin-left: 0px;'\n",
    "        'margin-right: 0px;'\n",
    "    '}'\n",
    "    '.note_body {'\n",
    "        'line-height: 1.25;'\n",
    "        'font-family: timesnewroman;'\n",
    "        'font-size: 18px;'\n",
    "        'margin-left: 0px;'\n",
    "        'margin-right: 0px;'\n",
    "        'color: #696969'\n",
    "    '}'\n",
    "    '.figure_caption {'\n",
    "        'line-height: 1.5;'\n",
    "        'font-family: timesnewroman;'\n",
    "        'font-size: 16px;'\n",
    "        'margin-left: 0px;'\n",
    "        'margin-right: 0px'\n",
    "    '</style>'\n",
    "))"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "SoS",
   "language": "sos",
   "name": "sos"
  },
  "language_info": {
   "codemirror_mode": "sos",
   "file_extension": ".sos",
   "mimetype": "text/x-sos",
   "name": "sos",
   "nbconvert_exporter": "sos_notebook.converter.SoS_Exporter",
   "pygments_lexer": "sos"
  },
  "sos": {
   "kernels": [
    [
     "Octave",
     "octave",
     "Octave",
     "#dff8fb"
    ],
    [
     "Python3",
     "python3",
     "Python3",
     "#FFD91A"
    ]
   ],
   "panel": {
    "displayed": true,
    "height": 0,
    "style": "side"
   },
   "version": "0.17.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}