{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Applying image-based deformations to streamlines\n\nThis example shows how to register streamlines into a template space by\napplying non-rigid deformations.\n\nAt times we will be interested in bringing a set of streamlines into some\ncommon, reference space to compute statistics out of the registered\nstreamlines. For a discussion on the effects of spatial normalization\napproaches on tractography the work by Green et al. [Greene17]_ can be read.\n\nFor brevity, we will include in this example only streamlines going through\nthe corpus callosum connecting left to right superior frontal cortex. The\nprocess of tracking and finding these streamlines is fully demonstrated in\nthe `streamline_tools` example. If this example has been run, we can read\nthe streamlines from file. Otherwise, we'll run that example first, by\nimporting it. This provides us with all of the variables that were created in\nthat example.\n\nIn order to get the deformation field, we will first use two b0 volumes. Both\nmoving and static images are assumed to be in RAS. The first one will be the\nb0 from the Stanford HARDI dataset:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from os.path import join as pjoin\n\nimport numpy as np\nimport nibabel as nib\n\nfrom dipy.core.gradients import gradient_table\nfrom dipy.data import get_fnames\nfrom dipy.io.gradients import read_bvals_bvecs\nfrom dipy.io.image import load_nifti_data, load_nifti, save_nifti\n\nhardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')\n\ndata, affine, hardi_img = load_nifti(hardi_fname, return_img=True)\nvox_size = hardi_img.header.get_zooms()[0]\nbvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)\ngtab = gradient_table(bvals, bvecs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The second one will be the T2-contrast MNI template image, which we'll need to\nreslice to 2x2x2 mm isotropic voxel resolution to match the HARDI data.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.data.fetcher import (fetch_mni_template, read_mni_template)\nfrom dipy.align.reslice import reslice\n\nfetch_mni_template()\nimg_t2_mni = read_mni_template(version=\"a\", contrast=\"T2\")\n\nnew_zooms = (2., 2., 2.)\ndata2, affine2 = reslice(np.asarray(img_t2_mni.dataobj), img_t2_mni.affine,\n img_t2_mni.header.get_zooms(), new_zooms)\nimg_t2_mni = nib.Nifti1Image(data2, affine=affine2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We filter the diffusion data from the Stanford HARDI dataset to find all the b0\nimages.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "b0_idx_stanford = np.where(gtab.b0s_mask)[0]\nb0_data_stanford = data[..., b0_idx_stanford]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We then remove the skull from them:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.segment.mask import median_otsu\n\nb0_masked_stanford, _ = median_otsu(b0_data_stanford,\n vol_idx=list(range(b0_data_stanford.shape[-1])),\n median_radius=4, numpass=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And go on to compute the Stanford HARDI dataset mean b0 image.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "mean_b0_masked_stanford = np.mean(b0_masked_stanford, axis=3,\n dtype=data.dtype)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will register the mean b0 to the MNI T2 image template non-rigidly to\nobtain the deformation field that will be applied to the streamlines. This is\njust one of the strategies that can be used to obtain an appropriate\ndeformation field. Other strategies include computing an FA template map as\nthe static image, and registering the FA map of the moving image to it. This\nmay may eventually lead to results with improved accuracy, since a T2-contrast\ntemplate image as the target for normalization does not provide optimal tissue\ncontrast for maximal SyN performance.\n\nWe will first perform an affine registration to roughly align the two volumes:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.align.imaffine import (MutualInformationMetric, AffineRegistration,\n transform_origins)\nfrom dipy.align.transforms import (TranslationTransform3D, RigidTransform3D,\n AffineTransform3D)\n\nstatic = np.asarray(img_t2_mni.dataobj)\nstatic_affine = img_t2_mni.affine\nmoving = mean_b0_masked_stanford\nmoving_affine = hardi_img.affine" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We estimate an affine that maps the origin of the moving image to that of the\nstatic image. We can then use this later to account for the offsets of each\nimage.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "affine_map = transform_origins(static, static_affine, moving, moving_affine)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We specify the mismatch metric:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "nbins = 32\nsampling_prop = None\nmetric = MutualInformationMetric(nbins, sampling_prop)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As well as the optimization strategy:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "level_iters = [10, 10, 5]\nsigmas = [3.0, 1.0, 0.0]\nfactors = [4, 2, 1]\naffine_reg = AffineRegistration(metric=metric, level_iters=level_iters,\n sigmas=sigmas, factors=factors)\ntransform = TranslationTransform3D()\n\nparams0 = None\ntranslation = affine_reg.optimize(static, moving, transform, params0,\n static_affine, moving_affine)\ntransformed = translation.transform(moving)\ntransform = RigidTransform3D()\n\nrigid_map = affine_reg.optimize(static, moving, transform, params0,\n static_affine, moving_affine,\n starting_affine=translation.affine)\ntransformed = rigid_map.transform(moving)\ntransform = AffineTransform3D()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We bump up the iterations to get a more exact fit:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "affine_reg.level_iters = [1000, 1000, 100]\nhighres_map = affine_reg.optimize(static, moving, transform, params0,\n static_affine, moving_affine,\n starting_affine=rigid_map.affine)\ntransformed = highres_map.transform(moving)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now perform the non-rigid deformation using the Symmetric Diffeomorphic\nRegistration (SyN) Algorithm proposed by Avants et al. [Avants09]_ (also\nimplemented in the ANTs software [Avants11]_):\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.align.imwarp import SymmetricDiffeomorphicRegistration\nfrom dipy.align.metrics import CCMetric\n\nmetric = CCMetric(3)\nlevel_iters = [10, 10, 5]\nsdr = SymmetricDiffeomorphicRegistration(metric, level_iters)\n\nmapping = sdr.optimize(static, moving, static_affine, moving_affine,\n highres_map.affine)\nwarped_moving = mapping.transform(moving)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We show the registration result with:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.viz import regtools\n\nregtools.overlay_slices(static, warped_moving, None, 0, 'Static', 'Moving',\n 'transformed_sagittal.png')\nregtools.overlay_slices(static, warped_moving, None, 1, 'Static', 'Moving',\n 'transformed_coronal.png')\nregtools.overlay_slices(static, warped_moving, None, 2, 'Static', 'Moving',\n 'transformed_axial.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ".. figure:: transformed_sagittal.png\n :align: center\n.. figure:: transformed_coronal.png\n :align: center\n.. figure:: transformed_axial.png\n :align: center\n\n Deformable registration result.\n\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now fetch a set of streamlines from the Stanford HARDI dataset.\nThose streamlines were generated during the `streamline_tools` example.\n\nWe read the streamlines from file in voxel space:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.data import fetch_stanford_tracks\nfrom dipy.io.streamline import load_tractogram\n\n\nstreamlines_files = fetch_stanford_tracks()\nlr_superiorfrontal_path = pjoin(streamlines_files[1],\n 'hardi-lr-superiorfrontal.trk')\n\nsft = load_tractogram(lr_superiorfrontal_path, 'same')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We then apply the obtained deformation field to the streamlines. Note that the\nprocess can be sensitive to image orientation and voxel resolution. Thus, we\nfirst apply the non-rigid warping and simultaneously apply a computed rigid\naffine transformation whose extents must be corrected to account for the\ndifferent voxel grids of the moving and static images.\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.tracking.streamline import deform_streamlines\n\n# Create an isocentered affine\ntarget_isocenter = np.diag(np.array([-vox_size, vox_size, vox_size, 1]))\n\n# Take the off-origin affine capturing the extent contrast between the mean b0\n# image and the template\norigin_affine = affine_map.affine.copy()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to align the FOV of the template and the mirror image of the\nstreamlines, we first need to flip the sign on the x-offset and y-offset so\nthat we get the mirror image of the forward deformation field.\n\nWe need to use the information about the origin offsets (i.e. between the\nstatic and moving images) that we obtained using :meth:`transform_origins`:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "origin_affine[0][3] = -origin_affine[0][3]\norigin_affine[1][3] = -origin_affine[1][3]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ":meth:`transform_origins` returns this affine transformation with (1, 1, 1)\nzooms and not (2, 2, 2), which means that the offsets need to be scaled by 2.\nThus, we scale z by the voxel size:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "origin_affine[2][3] = origin_affine[2][3]/vox_size" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But when scaling the z-offset, we are also implicitly scaling the y-offset as\nwell (by 1/2).Thus we need to correct for this by only scaling the y by the\nsquare of the voxel size (1/4, and not 1/2):\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "origin_affine[1][3] = origin_affine[1][3]/vox_size**2\n\n# Apply the deformation and correct for the extents\nmni_streamlines = deform_streamlines(\n sft.streamlines, deform_field=mapping.get_forward_field(),\n stream_to_current_grid=target_isocenter,\n current_grid_to_world=origin_affine, stream_to_ref_grid=target_isocenter,\n ref_grid_to_world=np.eye(4))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We display the original streamlines and the registered streamlines:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.viz import has_fury\n\n\ndef show_template_bundles(bundles, show=True, fname=None):\n\n scene = window.Scene()\n template_actor = actor.slicer(static)\n scene.add(template_actor)\n\n lines_actor = actor.streamtube(bundles, window.colors.orange,\n linewidth=0.3)\n scene.add(lines_actor)\n\n if show:\n window.show(scene)\n if fname is not None:\n window.record(scene, n_frames=1, out_path=fname, size=(900, 900))\n\n\nif has_fury:\n\n from fury import actor, window\n\n show_template_bundles(mni_streamlines, show=False,\n fname='streamlines_DSN_MNI.png')\n\n \"\"\"\n .. figure:: streamlines_DSN_MNI.png\n :align: center\n\n Streamlines before and after registration.\n\n The corpus callosum bundles have been deformed to adapt to the MNI\n template space.\n\n \"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we save the registered streamlines:\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from dipy.io.stateful_tractogram import Space, StatefulTractogram\nfrom dipy.io.streamline import save_tractogram\n\nsft = StatefulTractogram(mni_streamlines, img_t2_mni, Space.RASMM)\n\nsave_tractogram(sft, 'mni-lr-superiorfrontal.trk', bbox_valid_check=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n\n.. [Avants09] Avants, B. B., Epstein, C. L., Grossman, M., & Gee, J. C. (2009).\n Symmetric Diffeomorphic Image Registration with Cross-Correlation:\n Evaluating Automated Labeling of Elderly and Neurodegenerative Brain, 12(1),\n 26-41.\n\n.. [Avants11] Avants, B. B., Tustison, N., & Song, G. (2011). Advanced\n Normalization Tools (ANTS), 1-35.\n\n.. [Greene17] Greene, C., Cieslak, M., & Grafton, S. T. (2017). Effect of\n different spatial normalization approaches on tractography and structural\n brain networks. Network Neuroscience, 1-19.\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 }