{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ex 4: Parallel Transmission Pulse Design"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will wrap up the exercises for this session with a type of pulse design that visibly leverages some of the nicer abstractions within SigPy: parallel transmit pulse design.\n",
    "\n",
    "Specifically, in this exercise we will focus on a spatial-domain spiral parallel transmit pulse design.\n",
    "\n",
    "We start by writing our import statements one last time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "# typical sigpy and numpy imports\n",
    "import numpy as np\n",
    "import sigpy.mri as mr\n",
    "import sigpy.mri.rf as rf # importing our rf tools separately\n",
    "import sigpy.plot as pl\n",
    "\n",
    "# to assist with importing data, we will also import scipy\n",
    "import scipy.io as sio\n",
    "import scipy.ndimage.filters as filt\n",
    "\n",
    "import matplotlib.pyplot as mplib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this problem, use a small problem geometry: a 2D 32x32 design grid with 8 channels, to keep computation simple for Binder! However, if you are running these exercises locally and have more RAM available, larger designs are no problem. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = 32\n",
    "Nc = 8\n",
    "sens_shape = [Nc, dim, dim]\n",
    "img_shape = [dim, dim]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin, we will want some two-dimensional spatial magnetization pattern to design for. We provide one pattern in a matlab file for use in the data folder of this notebook, which we will load here. This is just to get you started - feel free to load other patterns or create your own! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat_struct = sio.loadmat('data/smallv.mat')\n",
    "d = mat_struct['d'].astype(np.single)\n",
    "pl.ImagePlot(d, title='Target Excitation Pattern')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an interesting pattern, but the edges are extremely sharp. To avoid introducing Gibbs ringing into our design from these sharp edges, we will perform a very slight blur of the pattern using scipy. Our scipy filter requires real-valued numbers but we will want to use complex values in the design later, so we will cast to complex after the blurring is performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = filt.gaussian_filter(d, 0.5)\n",
    "d = d.astype(np.complex)\n",
    "pl.ImagePlot(d, title='Blurred Target')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4a: Simulate our transmit coil sensitivities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we will break out the more general sigpy.mri submodule functions for the first time. Use [sigpy.mri.birdcage_maps()](https://sigpy.readthedocs.io/en/latest/generated/sigpy.mri.birdcage_maps.html#sigpy.mri.birdcage_maps) to simulate our B1+ transmit sensitivities. Note that these simulated profiles would be equally valid to use as receive sensitivities for a reconstruction problem, such as a SENSE recon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sens = mr.birdcage_maps(sens_shape)\n",
    "pl.ImagePlot(sens)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**If you haven't already, take a moment to play around with some of the hotkey features of sigpy.plot using the (complex-valued) sensitivities plot produced above. After clicking on the plot, try the following hotkeys:**\n",
    "* Up and down arrows - scroll through currently selected dimension\n",
    "* Left and right arrows - change which dimension is selected (indicated by [brackets])\n",
    "* m - magnitude plot\n",
    "* p - phase plot\n",
    "* r - real plot\n",
    "* i - imaginary plot\n",
    "\n",
    "You should be able to look at each of the coils' B1 sensitivities sequentially, and view their magnitude, phase, or real/imaginary components! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4b: Create a spiral trajectory.\n",
    "\n",
    "Before performing our design, we will also need to design time-varying gradient waveforms to produce a trajectory for our excitation.\n",
    "\n",
    "* Design a spiral-in trajectory using [*sigpy.mri.rf.spiral_arch()*](https://sigpy.readthedocs.io/en/latest/generated/sigpy.mri.rf.trajgrad.spiral_arch.html#sigpy.mri.rf.trajgrad.spiral_arch). Note that by default this trajectory designer produces a spiral-out trajectory - **you will need to flip the waveforms to spiral inward!** Numpy has functions that can assist with this, or you can use python built-ins (such as vector = vector[::-1])\n",
    "\n",
    "\n",
    "* Given: \n",
    "    * a FOV of 0.55 m\n",
    "    * a maximum slew rate of 150 mT/m/ms\n",
    "    * a maximum gradient amplitude of 30 mT/m\n",
    "    * dx = 0.025 m\n",
    "\n",
    "\n",
    "*Hint: to under- or oversample with the spiral use fov = fov/R*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fov = 0.55  # FOV in m\n",
    "gts = 6.4e-6  # hardware dwell time, s\n",
    "gslew = 150  # gradient slew rate in mT/m/ms\n",
    "gamp = 30 # maximum gradient amplitude in mT/m\n",
    "R = 1/2  # degree of undersampling of outer region of trajectory- let's oversample by a factor of 2\n",
    "dx = 0.025 # in m\n",
    "# construct a trajectory\n",
    "g, k, t, s = rf.spiral_arch(fov/R,dx,gts,gslew,gamp)\n",
    "\n",
    "#Note that this trajectory is a spiral-out trajectory. \n",
    "#Simply time-reverse it to create a spiral-in.\n",
    "k = k[::-1]\n",
    "g = g[::-1]\n",
    "\n",
    "\n",
    "mplib.figure()\n",
    "mplib.plot(k[:,0],k[:,1], color='orange')\n",
    "mplib.title('Constant density spiral-in trajectory kx and ky')\n",
    "mplib.figure()\n",
    "mplib.plot(g)\n",
    "mplib.title('Gradient waveforms')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4c: Design and plot the pulses\n",
    "\n",
    "Do this using the [sigpy.mri.rf.stspa() small tip spatial-domain pulse designer](https://sigpy.readthedocs.io/en/latest/generated/sigpy.mri.rf.ptx.stspa.html#sigpy.mri.rf.ptx.stspa). To keep the design fast, set explicit=False and perform ~ 10 iterations (explicit=True has not performed well in Binder). This should be sufficient to produce a good pattern."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Insert code for the pulse design here: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pulses = rf.stspa(d, sens, k, gts, explicit=False, max_iter=10)\n",
    "pl.LinePlot(pulses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4d: Check your result using a SigPy linop.\n",
    "\n",
    "One way of checking that our designed pulse is reasonable is by constructing the system matrix A used within the SmallTipSpatialPulseDesign app, and looking at $$ m = A b$$\n",
    "\n",
    "* Construct the system matrix A using the [sigpy.mri.linop.Sense() linop](https://sigpy.readthedocs.io/en/latest/generated/sigpy.mri.linop.Sense.html#sigpy.mri.linop.Sense). Recall that our A is the adjoint of the SENSE linear operator for a spatial domain pulse design.\n",
    "* Get and plot the magnetization m using A and the pulses you designed. This can be done in one line! \n",
    "\n",
    "Applying and transposing Linops is documented in detail on the [Linop documentation page](https://sigpy.readthedocs.io/en/latest/generated/sigpy.linop.Linop.html#sigpy.linop.Linop).\n",
    "\n",
    "The linear operator is converted to its adjoint with ```A = A.H```\n",
    "\n",
    "You can apply the operation y = A(x) with ```y = A * x```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = mr.linop.Sense(sens, coord=k, ishape=d.shape).H\n",
    "\n",
    "m = A * pulses\n",
    "\n",
    "pl.ImagePlot(m, title = 'Mxy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 4e: Use a bloch simulator to check the actual magnetization profile.\n",
    "\n",
    "A better way of verifying the accuracy of the pulses is by performing a bloch simulation. Do so using the [sigpy.mri.rf.sim.abrm_ptx() function](https://sigpy.readthedocs.io/en/latest/generated/sigpy.mri.rf.sim.abrm_ptx.html#sigpy.mri.rf.sim.abrm_ptx).\n",
    "\n",
    "Plot the Mxy and Mz profile once you're done. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make spatial variables, in units of m\n",
    "x, y = np.meshgrid(np.linspace(fov/2, -fov/2, dim), np.linspace(fov/2, -fov/2, dim))\n",
    "spatial = np.fliplr(np.concatenate((np.reshape(x, (dim*dim, 1)), np.reshape(y, (dim*dim, 1))), axis=1))\n",
    "gam = 42580 # Hz/mT\n",
    "a, b, m, mz = rf.abrm_ptx(pulses/33, spatial, g * gam * gts * 2 * np.pi, gts, fmap=None, sens=sens)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl.ImagePlot(np.reshape(2 * a * np.conj(b),img_shape), mode='m', title=('Mxy'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pl.ImagePlot(np.reshape(1-2*np.absolute(b)**2,img_shape), mode='m', title=('Mz'))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}