{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# R syntactic sugar\n",
    "\n",
    "SimpleITK in R attempts to make all the common R array operators and arithmetic work on images.\n",
    "\n",
    "First we load SimpleITK and set up a viewer for notebooks, using the jet colourmap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "library(SimpleITK)\n",
    "\n",
    "Dwidth <- grid::unit(10, \"cm\")\n",
    "## version of show using a jet colourmap.\n",
    "setMethod('show', '_p_itk__simple__Image',\n",
    "          function(object)\n",
    "          {\n",
    "            jet.colors <- colorRampPalette(c(\"#00007F\", \"blue\", \"#007FFF\", \"cyan\", \"#7FFF7F\", \"yellow\", \"#FF7F00\", \"red\", \"#7F0000\"))\n",
    "            a <- t(as.array(object))\n",
    "            rg <- range(a)\n",
    "            A <- (a-rg[1])/(rg[2]-rg[1])\n",
    "            dd <- dim(a)\n",
    "            sp <- object$GetSpacing()\n",
    "            sz <- object$GetSize()\n",
    "            worlddim <- sp*sz\n",
    "            worlddim <- worlddim/worlddim[1]\n",
    "            W <- Dwidth\n",
    "            H <- Dwidth * worlddim[2]\n",
    "            jc <- jet.colors(101)\n",
    "            A <- jc[A*100 + 1]               \n",
    "            dim(A) <- dim(a)\n",
    "            grid::grid.raster(A,default.units=\"mm\", width=W, height=H)\n",
    "\n",
    "          }\n",
    "\n",
    "          )\n",
    "## Make a smaller plot size\n",
    "options(repr.plot.height = 3, repr.plot.width=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img = GaussianSource(\"sitkUInt8\", size=c(64,64)) \n",
    "img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we note a problem with the swig bindings to the procedural interface - dispatching and defaults are not as nice as the python version. This means that we can't say\n",
    "```r\n",
    "GaborSource(\"sitkFloat32\", size=c(64,64), frequency=.03)\n",
    "```\n",
    "because there are arguments missing before _frequency_ and swig can't generate default arguments from c++ code.\n",
    "\n",
    "## Image slicing and indexing\n",
    "Lets start by accessing a single pixel - note that indexing starts from 1:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img[25,25]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cropping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img[17:49,]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Don't have the \"from the end\" slicing - negatives mean remove \n",
    "\n",
    "img[,16:(img$GetHeight() - 16)]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img = GaborSource(\"sitkFloat32\", size=c(64,64), sigma=c(16, 16), mean=c(32,32), frequency=.03)\n",
    "img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img[, -(1:16)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img[1:32, 1:32]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flipping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "corner=img[1:32, 1:32]\n",
    "corner[corner$GetWidth():1,]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tiling (using slicing)\n",
    "\n",
    "We cannot perform tiling via slicing directly on an image. The reason for this restriction is that an image is not equivalent to an array. Unlike an array, an image enforces the concept of physical spacing between pixels/voxels, this spacing has to be uniform for each of the axes. Thus a non-uniform slicing, e.g. (1,2,3,3,2,1), will generate an error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_expected": "spacing is not uniform"
   },
   "outputs": [],
   "source": [
    "corner[c(1:corner$GetWidth(), corner$GetWidth():1), c(1:corner$GetHeight(), corner$GetHeight():1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can still easily perform the tiling operation we want by slightly abusing the system, converting the image to an array, performing the tiling and then converting the array back to an image. Note that the new image will have unit spacings for each axis, the concept of physical spacing is lost once we convert to an array.\n",
    "\n",
    "The correct, maintaining the data as an image, way of performing this operation is to use a SimpleITK filter as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "arr = as.array(corner)\n",
    "as.image(arr[c(1:corner$GetWidth(), corner$GetWidth():1), c(1:corner$GetHeight(), corner$GetHeight():1)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tiling (with SimpleITK filters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Tile(corner, corner[corner$GetWidth():1,],corner[,corner$GetHeight():1],corner[corner$GetWidth():1, corner$GetHeight():1], c(2,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Slice extraction\n",
    "Extracting a single slice of a 3D image is a natural extension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "img3d = GaborSource(\"sitkFloat32\", size=rep(64,3), sigma=rep(16, 3), mean=rep(32,3), frequency=.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img3d[,,32]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img3d[16,,]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subsampling\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img3d[,seq(from=1, to=img3d$GetHeight(), by=3), 32]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mathematical Operators\n",
    "Most array operations are overloaded for images - meaning that an ITK filter is used to perform the pixel-wise operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## source the functions for loading data\n",
    "source(\"downloaddata.R\")\n",
    "cthead <- ReadImage(fetch_data(\"cthead1.png\"))\n",
    "img = Cast(cthead,\"sitkFloat32\")\n",
    "img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img[150,150]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "timg <- img**2\n",
    "timg\n",
    "timg[150,150]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "aimg <- img + 100\n",
    "aimg[150,150]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flip the image, x axis. This also changes the x coordinate of the origin. When we try to add the two images we will fail, because they do not occupy the same physical space. Again, SimpleITK images are not arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "simpleitk_error_expected": "Exception in SITK"
   },
   "outputs": [],
   "source": [
    "flipped_img <- img[img$GetWidth():1, ]\n",
    "img + flipped_img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we really want to add the intensity values of the two images, we can modify the meta-data of the flipped image to match that of the original. We know they have the same number of pixels on each axis, so we can do this safely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "flipped_img$CopyInformation(img)\n",
    "img + flipped_img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logic and comparison operators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img > 150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(img > 90) + (img > 150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m <- (img > 90) != (img > 150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img * Cast(m, \"sitkFloat32\") "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## It's about time\n",
    "\n",
    "SimpleITK's implementation of the mathematical and logical operators is computationally efficient. The following cells show that this is not an idle claim, we will time the same operations using SimpleITK images and R arrays.\n",
    "\n",
    "Create three images and equivalent arrays. In this case all images have default values for origin (0,0,0), spacing (1,1,1) and cosine matrix (I)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sz <- 512\n",
    "\n",
    "# Create a sz^3 image whose pixels are all set to one, then get the equivalent R array.\n",
    "im1 <- Image(sz, sz, sz, \"sitkFloat64\") + 1\n",
    "arr1 <- as.array(im1)\n",
    "\n",
    "# Create an sz^3 array whose pixels are set to the sequence 1..sz^3, then get the equivalent SimpleITK image.\n",
    "arr2 <- array(seq(1, sz^3, 1.0), c(sz, sz, sz))\n",
    "im2 <- as.image(arr2)\n",
    "\n",
    "# Create a sz^3 image whose pixels are all set to two, then get the equivalent R array.\n",
    "im3 <- Image(sz, sz, sz, \"sitkFloat64\") + 2\n",
    "arr3 <- as.array(im3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Time the same mathematical and logical operations using SimpleITK and R. Results are in seconds. Interpret the results according to the R documentation: \"The ‘user time’ is the CPU time charged for the execution of user instructions of the calling process. The ‘system time’ is the CPU time charged for execution by the system on behalf of the calling process.\"\n",
    "\n",
    "What we really care about is 'elapsed time' the wall clock. Note that the CPU usage for the SimpleITK code is higher than for the R code, due to the underlying multi-threaded implementation, but the elapsed time is lower.\n",
    "\n",
    "As an exercise, check that the two expressions below indeed return the same results. Note that the SimpleITK logical operator results are 0,1 , while the R results are TRUE, FALSE. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "system.time(((im1 * im2) / im3) > 24)\n",
    "system.time(((arr1 * arr2) / arr3) > 24)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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": 0
}