2D Demons Registration

demons2DDemo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      ": package ‘purrr’ was built under R version 3.2.5"
     ]
    }
   ],
   "source": [
    "library(SimpleITK)\n",
    "library(purrr)\n",
    "\n",
    "# Create the fixed image\n",
    "m1 <- matrix(0, 256, 256) \n",
    "m1[125:150, 130:170] <- 255\n",
    "fixed_image <- as.image(m1)\n",
    "\n",
    "# Create the moving image\n",
    "m2 <- matrix(0, 256, 256) \n",
    "m2[130:170, 130:170] <- 255\n",
    "moving_image <- as.image(m2)\n",
    "\n",
    "# Create the grid image to display the deformation, the sigma defines the \"width\" of the lines to avoid aliasing when displayed\n",
    "grid <- GridSource(outputPixelType=\"sitkUInt16\", size=c(256,256), sigma=c(0.3,0.3), gridSpacing=c(10.0,10.0))\n",
    "\n",
    "Show(fixed_image, \"fixed\")\n",
    "Show(moving_image, \"moving\")\n",
    "Show(grid, \"grid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#    \n",
    "# Args:\n",
    "#        image: The image we want to resample.\n",
    "#        shrink_factor: A number greater than one, such that the new image's size is original_size/shrink_factor.\n",
    "#        smoothing_sigma: Sigma for Gaussian smoothing, this is in physical (image spacing) units, not pixels.\n",
    "#    Return:\n",
    "#        Image which is a result of smoothing the input and then resampling it using the given sigma and shrink factor.\n",
    "#\n",
    "smooth_and_resample <- function(image, shrink_factor, smoothing_sigma)\n",
    "{\n",
    "    smoothed_image <- SmoothingRecursiveGaussian(image, smoothing_sigma)\n",
    "    \n",
    "    original_spacing <- image$GetSpacing()\n",
    "    original_size <- image$GetSize()\n",
    "    new_size <-  as.integer(round(original_size/shrink_factor))\n",
    "    new_spacing <- (original_size-1)*original_spacing/(new_size-1)\n",
    "\n",
    "    return(Resample(smoothed_image, new_size, Transform(), \n",
    "                    \"sitkLinear\", image$GetOrigin(),\n",
    "                    new_spacing, image$GetDirection(), 0.0, \n",
    "                    image$GetPixelID()))\n",
    "}\n",
    "\n",
    "#    \n",
    "# Run the given registration algorithm in a multiscale fashion. The original scale should not be given as input as the\n",
    "# original images are implicitly incorporated as the base of the pyramid.\n",
    "# Args:\n",
    "#   registration_algorithm: Any registration algorithm that has an Execute(fixed_image, moving_image, displacement_field_image)\n",
    "#                           method.\n",
    "#   fixed_image: Resulting transformation maps points from this image's spatial domain to the moving image spatial domain.\n",
    "#   moving_image: Resulting transformation maps points from the fixed_image's spatial domain to this image's spatial domain.\n",
    "#   initial_transform: Any SimpleITK transform, used to initialize the displacement field.\n",
    "#   shrink_factors: Shrink factors relative to the original image's size.\n",
    "#   smoothing_sigmas: Amount of smoothing which is done prior to resampling the image using the given shrink factor. These\n",
    "#                     are in physical (image spacing) units.\n",
    "# Returns: \n",
    "#    DisplacementFieldTransform\n",
    "#\n",
    "multiscale_demons <- function(registration_algorithm, fixed_image, moving_image, initial_transform = NULL, \n",
    "                              shrink_factors=NULL, smoothing_sigmas=NULL)\n",
    "{    \n",
    "    # Create image pyramids. \n",
    "    fixed_images <- c(fixed_image, \n",
    "                      if(!is.null(shrink_factors))\n",
    "                          map2(rev(shrink_factors), rev(smoothing_sigmas), \n",
    "                               ~smooth_and_resample(fixed_image, .x, .y))\n",
    "                      )\n",
    "    moving_images <- c(moving_image, \n",
    "                       if(!is.null(shrink_factors))\n",
    "                           map2(rev(shrink_factors), rev(smoothing_sigmas), \n",
    "                               ~smooth_and_resample(moving_image, .x, .y))\n",
    "                       )\n",
    "\n",
    "    # Uncomment the following two lines if you want to see your image pyramids.\n",
    "    #lapply(fixed_images, Show)\n",
    "    #lapply(moving_images, Show)\n",
    "    \n",
    "                              \n",
    "    # Create initial displacement field at lowest resolution. \n",
    "    # Currently, the pixel type is required to be sitkVectorFloat64 because of a constraint imposed by the Demons filters.\n",
    "    lastImage <- fixed_images[[length(fixed_images)]]\n",
    "    if(!is.null(initial_transform))\n",
    "    {\n",
    "        initial_displacement_field <- TransformToDisplacementField(initial_transform, \n",
    "                                                                  \"sitkVectorFloat64\",\n",
    "                                                                  lastImage$GetSize(),\n",
    "                                                                  lastImage$GetOrigin(),\n",
    "                                                                  lastImage$GetSpacing(),\n",
    "                                                                  lastImage$GetDirection())\n",
    "    }\n",
    "    else\n",
    "    {\n",
    "        initial_displacement_field <- Image(lastImage$GetWidth(), \n",
    "                                            lastImage$GetHeight(),                                            \n",
    "                                            \"sitkVectorFloat64\")\n",
    "        initial_displacement_field$CopyInformation(lastImage)\n",
    "    }\n",
    "    # Run the registration pyramid, run a registration at the top of the pyramid and then iterate: \n",
    "    # a. resampling previous deformation field onto higher resolution grid.\n",
    "    # b. register.\n",
    "    initial_displacement_field <- registration_algorithm$Execute(fixed_images[[length(fixed_images)]], \n",
    "                                                                moving_images[[length(moving_images)]], \n",
    "                                                                initial_displacement_field)\n",
    "    # This is a use case for a loop, because the operations depend on the previous step. Otherwise\n",
    "    # we need to mess around with tricky assignments to variables in different scopes\n",
    "    for (idx in seq(length(fixed_images)-1,1)) {\n",
    "        f_image <- fixed_images[[idx]]\n",
    "        m_image <- moving_images[[idx]]\n",
    "        initial_displacement_field <- Resample(initial_displacement_field, f_image)\n",
    "        initial_displacement_field <- registration_algorithm$Execute(f_image, m_image, initial_displacement_field)\n",
    "    }\n",
    "\n",
    "    return(DisplacementFieldTransform(initial_displacement_field))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Select a Demons filter and configure it.\n",
    "demons_filter <-  DiffeomorphicDemonsRegistrationFilter()\n",
    "demons_filter$SetNumberOfIterations(20)\n",
    "# Regularization (update field - viscous, total field - elastic).\n",
    "demons_filter$SetSmoothDisplacementField(TRUE)\n",
    "demons_filter$SetStandardDeviations(2.0)\n",
    "\n",
    "# Run the registration.\n",
    "tx <- multiscale_demons(demons_filter, \n",
    "                       fixed_image, \n",
    "                       moving_image,\n",
    "                       shrink_factors=c(4,2),\n",
    "                       smoothing_sigmas=c(8,4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "resampled_moving_image <- Resample(moving_image, fixed_image, tx)\n",
    "resampled_grid_image <- Resample(grid, fixed_image, tx)\n",
    "\n",
    "Show(resampled_moving_image, \"resampled moving\")\n",
    "Show(resampled_grid_image, \"resampled grid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}