image_analysis_tutorial_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Analysis with Python - <font color='orange'>Tutorial Pipeline Solutions</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. [Importing Modules & Packages](#import)\n",
    "3. [Loading & Handling Image Data](#load)\n",
    "4. [Preprocessing](#prepro)\n",
    "5. [Manual Thresholding & Threshold Detection](#thresh)\n",
    "6. [Adaptive Thresholding](#adaptive)\n",
    "7. [Improving Masks with Binary Morphology](#morpho)\n",
    "8. [Connected Components Labeling](#label)\n",
    "9. [Cell Segmentation by Seeding & Expansion](#seg)\n",
    "10. [Postprocessing: Removing Cells at the Image Border](#postpro)\n",
    "11. [Identifying Cell Edges](#edges)\n",
    "12. [Extracting Quantitative Measurements](#measure)\n",
    "13. [Simple Analysis & Visualization](#analysis)\n",
    "14. [Writing Output to Files](#write)\n",
    "15. [Batch Processing](#batch)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  About this Tutorial - <font color='orange'> Notes on the Solutions </font> <a id=about></a>\n",
    "\n",
    "These are the solutions to the tutorial pipeline (see `image_analysis_tutorial.py`). <font color=red>*For the best learning experience, it is highly recommended that you first try to implement a solution yourself!*</font> Only come here if you are totally stuck or if you have a working solution and would like to double-check it.\n",
    "\n",
    "Note that there are multiple ways of implementing any particular step in the pipeline, so if your solution is different from the solution here, it is not necessarily wrong. However, *some solutions are better than others* because they are...\n",
    "\n",
    "- **...more readable:**\n",
    "    - When reading the code, it is obvious and clear what the code is doing\n",
    "    - The code is clearly commented to help others (and your future self) understand it\n",
    "    \n",
    "    \n",
    "- **...more general:**\n",
    "    - The code still works if there are minor changes to the data (e.g. size of the image being processed)\n",
    "    - The code can easily be transformed into a solution for a similar problem\n",
    "    \n",
    "    \n",
    "- **...more computationally efficient:**\n",
    "    - No unnecessary copies of large datasets are made (memory efficiency)\n",
    "    - Faster algorithms are used, e.g. array operations instead of loops (CPU efficiency)\n",
    "    - No unnecessary loading and writing of data (io efficiency)\n",
    "\n",
    "\n",
    "It is up to you to decide if your solution is better, equally good, or not as good as the solution presented here. Either way, we hope you can learn something by looking at our solutions. **:)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importing Modules & Packages <a id=import></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The numerical array package numpy as np\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</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",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# The image processing package scipy.ndimage as ndi\n",
    "import scipy.ndimage as ndi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Side Note for Jupyter Notebook Users\n",
    "\n",
    "You can configure how the figures made by matplotlib are displayed."
   ]
  },
  {
   "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": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "We will now proceed to load one of the example images and verify that we get what we expect."
   ]
  },
  {
   "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",
    "filename = r'example_cells_1.tif'\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",
    "dirpath = r'example_data'  # Relative path\n",
    "#dirpath = r'/home/jack/data/python_image_analysis/example_data/example_cells_1.tif'  # Absolute path "
   ]
  },
  {
   "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",
    "from os.path import join\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",
    "filepath = join(dirpath, filename)\n",
    "print(filepath)"
   ]
  },
  {
   "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",
    "from skimage.io import imread\n",
    "\n",
    "# Load 'example_cells_1.tif' and store it in a variable.\n",
    "# Suggested name for the variable: img\n",
    "img = imread(filepath)"
   ]
  },
  {
   "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",
    "print(\"Loaded array is of type:\", type(img))\n",
    "\n",
    "# Print the shape of the array using the numpy-function 'shape'. \n",
    "# Make sure you understand the output!\n",
    "print(\"Loaded array has shape:\", img.shape)\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",
    "print(\"Loaded values are of type:\", img.dtype)\n",
    "\n",
    "# SOLUTION NOTE: The dtype should be 'uint8', because these are unsigned 8-bit integer images.\n",
    "#                This means that the intensity values range from 0 to 255 in steps of 1."
   ]
  },
  {
   "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 the 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",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(img, interpolation='none', cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing <a id=prepro></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Perform Gaussian smoothing and visualize the result."
   ]
  },
  {
   "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",
    "sigma = 3\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",
    "img_smooth = ndi.gaussian_filter(img, sigma)"
   ]
  },
  {
   "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",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "plt.show()\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",
    "plt.figure(figsize=(6,6))\n",
    "plt.imshow(img_smooth[400:600, 200:400], interpolation='none', cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "fig, ax = plt.subplots(1, 2, figsize=(10,7))\n",
    "ax[0].imshow(img, interpolation='none', cmap='gray')\n",
    "ax[1].imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "ax[0].set_title('Raw Image')\n",
    "ax[1].set_title('Smoothed Image')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Manual Thresholding & Threshold Detection <a id=thresh></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Try out manual thresholding and automated threshold detection."
   ]
  },
  {
   "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",
    "thresh = 70"
   ]
  },
  {
   "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",
    "mem = img_smooth > thresh\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",
    "print(mem.dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Visualize the result\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10,7))\n",
    "ax[0].imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "ax[1].imshow(mem, interpolation='none', cmap='gray')\n",
    "ax[0].set_title('Smoothed Image')\n",
    "ax[1].set_title('Thresholded Membranes')\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "\n",
    "# Prepare widget\n",
    "from ipywidgets import interact\n",
    "@interact(thresh=(10,250,10))\n",
    "def select_threshold(thresh=100):\n",
    "    \n",
    "    # Thresholding\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": "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",
    "\n",
    "# Import\n",
    "from skimage.filters.thresholding import threshold_otsu\n",
    "\n",
    "# Calculate and apply threshold\n",
    "thresh = threshold_otsu(img_smooth)\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": "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",
    "from skimage.filters.thresholding import try_all_threshold\n",
    "fig = try_all_threshold(img_smooth, figsize=(10,10), verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive Thresholding <a id=adaptive></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Implement the two steps of adaptive background subtraction:\n",
    "\n",
    "1. Use a strong \"mean filter\" with a circular `SE` to create the background image.\n",
    "\n",
    "2. Use the background image for thresholding."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1\n",
    "# ------\n",
    "\n",
    "# The expression below creates circular 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",
    "#\n",
    "#   - Is the resulting SE really circular?  \n",
    "#       >>> It is close enough for large i, not so much for small i.\n",
    "#\n",
    "#   - Can certain values of 'i' cause problems? If so, why?  \n",
    "#       >>> Even values create a slight asymmetry! It is best to use only odd values!\n",
    "#\n",
    "#   - What value of 'i' should be used for the SE?\n",
    "#       >>> My first guess was i=30, which is about 3x the membrane diameter. \n",
    "#           I tried out some other values but ultimately stuck with this.\n",
    "\n",
    "# Create SE\n",
    "i = 31\n",
    "SE = (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",
    "plt.imshow(SE, interpolation='none')\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "from skimage.filters import rank \n",
    "bg = rank.mean(img_smooth, footprint=SE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Visualize the resulting background image. Does what you get make sense?\n",
    "\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(bg, interpolation='none', cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "mem = img_smooth > bg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (v) Visualize and understand the output. \n",
    "\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(mem, interpolation='none', cmap='gray')\n",
    "plt.show()\n",
    "\n",
    "# Are you happy with this result as a membrane segmentation? \n",
    "# ->> It's pretty good, except for the speckles that litter the inside of the cells!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improving Masks with Binary Morphology <a id=morpho></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Improve the membrane segmentation from above with morphological operations."
   ]
  },
  {
   "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",
    "#mem_holefilled = np.logical_not(ndi.binary_fill_holes(np.logical_not(mem)))  # Long form\n",
    "mem_holefilled = ~ndi.binary_fill_holes(~mem)  # Short form"
   ]
  },
  {
   "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",
    "# artifacts at the image boundaries. Can you come up with a way of solving this?\n",
    "\n",
    "# New circular SE of appropriate size (determined by trial and error)\n",
    "i = 15\n",
    "SE = (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",
    "# One solution to the boundary artefact issue is padding with the reflection.\n",
    "# 'Padding' refers to the extension of the image at the boundaries, in this case using\n",
    "# a 'reflection' of the pixel values next to the boundary. If morphological operations \n",
    "# are done on the padded image, the boundary artefacts will occur in the padded region\n",
    "# outside the original image, which can simply be cropped out again at the end.\n",
    "pad_size = i+1\n",
    "mem_padded = np.pad(mem_holefilled, pad_size, mode='reflect')\n",
    "\n",
    "# Binary closing works well to round off the membranes and close gaps\n",
    "mem_final = ndi.binary_closing(mem_padded, structure=SE)\n",
    "\n",
    "# This slicing operation crops the padded image back to the original size\n",
    "mem_final = mem_final[pad_size:-pad_size, pad_size:-pad_size]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Visualize the final result\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10,7))\n",
    "ax[0].imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "ax[1].imshow(mem_final, interpolation='none', cmap='gray')\n",
    "ax[0].set_title('Smoothed Image')\n",
    "ax[1].set_title('Final Membrane Mask')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connected Components Labeling <a id=label></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Use your membrane segmentation for connected components labeling."
   ]
  },
  {
   "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",
    "cell_labels, _ = ndi.label(~mem_final)\n",
    "\n",
    "# Solution note: For functions with multiple outputs (here the labeled image on\n",
    "# the one hand and the number of detected objects on the other), it is a python\n",
    "# convention to unpack those outputs that will not be used in the remainder of \n",
    "# code into the variable '_' (underscore). This makes it clear to those reading\n",
    "# the code that the function returns multiple things but some of them are not\n",
    "# important in this case."
   ]
  },
  {
   "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",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(cell_labels, interpolation='none', cmap='inferno')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cell Segmentation by Seeding & Expansion <a id=seg></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Find seeds using the distance transform approach."
   ]
  },
  {
   "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",
    "dist_trans = ndi.distance_transform_edt(~mem_final)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Visualize the output and understand what you are seeing.\n",
    "\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(dist_trans, interpolation='none', cmap='viridis')\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "\n",
    "# Applying the filter\n",
    "dist_trans_smooth = ndi.gaussian_filter(dist_trans, sigma=5)\n",
    "\n",
    "# Visualizing\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(dist_trans_smooth, interpolation='none', cmap='viridis')\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "from skimage.feature import peak_local_max\n",
    "seed_coords = peak_local_max(dist_trans_smooth, min_distance=10)\n",
    "\n",
    "# Convert coords to mask as per skimage documentation\n",
    "seeds = np.zeros_like(dist_trans_smooth, dtype=bool)\n",
    "seeds[tuple(seed_coords.T)] = True"
   ]
  },
  {
   "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 to the visualization, use 'ndi.maximum_filter' to dilate the \n",
    "# seeds a little bit, making them bigger and thus better visible.\n",
    "\n",
    "# Dilate seeds\n",
    "seeds_dil = ndi.maximum_filter(seeds, size=10)\n",
    "\n",
    "# Create plot\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "plt.imshow(np.ma.array(seeds_dil, mask=seeds_dil==0), interpolation='none', cmap='autumn')\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "seeds_labeled = ndi.label(seeds)[0]\n",
    "\n",
    "# Visualize the final result (the labeled seeds) as an overlay on the raw (or smoothed) image\n",
    "seeds_labeled_dil = ndi.maximum_filter(seeds_labeled, size=10)  # Expand a bit for visualization\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "plt.imshow(np.ma.array(seeds_labeled_dil, mask=seeds_labeled_dil==0), interpolation='none', cmap='prism')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Expansion by Watershed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Expand your seeds by means of a watershed expansion."
   ]
  },
  {
   "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",
    "from skimage.segmentation import watershed\n",
    "ws = watershed(img_smooth, seeds_labeled)"
   ]
  },
  {
   "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",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "plt.imshow(ws, interpolation='none', cmap='prism', alpha=0.4)\n",
    "plt.show()"
   ]
  },
  {
   "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": [
    "#### <font color='orange'> Exercise Solution </font>\n",
    "\n",
    "Iterate through all the cells in your segmentation and remove those touching the image border."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (i) Create an image border mask\n",
    "\n",
    "border_mask = np.zeros(ws.shape, dtype=bool)\n",
    "border_mask = ndi.binary_dilation(border_mask, border_value=1)"
   ]
  },
  {
   "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",
    "clean_ws = np.copy(ws)\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",
    "for cell_ID in np.unique(ws):\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",
    "    cell_mask = ws==cell_ID \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",
    "    cell_border_overlap = np.logical_and(cell_mask, border_mask)  # Overlap of cell mask and boundary mask\n",
    "    total_overlap_pixels = np.sum(cell_border_overlap)            # Sum overlapping pixels\n",
    "\n",
    "    # If a cell touches the image boundary, delete it by setting its pixels in the segmentation to 0.\n",
    "    if total_overlap_pixels > 0: \n",
    "        clean_ws[cell_mask] = 0"
   ]
  },
  {
   "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",
    "for new_ID, cell_ID in enumerate(np.unique(clean_ws)[1:]):  # The [1:] excludes 0 from the list (background)!\n",
    "    clean_ws[clean_ws==cell_ID] = new_ID+1                  # The same here for the +1"
   ]
  },
  {
   "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",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "plt.imshow(np.ma.array(clean_ws, mask=clean_ws==0), interpolation='none', cmap='prism', alpha=0.4)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identifying Cell Edges <a id=edges></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Create a labeled mask of cell edges."
   ]
  },
  {
   "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",
    "edges = np.zeros_like(clean_ws)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (ii) Iterate over the cell IDs\n",
    "for cell_ID in np.unique(clean_ws)[1:]:\n",
    "\n",
    "    # (iii) Erode the cell's mask by 1 pixel\n",
    "    # Hint: 'ndi.binary_erode'\n",
    "    cell_mask = clean_ws==cell_ID\n",
    "    eroded_cell_mask = ndi.binary_erosion(cell_mask, iterations=1) # Increase iterations to make boundary wider!\n",
    "    \n",
    "    # (iv) Create the cell edge mask\n",
    "    # Hint: 'np.logical_xor'\n",
    "    edge_mask = np.logical_xor(cell_mask, eroded_cell_mask)\n",
    "    \n",
    "    # (v) Add the cell edge mask to the empty array generated above, labeling it with the cell's ID\n",
    "    edges[edge_mask] = cell_ID"
   ]
  },
  {
   "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",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(np.zeros_like(edges)[300:500, 300:500], cmap='gray', vmin=0, vmax=1)  # Simple black background\n",
    "plt.imshow(np.ma.array(edges, mask=edges==0)[300:500, 300:500], interpolation='none', cmap='prism')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting Quantitative Measurements <a id=measure></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</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 simplicty, we will here stick with a dictionary of lists."
   ]
  },
  {
   "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",
    "results = {\"cell_id\"      : [],\n",
    "           \"int_mean\"     : [],\n",
    "           \"int_mem_mean\" : [],\n",
    "           \"cell_area\"    : [],\n",
    "           \"cell_edge\"    : []}\n",
    "\n",
    "# Solution note: the spacing between the strings and colons doesn't matter for the code's\n",
    "# execution. It is used solely to make the code more readable!"
   ]
  },
  {
   "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",
    "# Iterate over cell IDs\n",
    "for cell_id in np.unique(clean_ws)[1:]:\n",
    "\n",
    "    # Mask the current cell and cell edge\n",
    "    cell_mask = clean_ws==cell_id\n",
    "    edge_mask = edges==cell_id\n",
    "    \n",
    "    # Get the measurements\n",
    "    results[\"cell_id\"].append(cell_id)\n",
    "    results[\"int_mean\"].append(np.mean(img[cell_mask]))\n",
    "    results[\"int_mem_mean\"].append(np.mean(img[edge_mask]))\n",
    "    results[\"cell_area\"].append(np.sum(cell_mask))\n",
    "    results[\"cell_edge\"].append(np.sum(edge_mask))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Print the results and check that they make sense\n",
    "\n",
    "for key in results.keys(): \n",
    "    print(key + \":\", results[key][:5], '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Analysis & Visualisation <a id=analysis></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Analyze and plot the extracted data in a variety of ways."
   ]
  },
  {
   "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",
    "# Custom function for nice printing of summary statistics.\n",
    "# Note the use of format strings for nice number padding.\n",
    "def print_summary(data):\n",
    "    print( \"  Mean:    {:7.2f}\".format(np.mean(data))   )\n",
    "    print( \"  Stdev:   {:7.2f}\".format(np.std(data))    )\n",
    "    print( \"  Max:     {:7.2f}\".format(np.max(data))    )\n",
    "    print( \"  Min:     {:7.2f}\".format(np.min(data))    )\n",
    "    print( \"  Median:  {:7.2f}\".format(np.median(data)) )\n",
    "\n",
    "# Calling the custom function for each dataset\n",
    "for key in results.keys():\n",
    "    print( '\\n'+key )\n",
    "    print_summary(results[key])\n",
    "    \n",
    "# There are also pre-made functions to get summary statistics,\n",
    "# for example 'scipy.stats.describe'.\n",
    "from scipy.stats import describe\n",
    "stat_summary = describe(results['int_mean'])\n",
    "\n",
    "print( '\\nscipy.stats.describe of int_mean' )\n",
    "for key in stat_summary._asdict().keys():\n",
    "    print( ' ', key+': ', stat_summary._asdict()[key] )"
   ]
  },
  {
   "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",
    "plt.figure(figsize=(3,6))\n",
    "plt.boxplot([results['int_mean'], results['int_mem_mean']], \n",
    "            labels=['int_mean', 'int_mem_mean'],\n",
    "            widths=0.6, notch=True)\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "\n",
    "plt.figure(figsize=(8,5))\n",
    "plt.scatter(results[\"cell_area\"], results[\"cell_edge\"],\n",
    "            edgecolor='k', s=30, alpha=0.5)\n",
    "plt.xlabel(\"cell area [pxl^2]\")\n",
    "plt.ylabel(\"cell edge length [pxl]\")\n",
    "\n",
    "# BONUS: Do you understand why you are seeing the pattern this produces? \n",
    "# ->> The curve reflects how circumference scales with area!\n",
    "\n",
    "# Can you generate a 'null model' curve that assumes all cells to be circular? \n",
    "cell_area_range = np.linspace(min(results[\"cell_area\"]), max(results[\"cell_area\"]), num=100)\n",
    "circle_circumference = 2 * np.pi * np.sqrt( cell_area_range / np.pi )\n",
    "plt.plot(cell_area_range, circle_circumference, color='r', alpha=0.8)\n",
    "plt.legend(['circles', 'data'], loc=2, fontsize=10)\n",
    "plt.show()\n",
    "\n",
    "# What is the result? Do you notice something odd about it? What could be the reason for\n",
    "# this and how could it be fixed?\n",
    "# ->> In general, the cells don't deviate all that much from the circular case.\n",
    "# ->> Strangely, some cells have a smaller outline than the circumference of a circle\n",
    "#     of equivalent area. This is mathematically impossible.\n",
    "# ->> A possible reason could be that the measures are taken in pixels, which leads\n",
    "#     to a so-called discretization error. It could be fixed by \"meshing\" the cell\n",
    "#     outline and interpolating a more accurate measurement of circumference."
   ]
  },
  {
   "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",
    "# Compute linear fit\n",
    "from scipy.stats import linregress\n",
    "linfit = linregress(results[\"cell_area\"], results[\"int_mem_mean\"])\n",
    "\n",
    "# Print all the results\n",
    "linprops = ['slope','interc','rval','pval','stderr']\n",
    "for index,prop in enumerate(linprops):\n",
    "    print( prop, '\\t', '{:4.2e}'.format(linfit[index]) )"
   ]
  },
  {
   "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",
    "#       ->> This p-value only means that, given a linear fit through this data, the slope of the\n",
    "#           fit is very unlikely to be zero. However, it does not make a statement on whether or\n",
    "#           not it makes sense to use a linear fit in the first place. Looking at the scatterplot\n",
    "#           below or at the correlation coefficient r, it is clear that a linear fit on this data\n",
    "#           is not meaningful.\n",
    "#       ->> Note also: With single-cell approaches, we quickly get to a large number of data points. \n",
    "#           This makes hypothesis testing in general less useful, as p-values tend to become very\n",
    "#           small even if the null hypothesis holds. It makes sense to instead report effect sizes.\n",
    "#           This is a tricky topic but well worth reading up on.\n",
    "#\n",
    "#   2) Could there be artifacts in our segmentation that bias this analysis?\n",
    "#       ->> Oversegmentation is an important source of bias here. If a cell is oversegmented,\n",
    "#           it will be considered as two or three cells. These will naturally have a lower\n",
    "#           cell area and will naturally have a lower membrane intensity because some of their\n",
    "#           edges are actually not on membranes. In other words, they will fall into the bottom\n",
    "#           left of the plot, distorting the data.\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",
    "#   ->> The max and min values in the data are a good choice.\n",
    "x_vals = [min(results[\"cell_area\"]), max(results[\"cell_area\"])]\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",
    "y_vals = [linfit[0] * x_vals[0] + linfit[1], linfit[0] * x_vals[1] + linfit[1]]\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",
    "plt.figure(figsize=(8,5))\n",
    "plt.scatter(results[\"cell_area\"], results[\"int_mem_mean\"], \n",
    "            edgecolor='k', s=30, alpha=0.5)\n",
    "plt.plot(x_vals, y_vals, color='red', lw=2, alpha=0.8)\n",
    "\n",
    "# Use 'plt.legend' to add information about the line to the plot.\n",
    "plt.legend([\"linear fit, Rsq={:4.2e}\".format(linfit[2]**2.0)], frameon=False, loc=4)\n",
    "\n",
    "# Label the plot and finally show it with 'plt.show'.\n",
    "plt.xlabel(\"cell area [pxl]\")\n",
    "plt.ylabel(\"Mean membrane intensity [a.u.]\")\n",
    "plt.title(\"Scatterplot with linear fit\")\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "areas_8bit = np.array(results[\"cell_area\"]) / max(results[\"cell_area\"]) * 255\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",
    "area_map = np.zeros_like(clean_ws, dtype=np.uint8) \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",
    "for index, cell_id in enumerate(results[\"cell_id\"]):\n",
    "    \n",
    "    # Mask the current cell and assign the cell's (re-scaled) area value to the cell's pixels.\n",
    "    area_map[clean_ws==cell_id] = areas_8bit[index]\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",
    "\n",
    "# Mask of outliers (the largest and smallest 5% of all cells)\n",
    "outlier_mask = np.logical_or(area_map > np.percentile(areas_8bit, 95),\n",
    "                             area_map < np.percentile(areas_8bit, 5))\n",
    "\n",
    "# Mask of all regions to leave blank (outliers + image boundary cells)\n",
    "full_mask = np.logical_or(area_map==0, outlier_mask)\n",
    "\n",
    "# Create the plot\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(img_smooth, interpolation='none', cmap='gray')\n",
    "plt.imshow(np.ma.array(area_map, mask=full_mask),\n",
    "           interpolation='none', cmap='viridis', alpha=0.6)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing Output to Files <a id=write></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Write the generated data into a variety of different output files."
   ]
  },
  {
   "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",
    "from skimage.io import imsave\n",
    "imsave(\"example_cells_1_edges.tif\", edges.astype(np.uint16))"
   ]
  },
  {
   "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",
    "# Create plot (but don't show)\n",
    "plt.scatter(results[\"cell_area\"], results[\"int_mem_mean\"], \n",
    "            edgecolor='k', s=30, alpha=0.5)\n",
    "plt.plot(x_vals, y_vals, color='red', lw=2, alpha=0.8)\n",
    "plt.legend([\"linear fit, Rsq={:4.2e}\".format(linfit[2]**2.0)], frameon=False, loc=4)\n",
    "plt.xlabel(\"cell area [pxl]\")\n",
    "plt.ylabel(\"Mean membrane intensity [a.u.]\")\n",
    "plt.title(\"Scatterplot with linear fit\")\n",
    "\n",
    "# Save as png and pdf\n",
    "plt.savefig('example_cells_1_scatterFit.png')\n",
    "plt.savefig('example_cells_1_scatterFit.pdf')\n",
    "plt.clf()  # Clear the figure buffer"
   ]
  },
  {
   "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",
    "np.save(\"example_cells_1_seg\", clean_ws)  # Save\n",
    "seg = np.load(\"example_cells_1_seg.npy\")  # Load\n",
    "print(clean_ws.shape, seg.shape)"
   ]
  },
  {
   "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 mode ('w'). It's best to do \n",
    "# so using the 'with'-statement (context manager) to make sure that the file object will be closed\n",
    "# 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",
    "import pickle\n",
    "with open('example_cells_1_results.pkl','wb') as outfile:\n",
    "    pickle.dump(results, outfile)\n",
    "\n",
    "# Note: Pickled files can be re-loaded again as follows:\n",
    "with open('example_cells_1_results.pkl', 'rb') as infile:\n",
    "    results_reloaded = pickle.load(infile)\n",
    "    print(results_reloaded.keys())"
   ]
  },
  {
   "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",
    "with open('example_cells_1_results.txt','w') as outfile:\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",
    "    header_string = '\\t'.join(results.keys()) + '\\n'\n",
    "    outfile.write(header_string)\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",
    "    for index in range(len(results['cell_id'])):\n",
    "        data_string = '\\t'.join([str(results[key][index]) for key in results.keys()]) + '\\n'\n",
    "        outfile.write(data_string)\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": [
    "#### <font color='orange'>Exercise Solution</font>\n",
    "\n",
    "Convert the pipeline to a single function in a python script.\n",
    "\n",
    "Import that function and write the code necessary to iterate over all input files in a folder, apply the function to each file, and collect the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solution Note\n",
    "\n",
    "The converted version of the pipeline can be found in `batch_processing_solution.py`. \n",
    "\n",
    "Note that the most of the exercise comments have been removed so it doesn't look too cluttered. However, this level of clean-up is probably a bit extreme; it is generally recommended to retain at least basic comments on the purpose of each code block. \n",
    "\n",
    "Also note that a doc string was added to the function definition (the string designated with three `\"\"\"` at the start and end, directly under the function definition). This is a very useful reference, since it is automatically recognized as a help message by Jupyter notebook (and other IDEs), so you can easily double-check what an imported function does, for example by typing `run_pipeline?` in a code cell."
   ]
  },
  {
   "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",
    "from batch_processing_solution import run_pipeline\n",
    "pip_seg, pip_results = run_pipeline(r\"example_data\", r'example_cells_1.tif')\n",
    "\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.imshow(np.zeros_like(pip_seg), interpolation='none', cmap='gray', vmax=1)  # Simple black background\n",
    "plt.imshow(np.ma.array(pip_seg, mask=pip_seg==0), interpolation='none', cmap='prism')\n",
    "plt.show()"
   ]
  },
  {
   "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",
    "# Get all files\n",
    "dirpath = r\"example_data\"\n",
    "from os import listdir\n",
    "filelist = listdir(dirpath)\n",
    "\n",
    "# Filter for target files: simple option\n",
    "# Note that this will use ALL files with a .tif ending, which in some circumstances\n",
    "# may include files that are not supposed to be used!\n",
    "target_files = []\n",
    "for fname in filelist:\n",
    "    if fname.endswith('.tif'):\n",
    "        target_files.append(fname)\n",
    "print(target_files)\n",
    "\n",
    "# Filter for target files: advanced option using regex and a list comprehension\n",
    "import re\n",
    "target_pattern = re.compile(\"^example_cells_\\d+\\.tif$\")\n",
    "target_files = [fname for fname in filelist if target_pattern.match(fname)]\n",
    "print(target_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (iii) Iterate over the relevant 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",
    "# Initialize empty dictionaries\n",
    "all_seg = {}\n",
    "all_results = {}\n",
    "\n",
    "# Iterate over files and run pipeline for each\n",
    "for fname in target_files:\n",
    "    pip_seg, pip_results = run_pipeline(dirpath, fname)\n",
    "    all_seg[fname] = pip_seg\n",
    "    all_results[fname] = pip_results"
   ]
  },
  {
   "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",
    "plt.figure(figsize=(8,5))\n",
    "\n",
    "colors = ['blue','red']\n",
    "for key, color in zip(sorted(all_results.keys()), colors):\n",
    "    plt.scatter(all_results[key][\"cell_area\"], all_results[key][\"cell_edge\"],\n",
    "                edgecolor='k', c=color, s=30, alpha=0.5, label=key)\n",
    "    \n",
    "plt.legend()\n",
    "plt.xlabel(\"cell area [pxl^2]\")\n",
    "plt.ylabel(\"cell edge length [pxl]\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='orange'>*This is the end of 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 [conda env:py39_tutorial_test]",
   "language": "python",
   "name": "conda-env-py39_tutorial_test-py"
  },
  "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}