Skip to content

Broader support of the SNIRF file format and enable reading GowerLab data #10555

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
agramfort merged 16 commits into from
May 5, 2022
Merged

Broader support of the SNIRF file format and enable reading GowerLab data #10555

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
6f57595
Add support for SNIRF timeunit
rob-luke Apr 21, 2022
08fc60f
Flake
rob-luke Apr 21, 2022
f4fb273
Fix landmark reading
rob-luke May 1, 2022
6b89d3e
Remove commented code
rob-luke May 1, 2022
3b78561
Merge branch 'main' into snirftimeunit
rob-luke May 1, 2022
ffa38db
More tests
rob-luke May 1, 2022
7e55e66
Merge branch 'snirftimeunit' of github.com:rob-luke/mne-python into s…
rob-luke May 1, 2022
7edec29
Report api change
rob-luke May 1, 2022
9fb918c
Update docs
rob-luke May 1, 2022
d28ecb7
Flake
rob-luke May 1, 2022
b1486d6
Update config.py
rob-luke May 1, 2022
5c0b057
Update config.py
rob-luke May 1, 2022
07c4c9b
Update test indicies
rob-luke May 4, 2022
df66ca5
Merge branch 'snirftimeunit' of github.com:rob-luke/mne-python into s…
rob-luke May 4, 2022
a8d6b91
More tests
rob-luke May 4, 2022
59a5c27
Update latest.inc
rob-luke May 5, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions doc/changes/latest.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,6 +51,8 @@ Enhancements

- Add ``'voronoi'`` as an option for the ``image_interp`` argument in :func:`mne.viz.plot_topomap` to plot a topomap without interpolation using a Voronoi parcelation (:gh:`10571` by `Alex Rockhill`_)

- Add support for reading data from Gowerlabs devices to :func:`mne.io.read_raw_snirf` (:gh:`10555` by `Robert Luke`_ and `Samuel Powell`_)

@rob-luke rob-luke May 1, 2022

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@samuelpowell are you ok with acknowledging your assistance here?

@samuelpowell samuelpowell May 1, 2022

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes that’s fine, thank you.

@larsoner larsoner May 3, 2022

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rob-luke this should be :newcontrib: and moved to the top of the new list


Bugs
~~~~
- Make ``color`` parameter check in in :func:`mne.viz.plot_evoked_topo` consistent (:gh:`10217` by :newcontrib:`T. Wang` and `Stefan Appelhoff`_)
Expand Down Expand Up @@ -82,3 +84,5 @@ Bugs
API and behavior changes
~~~~~~~~~~~~~~~~~~~~~~~~
- When creating BEM surfaces via :func:`mne.bem.make_watershed_bem` and :func:`mne.bem.make_flash_bem`, the ``copy`` parameter now defaults to ``True``. This means that instead of creating symbolic links inside the FreeSurfer subject's ``bem`` folder, we now create "actual" files. This should avoid troubles when sharing files across different operating systems and file systems (:gh:`10531` by `Richard Höchenberger`_)

- The ordering of channels returned by :func:`mne.io.read_raw_nirx` is now ordered by channel name, rather than the order provided by the manufacturer. This enables consistent ordering of channels across different file types(:gh:`10555` by `Robert Luke`_)
2 changes: 2 additions & 0 deletions doc/changes/names.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -443,3 +443,5 @@
.. _Matthias Dold: https://matthiasdold.de

.. _T. Wang: https://github.com/twang5

.. _Samuel Powell: https://github.com/samuelpowell
3 changes: 3 additions & 0 deletions mne/io/nirx/nirx.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -459,6 +459,9 @@ def __init__(self, fname, saturated, preload=False, verbose=None):
annot = Annotations(onset, duration, description, ch_names=ch_names)
self.set_annotations(annot)

sort_idx = np.argsort(self.ch_names)
self.pick(picks=sort_idx)

@agramfort agramfort May 5, 2022

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hum. This will make a copy of the data buffer but more importantly this is not consistent with the other readers which read the channels in the order they appear in the file on disk. I prefer to avoid doing magic in readers. Relevant functions down the road should deal with files with channels present in different orders, eg combine_evoked grand_average etc.

@larsoner larsoner May 5, 2022

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will make a copy of the data buffer

If preload=False it won't because we can pick without loaded data now. This means people who really need to save this 2x memory can do read_raw_nirx(...).load_data() instead. But really NIRS data tends to be pretty small so I don't think this will be an issue in practice.

more importantly this is not consistent with the other readers which read the channels in the order they appear in the file on disk.

I agree with this idea in principle, but it will require a big refactoring of our existing code, which assumes and requires that the frequency pairings are essentially properly "interlaced" to work correctly.

So to me I would like to add a comment # TODO: Remove this reordering, which would give us some months to properly craft order-agnostic code, which I think will be a bigger undertaking involving coordination between MNE-Python and MNE-NIRS and changes in probably 5-10 public functions...

@larsoner larsoner May 13, 2022 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rob-luke I realized this is problematic. Let's say we have this arrangement with the more-or-less standard "logical channel numbering":

S1 -- 1 -- D1 -- 2 -- S2
|          |           |
3          4           5
|          |           |
D3 -- 6 -- S3 -- 7 -- D4

The native file will be written as ['S1_D1 760', 'S1_D2_850', 'S2_D1_760', 'S2_D1_850', 'S1_D3_760', 'S1_D3_850', ...], i.e., a pair ordering of [S1_D1, S2_D1, S1_D3, S3_D1, S2_D4, S3_D3, S3_D4], and people will naturally expect these pairs to be associated with channel numbers 1 through 7. But this pick(argsort(ch_names)) will result in a reordered pair ordering of [S1_D1, S1_D3, S2_D1, S2_D4, S3_D1, S3_D3, S3_D4], so channels 1-7 from before now show up in the order [1, 3, 2, 5, 4, 6, 7].

Rather than this pick(argsort(...)) it seems like it would make sense to reorder the channels such that the original "logical channel" ordering is preserved, but such that pairs are still immediately adjacent. We could start by making sure that the first N channels are all fNIRS channels, and that they are paired. Non-fNIRS channels go after that in channels N+1:. Sound good?

Eventually like we said before we could/should relax this by making the channel logic smarter everywhere, but that's a bigger undertaking. All we need to do to fix this immediate/lpressing "logical channel" reordering problem I think is .pick(a_better_order) where a_better_order comes from a few lines of pairing logic.

@rob-luke rob-luke May 13, 2022

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I disagree that there is a standard logical channel ordering. Even for the trivial example which you illustrate I would assign the channel numbers differently to you.

You have assigned

  • Channel 1: s1 d1
  • Channel 2: s2 d1
  • Channel 3: s1 d3

I would have thought the logical ordering was

  • Channel 1: s1 d1
  • Channel 2: s1 d3
  • Channel 3: s2 d1

I also don't think that the channel number has any meaning. And it's the pair naming that should be used.

That's why I went for a simple sort. It will always be the same and can be easily described.

@larsoner larsoner May 14, 2022

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure it can differ by system or setup, but to me the acquisition system sets the "standard order" when it writes data to disk in some order. So the idea is just to respect the order as written to disk as much as possible until we can remove any ordering restrictions whatsoever.

FWIW we used to have this order problem with MEEG data, too, and removing it made the code cleaner in the end. So I'm looking forward to removing the paired-order restriction for fNIRS, as the cleanest long term solution is just to keep the order each manufacturer uses in the first place, or a user does with reorder_channels, etc. Currently reorder_chaanels is problematic for fNIRS and it shouldn't be!


def _read_segment_file(self, data, idx, fi, start, stop, cals, mult):
"""Read a segment of data from a file.

Expand Down
85 changes: 59 additions & 26 deletions mne/io/snirf/_snirf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,10 @@ def __init__(self, fname, optode_frame="unknown",
"MNE does not support this feature. "
"Only the first dataset will be processed.")

manafacturer = _get_metadata_str(dat, "ManufacturerName")
if (optode_frame == "unknown") & (manafacturer == "Gowerlabs"):
optode_frame = "head"

snirf_data_type = np.array(dat.get('nirs/data1/measurementList1'
'/dataType')).item()
if snirf_data_type not in [1, 99999]:
Expand All @@ -109,20 +113,8 @@ def __init__(self, fname, optode_frame="unknown",

last_samps = dat.get('/nirs/data1/dataTimeSeries').shape[0] - 1

samplingrate_raw = np.array(dat.get('nirs/data1/time'))
sampling_rate = 0
if samplingrate_raw.shape == (2, 1):
# specified as onset/samplerate
warn("Onset/sample rate SNIRF not yet supported.")
else:
# specified as time points
fs_diff = np.around(np.diff(samplingrate_raw), decimals=4)
if len(np.unique(fs_diff)) == 1:
# Uniformly sampled data
sampling_rate = 1. / np.unique(fs_diff)
else:
# print(np.unique(fs_diff))
warn("Non uniform sampled data not supported.")
sampling_rate = _extract_sampling_rate(dat)

if sampling_rate == 0:
warn("Unable to extract sample rate from SNIRF file.")

Expand Down Expand Up @@ -280,8 +272,7 @@ def natural_keys(text):
# Update info
info.update(subject_info=subject_info)

LengthUnit = np.array(dat.get('/nirs/metaDataTags/LengthUnit'))
LengthUnit = _correct_shape(LengthUnit)[0].decode('UTF-8')
LengthUnit = _get_metadata_str(dat, "LengthUnit")
scal = 1
if "cm" in LengthUnit:
scal = 100
Expand Down Expand Up @@ -331,15 +322,17 @@ def natural_keys(text):

if 'landmarkPos3D' in dat.get('nirs/probe/'):
diglocs = np.array(dat.get('/nirs/probe/landmarkPos3D'))
diglocs /= scal
digname = np.array(dat.get('/nirs/probe/landmarkLabels'))
nasion, lpa, rpa, hpi = None, None, None, None
extra_ps = dict()
for idx, dign in enumerate(digname):
if dign == b'LPA':
dign = dign.lower()
if dign in [b'lpa', b'al']:
lpa = diglocs[idx, :3]
elif dign == b'NASION':
elif dign in [b'nasion']:
nasion = diglocs[idx, :3]
elif dign == b'RPA':
elif dign in [b'rpa', b'ar']:
rpa = diglocs[idx, :3]
else:
extra_ps[f'EEG{len(extra_ps) + 1:03d}'] = \
Expand Down Expand Up @@ -420,14 +413,10 @@ def natural_keys(text):
annot.append(data[:, 0], 1.0, desc.decode('UTF-8'))
self.set_annotations(annot, emit_warning=False)

# Reorder channels to match expected ordering in MNE if required'
# MNE requires channels are paired as alternating wavelengths
if len(_validate_nirs_info(self.info, throw_errors=False)) == 0:
num_chans = len(self.ch_names)
chans = []
for idx in range(num_chans // 2):
chans.append(idx)
chans.append(idx + num_chans // 2)
self.pick(picks=chans)
sort_idx = np.argsort(self.ch_names)
self.pick(picks=sort_idx)

# Validate that the fNIRS info is correctly formatted
_validate_nirs_info(self.info)
Expand All @@ -447,3 +436,47 @@ def _correct_shape(arr):
if arr.shape == ():
arr = arr[np.newaxis]
return arr


def _get_timeunit_scaling(time_unit):
"""MNE expects time in seconds, return required scaling."""
scalings = {'ms': 1000, 's': 1, 'unknown': 1}
if time_unit in scalings:
return scalings[time_unit]
else:
raise RuntimeError(f'The time unit {time_unit} is not supported by '
'MNE. Please report this error as a GitHub'
'issue to inform the developers.')


def _extract_sampling_rate(dat):
"""Extract the sample rate from the time field."""
time_data = np.array(dat.get('nirs/data1/time'))
sampling_rate = 0
if len(time_data) == 2:
# specified as onset, samplerate
sampling_rate = 1. / (time_data[1] - time_data[0])
else:
# specified as time points
fs_diff = np.around(np.diff(time_data), decimals=4)
if len(np.unique(fs_diff)) == 1:
# Uniformly sampled data
sampling_rate = 1. / np.unique(fs_diff)
else:
warn("MNE does not currently support reading "
"SNIRF files with non-uniform sampled data.")

time_unit = _get_metadata_str(dat, "TimeUnit")
time_unit_scaling = _get_timeunit_scaling(time_unit)
sampling_rate *= time_unit_scaling

return sampling_rate


def _get_metadata_str(dat, field):
if field not in np.array(dat.get('nirs/metaDataTags')):
return None
data = dat.get(f'/nirs/metaDataTags/{field}')
data = _correct_shape(np.array(data))
data = str(data[0], 'utf-8')
return data
66 changes: 42 additions & 24 deletions mne/io/snirf/tests/test_snirf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,16 +50,25 @@
ft_od = op.join(testing_path, 'SNIRF', 'FieldTrip',
'220307_opticaldensity.snirf')

# GowerLabs
lumo110 = op.join(testing_path, 'SNIRF', 'GowerLabs', 'lumomat-1-1-0.snirf')


def _get_loc(raw, ch_name):
return raw.copy().pick(ch_name).info['chs'][0]['loc']


@requires_testing_data
@pytest.mark.filterwarnings('ignore:.*contains 2D location.*:')
@pytest.mark.filterwarnings('ignore:.*measurement date.*:')
@pytest.mark.parametrize('fname', ([sfnirs_homer_103_wShort,
nirx_nirsport2_103,
sfnirs_homer_103_153,
nirx_nirsport2_103,
nirx_nirsport2_103_2,
nirx_nirsport2_103_2,
kernel_hb
kernel_hb,
lumo110
]))
def test_basic_reading_and_min_process(fname):
"""Test reading SNIRF files and minimum typical processing."""
Expand All @@ -73,6 +82,18 @@ def test_basic_reading_and_min_process(fname):
assert 'hbr' in raw


@requires_testing_data
@pytest.mark.filterwarnings('ignore:.*measurement date.*:')
def test_snirf_gowerlabs():
"""Test reading SNIRF files."""
raw = read_raw_snirf(lumo110, preload=True)

assert raw._data.shape == (216, 274)
assert raw.info['dig'][0]['coord_frame'] == FIFF.FIFFV_COORD_HEAD
assert len(raw.ch_names) == 216
assert_allclose(raw.info['sfreq'], 10.0)


@requires_testing_data
def test_snirf_basic():
"""Test reading SNIRF files."""
Expand All @@ -85,33 +106,33 @@ def test_snirf_basic():
# Test channel naming
assert raw.info['ch_names'][:4] == ["S1_D1 760", "S1_D1 850",
"S1_D9 760", "S1_D9 850"]
assert raw.info['ch_names'][24:26] == ["S5_D13 760", "S5_D13 850"]
# assert raw.info['ch_names'][24:26] == ["S5_D8 760", "S5_D8 850"]

@larsoner larsoner May 3, 2022

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cruft, remove?

@rob-luke rob-luke May 4, 2022

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gross 😞 , thanks for finding this mess.


# Test frequency encoding
assert raw.info['chs'][0]['loc'][9] == 760
assert raw.info['chs'][1]['loc'][9] == 850

# Test source locations
assert_allclose([-8.6765 * 1e-2, 0.0049 * 1e-2, -2.6167 * 1e-2],
raw.info['chs'][0]['loc'][3:6], rtol=0.02)
_get_loc(raw, 'S1_D1 760')[3:6], rtol=0.02)
assert_allclose([7.9579 * 1e-2, -2.7571 * 1e-2, -2.2631 * 1e-2],
raw.info['chs'][4]['loc'][3:6], rtol=0.02)
_get_loc(raw, 'S2_D3 760')[3:6], rtol=0.02)
assert_allclose([-2.1387 * 1e-2, -8.8874 * 1e-2, 3.8393 * 1e-2],
raw.info['chs'][8]['loc'][3:6], rtol=0.02)
_get_loc(raw, 'S3_D2 760')[3:6], rtol=0.02)
assert_allclose([1.8602 * 1e-2, 9.7164 * 1e-2, 1.7539 * 1e-2],
raw.info['chs'][12]['loc'][3:6], rtol=0.02)
_get_loc(raw, 'S4_D4 760')[3:6], rtol=0.02)
assert_allclose([-0.1108 * 1e-2, 0.7066 * 1e-2, 8.9883 * 1e-2],
raw.info['chs'][16]['loc'][3:6], rtol=0.02)
_get_loc(raw, 'S5_D5 760')[3:6], rtol=0.02)

# Test detector locations
assert_allclose([-8.0409 * 1e-2, -2.9677 * 1e-2, -2.5415 * 1e-2],
raw.info['chs'][0]['loc'][6:9], rtol=0.02)
_get_loc(raw, 'S1_D1 760')[6:9], rtol=0.02)
assert_allclose([-8.7329 * 1e-2, 0.7577 * 1e-2, -2.7980 * 1e-2],
raw.info['chs'][3]['loc'][6:9], rtol=0.02)
_get_loc(raw, 'S1_D9 850')[6:9], rtol=0.02)
assert_allclose([9.2027 * 1e-2, 0.0161 * 1e-2, -2.8909 * 1e-2],
raw.info['chs'][5]['loc'][6:9], rtol=0.02)
_get_loc(raw, 'S2_D3 850')[6:9], rtol=0.02)
assert_allclose([7.7548 * 1e-2, -3.5901 * 1e-2, -2.3179 * 1e-2],
raw.info['chs'][7]['loc'][6:9], rtol=0.02)
_get_loc(raw, 'S2_D10 850')[6:9], rtol=0.02)

assert 'fnirs_cw_amplitude' in raw

Expand Down Expand Up @@ -186,9 +207,9 @@ def test_snirf_nirsport2():
assert_almost_equal(raw.info['sfreq'], 7.6, decimal=1)

# Test channel naming
assert raw.info['ch_names'][:4] == ['S1_D1 760', 'S1_D1 850',
'S1_D3 760', 'S1_D3 850']
assert raw.info['ch_names'][24:26] == ['S6_D4 760', 'S6_D4 850']
assert raw.info['ch_names'][:4] == ['S10_D3 760', 'S10_D3 850',
'S10_D9 760', 'S10_D9 850']
assert raw.info['ch_names'][24:26] == ['S15_D11 760', 'S15_D11 850']

# Test frequency encoding
assert raw.info['chs'][0]['loc'][9] == 760
Expand Down Expand Up @@ -226,7 +247,7 @@ def test_snirf_nirsport2_w_positions():
# Test channel naming
assert raw.info['ch_names'][:4] == ['S1_D1 760', 'S1_D1 850',
'S1_D6 760', 'S1_D6 850']
assert raw.info['ch_names'][24:26] == ['S6_D4 760', 'S6_D4 850']
assert raw.info['ch_names'][24:26] == ['S6_D14 760', 'S6_D14 850']

# Test frequency encoding
assert raw.info['chs'][0]['loc'][9] == 760
Expand All @@ -238,11 +259,12 @@ def test_snirf_nirsport2_w_positions():
# nirsite https://github.com/mne-tools/mne-testing-data/pull/86
# figure 3
allowed_distance_error = 0.005
distances = source_detector_distances(raw.info)
assert_allclose(distances[::2][:14],

@rob-luke rob-luke May 1, 2022

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As the ordering of the channels is now changed, I test specific source detector pairs, rather than an indexed channel.

[0.0304, 0.0411, 0.008, 0.0400, 0.008, 0.0310, 0.0411,
0.008, 0.0299, 0.008, 0.0370, 0.008, 0.0404, 0.008],
atol=allowed_distance_error)
assert_allclose(source_detector_distances(raw.copy().
pick("S1_D1 760").info),
[0.0304], atol=allowed_distance_error)
assert_allclose(source_detector_distances(raw.copy().
pick("S2_D2 760").info),
[0.0400], atol=allowed_distance_error)

# Test location of detectors
# The locations of detectors can be seen in the first
Expand All @@ -261,10 +283,6 @@ def test_snirf_nirsport2_w_positions():
assert_allclose(
mni_locs[2], [-0.0841, -0.0138, 0.0248], atol=allowed_dist_error)

assert raw.info['ch_names'][34][3:5] == 'D5'
assert_allclose(
mni_locs[34], [0.0845, -0.0451, -0.0123], atol=allowed_dist_error)

# Test location of sensors
# The locations of sensors can be seen in the second
# figure on this page...
Expand Down