{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Symmetric Diffeomorphic Registration in 3D\nThis example explains how to register 3D volumes using the Symmetric\nNormalization (SyN) algorithm proposed by Avants et al. [Avants09]_\n(also implemented in the ANTs software [Avants11]_)\n\nWe will register two 3D volumes from the same modality using SyN with the Cross\n-Correlation (CC) metric.\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 CCMetric\nfrom dipy.data import get_fnames\nfrom dipy.io.image import load_nifti\nfrom dipy.viz import regtools"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's fetch two b0 volumes, the first one will be the b0 from the Stanford\nHARDI dataset\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')\n\nstanford_b0, stanford_b0_affine = load_nifti(hardi_fname)\nstanford_b0 = np.squeeze(stanford_b0)[..., 0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The second one will be the same b0 we used for the 2D registration tutorial\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t1_fname, b0_fname = get_fnames('syn_data')\nsyn_b0, syn_b0_affine = load_nifti(b0_fname)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We first remove the skull from the b0's\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.segment.mask import median_otsu\nstanford_b0_masked, stanford_b0_mask = median_otsu(stanford_b0,\n                                                   median_radius=4,\n                                                   numpass=4)\nsyn_b0_masked, syn_b0_mask = median_otsu(syn_b0, median_radius=4, numpass=4)\n\nstatic = stanford_b0_masked\nstatic_affine = stanford_b0_affine\nmoving = syn_b0_masked\nmoving_affine = syn_b0_affine"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Suppose we have already done a linear registration to roughly align the two\nimages\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pre_align = np.array([[1.02783543e+00, -4.83019053e-02, -6.07735639e-02, -2.57654118e+00],\n                      [4.34051706e-03, 9.41918267e-01, -2.66525861e-01, 3.23579799e+01],\n                      [5.34288908e-02, 2.90262026e-01, 9.80820307e-01, -1.46216651e+01],\n                      [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As we did in the 2D example, we would like to visualize (some slices of) the\ntwo volumes by overlapping them over two channels of a color image. To do that\nwe need them to be sampled on the same grid, so let's first re-sample the\nmoving image on the static grid. We create an AffineMap to transform the moving\nimage towards the static image\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from dipy.align.imaffine import AffineMap\naffine_map = AffineMap(pre_align,\n                       static.shape, static_affine,\n                       moving.shape, moving_affine)\n\nresampled = affine_map.transform(moving)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "plot the overlapped middle slices of the volumes\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "regtools.overlay_slices(static, resampled, None, 1, 'Static', 'Moving',\n                        'input_3d.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: input_3d.png\n   :align: center\n\n   Static image in red on top of the pre-aligned moving image (in green).\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We want to find an invertible map that transforms the moving image into the\nstatic image. We will use the Cross-Correlation metric\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "metric = CCMetric(3)"
      ]
    },
    {
      "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 object 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 = [10, 10, 5]\nsdr = SymmetricDiffeomorphicRegistration(metric, level_iters)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Execute the optimization, which returns a DiffeomorphicMap object,\nthat can be used to register images back and forth between the static and\nmoving domains. We provide the pre-aligning matrix that brings the moving\nimage closer to the static image\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mapping = sdr.optimize(static, moving, static_affine, moving_affine, pre_align)"
      ]
    },
    {
      "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)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We plot the overlapped middle slices\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "regtools.overlay_slices(static, warped_moving, None, 1, 'Static',\n                        'Warped moving', 'warped_moving.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: warped_moving.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)\nregtools.overlay_slices(warped_static, moving, None, 1, 'Warped static',\n                        'Moving', 'warped_static.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: warped_static.png\n   :align: center\n\n   Static image transformed under the (inverse) transformation in red on top of\n   the moving image (in green). Note that the moving image has a lower\n   resolution.\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
}