ridges_er_stress.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e974507-2508-4827-8ce3-620ec97ed516",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import dependencies\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage.transform import rescale, resize\n",
    "from skimage.io import imread, imsave\n",
    "from skimage.filters import meijering, gaussian, threshold_otsu\n",
    "from skimage.morphology import opening, closing, disk, remove_small_objects\n",
    "from skimage.morphology import convex_hull_image, binary_closing\n",
    "from skimage.measure import label, regionprops\n",
    "from scipy import ndimage as ndi\n",
    "\n",
    "import math\n",
    "\n",
    "from glob import glob\n",
    "import os\n",
    "\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38787688-e63d-4564-ba1c-db6945a9539f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def filter_ridges(ridge, \n",
    "                  threshold_ridges=0.1, \n",
    "                  ridge_opening_element=disk(3)):\n",
    "    \n",
    "    # Threshold the ridge image to only retain values higher than a threshold, and to binarize the image\n",
    "    binary_ridge = ridge.copy()\n",
    "    binary_ridge[binary_ridge > threshold_ridges] = 1\n",
    "    binary_ridge[binary_ridge < 1] = 0\n",
    "\n",
    "    # Morphological opening and removing of small objects\n",
    "    opened_ridge = opening(binary_ridge, ridge_opening_element)\n",
    "    opened_ridge[~remove_small_objects(opened_ridge > 0, min_size=128)] = 0\n",
    "    \n",
    "    return binary_ridge, opened_ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36a3dccf-1596-4383-9e44-fa645603e011",
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_roi(img):\n",
    "    # Downscale the image\n",
    "    img_small = rescale(img, 0.25, anti_aliasing=False)\n",
    "    \n",
    "    # Threshold the image with a fixed threshold\n",
    "    img_smooth = img_small.copy()\n",
    "    \n",
    "    # Smooth the image with a Gaussian filter\n",
    "    smooth = gaussian(img_smooth, sigma=8)\n",
    "    smooth = gaussian(smooth, sigma=10)\n",
    "    \n",
    "    # Threshold with an Otsu threshold\n",
    "    thresh_value = threshold_otsu(smooth)\n",
    "    img_thresh_otsu = smooth > thresh_value\n",
    "\n",
    "    # Filling holes and convex hull\n",
    "    img_fill = ndi.binary_fill_holes(img_thresh_otsu)\n",
    "    img_fill1 = convex_hull_image(img_fill)\n",
    "    \n",
    "    # Create a label image\n",
    "    label_image = label(img_fill1)\n",
    "    min_area = 5000\n",
    "\n",
    "    filtered_labels = np.zeros(label_image.shape)\n",
    "    \n",
    "    img_regionprops = regionprops(label_image)\n",
    "\n",
    "    # Only keep regions with areas that are large enough\n",
    "    for index, region in enumerate(img_regionprops):\n",
    "        if region.area >= min_area:\n",
    "            print(index)\n",
    "            filtered_labels[label_image == (index + 1)] = 1\n",
    "    \n",
    "    largest_area = 0\n",
    "    largest_area_index = 0\n",
    "    \n",
    "    if np.max(filtered_labels) == 0:\n",
    "        for index, region in enumerate(img_regionprops):\n",
    "            # take regions with large enough areas\n",
    "            if region.area >= largest_area:\n",
    "                print(index)\n",
    "                largest_area_index = index\n",
    "                largest_area = region.area\n",
    "        filtered_labels[label_image == (largest_area_index + 1)] = 1\n",
    "    \n",
    "    filtered_labels_large = resize(filtered_labels, img.shape, anti_aliasing=False)\n",
    "    \n",
    "    return filtered_labels_large"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4dc057c-ec18-47e5-88f4-1dcb6f62d3e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract ROI and ridges\n",
    "\n",
    "img_path = '/path_to_data/'\n",
    "\n",
    "save_path = os.path.join('output_dir')\n",
    "\n",
    "img_files: list = glob(os.path.join(img_path, '*.tif'))\n",
    "\n",
    "ridge_densities = []\n",
    "img_name_list = []\n",
    "\n",
    "for img_file in img_files:\n",
    "    print(img_file)\n",
    "    img = imread(img_file)\n",
    "\n",
    "    if len(img.shape) == 3:\n",
    "        img = img[1, :, :]\n",
    "    \n",
    "    # Extract region of interest\n",
    "    filtered_labels_large = extract_roi(img)\n",
    "    \n",
    "    # Calculate ridges\n",
    "    ridge = meijering(img, sigmas=range(1, 10, 2), black_ridges=False)\n",
    "    \n",
    "    # Filter ridges\n",
    "    binary_ridge, opened_ridge = filter_ridges(ridge,\n",
    "                                               threshold_ridges=0.25,\n",
    "                                               ridge_opening_element=disk(2))\n",
    "    \n",
    "    # Filter the ridges with the ROI\n",
    "    filtered_ridge = opened_ridge.copy()\n",
    "    filtered_ridge[filtered_labels_large == 0] = 0\n",
    "    \n",
    "    # Calculate the ridge density in the ROI\n",
    "    ridge_density = sum(sum(filtered_ridge > 0.0001))/sum(sum(filtered_labels_large))\n",
    "    \n",
    "    # For display\n",
    "    filtered_ridge[filtered_ridge > 0.0001] = 1\n",
    "    filtered_ridge[filtered_ridge < 1] = 0\n",
    "    \n",
    "    filtered_ridge[filtered_labels_large > 0] += 0.3\n",
    "    \n",
    "    fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(16, 16))\n",
    "\n",
    "    ax[0, 0].imshow(ridge, cmap='gray')\n",
    "    ax[0, 0].set_title(r'Meijering ridges', fontsize=20)\n",
    "    ax[0, 1].imshow(binary_ridge, cmap='gray')\n",
    "    ax[0, 1].set_title(r'Thresholded ridges', fontsize=20)\n",
    "    ax[1, 0].imshow(opened_ridge, cmap='gray')\n",
    "    ax[1, 0].set_title(r'Opened ridges', fontsize=20)\n",
    "    ax[1, 1].imshow(filtered_ridge, cmap='gray')\n",
    "    ax[1, 1].set_title(r'Filtered ridges on filtered ROI', fontsize=20)\n",
    "    \n",
    "    for a in ax:\n",
    "        for a_2 in a:\n",
    "            a_2.set_axis_off()\n",
    "    \n",
    "    img_identifier = str(os.path.basename(img_file).split('.tif')[0])\n",
    "    plt.savefig(os.path.join(save_path, 'ridges_and_roi_steps_' + img_identifier))\n",
    "    \n",
    "    imsave(os.path.join(save_path, 'ridges_and_roi_' + img_identifier + '.png'), filtered_ridge)\n",
    "    \n",
    "    ridge_densities.append(ridge_density)\n",
    "    img_name_list.append(img_identifier)\n",
    "    \n",
    "density_dict = {'image': img_name_list,\n",
    "                'Density': ridge_densities}  \n",
    "       \n",
    "df = pd.DataFrame(density_dict) \n",
    "    \n",
    "# Saving the dataframe \n",
    "df.to_csv(os.path.join(save_path, 'densities.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3655067d-9d71-49d1-b94e-11d1b366ccac",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}