{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Diffeomorphic Registration with binary and fuzzy images\n\nThis example demonstrates registration of a binary and a fuzzy image.\nThis could be seen as aligning a fuzzy (sensed) image to a binary\n(e.g., template) image.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom skimage import draw, filters\nfrom dipy.align.imwarp import SymmetricDiffeomorphicRegistration\nfrom dipy.align.metrics import SSDMetric\nfrom dipy.viz import regtools"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's generate a sample template image as the combination of three ellipses.\nWe will generate the fuzzy (sensed) version of the image by smoothing\nthe reference image.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def draw_ellipse(img, center, axis):\n    rr, cc = draw.ellipse(center[0], center[1], axis[0], axis[1],\n                          shape=img.shape)\n    img[rr, cc] = 1\n    return img\n\n\nimg_ref = np.zeros((64, 64))\nimg_ref = draw_ellipse(img_ref, (25, 15), (10, 5))\nimg_ref = draw_ellipse(img_ref, (20, 45), (15, 10))\nimg_ref = draw_ellipse(img_ref, (50, 40), (7, 15))\n\nimg_in = filters.gaussian(img_ref, sigma=3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's define a small visualization function.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def show_images(img_ref, img_warp, fig_name):\n    fig, axarr = plt.subplots(ncols=2, figsize=(12, 5))\n    axarr[0].set_title('warped image & reference contour')\n    axarr[0].imshow(img_warp)\n    axarr[0].contour(img_ref, colors='r')\n    ssd = np.sum((img_warp - img_ref) ** 2)\n    axarr[1].set_title('difference, SSD=%.02f' % ssd)\n    im = axarr[1].imshow(img_warp - img_ref)\n    plt.colorbar(im)\n    fig.tight_layout()\n    fig.savefig(fig_name + '.png')\n\n\nshow_images(img_ref, img_in, 'input')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: input.png\n   :align: center\n\n   Input images before alignment.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's use the general Registration function with some naive parameters,\nsuch as set `step_length` as 1 assuming maximal step 1 pixel and a reasonably\nsmall number of iterations since the deformation with already aligned images\nshould be minimal.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sdr = SymmetricDiffeomorphicRegistration(metric=SSDMetric(img_ref.ndim),\n                                         step_length=1.0,\n                                         level_iters=[50, 100],\n                                         inv_iter=50,\n                                         ss_sigma_factor=0.1,\n                                         opt_tol=1.e-3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Perform the registration with equal images.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mapping = sdr.optimize(img_ref.astype(float), img_ref.astype(float))\nimg_warp = mapping.transform(img_ref, 'linear')\nshow_images(img_ref, img_warp, 'output-0')\nregtools.plot_2d_diffeomorphic_map(mapping, 5, 'map-0.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: output-0.png\n   :align: center\n.. figure:: map-0.png\n   :align: center\n\n   Registration results for default parameters and equal images.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Perform the registration with binary and fuzzy images.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mapping = sdr.optimize(img_ref.astype(float), img_in.astype(float))\nimg_warp = mapping.transform(img_in, 'linear')\nshow_images(img_ref, img_warp, 'output-1')\nregtools.plot_2d_diffeomorphic_map(mapping, 5, 'map-1.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: output-1.png\n   :align: center\n.. figure:: map-1.png\n   :align: center\n\n   Registration results for a naive parameter configuration.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Note, we are still using a multi-scale approach which makes `step_length`\nin the upper level multiplicatively larger.\nWhat happens if we set `step_length` to a rather small value?\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sdr.step_length = 0.1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Perform the registration and examine the output.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mapping = sdr.optimize(img_ref.astype(float), img_in.astype(float))\nimg_warp = mapping.transform(img_in, 'linear')\nshow_images(img_ref, img_warp, 'output-2')\nregtools.plot_2d_diffeomorphic_map(mapping, 5, 'map-2.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: output-2.png\n   :align: center\n.. figure:: map-2.png\n   :align: center\n\n   Registration results for decreased step size.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "An alternative scenario is to use just a single-scale level.\nEven though the warped image may look fine, the estimated deformations show\nthat it is off the mark.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sdr = SymmetricDiffeomorphicRegistration(metric=SSDMetric(img_ref.ndim),\n                                         step_length=1.0,\n                                         level_iters=[100],\n                                         inv_iter=50,\n                                         ss_sigma_factor=0.1,\n                                         opt_tol=1.e-3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Perform the registration.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mapping = sdr.optimize(img_ref.astype(float), img_in.astype(float))\nimg_warp = mapping.transform(img_in, 'linear')\nshow_images(img_ref, img_warp, 'output-3')\nregtools.plot_2d_diffeomorphic_map(mapping, 5, 'map-3.png')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. figure:: output-3.png\n   :align: center\n.. figure:: map-3.png\n   :align: center\n\n   Registration results for single level.\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
}