# Author: Robert Luke <mail@robertluke.net>
#
# License: BSD (3-clause)
import os.path as op
import pandas as pd
import numpy as np
import mne
from mne.transforms import apply_trans, _get_trans
from mne.utils import _validate_type, _check_fname, warn
from mne.io import BaseRaw
from mne.io.constants import FIFF
def _read_fold_xls(fname, atlas="Juelich"):
"""Read fOLD toolbox xls file.
The values are then manipulated in to a tidy dataframe.
Note the xls files are not included as no license is provided.
Parameters
----------
fname : str
Path to xls file.
atlas : str
Requested atlas.
"""
page_reference = {"AAL2": 2,
"AICHA": 5,
"Brodmann": 8,
"Juelich": 11,
"Loni": 14}
tbl = pd.read_excel(fname,
sheet_name=page_reference[atlas])
# Remove the spacing between rows
empty_rows = np.where(np.isnan(tbl["Specificity"]))[0]
tbl = tbl.drop(empty_rows).reset_index(drop=True)
# Empty values in the table mean its the same as above
for row_idx in range(1, tbl.shape[0]):
for col_idx, col in enumerate(tbl.columns):
if not isinstance(tbl[col][row_idx], str):
if np.isnan(tbl[col][row_idx]):
tbl.iloc[row_idx, col_idx] = \
tbl.iloc[row_idx - 1, col_idx]
tbl["Specificity"] = tbl["Specificity"] * 100
tbl["brainSens"] = tbl["brainSens"] * 100
return tbl
def _generate_montage_locations():
"""Get standard MNI montage locations in dataframe.
Data is returned in the same format as the eeg_positions library.
"""
# standard_1020 and standard_1005 are in MNI (fsaverage) space already,
# but we need to undo the scaling that head_scale will do
montage = mne.channels.make_standard_montage(
'standard_1005', head_size=0.09700884729534559)
for d in montage.dig:
d['coord_frame'] = FIFF.FIFFV_MNE_COORD_MNI_TAL
montage.dig[:] = montage.dig[3:]
montage.add_mni_fiducials() # now in fsaverage space
coords = pd.DataFrame.from_dict(
montage.get_positions()['ch_pos']).T
coords["label"] = coords.index
coords = coords.rename(columns={0: "x", 1: "y", 2: "z"})
return coords.reset_index(drop=True)
def _find_closest_standard_location(position, reference, *, out='label'):
"""Return closest montage label to coordinates.
Parameters
----------
position : array, shape (3,)
Coordinates.
reference : dataframe
As generated by _generate_montage_locations.
trans_pos : str
Apply a transformation to positions to specified frame.
Use None for no transformation.
"""
from scipy.spatial.distance import cdist
p0 = np.array(position)
p0.shape = (-1, 3)
head_mri_t, _ = _get_trans('fsaverage', 'head', 'mri')
p0 = apply_trans(head_mri_t, p0)
dists = cdist(p0, np.asarray(reference[['x', 'y', 'z']], float))
if out == 'label':
min_idx = np.argmin(dists)
return reference["label"][min_idx]
else:
assert out == 'dists'
return dists
[docs]
def fold_landmark_specificity(raw, landmark, fold_files=None,
atlas="Juelich", interpolate=False):
"""Return the specificity of each channel to a specified brain landmark.
Parameters
----------
raw : BaseRaw
The fNIRS data.
landmark : str
Landmark of interest. Must be present in fOLD toolbox data file.
fold_files : list | path-like | None
If None, will use the MNE_NIRS_FOLD_PATH config variable.
If path-like, should be a path to a directory containing '10-10.xls'
and '10-5.xls'. If list, should be paths to the fold toolbox files.
See the Notes section of :func:`~mne_nirs.io.fold_channel_specificity`
for details.
atlas : str
Brain atlas to use.
interpolate : bool
If the optimal source-detector pair is not found in the fOLD files
False (default) will yield no results for that pairing, whereas
True will use the next closest match. See Notes of
:func:`mne_nirs.io.fold_channel_specificity` for an example.
.. warning::
The sensitivity profile can differ substantially for nearest
neighbors, so use ``interpolate=True`` with caution.
Returns
-------
spec : array
Specificity values for each channel to brain landmark.
See Also
--------
fold_landmark_specificity
Notes
-----
Specificity values are provided by the fOLD toolbox
:footcite:`morais2018fnirs` excel files. See the Notes section of
:func:`~mne_nirs.io.fold_channel_specificity` for more details.
References
----------
.. footbibliography::
"""
_validate_type(landmark, str, 'landmark')
_validate_type(raw, BaseRaw, 'raw')
reference_locations = _generate_montage_locations()
fold_tbl = _check_load_fold(fold_files, atlas)
specificity = np.zeros(len(raw.ch_names))
for cidx in range(len(raw.ch_names)):
tbl = _source_detector_fold_table(
raw, cidx, reference_locations, fold_tbl, interpolate)
if len(tbl) > 0:
tbl["ContainsLmk"] = [landmark in la for la in tbl["Landmark"]]
tbl = tbl.query("ContainsLmk == True")["Specificity"]
if len(tbl) == 0:
continue
# print(f"No data for {src_name}-{det_name}")
elif len(tbl) == 1:
specificity[cidx] = tbl.values[0]
else:
raise RuntimeError("Multiple specificity values returned")
return np.array(specificity)
[docs]
def fold_channel_specificity(raw, fold_files=None, atlas="Juelich",
interpolate=False):
"""Return the landmarks and specificity a channel is sensitive to.
Parameters
----------
raw : BaseRaw
The fNIRS data.
fold_files : list | path-like | None
If None, will use the MNE_NIRS_FOLD_PATH config variable.
If path-like, should be a path to a directory containing '10-10.xls'
and '10-5.xls'. If list, should be paths to the fold toolbox files.
See Notes for details.
atlas : str
Brain atlas to use.
interpolate : bool
If the optimal source-detector pair is not found in the fOLD files
False (default) will yield no results for that pairing, whereas
True will use the next closest match. See Notes for an example.
.. warning:: The sensitivity profile can differ substantially for
nearest neighbors, so use ``interpolate=True`` with
caution.
Returns
-------
spec : list of DataFrame
List of dataframes, one for each channel.
See Also
--------
fold_landmark_specificity
Notes
-----
**fOLD Toolbox**
Specificity values are provided by the fOLD toolbox
:footcite:`morais2018fnirs` excel files.
For licensing reasons, these files are not distributed with MNE-NIRS.
You need to download them from
`the author's website <https://github.com/nirx/fOLD-public>`__.
To automatically utilize the ``MNE_NIRS_FOLD_PATH`` config for the
``fold_files`` parameter, you can download the entire ``fOLD-public``
repository `as a zip <https://github.com/nirx/fOLD-public/archive/refs/heads/master.zip>`__
and expand it to some suitable location like
``~/mne_data/fOLD/fOLD-public-master``, and then set the config value
on your machine by using :func:`mne:mne.set_config` like::
>>> mne.set_config('MNE_NIRS_FOLD_PATH', '~/mne_data/fOLD/fOLD-public-master/Supplementary')
From then on, :func:`~mne_nirs.io.fold_channel_specificity` and
:func:`~mne_nirs.io.fold_landmark_specificity` will automatically use this
directory to find the fOLD xls files when you pass ``fold_files=None``
(which is the default). We recommend following this procedure so that
the files can be reused automatically.
**Interpolation**
For an example of interpolation, consider the pairings P5-PO7 nor P6-PO8,
neither of which are listed in the fOLD toolbox:
- With ``interpolate=False``, the returned ``spec`` will not have entries
for these channels. When
- With ``interpolate=True``, entries like P7-PO7 and P8-PO8, respectively,
might be used instead. A warning is emitted if such substitutions are
made.
References
----------
.. footbibliography::
""" # noqa: E501
_validate_type(raw, BaseRaw, 'raw')
reference_locations = _generate_montage_locations()
fold_tbl = _check_load_fold(fold_files, atlas)
chan_spec = list()
for cidx in range(len(raw.ch_names)):
tbl = _source_detector_fold_table(
raw, cidx, reference_locations, fold_tbl, interpolate)
chan_spec.append(tbl.reset_index(drop=True))
return chan_spec
def _check_load_fold(fold_files, atlas):
_validate_type(fold_files, (list, 'path-like', None), 'fold_files')
if fold_files is None:
fold_files = mne.get_config('MNE_NIRS_FOLD_PATH')
if fold_files is None:
raise ValueError(
'MNE_NIRS_FOLD_PATH not set, either set it using '
'mne.set_config or pass fold_files as str or list')
if not isinstance(fold_files, list): # path-like
fold_files = _check_fname(
fold_files, overwrite='read', must_exist=True, name='fold_files',
need_dir=True)
fold_files = [op.join(fold_files, f'10-{x}.xls') for x in (5, 10)]
fold_tbl = pd.DataFrame()
for fi, fname in enumerate(fold_files):
fname = _check_fname(fname, overwrite='read', must_exist=True,
name=f'fold_files[{fi}]')
fold_tbl = pd.concat([fold_tbl, _read_fold_xls(fname, atlas=atlas)],
ignore_index=True)
return fold_tbl
def _source_detector_fold_table(raw, cidx, reference, fold_tbl, interpolate):
src = raw.info['chs'][cidx]['loc'][3:6]
det = raw.info['chs'][cidx]['loc'][6:9]
ref_lab = list(reference['label'])
dists = _find_closest_standard_location(
[src, det], reference, out='dists')
src_min, det_min = np.argmin(dists, axis=1)
src_name, det_name = ref_lab[src_min], ref_lab[det_min]
tbl = fold_tbl.query("Source == @src_name and Detector == @det_name")
dist = np.linalg.norm(dists[[0, 1], [src_min, det_min]])
# Try reversing source and detector
if len(tbl) == 0:
tbl = fold_tbl.query("Source == @det_name and Detector == @src_name")
if len(tbl) == 0 and interpolate:
# Try something hopefully not too terrible: pick the one with the
# smallest net distance
good = (np.in1d(fold_tbl['Source'], reference['label']) &
np.in1d(fold_tbl['Detector'], reference['label']))
assert good.any()
tbl = fold_tbl[good]
assert len(tbl)
src_idx = [ref_lab.index(src) for src in tbl['Source']]
det_idx = [ref_lab.index(det) for det in tbl['Detector']]
# Original
tot_dist = np.linalg.norm(
[dists[0, src_idx], dists[1, det_idx]], axis=0)
assert tot_dist.shape == (len(tbl),)
idx = np.argmin(tot_dist)
dist_1 = tot_dist[idx]
src_1, det_1 = ref_lab[src_idx[idx]], ref_lab[det_idx[idx]]
# And the reverse
tot_dist = np.linalg.norm(
[dists[0, det_idx], dists[1, src_idx]], axis=0
)
idx = np.argmin(tot_dist)
dist_2 = tot_dist[idx]
src_2, det_2 = ref_lab[det_idx[idx]], ref_lab[src_idx[idx]]
if dist_1 < dist_2:
new_dist, src_use, det_use = dist_1, src_1, det_1
else:
new_dist, src_use, det_use = dist_2, det_2, src_2
warn('No fOLD table entry for best matching source/detector pair '
f'{src_name}/{det_name} (RMS distance to the channel positions '
f'was {1000 * dist:0.1f} mm) for channel index {cidx}, '
'using next smallest available '
f'src/det pairing {src_use}/{det_use} (RMS distance '
f'{1000 * new_dist:0.1f} mm). Consider setting your channel '
'positions to standard 10-05 locations using raw.set_montage '
'if your pair does show up in the tables.',
module='mne_nirs', ignore_namespaces=('mne', 'mne_nirs'))
tbl = fold_tbl.query("Source == @src_use and Detector == @det_use")
tbl = tbl.copy()
tbl['BestSource'] = src_name
tbl['BestDetector'] = det_name
tbl['BestMatchDistance'] = dist
tbl['MatchDistance'] = new_dist
assert len(tbl)
else:
tbl = tbl.copy()
tbl['BestSource'] = src_name
tbl['BestDetector'] = det_name
tbl['BestMatchDistance'] = dist
tbl['MatchDistance'] = dist
tbl = tbl.copy() # don't get warnings about setting values later
return tbl