plot_project_overlapping_masks_to_atlas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project overlapping masks to atlas\n",
    "This notebook project a list of masks (generated from CONN for example or SPM), then projects to an atlas, then extract overlapping regions and plot these regions (and generate a mask so that you can use another application to generate a nice visualization).\n",
    "\n",
    "By Stephen Larroque from Coma Science Group, University of Liège, created on 2017-04-18.\n",
    "\n",
    "Version v1.6.6\n",
    "\n",
    "TODO:\n",
    "* generate maps for all x top overlaps from 1 to max if min_threshold is set to 'all' instead of int."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as pltcol\n",
    "import numpy as np\n",
    "import nibabel as nib\n",
    "from nilearn import image\n",
    "from nilearn import plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "voxel_threshold = 0.0001 # minimum value to be considered a non-background voxel signal (because background voxels not just 0.0, can be 0.0000000000001), 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",
    "dpi_resolution = 300\n",
    "# List of masks\n",
    "# NOTE: the first image will be used as the template to resample other masks!\n",
    "list_imgs = [\n",
    "            r'hyperc2019\\spmresults\\dex\\dmn-dex-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\doc-uws\\dmn-doc-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\doc-mcs\\dmn-docmcs-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\doc-emcs\\dmn-docemcs-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\doc-uws-acute\\dmn-docacute-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\doc-mcs-acute\\dmn-docacutemcs-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\ket\\dmn-ket-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\prop\\dmn-prop-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\sevo\\dmn-sevo-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\sleep\\dmn-sleep-hyperc.nii',\n",
    "            r'hyperc2019\\spmresults\\sleep-n2n3\\dmn-sleepn2n3-hyperc.nii',\n",
    "            ]\n",
    "#list_imgs = [\n",
    "#            r'hyperc2019\\spmresults\\dex-allsess\\dmn-dex-hypercgroup.nii',\n",
    "#            r'hyperc2019\\spmresults\\doc-allsess\\dmn-doc-hypercgroup.nii',\n",
    "#            r'hyperc2019\\spmresults\\doc-acute-allsess\\dmn-docacute-hypercgroup.nii',\n",
    "#            r'hyperc2019\\spmresults\\ket-allsess\\dmn-ket-hypercgroup.nii',\n",
    "#            r'hyperc2019\\spmresults\\prop-allsess\\dmn-prop-hypercgroup.nii',\n",
    "#            r'hyperc2019\\spmresults\\sevo-allsess\\dmn-sevo-hypercgroup.nii',\n",
    "#            r'hyperc2019\\spmresults\\sleep-allsess\\dmn-sleep-hypercgroup.nii',\n",
    "#            ]\n",
    "# Optional: name for each mask, will be used to compute the top overlapping regions table at the end\n",
    "list_imgs_names = [\n",
    "    'Dex',\n",
    "    'UWS',\n",
    "    'MCS',\n",
    "    'EMCS',\n",
    "    'UWS acute',\n",
    "    'MCS acute',\n",
    "    'Ketamine',\n",
    "    'Propofol',\n",
    "    'Sevo',\n",
    "    'Sleep N3-Awake',\n",
    "    'Sleep N3-N2',\n",
    "]\n",
    "#list_imgs_names = [\n",
    "#    'Dex',\n",
    "#    'DOC',\n",
    "#    'DOC acute',\n",
    "#    'Ketamine',\n",
    "#    'Propofol',\n",
    "#    'Sevo',\n",
    "#    'Sleep N3',\n",
    "#]\n",
    "\n",
    "# Load masks and resample to first\n",
    "imgs = []\n",
    "for img in list_imgs:\n",
    "    im = image.load_img(img)\n",
    "    if imgs:\n",
    "        if im.shape != imgs[0].shape:\n",
    "            im = image.resample_to_img(im, imgs[0])\n",
    "    im = image.threshold_img(im, voxel_threshold)\n",
    "    imgs.append(im)\n",
    "# PLot!\n",
    "plotting.plot_prob_atlas(imgs, view_type=\"filled_contours\",\n",
    "                    title=\"lala\", colorbar=True, cut_coords=(0,0,0), draw_cross=True, cmap=pltcol.ListedColormap(['b', 'g', 'r', 'c', 'm'], name='from_list', N=None))\n",
    "plotting.plot_roi(imgs[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import plotting, datasets\n",
    "# Atlas\n",
    "atlas_choice = 'aal2' # anatomytoolbox or aal2\n",
    "if atlas_choice == 'anatomytoolbox':\n",
    "    atlas_path = 'masks\\AnatomyToolbox_Atlas_Map.nii'  # TODO: build atlas variable with all infos and data (labels, indices, nib niftiimage with affine etc)\n",
    "else:\n",
    "    atlas = datasets.fetch_atlas_aal(version='SPM12', data_dir='atlas')\n",
    "# Load mask images\n",
    "imgs = []\n",
    "for img in list_imgs:\n",
    "    im = image.load_img(img)\n",
    "    if imgs: # if there is already at least one image loaded, resample subsequent images to the first one\n",
    "        if im.shape != imgs[0].shape:\n",
    "            im = image.resample_to_img(im, imgs[0])\n",
    "    im = image.threshold_img(im, voxel_threshold)\n",
    "    imgs.append(im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "atlas_im = image.load_img(atlas.maps)\n",
    "print('Atlas shape: %s' % str(atlas_im.shape))\n",
    "nb_regions = len(np.unique(atlas_im.get_data()))- 1\n",
    "print('%i regions in this atlas: %s' % (nb_regions, str(np.unique(atlas_im.get_data())))) # 48 regions because 0 is background\n",
    "print('%i labels' % len(atlas['labels']))\n",
    "print('%i indices: %s' % (len(atlas['indices']), atlas['indices']))\n",
    "atlas.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Number of non-zero voxels per map (after thresholding):')\n",
    "for im in imgs:\n",
    "    print(np.nonzero(im.get_data())[0].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just plot each map\n",
    "for im in imgs:\n",
    "    plotting.plot_stat_map(im)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_atlas_label(atlas, region_idx):\n",
    "    return atlas['labels'][atlas['indices'].index(str(region_idx))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Resample masks to atlas size\n",
    "imgs2 = []\n",
    "for img in imgs:\n",
    "    if img.shape != atlas_im.shape:\n",
    "        img = image.resample_to_img(img, atlas_im)\n",
    "    img = image.threshold_img(img, voxel_threshold)\n",
    "    imgs2.append(img)\n",
    "imgs = imgs2\n",
    "del imgs2\n",
    "imgs[0].shape\n",
    "\n",
    "# Extract activated atlas brain regions for each mask\n",
    "maps_regions = []\n",
    "maps_regions_idxs = []\n",
    "maps_regions_count = []\n",
    "for img in imgs:\n",
    "    # Extract only non zeros voxels indices from mask\n",
    "    im_data = img.get_data()\n",
    "    #np.extract(im_data>0, im_data)\n",
    "    vox_thres = np.nonzero(im_data)\n",
    "    # Compare with atlas regions to extract region indices\n",
    "    atlas_data = atlas_im.get_data()\n",
    "    region_indices = set()\n",
    "    region_count = {}\n",
    "    for x in zip(*vox_thres): # walk through all non zero voxels of mask\n",
    "        region_idx = atlas_data[x] # get equivalent voxel from atlas\n",
    "        if region_idx != 0: # if not background\n",
    "            # Append region index into the set (so that they are unique)\n",
    "            region_indices.add(region_idx)\n",
    "            # Increase the count of voxels activated in this region\n",
    "            region_label = get_atlas_label(atlas, region_idx)\n",
    "            if region_label not in region_count:\n",
    "                region_count[region_label] = 0\n",
    "            region_count[region_label] += 1\n",
    "    if 0 in region_indices:\n",
    "        region_indices.remove(0) # remove background, not part of the atlas labels\n",
    "    print('Atlas indices of brain regions activated in current mask: %s' % str(sorted(region_indices)))\n",
    "    # Extract brain region names from atlas that are present in this mask\n",
    "    matching_idxs = [int(idx) in region_indices for idx in atlas['indices']]\n",
    "    map_brain_regions = filter(None, [label if match else None for label, match in zip(atlas['labels'], matching_idxs)])\n",
    "    maps_regions.append(map_brain_regions)\n",
    "    maps_regions_idxs.append(region_indices)\n",
    "    maps_regions_count.append(region_count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "maps_regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display number of voxels for each region of each mask\n",
    "from collections import OrderedDict\n",
    "for i, region_count in enumerate(maps_regions_count):\n",
    "    print('== %i: %s' % (i, str(OrderedDict(sorted(region_count.items(), key=lambda x: x[1], reverse=True)))) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Compute count of voxels per atlas region\n",
    "atlas_regions_count = {}\n",
    "for coord in zip(*np.where(atlas_data)):\n",
    "    region_idx = atlas_data[coord]\n",
    "    region_label = get_atlas_label(atlas, region_idx)\n",
    "    if not region_label in atlas_regions_count:\n",
    "        atlas_regions_count[region_label] = 0\n",
    "    atlas_regions_count[region_label] += 1\n",
    "atlas_regions_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Compute percentage of voxels (maps_regions_count / atlas_regions_count)\n",
    "maps_regions_percent = []\n",
    "for regions_count in maps_regions_count:\n",
    "    maps_regions_percent.append({region_label: (float(region_count) / atlas_regions_count[region_label] * 100) for region_label, region_count in regions_count.items()})\n",
    "\n",
    "# Display percentage\n",
    "from collections import OrderedDict\n",
    "for i, region_percent in enumerate(maps_regions_percent):\n",
    "    print('== %i: %s' % (i+1, str(OrderedDict(sorted(region_percent.items(), key=lambda x: x[1], reverse=True)))) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rank how much each brain regions are overlapping\n",
    "overlap_min_thres = 1 # len(imgs)-1 # minimum overlapping threshold = minimum number of maps that need to have these regions activated\n",
    "overlap_min_voxels = 5 # minimum number of voxels required to consider a region really activated and not spurious activity\n",
    "overlap_min_percent = 1 # minimum percentage of voxels activated over total atlas region surface to consider as real activation and not spurious\n",
    "\n",
    "overlaps = {}\n",
    "for region in atlas['labels']:\n",
    "    overlaps[region] = 0\n",
    "    for regions_mask, regions_count, regions_percent in zip(maps_regions, maps_regions_count, maps_regions_percent):\n",
    "        if region in regions_mask and regions_count[region] >= overlap_min_voxels and regions_percent[region] >= overlap_min_percent:\n",
    "            overlaps[region] += 1\n",
    "overlaps_top = {k:v for k,v in overlaps.items() if v >= overlap_min_thres}\n",
    "print(len(overlaps_top))\n",
    "overlaps_top"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get atlas indices for top overlapping regions\n",
    "overlap_condition = [True if label in overlaps_top.keys() else False for label in atlas['labels']]\n",
    "overlap_idxs = [int(idx) for idx,cond in zip(atlas['indices'], overlap_condition) if cond]\n",
    "overlap_idxs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct atlas constrained to only top overlapping regions\n",
    "# Extract overlapping regions voxels\n",
    "ix = np.in1d(atlas_data, overlap_idxs).reshape(atlas_data.shape)\n",
    "vox_overlap_atlas = np.where(ix)\n",
    "if len(vox_overlap_atlas[0]) == 0:\n",
    "    print('Found no overlapping voxel, cannot generate a map!')\n",
    "else:\n",
    "    # Generate atlas with indices/labels\n",
    "    overlap_data = np.zeros(atlas_data.shape, dtype='uint8')\n",
    "    for coord in zip(*vox_overlap_atlas):\n",
    "        overlap_data[coord] = int(atlas_data[coord])\n",
    "print(overlap_data.dtype)\n",
    "np.where(overlap_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate a full brain atlas to show in other vis softwares (so that there is no missing brain)\n",
    "overlap_data_full = overlap_data.copy()\n",
    "for coord in zip(*np.where(atlas_data)):\n",
    "    if overlap_data_full[coord] == 0:\n",
    "        overlap_data_full[coord] = 7 # TODO: change value here if some areas are label 7, we want to use here a label that does not exist! 7 is grey in MRIcroGL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Convert numpy mask to nifti image using nibabel\n",
    "overlap_image = nib.Nifti1Image(overlap_data, affine=atlas_im.affine)\n",
    "overlap_image_full = nib.Nifti1Image(overlap_data_full, affine=atlas_im.affine)\n",
    "#overlap_image.header['regular'] = 'r'\n",
    "# Set the header intent_code to 1002 to signify it's an atlas (label regions) to ease visualization in other softwares\n",
    "# https://nifti.nimh.nih.gov/pub/dist/src/niftilib/nifti1.h\n",
    "# https://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/group__NIFTI1__INTENT__CODES.html#a26\n",
    "overlap_image.header['intent_code'] = 1002\n",
    "overlap_image_full.header['intent_code'] = 1002\n",
    "# Save overlap image as a mask\n",
    "nib.save(overlap_image, 'overlap_image.nii')  # just the overlap regions\n",
    "nib.save(overlap_image_full, 'overlap_image_full.nii')  # overlap regions + rest of the atlas regions as region 0 (black)\n",
    "\n",
    "# Plot!\n",
    "print('Overlapping brain regions with min overlap=%i: %s' % (overlap_min_thres, overlaps_top))\n",
    "#plotting.plot_prob_atlas(overlap_data, view_type=\"filled_contours\",\n",
    "                    #title=\"lala\", colorbar=True, cut_coords=(0,0,0), draw_cross=True, cmap=pltcol.ListedColormap(['b', 'g', 'r', 'c', 'm'], name='from_list', N=None))\n",
    "fig1 = plotting.plot_roi(overlap_image, title='Overlapping brain regions with min overlap=%i' % overlap_min_thres, cut_coords=[5, 10, 37], cmap=plt.cm.prism)\n",
    "fig1.savefig('overlap_image.png', dpi=dpi_resolution)\n",
    "fig2 = plotting.plot_glass_brain(overlap_image, title='Overlapping brain regions with min overlap=%i' % overlap_min_thres, cmap=plt.cm.prism)\n",
    "fig2.savefig('overlap_image_glass.png', dpi=dpi_resolution)\n",
    "print('Saved results in overlap_image.nii, overlap_image_full.nii, overlap_image.png and overlap_image_glass.png')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute overlap table, which is a summary of whether each region is present or not in each map\n",
    "import pandas as pd\n",
    "disable_overlap_min_checks = False  # can disable the minimum amount of voxels required to count the region in\n",
    "overlap_table = pd.DataFrame(columns=list_imgs_names + ['total', 'real_total'], dtype='uint8')\n",
    "\n",
    "for region_label, total in overlaps_top.items():\n",
    "    overlap_table.loc[region_label, :] = False\n",
    "    for map_id, regions in enumerate(maps_regions):\n",
    "        if region_label in regions:\n",
    "            if disable_overlap_min_checks or (maps_regions_count[map_id][region_label] >= overlap_min_voxels and maps_regions_percent[map_id][region_label] >= overlap_min_percent):\n",
    "                overlap_table.loc[region_label, overlap_table.columns[map_id]] = True\n",
    "    overlap_table.loc[region_label, 'total'] = total\n",
    "    overlap_table.loc[region_label, 'real_total'] = overlap_table.loc[region_label, list_imgs_names].astype(int).sum()\n",
    "\n",
    "# Display results\n",
    "# Better table formatting + enable hovering line by line\n",
    "from IPython.display import HTML\n",
    "def hover(hover_color=\"#ffff99\"):\n",
    "    return dict(selector=\"tr:hover\",\n",
    "                props=[(\"background-color\", \"%s\" % hover_color)])\n",
    "styles = [\n",
    "    hover(),\n",
    "    dict(selector=\"table\", props=[(\"border\", \"none\"),\n",
    "                                  (\"border-collapse\", \"collapse\")]),\n",
    "    dict(selector=\"th\", props=[(\"font-size\", \"110%\"),\n",
    "                               (\"text-align\", \"center\"),\n",
    "                               (\"border\", \"none\")]),\n",
    "    dict(selector=\"caption\", props=[(\"caption-side\", \"top\"),\n",
    "                                    (\"text-align\", \"center\")]),\n",
    "    dict(selector=\"td\", props=[(\"text-align\", \"center\"),\n",
    "                               (\"border\", \"none\"),\n",
    "                               (\"font-size\", \"110%\")]),\n",
    "    dict(selector=\"tr\", props=[(\"border\", \"none\")])\n",
    "]\n",
    "\n",
    "def green(val):\n",
    "    # From https://stackoverflow.com/questions/41555678/highlighting-multiple-cells-in-different-colors-with-pandas#\n",
    "    if val is True:\n",
    "        color = 'lightgreen'\n",
    "        return 'background-color: %s' % color\n",
    "    else:\n",
    "        return ''  # need to return a string in any case\n",
    "\n",
    "def green_image(val):\n",
    "    # If true, hide the text and show a picture inplace\n",
    "    if val is True:\n",
    "        return 'width:32px;height:27px;background:url(\"checkgreen_small.png\");background-repeat:no-repeat;background-position:center;text-indent:-9999px'\n",
    "    else:\n",
    "        return ''  # need to return a string in any case\n",
    "\n",
    "overlap_table.sort_values(['real_total'], ascending=False, inplace=True)  # sort\n",
    "overlap_table = overlap_table.replace(False, '')  # empty false (cleaner table, more legible)\n",
    "html = (overlap_table.style.set_caption('Overlapping hyperconnectivity regions summary table')\\\n",
    "                .set_table_styles(styles)\\\n",
    "                .applymap(green_image))  # highlight True cells in green\n",
    "\n",
    "html  # display!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the table to a csv file\n",
    "overlap_table.to_csv('overlap_table.csv', sep=';')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----------------------------------------------\n",
    "## Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "overlap_image.header.values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aal = nib.load('aal.nii.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "aal.header['sizeof_hdr']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "overlap_image.header['regular'] = 'r'\n",
    "overlap_image.header['intent_code'] = 1002"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in aal.header.keys():\n",
    "    #if aal.header[k] != overlap_image.header[k]:\n",
    "    print('%s=%s OR %s' % (k, aal.header[k], overlap_image.header[k]))\n",
    "    #elif k not in overlap_image.header.keys():\n",
    "    #    print('MISSING: %s=%s' % (k, aal.header[k]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "overlaps"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}