MRG: Add nirx file reader

[rob-luke]

Add support for reading NIRX fNIRS files.

Reference issue

As discussed in #4257 and in person with @larsoner

What does this implement/fix?

Add mne/io/nirx. which allows reading of data collected from NIRX machine

Additional information

A rough list of things to be completed:

• Add basic reader that stores fNIRS (2 x channels x samples) data in practical way
• Add reader for source-detector locations and store in Raw (useful place to start?)
• Finalise FIFF decisions
• Add tests
• Add docs
• Ensure compliance with style guidelines and naming convention of mne
• Record short example measurement and add to mne test data set
• Update doc/_includes/data_formats.rst
• Reduce dependancies (scipy, pandas)
• Update doc/whats_new.rst

Demo

Example of functionality in use on real data can be found at https://github.com/rob-luke/mne-fnirs-demo/blob/master/NIRx%20Reader.ipynb

[larsoner]

What is the appropriate channel type to use in create_info. Currently I am using misc as hbo/hbr is for post processed measurements

What are the types / units of the data for these channels? We can add new channel types to the FIF standard if necessary.

[rob-luke]

What are the types / units of the data for these channels?

The NIRX devices are continuous wave, single phase machines. The signal measured and stored to disk is light intensity. That is what I am reading in this pull request.

The next processing step (beyond the scope of this pull request, but maybe relevant if we are discussing adding types) is to convert from intensity to optical density (script example1:nirstoolbox, example2:homer2) which is unitless.

The final preprocessing step is to convert optical density to concentration changes of HbO and HbR usually using the modified Beer Lambert law (script example 1:nirstoolbox, example 2:homer2). The units are usually change in μMol.

There is a nice overview of this in [1] that includes this nice graphic:

And a more in depth discussion can be found in [2] from which I have taken the following extract:

[1] Pinti, Paola, et al. "Current status and issues regarding pre-processing of fNIRS neuroimaging data: An investigation of diverse signal filtering methods within a General Linear Model framework." Frontiers in human neuroscience 12 (2018): 505.
[2] Scholkmann, Felix, et al. "A review on continuous wave functional near-infrared spectroscopy and imaging instrumentation and methodology." Neuroimage 85 (2014): 6-27.

[codecov]

Codecov Report

Merging #6674 into master will increase coverage by 0.01%.
The diff coverage is 93.77%.

@@            Coverage Diff             @@
##           master    #6674      +/-   ##
==========================================
+ Coverage   89.65%   89.67%   +0.01%     
==========================================
  Files         422      425       +3     
  Lines       76690    76898     +208     
  Branches    12546    12579      +33     
==========================================
+ Hits        68759    68956     +197     
- Misses       5116     5121       +5     
- Partials     2815     2821       +6

[larsoner]

Add basic reader that stores fNIRS (2 x channels x samples) data in practical way

I would store these as interleaved like ch_1_a, ch_1_b, ch_2_a, ch_2_b etc. Depending on what a and b actually store (both light intensities?) then they can be named something reasonable.

source-detector locations and store in Raw

info['chs'][ii]['loc'] is a 12-element float array. For EEG channels usually ['loc'][:3] is the location of the electrode and ['loc'][3:6] is the location of the reference electrode. Sounds like you could do something similar for fNIRS.

Record short example measurement and add to mne test data set

I would actually put this toward the top. You already have code in this PR that you could use to make travis green. It's nice to get it green and keep it green.

[larsoner]

Looks like it's shaping up, a few ideas for how to move forward

[rob-luke]

@larsoner I have hit a bit of a road block, could you please help me understand why the circleci checks are failing. Am I missing some docs for the class RawNIRX? I think the relevant line in the logs is:

/home/circleci/project/mne/io/nirx/nirx.py:docstring of mne.io.read_raw_nirx:30: WARNING: py:obj reference target not found: RawNIRX

[agramfort]

try editing https://github.com/mne-tools/mne-python/blob/master/doc/conf.py#L585

[rob-luke]

It's nice to get it green and keep it green.

Thanks @larsoner and @agramfort. All green now!

I will keep working on functionality and improving the code, and ping you when its ready for a review.

[agramfort]

Can you update https://github.com/mne-tools/mne-python/blob/master/doc/manual/io.rst?

[rob-luke]

Check location info is stored in accessible manner

Locations are now stored in info['chs'][ii]['loc'] with:

• first 3 entries as the x,y,z coordinates of the channel
• second 3 entries as the x,y,z coordinates of the source
• third 3 entries as the x,y,z coordinates of the detector

The locations are in MNI coordinates.

Plotting them using following script produces...

for idx in range(len(raw.info['chs'])):

    if idx == 1: label = 'Channel'; 
    else: label = ''
    ax.scatter(raw.info['chs'][idx]['loc'][0], 
               raw.info['chs'][idx]['loc'][1], 
               raw.info['chs'][idx]['loc'][2],
               marker = 'x', edgecolors='k', color='k', label=label)
    
    if idx == 1: label = 'Source'; 
    else: label = ''
    ax.scatter(raw.info['chs'][idx]['loc'][3], 
               raw.info['chs'][idx]['loc'][4], 
               raw.info['chs'][idx]['loc'][5],
               marker = 'o', edgecolors='b', color='b', label=label)
    
    if idx == 1: label = 'Detector'; 
    else: label = ''
    ax.scatter(raw.info['chs'][idx]['loc'][6], 
               raw.info['chs'][idx]['loc'][7], 
               raw.info['chs'][idx]['loc'][8],
               marker = 'o', edgecolors='r', color='r', label=label)

which matches nicely with the figure from mne-tools/mne-testing-data#51

If I change the type to eeg I can run raw.plot_sensors() ...

raw.plot_sensors(ch_type='eeg') # Force to treat as eeg so uses channel info stored in first 3 loc values

[agramfort]

@rob-luke we should see how to handle the different channel type for fnirs_raw
It's not that easy to do for example just because the FIFF constants are part of a
standard shared by different packages.

there is no way you can use an existing channel type?

[rob-luke]

Thanks for the feedback @agramfort. I am new to MNE and FIFF, so don’t know anything about channel types. Is there a type you would recommend I use?

Most fNIRS machines are continuous wave two wavelength devices. The raw data saved to file by the device is a measure of light intensity per wavelength. This needs to be converted to optical density. Finally the optical density is converted to chromophore concentration concentration. So if I wrote my own software from scratch I would add fnirs_raw, fnirs_od, and fnirs_chroma as types.

I see there is already a HbO and HbR channel type. These are chroma and only relevant after the conversion steps.

But to keep things simple could we store the channel type as fnirs and put a variable somewhere else that indicates if the data is raw (V), optical density (unitless), or chomophore (Mol)?

I am thinking out loud here, so please let me know if you have any suggestions.

[larsoner]

It's not that easy to do for example just because the FIFF constants are part of a
standard shared by different packages.

there is no way you can use an existing channel type?

In short, no, there doesn't seem to be. The existing "coil types" are actually derived from an earlier actual measured quantity (which this reader reads) so we need a new coil type for it. Thus I proposed working with mne-tools/fiff-constants to add a new FIFF constant (or multiple if we need units, etc.) in the discussion above.

[agramfort]

we should sort out the FIFF constant first then. @larsoner ok for you?

…

[larsoner]

I would let @rob-luke get this PR mostly ready to go to make sure we have all the constants that we need, then ping the MEGIN folks. I'd rather bother them once then multiple times.

[agramfort]

ok but we need to agree of the constants and units then. I don't have any strong opinion. Note also that as you introduce a new data type every pick_types function should be updated.

…

[larsoner]

also that as you introduce a new data type every pick_types function should
be updated.

We still only have FNIRS as the channel type, so no need for a new kwarg. But we do need to modify the fnirs=str options to have new options in addition to just 'hbo', 'hbr'. @rob-luke can you look and see everywhere we use e.g. 'hbo' and try to add your new raw optode measurement types?

FYI this PR is harder than most "add a reader" PR because we're adding new channel types, but it seems to be converging to something usable!

[larsoner]

BTW ideally we would with your reader add some new example/tutorial for dealing with fNIRS data showing that it works. (It also helps reviewers quickly test your new functionality.) If you have some simple paradigm (N100?) that you have or could record quickly to show that "things are working" that would be good as a new MNE dataset. Basically it can be a shorter version of #4257 (no need for info creation boilerplate, etc.), i.e. upload a new dataset to osf.io, add a downloader for it (mne.datasets.fnirs_demo.data_path()), and document it (doc/manual/datasets.rst). That would be a good follow-up PR if you're up for it.

That way, if/once the conversion steps are added, they can basically just be appended to this tutorial, and at the end of the day the fNIRS pipeline is be nicely documented for people.

[larsoner]

Hopefully you can just edit

https://github.com/mne-tools/mne-python/blob/master/mne/io/pick.py

And update the existing fnirs-picking tests like

https://github.com/mne-tools/mne-python/blob/master/mne/io/tests/test_pick.py#L260
https://github.com/mne-tools/mne-python/blob/master/mne/io/tests/test_pick.py#L487-L494

[rob-luke]

BTW ideally we would with your reader add some new example/tutorial for dealing with fNIRS data showing that it works. (It also helps reviewers quickly test your new functionality.) If you have some simple paradigm (N100?) that you have or could record quickly to show that "things are working" that would be good as a new MNE dataset. Basically it can be a shorter version of #4257 (no need for info creation boilerplate, etc.), i.e. upload a new dataset to osf.io, add a downloader for it (mne.datasets.fnirs_demo.data_path()), and document it (doc/manual/datasets.rst). That would be a good follow-up PR if you're up for it.

That way, if/once the conversion steps are added, they can basically just be appended to this tutorial, and at the end of the day the fNIRS pipeline is be nicely documented for people.

Already started something with the same idea here: https://github.com/rob-luke/mne-fnirs-demo/blob/master/NIRx%20Reader.ipynb
Will make a new PR once this is complete

[rob-luke]

Thanks for the help @larsoner and @agramfort . How about this plan?

• The physical measurement device is always an optode. So add FIFF.FIFFV_COIL_FNIRS_OPTODE and this is used for all stages of processing.
  Requires modification of fiff-constants and MNE.

  • Currently there is fnirs_hbo and fnirs_hbr as coil types. Maybe I am taking this too literally, but these are derived quantities, not a physical device. All the other coils listed seem to be physical. I propose depreciating this in a future PR.

• The type of data can be either raw data, optical density data, or chromophore data (explanation post). So add the channel types FIFF.FIFFV_FNIRS_RAW_CH, FIFF.FIFFV_FNIRS_OD_CH, and FIFF.FIFFV_FNIRS_CHROMA_CH.
  Requires modification to fiff-constants and MNE.

  • Currently there is FIFF.FIFFV_FNIRS_CH = 1100. This doesn't inform the user what type of data they are actually dealing with. I propose depreciating this in a future PR.

• Define the following get_channel_types()

  fnirs_raw=dict(coil_type=FIFF.FIFFV_COIL_FNIRS_OPTODE,
                kind=FIFF.FIFFV_FNIRS_RAW_CH,
                unit=FIFF.FIFF.FIFF_UNIT_V)
  
  fnirs_od=dict(coil_type=FIFF.FIFFV_COIL_FNIRS_OPTODE,
               kind=FIFF.FIFFV_FNIRS_OD_CH,
               unit=FIFF.FIFF.FIFF_UNIT_NONE)
  
  fnirs_chroma=dict(coil_type=FIFF.FIFFV_COIL_FNIRS_OPTODE,
                   kind=FIFF.FIFFV_FNIRS_CHROMA_CH,
                   unit=FIFF.FIFF.FIFF_UNIT_MOL)

• Store the precise details of type in the channel name.
  For fnirs_raw store the wavelength. E.g. raw.info['ch_names'][3] == "S1-D9 850 (nm)"
  For fnirs_od store the wavelength. E.g. raw.info['ch_names'][3] == "S1-D9 850 (nm)"
  For fnirs_chroma store the chromophore. E.g. raw.info['ch_names'][3] == "S1-D9 HbO"

  • Will need to provide mechanisms to extract out the detailed type of interest, for example raw.pick_wavelength() or raw.pick_chroma()

Thoughts? Suggestions? Approval to proceed with the plan?

(@hamishib any thoughts, details on logic, terminology, units etc?)

---

Additional information that may help make this decision.

• Different fNIRS systems use different light wavelengths, knowing the exact wavelength is required for conversion from optical density to chromorphore. Some devices use more than two wavelengths.
• Similarly some NIRS techniques allow extraction of chromophores beyond hbo/hbr (example), so flexibility to expand the list of chromophore would be good, but probably not essential.

[larsoner]

The physical measurement device is always an optode. So add FIFF.FIFFV_COIL_FNIRS_OPTODE and this is used for all stages of processing.

... add the channel types FIFF.FIFFV_FNIRS_RAW_CH, FIFF.FIFFV_FNIRS_OD_CH, and FIFF.FIFFV_FNIRS_CHROMA_CH.

This proposal IIUC is basically:

• Have three channel types: raw, optical density, and chroma
• Have one coil type that is used for all three channel types

It is a bit unusual to have one coil type used for multiple channel types. Something that might be closer to the MNE way of doing things (and also ends up backward compatible) is to do it the other way:

• Have all fNIRS-related channels be of type FIFFV_FNIRS_CH
• Have different "coil types" to differentiate the raw, optical density, and chroma "sensors"

I think this is closer to the way MNE operates on these quantities. This is true even though the units for the different coil types will end up differing. For example, in MEG we have magnetometers and gradiometers with different units (T and T/m), but these are all FIFFV_MEG_CH type with different coil types associated with them. And in some systems, channels can be operated on by a linear denoising comp matrix, and applying this yields a changed coil_type (setting some higher order bits).

Ultimately we could go either way, but to me the simplicity and backward compatibility of the unified FIFFV_FNIRS_CH approach gets my vote.

Store the precise details of type in the channel name.

This seems reasonable to me. I think we are limited to 14 chars for ch names. But I can't think of a better place to store this information. We can drop the nm and just make it part of the spec that the units are always nm (since this makes sense for infrared wavelength specification).

[agramfort]

plan suggested by @larsoner is ok for me. Except that encoding information in channel name seems dangerous. Is it possible that the wavelength is channel dependent or can it assumed global like sfreq?

…

[rob-luke]

plan suggested by @larsoner is ok for me.

Great, me too.

Except that encoding information in channel name seems dangerous.

Agreed, encoding information in the channel name sounds dangerous. And doesn't work well with the limit on characters.

Is it possible that the wavelength is channel dependent or can it assumed
global like sfreq?

I'm still learning MNE terminology here, so excuse if I have misunderstood your comment. For each source-detector pair there is a measurement for each wavelength. The NIRx machine uses 760 and 850 nm, but other machines use different (and possibly more than two) wavelengths.
Currently I am interleaving the wavelength on @larsoner suggestion, such that the channel names look something like: S1-D1 760, S1-D1 850, S1-D2 760, S1-D2 850, etc. The plan was that after conversion to chroma the channel names would become : S1-D1 HbO, S1-D1 HbR, S1-D2 HbO, etc.

What's an alternative that we could use? Place a field like info['fnirs_measurement'] = [760, 850] for raw data and have the channel names repeating like: S1-D1, S1-D1, S1-D2, S1-D2, etc. Then after conversion we modify the info structure to be info['fnirs_measurement'] = ["HbO", "HbR"]. I'm open to totally different suggestions too. Happy to do this whatever however best fits the MNE way. @larsoner any suggestion?

Side note

Maybe it will help guide some of these software design decisions for me to illustrate how fNIRS data is usually displayed in publications. For an fNIRS analysis we want to plot the haemodynamic response function (below, taken from great short description here). For each channel (source-detector pair) we want to plot the change in chromophore concentration after a stimulus onset. It would be great if after reading the data, conversion to chromophore, and epoching, that we could call epochs.plot(channel_index) and it would plot the HbO in red, and HbR in blue (each chromophore has a specific color it should be plotted in).

[rob-luke]

Epoching and plotting now seems to work (example). Now back to getting those test ticks green.

[lgtm-com]

This pull request introduces 3 alerts when merging 2bd56e9 into d247e55 - view on LGTM.com

new alerts:

• 3 for Unhashable object hashed

[rob-luke]

@larsoner I think I have made all the changes you requested. Please let me know if I need to do anything else.

Could you also point me to a function that changes one coil type to another? So I can use this as the basis for the conversion to optical density.

[larsoner]

Could you also point me to a function that changes one coil type to another? So I can use this as the basis for the conversion to optical density.

There is no real magic it's just a dict assignment. mne.io.Raw.fix_mag_coil_types does it, as does mne.io.Raw.apply_gradient_compensation (under the hood).

FYI -- we try to stay close to Python builtin types in MNE where possible to make this stuff painless. For example Info is just a subclass of dict with a nice __repr__ and a couple of convenience methods. info['chs'][ii] as you have seen is also just a dict. It makes life a lot easier in cases like this.

Will look at code tomorrow and hopefully merge!

[larsoner]

So here are the things from above I'll try to restore:

1. switch from Pandas to NumPy (and remove related decorators)

2. remove logger.info call from test (usually we don't do this and it looked like debug leftover anyway)

3. remove scipy.spatial.distance function in favor of just np.linalg.norm, and removed TODO label since this should be plenty fast (and explicit enough) already

4. moved ..externals.pymatreader to the top of the method that uses it (makes it a bit easier to trace dependencies when nested imports are at the top)

5. fix indexing in _read_raw_segment -- it wasn't taking into account the requested idx when storing the data in data. In principle it could be made more efficient, but so could the original reading (which will probably be the bottleneck) but for now at least it should work.

6. This problem (along with one with subject['age']) was exposed by running our standard tests with _test_raw_reader. For now subject['age'] is commented out. It needs to be a julian calendar date, see:

   https://www.nmr.mgh.harvard.edu/mne/stable/generated/mne.Info.html

   birthday : tuple of int
   Birthday in (year, month, day) format.

   If you want, you can just ballpark this by taking the recording date and subtracting their age from the year (and saying in the Notes that this is as close as we get to storing the actual birthdate for this format).

I had deleted my branch locally and never pushed to my origin, but this SO thread had some tips. Nothing in my git reflog for some reason, but git reflog rob-luke/add-nirx-reader had five entries:

$ git reflog rob-luke/add-nirx-reader
b02600eba (rob-luke/add-nirx-reader) refs/remotes/rob-luke/add-nirx-reader@{0}: fetch rob-luke: forced-update
cfba17b4a refs/remotes/rob-luke/add-nirx-reader@{1}: update by push
7e2a7a436 refs/remotes/rob-luke/add-nirx-reader@{2}: update by push
111b3d659 refs/remotes/rob-luke/add-nirx-reader@{3}: update by push
935c4c850 refs/remotes/rob-luke/add-nirx-reader@{4}: fetch rob-luke: storing head

Then doing:

$ git checkout cfba17b4a
...
HEAD is now at cfba17b4a FIX: More efficient reading

Which is awesome because this is the last commit that I think I pushed. Then a simple git log showed me three commits I made, which I'm going to cherry-pick, squash, and push so you can look @rob-luke.

TL;DR: You somehow omitted up my commits during rebase (probably forgot to pull them before rebasing/adding commits, then rebased, then force-pushed) but git magic allowed me to recover them.

[larsoner]

@rob-luke let me know if you are okay with 5b1e6ef and we can merge (assuming CIs are green, if they are not I'll fix them since I broke them)

[rob-luke]

Thanks for catching my mistake @larsoner, I had incorrectly assumed your commits were in there. All good by me to merge.

[larsoner]

Awesome, go go go @rob-luke !

[rob-luke]

Thanks @larsoner and @agramfort for all your patience and help. Very excited to have merged my first PR.

[larsoner]

Store the precise details of type in the channel name.

This seems reasonable to me. I think we are limited to 14 chars for ch names. But I can't think of a better place to store this information.

@rob-luke sorry to come back to this after we have something working -- but it occurred to me that this can probably also just live in info['chs'][ii]['loc']. We have 12 floats to work with there and we are currently just using 9 of them (channel X/Y/Z, source X/Y/Z, detector X/Y/Z). So why don't we add one more entry to go from 9 to 10 elements stored/used in as:

• channel X, channel Y, channel Z (3)
• source X, source Y, source Z (3)
• detector X, detector Y, detector Z (3)
• wavelength in nm or even just meters (1)

? Then we don't have to put the wavelength in the channel name or worry about parsing it, lengths, etc.

[agramfort]

brillant idea !

…

[rob-luke]

Great idea. I will get on to it once I clear the current backlog of PRs.

[massich]

I like it too

…

On Mon, Oct 28, 2019, 22:36 Robert Luke ***@***.***> wrote: Great idea. I will get on to it once I clear the current backlog of PRs. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#6674?email_source=notifications&email_token=ABVX5YY3Z3UZ5WXHYBQKAV3QQ5LNLA5CNFSM4IMUZEE2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOECOPGBQ#issuecomment-547156742>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ABVX5Y5N7LE54KGLX6DZLQDQQ5LNLANCNFSM4IMUZEEQ> .

[rob-luke]

@larsoner I just took a first crack at this and I get an error that the channel names are not unique if I remove the frequency from them. This might not make much sense without the code changes that I made, but the error is below. Is having unique channel names something that is essential for some functions to work in MNE?

mne/io/nirx/nirx.py:39: in read_raw_nirx
    return RawNIRX(fname, preload, verbose)
</Repositories/mne-python/mne/externals/decorator.py:decorator-gen-193>:2: in __init__
    ???
mne/utils/_logging.py:90: in wrapper
    return function(*args, **kwargs)
mne/io/nirx/nirx.py:208: in __init__
    ch_types='fnirs_raw')
<//Repositories/mne-python/mne/externals/decorator.py:decorator-gen-27>:2: in create_info
    ???
mne/utils/_logging.py:90: in wrapper
    return function(*args, **kwargs)
mne/io/meas_info.py:1799: in create_info
    info._check_consistency()
mne/io/meas_info.py:615: in _check_consistency
    self['ch_names'] = _unique_channel_names(self['ch_names'])
mne/io/meas_info.py:146: in _unique_channel_names
    '%s. Applying running numbers for duplicates.' % dups)
_ _ _ _ _ _ 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

message = "Channel names are not unique, found duplicates for: {'S2-D3', 'S5-D5', 'S4-D4', 'S5-D8', 'S5-D6', 'S1-D1', 'S5-D7', 'S5-D13', 'S3-D11', 'S1-D9', 'S4-D12', 'S3-D2', 'S2-D10'}. Applying running numbers for duplicates."
category = <class 'RuntimeWarning'>, module = 'mne'

[agramfort]

arfff true. Channel names must be different with mne. I would keep the names different but don't parse the wavelength from the channel name but from loc

[rob-luke]

Makes sense, I'll give that a go. Thanks