s-track-figure6.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 dsiadapt as dsi\n",
    "from dipy.viz import fvtk\n",
    "import dipy.core.gradients as grad\n",
    "from dipy.data import get_sphere\n",
    "from dipy.tracking.local import LocalTracking\n",
    "from dipy.viz import fvtk\n",
    "from dipy.viz.colormap import line_colors\n",
    "from dipy.tracking import utils\n",
    "from dipy.tracking.local import ThresholdTissueClassifier\n",
    "from dipy.direction import peaks_from_model\n",
    "import time\n",
    "from dipy.io.trackvis import save_trk\n",
    "from dipy.io.pickles import save_pickle\n",
    "from dipy.io.pickles import load_pickle\n",
    "import glob\n",
    "from dipy.io.trackvis import save_trk\n",
    "\n",
    "%pylab inline\n",
    "np.set_printoptions(threshold=numpy.nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# ODF sphere directions\n",
    "sphere = get_sphere('symmetric724')\n",
    "#sphere = sphere.subdivide(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1633.97955608\n",
      "1684.98334193\n",
      "1687.14775681\n",
      "1620.3402009\n",
      "1264.95824313\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/anaconda/lib/python2.7/site-packages/dipy/viz/colormap.py:243: RuntimeWarning: invalid value encountered in divide\n",
      "  orient=np.abs(orient/np.linalg.norm(orient))\n"
     ]
    }
   ],
   "source": [
    "datapath = '/data/qytian/QTResearchResult/DSIQspace/DSIData'; # Whole brain data is only locally available \n",
    "wintypeArr = np.array(['none', 'hanning', 'hamming', 'blackman', 'none']);\n",
    "qgridsz = 17;\n",
    "\n",
    "for ii in np.arange(fnArr.shape[0]):\n",
    "    fn = fnArr[ii];\n",
    "    data = nib.load(datapath + '/' + fn + '/' + fn + '.nii.gz').get_data();\n",
    "    affine = np.identity(4);\n",
    "    wmmask = nib.load(datapath + '/' + fn + '/' + fn + '_wmmask.nii.gz').get_data();\n",
    "    brainmask = nib.load(datapath + '/' + fn + '/' + fn + '_masked_mask.nii.gz').get_data();\n",
    "    classifier = ThresholdTissueClassifier(brainmask, 0.5)\n",
    "\n",
    "    gtab = grad.gradient_table('data/' + fn + '_bvals.txt', 'data/' + fn + '_bvecs_dipy.txt');\n",
    "    seedmask = nib.load('data/' + fn + '_roi_mask.nii.gz').get_data();\n",
    "\n",
    "    mdd = np.loadtxt('data/' + fn + '_stats.txt')[0];\n",
    "    fov = np.loadtxt('data/' + fn + '_stats.txt')[1];\n",
    "    \n",
    "    intdistArr = np.array([fov/2, fov/2, fov/2, fov/2, mdd]); \n",
    "    \n",
    "    for jj in np.arange(wintypeArr.shape[0]):\n",
    "        wintype = wintypeArr[jj];\n",
    "        intdist = intdistArr[jj];   \n",
    "        rend = intdist / fov * qgridsz;\n",
    "        winwidth = 2. * dsi.create_qtable(gtab).max();\n",
    "        rend = intdist / fov * qgridsz;\n",
    "\n",
    "        dsimodel = dsi.DiffusionSpectrumModel(gtab, qgrid_size=qgridsz,\n",
    "                    filter_width=winwidth, filter_type=wintype, r_start=0, r_end=rend, r_step=0.2); \n",
    "        \n",
    "        peakname = 'figure6/' + fn + '_peak_qgridsz' + str(qgridsz) + '_odfdirnum' + str(sphere.vertices.shape[0]) + '_' + wintype + '_rend' + str(intdist) + '.pkl';\n",
    "        \n",
    "        if len(glob.glob(peakname)) == 0:\n",
    "            # Compute ODF peaks\n",
    "            tic = time.time()\n",
    "            dsipeaks = peaks_from_model(dsimodel, data, sphere,\n",
    "                         relative_peak_threshold=0.2,\n",
    "                         min_separation_angle=25, normalize_peaks=True, mask=brainmask);\n",
    "            toc = time.time() - tic\n",
    "            disp(toc)\n",
    "\n",
    "            save_pickle(peakname, dsipeaks);\n",
    "        else:\n",
    "            dsipeaks = load_pickle(peakname);\n",
    "        \n",
    "        # Track\n",
    "        dennum = 2;\n",
    "        seeds = utils.seeds_from_mask(seedmask, density=[dennum, dennum, dennum], affine=affine)\n",
    "\n",
    "        # Initialization of LocalTracking. The computation happens in the next step.\n",
    "        streamlines = LocalTracking(dsipeaks, classifier, seeds, affine, step_size=.5)\n",
    "\n",
    "        # Compute streamlines and store as a list.\n",
    "        streamlines = list(streamlines)\n",
    "\n",
    "        # Prepare the display objects.\n",
    "        color = line_colors(streamlines)\n",
    "        streamlines_actor = fvtk.line(streamlines, line_colors(streamlines))\n",
    "        \n",
    "        trackname = 'figure6/' + fn + '_track_qgridsz' + str(qgridsz) + '_odfdirnum' + str(sphere.vertices.shape[0]) + '_' + wintype + '_rend' + str(intdist) + '_density' + str(dennum) + '.trk';\n",
    "        save_trk(trackname, streamlines, affine, wmmask.shape)\n",
    "\n",
    "        # Create the 3d display.\n",
    "        # r = fvtk.ren()\n",
    "        # fvtk.add(r, streamlines_actor)\n",
    "        # fvtk.camera(r, pos=(0, 1, 0), viewup=(0, 0, 1))\n",
    "        # fvtk.show(r)"
   ]
  }
 ],
 "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
}