{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Tracking with Robust Unbiased Model-BAsed Spherical Deconvolution (RUMBA-SD)\n\nHere, we demonstrate fiber tracking using a probabilistic direction getter\nand RUMBA-SD, a model introduced in [CanalesRodriguez2015]_. This model adapts\nRichardson-Lucy deconvolution by assuming Rician or Noncentral Chi noise\ninstead of Gaussian, which more accurately reflects the noise from MRI\nscanners (see also `example_reconst_rumba`). This tracking tutorial is an\nextension on `example_tracking_probabilistic`.\n\nWe start by loading sample data and identifying a fiber response function.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from numpy.linalg import inv\n\nfrom dipy.core.gradients import gradient_table\nfrom dipy.data import get_fnames, small_sphere\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti, load_nifti_data\nfrom dipy.reconst.csdeconv import auto_response_ssst\nfrom dipy.tracking import utils\nfrom dipy.tracking.local_tracking import LocalTracking\nfrom dipy.tracking.streamline import Streamlines, transform_streamlines\nfrom dipy.tracking.stopping_criterion import ThresholdStoppingCriterion\nfrom dipy.viz import window, actor, colormap\nfrom dipy.reconst.rumba import RumbaSDModel\n\n# Enables/disables interactive visualization\ninteractive = False\n\nhardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')\nlabel_fname = get_fnames('stanford_labels')\nt1_fname = get_fnames('stanford_t1')\n\ndata, affine, hardi_img = load_nifti(hardi_fname, return_img=True)\nlabels = load_nifti_data(label_fname)\nt1_data, t1_aff = load_nifti(t1_fname)\nbvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)\n\nseed_mask = (labels == 2)\nwhite_matter = (labels == 1) | (labels == 2)\nseeds = utils.seeds_from_mask(seed_mask, affine, density=2)\n\nresponse, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)\n\nsphere = small_sphere"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now initialize a `RumbaSdModel` model and fit it globally by setting\n`voxelwise` to `False`. For this example, TV regularization (`use_tv`) will be\nturned off for efficiency, although its usage can provide more coherent results\nin fiber tracking. The fit will take about 5 minutes to complete.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rumba = RumbaSDModel(gtab, wm_response=response[0], n_iter=200,\n                     voxelwise=False, use_tv=False, sphere=sphere)\nrumba_fit = rumba.fit(data, mask=white_matter)\nodf = rumba_fit.odf()  # fODF\nf_wm = rumba_fit.f_wm  # white matter volume fractions"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To establish stopping criterion, a common technique is to use the Generalized\nFractional Anisotropy (GFA). One point of caution is that RUMBA-SD by default\nseparates the fODF from an isotropic compartment. This can bias GFA results\ncomputed on the fODF, although it will still generally be an effective\ncriterion.\n\nHowever, an alternative stopping criterion that takes advantage of this\nfeature is to use RUMBA-SD's white matter volume fraction map.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "stopping_criterion = ThresholdStoppingCriterion(f_wm, .25)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can visualize a slice of this mask.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n\nsli = f_wm.shape[2] // 2\nplt.figure()\n\nplt.subplot(1, 2, 1).set_axis_off()\nplt.imshow(f_wm[:, :, sli].T, cmap='gray', origin='lower')\n\nplt.subplot(1, 2, 2).set_axis_off()\nplt.imshow((f_wm[:, :, sli] > 0.25).T, cmap='gray', origin='lower')\n\nplt.savefig('f_wm_tracking_mask.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: f_wm_tracking_mask.png\n   :align: center\n\n   White matter volume fraction slice\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "These discrete fODFs can be used as a PMF in the `ProbabilisticDirectionGetter`\nfor sampling tracking directions. The PMF must be strictly non-negative;\nRUMBA-SD already adheres to this constraint so no further manipulation of the\nfODFs is necessary.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.direction import ProbabilisticDirectionGetter\nfrom dipy.io.stateful_tractogram import Space, StatefulTractogram\nfrom dipy.io.streamline import save_trk\n\nprob_dg = ProbabilisticDirectionGetter.from_pmf(odf, max_angle=30.,\n                                                sphere=sphere)\nstreamline_generator = LocalTracking(prob_dg, stopping_criterion, seeds,\n                                     affine, step_size=.5)\nstreamlines = Streamlines(streamline_generator)\n\ncolor = colormap.line_colors(streamlines)\nstreamlines_actor = actor.streamtube(\n    list(transform_streamlines(streamlines, inv(t1_aff))),\n    color, linewidth=0.1)\n\nvol_actor = actor.slicer(t1_data)\nvol_actor.display(x=40)\nvol_actor2 = vol_actor.copy()\nvol_actor2.display(z=35)\n\nscene = window.Scene()\nscene.add(vol_actor)\nscene.add(vol_actor2)\nscene.add(streamlines_actor)\nif interactive:\n    window.show(scene)\n\nwindow.record(scene, out_path='tractogram_probabilistic_rumba.png',\n              size=(800, 800))\n\nsft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)\nsave_trk(sft, \"tractogram_probabilistic_rumba.trk\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractogram_probabilistic_rumba.png\n   :align: center\n\n   RUMBA-SD tractogram\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n\n.. [CanalesRodriguez2015] Canales-Rodr\u00edguez, E. J., Daducci, A., Sotiropoulos,\n   S. N., Caruyer, E., Aja-Fern\u00e1ndez, S., Radua, J., Mendizabal, J. M. Y.,\n   Iturria-Medina, Y., Melie-Garc\u00eda, L., Alem\u00e1n-G\u00f3mez, Y., Thiran, J.-P.,\n   Sarr\u00f3, S., Pomarol-Clotet, E., & Salvador, R. (2015). Spherical\n   Deconvolution of Multichannel Diffusion MRI Data with Non-Gaussian Noise\n   Models and Spatial Regularization. PLOS ONE, 10(10), e0138910.\n   https://doi.org/10.1371/journal.pone.0138910\n\n.. include:: ../links_names.inc\n\n\n"
      ]
    }
  ],
  "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.9.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}