mrImaging.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MR Imaging"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Class:**     Psych 204a\n",
    "\n",
    "**Tutorial:**  MR Imaging\n",
    "\n",
    "**Author:**    Wandell\n",
    "\n",
    "**Date:**      03.15.04\n",
    "\n",
    "**Duration:**  90 minutes\n",
    "\n",
    "**Copyright:** Stanford University, Brian A. Wandell\n",
    "\n",
    "**Checked:**  \n",
    "\n",
    "- Oct 2007:   Rory Sayres\n",
    "- Sep 2009:   Jon Winawer\n",
    "\n",
    "Translated to Python by Michael Waskom, 10/2012"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial covers two related topics.  First, pulse sequence methods\n",
    "for measuring T1 or T2 are illustrated.  These are introduced by\n",
    "explaining two types of pulse sequences (inversion recovery and spin\n",
    "echo) that selectively emphasize the T1 or T2 tissue properties.\n",
    "\n",
    "Then, the tutorial continues with a very simple example of how one can\n",
    "form images of the relaxation constants (T1 or T2) at different positions\n",
    " \n",
    "**You should complete the tutorial \"mrTutMR\" prior to this one.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we set up the pylab environment and import a few display-relevant functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams[\"figure.figsize\"] = (8, 6)\n",
    "mpl.rcParams[\"axes.grid\"] = True\n",
    "from IPython.display import display, clear_output\n",
    "from time import sleep"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define two small functions to convert between coordinate systems.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def cart2pol(x, y):\n",
    "    theta = arctan2(y, x)\n",
    "    r = sqrt(x ** 2 + y ** 2)\n",
    "    return theta, r\n",
    "\n",
    "def pol2cart(theta, r):\n",
    "    x = r * cos(theta)\n",
    "    y = r * sin(theta)\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will define several functions that will be employed in the tutorial to illustrate concepts. The narrative tutorial begins below the function definitions. If you are curious about how to draw the plots, or to dig deeper into the code demonstrating the principles, feel free to read these. However, they are mostly for manipulating plots, so don't feel the need to concern yourself with the details.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def inversion_recovery(T1, f=None, ax=None):\n",
    "    \"\"\"Graphical illustration of the Inversion Recovery idea.\"\"\"\n",
    "    if not all([f, ax]):\n",
    "        f, ax = subplots(1, 1)\n",
    "        ax.set_aspect(\"equal\")\n",
    "        ax.set_xlim(-10, 10)\n",
    "        ax.set_ylim(-10, 10)   \n",
    "    \n",
    "    # Original net magnetization\n",
    "    m = [0, 10]\n",
    "    p0, r0 = cart2pol(*m)\n",
    "    arr = Arrow(0, 0, *m)\n",
    "    ax.add_patch(arr)\n",
    "    ax.set_title(\"Initial magnetization\")\n",
    "    display(f)\n",
    "    sleep(1.5)\n",
    "\n",
    "    # Flip the magnetization with a 180 pulse\n",
    "    ax.set_title(\"$180^o$ RF pulse\")\n",
    "    p, r = cart2pol(*m)\n",
    "    p = p0 - pi\n",
    "    m = pol2cart(p, r)\n",
    "    arr.remove()\n",
    "    arr = Arrow(0, 0, *m)\n",
    "    ax.add_patch(arr)\n",
    "    clear_output(True)\n",
    "    display(f)\n",
    "    sleep(.5)\n",
    "\n",
    "    # Let the T1 value decay for a while at the T1 rate.\n",
    "    ax.set_title(\"Recovery\")\n",
    "    for t in arange(0, 0.8, 0.15):  # Time in seconds\n",
    "        arr.remove()\n",
    "        r = r0 - 2 * r0 *exp(-t / T1)\n",
    "        p = p0 if r > 0 else p\n",
    "        m = pol2cart(p, abs(r))\n",
    "        arr = Arrow(0, 0, *m)\n",
    "        ax.add_patch(arr)\n",
    "        clear_output(True)\n",
    "        display(f)\n",
    "        sleep(.2)\n",
    "    \n",
    "    # Rotate the magnetization using a 90 deg pulse to measure it\n",
    "    ax.set_title(\"$90^o$ RF pulse\")\n",
    "    for i in range(5):\n",
    "        p, r = cart2pol(*m)\n",
    "        p = p - (pi / 2) / 5\n",
    "        m = pol2cart(p, r)\n",
    "        arr.remove()\n",
    "        arr = Arrow(0, 0, *m)\n",
    "        ax.add_patch(arr)\n",
    "        clear_output(True)\n",
    "        display(f)\n",
    "        sleep(.2)\n",
    "\n",
    "    ax.set_title(\"Measure RF now\")\n",
    "    clear_output(True)\n",
    "    display(f)\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def spin_echo(TE=16, f=None, ax=None):\n",
    "    \"\"\"Graphical illustration of the Spin Echo dephasing and echo formation.\"\"\"\n",
    "    if not all(f, ax):\n",
    "        f, ax = subplots(1, 1)\n",
    "        ax.set_aspect(\"equal\")\n",
    "        ax.set_xlim(-10, 10)\n",
    "        ax.set_ylim(-10, 10)\n",
    "        ax.set_xlabel(\"X axis\")\n",
    "        ax.set_ylabel(\"Y axis\")\n",
    "    \n",
    "    # Original net magnetization\n",
    "    ma = [10, 0]\n",
    "    p0, r0 = cart2pol(*ma)\n",
    "    arr_a = Arrow(0, 0, *ma, color=\"blue\")\n",
    "    ax.add_patch(arr_a)\n",
    "    mb = ma\n",
    "    arr_b = arrow(0, 0, *mb, color=\"green\")\n",
    "    ax.set_title(\"$90^o$ pulse and dephasing\")\n",
    "    display(f)\n",
    "    sleep(1.2)\n",
    "    \n",
    "    for i in range(TE):\n",
    "        arr_a.remove()\n",
    "        arr_b.remove()\n",
    "        pa, ra = cart2pol(*ma)\n",
    "        pa = pa - pi / 8\n",
    "        ma = pol2cart(pa, ra)\n",
    "        arr_a = Arrow(0, 0, *ma, color=\"blue\")\n",
    "        ax.add_patch(arr_a)\n",
    "        pb, rb = cart2pol(*mb)\n",
    "        pb = pb - pi / 10\n",
    "        mb = pol2cart(pb, rb)\n",
    "        arr_b = Arrow(0, 0, *mb, color=\"green\")\n",
    "        ax.add_patch(arr_b)\n",
    "        clear_output(True)\n",
    "        display(f)\n",
    "        sleep(.2)\n",
    "\n",
    "    # Apply a 180 deg pulse that rotates the spins around the x-axis\n",
    "    sleep(.8)\n",
    "    ax.set_title(\"Inverting ($180^o$ pulse)\")\n",
    "    clear_output(True)\n",
    "    display(f)\n",
    "    sleep(.5)\n",
    "    \n",
    "    arr_a.remove()\n",
    "    arr_b.remove()\n",
    "    pa, ra = cart2pol(*ma)\n",
    "    pa = -pa\n",
    "    ma = pol2cart(pa, ra)\n",
    "    arr_a = Arrow(0, 0, *ma, color=\"blue\")\n",
    "    ax.add_patch(arr_a)\n",
    "    pb, rb = cart2pol(*mb)\n",
    "    pb = -pb\n",
    "    mb = pol2cart(pb, rb)\n",
    "    arr_b = Arrow(0, 0, *mb, color=\"green\")\n",
    "    ax.add_patch(arr_b)\n",
    "    clear_output(True)\n",
    "    display(f)\n",
    "    \n",
    "    # Now keep going\n",
    "    ax.set_title(\"Catching up\")\n",
    "    clear_output(True)\n",
    "    display(f)\n",
    "    sleep(.5)\n",
    "\n",
    "    for i in range(TE):\n",
    "        arr_a.remove()\n",
    "        arr_b.remove()\n",
    "        pa, ra = cart2pol(*ma)\n",
    "        pa = pa - pi / 8\n",
    "        ma = pol2cart(pa, ra)\n",
    "        arr_a = Arrow(0, 0, *ma, color=\"blue\")\n",
    "        ax.add_patch(arr_a)\n",
    "        pb, rb = cart2pol(*mb)\n",
    "        pb = pb - pi / 10\n",
    "        mb = pol2cart(pb, rb)\n",
    "        arr_b = Arrow(0, 0, *mb, color=\"green\")\n",
    "        ax.add_patch(arr_b)\n",
    "        clear_output(True)\n",
    "        display(f)\n",
    "        sleep(.2)    \n",
    "    \n",
    "    ax.set_title(\"The echo arrives\")\n",
    "    clear_output(True)\n",
    "    display(f)\n",
    "    \n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def phase_encode(rate, spin_dir, n_steps=15):\n",
    "    \"\"\"Visualization of the phase-encoding.\"\"\"\n",
    "    f, axes = subplots(2, 2, figsize=(8, 8))\n",
    "    axes = axes.ravel()\n",
    "    for ax in axes:\n",
    "        ax.set_aspect(\"equal\")\n",
    "        ax.set_xlabel(\"X axis\")\n",
    "        ax.set_ylabel(\"Y axis\")\n",
    "        ax.set_xlim(-10, 10)\n",
    "        ax.set_ylim(-10, 10)\n",
    "    \n",
    "    a = empty(4, object)\n",
    "    for i in range(n_steps):\n",
    "        sleep(.2)\n",
    "        pa = zeros(4)\n",
    "        ra = zeros(4)\n",
    "        for j, ax in enumerate(axes):\n",
    "            pa[j], ra[j] = cart2pol(spin_dir[j, 0], spin_dir[j, 1]) \n",
    "            pa[j] = pa[j] - rate[j] \n",
    "            spin_dir[j, 0], spin_dir[j, 1] = pol2cart(pa[j], ra[j])\n",
    "            if i:\n",
    "                a[j].remove()\n",
    "            a[j] = Arrow(0, 0, *spin_dir[j, :])\n",
    "            ax.add_patch(a[j])\n",
    "        clear_output(True)\n",
    "        display(f)\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pulse Sequences for measuring T1 and T2 signals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### T1 Signals (used for anatomical images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Inversion-Recovery (IR)*\n",
    "\n",
    "Inversion-Recovery pulse sequences are a method for measuring  T1\n",
    "relaxation (spin-lattice). As the sequence name suggests, the pulse\n",
    "sequence first inverts the net magnetization ($180^o$ pulse). Then, the\n",
    "sequence simply pauses for a time, TI, to let the longitudinal\n",
    "magnetization recover towards steady state across the tissue. Then, a\n",
    "$90^o$ pulse is delivered that places the net magnetization in the\n",
    "transverse plane.  The transverse magnetization is measured right away,\n",
    "before significant dephasing, in order to estimate the T1 properties.\n",
    "\n",
    "To help you visualize the events, run this code a few times.  The first\n",
    "plot shows the 180-TI-90 sequence for a relatively slow T1 value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "inversion_recovery(T1=2.8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now suppose the tissue has a faster T1 relaxation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "inversion_recovery(T1=0.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 1:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we apply a 90-degree pulse at the time when exactly half the signal has\n",
    "recovered, what do you expect the transverse magnetization to be?\n",
    "\n",
    "*(Extra credit: is this the same as applying a pulse at T1? Why or why not?)*\n",
    "\n",
    "Comparing the two plots, you can see that the amplitude of the net\n",
    "transverse magnetization depends on the value of T1. The final $90^o$\n",
    "flip let us measure the size of the longitudinal magnetization.\n",
    "Because the measurement takes place immediately after this flip,\n",
    "there is not much time for spin-spin dephasing and we measure mainly\n",
    "properties of T1.  That is why such sequences are called 'T1-weighted'."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### T2 Signals (used for BOLD images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Spin Echo (Hahn)**\n",
    "\n",
    "In principle, to make a T2 measurement we need only to flip the net\n",
    "magnetization $90^o$ and then measure the signal decay as the spins\n",
    "dephase.  Because the T2 reduction is very fast compared to the T1\n",
    "reduction, the decay we see will be nearly all due to T2.\n",
    "\n",
    "The spin-spin dephasing occurs so quickly that it is almost impossible to\n",
    "obtain a T2 measurement soon enough after the 90 deg flip. The Hahn Spin\n",
    "Echo, and its partner the Gradient Echo, make it possible to separate the\n",
    "time between the pulse and the measurement.\n",
    "\n",
    "The next visualizataion shows two spins within a single voxel.  These\n",
    "spins are experiencing slightly differeent local magnetic fields.\n",
    "Consequently, one precesses around the origin slightly faster than the\n",
    "other.  After a little time the signals are well out of phase. Then, a\n",
    "180 deg inverting pulse is introduced.  This pulse rotates the spins\n",
    "around the horizontal axis (x-axis). Seen within the plane, this causes\n",
    "the two spins to reverse positions so that the leader becomes the\n",
    "follower. The leader is still in a higher local field, and so after a\n",
    "while it catches up.  At this moment the spin phasese come together to\n",
    "create an echo of the signal after the first 90 deg pulse.  This echo\n",
    "happens at a point in time that is well separated from the inverting\n",
    "pulse.\n",
    "\n",
    "The time until the inverse pulse determines when the echo will occur"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "spin_echo(TE=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MR Image Formation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now ask: How can we distinguish signals from different locations?  \n",
    "\n",
    "This is the key step in learning how to form an image of the MR time \n",
    "constant's properties. \n",
    "\n",
    "Consider a simple situation in which we have two beakers sitting next to\n",
    "one another on a table.  Both contain water.  To start, suppose the\n",
    "magnetic field is uniform, and suppose that we measure T1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to need these variables to develop the story"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Mo = 1                    # Net magnetization\n",
    "larmor_freq = [12, 18]    # Larmor frequency in MHz/10\n",
    "T1 = [1.5, 3]\n",
    "t = arange(0, 1, .005) * (3 * max(T1))  # Time samples in secs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def rf_signal(tau=1, net_mag=1, t_samples=None, larmor_freq=12, ph=0):\n",
    "    \"\"\"Estimate an rf signal based on the various parameters.\"\"\"\n",
    "    if t_samples is None:\n",
    "        t_samples = arange(0, 1, 0.005) * (4 * tau)\n",
    "\n",
    "    signal = exp_decay(tau, net_mag, t_samples) * cos(t_samples * larmor_freq + ph)\n",
    "    return signal\n",
    "\n",
    "def exp_decay(tau, Mo, t):\n",
    "    \"\"\"Create an exponential decay that will be used for either\n",
    "       longitudinal or transverse decay modeling.\"\"\"\n",
    "    return Mo * exp(-t / tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `rf_signal` produces the RF signal that we will measure given a particular time constant and Larmor frequency.\n",
    "\n",
    "The RF signal is periodic.  We can summarize the amplitude and\n",
    "frequency of this signal by plotting the Fourier Transform amplitude\n",
    "spectrum.  This plot measures the amplitude of each harmonic (temporal\n",
    "frequency) in the signal.\n",
    "\n",
    "We're going to be making quite a few plots with an RF signal over time and its spectral density.\n",
    "Let's write a general function so that we avoid typing the same thing over and over again.\n",
    "As a general rule of thumb for progamming, any time you find yourself using \"copy\" and \"paste\", it's time\n",
    "to write a function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def rf_plot(time, signal, title):\n",
    "    f, ax = subplots(1, 1, figsize=(6, 6))\n",
    "    ax.plot(time, signal)\n",
    "    ax.set_xlabel('Time (s)')\n",
    "    ax.set_ylabel('RF signal')\n",
    "    ax.set_title(title)\n",
    "\n",
    "def rf_kspace_plot(time, signal, titles):\n",
    "    f, (ax1, ax2) = subplots(1, 2, figsize=(14, 6))\n",
    "    ax1.plot(t, signal)\n",
    "    ax1.set_xlabel(\"Time (s)\")\n",
    "    ax1.set_ylabel(\"RF signal\")\n",
    "    ax1.set_title(titles[0])\n",
    "    ax2.psd(signal)\n",
    "    ax2.set_title(titles[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "signal = rf_signal(T1[0], Mo, t, larmor_freq[1])\n",
    "rf_kspace_plot(t, signal, [\"Beaker Signal\", \"Slice Selection\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, consider the signal as the relaxation constant (t_constant)\n",
    "increases.  Notice that over the same period of time, the rf_signal has a\n",
    "larger amplitude (there is less decay)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_constant = [1, 2, 4] \n",
    "f, axes = subplots(1, 3, figsize=(14, 5), sharey=True)\n",
    "for i, ax in enumerate(axes):\n",
    "    signal = rf_signal(t_constant[i], Mo, t, larmor_freq[0])\n",
    "    ax.set_ylim(-1, 1)\n",
    "    ax.plot(t, signal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 2:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you were to plot the Fourier Transform for each of the three \n",
    "subplots, how do you expect they would differ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_constant = [1, 2, 4]\n",
    "f, axes = subplots(1, 3, figsize=(14, 5), sharey=True)\n",
    "for i, ax in enumerate(axes):\n",
    "    signal = rf_signal(t_constant[i], Mo, t, larmor_freq[1])\n",
    "    ax.psd(signal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Changing magnetic field strength"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, suppose that we change the magnetic field strength.  Remember that the\n",
    "Larmor frequency is proportional to the magnetic field.  Consequently,\n",
    "the frequency of the RF signal will increase.  We can compare the signals\n",
    "at two different frequencies as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "larmor_freqs = [6, 124]\n",
    "f, axes = subplots(1, 2, figsize=figsize(13, 7), sharey=True)\n",
    "for i, freq in enumerate(larmor_freqs):\n",
    "    signal = rf_signal(T1[0], Mo, t, freq);\n",
    "    axes[i].plot(t, signal)\n",
    "    axes[i].set_title(\"Frequency: %d Hz\" % freq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is easiest to see the change in frequency if we compute the Fourier transform of the two signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f, axes = subplots(1, 2, figsize=figsize(13, 7), sharey=True)\n",
    "for i, freq in enumerate(larmor_freqs):\n",
    "    ax = axes[i]\n",
    "    signal = rf_signal(T1[0], Mo, t, freq)\n",
    "    ax.psd(signal)\n",
    "    ax.set_title('Fourier Spectrum for %d Hz Signal' % freq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This figure is important.  It shows that the frequency of the response\n",
    "from each beaker depends on the local magnetic field.  We have already\n",
    "seen that the amplitude at the response frequency depends on the time\n",
    "constant.  We can take advantage of these two observations by introducing\n",
    "a gradient into the magnetic field. \n",
    "\n",
    "By inserting a gradient, the two beakers experience slightly different\n",
    "magnetic fields and thus the signals have two different Larmor\n",
    "frequencies. The combined RF signal from the two beakers will be the\n",
    "sum of the signals from the two beakers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "signal =  rf_signal(T1[0], Mo, t, larmor_freqs[0]) + rf_signal(T1[0], Mo, t, larmor_freqs[1])\n",
    "rf_plot(t, signal, \"Two Beakers in Gradient\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the amplitude of the signals from each of the beakers by\n",
    "plotting the Fourier Transform spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rf_kspace_plot(t, signal, [\"Two beakers in gradient\", \"Fourier Space\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, if the two beakers represent substances with different time\n",
    "constants, we will be able to measure this by estimating the amplitudes\n",
    "of the two peaks in the spectrum.  Here, we create a signal in which the\n",
    "two beakers are in a gradient and the substances have slightly different\n",
    "time constants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "signal =  rf_signal(T1[0], Mo, t, larmor_freqs[0]) + rf_signal(T1[1], Mo, t, larmor_freqs[1])\n",
    "rf_plot(t, signal, \"Two beakers in gradient\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the amplitude of peaks in the signal change, reflecting\n",
    "the different time constants.  We can distinguish which amplitude is\n",
    "associated with which beaker because their responses are at different\n",
    "frequencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rf_kspace_plot(t, signal, [\"Two Beakers in Gradient\", \"Fourier Space\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 3:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When computing the RF signal in the last example, which variable(s)\n",
    "represented the magnetic field gradient?\n",
    "\n",
    "From answering Question 2, and understanding this simple case, you should\n",
    "understand how the gradient associates different RF signal frequencies\n",
    "with different spatial locations.\n",
    "\n",
    "This is why determining the frequency amplitude of the signal, using the\n",
    "Fourier Transform, is so important:  RF signal frequency is associated\n",
    "with spatial position.  The amplitude at that frequency can be\n",
    "interpreted as the decay rate. \n",
    "\n",
    "This simple example provides you with the basic principles of an\n",
    "important MR imaging term: *k-space.* Specifically, in this example the\n",
    "frequency axis of the Fourier Transform is analogous to k-space.\n",
    "\n",
    "There are important limitations to the method we have developed up this point.\n",
    "\n",
    "Mainly, this method works if we only need to make images of an array of\n",
    "beakers along a line. To make estimates of beakers spread across a table\n",
    "top (2D) or filling up a box (3D) we need to do more.  These methods will\n",
    "be explained below. \n",
    "\n",
    "The main difference between this example and general imaging is\n",
    "dimensionality. In MR imaging the stimuli don't fall along one dimension,\n",
    "so we can't simply use a one-dimensional frequency axis to assign\n",
    "position.  In general the position in k-space corresponds to a position\n",
    "in a two-dimensional plane that represents various spatial frequencies.\n",
    "\n",
    "(By the way, there are historical reasons why they call it k-space.  I\n",
    "forget these reasons.  If someone can remind me, I will give them a firm\n",
    "handshake and a warm smile.)\n",
    "\n",
    "We will start to extend the ideas on imaging in this first section to\n",
    "slightly more complex cases in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slice Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In nearly all fMRI imaging, the data are acquired using a series of\n",
    "measurements from different planar sections (slices) of the tissue. The\n",
    "measurements are made by selectively exciting the spins, one planar\n",
    "section at a time.  If only one plane is excited, then we can be\n",
    "confident that the signal we measure must arise from a voxel in that\n",
    "plane.  \n",
    " \n",
    "In the next part of this tutorial, we will reivew how we can excite\n",
    "the spins in one planar section.  Following that we will review how we\n",
    "can distinguish different positions within the excited plane.\n",
    " \n",
    "The principles used to understand slice selection are the same as the\n",
    "principles used to distinguish signals from beakers at two positions.\n",
    "The main difference is that in the previous example we used gradients to\n",
    "distinguish received signals.  In slice selection we use gradients to\n",
    "selectively deliver an excitation.\n",
    "\n",
    "The idea is this.  We introduce a magnetic field gradient across the \n",
    "sample changing the Larmor frequency across the sample.  \n",
    "When we deliver an RF pulse at a particular frequency, then,\n",
    "only the spins in one portion of the gradient field will be excited.\n",
    "\n",
    "What would such a pulse look like? Let's create some examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's perform some example calculations. We initialize the rf_pulse to zero, and then we initialize some parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n_time = 256.\n",
    "t = arange(n_time) / n_time - 0.5\n",
    "rf_pulse = zeros(int(n_time))\n",
    "\n",
    "pulse_duration = round(n_time / 2)\n",
    "pulse_start = n_time / 2 - pulse_duration / 2\n",
    "pulse_stop = pulse_start + pulse_duration - 1\n",
    "idx = slice(int(pulse_start), int(pulse_stop))\n",
    "\n",
    "pulse_t = t[idx]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose a pulse frequency.  This frequency will excite the planar\n",
    "section of the tissue that is at the Larmor frequency we wish to excite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pulse_freq = 25"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, create a sinusoid RF pulse at that frequency:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rf_pulse[idx] = sin(2 * pi * pulse_freq * pulse_t)\n",
    "rf_plot(t, rf_pulse, \"RF pulse\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the Fourier Transform spectrum of the RF pulse.\n",
    "The frequency content of the pulse determines which plane we excite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rf_kspace_plot(t, rf_pulse, [\"RF Pulse\", \"Slice Selection\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can control the position that is excited by adjusting the frequency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pulse_freq = 50\n",
    "rf_pulse[idx] = sin(2 * pi * pulse_freq * pulse_t)\n",
    "rf_kspace_plot(t, rf_pulse, [\"RF Pulse\", \"Slice Selection\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A second parameter we would like to control is the slice width.   There\n",
    "are two ways we can adjust the slice width.  One way, as described in the\n",
    "book, is to change the gradient.  A second way, illustrated here, is to\n",
    "change the timing of the RF pulse.\n",
    "\n",
    "In this example, we create an RF pulse that is the product of the\n",
    "sinusoid and a sinc function.  (Type `sinc?` to read about this\n",
    "important function). Each sinc function has its own frequency, too. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sinc_freq  = 20   # Try 10, 20 and 40.\n",
    "pulse_freq = 50  # Try 25, 50 and 75.\n",
    "\n",
    "rf_pulse[idx] = sin(2 * pi * pulse_freq * pulse_t)\n",
    "rf_pulse[idx] = rf_pulse[idx] * sinc(sinc_freq * pulse_t)\n",
    "\n",
    "rf_kspace_plot(t, rf_pulse, [\"RF Pulse\", \"Slice Selection\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 4:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the above code a few times, varying the pulse and sinc frequency \n",
    "values.  What effect does each parameter have on slice position and \n",
    "slice width?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Image Formation using Frequency Encoding and Phase Encoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [Hornak MRI book](http://www.cis.rit.edu/htbooks/mri/) has a very\n",
    "good discussion of imaging, including frequency and phase encoding.\n",
    "Please read Chapters 6 and 7 for a useful discussion of the principles\n",
    "further described here.  Also, the Huettel et al. book is very clear and\n",
    "useful in describing pulse sequences (Chapters 4 and 5).\n",
    " \n",
    "Earlier, we reviewed how to use a gradient field to associate different\n",
    "positions along a line with different RF signal frequencies.  That method\n",
    "is often called frequency encoding.  In this  section, we describe how to\n",
    "measure along the second spatial dimension.  This measurement is\n",
    "sometimes called the phase-encoding dimension.  Taken together, the\n",
    "methods described in this section are sometimes called Fourier Transform\n",
    "imaging.\n",
    " \n",
    "Consider the problem of identifying signals from 4 beakers placed at the\n",
    "corners of a square.  These four beakers are in the planar section and the\n",
    "spins were set in motion using the methods described in the previous\n",
    "section. \n",
    "\n",
    "First, we set up the basic parameters for the four beakers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "T1 = array([(1, 2), (3, 4)], float) / 2\n",
    "larmor_freq = [15., 50.]\n",
    "ph = [0, pi]\n",
    "t = arange(0, 5, .02)\n",
    "Mo = 1\n",
    "\n",
    "rate = [0., 0., 0., 0.]\n",
    "spin_dir = array([(10, 0), (10, 0), (10, 0), (10, 0)], float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the starting positions in the transverse plane:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "phase_encode(rate, spin_dir, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define a new plotting function to save some typing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_rf_2d(fs, ps, t, T1, Mo=1):\n",
    "    f, axes = subplots(2, 2, figsize=(8, 8))\n",
    "    signal = zeros((2, 2, len(t)))\n",
    "    freq_text = [\"F1\", \"F2\"]\n",
    "    fs = reshape(fs, (2, 2)).T\n",
    "    ps = reshape(ps, (2, 2)).T\n",
    "    for i in range(2):\n",
    "        for j in range(2):\n",
    "            signal[i, j, :] = rf_signal(T1[i, j], Mo, t, fs[i, j], ps[i, j])\n",
    "            axes[i, j].plot(t, signal[i, j, :])\n",
    "            axes[i, j].set_title(freq_text[j])\n",
    "    return signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we apply a gradient across the x-axis.  In the presence of this\n",
    "gradient, the signals from the two beakers on the left will differ from\n",
    "the RF signals emitted by the two beakers on the right. When we measure\n",
    "with only the x-gradient (Gx), we obtain the sum of the two beakers on\n",
    "the left in one frequency band and the sum of the two beakers on the\n",
    "right in a second frequency band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "freqs = [larmor_freq[0]] * 2 + [larmor_freq[1]] * 2\n",
    "phases = [ph[0]] * 4\n",
    "signal = plot_rf_2d(freqs, phases, t, T1, Mo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the total signal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s = (signal[0] + signal[1]).sum(axis=0)\n",
    "rf_kspace_plot(t, s, [\"Total RF\", \"Total RF Signal\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, suppose that we turn off the Gx gradient and introduce a gradient in\n",
    "the y-dimension (Gy).  This changes the Larmor frequency of the beakers\n",
    "at the top and bottom of the square.  Suppose the frequency for the\n",
    "bottom beakers is a little lower. Then the spins in the bottom beakers\n",
    "rotate more slowly and they end up pointing in a different direction from\n",
    "the spins at the top.  After applying Gy for a certain amount of time,\n",
    "the spins at the top and bottom will point in opposite directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rate = [pi / 8, pi / 8, pi / 16, pi / 16]\n",
    "spin_dir = array([(10, 0), (10, 0), (10, 0), (10, 0)], float)\n",
    "phase_encode(rate, spin_dir, 15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we switch off  Gy and turn on Gx.  As before we will measure the\n",
    "combination of the spins on the left at one frequency and the combination\n",
    "of the spins at the right at a different frequency.  Because the spins of\n",
    "the top and bottom beaker oppose one another, however, the total RF\n",
    "signal we obtain now is the difference between the top and bottom. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ps = [ph[0], ph[1]] * 2\n",
    "signal = plot_rf_2d(freqs, ps, t, T1, Mo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Total signal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "s = (signal[0] + signal[1]).sum(axis=0)\n",
    "rf_kspace_plot(t, s, [\"Total RF\", \"Total RF Signal\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These two measurements at the main frequencies provide an estimate of the\n",
    "sum and the difference of the time constants in the upper and lower\n",
    "beakers. Specifically, the two measurements are\n",
    "\n",
    "    m1 = a + b\n",
    "    m2 = a - b\n",
    "\n",
    "Measuring m1 and m2, we can solve for a and b.\n",
    "\n",
    "In this example, the frequency of the RF signal codes the x-dimension. \n",
    "This dimension is called the frequency dimension.  We keep adjusting the \n",
    "phase of the spins along the y-dimension. Hence, this dimension is called\n",
    "the phase-encode dimension.  People keep track of these two dimensions \n",
    "because the phase-encode direction is subject to more measurement error \n",
    "than the frequency encode direction. "
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}