3_-_Peakfinding_for_ring_widths.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ecfab0d4",
   "metadata": {},
   "source": [
    "## Find the peaks of the extracted growth line pixels\n",
    "\n",
    "The final stage of the process is to extract the peaks of the pixel intensities\n",
    "and measure the distance between them, and, because the peaks\n",
    "should correspond to the segmented rings, this is analogous to hand-measuring\n",
    "the distance between rings. The peaks are found via the use of the\n",
    "[`scipy.signal.find_peaks`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html)\n",
    "function. It relies on the notion of\n",
    "[peak prominence](https://en.wikipedia.org/wiki/Topographic_prominence) to\n",
    "determine which peaks in the signal, i.e., the vector of pixel intensities, are\n",
    "really *peaks* in the signal and which are just noise. Full details can be\n",
    "found in the project [README](./README.md)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0972f1ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import yaml\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.signal import find_peaks\n",
    "from shellai import peakfinding\n",
    "\n",
    "# load the config file\n",
    "with open(\"config.yaml\", \"r\") as fd:\n",
    "    cfg = yaml.safe_load(fd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6ef9cb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # #  General settings\n",
    "# base directory where we'll be saving things to\n",
    "dir_base = cfg['paths']['local']['base']\n",
    "\n",
    "# directory where the predictions are stored\n",
    "dir_predictions = \"output\"\n",
    "\n",
    "# directory where the ring widths and plots will be saved\n",
    "dir_saved_widths = \"ring_predictions\"\n",
    "\n",
    "# settings for both sets of predictions\n",
    "experiment_info = {\n",
    "    \"loo\": {\n",
    "        \"image_names\": cfg['images']['train_images'],\n",
    "        \"file_mask\": \"loo_{image_name:s}_all_predictions.npz\",\n",
    "        \"base_image_dir\": cfg['images']['train_dir']\n",
    "    },\n",
    "\n",
    "    \"final\": {\n",
    "        \"image_names\": cfg['images']['test_images'],\n",
    "        \"file_mask\": \"final_{image_name:s}_all_predictions.npz\",\n",
    "        \"base_image_dir\": cfg['images']['test_dir']\n",
    "    }\n",
    "}\n",
    "\n",
    "# location to a spreadsheet which contains the scale of the images\n",
    "# specifically, it contains the length in micrometers (μm) and its\n",
    "# corresponding length in pixels. see the included file for its structure:\n",
    "# pixel_absolute_conversions.xlsx\n",
    "conversion_spreadsheet_location = os.path.join(\n",
    "    dir_base,\n",
    "    \"pixel_absolute_conversions.xlsx\",\n",
    ")\n",
    "\n",
    "# whether to also plot the results\n",
    "create_plots = True\n",
    "\n",
    "# # # Peak finding settings\n",
    "# list of epochs we wish to get the peaks of the extracted growth line from\n",
    "epochs = [100, 150, 200]\n",
    "\n",
    "# for full details of how the peak finding works, and what the following\n",
    "# parameters do, please see:\n",
    "# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html\n",
    "\n",
    "# list of prominences to extract for\n",
    "peak_prominences = [0.05, 0.1]\n",
    "\n",
    "# minimum width of peaks\n",
    "peak_width = 1\n",
    "\n",
    "# minimum distance between peaks (in pixels). this value was selected by\n",
    "# looking at the minimum distance between the peaks training images (~30) and\n",
    "# choosing a smaller value (25)\n",
    "peak_min_dist = 25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fd43d3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# first, we extract the scaling from pixel length to micrometer\n",
    "# this uses the first row of the spreadsheet as column names in the dataframe\n",
    "df = pd.read_excel(conversion_spreadsheet_location)\n",
    "\n",
    "pix_to_um = {}\n",
    "\n",
    "# iterate through each row of the dataframe\n",
    "for _, row in df.iterrows():\n",
    "\n",
    "    # get the image name\n",
    "    image_name = str(row[\"Source\"])\n",
    "\n",
    "    # extract the value by which to multiply pixel distances to convert to um\n",
    "    scale = row[\"Length (um)\"] / row[\"Length (pixels)\"] # type:ignore\n",
    "    pix_to_um[image_name] = scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6438c39",
   "metadata": {},
   "outputs": [],
   "source": [
    "# gather all the experiments to extract the rings for\n",
    "exps = [(e, p) for e in epochs for p in peak_prominences]\n",
    "n_exps = len(exps)\n",
    "\n",
    "for experiment_name in experiment_info:\n",
    "    image_names = experiment_info[experiment_name][\"image_names\"]\n",
    "    file_mask = experiment_info[experiment_name][\"file_mask\"]\n",
    "    base_image_dir = experiment_info[experiment_name][\"base_image_dir\"]\n",
    "\n",
    "    for image_name in image_names:\n",
    "        pred_path = os.path.join(\n",
    "            dir_base, \n",
    "            dir_predictions,\n",
    "            file_mask.format(image_name=image_name)\n",
    "        )\n",
    "\n",
    "        # get the measurements along the growth line (ring_data_ and their\n",
    "        # corresponding pixel coordinates (flc)\n",
    "        with np.load(pred_path, allow_pickle=True) as fd:\n",
    "            # dictionary keys are ['preds', 'gt']\n",
    "            ring_data = fd[\"ring_data\"].item()\n",
    "\n",
    "            # coordinates of the line along which the ring data was extracted\n",
    "            # this should be shape (n, 2)\n",
    "            flc = fd[\"flc\"]\n",
    "\n",
    "        # scaling to convert px to um\n",
    "        px_to_um_scaling = pix_to_um[image_name]\n",
    "\n",
    "        base_image_dir = experiment_info[experiment_name][\"base_image_dir\"]\n",
    "        image_dir = os.path.join(base_image_dir, image_name)\n",
    "\n",
    "        # ground-truth distances in um. these are hand-measured\n",
    "        gt_ring_distances_um = peakfinding.parse_measurement_csv_from_directory(\n",
    "            image_dir\n",
    "        )\n",
    "\n",
    "        # storage for the extracted ring widths and ground truth\n",
    "        results = {}\n",
    "        results['original'] = gt_ring_distances_um\n",
    "\n",
    "        # for combination of epoch and prominence value\n",
    "        for (epoch, prominence) in exps:\n",
    "            # get the pixel values along the axis\n",
    "            r = ring_data['preds'][epoch]\n",
    "\n",
    "            # normalise them such that the largest value is one and the\n",
    "            # smallest is zero. note that this is done so that the prominence\n",
    "            # value, has as close to the same meaning for each set of pixels\n",
    "            rnormed = peakfinding.normalize_vector(r)\n",
    "\n",
    "            # find the peaks\n",
    "            r_peaks, _ = find_peaks(\n",
    "                rnormed, \n",
    "                prominence=prominence, \n",
    "                width=peak_width,\n",
    "                distance=peak_min_dist,\n",
    "            )\n",
    "\n",
    "            # calculate the distances between the extracted peaks,\n",
    "            # rescale them, and reverse the order. the latter is carried out\n",
    "            # to match the ground truth\n",
    "            predicted_ring_distances_um = peakfinding.get_predicted_ring_distances(\n",
    "                flc, r_peaks, px_to_um_scaling, reversed=True\n",
    "            )\n",
    "            \n",
    "            # store the resulting vector of ring widths\n",
    "            results[(epoch, prominence)] = predicted_ring_distances_um\n",
    "\n",
    "        # create a pandas dataframe containing the ring widths by starting\n",
    "        # with the ground truth in one column\n",
    "        df = pd.DataFrame(results['original'], columns=['Ground Truth'])\n",
    "\n",
    "        # and adding columns corresponding to each prominence and epoch combo\n",
    "        for key in exps:\n",
    "            df = pd.concat([df, pd.DataFrame(results[key], columns=[key])], axis=1)\n",
    "\n",
    "        # next, rename the index appropriately\n",
    "        df = df.rename_axis('Ring Number')\n",
    "\n",
    "        # and finally, save the csv file\n",
    "        csv_file_path = os.path.join(\n",
    "            dir_base, \n",
    "            dir_predictions, \n",
    "            dir_saved_widths,\n",
    "            f\"{experiment_name}_{image_name}_rings.csv\", \n",
    "        )\n",
    "        df.to_csv(csv_file_path, float_format=\"%0.2f\")\n",
    "\n",
    "        print(f'Saved: {csv_file_path}')\n",
    "\n",
    "        # if set to True, we create a plot of each experiment for a given image\n",
    "        if create_plots:\n",
    "            image_file_path = os.path.join(\n",
    "                dir_base, \n",
    "                dir_predictions, \n",
    "                dir_saved_widths,\n",
    "                f\"{experiment_name}_{image_name}_rings.pdf\", \n",
    "            )\n",
    "\n",
    "            gt = results['original']\n",
    "            gtmax = np.max(results['original'])\n",
    "            maxlen = np.max([len(x) for x in results.values()])\n",
    "\n",
    "            fig, axes = plt.subplots(\n",
    "                n_exps, 1, figsize=(15, 3 * n_exps), sharex=True, dpi=200\n",
    "            )\n",
    "            for i, (epoch, prominence) in enumerate(exps):\n",
    "                ax = axes[i]\n",
    "                lbl = f\"Epochs={epoch}, Prominence={prominence}\"\n",
    "                pred = results[(epoch, prominence)]\n",
    "\n",
    "                gtx = np.arange(maxlen - gt.size, maxlen)\n",
    "                gtp = np.arange(maxlen - pred.size, maxlen)\n",
    "\n",
    "                ax.plot(gtx, gt, 'o-', label='Ground truth')\n",
    "                ax.plot(gtp, pred, 'x-', label=lbl)\n",
    "                ax.legend(loc=\"upper right\")\n",
    "\n",
    "                ax.set_ylim([0, gtmax+ 50])\n",
    "\n",
    "            plt.savefig(image_file_path, bbox_inches='tight')\n",
    "            plt.close()\n",
    "\n",
    "            print(f'Saved: {image_file_path}')"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}