{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# Visualizing AFQ derivatives\n", "\n", "Visualizing the results of a pyAFQ analysis is useful because it allows us to\n", "inspect the results of the analysis and to communicate the results to others.\n", "The pyAFQ pipeline produces a number of different kinds of outputs, including\n", "visualizations that can be used for quality control and for quick examination\n", "of the results of the analysis.\n", "\n", "However, when communicating the results of pyAFQ analysis, it is often useful\n", "to have more specific control over the visualization that is produced. In\n", "addition, it is often useful to have visualizations that are visually appealing\n", "and striking. In this tutorial, we will use the [FURY](https://fury.gl/)\n", "library [@Garyfallidis2021] to visualize outputs of pyAFQ as publication-ready figures.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "import os\n", "import os.path as op\n", "import nibabel as nib\n", "import numpy as np\n", "from paths import afq_home\n", "\n", "from dipy.io.streamline import load_trk\n", "from dipy.tracking.streamline import transform_streamlines\n", "\n", "from fury import actor, window\n", "from fury.colormap import create_colormap\n", "\n", "import AFQ.data.fetch as afd\n", "from AFQ.viz.utils import PanelFigure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Get some data from HBN POD2\n", "\n", "The Healthy Brain Network Preprocessed Open Diffusion Derivatives (HBN POD2)\n", "is a collection of resources based on the Healthy Brain Network dataset [](https://doi.org/10.1038/s41597-022-01695-7). \n", "HBN POD2 includes data derivatives - including pyAFQ derivatives -\n", "from more than 2,000 subjects. \n", "\n", ":::{seealso}\n", "The data and the derivatives can be browsed at https://fcp-indi.s3.amazonaws.com/index.html#data/Projects/HBN/BIDS_curated/\n", ":::\n", "\n", "Here, we will visualize the results from one subject, together with their\n", "anatomy and using several variations. We start by downloading their\n", "pyAFQ-processed data using fetcher functions that download both the\n", "preprocessed data, as well as the pyAFQ-processed data [](https://doi.org/10.1038/sdata.2017.181).\n", "\n", ":::{warning}\n", "This process will take up about 1.75 GB of disk space.\n", ":::\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "afd.fetch_hbn_preproc([\"NDARAA948VFH\"])\n", "study_path = afd.fetch_hbn_afq([\"NDARAA948VFH\"])[1]\n", "\n", "deriv_path = op.join(\n", " study_path, \"derivatives\")\n", "\n", "afq_path = op.join(\n", " deriv_path,\n", " 'afq',\n", " 'sub-NDARAA948VFH',\n", " 'ses-HBNsiteRU')\n", "\n", "bundle_path = op.join(afq_path, 'bundles')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read data into memory\n", "\n", "The bundle coordinates from pyAFQ are always saved in the reference frame of\n", "the diffusion data from which they are generated, so we need an image file\n", "with the dMRI coordinates as a reference for loading the data (we could also\n", "use `\"same\"` here).\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": 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'))\n", "fa = fa_img.get_fdata()\n", "sft_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)\n", "sft_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\n", "\n", "Our first goal is to visualize the bundles with a background of the\n", "T1-weighted image, which provides anatomical context. We read in this data and\n", "transform the bundle coordinates, first into the RASMM common coordinate frame\n", "and then subsequently into the coordinate frame of the T1-weighted data (if\n", "you 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, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "t1w_img = nib.load(op.join(deriv_path,\n", " 'qsiprep/sub-NDARAA948VFH/anat/sub-NDARAA948VFH_desc-preproc_T1w.nii.gz'))\n", "t1w = t1w_img.get_fdata()\n", "sft_arc.to_rasmm()\n", "sft_cst.to_rasmm()\n", "arc_t1w = transform_streamlines(sft_arc.streamlines,\n", " np.linalg.inv(t1w_img.affine))\n", "cst_t1w = transform_streamlines(sft_cst.streamlines,\n", " np.linalg.inv(t1w_img.affine))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":::{note}\n", "A virtual frame buffer is needed if you are running this example on a machine that is not connected to a display (\"headless\"). If this is the case, you can either set an environment variable called `XVFB` to `1` or you can run the following code to initialize the virtual frame buffer:\n", "\n", "```python\n", "import xvfbwrapper\n", "from xvfbwrapper import Xvfb\n", "\n", "vdisplay = Xvfb()\n", "vdisplay.start()\n", "```\n", "\n", "And remember to stop it when you are done:\n", "\n", "```python\n", "vdisplay.stop()\n", "```\n", "\n", ":::" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualizing bundles with principal direction coloring\n", "\n", "The first visualization we will create will have the streamlines colored\n", "according to their direction. The color of each streamline will be RGB encoded\n", "according to the RAS of the average orientation of its segments. This is the\n", "default behavior in fury.\n", "\n", "Fury uses \"actors\" to render different kind of graphics. For the bundles, we\n", "will use the `line` actor. These objects are wrappers around the [vtkActor\n", "class](https://vtk.org/doc/nightly/html/classvtkActor.html), so methods of\n", "that class (like `GetProperty()`) are available to use. We like to set the\n", "aesthetics of the streamlines, so that they are rendered as tubes and with\n", "slightly thicker line-width than the default. We create a function that sets\n", "these properties of the line actor via the `GetProperty` method. We will\n", "reuse this function later on, also setting the key-word arguments to the call\n", "to `actor.line`, but for now we use the default setting, which colors each\n", "streamline based on the RAS orientation, and we set the line width to 8.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "def lines_as_tubes(sl, line_width, **kwargs):\n", " line_actor = actor.line(sl, **kwargs)\n", " line_actor.GetProperty().SetRenderLinesAsTubes(1)\n", " line_actor.GetProperty().SetLineWidth(line_width)\n", " return line_actor\n", "\n", "\n", "arc_actor = lines_as_tubes(arc_t1w, 8)\n", "cst_actor = lines_as_tubes(cst_t1w, 8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Slicer actors\n", "\n", "The anatomical image is rendered using `slicer` actors. These are actors that\n", "visualize one slice of a three dimensional volume. Again, we create a helper\n", "function that will slice a volume along the x, y, and z dimensions. This\n", "function returns a list of the slicers we want to include in our\n", "visualization. This can be one, two, or three slicers, depending on how many\n", "of {x,y,z} are set. If you are curious to understand what is going on in this\n", "function, take a look at the documentation for the\n", ":met:`actor.slicer.display_extent` method. \n", "\n", ":::{hint}\n", "For every dimension you select on, you want the full extent of the image on the two *other* two dimensions. We call the function on the T1-weighted data, selecting the `#` x-slice that is half-way through the x-dimension of the image (`shape[0]`) and the z slice that is a third of a way through that x dimension of the image (`shape[-1]`).\n", ":::\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "def slice_volume(data, x=None, y=None, z=None):\n", " slicer_actors = []\n", " slicer_actor_z = actor.slicer(data)\n", " if z is not None:\n", " slicer_actor_z.display_extent(\n", " 0, data.shape[0] - 1,\n", " 0, data.shape[1] - 1,\n", " z, z)\n", " slicer_actors.append(slicer_actor_z)\n", " if y is not None:\n", " slicer_actor_y = slicer_actor_z.copy()\n", " slicer_actor_y.display_extent(\n", " 0, data.shape[0] - 1,\n", " y, y,\n", " 0, data.shape[2] - 1)\n", " slicer_actors.append(slicer_actor_y)\n", " if x is not None:\n", " slicer_actor_x = slicer_actor_z.copy()\n", " slicer_actor_x.display_extent(\n", " x, x,\n", " 0, data.shape[1] - 1,\n", " 0, data.shape[2] - 1)\n", " slicer_actors.append(slicer_actor_x)\n", "\n", " return slicer_actors\n", "\n", "\n", "slicers = slice_volume(t1w, x=t1w.shape[0] // 2, z=t1w.shape[-1] // 3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Making a `scene`\n", "\n", "The next kind of fury object we will be working with is a `window.Scene`\n", "object. This is the (3D!) canvas on which we are drawing the actors. We\n", "initialize this object and call the `scene.add` method to add the actors.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "scene = window.Scene()\n", "\n", "scene.add(arc_actor)\n", "scene.add(cst_actor)\n", "for slicer in slicers:\n", " scene.add(slicer)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Showing the visualization\n", "\n", "If you are working in an interactive session, you can call:\n", "\n", "```python\n", "window.show(scene, size=(1200, 1200), reset_camera=False)\n", "```\n", "\n", "to see what the visualization looks like. This would pop up a window that will\n", "show you the visualization as it is now. You can interact with this\n", "visualization using your mouse to rotate the image in 3D, and mouse+ctrl or\n", "mouse+shift to pan and rotate it in plane, respectively. Use the scroll up and\n", "scroll down in your mouse to zoom in and out. Once you have found a view of\n", "the data that you like, you can close the window (as long as its open, it is\n", "blocking execution of any further commands in the Python interpreter!) and\n", "then you can query your scene for the \"camera settings\" by calling:\n", "\n", "```python\n", "scene.camera_info()\n", "```\n", "\n", "This will print out to the screen something like this:\n", "\n", "```\n", "# Active Camera\n", "Position (238.04, 174.48, 143.04)\n", "Focal Point (96.32, 110.34, 84.48)\n", "View Up (-0.33, -0.12, 0.94)\n", "```\n", "\n", "We can use the information we have gleaned to set the camera on subsequent visualization that use this scene object." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "scene.set_camera(position=(238.04, 174.48, 143.04),\n", " focal_point=(96.32, 110.34, 84.48),\n", " view_up=(-0.33, -0.12, 0.94))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Record the visualization\n", "\n", "If you are not working in an interactive session, or you have already figured\n", "out the camera settings that work for you, you can go ahead and record the\n", "image into a png file. We use a pretty high resolution here (2400 by 2400) so\n", "that 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, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "out_folder = op.join(afq_home, \"viz_results\")\n", "\n", "if not op.exists(op.join(out_folder, 'arc_cst1.png')):\n", " window.record(\n", " scene=scene,\n", " out_path=op.join(out_folder, 'arc_cst1.png'),\n", " size=(2400, 2400))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting bundle colors\n", "\n", "The default behavior produces a nice image! But we often want to be able to\n", "differentiate streamlines based on the bundle to which they belong. For this,\n", "we will set the color of each streamline, based on the bundle. We often use\n", "the 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\n", "options to pass here is an RGB triplet, which will determine the color of all\n", "of the streamlines in that actor. We can get the Tableau20 RGB triplets from\n", "Matplotlib's colormap module. We initialize our bundle actors again and clear\n", "the scene before adding these new actors and adding back the slice actors. We\n", "then call `record` to create this new figure.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "from matplotlib.cm import tab20\n", "color_arc = tab20.colors[18]\n", "color_cst = tab20.colors[2]\n", "\n", "arc_actor = lines_as_tubes(arc_t1w, 8, colors=color_arc)\n", "cst_actor = lines_as_tubes(cst_t1w, 8, colors=color_cst)\n", "\n", "scene.clear()\n", "\n", "scene.add(arc_actor)\n", "scene.add(cst_actor)\n", "for slicer in slicers:\n", " scene.add(slicer)\n", "\n", "if not op.exists(op.join(out_folder, 'arc_cst2.png')):\n", " window.record(\n", " scene=scene,\n", " out_path=op.join(out_folder, 'arc_cst2.png'),\n", " size=(2400, 2400))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding core bundles with tract profiles\n", "\n", "Next, we'd like to add a representation of the core of each bundle and plot\n", "the profile of some property along the length of the bundle. Another option\n", "for setting streamline colors is for each point along a streamline to have a a\n", "different color. Here, we will do this with only one streamline, which\n", "represents the core of the fiber -- it's median trajectory, and with a profile\n", "of the FA. This requires a few distinct steps. The first is adding another\n", "actor for each of the core fibers. We determine the core bundle as the median\n", "of the coordinates of the streamlines after each streamline is resampled to\n", "100 equi-distant points. In the next step, we extract the FA profiles for each\n", "of the bundles, using the `afq_profile` function. Before we do this, we make\n", "sure that the streamlines are back in the voxel coordinate frame of the FA\n", "volume. The core bundle and FA profile are put together into a new line actor,\n", "which we initialize as a thick tube (line width of 40) and with a call to\n", "create a colormap from the FA profiles (we can choose which color-map to use\n", "here. In this case, we chose `'viridis'`). Note that because we are passing\n", "only one streamline into the line actor, we have to pass it inside of square\n", "brackets (`[]`). This is because the `actor.line` initializer expects a\n", "sequence of streamlines as input. Before plotting things again, we initialize\n", "our bundle line actors again. This time, we set the opacity of the lines to\n", "0.1, so that they do not occlude the view of the core bundle visualizations.\n", "\n", ":::{note}\n", "In this case, the profile is FA, but you can use any sequence of 100 numbers in place of the FA profiles: group differences, p-values, etc.\n", "::" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "from dipy.tracking.streamline import set_number_of_points\n", "core_arc = np.median(np.asarray(set_number_of_points(arc_t1w, 100)), axis=0)\n", "core_cst = np.median(np.asarray(set_number_of_points(cst_t1w, 100)), axis=0)\n", "\n", "\n", "from dipy.stats.analysis import afq_profile\n", "sft_arc.to_vox()\n", "sft_cst.to_vox()\n", "arc_profile = afq_profile(fa, sft_arc.streamlines, affine=np.eye(4))\n", "cst_profile = afq_profile(fa, sft_cst.streamlines, affine=np.eye(4))\n", "\n", "core_arc_actor = lines_as_tubes(\n", " [core_arc],\n", " 40,\n", " colors=create_colormap(arc_profile, 'viridis')\n", ")\n", "\n", "core_cst_actor = lines_as_tubes(\n", " [core_cst],\n", " 40,\n", " colors=create_colormap(cst_profile, 'viridis')\n", ")\n", "\n", "scene.clear()\n", "\n", "arc_actor = lines_as_tubes(arc_t1w, 8, colors=color_arc, opacity=0.1)\n", "cst_actor = lines_as_tubes(cst_t1w, 8, colors=color_cst, opacity=0.1)\n", "\n", "scene.add(arc_actor)\n", "scene.add(cst_actor)\n", "for slicer in slicers:\n", " scene.add(slicer)\n", "scene.add(core_arc_actor)\n", "scene.add(core_cst_actor)\n", "\n", "if not op.exists(op.join(out_folder, 'arc_cst3.png')):\n", " window.record(\n", " scene=scene,\n", " out_path=op.join(out_folder, 'arc_cst3.png'),\n", " size=(2400, 2400))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding ROIs\n", "Another element that we can add into the visualization are volume renderings\n", "of regions of interest. For example, the waypoint ROIs that were used to\n", "define the bundles. Here, we will add the waypoints that were used to define\n", "the arcuate fasciculus in this subject. We start by reading this data in.\n", "Because it is represented in the space of the diffusion data, it too needs to\n", "be resampled into the T1-weighted space, before being visualized. After\n", "resampling, we might have some values between 0 and 1 because of the\n", "interpolation from the low resolution of the diffusion into the high\n", "resolution of the T1-weighted. We will include in the volume rendering any\n", "values larger than 0. The main change from the previous visualizations is the\n", "adition of a `contour_from_roi` actor for each of the ROIs. We select another\n", "color from the Tableau 20 palette to represent this, and use an opacity of\n", "0.5.\n", "\n", ":::{note}\n", "In this case, the surface rendered is of the waypoint ROIs, but very similar code can be used to render other surfaces, provided a volume that contains that surface. For example, the gray-white matter boundary could be visualized provided an array with a binary representation of the volume `#` enclosed in this boundary.\n", ":::\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "from dipy.align import resample\n", "\n", "waypoint1 = 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", "\n", "waypoint2 = 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", "\n", "waypoint1_xform = resample(waypoint1, t1w_img)\n", "waypoint2_xform = resample(waypoint2, t1w_img)\n", "waypoint1_data = waypoint1_xform.get_fdata() > 0\n", "waypoint2_data = waypoint2_xform.get_fdata() > 0\n", "\n", "scene.clear()\n", "\n", "arc_actor = lines_as_tubes(arc_t1w, 8, colors=color_arc)\n", "cst_actor = lines_as_tubes(cst_t1w, 8, colors=color_cst)\n", "\n", "scene.add(arc_actor)\n", "scene.add(cst_actor)\n", "for slicer in slicers:\n", " scene.add(slicer)\n", "\n", "surface_color = tab20.colors[0]\n", "\n", "waypoint1_actor = actor.contour_from_roi(waypoint1_data,\n", " color=surface_color,\n", " opacity=0.5)\n", "\n", "waypoint2_actor = actor.contour_from_roi(waypoint2_data,\n", " color=surface_color,\n", " opacity=0.5)\n", "\n", "\n", "scene.add(waypoint1_actor)\n", "scene.add(waypoint2_actor)\n", "\n", "\n", "if not op.exists(op.join(out_folder, 'arc_cst4.png')):\n", " window.record(\n", " scene,\n", " out_path=op.join(out_folder, 'arc_cst4.png'),\n", " size=(2400, 2400))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Making a Figure out of many fury panels\n", "\n", "We can also make a figure that contains multiple panels, each of which\n", "contains a different visualization. This is useful for communicating the\n", "results of an analysis. Here, we will make a figure with four panels, using\n", "some of the visualizations we have already created. We will use some\n", "convenient methods from pyAFQ.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "pf = PanelFigure(3, 2, 6, 9)\n", "pf.add_img(op.join(out_folder, 'arc_cst1.png'), 0, 0)\n", "pf.add_img(op.join(out_folder, 'arc_cst2.png'), 1, 0)\n", "pf.add_img(op.join(out_folder, 'arc_cst3.png'), 0, 1)\n", "pf.add_img(op.join(out_folder, 'arc_cst4.png'), 1, 1)\n", "if not op.exists(op.join(out_folder, 'arc_cst_fig.png')):\n", " pf.format_and_save_figure(op.join(out_folder, 'arc_cst_fig.png'))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.12.4" } }, "nbformat": 4, "nbformat_minor": 4 }