image_analysis_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Analysis with Python - <font color='teal'>Tutorial Pipeline</font>\n",
    "\n",
    "*originally created in 2016*<br>\n",
    "*updated and converted to a Jupyter notebook in 2017*<br>\n",
    "*updated and converted to python 3 in 2018*<br>\n",
    "*by Jonas Hartmann (Gilmour group, EMBL Heidelberg)*<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Table of Contents\n",
    "\n",
    "1. [About this Tutorial](#about)\n",
    "2. [Setup](#setup)\n",
    "3. [Importing Modules & Packages](#import)\n",
    "4. [Loading & Handling Image Data](#load)\n",
    "5. [Preprocessing](#prepro)\n",
    "6. [Manual Thresholding & Threshold Detection](#thresh)\n",
    "7. [Adaptive Thresholding](#adaptive)\n",
    "8. [Improving Masks with Binary Morphology](#morpho)\n",
    "9. [Connected Components Labeling](#label)\n",
    "10. [Cell Segmentation by Seeding & Expansion](#seg)\n",
    "11. [Postprocessing: Removing Cells at the Image Border](#postpro)\n",
    "12. [Identifying Cell Edges](#edges)\n",
    "13. [Extracting Quantitative Measurements](#measure)\n",
    "14. [Simple Analysis & Visualization](#analysis)\n",
    "15. [Writing Output to Files](#write)\n",
    "16. [Batch Processing](#batch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  About this Tutorial <a id=about></a>\n",
    "\n",
    "*This tutorial aims to teach the basics of (bio-)image processing with python, in particular the analysis of fluorescence microscopy data. The tutorial is self-explanatory and follows a \"learning by doing\" philosophy. It consists of step by step instructions that guide students through the construction of a 2D single-cell segmentation pipeline.*\n",
    "\n",
    "\n",
    "#### Instructions\n",
    "\n",
    "- This notebook contains detailed instructions on how to program a pipeline that can segment cells in 2D fluorescence microscopy images of a developing tissue.\n",
    "\n",
    "\n",
    "- Simply go through the instructions step by step and try to implement each step as best as you can.\n",
    "    - By doing so, you will learn about key concepts of (bio-)image processing with python!\n",
    "\n",
    "\n",
    "- If you are stuck...\n",
    "    - ...first think some more about the problem and see if you can get yourself unstuck.\n",
    "    - ...search the internet (in particular Stack Overflow) for a solution to your problem.\n",
    "    - ...if you are working through this tutorial in class, ask one of the tutors for help.\n",
    "    - ...if nothing else helps, you can have a look at the solutions (`image_analysis_tutorial_solutions.py`) for inspiration.\n",
    "        - **Beware:** If you simply copy a solution, you will learn *nothing* in the process! Use the solutions only as an inspiration for implementing your own version and make sure you fully understand how the solution works!\n",
    "        \n",
    "        \n",
    "- Some exercises are labeled as 'BONUS' exercises. It is up to you whether you'd like to complete them. \n",
    "    - If things are going slowly, you may want to skip them and come back to them later.\n",
    "\n",
    "\n",
    "#### Background \n",
    "\n",
    "The aim is to construct a pipeline for the identification and segmentation of cells in 2D confocal fluorescence microscopy images of a tissue with labeled membranes. This is among the most common tasks in bio-image analysis and is often essential for the extraction of useful quantitative information from microscopy data.\n",
    "\n",
    "The pipeline is constructed based on provided example images (see `example_data` directory), which are single-color spinning-disk confocal micrographs (objective: `40X 1.2NA W`) of cells in live zebrafish embryos in early development (~10h post fertilization), fluorescently labeled with a membrane-localized fusion protein (`mNeonGreen:Ggamma9`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Setup <a id=setup></a>\n",
    "\n",
    "#### Development Environment\n",
    "\n",
    "- Depending on your preferences, there are several alternatives for implementing this pipeline:\n",
    "    1. Implement it directly in this Jupyter notebook *[recommended]*.\n",
    "    2. Implement it in your favorite Integrated Development Environment (IDE) like PyCharm.\n",
    "    3. Old School: Implement it in a text editor (like Vim...) and run it on the terminal.\n",
    "\n",
    "\n",
    "#### Example Data\n",
    "\n",
    "Make sure that you have downloaded the example image data (`example_cells_1.tif` and `example_cells_2.tif`), which should be located in the directory `example_data`, which in turn should be in the same directory as this notebook.\n",
    "\n",
    "\n",
    "#### Python version\n",
    "\n",
    "- The solutions were last tested with **python 3.9**, but any version of python 3 should work without changes\n",
    "- To install python, jupyter notebook and the required modules, we recommend the **Anaconda Distribution**. The installer can be downloaded [here](https://www.anaconda.com/download/).\n",
    "\n",
    "\n",
    "#### Required modules\n",
    "\n",
    "- Make sure the following modules are installed before you get started:\n",
    "    - numpy\n",
    "    - scipy\n",
    "    - matplotlib\n",
    "    - scikit-image\n",
    "    - ipywidgets\n",
    "- All required modules come pre-installed if you are using the **[Anaconda distribution](https://www.anaconda.com/download/)** of python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Everything ready? Yes? Then let's get started!*\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing Modules & Packages <a id=import></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by importing the package NumPy, which enables the manipulation of numerical arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*<font color=red>Important note:</font> If you are not at all familiar with arrays and NumPy, we strongly recommend that you first complete an introductory tutorial on this topic before carrying on!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that, once imported, we can use functions/modules from the package, for example to create an array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.array([1, 2, 3])\n",
    "\n",
    "print(a)\n",
    "print(type(a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the package is imported under a variable name (here `np`). You can freely choose this name yourself. For example, it would be just as valid (but not as convenient) to write:\n",
    "\n",
    "```python\n",
    "import numpy as lovelyArrayTool\n",
    "a = lovelyArrayTool.array([1,2,3])\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "Using the import command as above, follow the instructions in the comments below to import two additional modules that we will be using frequently in this pipeline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The plotting module matplotlib.pyplot as plt\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# The image processing module scipy.ndimage as ndi\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Side Note for Jupyter Notebook Users\n",
    "\n",
    "You can configure how the figures made by matplotlib are displayed.\n",
    "\n",
    "The most common options are the following:\n",
    "\n",
    "- **inline**:   displays as static figure in code cell output\n",
    "- **notebook**: displays as interactive figure in code cell output \n",
    "- **qt**:       displays as interactive figure in a separate window\n",
    "\n",
    "Feel free to test them out on one of the figures you will generate later on in the tutorial. The code cell below shows how to set the different options. Note that combinations of different options in the same notebook do not always work well, so it is best to decide for one and use it throughout. You may need to restart the kernel (`Kernel > Restart`) when you change from one option to another."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set matplotlib backend\n",
    "%matplotlib inline \n",
    "#%matplotlib notebook \n",
    "#%matplotlib qt "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading & Handling Image Data <a id=load></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "Images are essentially just numbers (representing intensity) in an ordered grid of pixels. Image processing is simply to carry out mathematical operations on these numbers.\n",
    "\n",
    "The ideal object for storing and manipulating ordered grids of numbers is the **array**. Many mathematical operations are well defined on arrays and can be computed quickly by vector-based computation.\n",
    "\n",
    "Arrays can have any number of dimensions (or \"axes\"). For example, a 2D array could represent the x and y axis of a grayscale image (xy), a 3D array could contain a z-stack (zyx), a 4D array could also have multiple channels for each image (czyx) and a 5D array could have time on top of that (tczyx)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "We will now proceed to load one of the example images and verify that we get what we expect. \n",
    "\n",
    "Note: Before starting, it always makes sense to have a quick look at the data in Fiji/ImageJ so you know what you are working with!\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Specify the directory path and file name\n",
    "\n",
    "# Create a string variable with the name of the file you'd like to load (here: 'example_cells_1.tif').\n",
    "# Suggested name for the variable: filename\n",
    "# Note: Paths and filenames can contain slashes, empty spaces and other special symbols, which can cause \n",
    "#       trouble for programming languages under certain circumstances. To circumvent such trouble, add \n",
    "#       the letter r before your string definition to create a so-called 'raw string', which is not\n",
    "#       affected by these problems (e.g. `my_raw_string = r\"some string with funny symbols: \\\\\\!/~***!\"`).\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# If the file is not in the current working directory, you must also have a way of specifying the path\n",
    "# to the directory where the file is stored. Most likely, your example images are stored in a directory\n",
    "# called 'example_data' in the same folder as this notebook. Note that you can use either the full path\n",
    "#  - something like r\"/home/jack/data/python_image_analysis/example_data\"\n",
    "# or the relative path, starting from the current working directory\n",
    "#   - here that would just be r\"example_data\"\n",
    "# Create a string variable with the path to the directory that contains the file you'd like to load.\n",
    "# Suggested name for the variable: dirpath\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Combine the directory path and file name into one variable, the file path\n",
    "\n",
    "# Import the function 'join' from the module 'os.path'\n",
    "# This function automatically takes care of the slashes that need to be added when combining two paths.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Use the 'join' function to combine the directory path with the file name and create a new variable.\n",
    "# Print the result to see that everything is correct (this is always a good idea!)\n",
    "# Suggested name for the variable: filepath\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Load the image\n",
    "\n",
    "# Import the function 'imread' from the module 'skimage.io'.\n",
    "# (Note: If this gives you an error, please refer to the note below!)\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Load 'example_cells_1.tif' and store it in a variable.\n",
    "# Suggested name for the variable: img\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "*Important note for those who get an error when trying to import `imread` from `skimage.io`:*\n",
    "\n",
    "Some users have been experiencing problems with this module, even though the rest of skimage is installed correctly (running `import skimage` does not given an error). This may have something to do with operating system preferences. The easiest solution in this case is to install the module `tifffile` (with three`f`) and use the function `imread` from that module (it is identical to the `imread` function of `skimage.io` when reading `tif` files). \n",
    "\n",
    "The `tifffile` module does not come with the Anaconda distribution, so it's likely that you don't have it installed. To install it, save and exit Jupyter notebook, then go to a terminal and type `conda install -c conda-forge tifffile`. After the installation is complete, restart Jupyter notebook, come back here and import `imread` from `tifffile`. This should now hopefully work.\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iv) Check that everything is in order\n",
    "\n",
    "# Check that 'img' is a variable of type 'ndarray' - use Python's built-in function 'type'.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Print the shape of the array by looking at its 'shape' attribute. \n",
    "# Make sure you understand the output!\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Check the datatype of the individual numbers in the array. You can use the array attribute 'dtype' to do so.\n",
    "# Make sure you understand the output!\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (v) Look at the image to confirm that everything worked as intended\n",
    "\n",
    "# To plot an array as an image, use pyplot's functions 'plt.imshow' followed by 'plt.show'. \n",
    "# Check the documentation for 'plt.imshow' and note the parameters that can be specified, such as colormap (cmap)\n",
    "# and interpolation. Since you are working with scientific data, interpolation is unwelcome, so you should set it \n",
    "# to \"none\". The most common cmap for grayscale images is naturally \"gray\".\n",
    "# You may also want to adjust the size of the figure. You can do this by preparing the figure canvas with\n",
    "# the function 'plt.figure' before calling 'plt.imshow'. The canvas size is adjusted using the keyword argument\n",
    "# 'figsize' when calling 'plt.figure'.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing <a id=prepro></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "The goal of image preprocessing is to prepare or optimize the images to make further analysis easier. Usually, this boils down to increasing the signal-to-noise ratio by removing noise and background and by enhancing structures of interest.\n",
    "\n",
    "The specific preprocessing steps used in a pipeline depend on the type of sample, the microscopy technique used, the image quality, and the desired downstream analysis. \n",
    "\n",
    "The most common operations include:\n",
    "\n",
    "\n",
    "- Deconvolution\n",
    "    - Image reconstruction based on information about the PSF of the microscope\n",
    "    - These days deconvolution is often included with microscope software\n",
    "    - *Our example images are not deconvolved, but will do just fine regardless*\n",
    "\n",
    "\n",
    "- Conversion to 8-bit images to save memory / computational time\n",
    "    - *Our example images are already 8-bit*\n",
    "\n",
    "\n",
    "- Cropping of images to an interesting region\n",
    "    - *The field of view in our example images is fine as it is*\n",
    "\n",
    "\n",
    "- Smoothing of technical noise\n",
    "    - This is a very common step and usually helps to improve almost any type of downstream analysis\n",
    "    - Commonly used filters are the `Gaussian filter` and the `median filter`\n",
    "    - *Here we will be using a Gaussian filter.*\n",
    "\n",
    "\n",
    "- Corrections of technical artifacts\n",
    "    - Common examples are uneven illumination and multi-channel bleed-through\n",
    "    - *Here we will deal with uneven signal by adaptive/local thresholding*\n",
    "\n",
    "\n",
    "- Background subtraction\n",
    "    - There are various ways of sutracting background signal from an image\n",
    "    - Two different types are commonly distinguished:\n",
    "        - `uniform background subtraction` treats all regions of the image the same\n",
    "        - `adaptive or local background subtraction` automatically accounts for differences between regions of the image\n",
    "    - *Here we will do something similar to adaptive background subtraction when we do adaptive thresholding*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gaussian Smoothing\n",
    "\n",
    "A Gaussian filter smoothens an image by convolving it with a Gaussian-shaped kernel. In the case of a 2D image, the Gaussian kernel is also 2D and will look something like this:\n",
    "\n",
    "<img src=\"ipynb_images\\gaussian_kernel_grid.png\" alt=\"Gaussian Kernel Figure\" style=\"width: 300px;\"/>\n",
    "\n",
    "How much the image is smoothed by a Gaussian kernel is determined by the standard deviation  of the Gaussian distribution, usually referred to as **sigma** ($\\sigma$). A higher $\\sigma$ means a broader distribution and thus more smoothing.\n",
    "\n",
    "**How to choose the correct value of $\\sigma$?**\n",
    "\n",
    "This depends a lot on your images, in particular on the pixel size. In general, the chosen $\\sigma$ should be large enough to blur out noise but small enough so the \"structures of interest\" do not get blurred too much. Usually, the best value for $\\sigma$ is simply found by trying out some different options and looking at the result. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "Perform Gaussian smoothing and visualize the result.\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Create a variable for the smoothing factor sigma, which should be an integer value\n",
    "\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# After implementing the Gaussian smoothing function below, you can modify this variable \n",
    "# to find the ideal value of sigma."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Perform the smoothing on the image\n",
    "\n",
    "# To do so, use the Gaussian filter function 'ndi.gaussian_filter' from the \n",
    "# image processing module 'scipy.ndimage', which was imported at the start of the tutorial. \n",
    "# Check out the documentation of scipy to see how to use this function. \n",
    "# Allocate the output to a new variable.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Visualize the result using 'plt.imshow'\n",
    "\n",
    "# Compare with the original image visualized above. \n",
    "# Does the output make sense? Is this what you expected? \n",
    "# Can you optimize sigma such that the image looks smooth without blurring the membranes too much?\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# To have a closer look at a specific region of the image, crop that region out and show it in a \n",
    "# separate plot. Remember that you can crop arrays by \"indexing\" or \"slicing\" them similar to lists.\n",
    "# Use such \"zoomed-in\" views throughout this tutorial to take a closer look at your intermediate \n",
    "# results when necessary.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iv) BONUS: Show the raw and smoothed images side by side using 'plt.subplots'\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manual Thresholding & Threshold Detection <a id=thresh></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "The easiest way to distinguish foreground objects (here: membranes) from the image background is to threshold the image, meaning all pixels with an intensity above a certain threshold are accepted as foreground, all others are set as background.\n",
    "\n",
    "To find the best threshold for a given image, one option is to simply try out different thresholds manually. Alternatively, one of many algorithms for automated 'threshold detection' can be used. These algorithms use information about the image (such as the histogram) to automatically find a suitable threshold value, often under the assumption that the background and foreground pixels in an image belong to two clearly distinct populations in terms of their intensity. \n",
    "\n",
    "There are many different algorithms for threshold detection and it is often hard to predict which one will produce the nicest and most robust result for a particular dataset. It therefore makes sense to try out a bunch of different options.\n",
    "\n",
    "For this pipeline, we will ultimately use a more advanced thresholding approach, which also accounts (to some extent) for variations in signal across the field of view: adaptive thresholding. \n",
    "\n",
    "But first, let's experiment a bit with threshold detection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "Try out manual thresholding and automated threshold detection.\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Create a variable for a manually set threshold, which should be an integer\n",
    "\n",
    "# This can be changed later to find a suitable value.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Perform thresholding on the smoothed image\n",
    "\n",
    "# Remember that you can use relational (Boolean) expressions such as 'smaller' (<), 'equal' (==)\n",
    "# or 'greater or equal' (>=) with numpy arrays - and you can directly assign the result to a new\n",
    "# variable.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Check the dtype of your thresholded image\n",
    "# You should see that the dtype is 'bool', which stands for 'Boolean' and means the array\n",
    "# is now simply filled with 'True' and 'False', where 'True' is the foreground (the regions\n",
    "# above the threshold) and 'False' is the background.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Visualize the result\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iv) Try out different thresholds to find the best one\n",
    "\n",
    "# If you are using jupyter notebook, you can adapt the code below to\n",
    "# interactively change the threshold and look for the best one. These\n",
    "# kinds of interactive functions are called 'widgets' and are very \n",
    "# useful in exploratory data analysis to create greatly simplified\n",
    "# 'User Interfaces' (UIs) on the fly.\n",
    "# As a BONUS exercise, try to understand or look up how the widget works\n",
    "# and play around with it a bit!\n",
    "# (Note: If this just displays a static image without a slider to adjust\n",
    "#        the threshold or if it displays a text warning about activating\n",
    "#        the 'widgetsnbextension', check out the note below!)\n",
    "\n",
    "# Prepare widget\n",
    "from ipywidgets import interact\n",
    "@interact(thresh=(10,250,10))\n",
    "def select_threshold(thresh=100):\n",
    "    \n",
    "    # Thresholding\n",
    "    ### ADAPT THIS: Change 'img_smooth' into the variable you stored the smoothed image in!\n",
    "    mem = img_smooth > thresh\n",
    "    \n",
    "    # Visualization\n",
    "    plt.figure(figsize=(7,7))\n",
    "    plt.imshow(mem, interpolation='none', cmap='gray')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "\n",
    "*Important note for those who get a static image (no slider) or a text warning:*\n",
    "\n",
    "For some users, it is necessary to specifically activate the widgets plugin for Jupyter notebook. To do so, save and exit Jupyter notebook, then go to a terminal (or Anaconda prompt) and write `jupyter nbextension enable --py --sys-prefix widgetsnbextension`. After this, you should be able to restart Jupyter notebook and the widget should display correctly. \n",
    "\n",
    "If it still doesn't work, you may instead have to type `jupyter nbextension enable --py widgetsnbextension` in the terminal (or Anaconda prompt). However, note that this implies that your installation of Conda/Jupyter is not optimally configured (see [this GitHub issue](https://github.com/jupyter-widgets/ipywidgets/issues/541) for more information, although this is not something you necessarily need to worry about in the context of this course).\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (v) Perfom automated threshold detection with Otsu's method\n",
    "\n",
    "# The scikit-image module 'skimage.filters.thresholding' provides\n",
    "# several threshold detection algorithms. The most popular one \n",
    "# among them is Otsu's method. Using what you've learned so far,\n",
    "# import the 'threshold_otsu' function, use it to automatically \n",
    "# determine a threshold for the smoothed image, apply the threshold,\n",
    "# and visualize the result.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (vi) BONUS: Did you notice the 'try_all_threshold' function?\n",
    "\n",
    "# That's convenient! Use it to automatically test the threshold detection\n",
    "# functions in 'skimage.filters.thresholding'. Don't forget to adjust the\n",
    "# 'figsize' parameter so the resulting images are clearly visible.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive Thresholding <a id=adaptive></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "Simply applying a fixed intensity threshold does not always produce a foreground mask of sufficiently high quality, since background and foreground intensities often vary across the image. In our example image, for instance, the intensity drops at the image boundaries - a problem that cannot be resolved just by changing the threshold value.\n",
    "\n",
    "One way of addressing this issue is to use an *adaptive thresholding* algorithm, which adjusts the threshold locally in different regions of the image to account for varying intensities.\n",
    "\n",
    "Although `scikit-image` provides a function for adaptive thresholding (called `threshold_local`), we will here implement our own version, which is slightly different and will hopefully make the concept of adaptive thresholding very clear.\n",
    "\n",
    "Our approach to adaptive tresholding works in two steps:\n",
    "\n",
    "1. Generation of a \"background image\"\n",
    "\n",
    " This image should - across the entire image - always have higher intensities than the local background but lower intensities than the local foreground. This can be achieved by strong blurring/smoothing of the image, as illustrated in this 1D example:\n",
    "\n",
    " <img src=\"ipynb_images\\adaptive_bg_1D.png\" alt=\"Adaptive Background Figure\" style=\"width: 400px;\"/><br>\n",
    "    \n",
    "2. Thresholding of the original image with the background\n",
    "\n",
    " Instead of thresholding with a single value, every pixel in the image is thresholded with the corresponding pixel of the \"background image\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'> Exercise </font>\n",
    "\n",
    "Implement the two steps of the adaptive background subtraction:\n",
    "\n",
    "1. Use a strong \"mean filter\" (aka \"uniform filter\") to create the background image. This simply assigns each pixel the average value of its local neighborhood. Just like the Gaussian blur, this can be done by convolution, but this time using a \"uniform kernel\" like this one:\n",
    "\n",
    " <img src=\"ipynb_images\\uniform_filter_SE.png\" alt=\"Uniform Filter SE Figure\" style=\"width: 300px;\"/>\n",
    "    \n",
    " To define which pixels should be considered as the local neighborhood of a given pixel, a `structuring element` (`SE`) is    used. This is a small binary image where all pixels set to `1` will be considered as part of the neighborhood and all pixels set to `0` will not be considered. Here, we use a disc-shaped `SE`, as this reduces artifacts compared to a square `SE`.\n",
    "  \n",
    " *Side note:* A strong Gaussian blur would also work to create the background mask. For the Gaussian blur, the analogy to the `SE` is the `sigma` value, which in a way also determines the size of the local neighborhood.<br><br>\n",
    "\n",
    "2. Use the background image for thresholding. In practical terms, this works in exactly the same way as thresholding with a single value, since numpy arrays will automatically perform element-wise (pixel-by-pixel) comparisons when compared to other arrays of the same shape by a relational (Boolean) expression.\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1\n",
    "# ------\n",
    "\n",
    "# (i) Create a disk-shaped structuring element and asign it to a new variable.\n",
    "\n",
    "# Structuring elements are small binary images that indicate which pixels \n",
    "# should be considered as the 'neighborhood' of the central pixel. \n",
    "#\n",
    "# An example of a small disk-shaped SE would be this:\n",
    "#   0 0 1 0 0\n",
    "#   0 1 1 1 0\n",
    "#   1 1 1 1 1\n",
    "#   0 1 1 1 0\n",
    "#   0 0 1 0 0\n",
    "#\n",
    "# The expression below creates such structuring elements. \n",
    "# It is an elegant but complicated piece of code and at the moment it is not \n",
    "# necessary for you to understand it in detail. Use it to create structuring \n",
    "# elements of different sizes (by changing 'i') and find a way to visualize \n",
    "# the result (remember that the SE is just a small 'image').\n",
    "# \n",
    "# Try to answer the following questions: \n",
    "#   - Is the resulting SE really circular? \n",
    "#   - Could certain values of 'i' cause problems? If so, why?\n",
    "#   - What value of 'i' should be used for the SE?\n",
    "#     Note that, similar to the sigma in Gaussian smoothing, the size of the SE\n",
    "#     is first estimated based on the images and by thinking about what would \n",
    "#     make sense. Later it can be optimized by trial and error.\n",
    "\n",
    "# Create SE\n",
    "i = ???\n",
    "struct = (np.mgrid[:i,:i][0] - np.floor(i/2))**2 + (np.mgrid[:i,:i][1] - np.floor(i/2))**2 <= np.floor(i/2)**2\n",
    "\n",
    "# Visualize the result\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Create the background\n",
    "\n",
    "# Run a mean filter over the image using the disc SE and assign the output to a new variable.\n",
    "# Use the function 'skimage.filters.rank.mean'.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Visualize the resulting background image. Does what you get make sense?\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 2\n",
    "# ------\n",
    "\n",
    "# (iv) Threshold the Gaussian-smoothed original image against the background image created in step 1 \n",
    "#      using a relational expression\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (v) Visualize and understand the output. \n",
    "\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# What do you observe? \n",
    "# Are you happy with this result as a membrane segmentation? \n",
    "# Adapt the size of the circular SE to optimize the result!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improving Masks with Binary Morphology <a id=morpho></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "Morphological operations such as `erosion`, `dilation`, `closing` and `opening` are common tools used to improve masks after they are generated by thresholding. They can be used to fill small holes, remove noise, increase or decrease the size of an object, or smoothen mask outlines.\n",
    "\n",
    "Most morphological operations are once again simple kernel functions that are applied at each pixel of the image based on their neighborhood as defined by a `structuring element` (`SE`). For example, `dilation` simply assigns to the central pixel the maximum pixel value within the neighborhood; it is a maximum filter. Conversely, `erosion` is a minimum filter. Additional options emerge from combining the two: `morphological closing`, for example, is a `dilation` followed by an `erosion`. This is used to fill in gaps and holes or smoothing mask outlines without significantly changing the mask's area. Finally, there are also some more complicated morphological operations, such as `hole filling`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "Improve the membrane segmentation from above with morphological operations.\n",
    "\n",
    "Specifically, use `binary hole filling` to get rid of the speckles of foreground pixels that litter the insides of the cells. Furthermore, try different other types of morphological filtering to see how they change the image and to see if you can improve the membrane mask even more, e.g. by filling in gaps.\n",
    "\n",
    "Follow the instructions in the comments below. Visualize all intermediate results of your work and remember to \"zoom in\" to get a closer look by slicing out and then plotting a subsection of the image array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Get rid of speckles using binary hole filling\n",
    "\n",
    "# Use the function 'ndi.binary_fill_holes' for this. Be sure to check the docs to\n",
    "# understand exactly what it does. For this to work as intended, you will have to \n",
    "# invert the mask, which you can do using the function `np.logical_not` or the\n",
    "# corresponding operator '~'. Again, be sure to understand why this has to be done\n",
    "# and don't forget to revert the result back.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Try out other morphological operations to further improve the membrane mask\n",
    "\n",
    "# The various operations are available in the ndimage module, for example 'ndi.binary_closing'.\n",
    "# Play around and see how the different functions affect the mask. Can you optimize the mask, \n",
    "# for example by closing gaps?\n",
    "# Note that the default SE for these functions is a square. Feel free to create another disc-\n",
    "# shaped SE and see how that changes the outcome.\n",
    "# BONUS: If you pay close attention, you will notice that some of these operations introduce \n",
    "# artefacts at the image boundaries. Can you come up with a way of solving this? (Hint: 'np.pad')\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Visualize the final result\n",
    "\n",
    "### YOUR CODE HERE\n",
    "\n",
    "# At this point you should have a pretty neat membrane mask.\n",
    "# If you are not satisfied with the quality your membrane segmentation, you should go back \n",
    "# and fine-tune the size of the SE in the adaptive thresholding section and also optimize\n",
    "# the morphological cleaning operations.\n",
    "# Note that the quality of the membrane segmentation will have a significant impact on the \n",
    "# cell segmentation we will perform next."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connected Components Labeling <a id=label></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "Based on the membrane segmentation, we can get a preliminary segmentation of the cells in the image by considering each background region surrounded by membranes as a cell. This can already be good enough for many simple measurements.\n",
    "\n",
    "The only thing we still need to do in order to get there is to label each cell individually. Only if each separate cell has a unique number (an `ID`) assigned, values such as the mean intensity can be measured and analyzed at the single-cell level.\n",
    "\n",
    "The approach used to achieve this is called `connected components labeling`. It gives every connected group of foreground pixels a unique `ID` number."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "Use your membrane segmentation for connected components labeling.\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Label connected components\n",
    "\n",
    "# Use the function 'ndi.label' from the 'ndimage' module. \n",
    "# Note that this function labels foreground pixels (1s, not 0s), so you may need \n",
    "# to invert your membrane mask just as for hole filling above.\n",
    "# Also, note that 'ndi.label' returns another result in addition to the labeled \n",
    "# image. Read up on this in the function's documentation and make sure you don't\n",
    "# mix up the two outputs!\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Visualize the output\n",
    "\n",
    "# Here, it is no longer ideal to use a 'gray' colormap, since we want to visualize that each\n",
    "# cell has a unique ID. Play around with different colormaps (check the docs to see what\n",
    "# types of colormaps are available) and choose one that you are happy with.\n",
    "\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Take a close look at the picture and note mistakes in the segmentation. Depending on the\n",
    "# quality of your membrane mask, there will most likely be some cells that are 'fused', meaning \n",
    "# two or more cells are labeled as the same cell; this is called \"under-segmentation\". \n",
    "# We will resolve this issue in the next step. Note that our downstream pipeline does not involve \n",
    "# any steps to resolve \"over-segmentation\" (i.e. a cell being wrongly split into multiple labeled\n",
    "# areas), so you should tune your membrane mask such that this is not a common problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cell Segmentation by Seeding & Expansion <a id=seg></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "The segmentation we achieved by membrane masking and connected components labeling is a good start. We could for example use it to measure the fluorescence intensity in each cell's cytoplasm. However, we cannot use it to measure intensities at the membrane of the cells, nor can we use it to accurately measure features like cell shape or size.\n",
    "\n",
    "To improve this (and to resolve cases of under-segmentation), we can use a \"seeding & expansion\" strategy. Expansion algorithms such as the `watershed` start from a small `seed` and \"grow outward\" until they touch the boundaries of neighboring cells, which are themselves growing outward from neighboring seeds. Since the \"growth rate\" at the edge of the growing areas is dependent on image intensity (higher intensity means slower expansion), these expansion methods end up tracing the cells' outlines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seeding by Distance Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "A `seed image` contains a few pixels at the center of each cell labeled by a unique `ID` number and surrounded by zeros. The expansion algorithm will start from these central pixels and grow outward until all zeros are overwritten by an `ID` label. In the case of `watershed` expansion, one can imagine the `seeds` as the sources from which water pours into the cells and starts filling them up.\n",
    "\n",
    "For multi-channel images that contain a nuclear label, it is common practice to mask the nuclei by thresholding and use an eroded version of the nuclei as seeds for cell segmentation. However, there are good alternative seeding approaches for cases where nuclei are not available or not nicely separable by thresholding.\n",
    "\n",
    "Here, we will use a `distance transform` for seeding. In a `distance transform`, each pixel in the foreground (here the cells) is assigned a value corresponding to its distance from the closest background pixel (here the membrane segmentation). In other words, we encode within the image how far each pixel of a cell is away from the membrane (see figure below). The pixels furthest away from the membrane will be at the center of the cells and will have the highest values. Using a function to detect `local maxima`, we will find these high-value peaks and use them as seeds for our segmentation.\n",
    "\n",
    "<img src=\"ipynb_images\\distance_transform.png\" alt=\"Distance Transform Figure\" style=\"width: 900px;\"/>\n",
    "\n",
    "One big advantage of this approach is that it will create two separate seeds even if two cells are connected by a hole in the membrane segmentation. Thus, under-segmentation artifacts will be reduced."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'> Exercise </font>\n",
    "\n",
    "Find seeds using the distance transform approach.\n",
    "\n",
    "This involves the following three steps:\n",
    "\n",
    "1. Run the distance transform on your membrane mask.\n",
    "\n",
    "2. Due to irregularities in the membrane shape, the distance transform may have some smaller local maxima in addition to those at the center of the cells. This will lead to additional seeds, which will lead to over-segmentation. To resolve this problem, smoothen the distance transform using Gaussian smoothing. \n",
    "\n",
    "3. Find the seeds by detecting local maxima. Optimize the seeding by changing the amount of smoothing done in step 2, aiming to have exactly one seed for each cell (although this may not be perfectly achievable).\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Run a distance transform on the membrane mask\n",
    "\n",
    "# Use the function 'ndi.distance_transform_edt'.\n",
    "# You may need to invert your membrane mask so the distances are computed on\n",
    "# the cells, not on the membranes.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Visualize the output and understand what you are seeing.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Smoothen the distance transform\n",
    "\n",
    "# Use 'ndi.gaussian_filter' to do so.\n",
    "# You will have to optimize your choice of 'sigma' based on the outcome below.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iv) Retrieve the local maxima (the 'peaks') from the distance transform\n",
    "\n",
    "# Use the function 'peak_local_max' from the module 'skimage.feature'. By default, this function will return the\n",
    "# indices (aka coordinates) of the pixels where the local maxima are. However, we instead need a boolean mask of \n",
    "# the same shape as the original image, where all the local maximum pixels are labeled as `1` and everything else \n",
    "# as `0`. The documentation for `peak_local_max` shows how to convert the coordinates into a Boolean mask using\n",
    "# numpy (see the last of the examples in the docs). This is a bit of a technical detail, though, so feel free to\n",
    "# copy the conversion from the solutions.\n",
    "\n",
    "# Retrieve peak coordinates\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Convert coords to mask as per skimage documentation\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (v) Visualize the output as an overlay on the raw (or smoothed) image\n",
    "\n",
    "# If you just look at the local maxima image, it will simply look like a bunch of distributed dots.\n",
    "# To get an idea if the seeds are well-placed, you will need to overlay these dots onto the original image.\n",
    "\n",
    "# To do this, it is important to first understand a key point about how the 'pyplot' module works: \n",
    "# every plotting command is slapped on top of the previous plotting commands, until everything is ultimately \n",
    "# shown when 'plt.show' is called. Hence, you can first plot the raw (or smoothed) input image and then\n",
    "# plot the seeds on top of it before showing both with 'plt.show'.\n",
    "\n",
    "# As you can see if you try this, you will not get the desired result because the zero values in seed array\n",
    "# are painted in black over the image you want in the background. To solve this problem, you need to mask \n",
    "# these zero values before plotting the seeds. You can do this by creating an appropriately masked array\n",
    "# using the function 'np.ma.array' with the keyword argument 'mask'. \n",
    "# Check the docs or Stack Overflow to figure out how to do this.\n",
    "\n",
    "# BONUS: As an additional improvement for the visualization, use 'ndi.maximum_filter' to dilate the \n",
    "# seeds a little bit, making them bigger and thus better visible.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (vi) Optimize the seeding\n",
    "\n",
    "# Ideally, there should be exactly one seed for each cell.\n",
    "# If you are not satisfied with your seeding, go back to the smoothing step above and optimize 'sigma'\n",
    "# to get rid of additional maxima. You can also try using the keyword argument 'min_distance' in \n",
    "# 'peak_local_max' to solve cases where there are multiple small seeds at the center of a cell. Note \n",
    "# that good seeding is essential for a good segmentation with an expansion algorithm. However, no \n",
    "# segmentation is perfect, so it's okay if a few cells end up being oversegmented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (vii) Label the seeds\n",
    "\n",
    "# Use connected component labeling to give each cell seed a unique ID number.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Visualize the final result (the labeled seeds) as an overlay on the raw (or smoothed) image\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Expansion by Watershed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "To achieve a cell segmentation, the `seeds` now need to be expanded outward until they follow the outline of the cell. The most commonly used expansion algorithm is the `watershed`.\n",
    "\n",
    "Imagine the intensity in the raw/smoothed image as a topographical height profile; high-intensity regions are peaks, low-intensity regions are valleys. In this representation, cells are deep valleys (with the seeds at the center), enclosed by mountains. As the name suggests, the `watershed` algorithm can be understood as the gradual filling of this landscape with water, starting from the seed. As the water level rises, the seed expands - until it finally reaches the 'crest' of the cell membrane 'mountain range'. Here, the water would flow over into the neighboring valley, but since that valley is itself filled up with water from the neighboring cell's seed, the two water surfaces touch and the expansion stops.\n",
    "\n",
    "<img src=\"ipynb_images\\watershed_illustration.png\" alt=\"Watershed Figure\" style=\"width: 900px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "Expand your seeds by means of a watershed expansion.\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Perform watershed\n",
    "\n",
    "# Use the function 'watershed' from the module 'skimage.segmentation'.\n",
    "# Use the labeled cell seeds and the smoothed membrane image as input.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Show the result as transparent overlay over the smoothed input image\n",
    "\n",
    "# Like the masked overlay of the seeds, this can be achieved by making two calls to 'imshow',\n",
    "# one for the background image and one for the segmentation. Instead of masking away background,\n",
    "# this time you simply make the segmentation image semi-transparent by adjusting the keyword \n",
    "# argument 'alpha' of the 'imshow' function, which specifies opacity.\n",
    "# Be sure to choose an appropriate colormap that allows you to distinguish the segmented cells\n",
    "# even if cells with a very similar ID are next to each other (I would recommend 'prism').\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *A Note on Segmentation Quality*\n",
    "\n",
    "This concludes the segmentation of the cells in the example image. Depending on the quality you achieved in each step along the way, the final segmentation may be of greater or lesser quality (in terms of over-/under-segmentation errors).\n",
    "\n",
    "It should be noted that the segmentation will likely *never* be perfect, as there is usually a trade-off between over- and undersegmentation.\n",
    "\n",
    "This raises an important question: ***When should I stop trying to optimize my segmentation?***\n",
    "\n",
    "There is no absolute answer to this question but the best answer is probably this: ***When you can use it to address your biological questions!***\n",
    "\n",
    "*Importantly, this implies that you should already have relatively clear questions in mind when you are working on the segmentation!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Postprocessing: Removing Cells at the Image Border <a id=postpro></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "Since segmentation is never perfect, it often makes sense to explicitely remove artifacts afterwards. For example, one could filter out objects that are too small, have a very strange shape, or very strange intensity values. \n",
    "\n",
    "**Warning:** Filtering out objects is equivalent to the *removal of outliers* in data analysis and *should only be done for good reason and with caution!*\n",
    "\n",
    "As an example of postprocessing, we will now filter out a particular group of problematic cells: those that are being cut off at the image border."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'> Exercise </font>\n",
    "\n",
    "Iterate through all the cells in your segmentation and remove those touching the image border.\n",
    "\n",
    "Follow the instructions in the comments below. Note that the instructions will get a little less specific from here on, so you need to figure out how to approach a problem yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Create an image border mask\n",
    "\n",
    "# We need some way to check if a cell is at the border. For this, we generate a 'mask' of the image border,\n",
    "# i.e. a Boolean array of the same size as the image where only the border pixels are set to `1` and all \n",
    "# others to `0`, like this:\n",
    "#   1 1 1 1 1\n",
    "#   1 0 0 0 1\n",
    "#   1 0 0 0 1\n",
    "#   1 0 0 0 1\n",
    "#   1 1 1 1 1\n",
    "# There are multiple ways of generating this mask, for example by erosion or by array indexing.\n",
    "# It is up to you to find a way to do it.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) 'Delete' the cells at the border\n",
    "\n",
    "# When modifying a segmentation (in this case by deleting some cells), it makes sense\n",
    "# to work on a copy of the array, not on the original. This avoids unexpected behaviors,\n",
    "# especially within jupyter notebooks. Use the function 'np.copy' to copy an array.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Iterate over the IDs of all the cells in the segmentation. Use a for-loop and the \n",
    "# function 'np.unique' (remember that each cell in our segmentation is labeled with a \n",
    "# different integer value).\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "    # Create a mask that contains only the 'current' cell of the iteration\n",
    "    # Hint: Remember that the comparison of an array with some number (array==number)\n",
    "    #       returns a Boolean mask of the pixels in 'array' whose value is 'number'.\n",
    "    ### YOUR CODE HERE!\n",
    "    \n",
    "    # Using the cell mask and the border mask from above, test if the cell has pixels touching \n",
    "    # the image border or not.\n",
    "    # Hint: 'np.logical_and'\n",
    "    ### YOUR CODE HERE!\n",
    "\n",
    "    # If a cell touches the image boundary, delete it by setting its pixels in the segmentation to 0.\n",
    "    ### YOUR CODE HERE!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OPTIONAL: re-label the remaining cells to keep the numbering consistent from 1 to N (with 0 as background).\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Visualize the result\n",
    "\n",
    "# Show the result as transparent overlay over the raw or smoothed image. \n",
    "# Here you have to combine alpha (to make cells transparent) and 'np.ma.array'\n",
    "# (to hide empty space where the border cells were deleted).\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identifying Cell Edges <a id=edges></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "With the final segmentation in hand, we can now start to think about measurements and data analysis. However, to extract interesting measurements from our cells, the segmentation on its own is often not enough: additional masks that identify sub-regions for each cell allow more precise and more biologically relevant measurements.\n",
    "\n",
    "The most useful example of this is an additional mask that identifies only the edge pixels of each cell. This is useful for a number of purposes, including:\n",
    "\n",
    "- Edge intensity is a good measure of membrane intensity, which is often a desired readout.\n",
    "- The intensity profile along the edge may contain information on cell polarity.\n",
    "- The length of the edge (relative to the cell area) is an informative feature about the cell shape. \n",
    "- Showing colored edges is a nice way of visualizing cell segmentations.\n",
    "\n",
    "There are many ways of identifying edge pixels in a fully labeled segmentation. Here, we will use a simple and relatively fast method based on erosion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'> Exercise </font>\n",
    "\n",
    "Create a labeled mask of cell edges by following these steps:\n",
    "\n",
    "\n",
    "- Create an array of the same size and data type as the segmentation but filled with only zeros\n",
    "    - This will be your final cell edge mask; you gradually add cell edges as you iterate over cells\n",
    "    \n",
    "\n",
    "- *For each cell...*\n",
    "    - Erode the cell's mask by 1 pixel\n",
    "    - Using the eroded mask and the original mask, create a new mask of only the cell's edge pixels\n",
    "    - Add the cell's edge pixels into the empty image generated above, labeling them with the cell's original ID number\n",
    "\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Create an array of the same size and data type as the segmentation but filled with only zeros\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Iterate over the cell IDs\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "    # (iii) Erode the cell's mask by 1 pixel\n",
    "    # Hint: 'ndi.binary_erode'\n",
    "    ### YOUR CODE HERE!\n",
    "    \n",
    "    # (iv) Create the cell edge mask\n",
    "    # Hint: 'np.logical_xor'\n",
    "    ### YOUR CODE HERE!\n",
    "    \n",
    "    # (v) Add the cell edge mask to the empty array generated above, labeling it with the cell's ID\n",
    "    ### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (vi) Visualize the result\n",
    "\n",
    "# Note: Because the lines are so thin (1pxl wide), they may not be displayed correctly in small figures.\n",
    "#       You can 'zoom in' by showing a sub-region of the image which is then rendered bigger. You can\n",
    "#       also go back to the edge identification code and make the edges multiple pixels wide (but keep \n",
    "#       in mind that this will have an effect on your quantification results!).\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting Quantitative Measurements <a id=measure></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "The ultimate goal of image segmentation is of course the extraction of quantitative measurements, in this case on a single-cell level. Measures of interest can be based on intensity (in different channels) or on the size and shape of the cells.\n",
    "\n",
    "To exemplify how different properties of cells can be measured, we will extract the following:\n",
    "\n",
    "- Cell ID (so all other measurements can be traced back to the cell that was measured)\n",
    "- Mean intensity of each cell\n",
    "- Mean intensity at the membrane of each cell\n",
    "- The cell area, i.e. the number of pixels that make up the cell\n",
    "- The cell outline length, i.e. the number of pixels that make up the cell edge\n",
    "\n",
    "*Note: It makes sense to use smoothed/filtered/background-subtracted images for segmentation. When it comes to measurements, however, it's best to get back to the raw data!*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "Extract the measurements listed above for each cell and collect them in a dictionary.\n",
    "\n",
    "Note: The ideal data structure for data like this is the `DataFrame` offered by the module `Pandas`. However, for the sake of simplicity, we will here stick with a dictionary of lists.\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Create a dictionary that contains a key-value pairing for each measurement\n",
    "\n",
    "# The keys should be strings describing the type of measurement (e.g. 'intensity_mean') and \n",
    "# the values should be empty lists. These empty lists will be filled with the results of the\n",
    "# measurements.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Record the measurements for each cell\n",
    "\n",
    "# Iterate over the segmented cells ('np.unique').\n",
    "# Inside the loop, create a mask for the current cell and use it to extract the measurements listed above. \n",
    "# Add them to the appropriate list in the dictionary using the 'append' method.\n",
    "# Hint: Remember that you can get out all the values within a masked area by indexing the image \n",
    "#       with the mask. For example, 'np.mean(image[cell_mask])' will return the mean of all the \n",
    "#       intensity values of 'image' that are masked by 'cell_mask'!\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Print the results and check that they make sense\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Analysis & Visualisation <a id=analysis></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "By extracting quantitative measurements from an image we cross over from 'image analysis' to 'data analysis'. \n",
    "\n",
    "This section briefly explains how to do basic data analysis and plotting, including boxplots, scatterplots and linear fits. It also showcases how to map data back onto the image, creating an \"image-based heatmap\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "Analyze and plot the extracted data in a variety of ways.\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Familiarize yourself with the data structure of the results dict and summarize the results\n",
    "\n",
    "# Recall that dictionaries are unordered; a dataset of interest is accessed through its key.\n",
    "# In our case, the datasets inside the dict are lists of values, ordered in the same order\n",
    "# as the cell IDs. \n",
    "\n",
    "# For each dataset in the results dict, print its name (the key) along with its mean, standard \n",
    "# deviation, maximum, minimum, and median. The appropriate numpy methods (e.g. 'np.median') work\n",
    "# with lists just as well as with arrays.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Create a box plot showing the mean cell and mean membrane intensities for both channels. \n",
    "\n",
    "# Use the function 'plt.boxplot'. Use the 'label' keyword of 'plt.boxplot' to label the x axis with \n",
    "# the corresponding key names. Feel free to play around with the various options of the boxplot \n",
    "# function to make your plot look nicer. Remember that you can first call 'plt.figure' to adjust \n",
    "# settings such as the size of the plot.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Create a scatter plot of cell outline length over cell area\n",
    "\n",
    "# Use the function 'plt.scatter' for this. Be sure to properly label the \n",
    "# plot using 'plt.xlabel' and 'plt.ylabel'.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# BONUS: Do you understand why you are seeing the pattern this produces? Can you\n",
    "# generate a 'null model' curve that assumes all cells to be circular? What is\n",
    "# the result? Do you notice something odd about it? What could be the reason for\n",
    "# this and how could it be fixed?\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iv) Perform a linear fit of membrane intensity over cell area\n",
    "\n",
    "# Use the function 'linregress' from the module 'scipy.stats'. Be sure to read the docs to\n",
    "# understand the output of this function. Print the output.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (v) Think about the result\n",
    "\n",
    "# Note that the fit seems to return a highly significant p-value but a very low correlation \n",
    "# coefficient (r-value). Based on prior knowledge, we would not expect a linear correlation of \n",
    "# this sort to be present in our data. \n",
    "#\n",
    "# This should prompt several questions:\n",
    "#   1) What does this p-value actually mean? Check the docs of 'linregress'!\n",
    "#   2) Could there be artifacts in our segmentation that bias this analysis?\n",
    "#\n",
    "# In general, it's always good to be very careful when doing any kind of data analysis. Make sure you \n",
    "# understand the functions you are using and always check for possible errors or sources of bias!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (vi) Overlay the linear fit onto a scatter plot\n",
    "\n",
    "# Recall that a linear function is defined by `y = slope * x + intercept`.\n",
    "\n",
    "# To define the line you'd like to plot, you need two values of x (the starting point and\n",
    "# and the end point of the line). What values of x make sense? Can you get them automatically?\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# When you have the x-values for the starting point and end point, get the corresponding y \n",
    "# values from the fit through the equation above.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Plot the line with 'plt.plot'. Adjust the line's properties so it is well visible.\n",
    "# Note: Remember that you have to create the scatterplot before plotting the line so that\n",
    "#       the line will be placed on top of the scatterplot.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Use 'plt.legend' to add information about the line to the plot.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Label the plot and finally show it with 'plt.show'.\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (vii) Map the cell area back onto the image as a 'heatmap'\n",
    "\n",
    "# Scale the cell area data to 8bit so that it can be used as pixel intensity values.\n",
    "# Hint: if the largest cell area should correspond to the value 255 in uint8, then \n",
    "#       the other cell areas correspond to 'cell_area * 255 / largest_cell_area'.\n",
    "# Hint: To perform an operation on all cell area values at once, convert the list \n",
    "#       of cell areas to a numpy array.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Initialize a new image array; all values should be zeros, the shape should be identical \n",
    "# to the images we worked with before and the dtype should be uint8.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "# Iterate over the segmented cells. In addition to the cell IDs, the for-loop should\n",
    "# also include a simple counter (starting from 0) with which the area measurement can be \n",
    "# accessed by indexing.\n",
    "### YOUR CODE HERE!\n",
    "    \n",
    "    # Mask the current cell and assign the cell's (re-scaled) area value to the cell's pixels.\n",
    "    ### YOUR CODE HERE!\n",
    "\n",
    "# Visualize the result as a colored semi-transparent overlay over the raw/smoothed original input image.\n",
    "# BONUS: See if you can exclude outliers to make the color mapping more informative!\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing Output to Files <a id=write></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "The final step of the pipeline shows how to write various outputs of the pipeline to files.\n",
    "\n",
    "Data can be saved to files in a human-readable format such as text files (e.g. to import into Excel), in a format readable for other programs such as tif-images (e.g. to view in Fiji) or in a language-specific file that makes it easy to reload the data into python in the future (e.g. for further analysis)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'> Exercise </font>\n",
    "\n",
    "Write the generated data into a variety of different output files.\n",
    "\n",
    "Follow the instructions in the comments below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Write one or more of the images you produced to a tif file\n",
    "\n",
    "# Use the function 'imsave' from the 'skimage.io' module. Make sure that the array you are \n",
    "# writing is of integer type. If necessary, you can use the method 'astype' for conversions, \n",
    "# e.g. 'some_array.astype(np.uint8)' or 'some_array.astype(np.uint16)'. Careful when \n",
    "# converting a segmentation to uint8; if there are more than 255 cells, the 8bit format\n",
    "# doesn't have sufficient bit-depth to represent all cell IDs!\n",
    "#\n",
    "# You can also try adding the segmentation to the original image, creating an image with\n",
    "# two channels, one of them being the segmentation. \n",
    "#\n",
    "# After writing the file, load it into Fiji and check that everything worked as intended.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Write a figure to a png or pdf\n",
    "\n",
    "# Recreate the scatter plot from above (with or without the regression line), then save the figure\n",
    "# as a png using 'plt.savefig'. Alternatively, you can also save it to a pdf, which will create a\n",
    "# vector graphic that can be imported into programs like Adobe Illustrator.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Save the segmentation as a numpy file\n",
    "\n",
    "# Numpy files allow fast storage and reloading of numpy arrays. Use the function 'np.save'\n",
    "# to save the array and reload it using 'np.load'.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iv) Save the result dictionary as a pickle file\n",
    "\n",
    "# Pickling is a way of generating generic files from almost any python object, which can easily\n",
    "# be reloaded into python at a later point in time.\n",
    "# You will need to open an empty file object using 'open' in write-bytes mode ('wb'). It's best to  \n",
    "# do so using the 'with'-statement (context manager) to make sure that the file object will be \n",
    "# closed automatically when you are done with it.\n",
    "# Use the function 'pickle.dump' from the 'pickle' module to write the results to the file.\n",
    "# Hint: Refer to the python documention for input and output to understand how file objects are\n",
    "#       handled in python in general.\n",
    "\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "## Note: Pickled files can be re-loaded again as follows:\n",
    "#with open('my_filename.pkl', 'rb') as infile:\n",
    "#    reloaded = pickle.load(infile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (v) Write a tab-separated text file of the results dict\n",
    "\n",
    "# The most generic way of saving numeric results is a simple text file. It can be imported into \n",
    "# pretty much any other program.\n",
    "\n",
    "# To write normal text files, open an empty file object in write mode ('w') using the 'with'-statement.\n",
    "### YOUR CODE HERE!\n",
    "\n",
    "    # Use the 'file_object.write(string)' method to write strings to the file, one line at a time,\n",
    "    # First, write the header of the data (the result dict keys), separated by tabs ('\\t'). \n",
    "    # It makes sense to first generate a complete string with all the headers and then write this \n",
    "    # string to the file as one line. Note that you will need to explicitly write 'newline' characters \n",
    "    # ('\\n') at the end of the line to switch to the next line.\n",
    "    # Hint: the string method 'join' is very useful here!\n",
    "    ### YOUR CODE HERE!\n",
    "\n",
    "    # After writing the headers, iterate over all the cells and write the result data to the file line\n",
    "    # by line, by creating strings similar to the header string.\n",
    "    ### YOUR CODE HERE!\n",
    "\n",
    "# After writing the data, have a look at the output file in a text editor or in a spreadsheet\n",
    "# program like Excel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Batch Processing <a id=batch></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Background\n",
    "\n",
    "In practice, we never work with just a single image, so we would like to make it possible to run our analysis pipeline for multiple images and then collect and analyze all the results. This final section of the tutorial shows how to do just that."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='teal'>Exercise</font>\n",
    "\n",
    "To run a pipeline multiple times, it needs to be packaged into a function or - even better - as a separate module. Jupyter notebook is not well suited for this, so if you're working in a notebook, first extract your code to a `.py` file (see instructions below). If you are not working in a notebook, create a copy of your pipeline; we will modify this copy into a function that can then be called repeatedly for different images.\n",
    "\n",
    "To export a jupyter notebook as a `.py` file, use `File > Download as > Python (.py)`, then save the file. Open the resulting python script in a text editor or in an IDE like PyCharm. \n",
    "\n",
    "\n",
    "#### Let's clean the script a bit:\n",
    "\n",
    "- Remove the line `%matplotlib [inline|notebook|qt]`. It is not valid python code outside of a Jupyter notebook.\n",
    "\n",
    "\n",
    "- Go through the script and comment out everything related to plotting; when running a pipeline for dozens or hundreds of images, we usually do not want to generate tons of plots. Similarly, it can make sense to remove some print statments if you have many of them.\n",
    "\n",
    "\n",
    "- Remove the sections `Manual Thresholding` and `Connected Components Labeling`; they are not used in the final segmentation.\n",
    "\n",
    "\n",
    "- Remove the sections `Simple Analysis and Visualization` and `Writing Output to Files`; we will collect the output for each image when running the pipeline in a loop. That way, everything can be analyzed at once at the end. \n",
    "    - Note that, even though we skip it here, it is often very useful to store every input file's corresponding outputs in new files. When doing so, the output files should use the name of the input file modified with an additional suffix. For example, the results extracted when analyzing `img_1.tif` might best be stored as `img_1_results.pkl`.\n",
    "    - You can implement this approach for saving the segmentations and/or the result dicts as a *bonus* exercise!\n",
    "\n",
    "\n",
    "- Feel free to delete some of the background information to make the script more concise.\n",
    "\n",
    "\n",
    "#### Converting the pipeline to a function:\n",
    "\n",
    "Convert the entire pipeline into a function that accepts a directory and a filename as input, runs everything, and returns the final segmentation and the results dictionary. To do this, you must:\n",
    "\n",
    "- Add the function definition statement at the beginning of the script (after the imports)\n",
    "- Replace the 'hard-coded' directory path and filename by variables that are accepted by the function\n",
    "- Indent all the code\n",
    "- Add a return statement at the end\n",
    "\n",
    "\n",
    "#### Importing the function and running it for multiple input files:\n",
    "\n",
    "To actually run the pipeline function for multiple input files, we need to do the following:\n",
    "\n",
    "- Import the pipeline function from the `.py` file\n",
    "- Iterate over all the filenames in a directory\n",
    "- For each filename, call the pipeline function\n",
    "- Collect the returned results\n",
    "\n",
    "Once you have converted your pipeline into a function as described above, you can import and run it according to the instructions below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Test if your pipeline function actually works\n",
    "\n",
    "# Import your function using the normal python syntax for imports, like this:\n",
    "#   from your_module import your_function\n",
    "# Run the function and visualize the resulting segmentation. Make sure everything works as intended.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Get all relevant filenames from the input directory\n",
    "\n",
    "# Use the function 'listdir' from the module 'os' to get a list of all the files\n",
    "# in a directory. Find a way to filter out only the relevant input files, namely\n",
    "# \"example_cells_1.tif\" and \"example_cells_2.tif\". Of course, one would usually\n",
    "# do this for many more images, otherwise it's not worth the effort.\n",
    "# Hint: Loop over the filenames and use if statements to decide which ones to \n",
    "#       keep and which ones to throw away.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Iterate over the input filenames and run the pipeline function\n",
    "\n",
    "# Be sure to collect the output of the pipeline function in a way that allows\n",
    "# you to trace it back to the file it came from. You could for example use a\n",
    "# dictionary with the filenames as keys.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iv) Recreate one of the scatterplots from above but this time with all the cells\n",
    "\n",
    "# You can color-code the dots to indicate which file they came from. Don't forget to\n",
    "# add a corresponding legend.\n",
    "\n",
    "### YOUR CODE HERE!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='teal'>*Congratulations! You have completed the tutorial!*</font>\n",
    "\n",
    "**We hope you enjoyed the ride and learned a lot!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Concluding Remarks\n",
    "\n",
    "It's important to remember that the phrase ***\"Use it or loose it!\"*** fully applies for the skills taught in this tutorial.\n",
    "\n",
    "If you now just go back to the lab and don't touch python or image analysis for the next half year, most of the things you have learned here will be lost.\n",
    "\n",
    "So, what can you do?\n",
    "\n",
    "\n",
    "- If possible, start applying what you have learned to your own work right away\n",
    "\n",
    "\n",
    "- Even if your current work doesn't absolutely *need* coding / image analysis (which to be honest is hard to believe! ;p), you can still use it at least to make some nice plots!\n",
    "\n",
    "\n",
    "- Another very good approach is to find yourself an interesting little side project you can play around with\n",
    "\n",
    "\n",
    "- Of course, there is still much more to learn and the internet happens to be full of excellent tutorials!\n",
    "    - As a starting point, have a look at [Bio-IT's curated list of tutorials](https://bio-it.embl.de/coding-club/curated-tutorials/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***We wish you the best of luck for all your coding endeavors!***"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}