{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Visualizing AFQ derivatives\n\nVisualizing the results of a pyAFQ analysis is useful because it allows us to\ninspect the results of the analysis and to communicate the results to others.\nThe pyAFQ pipeline produces a number of different kinds of outputs, including\nvisualizations that can be used for quality control and for quick examination\nof the results of the analysis.\n\nHowever, when communicating the results of pyAFQ analysis, it is often useful\nto have more specific control over the visualization that is produced. In\naddition, it is often useful to have visualizations that are visually appealing\nand striking. In this tutorial, we will use the [fury](https://fury.gl/)\nlibrary [1]_ to visualize outputs of pyAFQ as publication-ready figures.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import os\nimport os.path as op\nimport nibabel as nib\nimport numpy as np\nfrom math import radians\n\nfrom dipy.io.streamline import load_trk\nfrom dipy.tracking.streamline import transform_streamlines\n\nfrom fury import actor, window\nfrom fury.colormap import create_colormap\n\n\nimport AFQ.data.fetch as afd\nfrom AFQ.viz.utils import PanelFigure"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Get some data from HBN POD2\nThe Healthy Brain Network Preprocessed Open Diffusion Derivatives (HBN POD2)\nis a collection of resources based on the Healthy Brain Network dataset\n[2, 3]_. HBN POD2 includes data derivatives - including pyAFQ derivatives -\nfrom more than 2,000 subjects. The data and the derivatives can be browsed at\nhttps://fcp-indi.s3.amazonaws.com/index.html#data/Projects/HBN/BIDS_curated/\n\nHere, we will visualize the results from one subject, together with their\nanatomy and using several variations. We start by downloading their\npyAFQ-processed data using fetcher functions that download both the\npreprocessed data, as well as the pyAFQ-processed data (Note that this\nwill take up about 1.75 GB of disk space):\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "afd.fetch_hbn_preproc([\"NDARAA948VFH\"])\nstudy_path = afd.fetch_hbn_afq([\"NDARAA948VFH\"])[1]\n\nderiv_path = op.join(\n    study_path, \"derivatives\")\n\nafq_path = op.join(\n    deriv_path,\n    'afq',\n    'sub-NDARAA948VFH',\n    'ses-HBNsiteRU')\n\nbundle_path = op.join(afq_path,\n                      'bundles')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Read data into memory\nThe bundle coordinates from pyAFQ are always saved in the reference frame of\nthe diffusion data from which they are generated, so we need an image file\nwith the dMRI coordinates as a reference for loading the data (we could also\nuse `\"same\"` here).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fa_img = nib.load(op.join(afq_path,\n                          'sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_model-DKI_FA.nii.gz'))\nfa = fa_img.get_fdata()\nsft_arc = load_trk(op.join(bundle_path,\n                           'sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_space-RASMM_model-CSD_desc-prob-afq-ARC_L_tractography.trk'), fa_img)\nsft_cst = load_trk(op.join(bundle_path,\n                           'sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_space-RASMM_model-CSD_desc-prob-afq-CST_L_tractography.trk'), fa_img)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Transform into the T1w reference frame\nOur first goal is to visualize the bundles with a background of the\nT1-weighted image, which provides anatomical context. We read in this data and\ntransform the bundle coordinates, first into the RASMM common coordinate frame\nand then subsequently into the coordinate frame of the T1-weighted data (if\nyou find this confusing, you can brush up on this topic in the\n[nibabel documentation](https://nipy.org/nibabel/coordinate_systems.html)).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t1w_img = nib.load(op.join(deriv_path,\n                           'qsiprep/sub-NDARAA948VFH/anat/sub-NDARAA948VFH_desc-preproc_T1w.nii.gz'))\nt1w = t1w_img.get_fdata()\nsft_arc.to_rasmm()\nsft_cst.to_rasmm()\narc_t1w = transform_streamlines(sft_arc.streamlines,\n                                np.linalg.inv(t1w_img.affine))\ncst_t1w = transform_streamlines(sft_cst.streamlines,\n                                np.linalg.inv(t1w_img.affine))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualizing bundles with principal direction coloring\nThe first visualization we will create will have the streamlines colored\naccording to their direction. The color of each streamline will be RGB encoded\naccording to the RAS of the average orientation of its segments. This is the\ndefault behavior in fury.\n\nFury uses \"actors\" to render different kind of graphics. For the bundles, we\nwill use the `line` actor. These objects are wrappers around the [vtkActor\nclass](https://vtk.org/doc/nightly/html/classvtkActor.html), so methods of\nthat class (like `GetProperty()`) are available to use. We like to set the\naesthetics of the streamlines, so that they are rendered as tubes and with\nslightly thicker line-width than the default. We create a function that sets\nthese properties of the line actor via the `GetProperty` method. We will\nreuse this function later on, also setting the key-word arguments to the call\nto `actor.line`, but for now we use the default setting, which colors each\nstreamline based on the RAS orientation, and we set the line width to 8.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "arc_actor = actor.streamlines(arc_t1w, thickness=8)\ncst_actor = actor.streamlines(cst_t1w, thickness=8)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Slicer actors\nThe anatomical image is rendered using `slicer` actors. These are actors that\nvisualize one slice of a three dimensional volume. We call the function on the\nT1-weighted data, selecting the x slice that is half-way through the x\ndimension of the image (`shape[0]`) and the z slice that is a third of a\nway through that x dimension of the image (`shape[-1]`).\nWe set the visibility of the y slice to `False`\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "slicer = actor.data_slicer(t1w,\n                           visibility=(\n                               True, False, True),\n                           initial_slices=(\n                               t1w.shape[0] // 2,\n                               t1w.shape[1] // 2,\n                               t1w.shape[-1] // 3))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Making a `scene`\nThe next kind of fury object we will be working with is a `window.Scene`\nobject. This is the (3D!) canvas on which we are drawing the actors. We\ninitialize this object and call the `scene.add` method to add the actors.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "scene = window.Scene()\n\nscene.add(arc_actor)\nscene.add(cst_actor)\nscene.add(slicer)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Showing the visualization\nIf you are working in an interactive session, you can call::\n\n    window.show(scene, size=(1200, 1200))\n\nto see what the visualization looks like. This would pop up a window that will\nshow you the visualization as it is now. You can interact with this\nvisualization using your mouse to rotate the image in 3D, and mouse+ctrl or\nmouse+shift to pan and rotate it in plane, respectively. Use the scroll up and\nscroll down in your mouse to zoom in and out. Once you have found a view of\nthe data that you like, you can close the window (as long as its open, it is\nblocking execution of any further commands in the Python interpreter!)\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Record the visualization\nIf you are not working in an interactive session, or you have already figured\nout the camera settings that work for you, you can go ahead and record the\nimage into a png file. We use a pretty high resolution here (2400 by 2400) so\nthat we get a nice crisp image. That also means the file is pretty large.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "out_folder = op.join(afd.afq_home, \"VizExample\")\nos.makedirs(out_folder, exist_ok=True)\n\ndef save_png(scene, name):\n    \"\"\"Helper function to PNGs in this example.\"\"\"\n    show_m = window.ShowManager(\n        scene=scene, window_type=\"offscreen\",\n        size=(2400, 2400)\n    )\n    window.update_camera(show_m.screens[0].camera, None, scene)\n    show_m.screens[0].controller.rotate((0, radians(-90)), None)\n    show_m.render()\n    show_m.window.draw()\n    show_m.snapshot(op.join(out_folder, name))\n\nsave_png(scene, 'arc_cst1.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Setting bundle colors\nThe default behavior produces a nice image! But we often want to be able to\ndifferentiate streamlines based on the bundle to which they belong. For this,\nwe will set the color of each streamline, based on the bundle. We often use\nthe Tableau20 color palette to set the colors for the different bundles. The\n`actor.line` initializer takes `colors` as a keyword argument and one of the\noptions to pass here is an RGB triplet, which will determine the color of all\nof the streamlines in that actor. We can get the Tableau20 RGB triplets from\nMatplotlib's colormap module. We initialize our bundle actors again and clear\nthe scene before adding these new actors and adding back the slice actors. We\nthen call `record` to create this new figure.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from matplotlib.cm import tab20\ncolor_arc = tab20.colors[18]\ncolor_cst = tab20.colors[2]\n\narc_actor = actor.streamlines(arc_t1w, thickness=8, colors=color_arc)\ncst_actor = actor.streamlines(cst_t1w, thickness=8, colors=color_cst)\n\nscene.clear()\n\nscene.add(arc_actor)\nscene.add(cst_actor)\nscene.add(slicer)\n\nsave_png(scene, 'arc_cst2.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding core bundles with tract profiles\nNext, we'd like to add a representation of the core of each bundle and plot\nthe profile of some property along the length of the bundle. Another option\nfor setting streamline colors is for each point along a streamline to have a a\ndifferent color. Here, we will do this with only one streamline, which\nrepresents the core of the fiber -- it's median trajectory, and with a profile\nof the FA. This requires a few distinct steps. The first is adding another\nactor for each of the core fibers. We determine the core bundle as the median\nof the coordinates of the streamlines after each streamline is resampled to\n100 equi-distant points. In the next step, we extract the FA profiles for each\nof the bundles, using the `afq_profile` function. Before we do this, we make\nsure that the streamlines are back in the voxel coordinate frame of the FA\nvolume. The core bundle and FA profile are put together into a new line actor,\nwhich we initialize as a thick tube (line width of 40) and with a call to\ncreate a colormap from the FA profiles (we can choose which color-map to use\nhere. In this case, we chose `'viridis'`). Note that because we are passing\nonly one streamline into the line actor, we have to pass it inside of square\nbrackets (`[]`). This is because the `actor.line` initializer expects a\nsequence of streamlines as input. Before plotting things again, we initialize\nour bundle line actors again. This time, we set the opacity of the lines to\n0.1, so that they do not occlude the view of the core bundle visualizations.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>In this case, the profile is FA, but you can use any sequence of 100\n  numbers in place of the FA profiles: group differences, p-values, etc.</p></div>\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.tracking.streamline import set_number_of_points\ncore_arc = np.median(np.asarray(set_number_of_points(arc_t1w, 100)), axis=0)\ncore_cst = np.median(np.asarray(set_number_of_points(cst_t1w, 100)), axis=0)\n\n\nfrom dipy.stats.analysis import afq_profile\nsft_arc.to_vox()\nsft_cst.to_vox()\narc_profile = afq_profile(fa, sft_arc.streamlines, affine=np.eye(4))\ncst_profile = afq_profile(fa, sft_cst.streamlines, affine=np.eye(4))\n\ncore_arc_actor = actor.streamlines(\n    [core_arc],\n    thickness=40,\n    colors=create_colormap(arc_profile, name='viridis')\n)\n\ncore_cst_actor = actor.streamlines(\n    [core_cst],\n    thickness=40,\n    colors=create_colormap(cst_profile, name='viridis')\n)\n\nscene.clear()\n\narc_actor = actor.streamlines(arc_t1w, thickness=8, colors=color_arc, opacity=0.1)\ncst_actor = actor.streamlines(cst_t1w, thickness=8, colors=color_cst, opacity=0.1)\n\nscene.add(arc_actor)\nscene.add(cst_actor)\nscene.add(slicer)\nscene.add(core_arc_actor)\nscene.add(core_cst_actor)\n\nsave_png(scene, 'arc_cst3.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Adding ROIs\nAnother element that we can add into the visualization are volume renderings\nof regions of interest. For example, the waypoint ROIs that were used to\ndefine the bundles. Here, we will add the waypoints that were used to define\nthe arcuate fasciculus in this subject. We start by reading this data in.\nBecause it is represented in the space of the diffusion data, it too needs to\nbe resampled into the T1-weighted space, before being visualized. After\nresampling, we might have some values between 0 and 1 because of the\ninterpolation from the low resolution of the diffusion into the high\nresolution of the T1-weighted. We will include in the volume rendering any\nvalues larger than 0. The main change from the previous visualizations is the\naddition of a `contour_from_roi` actor for each of the ROIs. We select another\ncolor from the Tableau 20 palette to represent this, and use an opacity of\n0.5.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>In this case, the surface rendered is of the waypoint ROIs, but very\n  similar code can be used to render other surfaces, provided a volume that\n  contains that surface. For example, the gray-white matter boundary could\n  be visualized provided an array with a binary representation of the volume #   enclosed in this boundary</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.align import resample\n\nwaypoint1 = nib.load(\n    op.join(\n        afq_path,\n        \"ROIs\", \"sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_desc-ROI-ARC_L-1-include.nii.gz\"))\n\nwaypoint2 = nib.load(\n    op.join(\n        afq_path,\n        \"ROIs\", \"sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_desc-ROI-ARC_L-2-include.nii.gz\"))\n\nwaypoint1_xform = resample(waypoint1, t1w_img)\nwaypoint2_xform = resample(waypoint2, t1w_img)\nwaypoint1_data = waypoint1_xform.get_fdata() > 0\nwaypoint2_data = waypoint2_xform.get_fdata() > 0\n\nscene.clear()\n\narc_actor = actor.streamlines(arc_t1w, thickness=8, colors=color_arc)\ncst_actor = actor.streamlines(cst_t1w, thickness=8, colors=color_cst)\n\nscene.add(arc_actor)\nscene.add(cst_actor)\nscene.add(slicer)\n\nsurface_color = tab20.colors[0]\n\nwaypoint1_actor = actor.contour_from_roi(waypoint1_data,\n                                         color=surface_color,\n                                         opacity=0.5)\n\nwaypoint2_actor = actor.contour_from_roi(waypoint2_data,\n                                         color=surface_color,\n                                         opacity=0.5)\n\n\nscene.add(waypoint1_actor)\nscene.add(waypoint2_actor)\n\nsave_png(scene, 'arc_cst4.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualizing tracts and tract profiles with a \"glass brain\"\nDisplaying tracts together with slices of anatomical data is beautiful\nbut sometimes poses a challenge because of occlusion. Another option is\nto visualize the data with a \"glass brain\". That is, a faint contour\nshowing where the edges of the brain are, so that the tracts are visible\nand the anatomical context can be understood. One way to generate the\ncontour of the brain is to use a brain mask generated from anatomical data.\nWe set its display color to black (`[0, 0, 0]`) and its opacity to a very low\nvalue so that the tracts can easily be seen. This visualization is best with\na bright background, so we set the background to white (`[1, 1, 1]`).\nIn addition, we will demonstrate another way to show tract profile\ninformation (based on Luo et al. 2025 [4]_): here the tract profile is\ninterpolated onto each streamline in the tract and mapped to a colormap.\nThis vividly shows the variations in tract profile values over space.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "brain_mask_img = nib.load(op.join(deriv_path,\n    'qsiprep/sub-NDARAA948VFH/anat/sub-NDARAA948VFH_desc-brain_mask.nii.gz'))\n\nbrain_mask_data = brain_mask_img.get_fdata()\n\nscene.clear()\n\nbrain_actor = actor.contour_from_roi(brain_mask_data,\n                                     color=[0, 0, 0],\n                                     opacity=0.05)\n\nscene.add(brain_actor)\n\nfor sl in arc_t1w:\n    # interpolate the 100 tract profiles values to match the number of points\n    # in the streamline:\n    interpolated_values = np.interp(np.linspace(0, 1, len(sl)),\n                                    np.linspace(0, 1, len(arc_profile)),\n                                    arc_profile)\n    colors = create_colormap(interpolated_values, name='Spectral')\n    line_actor = actor.streamlines([sl], thickness=8, colors=colors)\n    scene.add(line_actor)\n\nscene.background = (1, 1, 1)\nsave_png(scene, 'arc_cst5.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Making a Figure out of many fury panels\nWe can also make a figure that contains multiple panels, each of which\ncontains a different visualization. This is useful for communicating the\nresults of an analysis. Here, we will make a figure with four panels, using\nsome of the visualizations we have already created. We will use some\nconvenient methods from pyAFQ.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pf = PanelFigure(3, 2, 6, 9)\npf.add_img(op.join(out_folder, 'arc_cst1.png'), 0, 0)\npf.add_img(op.join(out_folder, 'arc_cst2.png'), 1, 0)\npf.add_img(op.join(out_folder, 'arc_cst3.png'), 0, 1)\npf.add_img(op.join(out_folder, 'arc_cst5.png'), 1, 1)\npf.format_and_save_figure(f\"arc_cst_fig.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n.. [1] Garyfallidis et al., (2021). FURY: advanced scientific visualization.\n   Journal of Open Source Software, 6(64), 3384,\n   https://doi.org/10.21105/joss.03384\n\n.. [2] Alexander LM, Escalera J, Ai L, et al. An open resource for\n    transdiagnostic research in pediatric mental health and learning\n    disorders. Sci Data. 2017;4:170181.\n\n.. [3] Richie-Halford A, Cieslak M, Ai L, et al. An analysis-ready and\n    quality controlled resource for pediatric brain white-matter research.\n    Scientific Data. 2022;9(1):1-27.\n\n.. [4] Luo A, et al. Two Axes of White Matter Development. In prep.\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.13.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}