{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# N-dimensional image processing\n",
    "This tutorial covers how to access and manipulate image sample values directly using [ImgLib2](https://imagej.net/ImgLib2)'s primary image structures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Added new repo: scijava.public\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "method": "display_data"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9187a75e-394c-4293-aafc-76a53df0d192",
       "version_major": 2,
       "version_minor": 0
      },
      "method": "display_data"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "ImageJ v2.0.0-rc-71 is ready to go."
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%classpath config resolver scijava.public https://maven.scijava.org/content/groups/public\n",
    "%classpath add mvn net.imagej imagej 2.0.0-rc-71\n",
    "ij = new net.imagej.ImageJ()\n",
    "\"ImageJ v${ij.getVersion()} is ready to go.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Images in ImageJ are backed by data structures from the [ImgLib2](https://imagej.net/ImgLib2) library. There are several important interfaces, but the two most crucial to understand initially are `IterableInterval` and `RandomAccessibleInterval`.\n",
    "\n",
    "As you can infer from their names, both of them are `Interval`s, which means finite, discrete point samplings in\n",
    "$\\mathbb{Z}^{n}$\n",
    "bounded in each dimension. For example, a `uint8` image with dimensions $1024 \\times 768$ is an interval in $\\mathbb{Z}^{2}$ bounded in X (the first dimension) by $[0, 1023]$ and bounded in Y (the second dimension) by $[0, 767]$, with each integer coordinate inside the interval possessing some value $v \\in \\mathbb{Z} : 0 \\leq v \\leq 255$.\n",
    "\n",
    "* An `IterableInterval` is an image which can be [_iterated_](https://en.wikipedia.org/wiki/Iteration#Computing) in some order. That is, you can loop over its samples, although the _iteration order_ may vary depending on the type of image—that is, you cannot rely on a particular sample order _a priori_. The good news is that an `IterableInterval` _does_ know its dimensional position at each iteration. An `IterableInterval` is essentially a [stream](https://en.wikipedia.org/wiki/Stream_(computing)).\n",
    "\n",
    "* A `RandomAccessibleInterval` is an image which can be _inspected_ at will at arbitrary positions within the interval. In other words, it provides [_random access_](https://en.wikipedia.org/wiki/Random_access) to the image samples.\n",
    "\n",
    "If you are familiar with the `java.io` package of the Java standard library, `IterableInterval` is to `java.io.InputStream` as `RandomAccessibleInterval` is to `java.io.RandomAccessFile`.\n",
    "\n",
    "Let's dive into the API for each of these types of images!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Iterating images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a demonstration of image iteration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The image is a net.imglib2.img.array.ArrayImg\n",
      "image instanceof IterableInterval? true\n",
      "Each sample is a net.imglib2.type.numeric.real.DoubleType\n",
      "Sample values = 0.0 1.0 2.0 3.0 4.0 1.0 2.0 3.0 4.0 5.0 2.0 3.0 4.0 5.0 6.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "null"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import net.imglib2.IterableInterval\n",
    "\n",
    "// Create a tiny image.\n",
    "image = ij.op().run(\"create.img\", [5, 3])\n",
    "println(\"The image is a \" + image.getClass().getName())\n",
    "println(\"image instanceof IterableInterval? \" + (image instanceof IterableInterval))\n",
    "println(\"Each sample is a \" + image.firstElement().getClass().getName())\n",
    "\n",
    "// Populate it with a diagonal gradient.\n",
    "ij.op().image().equation(image, \"p[0]+p[1]\")\n",
    "\n",
    "// Iterate over the image samples!\n",
    "print(\"Sample values =\")\n",
    "// A foreach loop in Groovy uses \"in\" to seperate arguments whereas a Java foreach loop uses \":\" to seperate arguments\n",
    "for (v in image)\n",
    "  print(\" \" + v)\n",
    "println()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The actual type of the created image object is `ArrayImg`, which is an _image container_ backed by one big array—in this case, a `double[]` because the samples are of type `DoubleType`.\n",
    "\n",
    "More generally: an `Img` is an object which is both an `IterableInterval` and a `RandomAccessibleInterval`. Other kinds of image containers include `PlanarImg` (one array per 2D slice) and `CellImg` (one array per N-dimensional block).\n",
    "\n",
    "You might be wondering: what if I need to know the dimensional position during iteration? The solution is to use a _cursor_:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Samples by position:\n",
      "\t(0, 0) = 0.0\t(1, 0) = 1.0\t(2, 0) = 2.0\t(3, 0) = 3.0\t(4, 0) = 4.0\n",
      "\t(0, 1) = 1.0\t(1, 1) = 2.0\t(2, 1) = 3.0\t(3, 1) = 4.0\t(4, 1) = 5.0\n",
      "\t(0, 2) = 2.0\t(1, 2) = 3.0\t(2, 2) = 4.0\t(3, 2) = 5.0\t(4, 2) = 6.0"
     ]
    },
    {
     "data": {
      "text/plain": [
       "null"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(\"Samples by position:\")\n",
    "c = image.localizingCursor()\n",
    "while (c.hasNext()) {\n",
    "  v = c.next()\n",
    "  xPos = c.getLongPosition(0)\n",
    "  yPos = c.getLongPosition(1)\n",
    "  if (xPos == 0) println()\n",
    "  print(\"\\t(\" + xPos + \", \" + yPos + \") = \" + v)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each cursor is a pointer into the image somewhere, which knows its position. It is possible to run multiple cursors over an image simultaneously.\n",
    "\n",
    "Note that for performance, we use a _localizing cursor_ above by calling `localizingCursor()`, because we knew we would query the position every time. If we were going to query the position only rarely, `cursor()` would be better—a vanilla cursor still knows its position, but does less bookkeeping and hence is faster to iterate in cases where you don't query very often."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Iteration orders"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the positions above, we see that the iteration order here happens to be nice and organized, with X moving fastest, followed by Y. This so-called _flat iteration order_ is what you get with `ArrayImg`, but may differ with other image containers. Here is an example of how a _cell image_ with $2 \\times 2 \\times 2$ blocks differs from an _array image_ of the same size:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Array image iteration order =\n",
      "\t(0, 0, 0)\t(1, 0, 0)\t(2, 0, 0)\t(3, 0, 0)\n",
      "\t(0, 1, 0)\t(1, 1, 0)\t(2, 1, 0)\t(3, 1, 0)\n",
      "\t(0, 0, 1)\t(1, 0, 1)\t(2, 0, 1)\t(3, 0, 1)\n",
      "\t(0, 1, 1)\t(1, 1, 1)\t(2, 1, 1)\t(3, 1, 1)\n",
      "\n",
      "Cell image iteration order =\n",
      "\t(0, 0, 0)\t(1, 0, 0)\t(0, 1, 0)\t(1, 1, 0)\n",
      "\t(0, 0, 1)\t(1, 0, 1)\t(0, 1, 1)\t(1, 1, 1)\n",
      "\t(2, 0, 0)\t(3, 0, 0)\t(2, 1, 0)\t(3, 1, 0)\n",
      "\t(2, 0, 1)\t(3, 0, 1)\t(2, 1, 1)\t(3, 1, 1)\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "null"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "// Create an array image and a cell (2x2x2) image.\n",
    "import net.imglib2.img.array.ArrayImgs\n",
    "import net.imglib2.img.cell.CellImgFactory\n",
    "long[] dims = [4, 2, 2]\n",
    "arrayImg = ArrayImgs.unsignedBytes(dims)\n",
    "cellImg = new CellImgFactory(2).create(dims, arrayImg.firstElement())\n",
    "\n",
    "def printPositions(ii, width) {\n",
    "  c = ii.localizingCursor()\n",
    "  col = 0\n",
    "  while (c.hasNext()) {\n",
    "    v = c.next()\n",
    "    xPos = c.getLongPosition(0)\n",
    "    yPos = c.getLongPosition(1)\n",
    "    zPos = c.getLongPosition(2)\n",
    "    print(\"\\t(\" + xPos + \", \" + yPos + \", \" + zPos + \")\")\n",
    "    if (++col == width) { col = 0;  println() }\n",
    "  }\n",
    "  println()\n",
    "}\n",
    "\n",
    "println()\n",
    "println(\"Array image iteration order =\")\n",
    "printPositions(arrayImg, 4)\n",
    "\n",
    "println(\"Cell image iteration order =\")\n",
    "printPositions(cellImg, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Access Images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accessing image samples directly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes you need to access samples in an order other than the natural one. This is possible as long as the image implements the `RandomAccessible` interface.\n",
    "\n",
    "By using a `RandomAccess` accessor, we can get access to any pixel of the image in any order:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Array image instanceof RandomAccessible? true\n",
      "Cell image instanceof RandomAccessible? true\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class =\"scijava\"><tbody><tr><td class =\"rowLabel\">ArrayImg</td><td><img src=\"data:image/png;charset=utf-8;base64,iVBORw0KGgoAAAANSUhEUgAAAIAAAABgCAYAAADVenpJAAAA30lEQVR42u3csQkAIAxFwb//0rFxgwiiucBrLcxhaZKURucSABAAAkAACAABIAAEgAAQAAJAAAgA/VmZ0dN6AfYButiBHQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAJQAgAASAABAAAkAAaAQA+SBCAAgAASAABIAAEAACQAAIAAEgAASAABAAAkAACAABIAAEgAAQAAJAAAgAASAABIAAEAACQAD06wwAAkAACAABIAAEgAAQAAJAAAgAASAABIAebAFWCTGOgmqCsgAAAABJRU5ErkJggg==\" /></td></tr><tr><td class =\"rowLabel\">CellImg</td><td><img src=\"data:image/png;charset=utf-8;base64,iVBORw0KGgoAAAANSUhEUgAAAIAAAABgCAYAAADVenpJAAAA30lEQVR42u3csQkAIAxFwb//0rFxgwiiucBrLcxhaZKURucSABAAAkAACAABIAAEgAAQAAJAAAgA/VmZ0dN6AfYButiBHQAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAJQAgAASAABAAAkAAaAQA+SBCAAgAASAABIAAEAACQAAIAAEgAASAABAAAkAACAABIAAEgAAQAAJAAAgAASAABIAAEAACQAD06wwAAkAACAABIAAEgAAQAAJAAAgAASAABIAebAFWCTGOgmqCsgAAAABJRU5ErkJggg==\" /></td></tr></tbody></table>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import net.imglib2.RandomAccessibleInterval\n",
    "import net.imglib2.img.array.ArrayImgs\n",
    "import net.imglib2.img.cell.CellImgFactory\n",
    "import net.imglib2.type.numeric.integer.UnsignedByteType;\n",
    "import net.imglib2.type.numeric.real.DoubleType;\n",
    "\n",
    "// create array image and cell image\n",
    "long[] dims = [128,96]\n",
    "arrayImg = ArrayImgs.unsignedBytes(dims);\n",
    "cellImg = new CellImgFactory(1).create(dims, arrayImg.firstElement())\n",
    "\n",
    "// is the image RandomAccessible\n",
    "println(\"Array image instanceof RandomAccessible? \" + (arrayImg instanceof RandomAccessibleInterval))\n",
    "println(\"Cell image instanceof RandomAccessible? \" + (cellImg instanceof RandomAccessibleInterval))\n",
    "\n",
    "// start drawing\n",
    "def drawimage(image){\n",
    "    ra = image.randomAccess()\n",
    "    ra.setPosition(15,0)\n",
    "    for (y in 15..45) {\n",
    "        ra.setPosition(y, 1)\n",
    "        UnsignedByteType t = ra.get()\n",
    "        t.set(255)\n",
    "    }\n",
    "\n",
    "    ra.setPosition(15,1)\n",
    "    for (x in 0..cellImg.dimension(0) - 1){\n",
    "        ra.setPosition(x, 0)\n",
    "        UnsignedByteType t =  ra.get()\n",
    "        t.set(255)\n",
    "    }\n",
    "\n",
    "    ra.setPosition(112, 0)\n",
    "    for (y in 15..45) {\n",
    "        ra.setPosition(y, 1)\n",
    "        UnsignedByteType t = ra.get()\n",
    "        t.set(255)\n",
    "    }\n",
    "\n",
    "    ra.setPosition(80, 1)\n",
    "    for (x in 57..71){\n",
    "        ra.setPosition(x, 0)\n",
    "        UnsignedByteType t =  ra.get()\n",
    "        t.set(255)\n",
    "    }\n",
    "}\n",
    "\n",
    "drawimage(arrayImg)\n",
    "drawimage(cellImg)\n",
    "ij.notebook().display(['ArrayImg': arrayImg,'CellImg': cellImg])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TODO:\n",
    "* How to perform image processing operations such as convolution, deconvolution and Fourier transforms.\n",
    "* How to work with non-rectangular regions of interest (ROIs) and do computational geometry.\n",
    "* How to compute image features and use them with machine learning algorithms."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Groovy",
   "language": "groovy",
   "name": "groovy"
  },
  "language_info": {
   "codemirror_mode": "groovy",
   "file_extension": ".groovy",
   "mimetype": "",
   "name": "Groovy",
   "nbconverter_exporter": "",
   "version": "2.5.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}