Module2-ImageProcessing/M2A_Basic_Image_Processing.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zKSdkHsA_Wv_"
      },
      "source": [
        "Introduction to Image Processing - Basic Image Manipulations\n",
        "<html>\n",
        "    <summary></summary>\n",
        "         <div> <p></p> </div>\n",
        "         <div style=\"font-size: 20px; width: 800px;\"> \n",
        "              <h1>\n",
        "               <left>Basic Image Manipulation in Python</left>\n",
        "              </h1>\n",
        "              <p><left>============================================================================</left> </p>\n",
        "<pre>Course: UQ-Bio 2024\n",
        "Instructor: Luis Aguilera\n",
        "Authors: Luis Aguilera, Will Raymond, Brian Munsky\n",
        "Contact Info: luis.aguilera@colostate.edu\n",
        "</pre>\n",
        "         </div>\n",
        "    </p>\n",
        "\n",
        "</html>\n",
        "\n",
        "<details>\n",
        "  <summary>Copyright info</summary>\n",
        "\n",
        "```\n",
        "Copyright 2024 Brian Munsky\n",
        "\n",
        "Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:\n",
        "\n",
        "1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.\n",
        "\n",
        "2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.\n",
        "\n",
        "3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.\n",
        "\n",
        "THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\n",
        "```\n",
        "<details>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YmGRxtCzqWqs"
      },
      "source": [
        "# Abstract\n",
        "\n",
        "This notebook provides a list of procedures to analyze microscope images. The notebook describes what a digital image is, and how to extract relevant information from the image. At the end of the tutorial, you should have obtained the computational skills to implement the following list of objectives independently.\n",
        "\n",
        "## List of objectives\n",
        "\n",
        "\n",
        "1. Load the **python modules** commonly used to work with microscope data.\n",
        "2. Understand and explain what a **digital image** is in terms of matrices and tensors.\n",
        "3. Understand and explain what a **monochromatic image** is and a **color image** is.\n",
        "4. Select and **slice the dimensions** in a sequence of microscope images.\n",
        "5. Apply different **filters** to remove noise from an image using linear algebra operations.\n",
        "6. Perform **basic mathematic operations** involved in image processing, including rotation, translation, and scaling."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MabWJ6vkA5Dh"
      },
      "source": [
        "# Importing Libraries"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_XltNqo0PlDt"
      },
      "source": [
        "The following lines of code import and install some libraries. For more information, look at the library name on the  Python Package Index [(PyPI)](https://pypi.org/)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "oUp21jU285ef"
      },
      "outputs": [],
      "source": [
        "# Loading libraries for the course\n",
        "import random\n",
        "import matplotlib.pyplot as plt             # Library used for plotting\n",
        "from matplotlib.patches import Rectangle    # Module to plot a rectangle in the image\n",
        "import urllib.request                       # Library to download data\n",
        "import numpy as np                          # Library for array manipulation\n",
        "import seaborn as sn                        # Library for advanced plotting\n",
        "import pandas as pd                         # Library to manipulate data frames\n",
        "import tifffile                             # Library to store numpy arrays as TIFF\n",
        "import pathlib                              # Library to work with file paths\n",
        "from pathlib import Path                    # Library to work with file paths\n",
        "import skimage                              # Library for image manipulation. scikit-image\n",
        "from skimage.io import imread               # Module from skimage\n",
        "from matplotlib import animation            # Module to plot animations"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "K3jGYDf50dZW"
      },
      "outputs": [],
      "source": [
        "students = ['Daria', 'Wiktoria','Jinyi','Wei','Pilar','Camila','Jayden','Humberto',\n",
        "'Jakini','Varvara','Kanniese','Kristi','Eliud','Nam','Lakshay','Margaret',\n",
        "'Grace','Apurv','Yinfei','Samantha','Yingli','Rena','Emily','Pegah','Al-Bashir']"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "S7-k93i6FR4w"
      },
      "outputs": [],
      "source": [
        "random.choice(students)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Nxsu81OwX90z"
      },
      "source": [
        "A library is a collection of code to perform a specific task. For example, check the [scikit-image library webpage](https://scikit-image.org).\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H6LteMHOZhYL"
      },
      "source": [
        "Python is the most popular programming language in data science, and the most popular programming language in CS [tiobe-index](https://www.tiobe.com/tiobe-index/).\n",
        "\n",
        "\n",
        "# Python is an ecosistem with a comprensive list of well-mantanined libraries."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "APoMxFVJPuhj"
      },
      "source": [
        "<img src= https://numpy.org/images/content_images/ds-landscape.png alt=\"drawing\" width=\"1200\"/>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xF80TYlZarBB"
      },
      "source": [
        "# Working with images in python"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ms-rgCqfiyPE"
      },
      "source": [
        "![alt text](FigsA/Module_1_1/Slide2.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q5WNA54G_Sdd"
      },
      "source": [
        "Image modified from: Gonzalez, Rafael & Faisal, Zahraa. (2019). Digital Image Processing Second Edition."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9rw73zKSSzxS"
      },
      "source": [
        "## Downloading, opening, and visualizing images\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CryOMGJyTJKW"
      },
      "outputs": [],
      "source": [
        "# Downloading the image from figshare SupFig1c_BG_MAX_Cell04.tif\n",
        "urls = ['https://ndownloader.figshare.com/files/26751209','https://ndownloader.figshare.com/files/26751203','https://ndownloader.figshare.com/files/26751212','https://ndownloader.figshare.com/files/26751218']\n",
        "print('Downloading file...')\n",
        "urllib.request.urlretrieve(urls[1], './image_cell.tif') "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "H9iaLiiUZ_l9"
      },
      "outputs": [],
      "source": [
        "# Importing the image as variable img\n",
        "figName = './image_cell.tif'\n",
        "img = imread(str(figName))\n",
        "type(img)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "img.max()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "During this course we have been constantly using [NumPy](https://numpy.org), [Matplotlib](https://matplotlib.org),\n",
        "[Scipy](https://scipy.org), [Pandas](https://pandas.pydata.org), [Scikit-Image](https://scikit-image.org), [TensorFlow](https://www.tensorflow.org), etc."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# If an image is an array, how can we find the size of that array?\n",
        "random.choice(students)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# What should we do if we want to know the size of the image?\n",
        "np.shape(img)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "img.shape"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "If you need a refresher, please checkgo back over the tutorials in Module 1 for basic Python and numpy. Also check this link for [cheatsheets](https://www.datacamp.com/resources/data-science-and-analytics-cheatsheets).\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dw1h-pNreUMt"
      },
      "source": [
        "## Understanding digital images"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n_7ggMeAjrUZ"
      },
      "source": [
        "![alt text](FigsA/Module_1_1/Slide3.png)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5dUYyyKvSQbw"
      },
      "source": [
        "**Displaying** a section of the image. Notice that an image is only a matrix of numbers."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EbpGVc5afAH9"
      },
      "source": [
        "What is a digital image?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9UZmW3AXaQpg"
      },
      "outputs": [],
      "source": [
        "# What is img?\n",
        "print('Image type =', type(img))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2pu8D7shV8Vu"
      },
      "source": [
        "It is essential to understand that images in Python are stored as NumPy objects.  Need a refresher?  Check out the totorials in Module 1."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "r7RLPW07e9_8"
      },
      "source": [
        "\n",
        "What is the shape of the image?\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FhXrbP-6e7dk"
      },
      "outputs": [],
      "source": [
        "print('Image dimensions, Shape =', img.shape )   #[T,Y,X,C]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "IWcnOKXt6Eaj"
      },
      "outputs": [],
      "source": [
        "# Attributes in img object\n",
        "print([x for x in dir(img) if '__' not in x] ) #what functions (excluding internal functions, denoted by __xx__) does a list object have?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "img.max()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HFw7CDHXfJsM"
      },
      "source": [
        "**Displaying** a section of the image. Notice that an image is only a matrix of numbers."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9v0EcJ84QQ0g"
      },
      "source": [
        "From the [image's publication](https://www.biorxiv.org/content/10.1101/2020.04.03.024414v2) we can obtain the **metadata**. Indicating that the following information:\n",
        "\n",
        "Dimension  | Meaning |  Value\n",
        "---------|---------- |----------\n",
        "0   | Time        | 35 (frames)\n",
        "1   | Y-dimension | 512 pixels\n",
        "2   | X-dimension | 512 pixels\n",
        "3   | Color       | 3 color image (R,G,B)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2ElEneHLQiON"
      },
      "source": [
        "![alt text](FigsA/Module_1_1/Slide4.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lUitS42PQnse"
      },
      "source": [
        "High-order tensor are also refered as Hyperstacks in software like [imagej](https://imagej.net/software/fiji/)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kViaZsEhQ2fL"
      },
      "source": [
        "# Slicing a high-order tensor to display it as a 2D image."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7sdOtVMERCS0"
      },
      "source": [
        "The image is a numpy array with dimensions:\n",
        "```\n",
        "image[time, y-axis, x-axis, colors]\n",
        "```\n",
        "To display our tensor we can select a 2D subsection.\n",
        "\n",
        "```\n",
        "image[selected_time, : , : , selected_color]\n",
        "```\n",
        "\n",
        "\n",
        "* Notice that when we slice [start: end(not including this value) :step]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "wSpcwfvMTdv1"
      },
      "outputs": [],
      "source": [
        "# How would we show just a grayscale image corresponding to the 3rd channel, at the \n",
        "# second time point?\n",
        "random.choice(students)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "IRlgXhLCqePQ"
      },
      "outputs": [],
      "source": [
        "# Recall, the image shape is:\n",
        "print('Image dimensions, Shape =', img.shape )   #[T,Y,X,C]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZzeF30uxxCcO"
      },
      "outputs": [],
      "source": [
        "# Plotting an image\n",
        "plt.figure(figsize=(7,7))\n",
        "selected_frame = 0\n",
        "selected_color_channel = 2\n",
        "plt.imshow( img[selected_frame,100:300,100:300,selected_color_channel], cmap='gray') # Notice that only a time point and a color is plotted\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# find the maximum value of the image\n",
        "img.max()\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZFiHoZwdc2IZ"
      },
      "outputs": [],
      "source": [
        "# Notice the difference between plotting an image (using imshow) and ploting a time course\n",
        "# (using plot). The former is used to plot an image, the latter is used to plot a line.\n",
        "plt.plot([1,2,3,4])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "N51bwNQkIPC8"
      },
      "outputs": [],
      "source": [
        "# How would I find the brightest pixel in the image?\n",
        "random.choice(students)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Y_yUuV207STB"
      },
      "source": [
        "#**Intensity** values in the image"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1ph7AsyjgYVj"
      },
      "outputs": [],
      "source": [
        "# Minimum and maximum intensity values on the image\n",
        "max_intensity_value = np.max(img)\n",
        "min_intensity_value = np.min(img)\n",
        "quant_intensity_value = np.quantile(img, [0.9,0.8,0.5])    # 0.9 is equivalent to 90 percentile\n",
        "\n",
        "print('Maximum intensity : ', max_intensity_value)\n",
        "print('Minimum intensity : ', min_intensity_value)\n",
        "print('Quantile intensity: ', quant_intensity_value )\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ppQ0yOqbf-Gw"
      },
      "source": [
        "To understand the parameters that you need to pass to a method use the ```help``` function\n",
        "```\n",
        "help(method)\n",
        "```\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "37CZnceggM4c"
      },
      "outputs": [],
      "source": [
        "#help(np.quantile)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u2HVf-6R7WPz"
      },
      "source": [
        "Intensity distribution in the image"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PN5SazcD7c52"
      },
      "outputs": [],
      "source": [
        "# Plotting the intensity distribution for a specific time point and an specific channel\n",
        "plt.figure(figsize=(3,2))\n",
        "plt.hist(img[:,:,:,2].flatten(), bins=80,color='orangered')\n",
        "plt.xlabel('Value')\n",
        "plt.ylabel('Frequency')\n",
        "plt.title('Intensity Histogram')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qOGE49vGg51F"
      },
      "source": [
        "Summary of image properties:\n",
        "\n",
        "* 4 dimensional tensor [T, Y, X, C].\n",
        "* Numpy array\n",
        "* Intensity range (0, 7965)\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "sXp78njsJ_Ue"
      },
      "outputs": [],
      "source": [
        "#Here, we show how to visualize intensity values in 3D.\n",
        "\n",
        "# Pick a frame and a color channel\n",
        "selected_frame = 0\n",
        "selected_color_channel = 1\n",
        "\n",
        "# Select the image\n",
        "selected_image = img[selected_frame,:,:,selected_color_channel]\n",
        "\n",
        "# Create a grid of X and Y values the same size as the image\n",
        "space= np.arange(0, selected_image.shape[0], 1)\n",
        "xx, yy = np.meshgrid(space,space)\n",
        "\n",
        "# Create a figure\n",
        "fig = plt.figure(figsize=(15,7))\n",
        "\n",
        "# Set up the axes for the first plot and show the image\n",
        "ax = fig.add_subplot(1, 2, 1)\n",
        "ax.imshow(selected_image, cmap='gray') # Reds_r\n",
        "\n",
        "# Set up the axes for the second plot and show the 3D surface\n",
        "ax2 = fig.add_subplot(1, 2, 2, projection='3d')\n",
        "ax2.plot_surface(xx, yy , selected_image,  rstride=20, cstride=20, shade=False, cmap='gray')\n",
        "ax2.view_init(20, 45)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aUpbCqxnNgXe"
      },
      "source": [
        "## Intensity in images (Bit depth)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WSMC6pqVhmJd"
      },
      "source": [
        "Bit depth is the information stored on each pixel in the image.\n",
        "\n",
        "Bits  | Range of values: $2^n$\n",
        "---------|------------------\n",
        "1 bit    | 2\n",
        "8 bit    | 256\n",
        "12 bit   | 4096\n",
        "16 bit   | 65536\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "RR-DZx0HpJDv"
      },
      "outputs": [],
      "source": [
        "# Here is a function to convert bit depth in an image\n",
        "# Note - you will lose information if you convert to a lower bit depth\n",
        "def convert_bit_depth(img, target_type_min, target_type_max, target_type):\n",
        "    '''\n",
        "    This function is intended to normalize img and change the image to the specified target_type\n",
        "      img: numpy array\n",
        "      target_type_min: int\n",
        "      target_type_max: int\n",
        "      target_type: str, options are: np.uint\n",
        "    '''\n",
        "    imin = img.min()\n",
        "    imax = img.max()\n",
        "    a = (target_type_max - target_type_min) / (imax - imin)\n",
        "    b = target_type_max - a * imax\n",
        "    new_img = (a * img + b).astype(target_type)\n",
        "    return new_img"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xdsuzjVaN2Px"
      },
      "source": [
        "Check this [link](https://numpy.org/doc/stable/user/basics.types.html) for a complete list of numpy data types."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SHAqS4BZhlOx"
      },
      "outputs": [],
      "source": [
        "# Normalizing and converting images between different bit-depths\n",
        "\n",
        "#Convert an image to unsigned byte format, with values in [0, 1]\n",
        "img_int1 = convert_bit_depth(img, 0,1,target_type=np.bool_)\n",
        "\n",
        "#Convert an image to unsigned byte format, with values in [0, 8]\n",
        "img_int3 = convert_bit_depth(img, 0,8,target_type=np.uint8)\n",
        "\n",
        "#Convert an image to unsigned byte format, with values in [0, 255]\n",
        "img_int8 = convert_bit_depth(img, 0,255,target_type=np.uint8)\n",
        "\n",
        "#Convert an image to unsigned byte format, with values in [0, 255]\n",
        "img_int12 = convert_bit_depth(img, 0,4095,target_type=np.uint16)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zY6HplbNnxTw"
      },
      "outputs": [],
      "source": [
        "# Let's look at the range of values in each image\n",
        "print('Range in 1-bit image: [', np.amin(img_int1),',' ,np.amax(img_int1) , ']' )\n",
        "print('Range in 3-bit image: [', np.amin(img_int3),',' ,np.amax(img_int3) , ']' )\n",
        "print('Range in 8-bit image: [', np.amin(img_int8),',' ,np.amax(img_int8) , ']' )\n",
        "print('Range in 12-bit image: [', np.amin(img_int12),',' ,np.amax(img_int12) , ']' )\n",
        "print('Range in 16-bit image: [', np.amin(img),',' ,np.amax(img) , ']' ) # notice that the max value in this particular image is ~8K, but the maximum possible range in an int16 image is 65.5K"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "W_aPRFIXkqXu"
      },
      "outputs": [],
      "source": [
        "#Let's display the images for a Side-by-side comparison of the same image at different bit depths\n",
        "\n",
        "# Create a figure\n",
        "fig, ax = plt.subplots(1,5, figsize=(40, 20))\n",
        "\n",
        "# Show the image at 1 bit\n",
        "ax[0].imshow(img_int1[0,:,:,0],cmap='gray')\n",
        "ax[0].set_title('1bit',fontsize=24)\n",
        "\n",
        "# Show the image at 3 bit\n",
        "ax[1].imshow(img_int3[0,:,:,0],cmap='gray')\n",
        "ax[1].set_title('3bit', fontsize=24)\n",
        "\n",
        "# Show the image at 8 bit\n",
        "ax[2].imshow(img_int8[0,:,:,0],cmap='gray')\n",
        "ax[2].set_title('8bit', fontsize=24)\n",
        "\n",
        "# Show the image at 12 bit\n",
        "ax[3].imshow(img_int12[0,:,:,0],cmap='gray')\n",
        "ax[3].set_title('12bit', fontsize=24)\n",
        "\n",
        "# Show the image at 16 bit\n",
        "ax[4].imshow(img[0,:,:,0],cmap='gray')\n",
        "ax[4].set_title('16bit', fontsize=24)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qagrynMjV-gD"
      },
      "source": [
        "The human eye cannot distinguish between a 12-bit image and a 16-bit image."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "xfxVTg3qEKRU"
      },
      "outputs": [],
      "source": [
        "# Visualizing the intensity values in images with different bit depths\n",
        "\n",
        "# Let's zoom in to a 10x10 section of the image\n",
        "min_selection_area = 200\n",
        "max_selection_area = min_selection_area+10\n",
        "\n",
        "# Extracting the 10x10 section of the image as a data frame\n",
        "df_3bit = pd.DataFrame(img_int3[0,min_selection_area:max_selection_area,min_selection_area:max_selection_area,0] ) # Range in 3-bit image: [ 0 , 8 ]\n",
        "df_8bit = pd.DataFrame(img_int8[0,min_selection_area:max_selection_area,min_selection_area:max_selection_area,0] ) # Range in 8-bit image: [ 0 , 255 ]\n",
        "df_16bit = pd.DataFrame(img[0,min_selection_area:max_selection_area,min_selection_area:max_selection_area,0] ) # Range in 16-bit image: [ 0 , 65536 ]. In this particular image, the original maximum value is 6380\n",
        "\n",
        "# Plotting the intensity values in the 10x10 section of the image\n",
        "fig, ax = plt.subplots(1,3, figsize=(30, 7))\n",
        "\n",
        "# Plotting the heatmap of a section in the image\n",
        "sn.heatmap(df_3bit, annot=True,cmap=\"gray\",fmt='d', ax=ax[0])\n",
        "ax[0].set_title('3-bit image', fontsize=24)\n",
        "sn.heatmap(df_8bit, annot=True,cmap=\"gray\",fmt='d', ax=ax[1])\n",
        "ax[1].set_title('8-bit image', fontsize=24)\n",
        "sn.heatmap(df_16bit, annot=True,cmap=\"gray\",fmt='d', ax=ax[2])\n",
        "ax[2].set_title('16-bit image', fontsize=24)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ar0uiRhfOlGT"
      },
      "source": [
        "#### File size for two different data types and bit depths\n",
        "\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# How does the bit depth affect the file size? \n",
        "print(random.choice(students))\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KwSKSyPYINNa"
      },
      "outputs": [],
      "source": [
        "# Saving the images to disk\n",
        "tifffile.imwrite('temp_img_int8.tif', img_int8)\n",
        "tifffile.imwrite('temp_img_int16.tif', img)\n",
        "\n",
        "# Loading the images\n",
        "print(\"File size of the 8-bit image in Mb is: \", round(Path('temp_img_int8.tif').stat().st_size/1e6))\n",
        "\n",
        "# How big do you think the file size will be for the 16-bit image? \n",
        "# Uncomment the next line to see the answer\n",
        "print(\"File size of the 16-bit image in Mb is: \", round(Path('temp_img_int16.tif').stat().st_size/1e6))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZjonhgDtzNz2"
      },
      "source": [
        "# What is Color in an image?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MYl5yDAdzUYd"
      },
      "source": [
        "An example of a monochromatic system is the black and white television.\n",
        "\n",
        "<img src= https://media.wired.com/photos/5eacb6979fa72e5901a7bfa3/191:100/w_5100,h_2670,c_limit/Gear-TV-binge-10079749.jpg alt=\"drawing\" width=\"800\"/>\n",
        "\n",
        "Image source: https://www.wired.com/story/what-we-keep-rebinging/"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Here is how to download an image from the web and display it.  Just so you know we really can do it, \n",
        "# let's also flip the image upside down and change the color map to 'hot'.\n",
        "\n",
        "# Load an jpg image from an html link and create a numpy array\n",
        "url = 'https://media.wired.com/photos/5eacb6979fa72e5901a7bfa3/191:100/w_5100,h_2670,c_limit/Gear-TV-binge-10079749.jpg'\n",
        "\n",
        "# Download the image\n",
        "urllib.request.urlretrieve(url, './image.jpg')\n",
        "\n",
        "# Load the image\n",
        "img_jpg = imread('./image.jpg',as_gray=True) \n",
        "print(f'The size of this image is {img_jpg.shape}. The image is {img_jpg.dtype}.')\n",
        "\n",
        "plt.imshow(img_jpg, cmap='magma') # 'Blues'\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jG9zQZVczkUA"
      },
      "source": [
        "Monochromatic images are normally associated with Black and White colors, but we can associate a different color map. Notice that this will only map the color intensity to a predefined colormap."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "UqXI2160zBjU"
      },
      "outputs": [],
      "source": [
        "# There are lots of Colormaps\n",
        "\n",
        "cmaps = [('Perceptually Uniform Sequential', [\n",
        "            'viridis', 'plasma', 'inferno', 'magma', 'cividis']),\n",
        "         ('Sequential', [\n",
        "            'Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds',\n",
        "            'YlOrBr', 'YlOrRd', 'OrRd', 'PuRd', 'RdPu', 'BuPu',\n",
        "            'GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn']),\n",
        "         ('Sequential (2)', [\n",
        "            'binary', 'gist_yarg', 'gist_gray', 'gray', 'bone', 'pink',\n",
        "            'spring', 'summer', 'autumn', 'winter', 'cool', 'Wistia',\n",
        "            'hot', 'afmhot', 'gist_heat', 'copper']),\n",
        "         ('Diverging', [\n",
        "            'PiYG', 'PRGn', 'BrBG', 'PuOr', 'RdGy', 'RdBu',\n",
        "            'RdYlBu', 'RdYlGn', 'Spectral', 'coolwarm', 'bwr', 'seismic']),\n",
        "         ('Cyclic', ['twilight', 'twilight_shifted', 'hsv']),\n",
        "         ('Qualitative', [\n",
        "            'Pastel1', 'Pastel2', 'Paired', 'Accent',\n",
        "            'Dark2', 'Set1', 'Set2', 'Set3',\n",
        "            'tab10', 'tab20', 'tab20b', 'tab20c'])]\n",
        "\n",
        "\n",
        "gradient = np.linspace(0, 1, 256)\n",
        "gradient = np.vstack((gradient, gradient))\n",
        "\n",
        "def plot_color_gradients(cmap_category, cmap_list):\n",
        "    # Create figure and adjust figure height to number of colormaps\n",
        "    nrows = len(cmap_list)\n",
        "    figh = 0.35 + 0.15 + (nrows + (nrows-1)*0.1)*0.22\n",
        "    fig, axs = plt.subplots(nrows=nrows, figsize=(6.4, figh/2))\n",
        "    fig.subplots_adjust(top=1-.35/figh, bottom=.15/figh, left=0.2, right=0.99)\n",
        "\n",
        "    axs[0].set_title(cmap_category + ' colormaps', fontsize=14)\n",
        "\n",
        "    for ax, cmap_name in zip(axs, cmap_list):\n",
        "        ax.imshow(gradient, aspect='auto', cmap=cmap_name)\n",
        "        ax.text(-.01, .5, cmap_name, va='center', ha='right', fontsize=10,\n",
        "                transform=ax.transAxes)\n",
        "\n",
        "    # Turn off *all* ticks & spines, not just the ones with colormaps.\n",
        "    for ax in axs:\n",
        "        ax.set_axis_off()\n",
        "\n",
        "\n",
        "for cmap_category, cmap_list in cmaps:\n",
        "    plot_color_gradients(cmap_category, cmap_list)\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KoN6L5Ga0eQC"
      },
      "source": [
        "Using a colormap in a monochromatic image."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "kn6nZQZj3V-v"
      },
      "outputs": [],
      "source": [
        "# Pick a color map from the list above, and visualize \n",
        "# the image in 2D and 3D using that color map\n",
        "\n",
        "# Choose a color map\n",
        "selected_color_map ='Spectral' #'Spectral' # 'gray'\n",
        "\n",
        "# Make a grid of X and Y values the same size as the image\n",
        "space= np.arange(0, selected_image.shape[0], 1)\n",
        "xx, yy = np.meshgrid(space,space)\n",
        "fig = plt.figure(figsize=(15,7))\n",
        "\n",
        "# Set up the axes for the first plot\n",
        "ax = fig.add_subplot(1, 2, 1)\n",
        "ax.imshow(selected_image,cmap=selected_color_map) # coolwarm\n",
        "\n",
        "# Set up the axes for the second plot\n",
        "ax2 = fig.add_subplot(1, 2, 2, projection='3d')\n",
        "ax2.plot_surface(xx, yy , selected_image,  rstride=20, cstride=20, shade=False, cmap=selected_color_map) #\n",
        "ax2.view_init(20, 45)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o0mGFJ8dM8zl"
      },
      "source": [
        "## Color images.\n",
        "\n",
        "In a color image each pixel in the image is assigned a combination of intensity values for the red, green, and blue color channels, which together determine the color appearance of each pixel.\n",
        "\n",
        "Most image processing libraries use the \"consensus\" order RGB. A very important library for image processing ([OpenCV](https://opencv.org)) uses the order BGR!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "2IlD3VM7sGAL"
      },
      "outputs": [],
      "source": [
        "# Plotting a color image\n",
        "\n",
        "# For many images, the color channels are stored in the last dimension. \n",
        "# For microscopy images, the colors are usually related to different fluorophores.\n",
        "# For example, in this three color image, the first channel is red and corresponds\n",
        "# to phosphorylated RNAP, the second is green and corresponds to phosphorylated RNAP,\n",
        "# the third is blue, which corresponds to RNA.\n",
        "\n",
        "# Make a 3x3 grid of images\n",
        "fig, ax = plt.subplots(1,4, figsize=(20, 7))\n",
        "Red = img[0,:,:,0]\n",
        "\n",
        "# Plotting the color channels\n",
        "ax[0].imshow(Red,cmap='Reds_r')\n",
        "ax[0].set(title='Red Channel (Ch 0) - RNAP')\n",
        "Green = img[0,:,:,1]\n",
        "ax[1].imshow(Green,cmap='Greens_r')\n",
        "ax[1].set(title='Green Channel (Ch 1) - RNAP-Phos')\n",
        "Blue = img[0,:,:,2]\n",
        "ax[2].imshow(Blue,cmap='Blues_r')\n",
        "ax[2].set(title='Blue Channel (Ch 2) - RNA')\n",
        "\n",
        "# Now to create a merged color image. \n",
        "# We need to rescale the image to avoid saturation of the colors.\n",
        "number_channels =3\n",
        "rescaled_image = np.zeros_like(img_int8)\n",
        "# Rescaling image to avoid saturation\n",
        "color_image = img[0,:,:,:].copy()\n",
        "rescaled_image = np.zeros_like(color_image)\n",
        "for index_channels in range (number_channels):\n",
        "  normalized_image_temp = color_image[ :, :, index_channels]\n",
        "  # Pick a percentile value to avoid saturation\n",
        "  max_val = np.percentile(color_image[:,:,index_channels], 97)\n",
        "  # Set values above the percentile to the percentile value\n",
        "  normalized_image_temp [normalized_image_temp > max_val] = max_val\n",
        "  # Saved the thresholded image in the rescaled image\n",
        "  rescaled_image[:,:,index_channels] = normalized_image_temp\n",
        "\n",
        "# to show a color image in matplotlib, we need to convert it to uint8.\n",
        "img_int8_color = convert_bit_depth(rescaled_image, 0,255,target_type=np.uint8)\n",
        "\n",
        "# Now, let's plot the color image\n",
        "ax[3].imshow(img_int8_color)\n",
        "ax[3].set(title='Color image')\n",
        "ax[0].axis('off')\n",
        "ax[1].axis('off')\n",
        "ax[2].axis('off')\n",
        "ax[3].axis('off')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uf1706KGOOv3"
      },
      "source": [
        "## Working with images in Python"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MjUK0rpoOW6K"
      },
      "source": [
        "### Basic image manipulation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "it4IFSnMOqFp"
      },
      "source": [
        "#### Slicing"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-x_xa1Ow4kW1"
      },
      "source": [
        "In this section, we select parts of the image.\n",
        "\n",
        "The image is a numpy array with dimensions:\n",
        "```\n",
        "image[time, y-axis, x-axis, colors]\n",
        "```\n",
        "\n",
        "If we need to select the following elements:\n",
        "* timepoint (frame) 5\n",
        "* y-axis from 100 to 200 pixel\n",
        "* x-axis from 230 to 300 pixel\n",
        "* \"Green\" color (Color 1 in the standard format [R,G,B]),\n",
        "\n",
        "We would slice the numpy array as follows:\n",
        "\n",
        "```\n",
        "image[5, 99:200, 229:300, 1]\n",
        "```\n",
        "\n",
        "* Notice that when we slice [start: end(not including this value) :step]\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "wcJ0-aZlOj5S"
      },
      "outputs": [],
      "source": [
        "# How would you crop the image to only show the middle 3rd of the image in both directions?\n",
        "random.choice(students)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-s8VySLIa5cO"
      },
      "outputs": [],
      "source": [
        "# Recall that the image shape is:\n",
        "print('Image shape: ', img_int8.shape)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FkxEwsc54BYp"
      },
      "outputs": [],
      "source": [
        "# Plotting a subsection of the image\n",
        "# Time point: 0\n",
        "# Y-range: [99:200]\n",
        "# X-range: [229:300]\n",
        "# Channel: Green (1)\n",
        "plt.imshow(img_int8[0, 99:200, 229:300, 1], cmap='Greens')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "plt.imshow(img_int8[0, 99:200, 229:300, 1], cmap='Greens_r')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RRwow-M54eHs"
      },
      "outputs": [],
      "source": [
        "# Ploting a subsection of the image\n",
        "# Time point: 22\n",
        "# Y-range: [230:300]\n",
        "# X-range: [155:350]\n",
        "# Channel: Blue (2)\n",
        "plt.figure(figsize=(7,7))\n",
        "plt.imshow(img[22,230:300,155:350,2],cmap='gray') \n",
        "# Notice that only one time point and one color is plotted\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wqf1CPcqOv_x"
      },
      "source": [
        "#### Thresholding"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# dont do this!\n",
        "img_copy = img.copy()\n",
        "img_copy[:]= img_copy[:]*0\n",
        "#print(img)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "BW2FvcR-Psk1"
      },
      "outputs": [],
      "source": [
        "# How do you make a copy of an image?\n",
        "img_copy = img.copy() # Making a copy of our original image\n",
        "img_copy = 1.0*img # Making a copy of our original image\n",
        "# img_copy = img # DO NOT DO THIS. This is not making a copy of the image, but rather creating a new reference to the same image.\n",
        "plt.imshow(img_copy[0,:,:,2],cmap='gray') "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3i09YvCR692X"
      },
      "source": [
        "## Applying thresholds to images.\n",
        "How would you set values less than the some threshold equal to zero?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "random.choice(students)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "img_copy = img.copy() # Making a copy of our img\n",
        "img_section = img_copy[0,:,:,2] # Selecting a time point and color channel # test channel 2\n",
        "threshold_value = 1000\n",
        "img_section=img_section*(img_section >= threshold_value )  # Thresholding image values larger than 'threshold_value' and replacing with 'threshold_value'\n",
        "plt.imshow(img_section, cmap='gray') \n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# How would you set values more than the the image mean equal to zero?\n",
        "# img_section[??? ]=???  # Thresholding image values larger than the mean equal to the mean\n",
        "plt.imshow(img_section,cmap='gray') "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Au0WrfoG0CoI"
      },
      "outputs": [],
      "source": [
        "# Making a binary image.\n",
        "\n",
        "img_copy = img.copy() # Making a copy of our img\n",
        "img_section = img_copy[34,:,:,2] # Selecting a time point and color channel # test channel 2\n",
        "\n",
        "threshold_value = 1000\n",
        "img_binary = img_section>threshold_value\n",
        "print(np.min(img_binary) , np.max(img_binary))\n",
        "plt.imshow(img_binary,cmap='gray')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "Ppdq8f4DvBSo"
      },
      "outputs": [],
      "source": [
        "# Plotting the image as the 3d dimension figure. Thresholding.\n",
        "\n",
        "# Plotting the image as the 3d dimension figure.\n",
        "space= np.arange(0, img_section.shape[0], 1)\n",
        "xx, yy = np.meshgrid(space,space)\n",
        "color_map = 'plasma' #'Greys_r'\n",
        "\n",
        "#create a figure with 2x2 subplots\n",
        "fig, ax = plt.subplots(2,2, figsize=(15, 15))\n",
        "# show the image in 2D\n",
        "ax[0][0].imshow(img_section,cmap=color_map) # coolwarm\n",
        "ax[0][0].set_title('Original Image')\n",
        "\n",
        "# show the image in 3D\n",
        "ax[0][1] = fig.add_subplot(2, 2, 2, projection='3d')\n",
        "ax[0][1].plot_surface(xx, yy , img_section,  rstride=20, cstride=20, shade=False, cmap=color_map) #\n",
        "ax[0][1].view_init(20, 45)\n",
        "\n",
        "# show the binary image in 2D\n",
        "ax[1][0].imshow(img_binary,cmap=color_map) # coolwarm\n",
        "ax[1][0].set_title('Binary Image')\n",
        "\n",
        "# show the binary image in 3D\n",
        "ax[1][1] = fig.add_subplot(2, 2, 4, projection='3d')\n",
        "ax[1][1].plot_surface(xx, yy , img_binary,  rstride=20, cstride=20, shade=False, cmap=color_map) #\n",
        "ax[1][1].view_init(20, 45)\n",
        "\n",
        "plt.show()\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "np.max(img_binary)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xmfGuqIbOlwY"
      },
      "source": [
        "## Filters"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JuLkk8wFLJvD"
      },
      "source": [
        "[Filters](https://ai.stanford.edu/~syyeung/cvweb/tutorial1.html) are used for:\n",
        "\n",
        "*   Noise reduction\n",
        "*   Edge detection\n",
        "*   Sharpening\n",
        "*   Blurring\n",
        "\n",
        "The mathematical operation is a 2D convolution. This convolution involves defining a smaller kernel matrix and applying the same mathematical operation to each pixel in the entire image. A more complete explanation can be found in this amizing [video](https://youtu.be/8rrHTtUzyZA?t=72) by 3Blue1Brown.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2gbyFSTKkig2"
      },
      "source": [
        "![alt text](FigsA/Module_1_1/Slide5.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xoq0pDWRgago"
      },
      "source": [
        "##### Gaussian Filter. Noise reduction and blurring."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mqPRfQwVgfNX"
      },
      "source": [
        "$G_\\sigma = \\frac{1}{2\\pi\\sigma^2}e{\\frac{x^2+y^2}{2\\sigma^2}}$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "ol0Hsh3RhCRn"
      },
      "outputs": [],
      "source": [
        "#@title Section that creates the Gaussian Kernel Matrix\n",
        "def gaussian_kernel (size_matrix,sigma):\n",
        "  '''\n",
        "  This function returns a normalized gaussian kernel matrix\n",
        "  size_matrix : int\n",
        "  sigma: float\n",
        "  '''\n",
        "  ax = np.linspace(-(size_matrix - 1) / 2., (size_matrix - 1) / 2., size_matrix)\n",
        "  xx, yy = np.meshgrid(ax, ax)\n",
        "  kernel = np.exp(-0.5 * (np.square(xx) + np.square(yy)) / np.square(sigma))\n",
        "  kernel = kernel/kernel.sum() # Normalizing to the sum\n",
        "  return kernel\n",
        "\n",
        "# Gaussian kernel matrix for different sigmas\n",
        "kernel_gaussian_sigma_3 = gaussian_kernel (size_matrix=20,sigma=3)\n",
        "kernel_gaussian_sigma_5 = gaussian_kernel (size_matrix=20,sigma=5)\n",
        "kernel_gaussian_sigma_10 = gaussian_kernel (size_matrix=20,sigma=10)\n",
        "\n",
        "\n",
        "print(sum(kernel_gaussian_sigma_3.flatten()))\n",
        "print(sum(kernel_gaussian_sigma_5.flatten()))\n",
        "print(sum(kernel_gaussian_sigma_10.flatten()))\n",
        "\n",
        "\n",
        "# Side-by-side comparison\n",
        "fig, ax = plt.subplots(1,3, figsize=(20, 10))\n",
        "ax[0].imshow(kernel_gaussian_sigma_3,cmap='gray')\n",
        "ax[0].set(title='Gaussian kernel $\\sigma$ =3')\n",
        "ax[1].imshow(kernel_gaussian_sigma_5,cmap='gray')\n",
        "ax[1].set(title='Gaussian kernel $\\sigma$ =5')\n",
        "ax[2].imshow(kernel_gaussian_sigma_10,cmap='gray')\n",
        "ax[2].set(title='Gaussian kernel $\\sigma$ =10')\n",
        "\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "qqdWFAiwsBad"
      },
      "outputs": [],
      "source": [
        "# Visualizing a Gaussian Kernel in 3D\n",
        "size_spot = 20\n",
        "spot_sigma = 3\n",
        "space = np.linspace(-(size_spot - 1) / 2., (size_spot - 1) / 2., size_spot)\n",
        "xx, yy = np.meshgrid(space, space)\n",
        "kernel = np.exp(-0.5 * (np.square(xx) + np.square(yy)) / np.square(spot_sigma))\n",
        "#kernel = kernel / np.amax(kernel) * 255  # Normalizing with respect to max and changing the range to [0,255]\n",
        "kernel = kernel/kernel.sum()\n",
        "print('Sum of the intensity in the kernel is: ',sum(kernel.flatten()))\n",
        "# Plotting\n",
        "fig = plt.figure(figsize=plt.figaspect(0.3))\n",
        "# Set up the axes for the first plot\n",
        "ax = fig.add_subplot(1, 2, 1)\n",
        "ax.imshow(kernel,cmap='gray') # Reds_r\n",
        "# Set up the axes for the second plot\n",
        "ax = fig.add_subplot(1, 2, 2, projection='3d')\n",
        "ax.plot_surface(xx, yy, kernel, cmap='gray')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MSNOiAivMRgj"
      },
      "source": [
        "Example using [gaussian filter scipy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html). For a complete  list of filters in scipy, use the following [link](https://docs.scipy.org/doc/scipy/reference/ndimage.html)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "wCvsao4vLWN9"
      },
      "outputs": [],
      "source": [
        "# Importing the library with the filter modules\n",
        "from scipy.ndimage import gaussian_filter\n",
        "\n",
        "img_copy = img.copy() # Making a copy of our img\n",
        "img_section = img_copy[0,:,:,0] # Selecting a time point and color channel\n",
        "\n",
        "# Applying the filter\n",
        "img_gaussian_filter_simga_1 = gaussian_filter(img_section, sigma=1)\n",
        "img_gaussian_filter_simga_10 = gaussian_filter(img_section, sigma=10)\n",
        "\n",
        "# Side-by-side comparison\n",
        "fig, ax = plt.subplots(1,3, figsize=(30, 10))\n",
        "ax[0].imshow(img_section,cmap='gray')\n",
        "ax[0].set_title('Original',fontsize=24)\n",
        "\n",
        "# Noise reduction\n",
        "ax[1].imshow(img_gaussian_filter_simga_1,cmap='gray')\n",
        "ax[1].set_title('Gaussian Filter $\\sigma$ =1 Noise reduction',fontsize=24)\n",
        "\n",
        "# Blurring\n",
        "ax[2].imshow(img_gaussian_filter_simga_10,cmap='gray')\n",
        "ax[2].set_title('Gaussian Filter $\\sigma$ =10 Image Blurring',fontsize=24)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qt7m9Z77NDK9"
      },
      "source": [
        "Filters in scikit-image ([difference of gaussians](https://scikit-image.org/docs/stable/api/skimage.filters.html#skimage.filters.difference_of_gaussians)).\n",
        "\n",
        "This filter is used to locate elements between a low and a high value.\n",
        "\n",
        " For a complete list of filters in scikit-image, use the following [link](https://scikit-image.org/docs/stable/api/skimage.filters.html)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RmrdlxfqN4vf"
      },
      "outputs": [],
      "source": [
        "# Importing skimage filters module\n",
        "from skimage.filters import difference_of_gaussians\n",
        "\n",
        "img_copy = img.copy() # Making a copy of our img\n",
        "img_section = img_copy[0,:,:,0] # Selecting a time point and color channel\n",
        "\n",
        "# Applying the filter to our image\n",
        "img_diff_gaussians = difference_of_gaussians(img_section,low_sigma=1, high_sigma=10)\n",
        "#img_diff_gaussians = difference_of_gaussians(img_section,low_sigma=5, high_sigma=10)\n",
        "\n",
        "# Side-by-side comparison\n",
        "fig, ax = plt.subplots(1,2, figsize=(20, 10))\n",
        "ax[0].imshow(img_section,cmap='gray')\n",
        "ax[0].set_title('Original',fontsize=24)\n",
        "ax[1].imshow(img_diff_gaussians,cmap='gray')\n",
        "ax[1].set_title('Difference of Gaussians', fontsize=24)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xu1xikL-O3iP"
      },
      "source": [
        "#### Rotation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "aV3x-0XvWdgR"
      },
      "source": [
        "\n",
        "Library ([rotate scipy](https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.rotate.html#scipy.ndimage.rotate))."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bYkFWx7_LUc2"
      },
      "outputs": [],
      "source": [
        "# Importing skimage rotation module\n",
        "from scipy import ndimage as nd\n",
        "\n",
        "img_copy = img.copy() # Making a copy of our img\n",
        "img_section = img_copy[0,:,:,0] # Selecting a time point and color channel\n",
        "\n",
        "# Rotating the image to a given right angle. \n",
        "selected_angle = 90\n",
        "img_rotation = nd.rotate(img_section, angle=selected_angle)\n",
        "\n",
        "# Side-by-side comparison\n",
        "fig, ax = plt.subplots(1,2, figsize=(20, 10))\n",
        "ax[0].imshow(img_section,cmap='gray')\n",
        "ax[0].set_title('Original',fontsize=24)\n",
        "ax[1].imshow(img_rotation,cmap='gray')\n",
        "ax[1].set_title( 'Image rotated by '+str(selected_angle)+ ' degrees' ,fontsize=24)\n",
        "plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dVqFESd_PFNk"
      },
      "source": [
        "#### Image transformation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2I99VaYAJ5pE"
      },
      "source": [
        "Consists of applying rotation, scaling, and translation processes to the image."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i-dok6uhYYp-"
      },
      "source": [
        "List of available [transformations in skimage](https://scikit-image.org/docs/stable/auto_examples/transform/plot_transform_types.html). Blog with more information about [applying transformations to images](https://towardsdatascience.com/image-processing-with-python-applying-homography-for-image-warping-84cd87d2108f)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "df60Z4LwYeen"
      },
      "outputs": [],
      "source": [
        "# Importing skimage transformation module\n",
        "from skimage import transform\n",
        "\n",
        "img_copy = img.copy() # Making a copy of our img\n",
        "img_section = img_copy[0,:,:,0] # Selecting a time point and color channel\n",
        "\n",
        "#  Transformation matrix\n",
        "tform = transform.SimilarityTransform(\n",
        "    scale = 0.95,                  # float, scaling value\n",
        "    rotation = np.pi/90,           # Rotation angle in counter-clockwise direction as radians. pi/180 rad = 1 degrees\n",
        "    translation=(100, 1))          # (x, y) values for translation\n",
        "print('Transformation matrix : \\n', tform.params , '\\n')\n",
        "\n",
        "# Applying the transformation\n",
        "tf_img = transform.warp(img_section, tform.inverse)\n",
        "\n",
        "# Side-by-side comparison\n",
        "fig, ax = plt.subplots(1,2, figsize=(20, 10))\n",
        "ax[0].imshow(img_section,cmap='gray')\n",
        "ax[0].set_title('Original',fontsize=24)\n",
        "ax[1].imshow(tf_img,cmap='gray')\n",
        "ax[1].set_title('Transformation',fontsize=24)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9ERlWkn7xnu3"
      },
      "source": [
        "This process is used before merging images taken from multiple cameras. Image registration [Check this video](https://www.youtube.com/watch?v=nNZJw0kgzdg&list=LL&index=7)\n",
        "\n",
        "![alt text](FigsA/Module_1_1/Slide7.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ff6-FeqUA9he"
      },
      "source": [
        "## Working with a sequence of images"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jZfytq3ss_vI"
      },
      "source": [
        "### Video"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "\n",
        "# !conda install -c conda-forge ffmpeg"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "R3dHENYartCw"
      },
      "outputs": [],
      "source": [
        "#  %conda install -c conda-forge ffmpeg\n",
        "\n",
        "# Visualizing all frames (time points).\n",
        "# blit=True re-draws only the parts that have changed\n",
        "fig,axes = plt.subplots(1,3,dpi=120,figsize=(8,3))\n",
        "i=0\n",
        "# Define inital frames\n",
        "Red = img[i,:,:,0]\n",
        "im1 = axes[0].imshow(Red,cmap='Reds_r')\n",
        "Green = img[i,:,:,1]\n",
        "im2 = axes[1].imshow(Green,cmap='Greens_r')\n",
        "Blue = img[i,:,:,2]\n",
        "im3 =  axes[2].imshow(Blue,cmap='Blues_r')\n",
        "axes[0].axis('off')\n",
        "axes[1].axis('off')\n",
        "axes[2].axis('off')\n",
        "\n",
        "def movieFrame(i):\n",
        "  Red = img[i,:,:,0]\n",
        "  Green = img[i,:,:,1]\n",
        "  Blue = img[i,:,:,2]\n",
        "  images = [Red,Green,Blue]\n",
        "  image_handles = [im1,im2,im3]\n",
        "  for k,image_n in enumerate(images):\n",
        "    image_handles[k].set_array(images[k])\n",
        "  return image_handles\n",
        "\n",
        "plt.close()\n",
        "anim = animation.FuncAnimation(fig, movieFrame, frames=img.shape[0], interval=20, blit=True)\n",
        "from IPython.display import HTML\n",
        "HTML(anim.to_html5_video())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UJJ1bUOpqH-M"
      },
      "source": [
        "### Images with 3-dimensional space, fluorescence in situ hybridization (FISH) images."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "rZf1xgAxqORF"
      },
      "outputs": [],
      "source": [
        "# # Images with a 3D space\n",
        "\n",
        "# Downloading the image to our local drive\n",
        "# You can skip this, since I provide the files over GitHub\n",
        "# found_files = list(Path('.').rglob('FISH_example.zip'))\n",
        "# if len(found_files) != 0:\n",
        "#   print(f\"File already downloaded and can be found in {found_files[0]}.\")\n",
        "# else:\n",
        "#   !wget --no-check-certificate 'https://www.dropbox.com/s/6y7yqlmlnnxm7rs/FISH_example.zip?dl=0' -r -A 'uc*' -e robots=off -nd -O 'FISH_example.zip'\n",
        "#   # !wget --no-check-certificate 'https://www.dropbox.com/s/1ZR6nY9kscgLEefZFwBPwsbog-JlVwE2B/FISH_example.zip?dl=0' -r -A 'uc*' -e robots=off -nd -O 'FISH_example.zip'\n",
        "#   !unzip Module4-ImageProcessing/FISH_example.zip\n",
        "\n",
        "# Importing the image as variable img.\n",
        "# You might need to change the path to the file depending on your local setup.\n",
        "figName_FISH = 'FISH_example.tif'\n",
        "img_FISH = imread(figName_FISH)\n",
        "\n",
        "# This image has dimension [Z,Y,X]\n",
        "print(f\"The shape of the image is {img_FISH.shape}\")\n",
        "\n",
        "#Removing outliers\n",
        "max_val = np.percentile(img_FISH, 99)\n",
        "img_FISH [img_FISH> max_val] = max_val\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "I3xubFa5AwsB"
      },
      "outputs": [],
      "source": [
        "# Plotting the FISH image\n",
        "number_columns =4\n",
        "number_plots = img_FISH.shape[0]\n",
        "number_rows = (number_plots + number_columns - 1) // number_columns  \n",
        "\n",
        "fig, ax = plt.subplots(number_rows, number_columns, figsize=(12, 15))\n",
        "ax = ax.flatten()  # Flatten in case of a 2D array\n",
        "for i in range(number_plots):\n",
        "    ax[i].imshow(img_FISH[i, :, :], cmap='gray')\n",
        "    ax[i].set(title='z = ' + str(i))\n",
        "    ax[i].axis('off')\n",
        "for i in range(number_plots, len(ax)):\n",
        "    ax[i].axis('off')  # Turn off unused subplots\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "PjAK5e47yHaP"
      },
      "outputs": [],
      "source": [
        "# FISH visualizer\n",
        "def FISH_viewer(z_value=1):\n",
        "    '''\n",
        "    This function is intended to display an image from an array of images (specifically, video: img_int8). \n",
        "         img_int8 is a numpy array with dimension [T,Y,X,C].\n",
        "    drop_channel : str with options 'Ch_0', 'Ch_1', 'Ch_2', 'All'\n",
        "    time: int with range 0 to the number of frames in video.\n",
        "    '''\n",
        "    plt.figure(1)\n",
        "    temp_FISH_image = img_FISH[z_value,:,:]\n",
        "    plt.imshow(temp_FISH_image,cmap='gray')\n",
        "    plt.show()\n",
        "\n",
        "# Defining an interactive plot\n",
        "# %pip install ipywidgets\n",
        "from ipywidgets import interactive, HBox, VBox, Layout\n",
        "import ipywidgets as widgets\n",
        "\n",
        "interactive_plot = interactive(FISH_viewer,z_value = widgets.IntSlider(min=0,max=img_FISH.shape[0]-1,step=1,value=0,description='z-value'))       \n",
        "# time slider parameters\n",
        "# Creates the controls\n",
        "controls = HBox(interactive_plot.children[:-1], layout = Layout(flex_flow='row wrap'))\n",
        "# Creates the outputs\n",
        "output = interactive_plot.children[-1]\n",
        "\n",
        "# Display the controls and output as an interactive widget\n",
        "display(VBox([controls, output]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_FLGKKdKmR_T"
      },
      "source": [
        "## Operations on multiple images"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wYkIt7-Okq9U"
      },
      "source": [
        "![alt text](FigsA/Module_1_1/Slide6.png)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kVOEK4mKzCul"
      },
      "source": [
        "Maximum projections"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "aq7o2XFay0bC"
      },
      "outputs": [],
      "source": [
        "# Making a copy of our sequence of images\n",
        "img_FISH_copy = img_FISH.copy() # Making a copy of our img\n",
        "\n",
        "# Applying a maximum projection\n",
        "#img_max_z_projection = np.max(img_FISH, axis=0)\n",
        "img_max_z_projection = np.mean(img_FISH, axis=0)\n",
        "\n",
        "# Plotting\n",
        "plt.figure(figsize=(7,7))\n",
        "plt.imshow(img_max_z_projection,cmap='Greys_r')\n",
        "plt.axis('off')\n",
        "plt.show()\n",
        "\n",
        "# Printing results\n",
        "print('Dimensions on the original sequence of images :', img_FISH.shape, '\\n')\n",
        "print('Dimensions on the maximum projection :', img_max_z_projection.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3LKgKz3UdL1c"
      },
      "source": [
        "Normalizing intensity for every channel and time point"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5tQxaYCN2X0M"
      },
      "source": [
        "Max-Min Normalization\n",
        "\n",
        "$img_{norm} = \\frac{img-min(img)}{max(img)-min(img)}$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "X873jFGndLCK"
      },
      "outputs": [],
      "source": [
        "img_normalized = np.zeros_like(img)   # Preallocating memory\n",
        "number_timepoints, y_dim, x_dim, number_channels = img.shape[0], img.shape[1], img.shape[2], img.shape[3] # Obtaining the dimensions size\n",
        "\n",
        "# Normalization using a nested for-loop\n",
        "for index_channels in range (number_channels): # Iteration for every channel\n",
        "    for index_time in range (number_timepoints): # Iterating for every time point\n",
        "        max_val = np.amax(img[index_time,:,:,index_channels])\n",
        "        min_val = np.amin(img[index_time,:,:,index_channels])\n",
        "        img_normalized[index_time,:,:,index_channels] = (img[index_time,:,:,index_channels]-min_val) / (max_val-min_val) # Normalization\n",
        "\n",
        "# Printing the output\n",
        "print('Range values in the original sequence of images: (' , np.amin(img) ,',', np.amax(img) ,')\\n' )\n",
        "print('Range values in the normalized sequence of images: (' , np.amin(img_normalized) ,',', np.amax(img_normalized) ,')\\n' )\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Oh_Ubv6ARPFM"
      },
      "source": [
        "Transposing dimensions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vpN1ZaXqRRcp"
      },
      "outputs": [],
      "source": [
        "# Making a copy of our sequence of images\n",
        "img_int8_copy = img_int8.copy() # Making a copy of our img [T, Y, X, C]\n",
        "\n",
        "# Reshaping the video. Changing the Time position (0) to the last place (3).  [C, Y, X, T]\n",
        "img_transposed = np.transpose(img_int8_copy, (3, 1, 2, 0))\n",
        "\n",
        "# Printing results\n",
        "print('Dimensions on the original sequence of images :', img_int8_copy.shape, '\\n')\n",
        "print('Dimensions on the maximum projection :', img_transposed.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "UFBJIdmvbTKM"
      },
      "source": [
        "# References"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_SzkyFx2bYPG"
      },
      "source": [
        "* Images downloaded from https://figshare.com from publication: \"Forero-Quintero, L.S., Raymond, W., Handa, T. et al. Live-cell imaging reveals the spatiotemporal organization of endogenous RNA polymerase II phosphorylation at a single gene. Nat Commun 12, 3158 (2021). https://doi.org/10.1038/s41467-021-23417-0\"\n",
        "\n",
        "* Gonzalez, Rafael & Faisal, Zahraa. (2019). Digital Image Processing Second Edition.\n",
        ">\n"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "private_outputs": true,
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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.10.14"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}