{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Automatic Fiber Bundle Extraction with RecoBundles\n\nThis example explains how we can use RecoBundles [Garyfallidis17]_ to extract\nbundles from tractograms.\n\nFirst import the necessary modules.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\nfrom dipy.align.streamlinear import whole_brain_slr\nfrom dipy.data import get_two_hcp842_bundles\nfrom dipy.data import (fetch_target_tractogram_hcp,\n                       fetch_bundle_atlas_hcp842,\n                       get_bundle_atlas_hcp842,\n                       get_target_tractogram_hcp)\nfrom dipy.io.stateful_tractogram import Space, StatefulTractogram\nfrom dipy.io.streamline import load_trk, save_trk\nfrom dipy.io.utils import create_tractogram_header\nfrom dipy.segment.bundles import RecoBundles\nfrom dipy.viz import actor, window"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Download and read data for this tutorial\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "target_file, target_folder = fetch_target_tractogram_hcp()\natlas_file, atlas_folder = fetch_bundle_atlas_hcp842()\n\natlas_file, all_bundles_files = get_bundle_atlas_hcp842()\ntarget_file = get_target_tractogram_hcp()\n\nsft_atlas = load_trk(atlas_file, \"same\", bbox_valid_check=False)\natlas = sft_atlas.streamlines\natlas_header = create_tractogram_header(atlas_file,\n                                        *sft_atlas.space_attributes)\n\nsft_target = load_trk(target_file, \"same\", bbox_valid_check=False)\ntarget = sft_target.streamlines\ntarget_header = create_tractogram_header(target_file,\n                                         *sft_target.space_attributes)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize atlas tractogram and target tractogram before registration\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(atlas, colors=(1, 0, 1)))\nscene.add(actor.line(target, colors=(1, 1, 0)))\nwindow.record(scene, out_path='tractograms_initial.png', size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractograms_initial.png\n   :align: center\n\n   Atlas and target before registration.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will register target tractogram to model atlas' space using streamlinear\nregistration (SLR) [Garyfallidis15]_\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "moved, transform, qb_centroids1, qb_centroids2 = whole_brain_slr(\n    atlas, target, x0='affine', verbose=True, progressive=True,\n    rng=np.random.RandomState(1984))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We save the transform generated in this registration, so that we can use\nit in the bundle profiles example\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "np.save(\"slr_transform.npy\", transform)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize atlas tractogram and target tractogram after registration\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(atlas, colors=(1, 0, 1)))\nscene.add(actor.line(moved, colors=(1, 1, 0)))\nwindow.record(scene, out_path='tractograms_after_registration.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: tractograms_after_registration.png\n   :align: center\n\n   Atlas and target after registration.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Extracting bundles using RecoBundles [Garyfallidis17]_\n\nRecoBundles requires a model (reference) bundle and tries to extract similar\nlooking bundle from the input tractogram. There are some key parameters that\nusers can set as per requirements. Here are the key threshold parameters\nmeasured in millimeters and their function in Recobundles:\n\n    - model_clust_thr : It will use QuickBundles to get the centroids of the\n    model bundle and work with centroids instead of all streamlines. This helps\n    to make RecoBundles faster. The larger the value of the threshold, the\n    fewer centroids will be, the and smaller the threshold value, the more\n    centroids will be. If you prefer to use all the streamlines of the model\n    bundle, you can set this threshold to 0.01 mm.\n    Recommended range of the model_clust_thr is 0.01 - 3.0 mm.\n\n    - reduction_thr : This threshold will be used to reduce the search space\n    for finding the streamlines that match model bundle streamlines in shape.\n    Instead of looking at the entire tractogram, now we will be looking at\n    neighboring region of a model bundle in the tractogram. Increase the\n    threshold to increase the search space.\n    Recommended range of the reduction_thr is 15 - 30 mm.\n\n    - pruning_thr : This threshold will filter the streamlines for which the\n    distance to the model bundle is greater than the pruning_thr.\n    This serves to filter the neighborhood area (search space) to get\n    streamlines that are like the model bundle.\n    Recommended range of the pruning_thr is 8 - 12 mm.\n\n    - reduction_distance and pruning_distance : Distance method used internally.\n    Minimum Diferect Flip distance (mdf) or Mean Average Minimum (mam).\n    Default is set to mdf.\n\n    - slr : If slr flag is set to True, local registration of model bundle with\n    neighbouring area will be performed. Default and recommended is True.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read Arcuate Fasciculus Left and Corticospinal Tract Left bundles from already\nfetched atlas data to use them as model bundle. Let's visualize the\nArcuate Fasciculus Left model bundle.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model_af_l_file, model_cst_l_file = get_two_hcp842_bundles()\nsft_af_l = load_trk(model_af_l_file, \"same\", bbox_valid_check=False)\nmodel_af_l = sft_af_l.streamlines\n\ninteractive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(model_af_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='AF_L_model_bundle.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: AF_L_model_bundle.png\n   :align: center\n\n   Model Arcuate Fasciculus Left bundle\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "rb = RecoBundles(moved, verbose=True, rng=np.random.RandomState(2001))\n\nrecognized_af_l, af_l_labels = rb.recognize(model_bundle=model_af_l,\n                                            model_clust_thr=0.1,\n                                            reduction_thr=15,\n                                            pruning_thr=7,\n                                            reduction_distance='mdf',\n                                            pruning_distance='mdf',\n                                            slr=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize extracted Arcuate Fasciculus Left bundle\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(recognized_af_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='AF_L_recognized_bundle.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: AF_L_recognized_bundle.png\n   :align: center\n\n   Extracted Arcuate Fasciculus Left bundle\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Save the bundle as a trk file. Let's save the recognized bundle in the\ncommon space (atlas space), in this case, MNI space.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_af_l = StatefulTractogram(recognized_af_l, atlas_header,\n                               Space.RASMM)\nsave_trk(reco_af_l, \"AF_L_rec_1.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's save the recognized bundle in the original space of the subject anatomy.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_af_l = StatefulTractogram(target[af_l_labels], target_header,\n                               Space.RASMM)\nsave_trk(reco_af_l, \"AF_L_org_1.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, let's increase the reduction_thr and pruning_thr values.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "recognized_af_l, af_l_labels = rb.recognize(model_bundle=model_af_l,\n                                            model_clust_thr=0.1,\n                                            reduction_thr=20,\n                                            pruning_thr=10,\n                                            reduction_distance='mdf',\n                                            pruning_distance='mdf',\n                                            slr=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize extracted Arcuate Fasciculus Left bundle.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(recognized_af_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='AF_L_recognized_bundle2.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: AF_L_recognized_bundle2.png\n   :align: center\n\n   Extracted Arcuate Fasciculus Left bundle\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Save the bundle as a trk file. Let's save the recognized bundle in the\ncommon space (atlas space), in this case, MNI space.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_af_l = StatefulTractogram(recognized_af_l, atlas_header,\n                               Space.RASMM)\nsave_trk(reco_af_l, \"AF_L_rec_2.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's save the recognized bundle in the original space of the subject anatomy.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_af_l = StatefulTractogram(target[af_l_labels], target_header,\n                               Space.RASMM)\nsave_trk(reco_af_l, \"AF_L_org_2.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, let's increase the reduction_thr and pruning_thr values further.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "recognized_af_l, af_l_labels = rb.recognize(model_bundle=model_af_l,\n                                            model_clust_thr=0.1,\n                                            reduction_thr=25,\n                                            pruning_thr=12,\n                                            reduction_distance='mdf',\n                                            pruning_distance='mdf',\n                                            slr=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize extracted Arcuate Fasciculus Left bundle.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(recognized_af_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='AF_L_recognized_bundle3.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: AF_L_recognized_bundle3.png\n   :align: center\n\n   Extracted Arcuate Fasciculus Left bundle\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Save the bundle as a trk file. Let's save the recognized bundle in the\ncommon space (atlas space), in this case, MNI space.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_af_l = StatefulTractogram(recognized_af_l, atlas_header,\n                               Space.RASMM)\nsave_trk(reco_af_l, \"AF_L_rec_3.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's save the recognized bundle in the original space of the subject anatomy.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_af_l = StatefulTractogram(target[af_l_labels], target_header,\n                               Space.RASMM)\nsave_trk(reco_af_l, \"AF_L_org_3.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's apply auto-calibrated RecoBundles on the output of standard RecoBundles.\nThis step will filter out the outlier streamlines. This time, the RecoBundles'\nextracted bundle will serve as a model bundle. As a rule of thumb, provide\nlarger threshold values in standard RecoBundles function and smaller values in\nthe auto-calibrated RecoBundles (refinement) step.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "r_recognized_af_l, r_af_l_labels = rb.refine(model_bundle=model_af_l,\n                                             pruned_streamlines=recognized_af_l,\n                                             model_clust_thr=0.1,\n                                             reduction_thr=15,\n                                             pruning_thr=6,\n                                             reduction_distance='mdf',\n                                             pruning_distance='mdf',\n                                             slr=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize extracted refined Arcuate Fasciculus Left bundle.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(r_recognized_af_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='AF_L_refine_recognized_bundle.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: AF_L_refine_recognized_bundle.png\n   :align: center\n\n   Extracted Arcuate Fasciculus Left bundle\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Save the bundle as a trk file. Let's save the recognized bundle in the\ncommon space (atlas space), in this case, MNI space.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_af_l = StatefulTractogram(r_recognized_af_l, atlas_header,\n                               Space.RASMM)\nsave_trk(reco_af_l, \"AF_L_rec_refine.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's save the recognized bundle in the original space of the subject anatomy.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_af_l = StatefulTractogram(target[r_af_l_labels], target_header,\n                               Space.RASMM)\nsave_trk(reco_af_l, \"AF_L_org_refine.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's load Corticospinal Tract Left model bundle and visualize it.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sft_cst_l = load_trk(model_cst_l_file, \"same\", bbox_valid_check=False)\nmodel_cst_l = sft_cst_l.streamlines\n\ninteractive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(model_cst_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='CST_L_model_bundle.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: CST_L_model_bundle.png\n   :align: center\n\n   The Corticospinal tract model bundle\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "recognized_cst_l, cst_l_labels = rb.recognize(model_bundle=model_cst_l,\n                                              model_clust_thr=0.1,\n                                              reduction_thr=15,\n                                              pruning_thr=7,\n                                              reduction_distance='mdf',\n                                              pruning_distance='mdf',\n                                              slr=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize extracted Corticospinal tract Left bundle\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(recognized_cst_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='CST_L_recognized_bundle.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: CST_L_recognized_bundle.png\n   :align: center\n\n   Extracted Corticospinal tract Left bundle\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Save the bundle as a trk file. Let's save the recognized bundle in the\ncommon space (atlas space), in this case, MNI space.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_cst_l = StatefulTractogram(recognized_cst_l, atlas_header,\n                                Space.RASMM)\nsave_trk(reco_cst_l, \"CST_L_rec_1.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's save the recognized bundle in the original space of the subject anatomy.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_cst_l = StatefulTractogram(target[cst_l_labels], target_header,\n                                Space.RASMM)\nsave_trk(reco_cst_l, \"CST_L_org_1.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, let's increase the reduction_thr and pruning_thr values.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "recognized_cst_l, cst_l_labels = rb.recognize(model_bundle=model_cst_l,\n                                              model_clust_thr=0.1,\n                                              reduction_thr=20,\n                                              pruning_thr=10,\n                                              reduction_distance='mdf',\n                                              pruning_distance='mdf',\n                                              slr=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize extracted Corticospinal tract Left bundle.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(recognized_cst_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='CST_L_recognized_bundle2.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: CST_L_recognized_bundle2.png\n   :align: center\n\n   Extracted Corticospinal tract Left bundle\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Save the bundle as a trk file. Let's save the recognized bundle in the\ncommon space (atlas space), in this case, MNI space.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_cst_l = StatefulTractogram(recognized_cst_l, atlas_header,\n                                Space.RASMM)\nsave_trk(reco_cst_l, \"CST_L_rec_2.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's save the recognized bundle in the original space of the subject anatomy.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_cst_l = StatefulTractogram(target[cst_l_labels], target_header,\n                                Space.RASMM)\nsave_trk(reco_cst_l, \"CST_L_org_2.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, let's increase the reduction_thr and pruning_thr values further.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "recognized_cst_l, cst_l_labels = rb.recognize(model_bundle=model_cst_l,\n                                              model_clust_thr=0.1,\n                                              reduction_thr=25,\n                                              pruning_thr=12,\n                                              reduction_distance='mdf',\n                                              pruning_distance='mdf',\n                                              slr=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize extracted Corticospinal tract Left bundle.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(recognized_cst_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='CST_L_recognized_bundle3.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: CST_L_recognized_bundle3.png\n   :align: center\n\n   Extracted Corticospinal tract Left bundle\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Save the bundle as a trk file. Let's save the recognized bundle in the\ncommon space (atlas space), in this case, MNI space.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_cst_l = StatefulTractogram(recognized_cst_l, atlas_header,\n                                Space.RASMM)\nsave_trk(reco_cst_l, \"CST_L_rec_3.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's save the recognized bundle in the original space of the subject anatomy.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_cst_l = StatefulTractogram(target[cst_l_labels], target_header,\n                                Space.RASMM)\nsave_trk(reco_cst_l, \"CST_L_org_3.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's apply auto-calibrated RecoBundles on the output of standard RecoBundles.\nThis step will filter out the outlier streamlines. This time, the RecoBundles'\nextracted bundle will serve as a model bundle. As a rule of thumb, provide\nlarger threshold values in standard RecoBundles function and smaller values in\nthe auto-calibrated RecoBundles (refinement) step.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "r_recognized_cst_l, r_cst_l_labels = rb.refine(\n                                     model_bundle=model_cst_l,\n                                     pruned_streamlines=recognized_cst_l,\n                                     model_clust_thr=0.1,\n                                     reduction_thr=15,\n                                     pruning_thr=6,\n                                     reduction_distance='mdf',\n                                     pruning_distance='mdf',\n                                     slr=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "let's visualize extracted refined Corticospinal tract Left bundle.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "interactive = False\n\nscene = window.Scene()\nscene.SetBackground(1, 1, 1)\nscene.add(actor.line(r_recognized_cst_l))\nscene.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),\n                 position=(-360.11, -30.46, -40.44),\n                 view_up=(-0.03, 0.028, 0.89))\nwindow.record(scene, out_path='CST_L_refine_recognized_bundle.png',\n              size=(600, 600))\nif interactive:\n    window.show(scene)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: CST_L_refine_recognized_bundle.png\n   :align: center\n\n   Extracted refined Corticospinal tract Left bundle\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Save the bundle as a trk file. Let's save the recognized bundle in the\ncommon space (atlas space), in this case, MNI space.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_cst_l = StatefulTractogram(r_recognized_cst_l, atlas_header,\n                                Space.RASMM)\nsave_trk(reco_cst_l, \"CST_L_rec_refine.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's save the recognized bundle in the original space of the subject anatomy.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reco_cst_l = StatefulTractogram(target[r_cst_l_labels], target_header,\n                                Space.RASMM)\nsave_trk(reco_cst_l, \"CST_L_org_refine.trk\", bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This example shows how changing different threshold parameters can change the\noutput. It is up to the user to decide what is the desired output and select\nparameters accordingly.\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n\n.. [Garyfallidis17] Garyfallidis et al. Recognition of white matter\n        bundles using local and global streamline-based registration\n        and clustering, Neuroimage, 2017.\n\n.. [Chandio2020] Chandio, B.Q., Risacher, S.L., Pestilli, F.,\n        Bullock, D., Yeh, FC., Koudoro, S., Rokem, A., Harezlak, J., and\n        Garyfallidis, E. Bundle analytics, a computational framework for\n        investigating the shapes and profiles of brain pathways across\n        populations. Sci Rep 10, 17149 (2020)\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
}