{
 "metadata": {
  "name": "phantom"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Phantom Simulation\n",
      "\n",
      "With a given fiber configuration, simulate the signal measured in S-space.\n",
      "\n",
      "The volume, ``V``, generated has shape ``(M, N, P, B)``, where ``B`` is the number of b-vectors.\n",
      "Each measurement ``V[x, y, z, ...]`` is the integral of the diffusion probability density\n",
      "function along the ``B`` different measurement directions.  ``V[x, y, z, 0]`` is the reference\n",
      "measurement, $S_0$, where no gradient is applied.\n",
      "\n",
      "The signal is generated by tracing the fiber contour, determining its gradient\n",
      "at each position, and then adding the response of a single fiber along that direction.\n",
      "The response of a single fibre is given by\n",
      "\n",
      "$$\\Delta S = S_0 \\exp(-b \\mathbf{b}^\\mathrm{T} R \\Lambda R^\\mathrm{T} \\mathbf{b} ) $$\n",
      "\n",
      "where $S_0$ is the reference signal of the fibre, $\\Lambda$ the diffusion tensor, and $R$\n",
      "the rotation matrix to align the diffusion tensor with the fiber direction. For\n",
      "white matter, we typically take\n",
      "\n",
      "$$\\Lambda = \\left(\\begin{array}{ccc}\n",
      "\\lambda_{\\parallel} & 0 & 0\\\\\\\\\n",
      "0 & \\lambda_{\\perp} & 0\\\\\\\\\n",
      "0 & 0 & \\lambda_{\\perp}\n",
      "\\end{array}\\right) = \n",
      "10^{-3}\\left(\\begin{array}{ccc}\n",
      "1.4 & 0 & 0\\\\\\\\\n",
      "0 & 0.35 & 0\\\\\\\\\n",
      "0 & 0 & 0.35\n",
      "\\end{array}\\right) \\mathrm{mm}^2 \\cdot \\mathrm{s}^{-1}.$$\n",
      "\n",
      "Units:\n",
      "\n",
      "- b: $\\mathrm{s} / \\mathrm{mm}^2$\n",
      "- $\\Lambda$: $\\mathrm{mm}^2 / \\mathrm{s}$\n",
      "- Note that, with the above units, the exponent is dimensionless.\n",
      "- $1\\ \\mathrm{T} = \\frac{\\mathrm{N} \\cdot \\mathrm{s}}{\\mathrm{C} \\cdot \\mathrm{m}} =\n",
      "\\frac{\\mathrm{kg}}{\\mathrm{C}\\mathrm{s}}$\n",
      "\n",
      "Also:\n",
      "\n",
      "- $b = \\gamma^2G \\delta^2(\\Delta - \\delta / 3)$\n",
      "- $\\gamma = \\frac{e}{2 m_p} g = 4257\\ \\mathrm{Hz}/\\mathrm{gauss}$ : Gyromagnetic ratio of the proton\n",
      "- $g \\approx 5.586$: g-factor of the nucleus (dimensionless number)\n",
      "- $m_p$: Mass of the proton\n",
      "- $G$: field gradient amplitude\n",
      "- $\\Delta$, $\\delta$: time constants related to the pulse sequence\n",
      "\n",
      "<img src=\"files/phantom/dmri_process.svg\" width=\"40%\"/>"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "from dipy.sims.phantom import orbital_phantom, add_rician_noise\n",
      "from dipy.viz import fvtk"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def f1(t):\n",
      "    t = np.linspace(-1, 1, len(t)) * 0.95\n",
      "    x = t\n",
      "    y = 0 * t\n",
      "    z = t ** 2 + 0.05\n",
      "    return x, y, z"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def f2(t):\n",
      "    t = np.linspace(-1, 1, len(t)) * 0.95\n",
      "    #x = (1 - np.cos(t)) ** 2\n",
      "    x = t ** 2 - 0.05\n",
      "    y = t\n",
      "    z = 0 * x\n",
      "    return x, y, z"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "def f3(t):\n",
      "    t = np.linspace(-1, 1, len(t)) * 0.95\n",
      "    x = -(t ** 2 + 0.05)\n",
      "    y = t\n",
      "    z = 0 * x\n",
      "    return x, y, z"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "t = np.linspace(0, 2 * np.pi, 100)\n",
      "data_shape = (33, 33, 33, 65)\n",
      "origin = np.array([16, 16, 16])\n",
      "scale = 0.9 * origin\n",
      "\n",
      "vol = orbital_phantom(func=f1, t=t, datashape=data_shape, origin=origin, scale=scale)\n",
      "vol += orbital_phantom(func=f2, t=t, datashape=data_shape, origin=origin, scale=scale)\n",
      "vol += orbital_phantom(func=f3, t=t, datashape=data_shape, origin=origin, scale=scale)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "r=fvtk.ren()\n",
      "fvtk.add(r, fvtk.volume(vol[..., 0]))\n",
      "fvtk.show(r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "vol.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "(33, 33, 33, 65)"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}