Virtual_Lesion_Step2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Virtual Lesion\n",
    "This series of notebooks is used to determine the evidence for a pathway connecting two ROIs following those described in Pestilli et al.<sup>[1]</sup> and Leong et al.<sup>[2]</sup>.\n",
    "\n",
    "Calculation steps done in this notebook:\n",
    "\n",
    "<ol>\n",
    "    <li>Calculate the neighborhood of the candidate streamlines (Voxels where the streamlines pass through)</li>\n",
    "    <li>Find fibers in the whole brain connectome that pass through the neibghorhood</li>\n",
    "    <li>Combine whole brain connectome with candidate streamlines into one steamline set</li>\n",
    "    <li>Run the LiFE optimization on this combined streamline set and retrieve RMSE between the original data set and the LiFE prediction</li>\n",
    "    <li>Run the LiFE otimization and RMSE calculation again on the whole brain, but this time without the candidate streamlines</li>\n",
    "</ol>\n",
    "\n",
    "**<small>Do not run the parallization on all cores. This notebook is parallalized by subjects. The limiting factor is the memory usage per process.</small>\n",
    "\n",
    "<sup>[1]</sup> <i>PMID: 25194848</i> <br/>\n",
    "<sup>[2]</sup> <i>PMID: 26748088</i> <br/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from dipy.tracking.vox2track import streamline_mapping\n",
    "import os\n",
    "import numpy as np\n",
    "from nibabel import trackvis as tv\n",
    "from utilities import *\n",
    "import nibabel as nib\n",
    "import dipy.tracking.life as life\n",
    "from dipy.tracking import utils\n",
    "import copy\n",
    "from parallelization import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Library of Files\n",
    "path = '/hcp/'\n",
    "path_saveing = '/data/hcp/data/'\n",
    "\n",
    "subjects = os.listdir(path_saveing)\n",
    "\n",
    "subjects_sorted = sorted(subjects)\n",
    "subjects_sorted.remove('.nii.gz')  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_LiFE(subject):\n",
    "    print 'Process subject ' + subject\n",
    "\n",
    "    if os.path.isfile(os.path.join(path_saveing, subject, 'RVLPFC2FIRSTamyg_bigRight_LIFE_started2.txt')) == False:\n",
    "        print \"LiFE Files do not exist for this subject, start calculation.\"\n",
    "\n",
    "        if os.path.isfile(os.path.join(path_saveing, subject, 'RVLPFC2FIRSTamyg_bigRight_clustered.trk')) == True and os.path.isfile(os.path.join(path_saveing, subject, '2M_SIFT.trk')) == True:\n",
    "            print \"All neccessary files there, continue ...\"\n",
    "\n",
    "            print \"Show other processes that this subject is processed\"\n",
    "            done = np.array([1])\n",
    "            np.savetxt(os.path.join(path_saveing, subject, 'RVLPFC2FIRSTamyg_bigRight_LIFE_started2.txt'), done, delimiter=',')\n",
    "            \n",
    "            try: \n",
    "                directory_output = os.path.join(path_saveing, subject)\n",
    "\n",
    "                print \"Start calculation for subject %s\" % subject\n",
    "                f_streamlines = os.path.join(path_saveing, subject, 'RVLPFC2FIRSTamyg_bigRight_clustered.trk')\n",
    "                f_in_nifti = os.path.join(path, subject, 'T1w/Diffusion/data.nii.gz')\n",
    "\n",
    "                streams, hdr = tv.read(f_streamlines, points_space='voxel')\n",
    "                streamlines = [i[0] for i in streams]\n",
    "\n",
    "                data, affine, gtab, header, shell_mask = load_hcp_data(path, subject)\n",
    "                dim = header['dim'][1:4]\n",
    "                \n",
    "                # Otherwise all weights are NaN\n",
    "                data[data <= 0.0] = 1.0\n",
    "\n",
    "                print \"Calculating neighborhood with LiFE\"\n",
    "                fiber_model_neighborhood = life.FiberModel(gtab)\n",
    "                fiber_fit_neighborhood = fiber_model_neighborhood.fit(data, streamlines, affine=np.eye(4))\n",
    "                indices_neighborhood = fiber_fit_neighborhood.vox_coords\n",
    "\n",
    "                neighborhood  = np.zeros(dim, dtype=bool)\n",
    "                for i in range(indices_neighborhood.shape[0]):\n",
    "                    neighborhood[indices_neighborhood[i][0], indices_neighborhood[i][1], indices_neighborhood[i][2]] = 1\n",
    "\n",
    "                save_as_nifti(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_neighborhood\", neighborhood.astype(np.int), affine)\n",
    "                \n",
    "                print 'Find fiber that pass through neighborhood'\n",
    "                f_streamlines_whole_brain = path_saveing + subject + \"/2M_SIFT.trk\"\n",
    "                streams_whole_brain, hdr_whole_brain = tv.read(f_streamlines_whole_brain, points_space='voxel')\n",
    "                streamlines_whole_brain = [i[0] for i in streams_whole_brain]\n",
    "                neighborhood_streamlines = utils.target(streamlines_whole_brain, neighborhood, affine=np.eye(4))\n",
    "\n",
    "                neighborhood_streamlines = list(neighborhood_streamlines)\n",
    "\n",
    "                strm = ((sl, None, None) for sl in neighborhood_streamlines)\n",
    "                tv.write(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_2M_SIFT_without_path.trk\", strm,  hdr_mapping=hdr_whole_brain, points_space='voxel')\n",
    "\n",
    "                print \"Combine streamlines\"\n",
    "                streamlines_together = neighborhood_streamlines + streamlines\n",
    "                strm_together = ((sl, None, None) for sl in streamlines_together)\n",
    "                tv.write(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_2M_SIFT_with_path.trk\", strm_together,  hdr_mapping=hdr_whole_brain, points_space='voxel')\n",
    "\n",
    "                print \"Start LiFE optimization with new path\"\n",
    "                fiber_model_together = life.FiberModel(gtab)\n",
    "                fiber_fit_together = fiber_model_together.fit(data, streamlines_together, affine=np.eye(4))\n",
    "                model_predict_together = fiber_fit_together.predict()\n",
    "                indices_together = fiber_fit_together.vox_coords\n",
    "\n",
    "                mask_with  = np.zeros(dim, dtype=bool)\n",
    "                whole_brain_together = np.zeros(header['dim'][1:5])\n",
    "                for i in range(indices_together.shape[0]):\n",
    "                    whole_brain_together[indices_together[i][0], indices_together[i][1], indices_together[i][2]] = model_predict_together[i]\n",
    "                    mask_with[indices_together[i][0], indices_together[i][1], indices_together[i][2]] = 1\n",
    "\n",
    "                save_as_nifti(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_LiFE_prediction_with_path\", whole_brain_together, affine)\n",
    "                save_as_matlab_file(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_LiFE_betas_with_path\", beta = fiber_fit_together.beta)\n",
    "                save_as_nifti(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_LiFE_mask_with_path\", mask_with.astype(np.int), affine)\n",
    "\n",
    "                print \"Calculate RMSE with\"\n",
    "                model_error_together = model_predict_together - fiber_fit_together.data\n",
    "                model_rmse_together = np.sqrt(np.mean(model_error_together[..., ~gtab.b0s_mask] ** 2, -1))\n",
    "\n",
    "                whole_brain_rmse_together = np.zeros(dim)\n",
    "                for i in range(indices_together.shape[0]):\n",
    "                    whole_brain_rmse_together[indices_together[i][0], indices_together[i][1], indices_together[i][2]] = model_rmse_together[i]\n",
    "\n",
    "                save_as_nifti(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_LiFE_rmse_with_path\", whole_brain_rmse_together, affine)\n",
    "                \n",
    "                print \"Start LiFE optimization without new path\"\n",
    "                fiber_fit = copy.deepcopy(fiber_fit_together)\n",
    "                fiber_fit.beta[-len(streamlines):] = 0\n",
    "                model_predict = fiber_fit.predict()\n",
    "                indices = fiber_fit.vox_coords\n",
    "\n",
    "                whole_brain = np.zeros(header['dim'][1:5])\n",
    "                mask_without  = np.zeros(dim, dtype=bool)\n",
    "                for i in range(indices.shape[0]):\n",
    "                    whole_brain[indices[i][0], indices[i][1], indices[i][2]] = model_predict[i]\n",
    "                    mask_without[indices[i][0], indices[i][1], indices[i][2]] = 1\n",
    "\n",
    "                save_as_nifti(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_LiFE_prediction_without_path\", whole_brain, affine)\n",
    "                save_as_matlab_file(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_LiFE_betas_without_path\", beta = fiber_fit.beta)\n",
    "                save_as_nifti(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_LiFE_mask_without_path\", mask_without.astype(np.int), affine)\n",
    "\n",
    "                print \"Calculate RMSE without\"\n",
    "                model_error = model_predict - fiber_fit.data\n",
    "                model_rmse = np.sqrt(np.mean(model_error[..., ~gtab.b0s_mask] ** 2, -1))\n",
    "\n",
    "                whole_brain_rmse = np.zeros(dim)\n",
    "                for i in range(indices.shape[0]):\n",
    "                    whole_brain_rmse[indices[i][0], indices[i][1], indices[i][2]] = model_rmse[i]\n",
    "\n",
    "                save_as_nifti(path_saveing + subject + \"/RVLPFC2FIRSTamyg_bigRight_LiFE_rmse_without_path\", whole_brain_rmse, affine)\n",
    "                \n",
    "                print \"All done\"\n",
    "            except:\n",
    "                print \"An error occured while computing LiFE. Skip this subject.\"\n",
    "        else:\n",
    "            print \"Some input files are missing, skip this subject.\"\n",
    "    else:\n",
    "        print \"LiFE Files exist already for this subject, skip calculation.\"\n",
    "    \n",
    "    return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "subjects_for_life = subjects_sorted\n",
    "p = parallelization(maximum_number_of_cores=7, display=False)\n",
    "return_values = p.start(run_LiFE, len(subjects_for_life), subjects_for_life)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}