{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Symmetric Diffeomorphic Registration in 2D\nThis example explains how to register 2D images using the Symmetric\nNormalization (SyN) algorithm proposed by Avants et al. [Avants09]_\n(also implemented in the ANTs software [Avants11]_)\n\nWe will perform the classic Circle-To-C experiment for diffeomorphic\nregistration\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nfrom dipy.align.imwarp import SymmetricDiffeomorphicRegistration\nfrom dipy.align.metrics import SSDMetric, CCMetric\nimport dipy.align.imwarp as imwarp\nfrom dipy.data import get_fnames\nfrom dipy.io.image import load_nifti_data\nfrom dipy.segment.mask import median_otsu\nfrom dipy.viz import regtools\n\n\nfname_moving = get_fnames('reg_o')\nfname_static = get_fnames('reg_c')\n\nmoving = np.load(fname_moving)\nstatic = np.load(fname_static)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To visually check the overlap of the static image with the transformed moving\nimage, we can plot them on top of each other with different channels to see\nwhere the differences are located\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "regtools.overlay_images(static, moving, 'Static', 'Overlay', 'Moving',\n                        'input_images.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: input_images.png\n   :align: center\n\n   Input images.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We want to find an invertible map that transforms the moving image (circle)\ninto the static image (the C letter).\n\nThe first decision we need to make is what similarity metric is appropriate\nfor our problem. In this example we are using two binary images, so the Sum\nof Squared Differences (SSD) is a good choice.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dim = static.ndim\nmetric = SSDMetric(dim)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we define an instance of the registration class. The SyN algorithm uses\na multi-resolution approach by building a Gaussian Pyramid. We instruct the\nregistration instance to perform at most $[n_0, n_1, ..., n_k]$ iterations at\neach level of the pyramid. The 0-th level corresponds to the finest resolution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "level_iters = [200, 100, 50, 25]\n\nsdr = SymmetricDiffeomorphicRegistration(metric, level_iters, inv_iter=50)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we execute the optimization, which returns a DiffeomorphicMap object,\nthat can be used to register images back and forth between the static and\nmoving domains\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mapping = sdr.optimize(static, moving)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It is a good idea to visualize the resulting deformation map to make sure the\nresult is reasonable (at least, visually)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "regtools.plot_2d_diffeomorphic_map(mapping, 10, 'diffeomorphic_map.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: diffeomorphic_map.png\n   :align: center\n\n   Deformed lattice under the resulting diffeomorphic map.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's warp the moving image and see if it gets similar to the static image\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "warped_moving = mapping.transform(moving, 'linear')\nregtools.overlay_images(static, warped_moving, 'Static', 'Overlay',\n                        'Warped moving', 'direct_warp_result.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: direct_warp_result.png\n   :align: center\n\n   Moving image transformed under the (direct) transformation in green on top\n   of the static image (in red).\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And we can also apply the inverse mapping to verify that the warped static\nimage is similar to the moving image\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "warped_static = mapping.transform_inverse(static, 'linear')\nregtools.overlay_images(warped_static, moving, 'Warped static', 'Overlay',\n                        'Moving', 'inverse_warp_result.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: inverse_warp_result.png\n   :align: center\n\n   Static image transformed under the (inverse) transformation in red on top\n   of the moving image (in green).\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's register a couple of slices from a b0 image using the Cross\nCorrelation metric. Also, let's inspect the evolution of the registration.\nTo do this we will define a function that will be called by the registration\nobject at each stage of the optimization process. We will draw the current\nwarped images after finishing each resolution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def callback_CC(sdr, status):\n    # Status indicates at which stage of the optimization we currently are\n    # For now, we will only react at the end of each resolution of the scale\n    # space\n    if status == imwarp.RegistrationStages.SCALE_END:\n        # get the current images from the metric\n        wmoving = sdr.metric.moving_image\n        wstatic = sdr.metric.static_image\n        # draw the images on top of each other with different colors\n        regtools.overlay_images(wmoving, wstatic, 'Warped moving', 'Overlay',\n                                'Warped static')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we are ready to configure and run the registration. First load the data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t1_name, b0_name = get_fnames('syn_data')\ndata = load_nifti_data(b0_name)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We first remove the skull from the b0 volume\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "b0_mask, mask = median_otsu(data, median_radius=4, numpass=4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And select two slices to try the 2D registration\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "static = b0_mask[:, :, 40]\nmoving = b0_mask[:, :, 38]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After loading the data, we instantiate the Cross-Correlation metric. The metric\nreceives three parameters: the dimension of the input images, the standard\ndeviation of the Gaussian Kernel to be used to regularize the gradient and the\nradius of the window to be used for evaluating the local normalized cross\ncorrelation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sigma_diff = 3.0\nradius = 4\nmetric = CCMetric(2, sigma_diff, radius)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's use a scale space of 3 levels\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "level_iters = [100, 50, 25]\nsdr = SymmetricDiffeomorphicRegistration(metric, level_iters)\nsdr.callback = callback_CC"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And execute the optimization\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mapping = sdr.optimize(static, moving)\n\nwarped = mapping.transform(moving)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can see the effect of the warping by switching between the images before and\nafter registration\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "regtools.overlay_images(static, moving, 'Static', 'Overlay', 'Moving',\n                        't1_slices_input.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: t1_slices_input.png\n   :align: center\n\n   Input images.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "regtools.overlay_images(static, warped, 'Static', 'Overlay', 'Warped moving',\n                        't1_slices_res.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: t1_slices_res.png\n   :align: center\n\n   Moving image transformed under the (direct) transformation in green on top\n   of the static image (in red).\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And we can apply the inverse warping too\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "inv_warped = mapping.transform_inverse(static)\nregtools.overlay_images(inv_warped, moving, 'Warped static', 'Overlay',\n                        'moving', 't1_slices_res2.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: t1_slices_res2.png\n   :align: center\n\n   Static image transformed under the (inverse) transformation in red on top\n   of the moving image (in green).\n\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, let's see the deformation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "regtools.plot_2d_diffeomorphic_map(mapping, 5, 'diffeomorphic_map_b0s.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: diffeomorphic_map_b0s.png\n   :align: center\n\n   Deformed lattice under the resulting diffeomorphic map.\n\n## 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.. 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
}