1 - Handling Image Data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1: Handling Image Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, we will learn how to manipulate image data as arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import imageio.v2 as imageio\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "\n",
    "plt.rcParams['figure.dpi'] = 200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import skimage as sk\n",
    "import skimage.filters as skf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Loading and visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Throughout this practical, we will be using two image sets from the Broad Bioimage Benchmark Collection (BBBC, https://bbbc.broadinstitute.org/):\n",
    "- BBBC010 (https://bbbc.broadinstitute.org/BBBC010/), a dataset of brightfield microscopy images of a C. elegans live/dead assay;\n",
    "- BBBC020 (https://bbbc.broadinstitute.org/BBBC020/), a dataset of fluorescence microscopy images of the cell membrane and nucleus of murine bone-marrow derived macrophages."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Load \n",
    "\n",
    "Load sample images from each datasets by running the lines below. \n",
    "\n",
    "*There exists plenty of different libraries to load image data, but here we will rely on imageio (https://pypi.org/project/imageio/) as it can handle a wide range of formats that are relevant to microscopy.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BBBC010\n",
    "# Live worms example\n",
    "bbbc010_live = imageio.imread('data/Part_1/BBBC010/live.tif')\n",
    "# Dead worms example\n",
    "bbbc010_dead = imageio.imread('data/Part_1/BBBC010/dead.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BBBC020\n",
    "# Cell membranes example\n",
    "bbbc020_cell = imageio.imread('data/Part_1/BBBC020/membranes.tif')\n",
    "# Nuclei example\n",
    "bbbc020_nuclei = imageio.imread('data/Part_1/BBBC020/nuclei.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Visualize\n",
    "\n",
    "Visualize the images with matplotlib (https://matplotlib.org/) by running the lines below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BBBC010\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "\n",
    "axes[0].imshow(bbbc010_live)\n",
    "axes[0].axis('off')\n",
    "axes[0].set_title('Live C. elegans')\n",
    "\n",
    "axes[1].imshow(bbbc010_dead)\n",
    "axes[1].axis('off')\n",
    "axes[1].set_title('Dead C. elegans')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# BBBC020\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "\n",
    "axes[0].imshow(bbbc020_cell)\n",
    "axes[0].axis('off')\n",
    "axes[0].set_title('Cell membrane marker')\n",
    "\n",
    "axes[1].imshow(bbbc020_nuclei)\n",
    "axes[1].axis('off')\n",
    "axes[1].set_title('Nuclear marker')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Colormaps\n",
    "\n",
    "As the image data we are working with are grayscale, what we see in 1.2 is pseudo-color. We can change that for any other colormap that matplotlib offers (see https://matplotlib.org/stable/tutorials/colors/colormaps.html). Run the code below to first visualize the BBBC020 images first in grayscale, and then adapt it to use your favourite colormap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grayscale\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "\n",
    "axes[0].imshow(bbbc020_cell, cmap='gray')\n",
    "axes[0].axis('off')\n",
    "axes[0].set_title('Cell membrane marker')\n",
    "\n",
    "axes[1].imshow(bbbc020_nuclei, cmap='gray')\n",
    "axes[1].axis('off')\n",
    "axes[1].set_title('Nuclear marker')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your favourite colormap\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Cropping\n",
    "\n",
    "Since images are stored as arrays, we can vizualize regions of interest by specifying a range of image coordinates to be displayed. Run the code lines below to visualize a crop of one of the BBBC010 images. Can you modify it to display a different part of the image? \n",
    "\n",
    "*Bonus challenge*: can you adapt the code so that the origin of the crop and its size are expressed relatively to the image dimensions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initial crop\n",
    "origin = (100,250)\n",
    "size = (200,200)\n",
    "\n",
    "plt.imshow(bbbc010_live[origin[0]:origin[0]+size[0],origin[1]:origin[1]+size[1]], cmap='gray')\n",
    "plt.axis('off')\n",
    "plt.title('BBBC010 - RoI')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Crop around different area\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bonus challenge: crop in relative coordinates\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5 Roi editing\n",
    "Although we used it only to visualize images up to now, matplotlib is a general plotting library that can do much more. One relevant feature in the context of this tutorial is the ability to combine classical plots with images. Run the code below to see one such example in which a red rectangle is overlaid to surround a region of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.patches as patches\n",
    "\n",
    "# Crop from 1.3\n",
    "origin = (100,250)\n",
    "size = (200,200)\n",
    "\n",
    "# Note that the x and y axes are swapped for images and plots\n",
    "rect = patches.Rectangle(origin[::-1], size[0], size[1], linewidth=1, edgecolor='r', facecolor='none')\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "\n",
    "axes[0].imshow(bbbc010_live, cmap='gray')\n",
    "axes[0].axis('off')\n",
    "axes[0].add_patch(rect)\n",
    "axes[0].set_title('Highlighted RoI')\n",
    "\n",
    "axes[1].imshow(bbbc010_live[origin[0]:origin[0]+size[0],origin[1]:origin[1]+size[1]], cmap='gray')\n",
    "axes[1].axis('off')\n",
    "axes[1].set_title('RoI')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Image modifications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know how to load and visualize image data, let's have a look at basic operations we can do to modify them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Masking\n",
    "\n",
    "In the same way that we were able to visualize a crop of the image, we can modify selected portion of the data. By adapting the code from 1.4, \"erase\" a patch of a copy of one of the images we are working with by setting to zero all pixels in a range of coordinates. This operation is usually referred to as *masking*. To check the result, visualize the modified image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do not forget to work on a copy of the image!\n",
    "original_img = bbbc010_live\n",
    "masked_img = copy.deepcopy(original_img)\n",
    "\n",
    "# Region to be masked\n",
    "(width, height) = masked_img.shape\n",
    "origin = (int(0.5*width),int(0.5*height))\n",
    "size = (int(0.2*width),int(0.2*height))\n",
    "\n",
    "# Masking\n",
    "masked_img[origin[0]:origin[0]+size[0],origin[1]:origin[1]+size[1]] = 0\n",
    "\n",
    "# Visualize result\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Inverting\n",
    "\n",
    "Further to setting pixels to whichever constant value we chose (as we did when masking in 2.1), we can also apply modifications that depend on the original image values. A simple example is *image inversion*, also called *image negative*. The grayscale images we are working with are all in the [0,255] range. To compute the inverted (or negative) version of the image of your choice, you will therefore need to swap values so that what used to be 0 is mapped to 255, what used to 255 is mapped to 0, and so on for all values in between. Implement image inversion and vizualise the result below.\n",
    "\n",
    "*Bonus*: if you only want the negative of the image for esthetic purposes (i.e., it looks better, but this is not needed for further processing), there is no need to modify the data directly. Instead, it is best to play on the visualization directly with colormaps. By adapting what you did in 1.3, display the inverse of the original image without modifying its data.\n",
    "\n",
    "*Bonus challenge*: the [0,255] range is the realm of 8-bit images, but there are many more formats out there. If you are up for the challenge, adapt your code to invert an image from each dataset in its raw format (16-bit and RGB). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "original_img = bbbc020_cell\n",
    "inverted_img = copy.deepcopy(original_img)\n",
    "\n",
    "# Inversion\n",
    "##### Add you code here ##### \n",
    "        \n",
    "# Visualize result\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "\n",
    "axes[0].imshow(original_img, cmap='gray')\n",
    "axes[0].axis('off')\n",
    "axes[0].set_title('Original image')\n",
    "\n",
    "axes[1].imshow(inverted_img, cmap='gray')\n",
    "axes[1].axis('off')\n",
    "axes[1].set_title('Inverted image')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bonus: play on colormap (no image modification)\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bonus challenge: invert a 16-bit image\n",
    "original_img = imageio.imread('data/Part_1/BBBC010/dead_raw.tif')\n",
    "inverted_img = copy.deepcopy(original_img)\n",
    "\n",
    "# Inversion\n",
    "##### Add you code here ##### \n",
    "        \n",
    "# Visualize result\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "\n",
    "axes[0].imshow(original_img, cmap='gray')\n",
    "axes[0].axis('off')\n",
    "axes[0].set_title('Original image')\n",
    "\n",
    "axes[1].imshow(inverted_img, cmap='gray')\n",
    "axes[1].axis('off')\n",
    "axes[1].set_title('Inverted image')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bonus challenge: invert an RGB image\n",
    "original_img = imageio.imread('data/Part_1/BBBC020/membranes_raw.tif')\n",
    "inverted_img = copy.deepcopy(original_img)\n",
    "\n",
    "# Inversion\n",
    "##### Add you code here ##### \n",
    "        \n",
    "# Visualize result\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "\n",
    "axes[0].imshow(original_img)\n",
    "axes[0].axis('off')\n",
    "axes[0].set_title('Original image')\n",
    "\n",
    "axes[1].imshow(inverted_img)\n",
    "axes[1].axis('off')\n",
    "axes[1].set_title('Inverted image')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Filtering\n",
    "\n",
    "A general way to modify an image based on its original values is *filtering*. Filtering consists of applying a so-called convolution kernel to an input image (https://en.wikipedia.org/wiki/Kernel_(image_processing). Although the technical details of image filtering are beyond the scope of this tutorial, we can get a sense of what filtering means by considering some examples provided below. Run the code below to see how 1) a Gaussian kernel, 2) the Prewitt filter, and 3) the Frangi filter process the input image. What are their effects? For the Gaussian kernel and Frangi filters, try to modify the value of the parameter called sigma. How does it affect the result? For the Gaussian kernel, do you intuitively understand why?\n",
    "\n",
    "*The scikit-image library offers many more filters for you to try out and play with if you are interested: https://scikit-image.org/docs/dev/api/skimage.filters.html.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "original_img = bbbc010_live\n",
    "\n",
    "# Important step: make data floating-point-valued for further processing\n",
    "original_img = sk.img_as_float(original_img)\n",
    "\n",
    "# Apply Gaussian kernel\n",
    "filtered_img_gaussian = skf.gaussian(original_img, sigma=10)\n",
    "\n",
    "# Apply Prewitt filter\n",
    "filtered_img_prewitt = skf.prewitt(original_img)\n",
    "\n",
    "# Apply Frangi filter\n",
    "filtered_img_frangi = skf.frangi(original_img, sigmas=range(1, 10, 2))\n",
    "\n",
    "# Visualize result\n",
    "fig, axes = plt.subplots(nrows=1, ncols=4)\n",
    "\n",
    "axes[0].imshow(original_img)\n",
    "axes[0].axis('off')\n",
    "axes[0].set_title('Original image', fontsize=5)\n",
    "\n",
    "axes[1].imshow(filtered_img_gaussian)\n",
    "axes[1].axis('off')\n",
    "axes[1].set_title('Filtered image (Gaussian)', fontsize=5)\n",
    "\n",
    "axes[2].imshow(filtered_img_prewitt)\n",
    "axes[2].axis('off')\n",
    "axes[2].set_title('Filtered image (Prewitt)', fontsize=5)\n",
    "\n",
    "axes[3].imshow(filtered_img_frangi)\n",
    "axes[3].axis('off')\n",
    "axes[3].set_title('Filtered image (Frangi)', fontsize=5)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Modify the parameter called sigma and look at the result\n",
    "##### Add you code here ##### "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Image processing\n",
    "\n",
    "Filtering is at the core of many image processing methods. Run the code below to see what happens when taking the difference between one of the BBBC010 images and its Gaussian-blurred version. What do you notice? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "original_img = bbbc010_live\n",
    "\n",
    "# Important step: make data floating-point-valued for further processing\n",
    "original_img = sk.img_as_float(original_img)\n",
    "\n",
    "# Apply Gaussian kernel with large-enough sigma to blur the objects\n",
    "filtered_img = skf.gaussian(original_img, sigma=20)\n",
    "\n",
    "# Compute the difference between the original and filtered image\n",
    "diff_img = original_img - filtered_img\n",
    "\n",
    "# Visualize result\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3)\n",
    "\n",
    "axes[0].imshow(original_img)\n",
    "axes[0].axis('off')\n",
    "axes[0].set_title('Original image', fontsize=5)\n",
    "\n",
    "axes[1].imshow(filtered_img)\n",
    "axes[1].axis('off')\n",
    "axes[1].set_title('Filtered image', fontsize=5)\n",
    "\n",
    "axes[2].imshow(diff_img)\n",
    "axes[2].axis('off')\n",
    "axes[2].set_title('Difference image', fontsize=5) # What do you notice?\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Export"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The images we created or modified in this exercise are Python array that only live within this notebook and will disappear whenever it stops running. In order to save them for later, they need to be exported into files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 tiff export\n",
    "\n",
    "Run the code below to save the image from 2.1 as a tiff file. You can then adapt it to save the other images you generated, or try out different file formats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imageio.imwrite('data/Part_1/BBBC010/my_masked_image.tif', masked_img)"
   ]
  }
 ],
 "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}