Skip to content
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

Jump to bottom

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

Merged
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 all 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.
Jump to
Jump to file
Failed to load files.
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 @@ -25,6 +25,8 @@ Enhancements
~~~~~~~~~~~~
- Add support for ahdr files in :func:`mne.io.read_raw_brainvision` (:gh:`10515` by :newcontrib:`Alessandro Tonin`)

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

- Add support for ``overview_mode`` in :meth:`raw.plot() <mne.io.Raw.plot>` and related functions/methods (:gh:`10501` by `Eric Larson`_)

- Add :meth:`mne.io.Raw.crop_by_annotations` method to get chunks of Raw data based on :class:`mne.Annotations`. (:gh:`10460` by `Alex Gramfort`_)
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
2 changes: 1 addition & 1 deletion mne/channels/tests/test_interpolation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -302,6 +302,6 @@ def test_interpolation_nirs():
assert bad_0_std_pre_interp > np.std(raw_od._data[bad_0])
raw_haemo = beer_lambert_law(raw_od, ppf=6)
raw_haemo.info['bads'] = raw_haemo.ch_names[2:4]
assert raw_haemo.info['bads'] == ['S1_D2 hbo', 'S1_D2 hbr']
assert raw_haemo.info['bads'] == ['S10_D11 hbo', 'S10_D11 hbr']
raw_haemo.interpolate_bads()
assert raw_haemo.info['bads'] == []
4 changes: 2 additions & 2 deletions mne/datasets/config.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@
# respective repos, and make a new release of the dataset on GitHub. Then
# update the checksum in the MNE_DATASETS dict below, and change version
# here: ↓↓↓↓↓ ↓↓↓
RELEASES = dict(testing='0.134', misc='0.23')
RELEASES = dict(testing='0.135', misc='0.23')
TESTING_VERSIONED = f'mne-testing-data-{RELEASES["testing"]}'
MISC_VERSIONED = f'mne-misc-data-{RELEASES["misc"]}'

Expand All @@ -111,7 +111,7 @@
# Testing and misc are at the top as they're updated most often
MNE_DATASETS['testing'] = dict(
archive_name=f'{TESTING_VERSIONED}.tar.gz', # 'mne-testing-data',
hash='md5:9d234ff0156e6feccbbbc01469a861eb',
hash='md5:db40791e786fa776e8717cf50b43b0ec',
url=('https://codeload.github.com/mne-tools/mne-testing-data/'
f'tar.gz/{RELEASES["testing"]}'),
folder_name='MNE-testing-data',
Expand Down
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)
Copy link
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.

Copy link
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...

Copy link
Member

@larsoner larsoner May 13, 2022 •
edited

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.

Copy link
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.

Copy link
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
130 changes: 68 additions & 62 deletions mne/io/nirx/tests/test_nirx.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,11 +95,12 @@ def test_nirsport_v2():
# 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],
[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 @@ -118,9 +119,9 @@ def test_nirsport_v2():
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 raw.info['ch_names'][13][3:5] == 'D5'
assert_allclose(
mni_locs[34], [0.0845, -0.0451, -0.0123], atol=allowed_dist_error)
mni_locs[13], [0.0845, -0.0451, -0.0123], atol=allowed_dist_error)

# Test location of sensors
# The locations of sensors can be seen in the second
Expand All @@ -138,9 +139,9 @@ def test_nirsport_v2():
assert_allclose(
mni_locs[9], [-0.0, -0.1195, 0.0142], atol=allowed_dist_error)

assert raw.info['ch_names'][34][:2] == 'S8'
assert raw.info['ch_names'][39][:2] == 'S8'
assert_allclose(
mni_locs[34], [0.0828, -0.046, 0.0285], atol=allowed_dist_error)
mni_locs[39], [0.0828, -0.046, 0.0285], atol=allowed_dist_error)

assert len(raw.annotations) == 3
assert raw.annotations.description[0] == '1.0'
Expand Down Expand Up @@ -290,7 +291,7 @@ def test_nirx_15_2_short():
# 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"]

# Test frequency encoding
assert raw.info['chs'][0]['loc'][9] == 760
Expand All @@ -307,17 +308,18 @@ def test_nirx_15_2_short():
# nirsite https://github.com/mne-tools/mne-testing-data/pull/51
# step 4 figure 2
allowed_distance_error = 0.0002
distances = source_detector_distances(raw.info)
assert_allclose(distances[::2], [
0.0304, 0.0078, 0.0310, 0.0086, 0.0416,
0.0072, 0.0389, 0.0075, 0.0558, 0.0562,
0.0561, 0.0565, 0.0077], 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_D10 760").info),
[0.0086], atol=allowed_distance_error)

# Test which channels are short
# These are the ones marked as red at
# https://github.com/mne-tools/mne-testing-data/pull/51 step 4 figure 2
is_short = short_channels(raw.info)
assert_array_equal(is_short[:9:2], [False, True, False, True, False])
assert_array_equal(is_short[:9:2], [False, True, True, False, True])
is_short = short_channels(raw.info, threshold=0.003)
assert_array_equal(is_short[:3:2], [False, False])
is_short = short_channels(raw.info, threshold=50)
Expand All @@ -344,29 +346,29 @@ def test_nirx_15_2_short():
assert_allclose(
mni_locs[0], [-0.0841, -0.0464, -0.0129], atol=allowed_dist_error)

assert raw.info['ch_names'][4][3:5] == 'D3'
assert raw.info['ch_names'][6][3:5] == 'D3'
assert_allclose(
mni_locs[4], [0.0846, -0.0142, -0.0156], atol=allowed_dist_error)
mni_locs[6], [0.0846, -0.0142, -0.0156], atol=allowed_dist_error)

assert raw.info['ch_names'][8][3:5] == 'D2'
assert raw.info['ch_names'][10][3:5] == 'D2'
assert_allclose(
mni_locs[8], [0.0207, -0.1062, 0.0484], atol=allowed_dist_error)
mni_locs[10], [0.0207, -0.1062, 0.0484], atol=allowed_dist_error)

assert raw.info['ch_names'][12][3:5] == 'D4'
assert raw.info['ch_names'][14][3:5] == 'D4'
assert_allclose(
mni_locs[12], [-0.0196, 0.0821, 0.0275], atol=allowed_dist_error)
mni_locs[14], [-0.0196, 0.0821, 0.0275], atol=allowed_dist_error)

assert raw.info['ch_names'][16][3:5] == 'D5'
assert raw.info['ch_names'][18][3:5] == 'D5'
assert_allclose(
mni_locs[16], [-0.0360, 0.0276, 0.0778], atol=allowed_dist_error)
mni_locs[18], [-0.0360, 0.0276, 0.0778], atol=allowed_dist_error)

assert raw.info['ch_names'][19][3:5] == 'D6'
assert raw.info['ch_names'][20][3:5] == 'D6'
assert_allclose(
mni_locs[19], [0.0352, 0.0283, 0.0780], atol=allowed_dist_error)
mni_locs[20], [0.0352, 0.0283, 0.0780], atol=allowed_dist_error)

assert raw.info['ch_names'][21][3:5] == 'D7'
assert raw.info['ch_names'][23][3:5] == 'D7'
assert_allclose(
mni_locs[21], [0.0388, -0.0477, 0.0932], atol=allowed_dist_error)
mni_locs[23], [0.0388, -0.0477, 0.0932], atol=allowed_dist_error)


@requires_testing_data
Expand All @@ -381,7 +383,7 @@ def test_nirx_15_3_short():
# Test channel naming
assert raw.info['ch_names'][:4] == ["S1_D2 760", "S1_D2 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"]

# Test frequency encoding
assert raw.info['chs'][0]['loc'][9] == 760
Expand All @@ -398,17 +400,18 @@ def test_nirx_15_3_short():
# Test distance between optodes matches values from
# https://github.com/mne-tools/mne-testing-data/pull/72
allowed_distance_error = 0.001
distances = source_detector_distances(raw.info)
assert_allclose(distances[::2], [
0.0304, 0.0078, 0.0310, 0.0086, 0.0416,
0.0072, 0.0389, 0.0075, 0.0558, 0.0562,
0.0561, 0.0565, 0.0077], atol=allowed_distance_error)
assert_allclose(source_detector_distances(raw.copy().
pick("S1_D2 760").info),
[0.0304], atol=allowed_distance_error)
assert_allclose(source_detector_distances(raw.copy().
pick("S5_D13 760").info),
[0.0076], atol=allowed_distance_error)

# Test which channels are short
# These are the ones marked as red at
# https://github.com/mne-tools/mne-testing-data/pull/72
is_short = short_channels(raw.info)
assert_array_equal(is_short[:9:2], [False, True, False, True, False])
assert_array_equal(is_short[:9:2], [False, True, False, True, True])
is_short = short_channels(raw.info, threshold=0.003)
assert_array_equal(is_short[:3:2], [False, False])
is_short = short_channels(raw.info, threshold=50)
Expand All @@ -435,25 +438,25 @@ def test_nirx_15_3_short():
assert_allclose(
mni_locs[4], [0.0846, -0.0142, -0.0156], atol=allowed_dist_error)

assert raw.info['ch_names'][8][3:5] == 'D3'
assert raw.info['ch_names'][10][3:5] == 'D3'
assert_allclose(
mni_locs[8], [0.0207, -0.1062, 0.0484], atol=allowed_dist_error)
mni_locs[10], [0.0207, -0.1062, 0.0484], atol=allowed_dist_error)

assert raw.info['ch_names'][12][3:5] == 'D4'
assert raw.info['ch_names'][14][3:5] == 'D4'
assert_allclose(
mni_locs[12], [-0.0196, 0.0821, 0.0275], atol=allowed_dist_error)
mni_locs[14], [-0.0196, 0.0821, 0.0275], atol=allowed_dist_error)

assert raw.info['ch_names'][16][3:5] == 'D5'
assert raw.info['ch_names'][18][3:5] == 'D5'
assert_allclose(
mni_locs[16], [-0.0360, 0.0276, 0.0778], atol=allowed_dist_error)
mni_locs[18], [-0.0360, 0.0276, 0.0778], atol=allowed_dist_error)

assert raw.info['ch_names'][19][3:5] == 'D6'
assert raw.info['ch_names'][20][3:5] == 'D6'
assert_allclose(
mni_locs[19], [0.0388, -0.0477, 0.0932], atol=allowed_dist_error)
mni_locs[20], [0.0388, -0.0477, 0.0932], atol=allowed_dist_error)

assert raw.info['ch_names'][21][3:5] == 'D7'
assert raw.info['ch_names'][22][3:5] == 'D7'
assert_allclose(
mni_locs[21], [-0.0394, -0.0483, 0.0928], atol=allowed_dist_error)
mni_locs[22], [-0.0394, -0.0483, 0.0928], atol=allowed_dist_error)


@requires_testing_data
Expand Down Expand Up @@ -501,8 +504,8 @@ def test_nirx_15_2():
tzinfo=dt.timezone.utc)

# Test channel naming
assert raw.info['ch_names'][:4] == ["S1_D1 760", "S1_D1 850",
"S1_D10 760", "S1_D10 850"]
assert raw.info['ch_names'][:4] == ["S10_D10 760", "S10_D10 850",
"S10_D9 760", "S10_D9 850"]

# Test info import
assert raw.info['subject_info'] == dict(sex=1, first_name="TestRecording",
Expand All @@ -519,13 +522,13 @@ def test_nirx_15_2():
head_mri_t, _ = _get_trans('fsaverage', 'head', 'mri')
mni_locs = apply_trans(head_mri_t, locs)

assert raw.info['ch_names'][0][3:5] == 'D1'
assert raw.info['ch_names'][28][3:5] == 'D1'
assert_allclose(
mni_locs[0], [-0.0292, 0.0852, -0.0142], atol=allowed_dist_error)
mni_locs[28], [-0.0292, 0.0852, -0.0142], atol=allowed_dist_error)

assert raw.info['ch_names'][15][3:5] == 'D4'
assert raw.info['ch_names'][42][3:5] == 'D4'
assert_allclose(
mni_locs[15], [-0.0739, -0.0756, -0.0075], atol=allowed_dist_error)
mni_locs[42], [-0.0739, -0.0756, -0.0075], atol=allowed_dist_error)

# Old name aliases for backward compat
assert 'fnirs_cw_amplitude' in raw
Expand All @@ -549,12 +552,12 @@ def test_nirx_15_0():
tzinfo=dt.timezone.utc)

# Test channel naming
assert raw.info['ch_names'][:12] == ["S1_D1 760", "S1_D1 850",
assert raw.info['ch_names'][:12] == ["S10_D10 760", "S10_D10 850",
"S1_D1 760", "S1_D1 850",
"S2_D2 760", "S2_D2 850",
"S3_D3 760", "S3_D3 850",
"S4_D4 760", "S4_D4 850",
"S5_D5 760", "S5_D5 850",
"S6_D6 760", "S6_D6 850"]
"S5_D5 760", "S5_D5 850"]

# Test info import
assert raw.info['subject_info'] == {'birthday': (2004, 10, 27),
Expand All @@ -572,20 +575,23 @@ def test_nirx_15_0():
head_mri_t, _ = _get_trans('fsaverage', 'head', 'mri')
mni_locs = apply_trans(head_mri_t, locs)

assert raw.info['ch_names'][0][3:5] == 'D1'
assert raw.info['ch_names'][2][3:5] == 'D1'
assert_allclose(
mni_locs[0], [0.0287, -0.1143, -0.0332], atol=allowed_dist_error)
mni_locs[2], [0.0287, -0.1143, -0.0332], atol=allowed_dist_error)

assert raw.info['ch_names'][15][3:5] == 'D8'
assert raw.info['ch_names'][17][3:5] == 'D8'
assert_allclose(
mni_locs[15], [-0.0693, -0.0480, 0.0657], atol=allowed_dist_error)
mni_locs[17], [-0.0693, -0.0480, 0.0657], atol=allowed_dist_error)

# Test distance between optodes matches values from
allowed_distance_error = 0.0002
distances = source_detector_distances(raw.info)
assert_allclose(distances[::2], [
0.0301, 0.0315, 0.0343, 0.0368, 0.0408,
0.0399, 0.0393, 0.0367, 0.0336, 0.0447], atol=allowed_distance_error)

assert_allclose(source_detector_distances(raw.copy().
pick("S1_D1 760").info),
[0.0300], atol=allowed_distance_error)
assert_allclose(source_detector_distances(raw.copy().
pick("S7_D7 760").info),
[0.0392], atol=allowed_distance_error)


@requires_testing_data
Expand Down