s-stats.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import nibabel as nib\n",
    "import dipy.reconst.dti as dti\n",
    "import dsiadapt as dsi\n",
    "import dipy.core.gradients as grad\n",
    "import matplotlib \n",
    "%pylab inline\n",
    "np.set_printoptions(threshold=numpy.nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "123\n",
      "4250.0\n",
      "57\n",
      "2000.0\n",
      "81\n",
      "1680.0\n"
     ]
    }
   ],
   "source": [
    "# Estimate diffusion rates using DTI fitting\n",
    "\n",
    "fnArr = np.array(['DSI17_exvivo', 'DSI11_invivo_b10k', 'DSI11_invivo_b7k'])\n",
    "bmaxthrArr = np.array([4500, 2000, 2000]);\n",
    "adArr = np.zeros(fnArr.shape[0]); # axial diffusion rate\n",
    "\n",
    "for ii in np.arange(fnArr.shape[0]):\n",
    "    fn = fnArr[ii];\n",
    "    data = nib.load('data/' + fn + '_cc.nii.gz').get_data();\n",
    "    gtab = grad.gradient_table('data/' + fn + '_bvals.txt', 'data/' + fn + '_bvecs.txt');\n",
    "    mask = gtab.bvals <= bmaxthrArr[ii]; # Use only data with bvals less than 3000\n",
    "    datamask = data[:, :, :, mask]\n",
    "    gtabmask = grad.gradient_table_from_bvals_bvecs(gtab.bvals[mask], gtab.bvecs[mask, :])\n",
    "    disp(gtabmask.bvals.shape) # Number of measurements used for DTI fitting\n",
    "    disp(gtabmask.bvals.max())\n",
    "    \n",
    "    tensormodel = dti.TensorModel(gtabmask)\n",
    "    tensorfit = tensormodel.fit(datamask)\n",
    "    adArr[ii] = tensorfit.ad.mean()\n",
    "\n",
    "adArr = np.concatenate((adArr[0:1], adArr[0:1], adArr)) # Use the ad from ex vivo DSI17 for ex vivo DSI15 and ex vivo DSI11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.00017951  0.00017951  0.00017951  0.00142466  0.00163717]\n"
     ]
    }
   ],
   "source": [
    "disp(adArr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[[[-0.96167792,  0.20693869, -0.17986648],\n",
       "          [ 0.26721999,  0.56050859, -0.78385177],\n",
       "          [ 0.06139255,  0.80187686,  0.5943269 ]],\n",
       "\n",
       "         [[-0.97242061,  0.13196121, -0.1923133 ],\n",
       "          [ 0.23205059,  0.6303544 , -0.74081432],\n",
       "          [-0.02346678,  0.76500952,  0.64359128]]]],\n",
       "\n",
       "\n",
       "\n",
       "       [[[[-0.98296716,  0.16020355, -0.0900577 ],\n",
       "          [ 0.15738741,  0.4807815 , -0.86259977],\n",
       "          [ 0.09489347,  0.86208119,  0.49780644]],\n",
       "\n",
       "         [[-0.99204213,  0.09891501, -0.0778989 ],\n",
       "          [ 0.11463644,  0.45375531, -0.88372202],\n",
       "          [ 0.05206634,  0.88561952,  0.46148364]]]],\n",
       "\n",
       "\n",
       "\n",
       "       [[[[-0.99974284, -0.02043197,  0.00983773],\n",
       "          [-0.01676935,  0.37407194, -0.92724806],\n",
       "          [-0.01526549,  0.92717459,  0.37431837]],\n",
       "\n",
       "         [[-0.99896407,  0.04293378, -0.01508249],\n",
       "          [ 0.03005452,  0.37360834, -0.92709953],\n",
       "          [ 0.03416894,  0.92659241,  0.37451167]]]],\n",
       "\n",
       "\n",
       "\n",
       "       [[[[-0.98087337,  0.10256257, -0.16543385],\n",
       "          [-0.1610008 ,  0.05014837,  0.9856794 ],\n",
       "          [ 0.10939005,  0.99346165, -0.03267654]],\n",
       "\n",
       "         [[-0.95518342,  0.17954174, -0.23534952],\n",
       "          [-0.22104474,  0.09619154,  0.97050833],\n",
       "          [ 0.19688539,  0.97903624, -0.05219381]]]]])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tensorfit.evecs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Estimate mean diffusion distance\n",
    "fnArr = np.array(['DSI11_exvivo', 'DSI15_exvivo', 'DSI17_exvivo', 'DSI11_invivo_b10k', 'DSI11_invivo_b7k'])\n",
    "DELTAArr = np.array([29.4 * 10**-3, 29.4 * 10**-3, 29.4 * 10**-3, 20.9 * 10**-3, 49.2 * 10**-3]); # in second\n",
    "deltaArr = np.array([16.7 * 10**-3, 16.7 * 10**-3, 16.7 * 10**-3, 12.9 * 10**-3, 42.3 * 10**-3]); # in second\n",
    "fovArr = np.zeros(fnArr.shape[0]);\n",
    "mddArr = np.zeros(fnArr.shape[0]);\n",
    "\n",
    "for ii in np.arange(fnArr.shape[0]):\n",
    "    fn = fnArr[ii];\n",
    "    gtab = grad.gradient_table('data/' + fn + '_bvals.txt', 'data/' + fn + '_bvecs.txt');\n",
    "    DELTA = DELTAArr[ii]; # second\n",
    "    delta = deltaArr[ii]; # second\n",
    "    ad = adArr[ii]; # mm2/s\n",
    "    \n",
    "    # Compute FOV\n",
    "    bmax = gtab.bvals.max(); # s/mm2\n",
    "    qmax = np.sqrt(bmax / (DELTA - delta / 3)) / (2 * np.pi); # mm-1\n",
    "    deltaq = qmax / dsi.create_qtable(gtab).max(); # mm-1\n",
    "    fovArr[ii] = 1 / deltaq; # mm\n",
    "    \n",
    "    # Compute MDD\n",
    "    mddArr[ii] = np.sqrt(6 * ad * (DELTA - delta / 3)); # mm\n",
    "    \n",
    "    np.savetxt('data/' + fn + '_stats.txt', (mddArr[ii], fovArr[ii]));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.00506661  0.00506661  0.00506661  0.01191204  0.01856845]\n",
      "[ 0.02797822  0.0391695   0.04476515  0.04047659  0.07034843]\n"
     ]
    }
   ],
   "source": [
    "disp(mddArr)\n",
    "disp(fovArr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.          0.9682532   0.88987076  0.78094232  0.7204265   0.64146775\n",
      "  0.63224626  0.67720091  0.56481433  0.50314587  0.50718021  0.42072907\n",
      "  0.42533979  0.46971807  0.41611832  0.3290908   0.45127514]\n",
      "[     0.    450.    450.   1900.   1900.   4200.   4200.   7500.   7500.\n",
      "  11700.  11700.  16900.  16900.  22950.  22950.  30000.  30000.]\n",
      "[ 1.          0.82066344  0.92743607  0.61575674  0.65134762  0.46060815\n",
      "  0.4975812   0.38908086  0.40255701  0.2950933   0.34208708]\n",
      "[     0.    400.    400.   1600.   1600.   3600.   3600.   6400.   6400.\n",
      "  10000.  10000.]\n",
      "[ 1.          0.86217952  0.72875309  0.71066624  0.67684025  0.53699291\n",
      "  0.54359198  0.33684531  0.44082347  0.41415536  0.39524052]\n",
      "[    0.   280.   280.  1120.  1120.  2520.  2520.  4480.  4480.  7000.\n",
      "  7000.]\n"
     ]
    }
   ],
   "source": [
    "# Estimate signal at the q-space edge\n",
    "fnArr = np.array(['DSI17_exvivo', 'DSI11_invivo_b10k', 'DSI11_invivo_b7k'])\n",
    "bmaxthrArr = np.array([4500, 2000, 2000]);\n",
    "adArr = np.zeros(fnArr.shape[0]); # axial diffusion rate\n",
    "\n",
    "for ii in np.arange(fnArr.shape[0]):\n",
    "    fn = fnArr[ii];\n",
    "    data = nib.load('data/' + fn + '_cc.nii.gz').get_data();\n",
    "    gtab = grad.gradient_table('data/' + fn + '_bvals.txt', 'data/' + fn + '_bvecs.txt');\n",
    "    ori = np.array([0, 1, 0]);\n",
    "    mask = np.logical_or(np.dot(gtab.bvecs, ori) == 1, np.dot(gtab.bvecs, ori) == -1);\n",
    "    mask[0] = True\n",
    "    oribvals = gtab.bvals[mask];\n",
    "    oridata = data[:, :, :, mask] * 1.0;\n",
    "    oridatasum = oridata.sum(axis=0)\n",
    "    oridatasum = np.squeeze(oridatasum.sum(axis=1))\n",
    "    oridatasum = oridatasum / oridatasum[0];\n",
    "    disp(oridatasum)\n",
    "    disp(oribvals)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.32058978000000005"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(0.30603716 + 0.29278132 + 0.39364102 + 0.28989962) / 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0552868"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(0.06185211 + 0.04872149) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.112747245"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(0.10041282 + 0.12508167) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4165505825"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(0.46971807 + 0.41611832 + 0.3290908 + 0.45127514) / 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.31859019"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(0.2950933 + 0.34208708) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.40469794"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(0.41415536 + 0.39524052) / 2"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}