# Reading and processing NIRS data with MNE-Python

This tutorial demonstrates how to read and process NIRS data using
MNE-Python and `almirah`.

## Setup

First, we'll import the necessary libraries and set up the logging and
warning configurations.

```python
import mne
import warnings
import numpy as np
import matplotlib.pyplot as plt
from itertools import compress

from almirah import Layout

mne.set_log_level(False)
warnings.filterwarnings('ignore')
```

## Loading the Data

Next, we'll set up the layout to access the NIRS data.

```python
lay = Layout(root="/path/to/data", specification_name="bids")
print(lay)
```

This should output:

    <Layout root: '/path/to/data'>

We can query the layout to find all NIRS files with the `.snirf` extension:

```python
files = lay.query(datatype="nirs", extension=".snirf")
len(files)
```

This gives the total number of NIRS files:

    1315

## Querying a Specific File

To query a specific file, we can filter by subject, task, datatype, and extension:

```python
file = lay.query(subject="D0019", task="rest", datatype="nirs", extension=".snirf")[0]
print(file.rel_path)
```

This should output the relative path of the file:

    sub-D0019/ses-111/nirs/sub-D0019_ses-111_task-rest_run-01_nirs.snirf

We can then download the NIRS file:

```python
file.download()
```

## Reading Raw Data

Next, we read the raw NIRS data:

```python
raw = mne.io.read_raw_snirf(file.path)
raw.load_data()
```

<details open>
    <summary><strong>General</strong></summary>
    <table class="table table-hover table-striped table-sm table-responsive small">
        <tr>
            <th>Measurement date</th>
            <td>November 12, 1917  00:00:00 GMT</td>
        </tr>
        <tr>
            <th>Experimenter</th>
            <td>Unknown</td>
        </tr>
        <tr>
            <th>Participant</th>
            <td>mne_anonymize</td>
        </tr>
    </table>
    </details>
    <details open>
        <summary><strong>Channels</strong></summary>
        <table class="table table-hover table-striped table-sm table-responsive small">
            <tr>
                <th>Digitized points</th>
                <td>21 points</td>
            </tr>
            <tr>
                <th>Good channels</th>
                <td>36 fNIRS (CW amplitude)</td>
            </tr>
            <tr>
                <th>Bad channels</th>
                <td>None</td>
            </tr>
            <tr>
                <th>EOG channels</th>
                <td>Not available</td>
            </tr>
            <tr>
                <th>ECG channels</th>
                <td>Not available</td>
            </tr>
        </table>
        </details>
        <details open>
            <summary><strong>Data</strong></summary>
            <table class="table table-hover table-striped table-sm table-responsive small">
                <tr>
                    <th>Sampling frequency</th>
                    <td>15.62 Hz</td>
                </tr>
                <tr>
                    <th>Highpass</th>
                    <td>0.00 Hz</td>
                </tr>
                <tr>
                    <th>Lowpass</th>
                    <td>7.81 Hz</td>
                </tr>
                <tr>
                    <th>Filenames</th>
                    <td>sub-D0019_ses-111_task-rest_run-01_nirs.snirf</td>
                </tr>
                <tr>
                    <th>Duration</th>
                    <td>00:10:03 (HH:MM:SS)</td>
                </tr>
            </table>
            </details>

```python
print(raw.info)
```

    <Info | 9 non-empty values
     bads: []
     ch_names: S1_D1 760, S1_D1 850, S1_D2 760, S1_D2 850, S2_D2 760, S2_D2 ...
     chs: 36 fNIRS (CW amplitude)
     custom_ref_applied: False
     dig: 21 items (3 Cardinal, 18 EEG)
     highpass: 0.0 Hz
     lowpass: 7.8 Hz
     meas_date: 1917-11-12 00:00:00 UTC
     nchan: 36
     projs: []
     sfreq: 15.6 Hz
     subject_info: 4 items (dict)
    >

## Preprocessing

We pick the NIRS channels, compute the source-detector distances, and
filter out channels with a distance greater than 0.01:

```python
picks = mne.pick_types(raw.info, meg=False, fnirs=True)
dists = mne.preprocessing.nirs.source_detector_distances(raw.info, picks=picks)
raw.pick(picks[dists > 0.01])
raw.plot(n_channels=len(raw.ch_names), duration=500, show_scrollbars=False)
plt.show()
```
    
![png](../images/nirs/raw-plot.png)
    
We convert the raw data to optical density:

```python
raw_od = mne.preprocessing.nirs.optical_density(raw)
raw_od.plot(n_channels=len(raw_od.ch_names), duration=500, show_scrollbars=False)
plt.show()
```

![png](../images/nirs/optical-density-plot.png)
    
Next, we compute the Scalp Coupling Index (SCI) and plot the distribution:

```python
sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od)
fig, ax = plt.subplots()
ax.hist(sci)
ax.set(xlabel="Scalp Coupling Index", ylabel="Count", xlim=[0, 1])
plt.show()
```

![png](../images/nirs/scalp-coupling-index-plot.png)
    
We mark channels with a low SCI as bad:

```python
raw_od.info["bads"] = list(compress(raw_od.ch_names, sci < 0.5))
```

We convert the optical density data to hemoglobin concentration:

```python
raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od, ppf=0.1)
raw_haemo.plot(n_channels=len(raw_haemo.ch_names), duration=500, show_scrollbars=False)
plt.show()
```
    
![png](../images/nirs/raw-haemo-plot.png)
    
## Power Spectral Analysis

We compute and plot the Power Spectral Density (PSD) before and after filtering:

```python
fig = raw_haemo.compute_psd().plot(average=True, amplitude=False)
fig.suptitle("Before filtering", weight="bold", size="x-large")
raw_haemo = raw_haemo.filter(0.05, 0.7, h_trans_bandwidth=0.2, l_trans_bandwidth=0.02)
fig = raw_haemo.compute_psd().plot(average=True)
fig.suptitle("After filtering", weight="bold", size="x-large")
plt.show()
```

![png](../images/nirs/raw-haemo-psd-before-filtering.png)
        
![png](../images/nirs/raw-haemo-psd-after-filtering.png)
    
## Epoching

We create epochs from the continuous data:

```python
epochs = mne.make_fixed_length_epochs(raw_haemo, duration=30, preload=False)
```

Finally, we plot the epochs:

```python
epochs.plot_image(combine="mean", vmin=-30, vmax=30,
                             ts_args=dict(ylim=dict(hbo=[-15, 15],
                                                    hbr=[-15, 15])))
plt.show()
```
    
![png](../images/nirs/deoxyhemoglobin-plot.png)
        
![png](../images/nirs/oxyhemoglobin-plot.png)

This concludes the tutorial. You've learned how to read, and process
near-infrared spectroscopy data using MNE-Python and `almirah`.