{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Preprocessing functional near-infrared spectroscopy (fNIRS) data\n\nThis tutorial covers how to convert functional near-infrared spectroscopy (fNIRS) data\nfrom raw measurements to relative oxyhaemoglobin (HbO) and deoxyhaemoglobin (HbR)\nconcentration, view the average waveform, and topographic representation of the\nresponse.\n\nHere we will work with the `fNIRS motor data <fnirs-motor-dataset>`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: The MNE-Python contributors.\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from itertools import compress\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport mne\n\nfnirs_data_folder = mne.datasets.fnirs_motor.data_path()\nfnirs_cw_amplitude_dir = fnirs_data_folder / \"Participant-1\"\nraw_intensity = mne.io.read_raw_nirx(fnirs_cw_amplitude_dir, verbose=True)\nraw_intensity.load_data()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Providing more meaningful annotation information\n\nFirst, we attribute more meaningful names to the trigger codes which are\nstored as annotations. Second, we include information about the duration of\neach stimulus, which was 5 seconds for all conditions in this experiment.\nThird, we remove the trigger code 15, which signaled the start and end\nof the experiment and is not relevant to our analysis.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw_intensity.annotations.set_durations(5)\nraw_intensity.annotations.rename(\n    {\"1.0\": \"Control\", \"2.0\": \"Tapping/Left\", \"3.0\": \"Tapping/Right\"}\n)\nunwanted = np.nonzero(raw_intensity.annotations.description == \"15.0\")\nraw_intensity.annotations.delete(unwanted)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Viewing location of sensors over brain surface\n\nHere we validate that the location of sources-detector pairs and channels\nare in the expected locations. Source-detector pairs are shown as lines\nbetween the optodes, channels (the mid point of source-detector pairs) are\noptionally shown as orange dots. Source are optionally shown as red dots and\ndetectors as black.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "subjects_dir = mne.datasets.sample.data_path() / \"subjects\"\n\nbrain = mne.viz.Brain(\n    \"fsaverage\", subjects_dir=subjects_dir, background=\"w\", cortex=\"0.5\"\n)\nbrain.add_sensors(\n    raw_intensity.info,\n    trans=\"fsaverage\",\n    fnirs=[\"channels\", \"pairs\", \"sources\", \"detectors\"],\n)\nbrain.show_view(azimuth=20, elevation=60, distance=400)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Selecting channels appropriate for detecting neural responses\n\nFirst we remove channels that are too close together (short channels) to\ndetect a neural response (less than 1 cm distance between optodes).\nThese short channels can be seen in the figure above.\nTo achieve this we pick all the channels that are not considered to be short.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "picks = mne.pick_types(raw_intensity.info, meg=False, fnirs=True)\ndists = mne.preprocessing.nirs.source_detector_distances(\n    raw_intensity.info, picks=picks\n)\nraw_intensity.pick(picks[dists > 0.01])\nraw_intensity.plot(\n    n_channels=len(raw_intensity.ch_names), duration=500, show_scrollbars=False\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Converting from raw intensity to optical density\n\nThe raw intensity values are then converted to optical density.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)\nraw_od.plot(n_channels=len(raw_od.ch_names), duration=500, show_scrollbars=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Evaluating the quality of the data\n\nAt this stage we can quantify the quality of the coupling\nbetween the scalp and the optodes using the scalp coupling index. This\nmethod looks for the presence of a prominent synchronous signal in the\nfrequency range of cardiac signals across both photodetected signals.\n\nIn this example the data is clean and the coupling is good for all\nchannels, so we will not mark any channels as bad based on the scalp\ncoupling index.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od)\nfig, ax = plt.subplots(layout=\"constrained\")\nax.hist(sci)\nax.set(xlabel=\"Scalp Coupling Index\", ylabel=\"Count\", xlim=[0, 1])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this example we will mark all channels with a SCI less than 0.5 as bad\n(this dataset is quite clean, so no channels are marked as bad).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw_od.info[\"bads\"] = list(compress(raw_od.ch_names, sci < 0.5))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "At this stage it is appropriate to inspect your data\n(for instructions on how to use the interactive data visualisation tool\nsee `tut-visualize-raw`)\nto ensure that channels with poor scalp coupling have been removed.\nIf your data contains lots of artifacts you may decide to apply\nartifact reduction techniques as described in `ex-fnirs-artifacts`.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Converting from optical density to haemoglobin\n\nNext we convert the optical density data to haemoglobin concentration using\nthe modified Beer-Lambert law.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od, ppf=0.1)\nraw_haemo.plot(n_channels=len(raw_haemo.ch_names), duration=500, show_scrollbars=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Removing heart rate from signal\n\nThe haemodynamic response has frequency content predominantly below 0.5 Hz.\nAn increase in activity around 1 Hz can be seen in the data that is due to\nthe person's heart beat and is unwanted. So we use a low pass filter to\nremove this. A high pass filter is also included to remove slow drifts\nin the data.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw_haemo_unfiltered = raw_haemo.copy()\nraw_haemo.filter(0.05, 0.7, h_trans_bandwidth=0.2, l_trans_bandwidth=0.02)\nfor when, _raw in dict(Before=raw_haemo_unfiltered, After=raw_haemo).items():\n    fig = _raw.compute_psd().plot(\n        average=True, amplitude=False, picks=\"data\", exclude=\"bads\"\n    )\n    fig.suptitle(f\"{when} filtering\", weight=\"bold\", size=\"x-large\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Extract epochs\n\nNow that the signal has been converted to relative haemoglobin concentration,\nand the unwanted heart rate component has been removed, we can extract epochs\nrelated to each of the experimental conditions.\n\nFirst we extract the events of interest and visualize them to ensure they are\ncorrect.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "events, event_dict = mne.events_from_annotations(raw_haemo)\nfig = mne.viz.plot_events(events, event_id=event_dict, sfreq=raw_haemo.info[\"sfreq\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next we define the range of our epochs, the rejection criteria,\nbaseline correction, and extract the epochs. We visualize the log of which\nepochs were dropped.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reject_criteria = dict(hbo=80e-6)\ntmin, tmax = -5, 15\n\nepochs = mne.Epochs(\n    raw_haemo,\n    events,\n    event_id=event_dict,\n    tmin=tmin,\n    tmax=tmax,\n    reject=reject_criteria,\n    reject_by_annotation=True,\n    proj=True,\n    baseline=(None, 0),\n    preload=True,\n    detrend=None,\n    verbose=True,\n)\nepochs.plot_drop_log()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## View consistency of responses across trials\n\nNow we can view the haemodynamic response for our tapping condition.\nWe visualize the response for both the oxy- and deoxyhaemoglobin, and\nobserve the expected peak in HbO at around 6 seconds consistently across\ntrials, and the consistent dip in HbR that is slightly delayed relative to\nthe HbO peak.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epochs[\"Tapping\"].plot_image(\n    combine=\"mean\",\n    vmin=-30,\n    vmax=30,\n    ts_args=dict(ylim=dict(hbo=[-15, 15], hbr=[-15, 15])),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also view the epoched data for the control condition and observe\nthat it does not show the expected morphology.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epochs[\"Control\"].plot_image(\n    combine=\"mean\",\n    vmin=-30,\n    vmax=30,\n    ts_args=dict(ylim=dict(hbo=[-15, 15], hbr=[-15, 15])),\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## View consistency of responses across channels\n\nSimilarly we can view how consistent the response is across the optode\npairs that we selected. All the channels in this data are located over the\nmotor cortex, and all channels show a similar pattern in the data.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(15, 6), layout=\"constrained\")\nclims = dict(hbo=[-20, 20], hbr=[-20, 20])\nepochs[\"Control\"].average().plot_image(axes=axes[:, 0], clim=clims)\nepochs[\"Tapping\"].average().plot_image(axes=axes[:, 1], clim=clims)\nfor column, condition in enumerate([\"Control\", \"Tapping\"]):\n    for ax in axes[:, column]:\n        ax.set_title(f\"{condition}: {ax.get_title()}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot standard fNIRS response image\n\nNext we generate the most common visualisation of fNIRS data: plotting\nboth the HbO and HbR on the same figure to illustrate the relation between\nthe two signals.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "evoked_dict = {\n    \"Tapping/HbO\": epochs[\"Tapping\"].average(picks=\"hbo\"),\n    \"Tapping/HbR\": epochs[\"Tapping\"].average(picks=\"hbr\"),\n    \"Control/HbO\": epochs[\"Control\"].average(picks=\"hbo\"),\n    \"Control/HbR\": epochs[\"Control\"].average(picks=\"hbr\"),\n}\n\n# Rename channels until the encoding of frequency in ch_name is fixed\nfor condition in evoked_dict:\n    evoked_dict[condition].rename_channels(lambda x: x[:-4])\n\ncolor_dict = dict(HbO=\"#AA3377\", HbR=\"b\")\nstyles_dict = dict(Control=dict(linestyle=\"dashed\"))\n\nmne.viz.plot_compare_evokeds(\n    evoked_dict, combine=\"mean\", ci=0.95, colors=color_dict, styles=styles_dict\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## View topographic representation of activity\n\nNext we view how the topographic activity changes throughout the response.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "times = np.arange(-3.5, 13.2, 3.0)\ntopomap_args = dict(extrapolate=\"local\")\nepochs[\"Tapping\"].average(picks=\"hbo\").plot_joint(\n    times=times, topomap_args=topomap_args\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compare tapping of left and right hands\n\nFinally we generate topo maps for the left and right conditions to view\nthe location of activity. First we visualize the HbO activity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "times = np.arange(4.0, 11.0, 1.0)\nepochs[\"Tapping/Left\"].average(picks=\"hbo\").plot_topomap(times=times, **topomap_args)\nepochs[\"Tapping/Right\"].average(picks=\"hbo\").plot_topomap(times=times, **topomap_args)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And we also view the HbR activity for the two conditions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epochs[\"Tapping/Left\"].average(picks=\"hbr\").plot_topomap(times=times, **topomap_args)\nepochs[\"Tapping/Right\"].average(picks=\"hbr\").plot_topomap(times=times, **topomap_args)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And we can plot the comparison at a single time point for two conditions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(\n    nrows=2,\n    ncols=4,\n    figsize=(9, 5),\n    gridspec_kw=dict(width_ratios=[1, 1, 1, 0.1]),\n    layout=\"constrained\",\n)\nvlim = (-8, 8)\nts = 9.0\n\nevoked_left = epochs[\"Tapping/Left\"].average()\nevoked_right = epochs[\"Tapping/Right\"].average()\n\nevoked_left.plot_topomap(\n    ch_type=\"hbo\", times=ts, axes=axes[0, 0], vlim=vlim, colorbar=False, **topomap_args\n)\nevoked_left.plot_topomap(\n    ch_type=\"hbr\", times=ts, axes=axes[1, 0], vlim=vlim, colorbar=False, **topomap_args\n)\nevoked_right.plot_topomap(\n    ch_type=\"hbo\", times=ts, axes=axes[0, 1], vlim=vlim, colorbar=False, **topomap_args\n)\nevoked_right.plot_topomap(\n    ch_type=\"hbr\", times=ts, axes=axes[1, 1], vlim=vlim, colorbar=False, **topomap_args\n)\n\nevoked_diff = mne.combine_evoked([evoked_left, evoked_right], weights=[1, -1])\n\nevoked_diff.plot_topomap(\n    ch_type=\"hbo\", times=ts, axes=axes[0, 2:], vlim=vlim, colorbar=True, **topomap_args\n)\nevoked_diff.plot_topomap(\n    ch_type=\"hbr\", times=ts, axes=axes[1, 2:], vlim=vlim, colorbar=True, **topomap_args\n)\n\nfor column, condition in enumerate([\"Tapping Left\", \"Tapping Right\", \"Left-Right\"]):\n    for row, chroma in enumerate([\"HbO\", \"HbR\"]):\n        axes[row, column].set_title(f\"{chroma}: {condition}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lastly, we can also look at the individual waveforms to see what is\ndriving the topographic plot above.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 4), layout=\"constrained\")\nmne.viz.plot_evoked_topo(\n    epochs[\"Left\"].average(picks=\"hbo\"), color=\"b\", axes=axes, legend=False\n)\nmne.viz.plot_evoked_topo(\n    epochs[\"Right\"].average(picks=\"hbo\"), color=\"r\", axes=axes, legend=False\n)\n\n# Tidy the legend:\nleg_lines = [line for line in axes.lines if line.get_c() == \"b\"][:1]\nleg_lines.append([line for line in axes.lines if line.get_c() == \"r\"][0])\nfig.legend(leg_lines, [\"Left\", \"Right\"], loc=\"lower right\")"
      ]
    }
  ],
  "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.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}