{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Read/Write streamline files\n\n## Overview\n\nDIPY_ can read and write many different file formats. In this example\nwe give a short introduction on how to use it for loading or saving\nstreamlines. The new stateful tractogram class was made to reduce errors\ncaused by spatial transformation and complex file format convention.\n\nRead `faq`\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\n\nimport nibabel as nib\nimport numpy as np\nfrom dipy.io.stateful_tractogram import Space, StatefulTractogram\nfrom dipy.io.streamline import load_tractogram, save_tractogram\nfrom dipy.io.utils import (create_nifti_header, get_reference_info,\n                           is_header_compatible)\nfrom dipy.tracking.streamline import select_random_set_of_streamlines\nfrom dipy.tracking.utils import density_map\n\nfrom dipy.data.fetcher import (fetch_file_formats,\n                               get_file_formats)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First fetch the dataset that contains 5 tractography file of 5 file formats:\n\n- cc_m_sub.trk\n- laf_m_sub.tck\n- lpt_m_sub.fib\n- raf_m_sub.vtk\n- rpt_m_sub.dpy\n\nAnd their reference anatomy, common to all 5 files:\n\n- template0.nii.gz\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fetch_file_formats()\nbundles_filename, ref_anat_filename = get_file_formats()\nfor filename in bundles_filename:\n    print(os.path.basename(filename))\nreference_anatomy = nib.load(ref_anat_filename)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load tractogram will support 5 file formats, functions like load_trk or\nload_tck will simply be restricted to one file format\n\nTRK files contain their own header (when written properly), so they\ntechnically do not need a reference. (See how below)\n\n``cc_trk = load_tractogram(bundles_filename[0], 'same')``\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cc_sft = load_tractogram(bundles_filename[0], reference_anatomy)\nprint(cc_sft)\nlaf_sft = load_tractogram(bundles_filename[1], reference_anatomy)\nraf_sft = load_tractogram(bundles_filename[3], reference_anatomy)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "These files contain invalid streamlines (negative values once in voxel space)\nThis is not considered a valid tractography file, but it is possible to load\nit anyway.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "lpt_sft = load_tractogram(bundles_filename[2], reference_anatomy,\n                          bbox_valid_check=False)\nrpt_sft = load_tractogram(bundles_filename[4], reference_anatomy,\n                          bbox_valid_check=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The function ``load_tractogram`` requires a reference, any of the following\ninputs is considered valid (as long as they are in the same share space)\n- Nifti filename\n- Trk filename\n- nib.nifti1.Nifti1Image\n- nib.streamlines.trk.TrkFile\n- nib.nifti1.Nifti1Header\n- Trk header (dict)\n- Stateful Tractogram\n\nThe reason why this parameter is required is to guarantee all information\nrelated to space attributes is always present.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "affine, dimensions, voxel_sizes, voxel_order = get_reference_info(\n    reference_anatomy)\nprint(affine)\nprint(dimensions)\nprint(voxel_sizes)\nprint(voxel_order)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If you have a Trk file that was generated using a particular anatomy,\nto be considered valid all fields must correspond between the headers.\nIt can be easily verified using this function, which also accept\nthe same variety of input as ``get_reference_info``\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(is_header_compatible(reference_anatomy, bundles_filename[0]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If a TRK was generated with a valid header, but the reference NIFTI was lost\na header can be generated to then generate a fake NIFTI file.\n\nIf you wish to manually save Trk and Tck file using nibabel streamlines\nAPI for more freedom of action (not recommended for beginners) you can\ncreate a valid header using create_tractogram_header\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nifti_header = create_nifti_header(affine, dimensions, voxel_sizes)\nnib.save(nib.Nifti1Image(np.zeros(dimensions), affine, nifti_header),\n         'fake.nii.gz')\nnib.save(reference_anatomy, os.path.basename(ref_anat_filename))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Once loaded, no matter the original file format, the stateful tractogram is\nself-contained and maintains a valid state. By requiring a reference the\ntractogram's spatial transformation can be easily manipulated.\n\nLet's save all files as TRK to visualize in TrackVis for example.\nHowever, when loaded the lpt and rpt files contain invalid streamlines and\nfor particular operations/tools/functions it is safer to remove them\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "save_tractogram(cc_sft, 'cc.trk')\nsave_tractogram(laf_sft, 'laf.trk')\nsave_tractogram(raf_sft, 'raf.trk')\n\nprint(lpt_sft.is_bbox_in_vox_valid())\nlpt_sft.remove_invalid_streamlines()\nprint(lpt_sft.is_bbox_in_vox_valid())\nsave_tractogram(lpt_sft, 'lpt.trk')\n\nprint(rpt_sft.is_bbox_in_vox_valid())\nrpt_sft.remove_invalid_streamlines()\nprint(rpt_sft.is_bbox_in_vox_valid())\nsave_tractogram(rpt_sft, 'rpt.trk')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Some functions in DIPY require streamlines to be in voxel space so computation\ncan be performed on a grid (connectivity matrix, ROIs masking, density map).\nThe stateful tractogram class provides safe functions for such manipulation.\nThese functions can be called safely over and over, by knowing in which state\nthe tractogram is operating, and compute only necessary transformations\n\nNo matter the state, functions such as ``save_tractogram`` or\n``removing_invalid_coordinates`` can be called safely and the transformations\nare handled internally when needed.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cc_sft.to_voxmm()\nprint(cc_sft.space)\ncc_sft.to_rasmm()\nprint(cc_sft.space)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's move them all to voxel space, subsample them to 100 streamlines,\ncompute a density map and save everything for visualisation in another\nsoftware such as Trackvis or MI-Brain.\n\nTo access volume information in a grid, the corner of the voxel must be\nconsidered the origin in order to prevent negative values.\nAny operation doing interpolation or accessing a grid must use the\nfunction 'to_vox()' and 'to_corner()'\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cc_sft.to_vox()\nlaf_sft.to_vox()\nraf_sft.to_vox()\nlpt_sft.to_vox()\nrpt_sft.to_vox()\n\ncc_sft.to_corner()\nlaf_sft.to_corner()\nraf_sft.to_corner()\nlpt_sft.to_corner()\nrpt_sft.to_corner()\n\ncc_streamlines_vox = select_random_set_of_streamlines(cc_sft.streamlines,\n                                                      1000)\nlaf_streamlines_vox = select_random_set_of_streamlines(laf_sft.streamlines,\n                                                       1000)\nraf_streamlines_vox = select_random_set_of_streamlines(raf_sft.streamlines,\n                                                       1000)\nlpt_streamlines_vox = select_random_set_of_streamlines(lpt_sft.streamlines,\n                                                       1000)\nrpt_streamlines_vox = select_random_set_of_streamlines(rpt_sft.streamlines,\n                                                       1000)\n\n# Same dimensions for every stateful tractogram, can be re-use\naffine, dimensions, voxel_sizes, voxel_order = cc_sft.space_attributes\ncc_density = density_map(cc_streamlines_vox, np.eye(4), dimensions)\nlaf_density = density_map(laf_streamlines_vox, np.eye(4), dimensions)\nraf_density = density_map(raf_streamlines_vox, np.eye(4), dimensions)\nlpt_density = density_map(lpt_streamlines_vox, np.eye(4), dimensions)\nrpt_density = density_map(rpt_streamlines_vox, np.eye(4), dimensions)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Replacing streamlines is possible, but if the state was modified between\noperations such as this one is not recommended:\n-> cc_sft.streamlines = cc_streamlines_vox\n\nIt is recommended to re-create a new StatefulTractogram object and\nexplicitly specify in which space the streamlines are. Be careful to follow\nthe order of operations.\n\nIf the tractogram was from a Trk file with metadata, this will be lost.\nIf you wish to keep metadata while manipulating the number or the order\nlook at the function StatefulTractogram.remove_invalid_streamlines() for more\ndetails\n\nIt is important to mention that once the object is created in a consistent state\nthe ``save_tractogram`` function will save a valid file. And then the function\n``load_tractogram`` will load them in a valid state.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cc_sft = StatefulTractogram(cc_streamlines_vox, reference_anatomy, Space.VOX)\nlaf_sft = StatefulTractogram(laf_streamlines_vox, reference_anatomy, Space.VOX)\nraf_sft = StatefulTractogram(raf_streamlines_vox, reference_anatomy, Space.VOX)\nlpt_sft = StatefulTractogram(lpt_streamlines_vox, reference_anatomy, Space.VOX)\nrpt_sft = StatefulTractogram(rpt_streamlines_vox, reference_anatomy, Space.VOX)\n\nprint(len(cc_sft), len(laf_sft), len(raf_sft), len(lpt_sft), len(rpt_sft))\nsave_tractogram(cc_sft, 'cc_1000.trk')\nsave_tractogram(laf_sft, 'laf_1000.trk')\nsave_tractogram(raf_sft, 'raf_1000.trk')\nsave_tractogram(lpt_sft, 'lpt_1000.trk')\nsave_tractogram(rpt_sft, 'rpt_1000.trk')\n\nnib.save(nib.Nifti1Image(cc_density, affine, nifti_header),\n         'cc_density.nii.gz')\nnib.save(nib.Nifti1Image(laf_density, affine, nifti_header),\n         'laf_density.nii.gz')\nnib.save(nib.Nifti1Image(raf_density, affine, nifti_header),\n         'raf_density.nii.gz')\nnib.save(nib.Nifti1Image(lpt_density, affine, nifti_header),\n         'lpt_density.nii.gz')\nnib.save(nib.Nifti1Image(rpt_density, affine, nifti_header),\n         'rpt_density.nii.gz')"
      ]
    }
  ],
  "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
}