02 - Figure 3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/oesteban/workspace/niworkflows/niworkflows/__init__.py:24: UserWarning: \n",
      "This call to matplotlib.use() has no effect because the backend has already\n",
      "been chosen; matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n",
      "or matplotlib.backends is imported for the first time.\n",
      "\n",
      "The backend was *originally* set to 'pgf' by the following code:\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/runpy.py\", line 193, in _run_module_as_main\n",
      "    \"__main__\", mod_spec)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/runpy.py\", line 85, in _run_code\n",
      "    exec(code, run_globals)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py\", line 16, in <module>\n",
      "    app.launch_new_instance()\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/traitlets/config/application.py\", line 658, in launch_instance\n",
      "    app.start()\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/kernelapp.py\", line 477, in start\n",
      "    ioloop.IOLoop.instance().start()\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/zmq/eventloop/ioloop.py\", line 177, in start\n",
      "    super(ZMQIOLoop, self).start()\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/tornado/ioloop.py\", line 888, in start\n",
      "    handler_func(fd_obj, events)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/tornado/stack_context.py\", line 277, in null_wrapper\n",
      "    return fn(*args, **kwargs)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py\", line 440, in _handle_events\n",
      "    self._handle_recv()\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py\", line 472, in _handle_recv\n",
      "    self._run_callback(callback, msg)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/zmq/eventloop/zmqstream.py\", line 414, in _run_callback\n",
      "    callback(*args, **kwargs)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/tornado/stack_context.py\", line 277, in null_wrapper\n",
      "    return fn(*args, **kwargs)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 283, in dispatcher\n",
      "    return self.dispatch_shell(stream, msg)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 235, in dispatch_shell\n",
      "    handler(stream, idents, msg)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 399, in execute_request\n",
      "    user_expressions, allow_stdin)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/ipkernel.py\", line 196, in do_execute\n",
      "    res = shell.run_cell(code, store_history=store_history, silent=silent)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/ipykernel/zmqshell.py\", line 533, in run_cell\n",
      "    return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2717, in run_cell\n",
      "    interactivity=interactivity, compiler=compiler, result=result)\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2821, in run_ast_nodes\n",
      "    if self.run_code(code, result):\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2881, in run_code\n",
      "    exec(code_obj, self.user_global_ns, self.user_ns)\n",
      "  File \"<ipython-input-1-d90d8da8b2c1>\", line 14, in <module>\n",
      "    from matplotlib import pyplot as plt\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/matplotlib/pyplot.py\", line 71, in <module>\n",
      "    from matplotlib.backends import pylab_setup\n",
      "  File \"/home/oesteban/.anaconda3/lib/python3.6/site-packages/matplotlib/backends/__init__.py\", line 16, in <module>\n",
      "    line for line in traceback.format_stack()\n",
      "\n",
      "\n",
      "  matplotlib.use('Agg')\n",
      "/home/oesteban/.anaconda3/lib/python3.6/importlib/_bootstrap.py:219: ImportWarning: can't resolve package from __spec__ or __package__, falling back on __name__ and __path__\n",
      "  return f(*args, **kwds)\n"
     ]
    }
   ],
   "source": [
    "# %load_ext autoreload\n",
    "# %autoreload 2\n",
    "# %matplotlib inline\n",
    "import os\n",
    "from pathlib import Path\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import nibabel as nb\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib as mpl\n",
    "mpl.use('pgf')\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import gridspec\n",
    "import seaborn as sn\n",
    "import palettable\n",
    "\n",
    "from niworkflows.data import get_template\n",
    "\n",
    "from nilearn.image import concat_imgs, mean_img\n",
    "from nilearn import plotting\n",
    "\n",
    "warnings.simplefilter('ignore')\n",
    "\n",
    "DATA_HOME = Path(os.getenv('FMRIPREP_DATA_HOME', os.getcwd())).resolve()\n",
    "DS030_HOME = DATA_HOME / 'ds000030' / '1.0.3'\n",
    "DERIVS_HOME = DS030_HOME / 'derivatives'\n",
    "ATLAS_HOME = get_template('MNI152NLin2009cAsym')\n",
    "ANALYSIS_HOME = DERIVS_HOME / 'fmriprep_vs_feat_2.0-oe'\n",
    "\n",
    "fprep_home = DERIVS_HOME / 'fmriprep_1.0.8' / 'fmriprep'\n",
    "feat_home = DERIVS_HOME / 'fslfeat_5.0.10' / 'featbids'\n",
    "\n",
    "out_folder = Path(os.getenv('FMRIPREP_OUTPUTS') or '').resolve()\n",
    "\n",
    "# Load MNI152 nonlinear, asymmetric 2009c atlas\n",
    "atlas = nb.load(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_T1w.nii.gz'))\n",
    "mask1mm = nb.load(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_brainmask.nii.gz')).get_data() > 0\n",
    "mask2mm = nb.load(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-02_brainmask.nii.gz')).get_data() > 0\n",
    "data = atlas.get_data()\n",
    "data[~mask1mm] = data[~mask1mm].max()\n",
    "atlas = nb.Nifti1Image(data, atlas.affine, atlas.header)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# sn.set_style(\"whitegrid\", {\n",
    "#     'ytick.major.size': 5,\n",
    "#     'xtick.major.size': 5,\n",
    "# })\n",
    "# sn.set_context(\"notebook\", font_scale=1.5)\n",
    "\n",
    "# pgf_with_custom_preamble = {\n",
    "#     'ytick.major.size': 0,\n",
    "#     'xtick.major.size': 0,\n",
    "#     'font.size': 30,\n",
    "#     'font.sans-serif': ['HelveticaLTStd-Light'],\n",
    "#     'font.family': 'sans-serif', # use serif/main font for text elements\n",
    "#     'text.usetex': False,    # use inline math for ticks\n",
    "# }\n",
    "# mpl.rcParams.update(pgf_with_custom_preamble)\n",
    "\n",
    "\n",
    "pgf_with_custom_preamble = {\n",
    "    'text.usetex': True,    # use inline math for ticks\n",
    "    'pgf.rcfonts': False,   # don't setup fonts from rc parameters\n",
    "    'pgf.texsystem': 'xelatex',\n",
    "    'verbose.level': 'debug-annoying',\n",
    "    \"pgf.preamble\": [\n",
    "        r\"\"\"\\usepackage{fontspec}\n",
    "\\setsansfont{HelveticaLTStd-Light}[\n",
    "Extension=.otf,\n",
    "BoldFont=HelveticaLTStd-Bold,\n",
    "ItalicFont=HelveticaLTStd-LightObl,\n",
    "BoldItalicFont=HelveticaLTStd-BoldObl,\n",
    "]\n",
    "\\setmainfont{HelveticaLTStd-Light}[\n",
    "Extension=.otf,\n",
    "BoldFont=HelveticaLTStd-Bold,\n",
    "ItalicFont=HelveticaLTStd-LightObl,\n",
    "BoldItalicFont=HelveticaLTStd-BoldObl,\n",
    "]\n",
    "\\setmonofont{Inconsolata-dz}\n",
    "\"\"\",\n",
    "        r'\\renewcommand\\familydefault{\\sfdefault}',\n",
    "#         r'\\setsansfont[Extension=.otf]{Helvetica-LightOblique}',\n",
    "#         r'\\setmainfont[Extension=.ttf]{DejaVuSansCondensed}',\n",
    "#         r'\\setmainfont[Extension=.otf]{FiraSans-Light}',\n",
    "#         r'\\setsansfont[Extension=.otf]{FiraSans-Light}',\n",
    "    ]\n",
    "}\n",
    "mpl.rcParams.update(pgf_with_custom_preamble)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mean_std_map(pipe_home, meanmask, force=False, lazy=False, maskval=1000):\n",
    "    pipe_std = pipe_home / 'summary_stdev.nii.gz'\n",
    "    pipe_mean = pipe_home / 'summary_means.nii.gz'\n",
    "\n",
    "    if force or not pipe_mean.is_file():\n",
    "        print('Forced or %s not found' % pipe_mean)\n",
    "        all_mus = []\n",
    "        if lazy:\n",
    "            all_mus = [nb.load(str(f)) for f in pipe_home.glob(\n",
    "                'sub-*/func/sub-*_task-stopsignal_bold_space-MNI152NLin2009cAsym_avgpreproc.nii.gz')]\n",
    "        \n",
    "        if not all_mus:\n",
    "            print('Generating means file')\n",
    "            pipe_files = list(pipe_home.glob(\n",
    "                'sub-*/func/sub-*_task-stopsignal_bold_space-MNI152NLin2009cAsym_preproc.nii.gz'))\n",
    "            all_mus = []\n",
    "            for f in pipe_files:\n",
    "                mean = mean_img(str(f))\n",
    "                data = mean.get_data()\n",
    "                sigma = np.percentile(data[meanmask], 50) / maskval\n",
    "                data /= sigma\n",
    "                all_mus.append(nb.Nifti1Image(data, mean.affine, mean.header))\n",
    "                \n",
    "        meannii = concat_imgs(all_mus, auto_resample=False)\n",
    "        meannii.to_filename(str(pipe_mean))\n",
    "        force = True\n",
    "\n",
    "    if force or not pipe_std.is_file():\n",
    "        print('Generating standard deviation map')\n",
    "        meannii = nb.load(str(pipe_mean))\n",
    "        nb.Nifti1Image(meannii.get_data().std(3), meannii.affine, meannii.header).to_filename(str(pipe_std))\n",
    "        \n",
    "    return pipe_mean, pipe_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Use the WM mask to normalize intensities of EPI means\n",
    "meanmask = nb.load(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-02_class-WM_probtissue.nii.gz')).get_data() > 0.9\n",
    "\n",
    "# Calculate average and std\n",
    "fprep_mean, fprep_std = mean_std_map(fprep_home, meanmask)\n",
    "feat_mean, feat_std = mean_std_map(feat_home, meanmask)\n",
    "\n",
    "# Trick to avoid nilearn zooming in\n",
    "feat_nii = nb.load(str(feat_std))\n",
    "\n",
    "fd = feat_nii.get_data()\n",
    "fd[0, 0, :] = 50\n",
    "fd[0, -1, :] = 50\n",
    "fd[-1, 0, :] = 50\n",
    "fd[-1, -1, :] = 50\n",
    "\n",
    "nb.Nifti1Image(fd, feat_nii.affine, feat_nii.header).to_filename('newfeat.nii.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(str(ANALYSIS_HOME / 'smoothness.csv'))\n",
    "\n",
    "plt.clf()\n",
    "fig = plt.gcf()\n",
    "_ = fig.set_size_inches(15, 2 * 3.1)\n",
    "# gs = gridspec.GridSpec(2, 4, width_ratios=[38, 7, 60, 10], height_ratios=[1, 1],  hspace=0.0, wspace=0.03)\n",
    "gs = gridspec.GridSpec(2, 3, width_ratios=[42, 9, 64], height_ratios=[1, 1],  hspace=0.0, wspace=0.03)\n",
    "\n",
    "a_ax1 = plt.subplot(gs[0, 0])\n",
    "a_ax2 = plt.subplot(gs[1, 0])\n",
    "\n",
    "\n",
    "fmriprep_smooth = df[df.pipeline.str.contains('fmriprep')][['fwhm_pre', 'fwhm_post']]\n",
    "feat_smooth = df[df.pipeline.str.contains('feat')][['fwhm_pre', 'fwhm_post']]\n",
    "\n",
    "cols = palettable.tableau.ColorBlind_10.hex_colors\n",
    "sn.distplot(fmriprep_smooth.fwhm_post, color=cols[0], ax=a_ax2,\n",
    "            axlabel='Smoothing', label='fMRIPrep')\n",
    "sn.distplot(feat_smooth.fwhm_post, color=cols[1], ax=a_ax2,\n",
    "            axlabel='Smoothing', label=r'\\texttt{feat}')\n",
    "sn.distplot(fmriprep_smooth.fwhm_pre, color=cols[0], ax=a_ax1,\n",
    "            axlabel='Smoothing', label='fMRIPrep')\n",
    "sn.distplot(feat_smooth.fwhm_pre, color=cols[1], ax=a_ax1,\n",
    "            axlabel='Smoothing', label=r'\\texttt{feat}')\n",
    "\n",
    "a_ax2.set_xlim([3, 8.8])\n",
    "a_ax2.set_xticks([])\n",
    "a_ax2.set_xticklabels([])\n",
    "a_ax2.xaxis.tick_bottom()\n",
    "a_ax2.grid(False)\n",
    "a_ax2.set_xlabel('')\n",
    "\n",
    "a_ax2.set_ylim([-1.1, 9.9])\n",
    "a_ax2.set_yticks([])\n",
    "a_ax2.spines['left'].set_visible(False)\n",
    "\n",
    "a_ax1.set_ylabel(r'\\noindent\\parbox{4.8cm}{\\centering\\textbf{Before smoothing} fraction of images}',\n",
    "                 size=13)\n",
    "a_ax1.yaxis.set_label_coords(-0.1, 0.4)\n",
    "a_ax2.set_ylabel(r'\\noindent\\parbox{4.8cm}{\\centering\\textbf{After smoothing} fraction of images}',\n",
    "                 size=13)\n",
    "a_ax2.yaxis.set_label_coords(-0.1, 0.6)\n",
    "\n",
    "# ax4.spines['bottom'].set_position(('outward', 20))\n",
    "a_ax2.invert_yaxis()\n",
    "a_ax2.spines['top'].set_visible(False)\n",
    "a_ax2.spines['bottom'].set_visible(False)\n",
    "a_ax2.spines['left'].set_visible(False)\n",
    "a_ax2.spines['right'].set_visible(False)\n",
    "\n",
    "a_ax1.set_xlim([3, 8.8])\n",
    "a_ax1.set_ylim([-0.6, 10.4])\n",
    "a_ax1.grid(False)\n",
    "\n",
    "a_ax1.set_xlabel('(mm)')\n",
    "a_ax1.xaxis.set_label_coords(0.95, 0.1)\n",
    "a_ax1.set_yticks([])\n",
    "a_ax1.set_yticklabels([])\n",
    "a_ax1.set_xticks([3, 4, 5, 6 , 7, 8])\n",
    "a_ax1.set_xticklabels([3, 4, 5, 6 , 7, 8])\n",
    "a_ax1.tick_params(axis='x', zorder=100, direction='inout')\n",
    "a_ax1.spines['left'].set_visible(False)\n",
    "a_ax1.spines['right'].set_visible(False)\n",
    "a_ax1.spines['top'].set_visible(False)\n",
    "a_ax1.spines['bottom'].set_visible(True)\n",
    "a_ax1.zorder = 100\n",
    "\n",
    "# a_ax2.xaxis.set_label_position('top')\n",
    "# a_ax2.xaxis.set_label_coords(0.45, 0.95)\n",
    "a_ax1.annotate(\n",
    "    r'\\noindent\\parbox{6.8cm}{\\centering\\textbf{Estimated smoothness} full~width~half~maximum~(mm)}',\n",
    "    xy=(0.15, 0.8), xycoords='axes fraction', xytext=(.0, .0),\n",
    "    textcoords='offset points', va='center', color='k', size=13,\n",
    ")\n",
    "\n",
    "legend = a_ax2.legend(ncol=2, loc='upper center', bbox_to_anchor=(0.5, 0.45), prop={'size': 15})\n",
    "legend.get_frame().set_facecolor('w')\n",
    "legend.get_frame().set_edgecolor('none')\n",
    "\n",
    "# a_ax2.annotate(\n",
    "#     r'\\noindent\\parbox{15cm}{Panels A, B present statistics derived from N=257 biologically independent participants}',\n",
    "#     xy=(-0.2, 0.03), xycoords='axes fraction', xytext=(.0, .0),\n",
    "#     textcoords='offset points', va='center', color='k', size=11,\n",
    "# )\n",
    "\n",
    "###### PANEL B\n",
    "\n",
    "b_ax1 = plt.subplot(gs[0, 2])\n",
    "b_ax2 = plt.subplot(gs[1, 2])\n",
    "\n",
    "\n",
    "thres = 20\n",
    "vmin = 50\n",
    "vmax = 200\n",
    "\n",
    "disp = plotting.plot_anat(str(fprep_std), display_mode='z', annotate=False,\n",
    "                          cut_coords=[-5, 10, 20], cmap='cividis', threshold=thres, vmin=vmin, vmax=vmax,\n",
    "                          axes=b_ax1)\n",
    "disp.annotate(size=12, left_right=True, positions=True)\n",
    "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_class-CSF_probtissue.nii.gz'), colors=['k'], levels=[0.8])\n",
    "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_class-WM_probtissue.nii.gz'), colors=['w'], levels=[0.8], linewidths=[1], alpha=0.7)\n",
    "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_brainmask.nii.gz'), colors=['k'], levels=[0.8], linewidths=[3], alpha=.7)\n",
    "\n",
    "\n",
    "disp = plotting.plot_anat('newfeat.nii.gz', display_mode='z', annotate=False,\n",
    "                          cut_coords=[-5, 10, 20], cmap='cividis', threshold=thres, vmin=vmin, vmax=vmax,\n",
    "                          axes=b_ax2)\n",
    "disp.annotate(size=12, left_right=False, positions=False)\n",
    "disp.annotate(size=12, left_right=False, positions=False, scalebar=True,\n",
    "              loc=3, size_vertical=2, label_top=False, frameon=True, borderpad=0.1)\n",
    "\n",
    "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_class-CSF_probtissue.nii.gz'), colors=['k'], levels=[0.8])\n",
    "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_class-WM_probtissue.nii.gz'), colors=['w'], levels=[0.8], linewidths=[1], alpha=0.7)\n",
    "disp.add_contours(str(ATLAS_HOME / 'tpl-MNI152NLin2009cAsym_space-MNI_res-01_brainmask.nii.gz'), colors=['k'], levels=[0.8], linewidths=[3], alpha=.7)\n",
    "\n",
    "b_ax1.annotate(\n",
    "    'fMRIPrep',\n",
    "    xy=(0., .5), xycoords='axes fraction', xytext=(-20, .0),\n",
    "    textcoords='offset points', va='center', color='k', size=15,\n",
    "    rotation=90)\n",
    "\n",
    "b_ax2.annotate(\n",
    "    r'\\texttt{feat}',\n",
    "    xy=(0., .5), xycoords='axes fraction', xytext=(-20, .0),\n",
    "    textcoords='offset points', va='center', color='k', size=12,\n",
    "    rotation=90)\n",
    "\n",
    "\n",
    "# inner_grid = gridspec.GridSpecFromSubplotSpec(1, 2, width_ratios=[1, 15],\n",
    "#                                               subplot_spec=gs[:, -1], wspace=0.01)\n",
    "\n",
    "# b_ax3 = fig.add_subplot(inner_grid[0])\n",
    "# gradient = np.hstack((np.zeros((50,)), np.linspace(0, 1, 120), np.ones((130,))))[::-1]\n",
    "# gradient = np.vstack((gradient, gradient))\n",
    "# b_ax3.imshow(gradient.T, aspect='auto', cmap=plt.get_cmap('cividis'))\n",
    "# b_ax3.xaxis.set_ticklabels([])\n",
    "# b_ax3.xaxis.set_ticks([])\n",
    "# b_ax3.yaxis.set_ticklabels([])\n",
    "# b_ax3.yaxis.set_ticks([])\n",
    "\n",
    "# b_ax4 = fig.add_subplot(inner_grid[1])\n",
    "# sn.distplot(nb.load(str(fprep_std)).get_data()[mask2mm], label='fMRIPrep',\n",
    "#             vertical=True, ax=b_ax4, kde=False, norm_hist=True)\n",
    "# sn.distplot(nb.load(str(feat_std)).get_data()[mask2mm], label=r'\\texttt{feat}', vertical=True,\n",
    "#             color='darkorange', ax=b_ax4, kde=False, norm_hist=True)\n",
    "\n",
    "# # plt.gca().set_ylim((0, 300))\n",
    "# plt.legend(prop={'size': 15}, edgecolor='none')\n",
    "\n",
    "# b_ax4.xaxis.set_ticklabels([])\n",
    "# b_ax4.xaxis.set_ticks([])\n",
    "# b_ax4.yaxis.set_ticklabels([])\n",
    "# b_ax4.yaxis.set_ticks([])\n",
    "\n",
    "# plt.axis('off')\n",
    "# b_ax3.axis('off')\n",
    "# b_ax4.axis('off')\n",
    "\n",
    "\n",
    "a_ax1.set_title('A', fontdict={'fontsize': 24}, loc='left', x=-0.2);\n",
    "b_ax1.set_title('B', fontdict={'fontsize': 24}, loc='left');\n",
    "\n",
    "plt.savefig(str(out_folder / 'figure03.pdf'),\n",
    "            format='pdf', bbox_inches='tight', pad_inches=0.2, dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coords = [-27, 0, 7]\n",
    "\n",
    "thres = 20\n",
    "vmin = 50\n",
    "vmax = 200\n",
    "\n",
    "\n",
    "# Plot\n",
    "plt.clf()\n",
    "fig = plt.figure(figsize=(20,10))\n",
    "\n",
    "plotting.plot_anat('newfeat.nii.gz', cut_coords=coords, colorbar=True, cmap='cividis',\n",
    "                   threshold=thres, vmin=vmin, vmax=vmax, title='feat',\n",
    "                   axes=plt.subplot(2,2,1)\n",
    ");\n",
    "plotting.plot_anat(str(fprep_std), cut_coords=coords, colorbar=True, cmap='cividis',\n",
    "                   threshold=thres, vmin=vmin, vmax=vmax, title='fmriprep',\n",
    "                   axes=plt.subplot(2,2,3)\n",
    ");\n",
    "plotting.plot_glass_brain(str(feat_std), threshold=200, colorbar=True, title='feat',\n",
    "                          axes=plt.subplot(2,2,2));\n",
    "plotting.plot_glass_brain(str(fprep_std), threshold=200, colorbar=True, title='fmriprep',\n",
    "                          axes=plt.subplot(2,2,4));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}