2d-spot-detection-ij.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 2D Spot Detection with ImageJ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Discipline: Spot / object counting (SptCnt)\n",
    "\n",
    "Task: Estimate the number of objects\n",
    "\n",
    "This project illustrates the 2D counting of small vesicle like objects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The input images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The images were generated by [SIMCEP](http://www.cs.tut.fi/sgn/csb/simcep/tool.html), a widefield fluorescence microscopy biological images simulator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "First set paths to the input and ground-truth folders and create the output folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "IN_FOLDER = \"/home/jovyan/data/in/\"\n",
    "GT_FOLDER = \"/home/jovyan/data/gt/\"\n",
    "OUT_FOLDER = \"/home/jovyan/data/out/\"\n",
    "\n",
    "!mkdir -p $OUT_FOLDER"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The input images are in the [ome-tiff](https://docs.openmicroscopy.org/ome-model/6.0.1/ome-tiff/) format. We display the input images. If after execution of the cell the images do not show, run the cell a second time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import tifffile\n",
    "import matplotlib.pyplot\n",
    "files = [i.path for i in os.scandir(IN_FOLDER) if i.is_file()]\n",
    "in_paths_sources = [i for i in files if ('.tif' in i) ]\n",
    "sources = []\n",
    "for path in in_paths_sources:\n",
    "    source = tifffile.imread(path)\n",
    "    tifffile.imshow(source, cmap='gray')\n",
    "    sources.append(source)\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The image analysis workflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The ImageJ macro uses the Laplacian of Gaussian filter from [FeatureJ](https://imagej.net/FeatureJ) with the given standard deviation of the Gaussian derivative kernels and detects the minima with the given noise tolerance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "%%writefile  /fiji/macros/macro.ijm\n",
    "// Author: Sébastien Tosi (IRB Barcelona) \n",
    "// Version: 1.0\n",
    "// Date: 21/04/2017\n",
    "\n",
    "// The default input and output folder\n",
    "inputDir = \"/dockershare/667/in/\";\n",
    "outputDir = \"/dockershare/667/out/\";\n",
    "\n",
    "// Functional parameters\n",
    "LapRad = 2;\n",
    "NoiseTol = 2.5;\n",
    "\n",
    "arg = getArgument();\n",
    "parts = split(arg, \",\");\n",
    "\n",
    "setBatchMode(true);\n",
    "for(i=0; i<parts.length; i++) {\n",
    "\tnameAndValue = split(parts[i], \"=\");\n",
    "\tif (indexOf(nameAndValue[0], \"input\")>-1) inputDir=nameAndValue[1];\n",
    "\tif (indexOf(nameAndValue[0], \"output\")>-1) outputDir=nameAndValue[1];\n",
    "\tif (indexOf(nameAndValue[0], \"radius\")>-1) LapRad=nameAndValue[1];\n",
    "\tif (indexOf(nameAndValue[0], \"noise\")>-1) NoiseTol=nameAndValue[1];\n",
    "}\n",
    "\n",
    "images = getFileList(inputDir);\n",
    "\n",
    "for(i=0; i<images.length; i++) {\n",
    "\timage = images[i];\n",
    "\tif (endsWith(image, \".tif\")) {\n",
    "\t\t// Open image\n",
    "\t\topen(inputDir + \"/\" + image);\n",
    "\t\twidth = getWidth();\n",
    "\t\theight = getHeight();\n",
    "\t\t\n",
    "\t\t// Processing\n",
    "\t\trun(\"Clear Results\", \"\");\n",
    "\t\trun(\"FeatureJ Laplacian\", \"compute smoothing=\"+d2s(LapRad,2));\n",
    "\t\trun(\"Find Maxima...\", \"noise=\"+d2s(NoiseTol,2)+\" output=List light\");\n",
    "\t\t\n",
    "\t\t// Export results\n",
    "\t\tnewImage(\"Mask\", \"16-bit black\", width, height, 1);\n",
    "\t\tfor(r=0;r<nResults;r++)\n",
    "\t\t{\n",
    "\t\t\tXPos = getResult(\"X\",r);\n",
    "\t\t\tYPos = getResult(\"Y\",r);\n",
    "\t\t\tsetPixel(XPos,YPos,65535);\t\t\t\n",
    "\t\t}\n",
    "\t\tsave(outputDir + \"/\" + image);\n",
    "\t\t// Cleanup\n",
    "\t\trun(\"Close All\");\n",
    "\t}\n",
    "}\n",
    "run(\"Quit\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Running the workflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The parameters for running the workflow are:\n",
    "* --infolder: the folder containing the input images\n",
    "* --gtfolder: the folder containing the ground truth data, only needed for benchmarking\n",
    "* --outfolder: the folder into which the result images will be written\n",
    "\n",
    "It has 4 switches that control the usage of the BIAflows server:\n",
    "* --no_download (-nd): the images are not downloaded from the server, they must already be in the input folder\n",
    "* --no_annotations_upload (-nau): the workflow results are only written into the output folder, but not uploaded to the server\n",
    "* --no_metrics_computation (-nmc): the metrics are neither computed nor uploaded to the server\n",
    "* --no_metrics_upload (-nmu): the metrics are computed but not uploaded to the server\n",
    "\n",
    "The switch --local (-l) is equivalent to using the three switches --no_download, --no_annotations_upload and --no_metrics_upload."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "If the analysis workflow has parameters these must also be passed. Their names can be found at the end of the file [descriptor.json](https://github.com/Neubias-WG5/W_SpotDetection-IJ/blob/master/descriptor.json).\n",
    "\n",
    "The workflow parameters are:\n",
    "\n",
    "* STD_DEV: The standard deviation of the Gaussian derivative kernels used for computing the second-order partial derivatives of the Laplacian. Must be larger than zero. The value is realated to the size of the spots that should be detected. The bigger the spots are the bigger the STD_DEV should be. See also the [documentation of FeatureJ](https://imagescience.org/meijering/software/featurej/laplacian/).\n",
    "\n",
    "* NOISE: The noise tolerance for detecting the minima in the filtered image. Minima are ignored if they do not stand out from the surroundings by more than this value. See also the [ImageJ user guide](https://imagej.nih.gov/ij/docs/guide/146-29.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "STD_DEV = 2\n",
    "NOISE = 2.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "%%capture cap_out --no-stderr\n",
    "!python /app/wrapper.py  --ij_radius $STD_DEV --ij_noise $NOISE --infolder $IN_FOLDER --gtfolder $GT_FOLDER --outfolder $OUT_FOLDER --local"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case the metric is the relative error of the spot count for each image. The metric is written to standard-out by BIAFlows. We get it from standard-out and store it in a dictionary which has the name of the images has keys and the relative errors of the spot counts as values.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output = cap_out.stdout\n",
    "lines = output.split(\"\\n\")\n",
    "relativeError = {}\n",
    "for line in lines:\n",
    "    if len(line)>0 and line[0]==\">\":\n",
    "        line = line.replace(\">\", \"\")\n",
    "        key = line.split(\":\")[0].strip()\n",
    "        value = float(line.split(\":\")[2].replace(\"]\", \"\").strip())\n",
    "        relativeError[key] = value\n",
    "print(relativeError)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Displaying the results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We display the detected spots as red circles on the input images and count the number of spots per image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "from skimage.measure import label, regionprops\n",
    "from skimage.color import label2rgb\n",
    "\n",
    "files = [i.path for i in os.scandir(OUT_FOLDER) if i.is_file()]\n",
    "paths_sources = [i for i in files if ('.tif' in i) ]\n",
    "pointImages = []\n",
    "index = 0\n",
    "nrOfSpots = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "for path in paths_sources:\n",
    "    fig, ax = plt.subplots(figsize=(20, 12))\n",
    "    pointImage = tifffile.imread(path)\n",
    "    label_image = label(pointImage)\n",
    "    image_label_overlay = label2rgb(label_image, image=sources[index])\n",
    "    ax.imshow(sources[index], cmap=\"gray\")\n",
    "    spotCount = 0\n",
    "    for region in regionprops(label_image):\n",
    "        (x,y) = region.centroid\n",
    "        radius = region.equivalent_diameter/2.0\n",
    "        rect = mpatches.Circle((y,x), 5, fill=False, edgecolor='red')\n",
    "        ax.add_patch(rect)\n",
    "        spotCount = spotCount + 1\n",
    "    nrOfSpots.append(spotCount)\n",
    "        \n",
    "    index = index + 1\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    print(\"number of spots: \" + str(spotCount))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Finally we display the number of spots for each image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import statistics\n",
    "index = 0\n",
    "for path in in_paths_sources:\n",
    "    print(path + \"\\t\" + str(nrOfSpots[index]))\n",
    "    index = index + 1\n",
    "    \n",
    "print(\"Mean\" + \"\\t\" + str(statistics.mean(nrOfSpots)))\n",
    "print(\"StdDev\" + \"\\t\" + str(statistics.stdev(nrOfSpots)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And we calculate and display a boxplot of the distribution of the number of cells per image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.boxplot(nrOfSpots)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary of the metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have the relative error of the spot count per image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pprint\n",
    "pprint.pprint(relativeError)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculate the mean and standard-deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import statistics\n",
    "print(\"Relative Error\")\n",
    "print(\"Mean:\" + \"\\t\" + str(statistics.mean(relativeError.values())))\n",
    "print(\"StdDev:\" + \"\\t\" + str(statistics.stdev(relativeError.values())))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.boxplot(relativeError.values());"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}