plot-residual-errors-recursive-checker.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Residual error automatic checker\n",
    "This notebook recursively find all residual errors files (default: ResMS.hdr|nii) and computes various descriptive stats measures and raise an alert if a measure is reaching above a preset threshold.\n",
    "\n",
    "Note that all stats and alerts are on non-zeros values only, in other words you should not count the quantity of 0 values in your calculations (also don't be afraid if you don't see any 0 value in histograms, that's normal, the count of 0 values are provided in descriptive stats for your information).\n",
    "\n",
    "By Stephen Larroque from Coma Science Group, University of Liège, created on 2018-01-27.\n",
    "\n",
    "Version v1.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "# BEWARE: autoreload works on functions and on general code, but NOT on new class methods:\n",
    "# if you add or change the name of a method, you have to reload the kernel!\n",
    "# also it will fail if you use super() calls in the classes you change\n",
    "# ALSO AUTORELOAD SHOULD BE THE FIRST LINE EVER EXECUTED IN YOUR IPYTHON NOTEBOOK!!!\n",
    "\n",
    "# Profilers:\n",
    "# http://pynash.org/2013/03/06/timing-and-profiling/\n",
    "# http://mortada.net/easily-profile-python-code-in-jupyter.html\n",
    "# use %lprun -m module func(*args, **kwargs)\n",
    "try:\n",
    "    %load_ext line_profiler\n",
    "    %load_ext memory_profiler\n",
    "except ImportError as exc:\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate figure inside IPython Notebook (must be called before any import of matplotlib, direct or indirect!)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pprint\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import nibabel as nib\n",
    "import os\n",
    "import re\n",
    "import scipy.stats\n",
    "from nilearn import image\n",
    "from nilearn import plotting\n",
    "#try:\n",
    "#    from niwidgets import NiftiWidget\n",
    "#    niwidgets_available = True\n",
    "#except Exception:\n",
    "#    niwidgets_available = False\n",
    "#    pass\n",
    "\n",
    "try:\n",
    "    from scandir import walk # use the faster scandir module if available (Python >= 3.5), see https://github.com/benhoyt/scandir\n",
    "except ImportError as exc:\n",
    "    from os import walk # else, default to os.walk()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PARAMETERS - EDIT ME\n",
    "rootpath = r'Z:\\some_folder'  # path from where to start recursively looking for residual error files\n",
    "RE_resms = 'ResMS\\.(img|nii)'  # regular expression to match residual error filenames\n",
    "alert_thresholds = {'max': 0.2, 'mean': 0.05, 'pct99': 0.1}  # alert thresholds, can use any descriptive stat calculated below. Must be formatted as a dictionary with 'parameter': threshold-value\n",
    "verbose = False  # in verbose mode, descriptive stats and histograms of all files will be displayed, whether it raises an alert or not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def recwalk(inputpath, sorting=True, folders=False, topdown=True, filetype=None, regex=None):\n",
    "    '''Recursively walk through a folder. This provides a mean to flatten out the files restitution (necessary to show a progress bar). This is a generator.'''\n",
    "    if filetype and isinstance(filetype, list):\n",
    "        filetype = tuple(filetype)  # str.endswith() only accepts a tuple, not a list\n",
    "    if regex:\n",
    "        RE_precomp = re.compile(regex, re.I)\n",
    "    # If it's only a single file, return this single file\n",
    "    if os.path.isfile(inputpath):\n",
    "        abs_path = fullpath(inputpath)\n",
    "        yield os.path.dirname(abs_path), os.path.basename(abs_path)\n",
    "    # Else if it's a folder, walk recursively and return every files\n",
    "    else:\n",
    "        for dirpath, dirs, files in walk(inputpath, topdown=topdown):\n",
    "            if sorting:\n",
    "                files.sort()\n",
    "                dirs.sort()  # sort directories in-place for ordered recursive walking\n",
    "            # return each file\n",
    "            for filename in files:\n",
    "                if (not filetype or filename.endswith(filetype)) and (not regex or RE_precomp.search(filename)):\n",
    "                    yield (dirpath, filename)  # return directory (full path) and filename\n",
    "            # return each directory\n",
    "            if folders:\n",
    "                for folder in dirs:\n",
    "                    yield (dirpath, folder)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "voxel_threshold = 0.0001 # minimum threshold to consider as a voxel and not just background noise (because background voxels can be 0.000001 for example), can be float or str ('1%' to give a percentage). TODO: autodetect minimum value (can be -4, 0.02, etc) as the background and use it as the threshold value.\n",
    "for dirpath, filename in recwalk(rootpath, regex=RE_resms):\n",
    "    # Get fullpath of image\n",
    "    fpath = os.path.join(dirpath, filename)\n",
    "    # Load image data\n",
    "    im = image.load_img(fpath)\n",
    "    # Remove background noise voxels    \n",
    "    im = image.threshold_img(im, voxel_threshold)\n",
    "    imdata = im.get_data()\n",
    "    # Calculate descriptive stats\n",
    "    imdata_nz = imdata[imdata != 0]  # only on non-zeros values\n",
    "    desc_stats = {'mean': imdata_nz.mean(),\n",
    "                  'median': np.median(imdata_nz),\n",
    "                  'min': imdata_nz.min(),\n",
    "                  'max': imdata_nz.max(),\n",
    "                  'range': np.ptp(imdata_nz),\n",
    "                  'pct25': np.percentile(imdata_nz,25,interpolation='lower'),\n",
    "                  'pct50': np.percentile(imdata_nz,50,interpolation='lower'),\n",
    "                  'pct75': np.percentile(imdata_nz,75,interpolation='lower'),\n",
    "                  'pct90': np.percentile(imdata_nz,90,interpolation='lower'),\n",
    "                  'pct99': np.percentile(imdata_nz,99,interpolation='lower'),\n",
    "                  'iqr': scipy.stats.iqr(imdata_nz, rng=(25, 75), interpolation='lower'),\n",
    "                  'var': np.var(imdata_nz),\n",
    "                  'stdev': np.std(imdata_nz),\n",
    "                  'skew': scipy.stats.skew(imdata_nz),\n",
    "                  'kurtosis': scipy.stats.kurtosis(imdata_nz, fisher=True),\n",
    "                  'count_nonzero': np.size(imdata_nz),\n",
    "                  'count': np.size(imdata),\n",
    "                  'count_zero': np.size(imdata[imdata == 0]),\n",
    "                  }\n",
    "    # Check if we reached aboved alert thresholds on any descriptive stats feature\n",
    "    alert = False\n",
    "    if alert_thresholds is not None:\n",
    "        for feature, threshold in alert_thresholds.items():\n",
    "            if desc_stats[feature] > threshold:\n",
    "                alert = True\n",
    "                break\n",
    "    # If the alert is raised or we are in verbose mode, display infos about this file\n",
    "    if alert or verbose:\n",
    "        # Display which file raised the alert (useful for user to find the incriminated file)\n",
    "        if alert:\n",
    "            print('Alert! Residual errors are above thresholds in file: %s' % fpath)\n",
    "        else:\n",
    "            print('Displaying stats for residual error file: %s' % fpath)\n",
    "        # Print descriptive stats\n",
    "        pprint.pprint(desc_stats)\n",
    "        # Show histogram\n",
    "        frq, edges = np.histogram(imdata_nz, bins=100)\n",
    "        fig, ax = plt.subplots()\n",
    "        ax.bar(edges[:-1], frq, width=np.diff(edges), ec=\"k\", align=\"edge\")\n",
    "        plt.xlabel('Residual error value voxel-wise')\n",
    "        plt.ylabel('Occurrences')\n",
    "        plt.show()\n",
    "        # Show spatial localization on brain of the residual errors\n",
    "        # If niwidgets is available, show an interactive visualization, else just show nilearn static visualization of the max value\n",
    "        #if niwidgets_available:\n",
    "            #my_widget = NiftiWidget(fpath)\n",
    "            #my_widget.nifti_plotter()\n",
    "        #else:\n",
    "        fig2 = plotting.plot_glass_brain(fpath, colorbar=True)\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}