{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Fiber to bundle coherence measures\n\nThis demo presents the fiber to bundle coherence (FBC) quantitative\nmeasure of the alignment of each fiber with the surrounding fiber bundles\n[Meesters2016]_. These measures are useful in 'cleaning' the results of\ntractography algorithms, since low FBCs indicate which fibers are isolated and\npoorly aligned with their neighbors, as shown in the figure below.\n\n\n.. figure:: _static/fbc_illustration.png\n   :scale: 60 %\n   :align: center\n\n   On the left this figure illustrates (in 2D) the contribution of two fiber\n   points to the kernel density estimator. The kernel density estimator is the\n   sum over all such locally aligned kernels. The local fiber to bundle\n   coherence, shown on the right, color-coded for each fiber, is obtained by\n   evaluating the kernel density estimator along the fibers. One spurious\n   fiber is present which is isolated and badly aligned with the other fibers,\n   and can be identified by a low LFBC value in the region where it deviates\n   from the bundle. Figure adapted from [Portegies2015]_.\n\nHere we implement FBC measures based on kernel density estimation in the\nnon-flat 5D position-orientation domain. First we compute the kernel density\nestimator induced by the full lifted output (defined in the space of positions\nand orientations) of the tractography. Then, the Local FBC (LFBC) is the\nresult of evaluating the estimator along each element of the lifted fiber.\nA whole fiber measure, the relative FBC (RFBC), is calculated\nby the minimum of the moving average LFBC along the fiber.\nDetails of the computation of FBC can be found in [Portegies2015]_.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The FBC measures are evaluated on the Stanford HARDI dataset\n(150 orientations, b=2000 $s/mm^2$) which is one of the standard example\ndatasets in DIPY_.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Enables/disables interactive visualization\ninteractive = False\n\nimport numpy as np\nfrom dipy.core.gradients import gradient_table\nfrom dipy.data import get_fnames\nfrom dipy.io.image import load_nifti_data, load_nifti\nfrom dipy.io.gradients import read_bvals_bvecs\n\n# Fix seed\nnp.random.seed(1)\n\n# Read data\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 = load_nifti(hardi_fname)\nlabels = load_nifti_data(label_fname)\nt1_data = load_nifti_data(t1_fname)\nbvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)\n\n\n\n# Select a relevant part of the data (left hemisphere)\n# Coordinates given in x bounds, y bounds, z bounds\ndshape = data.shape[:-1]\nxa, xb, ya, yb, za, zb = [15, 42, 10, 65, 18, 65]\ndata_small = data[xa:xb, ya:yb, za:zb]\nselectionmask = np.zeros(dshape, 'bool')\nselectionmask[xa:xb, ya:yb, za:zb] = True"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The data is first fitted to the Constant Solid Angle (CDA) ODF Model. CSA is a\ngood choice to estimate general fractional anisotropy (GFA), which the stopping\ncriterion can use to restrict fiber tracking to those areas where the ODF\nshows significant restricted diffusion, thus creating a region-of-interest in\nwhich the computations are done.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Perform CSA\nfrom dipy.reconst.shm import CsaOdfModel\nfrom dipy.data import default_sphere\nfrom dipy.direction import peaks_from_model\n\ncsa_model = CsaOdfModel(gtab, sh_order=6)\ncsa_peaks = peaks_from_model(csa_model, data, default_sphere,\n                             relative_peak_threshold=.6,\n                             min_separation_angle=45,\n                             mask=selectionmask)\n\n# Stopping Criterion\nfrom dipy.tracking.stopping_criterion import ThresholdStoppingCriterion\n\nstopping_criterion = ThresholdStoppingCriterion(csa_peaks.gfa, 0.25)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In order to perform probabilistic fiber tracking we first fit the data to the\nConstrained Spherical Deconvolution (CSD) model in DIPY. This model represents\neach voxel in the data set as a collection of small white matter fibers with\ndifferent orientations. The density of fibers along each orientation is known\nas the Fiber Orientation Distribution (FOD), used in the fiber tracking.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Perform CSD on the original data\nfrom dipy.reconst.csdeconv import auto_response_ssst\nfrom dipy.reconst.csdeconv import ConstrainedSphericalDeconvModel\n\nresponse, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)\ncsd_model = ConstrainedSphericalDeconvModel(gtab, response)\ncsd_fit = csd_model.fit(data_small)\ncsd_fit_shm = np.lib.pad(csd_fit.shm_coeff, ((xa, dshape[0]-xb),\n                                             (ya, dshape[1]-yb),\n                                             (za, dshape[2]-zb),\n                                             (0, 0)), 'constant')\n\n# Probabilistic direction getting for fiber tracking\nfrom dipy.direction import ProbabilisticDirectionGetter\n\nprob_dg = ProbabilisticDirectionGetter.from_shcoeff(csd_fit_shm,\n                                                    max_angle=30.,\n                                                    sphere=default_sphere)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The optic radiation is reconstructed by tracking fibers from the calcarine\nsulcus (visual cortex V1) to the lateral geniculate nucleus (LGN). We seed\nfrom the calcarine sulcus by selecting a region-of-interest (ROI) cube of\ndimensions 3x3x3 voxels.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set a seed region region for tractography.\nfrom dipy.tracking import utils\n\nmask = np.zeros(data.shape[:-1], 'bool')\nrad = 3\nmask[26-rad:26+rad, 29-rad:29+rad, 31-rad:31+rad] = True\nseeds = utils.seeds_from_mask(mask, affine, density=[4, 4, 4])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Local Tracking is used for probabilistic tractography which takes the\ndirection getter along with the stopping criterion and seeds as input.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Perform tracking using Local Tracking\nfrom dipy.tracking.local_tracking import LocalTracking\n\nstreamlines_generator = LocalTracking(prob_dg, stopping_criterion, seeds,\n                                      affine, step_size=.5)\n\n# Compute streamlines.\nfrom dipy.tracking.streamline import Streamlines\nstreamlines = Streamlines(streamlines_generator)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In order to select only the fibers that enter into the LGN, another ROI is\ncreated from a cube of size 5x5x5 voxels. The near_roi command is used to find\nthe fibers that traverse through this ROI.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set a mask for the lateral geniculate nucleus (LGN)\nmask_lgn = np.zeros(data.shape[:-1], 'bool')\nrad = 5\nmask_lgn[35-rad:35+rad, 42-rad:42+rad, 28-rad:28+rad] = True\n\n# Select all the fibers that enter the LGN and discard all others\nfiltered_fibers2 = utils.near_roi(streamlines, affine, mask_lgn, tol=1.8)\n\nsfil = []\nfor i in range(len(streamlines)):\n    if filtered_fibers2[i]:\n        sfil.append(streamlines[i])\nstreamlines = Streamlines(sfil)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Inspired by [Rodrigues2010]_, a lookup-table is created, containing rotated\nversions of the fiber propagation kernel $P_t$ [DuitsAndFranken2011]_\nrotated over a discrete set of orientations. See the\n[Contextual enhancement example](https://dipy.org/documentation/latest/examples_built/contextual_enhancement/#example-contextual-enhancement)\nfor more details regarding the kernel. In order to ensure rotationally\ninvariant processing, the discrete orientations are required to be equally\ndistributed over a sphere. By default, a sphere with 100 directions is obtained\nfrom electrostatic repulsion in DIPY.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Compute lookup table\nfrom dipy.denoise.enhancement_kernel import EnhancementKernel\n\nD33 = 1.0\nD44 = 0.02\nt = 1\nk = EnhancementKernel(D33, D44, t)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The FBC measures are now computed, taking the tractography results and the\nlookup tables as input.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Apply FBC measures\nfrom dipy.tracking.fbcmeasures import FBCMeasures\n\nfbc = FBCMeasures(streamlines, k)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After calculating the FBC measures, a threshold can be chosen on the relative\nFBC (RFBC) in order to remove spurious fibers. Recall that the relative FBC\n(RFBC) is calculated by the minimum of the moving average LFBC along the fiber.\nIn this example we show the results for threshold 0 (i.e. all fibers are\nincluded) and 0.2 (removing the 20 percent most spurious fibers).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Calculate LFBC for original fibers\nfbc_sl_orig, clrs_orig, rfbc_orig = \\\n  fbc.get_points_rfbc_thresholded(0, emphasis=0.01)\n\n# Apply a threshold on the RFBC to remove spurious fibers\nfbc_sl_thres, clrs_thres, rfbc_thres = \\\n  fbc.get_points_rfbc_thresholded(0.125, emphasis=0.01)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The results of FBC measures are visualized, showing the original fibers\ncolored by LFBC (see `optic_radiation_before_cleaning`), and the fibers\nafter the cleaning procedure via RFBC thresholding (see\n`optic_radiation_after_cleaning`).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Visualize the results\nfrom dipy.viz import window, actor\n\n# Create scene\nscene = window.Scene()\n\n# Original lines colored by LFBC\nlineactor = actor.line(fbc_sl_orig, np.vstack(clrs_orig), linewidth=0.2)\nscene.add(lineactor)\n\n# Horizontal (axial) slice of T1 data\nvol_actor1 = actor.slicer(t1_data, affine=affine)\nvol_actor1.display(z=20)\nscene.add(vol_actor1)\n\n# Vertical (sagittal) slice of T1 data\nvol_actor2 = actor.slicer(t1_data, affine=affine)\nvol_actor2.display(x=35)\nscene.add(vol_actor2)\n\n# Show original fibers\nscene.set_camera(position=(-264, 285, 155),\n                 focal_point=(0, -14, 9),\n                 view_up=(0, 0, 1))\nwindow.record(scene, n_frames=1, out_path='OR_before.png', size=(900, 900))\nif interactive:\n    window.show(scene)\n\n# Show thresholded fibers\nscene.rm(lineactor)\nscene.add(actor.line(fbc_sl_thres, np.vstack(clrs_thres), linewidth=0.2))\nwindow.record(scene, n_frames=1, out_path='OR_after.png', size=(900, 900))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n.. figure:: OR_before.png\n   :align: center\n\n   The optic radiation obtained through probabilistic tractography colored by\n   local fiber to bundle coherence.\n\n\n.. figure:: OR_after.png\n   :align: center\n\n   The tractography result is cleaned (shown in bottom) by removing fibers\n   with a relative FBC (RFBC) lower than the threshold $\\tau = 0.2$.\n\n## Acknowledgments\nThe techniques are developed in close collaboration with Pauly Ossenblok of\nthe Academic Center of Epileptology Kempenhaeghe & Maastricht UMC+.\n\n## References\n\n.. [Meesters2016] S. Meesters, G. Sanguinetti, E. Garyfallidis, J. Portegies,\n   P. Ossenblok, R. Duits. (2016) Cleaning output of tractography via fiber to\n   bundle coherence, a new open source implementation. Human Brain Mapping\n   Conference 2016.\n\n.. [Portegies2015] J. Portegies, R. Fick, G. Sanguinetti, S. Meesters,\n   G.Girard, and R. Duits. (2015) Improving Fiber Alignment in HARDI by\n   Combining Contextual PDE flow with Constrained Spherical Deconvolution. PLoS\n   One.\n\n.. [DuitsAndFranken2011] R. Duits and E. Franken (2011) Left-invariant\n   diffusions on the space of positions and orientations and their application\n   to crossing-preserving smoothing of HARDI images. International Journal of\n   Computer Vision, 92:231-264.\n\n.. [Rodrigues2010] P. Rodrigues, R. Duits, B. Romeny, A. Vilanova (2010).\n   Accelerated Diffusion Operators for Enhancing DW-MRI. Eurographics Workshop\n   on Visual Computing for Biology and Medicine. The Eurographics Association.\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
}