{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Building RDME geometry from scratch: jLM\n",
    "\n",
    "Tyler Earnest\n",
    "\n",
    "\"Hands-On\" Workshop on Cell Scale simulations at Urbana, IL\n",
    "\n",
    "July 20, 2018"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# jLM Philosophy\n",
    "\n",
    "Work with binary masks:\n",
    "* 3D data type with same shape as RDME lattice\n",
    "* Boolean: each element is `True` or `False`\n",
    "* True means that site is in the region\n",
    "\n",
    "Site type lattice is constructed from multiple disjoint binary masks\n",
    "\n",
    "Using binary masks allows using\n",
    "* Set operations to build up regions\n",
    "* Morphological image processing to sculpt regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.simplefilter(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Standard imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.ndimage as spnd\n",
    "import ipywidgets as ipw\n",
    "import h5py\n",
    "import itertools\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from jLM.RegionBuilder import RegionBuilder\n",
    "from jLM.RDME import Sim as RDMESim\n",
    "import jLM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import seaborn\n",
    "seaborn.set_context(\"talk\", font_scale=1.4)\n",
    "seaborn.set_style(\"darkgrid\")\n",
    "warnings.simplefilter(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Create the simulation object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "32e-9*32"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "sim = RDMESim(\"E. coli\",\n",
    "              'tmp.lm',\n",
    "              [32,32,128],\n",
    "              32e-9,\n",
    "              \"extracellular\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Define the regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "ext = sim.region(\"extracellular\")\n",
    "mem = sim.region(\"membrane\")\n",
    "cyt = sim.region(\"cytoplasm\")\n",
    "nuc = sim.region(\"nucleoid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Create a `RegionBuilder`\n",
    "\n",
    "The region builder is a helper object which is used to construct site lattice geometry\n",
    "\n",
    "Instead of defining the geometry first then discretizing it to a lattice, we work in the discretized representation from the start."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "B = RegionBuilder(sim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Define the dimensions and capsule centers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "cell_length = 3.6e-6/sim.latticeSpacing\n",
    "cell_width = 0.8e-6/sim.latticeSpacing\n",
    "nucleoid_width = cell_width/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "cylinderHalfLength = 0.5*(cell_length-cell_width)\n",
    "p1 = B.center - [0,0,cylinderHalfLength]\n",
    "p2 = B.center + [0,0,cylinderHalfLength]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Create the capsule\n",
    "\n",
    "This region will be carved into the cytoplasm\n",
    "\n",
    "Instead of stacking layers, we sculpt the geometry and combine the non-overlapping regions to form the site lattice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "sphere1 = B.ellipsoid(radius=cell_width/2, center=p1)\n",
    "sphere2 = B.ellipsoid(radius=cell_width/2, center=p2)\n",
    "cylinder = B.cylinder(radius=cell_width/2, length=cell_length-cell_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Regions are combined using set operations\n",
    "\n",
    "* Union: `|`\n",
    "* Intersection: `&`\n",
    "* Complement `~`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "capsule = sphere1 | sphere2 | cylinder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Inspect the capsule\n",
    "Need to make sure we got the positioning and dimensions correct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "B.showStack(capsule, scl=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Nucleoid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "sphere1 = B.ellipsoid(radius=nucleoid_width/2, center=p1)\n",
    "sphere2 = B.ellipsoid(radius=nucleoid_width/2, center=p2)\n",
    "cylinder = B.cylinder(radius=nucleoid_width/2, length=cell_length-cell_width)\n",
    "nucleoid = sphere1 | sphere2 | cylinder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Cytoplasm\n",
    "\n",
    "The cytoplasm is outside the nucleoid. We take the initial capsule and subtract the nucleoid region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "cytoplasm = capsule & ~nucleoid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Membrane\n",
    "\n",
    "The membrane is created by growing the capsule by a lattice site, then subtracting the initial capsule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "How to modify the lattice in  this way?\n",
    "* Morphological image processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Binary morphology\n",
    "\n",
    "Basic idea\n",
    "* Probe an image with a pre-defined shape and generate a new image based on how the shape interacts with the image\n",
    "* The probe shape is called the *structuring element*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Binary morphology\n",
    "\n",
    "Standard Morphological operations\n",
    "* Erosion\n",
    "* Dilation\n",
    "* Opening\n",
    "* Closing\n",
    "\n",
    "These can all be described as using set notation, but it's easier to explain with examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Structuring element\n",
    "\n",
    "Let's define a structuring element with nearest neighbor connectivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "se = spnd.generate_binary_structure(2,1)\n",
    "im = np.pad(se,(4,4), \"constant\")\n",
    "plt.imshow(im, interpolation='nearest')\n",
    "plt.scatter([5],[5])\n",
    "_=plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Example binary image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "img = np.zeros((64,64), dtype=bool)\n",
    "img[10:30,10:30] = True\n",
    "img[50:53,10:13] = True\n",
    "img[50,50] = True\n",
    "plt.imshow(img, interpolation='nearest')\n",
    "_=plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Dilation\n",
    "\n",
    "The morphological operation \"dilation\" is written as\n",
    "\n",
    "$$ I_{\\mathrm{dil}} = I_0\\oplus S $$\n",
    "\n",
    "for an image $I$ and a structuring element $S$.\n",
    "\n",
    "Essentially\n",
    "* $S$ is placed over each pixel in $I$\n",
    "* If any of the pixels in $I$ overlap with $S$, then $I_{\\mathrm{dil}}$ is set to True at that position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "dil = spnd.binary_dilation(img, se)\n",
    "\n",
    "fig,axs=plt.subplots(ncols=2,figsize=(10,5))\n",
    "axs[0].imshow(img, interpolation='nearest')\n",
    "axs[0].set_title(\"Input image\")\n",
    "axs[0].axis(\"off\")\n",
    "axs[1].imshow(dil, interpolation='nearest')\n",
    "axs[1].set_title(\"Dilated image\")\n",
    "axs[1].axis(\"off\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plt.imshow(1.0*img+dil)\n",
    "plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Try a different SE\n",
    "\n",
    "Second nearest neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "se = spnd.generate_binary_structure(2,2)\n",
    "im = np.pad(se,(4,4), \"constant\")\n",
    "plt.imshow(im, interpolation='nearest')\n",
    "plt.scatter([5],[5])\n",
    "_=plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dil = spnd.binary_dilation(img, se)\n",
    "\n",
    "fig,axs=plt.subplots(ncols=2,figsize=(10,5))\n",
    "axs[0].imshow(img, interpolation='nearest')\n",
    "axs[0].set_title(\"Input image\")\n",
    "axs[0].axis(\"off\")\n",
    "axs[1].imshow(dil, interpolation='nearest')\n",
    "axs[1].set_title(\"Dilated image\")\n",
    "axs[1].axis(\"off\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plt.imshow(1.0*img+dil)\n",
    "plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Erosion\n",
    "\n",
    "Erosion is like the reverse of dilation:\n",
    "\n",
    "$$ I_{\\mathrm{ero}} = I_0\\ominus S $$\n",
    "\n",
    "\n",
    "* $S$ is placed over each pixel in $I$\n",
    "* If any of the pixels in $S$ overlap with the background, then $I_{\\mathrm{ero}}$ is set to False at that position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "se = spnd.generate_binary_structure(2,1)\n",
    "dil = spnd.binary_erosion(img, se)\n",
    "\n",
    "fig,axs=plt.subplots(ncols=2,figsize=(10,5))\n",
    "axs[0].imshow(img, interpolation='nearest')\n",
    "axs[0].set_title(\"Input image\")\n",
    "axs[0].axis(\"off\")\n",
    "axs[1].imshow(dil, interpolation='nearest')\n",
    "axs[1].set_title(\"Eroded image\")\n",
    "axs[1].axis(\"off\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Other operations\n",
    "\n",
    "Less useful for generating geometry, but highly useful for working with data\n",
    "\n",
    "Opening\n",
    "$$ I\\circ S = (I\\ominus S)\\oplus S $$\n",
    "* Erosion followed by dilation\n",
    "* Useful for removing small objects/noise\n",
    "\n",
    "Closing\n",
    "$$ I\\bullet S = (I\\oplus S)\\ominus S $$\n",
    "* Dilation followed by erosion\n",
    "* Useful for filling holes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Constructing boundaries\n",
    "\n",
    "Inside\n",
    "$$I_{in} = I \\setminus (I\\ominus S) $$\n",
    "or outside\n",
    "$$I_{out} = (I\\oplus S) \\setminus I $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "se = spnd.generate_binary_structure(2,2)\n",
    "inner = img & ~spnd.binary_erosion(img,se)\n",
    "outer = spnd.binary_dilation(img,se) & ~img\n",
    "\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(10,5))\n",
    "ax[0].imshow(inner, interpolation='nearest')\n",
    "ax[0].set_title(\"inner\")\n",
    "ax[1].imshow(outer, interpolation='nearest')\n",
    "ax[1].set_title(\"outer\")\n",
    "ax[0].axis(\"off\")\n",
    "ax[1].axis(\"off\")\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Back to membrane\n",
    "\n",
    "First we grow the capsule by applying a dilation with a structuring element which includes the first and second nearest neighbors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dil_capsule = B.dilate(capsule, se=B.se26)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Then we subtract the original capsule from its dilation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "membrane = dil_capsule & ~capsule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Check that the region is contiguous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labeled, n = spnd.label(membrane)\n",
    "print(\"disjoint regions\", n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Extracellular\n",
    "\n",
    "Since the dilated capsule contains the membrane, cytoplasm, and nucleoid, the extracellular region is just its complement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "extracellular = ~dil_capsule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Compose the site lattice\n",
    "\n",
    "Now that we have all regions defined, we create the site type lattice by calling `B.compose()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "B.compose(\n",
    "    (ext, extracellular),\n",
    "    (cyt, cytoplasm),\n",
    "    (nuc, nucleoid),\n",
    "    (mem, membrane) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Check our work"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "sim.showRegionStack(scl=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "sim.displayGeometry() "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}