slideshow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e8a7ba2b",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Quantifying solidification of metallic alloys with scikit-image\n",
    "\n",
    "<center><h2>C. Gus Becker</h2></center>\n",
    "<center>\n",
    "    <img src=\"images/csm_logo.png\" width=\"400\">\n",
    "    <img src=\"images/canfsa_logo.png\" width=\"300\">\n",
    "</center>\n",
    "<center><h3>&</h3></center>\n",
    "<center><h2>Marianne Corvellec</h2></center>\n",
    "<center>\n",
    "    <img src=\"images/skimage_logo.png\" width=\"300\">\n",
    "    <img src=\"images/igdore_logo.png\" width=\"300\">\n",
    "</center>\n",
    "\n",
    "#### BIDS ImageXD Conference — May 18, 2021 ― Session 4: Materials Science"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7f4ce47",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Context and Motivation\n",
    "<img src=\"images/melt_pool.jpg\">\n",
    "<h4><a href=\"https://commons.wikimedia.org/wiki/File:Selective_laser_melting_system_schematic.jpg\">Materialgeeza</a>, <a href=\"https://creativecommons.org/licenses/by-sa/3.0\">CC BY-SA 3.0</a> (Wikimedia Commons)</h4>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aea7cf24",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction to scikit-image\n",
    "<br />\n",
    "<img src=\"images/skimage_logo.png\">\n",
    "<br />\n",
    "https://scikit-image.org/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4867558",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import skimage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22a41e37",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60682f2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat = data.cat()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "807eab62",
   "metadata": {},
   "outputs": [],
   "source": [
    "type(cat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c861f89b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat.ndim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "404098e0",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import plotly.express as px"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d10368cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "px.imshow(cat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f58f631",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "retina = data.retina()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a80bd008",
   "metadata": {},
   "outputs": [],
   "source": [
    "px.imshow(retina)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b472ced2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import color, filters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf5da6b6",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "image = retina[:, :, 2]\n",
    "\n",
    "thresh = filters.threshold_multiotsu(image)\n",
    "regions = np.digitize(image, bins=thresh)\n",
    "\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(10, 5))\n",
    "ax[0].imshow(image)\n",
    "ax[0].set_title('Channel 2 in RGB')\n",
    "ax[1].imshow(regions)\n",
    "ax[1].set_title('Multi-Otsu thresholding')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b613289f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "conv_img = color.rgb2lab(retina)\n",
    "image = conv_img[:, :, 1]\n",
    "\n",
    "thresh = filters.threshold_multiotsu(image)\n",
    "regions = np.digitize(image, bins=thresh)\n",
    "\n",
    "fig, ax = plt.subplots(ncols=2, figsize=(10, 5))\n",
    "ax[0].imshow(image)\n",
    "ax[0].set_title('Channel 1 in CIELAB')\n",
    "ax[1].imshow(regions)\n",
    "ax[1].set_title('Multi-Otsu thresholding')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c33fdb8d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Track Solid-Liquid Interface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0867925",
   "metadata": {},
   "outputs": [],
   "source": [
    "import imageio\n",
    "\n",
    "collection = []\n",
    "reader = imageio.get_reader('data/nickel_solidification.tif')\n",
    "for img in reader:\n",
    "    collection.append(img)\n",
    "    \n",
    "len(collection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8296ddc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "collection[0].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bad6eeee",
   "metadata": {},
   "outputs": [],
   "source": [
    "full_frames = np.array(collection)\n",
    "full_frames.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f3edc1a",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig = px.imshow(full_frames, animation_frame=0, binary_string=True)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cfdddd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "frames = np.stack([img[60:110, 175:250] for img in collection])\n",
    "frames.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a06c3324",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig = px.imshow(frames, animation_frame=0, binary_string=True)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b8d95b5",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Import local modules containing plotting and processing functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1c968b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import plotting module from local package\n",
    "from solidification import plot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08cef2ed",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Show raw images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f01d4a22",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plot.plot_imgs(\n",
    "    [frames[1], frames[4], frames[7]],\n",
    "    row_title='Raw Images'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b69ab9d3",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Apply Gaussian filter and subtract current image from succeeding image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26314e5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import filters\n",
    "\n",
    "def sub_imgs(collection, img_i_list):\n",
    "    sub_img_list = []\n",
    "    for i in img_i_list:\n",
    "        sub_img_list.append(\n",
    "            filters.gaussian(collection[i + 1]) \n",
    "            - filters.gaussian(collection[i])\n",
    "        )\n",
    "    return sub_img_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df0f7a82",
   "metadata": {},
   "outputs": [],
   "source": [
    "sub_img_list = sub_imgs(frames, [1, 4, 7])\n",
    "\n",
    "fig, axes = plot.plot_imgs(\n",
    "    sub_img_list,\n",
    "    row_title='Subtracted from Succeeding'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c990008f",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Clip highest and lowest intensities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "545ced70",
   "metadata": {},
   "outputs": [],
   "source": [
    "def clip_imgs(img_list, low=5, high=95):\n",
    "    clip_img_list = []\n",
    "    for img in img_list:\n",
    "        # Clip top and bottom image intensities \n",
    "        # (assign low/high for all values below/above low/high)\n",
    "        p_low, p_high = np.percentile(img, [low, high])\n",
    "        img -= p_low\n",
    "        img[(img < 0.0)] = 0.0\n",
    "        img = img / p_high\n",
    "        img[(img > 1.0)] = 1.0\n",
    "        clip_img_list.append(img)\n",
    "    return clip_img_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f703dd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "clip_img_list = clip_imgs(sub_img_list)\n",
    "\n",
    "fig, axes = plot.plot_imgs(\n",
    "    clip_img_list, \n",
    "    row_title='Clip Intensities'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d1d061e",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Invert image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ddd788c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def invert_imgs(img_list):\n",
    "    invert_img_list = []\n",
    "    for img in img_list:\n",
    "        img = 1 - img\n",
    "        invert_img_list.append(img)\n",
    "    return invert_img_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec81093c",
   "metadata": {},
   "outputs": [],
   "source": [
    "invert_img_list = invert_imgs(clip_img_list)\n",
    "\n",
    "fig, axes = plot.plot_imgs(\n",
    "    invert_img_list,\n",
    "    row_title='Inverted'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1cbfcca",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Denoise using `skimage.restoration.denoise_tv_chambolle`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42803490",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import restoration\n",
    "\n",
    "def denoise_imgs(img_list):\n",
    "    denoise_img_list = []\n",
    "    for img in img_list:\n",
    "        img = restoration.denoise_tv_chambolle(img, weight=0.15)\n",
    "        denoise_img_list.append(img)\n",
    "    return denoise_img_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2d577cb",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "denoise_img_list = denoise_imgs(invert_img_list)\n",
    "\n",
    "fig, axes = plot.plot_imgs(\n",
    "    denoise_img_list, \n",
    "    row_title='Denoised'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a21f013",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Threshold and create binary image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9382c360",
   "metadata": {},
   "outputs": [],
   "source": [
    "def binarize_imgs(img_list, thresh_val=0.4):\n",
    "    binarize_img_list = []\n",
    "    for img in img_list:\n",
    "        img_mask = np.zeros(img.shape, dtype=int)\n",
    "        img_mask[(img > thresh_val)] = 1\n",
    "        binarize_img_list.append(img_mask)\n",
    "    return binarize_img_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de565e9d",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "binarize_img_list = binarize_imgs(denoise_img_list)\n",
    "\n",
    "fig, axes = plot.plot_imgs(\n",
    "    binarize_img_list, \n",
    "    row_title='Binarized'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39490b9d",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Filter minimum region size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "def77126",
   "metadata": {},
   "outputs": [],
   "source": [
    "px.imshow(binarize_img_list[0], color_continuous_scale='viridis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d51b557",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from skimage import measure\n",
    "\n",
    "label_img = measure.label(binarize_img_list[0])\n",
    "# We have been processing images 1, 4, and 7, so \n",
    "# binarize_img_list[0] corresponds to collection[1]\n",
    "img_label_overlay = color.label2rgb(\n",
    "    label_img, image=frames[1], bg_label=0\n",
    ")\n",
    "\n",
    "px.imshow(img_label_overlay, color_continuous_scale='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a54cb95",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "props = measure.regionprops_table(\n",
    "    label_img, properties=('label', 'area', 'bbox')\n",
    ")\n",
    "props_df = pd.DataFrame(props)\n",
    "props_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "979ebfb8",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "props_df.sort_values('area', ascending=False).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8e48c53",
   "metadata": {},
   "outputs": [],
   "source": [
    "regions = measure.regionprops(label_img)\n",
    "\n",
    "[region.area for region in regions if region.area > 100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eeafbf00",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from skimage import morphology\n",
    "\n",
    "def filter_size_masks(mask_list, min_size=100):\n",
    "    filter_size_img_list = []\n",
    "    for mask in mask_list:\n",
    "        labels = morphology.label(mask, connectivity=2)\n",
    "        regions = measure.regionprops(labels)\n",
    "        filtered_mask = np.zeros(labels.shape, dtype=int)\n",
    "        for region in regions:\n",
    "            if region.area > min_size:\n",
    "                filtered_mask[(labels == region.label)] = 1\n",
    "                filter_size_img_list.append(filtered_mask)\n",
    "    return filter_size_img_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5211b69",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "filter_size_img_list = filter_size_masks(binarize_img_list)\n",
    "len(filter_size_img_list)\n",
    "fig, axes = plot.plot_imgs(\n",
    "    filter_size_img_list, \n",
    "    row_title='Size Filtered'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f52712a3",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Overlay region bounding box on raw image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7130bb3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "props_df.sort_values('area', ascending=False).head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dfdaaf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "[region.bbox for region in regions if region.area > 100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f140e28e",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "[region.bbox for region in regions if region.area > 100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be0ac056",
   "metadata": {},
   "outputs": [],
   "source": [
    "px.imshow(img_label_overlay, color_continuous_scale='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7516a9e",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def gen_bbox_list(mask_list, min_size=100):\n",
    "    bbox_list = []\n",
    "    for mask in mask_list:\n",
    "        labels = morphology.label(mask, connectivity=2)\n",
    "        regions = measure.regionprops(labels)\n",
    "        filtered_mask = np.zeros(labels.shape, dtype=int)\n",
    "        # Iterate through regions in case there is more than one\n",
    "        # larger than the minimum size\n",
    "        minr_list, minc_list, maxr_list, maxc_list = [], [], [], []\n",
    "        for region in regions:\n",
    "            if region.area > min_size:\n",
    "                filtered_mask[(labels == region.label)] = 1\n",
    "                minr_list.append(region.bbox[0]) \n",
    "                minc_list.append(region.bbox[1]) \n",
    "                maxr_list.append(region.bbox[2]) \n",
    "                maxc_list.append(region.bbox[3])\n",
    "            # Find the true min and max bounding rows and cols\n",
    "            minr = min(minr_list)\n",
    "            minc = min(minc_list)\n",
    "            maxr = max(maxr_list)\n",
    "            maxc = max(maxc_list)\n",
    "        bbox_list.append((minr, minc, maxr, maxc))\n",
    "    return bbox_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f9f1042",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "bbox_list = gen_bbox_list(filter_size_img_list)\n",
    "\n",
    "fig, axes = plot.plot_bbox(\n",
    "    filter_size_img_list,\n",
    "    bbox_list, \n",
    "    row_title='Region bounds on mask'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b66b981c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plot.plot_bbox(\n",
    "    [frames[1], frames[4], frames[7]],\n",
    "    bbox_list, \n",
    "    row_title='Region bounds on raw data'\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}