MRG: Store NIRx location data in meters

doc/changes/latest.inc

@@ -14,6 +14,7 @@ Current (0.20.dev0)
 
 - Add reader for NIRx data in :func:`mne.io.read_raw_nirx` by `Robert Luke`_
 
+- Add support for plotting fNIRS channels in :func:`mne.viz.plot_alignment` by `Eric Larson`_
 
 Changelog
 ~~~~~~~~~
@@ -23,4 +24,3 @@ Bug
 
 API
 ~~~
-

mne/coreg.py

@@ -1209,7 +1209,8 @@ def _scale_xfm(subject_to, xfm_fname, mri_name, subject_from, scale,
     _write_fs_xfm(fname_to, T_ras_mni['trans'], kind)
 
 
-def get_mni_fiducials(subject, subjects_dir=None):
+@verbose
+def get_mni_fiducials(subject, subjects_dir=None, verbose=None):
     """Estimate fiducials for a subject.
 
     Parameters
@@ -1219,6 +1220,7 @@ def get_mni_fiducials(subject, subjects_dir=None):
     subjects_dir : None | str
         Override the SUBJECTS_DIR environment variable
         (sys.environ['SUBJECTS_DIR'])
+    %(verbose)s
 
     Returns
     -------
@@ -1249,6 +1251,8 @@ def get_mni_fiducials(subject, subjects_dir=None):
     # Read fsaverage fiducials file and subject Talairach.
     fids, coord_frame = read_fiducials(fname_fids_fs)
     assert coord_frame == FIFF.FIFFV_COORD_MRI
+    if subject == 'fsaverage':
+        return fids  # special short-circuit for fsaverage
     mni_mri_t = invert_transform(_read_talxfm(subject, subjects_dir))
 
     # Convert to mm since this is Freesurfer's unit.

mne/defaults.py

@@ -52,6 +52,7 @@
         eeg_scale=4e-3, eegp_scale=20e-3, eegp_height=0.1,
         ecog_scale=5e-3,
         seeg_scale=5e-3,
+        fnirs_scale=5e-3,
         hpi_scale=15e-3,
 
         head_color=(0.988, 0.89, 0.74),
@@ -60,6 +61,7 @@
         eeg_color=(1., 0.596, 0.588), eegp_color=(0.839, 0.15, 0.16),
         ecog_color=(1., 1., 1.),
         seeg_color=(1., 1., .3),
+        fnirs_color=(1., .4, .3),
         lpa_color=(1., 0., 0.),
         nasion_color=(0., 1., 0.),
         rpa_color=(0., 0., 1.),

mne/io/nirx/nirx.py

@@ -9,10 +9,11 @@
 import numpy as np
 
 from ..base import BaseRaw
-from ..meas_info import create_info
-from ...utils import logger, verbose, fill_doc
 from ..constants import FIFF
+from ..meas_info import create_info, _format_dig_points
 from ...annotations import Annotations
+from ...transforms import apply_trans, _get_trans
+from ...utils import logger, verbose, fill_doc
 
 
 @fill_doc
@@ -57,6 +58,7 @@ class RawNIRX(BaseRaw):
     @verbose
     def __init__(self, fname, preload=False, verbose=None):
         from ...externals.pymatreader import read_mat
+        from ...coreg import get_mni_fiducials  # avoid circular import prob
         logger.info('Loading %s' % fname)
 
         # Check if required files exist and store names for later use
@@ -149,9 +151,31 @@ def __init__(self, fname, preload=False, verbose=None):
         #   Channels are defined as the midpoint between source and detector
         mat_data = read_mat(files['probeInfo.mat'], uint16_codec=None)
         requested_channels = mat_data['probeInfo']['probes']['index_c']
-        src_locs = mat_data['probeInfo']['probes']['coords_s3'] * 10
-        det_locs = mat_data['probeInfo']['probes']['coords_d3'] * 10
-        ch_locs = mat_data['probeInfo']['probes']['coords_c3'] * 10
+        src_locs = mat_data['probeInfo']['probes']['coords_s3'] / 100.
+        det_locs = mat_data['probeInfo']['probes']['coords_d3'] / 100.
+        ch_locs = mat_data['probeInfo']['probes']['coords_c3'] / 100.
+
+        # These are all in MNI coordinates, so let's transform them to
+        # the Neuromag head coordinate frame
+        mri_head_t, _ = _get_trans('fsaverage', 'mri', 'head')
+        src_locs = apply_trans(mri_head_t, src_locs)
+        det_locs = apply_trans(mri_head_t, det_locs)
+        ch_locs = apply_trans(mri_head_t, ch_locs)
+
+        # Set up digitization
+        dig = get_mni_fiducials('fsaverage', verbose=False)
+        for fid in dig:
+            fid['r'] = apply_trans(mri_head_t, fid['r'])
+            fid['coord_frame'] = FIFF.FIFFV_COORD_HEAD
+        for ii, ch_loc in enumerate(ch_locs, 1):
+            dig.append(dict(
+                kind=FIFF.FIFFV_POINT_EEG,  # misnomer but probably okay
+                r=ch_loc,
+                ident=ii,
+                coord_frame=FIFF.FIFFV_COORD_HEAD,
+            ))
+        dig = _format_dig_points(dig)
+        del mri_head_t
 
         # Determine requested channel indices
         # The wl1 and wl2 files include all possible source - detector pairs.
@@ -180,7 +204,7 @@ def prepend(list, str):
         info = create_info(chnames,
                            samplingrate,
                            ch_types='fnirs_raw')
-        info.update({'subject_info': subject_info})
+        info.update(subject_info=subject_info, dig=dig)
 
         # Store channel, source, and detector locations
         # The channel location is stored in the first 3 entries of loc.
@@ -258,7 +282,7 @@ def _probe_distances(self):
                 for ch in self.info['chs']]
         return np.array(dist, float)
 
-    def _short_channels(self, threshold=10.0):
+    def _short_channels(self, threshold=0.01):
         """Return a vector indicating which channels are short.
 
         Channels with distance less than `threshold` are reported as short.

mne/io/nirx/tests/test_nirx.py

@@ -4,12 +4,13 @@
 
 import os.path as op
 
-import numpy as np
+from numpy.testing import assert_allclose, assert_array_equal
 
+from mne.datasets.testing import data_path, requires_testing_data
 from mne.io import read_raw_nirx
 from mne.io.tests.test_raw import _test_raw_reader
+from mne.transforms import apply_trans, _get_trans
 from mne.utils import run_tests_if_main
-from mne.datasets.testing import data_path, requires_testing_data
 
 fname_nirx = op.join(data_path(download=False),
                      'NIRx', 'nirx_15_2_recording_w_short')
@@ -25,58 +26,79 @@ def test_nirx():
     assert raw.info['sfreq'] == 12.5
 
     # Test channel naming
-    assert raw.info['ch_names'][0] == "S1-D1 760"
-    assert raw.info['ch_names'][1] == "S1-D1 850"
-    assert raw.info['ch_names'][2] == "S1-D9 760"
-    assert raw.info['ch_names'][3] == "S1-D9 850"
-    assert raw.info['ch_names'][24] == "S5-D13 760"
-    assert raw.info['ch_names'][25] == "S5-D13 850"
+    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"]
 
     # Test info import
-    assert raw.info['subject_info']['sex'] == 1
-    assert raw.info['subject_info']['first_name'] == "MNE"
-    assert raw.info['subject_info']['middle_name'] == "Test"
-    assert raw.info['subject_info']['last_name'] == "Recording"
+    assert raw.info['subject_info'] == dict(sex=1, first_name="MNE",
+                                            middle_name="Test",
+                                            last_name="Recording")
 
     # Test distance between optodes matches values from
     # nirsite https://github.com/mne-tools/mne-testing-data/pull/51
     # step 4 figure 2
-    allowed_distance_error = 0.2
+    allowed_distance_error = 0.0002
     distances = raw._probe_distances()
-    assert abs(distances[0] - 30.4) < allowed_distance_error
-    assert abs(distances[2] - 7.8) < allowed_distance_error
-    assert abs(distances[4] - 31.0) < allowed_distance_error
-    assert abs(distances[6] - 8.6) < allowed_distance_error
-    assert abs(distances[8] - 41.6) < allowed_distance_error
-    assert abs(distances[10] - 7.2) < allowed_distance_error
-    assert abs(distances[12] - 38.9) < allowed_distance_error
-    assert abs(distances[14] - 7.5) < allowed_distance_error
-    assert abs(distances[16] - 55.8) < allowed_distance_error
-    assert abs(distances[18] - 56.2) < allowed_distance_error
-    assert abs(distances[20] - 56.1) < allowed_distance_error
-    assert abs(distances[22] - 56.5) < allowed_distance_error
-    assert abs(distances[24] - 7.7) < allowed_distance_error
+    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)
 
     # 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
     short_channels = raw._short_channels()
-    assert short_channels[0] is np.False_
-    assert short_channels[2] is np.True_
-    assert short_channels[4] is np.False_
-    assert short_channels[6] is np.True_
-    assert short_channels[8] is np.False_
-    short_channels = raw._short_channels(threshold=3)
-    assert short_channels[0] is np.False_
-    assert short_channels[2] is np.False_
+    assert_array_equal(short_channels[:9:2], [False, True, False, True, False])
+    short_channels = raw._short_channels(threshold=0.003)
+    assert_array_equal(short_channels[:3:2], [False, False])
     short_channels = raw._short_channels(threshold=50)
-    assert short_channels[0] is np.True_
-    assert short_channels[2] is np.True_
+    assert_array_equal(short_channels[:3:2], [True, True])
 
     # Test trigger events
-    assert raw.annotations[0]['description'] == '3.0'
-    assert raw.annotations[1]['description'] == '2.0'
-    assert raw.annotations[2]['description'] == '1.0'
+    assert_array_equal(raw.annotations.description, ['3.0', '2.0', '1.0'])
+
+    # Test location of detectors
+    # The locations of detectors can be seen in the first
+    # figure on this page...
+    # https://github.com/mne-tools/mne-testing-data/pull/51
+    # And have been manually copied below
+    # These values were reported in mm, but according to this page...
+    # https://mne.tools/stable/auto_tutorials/intro/plot_40_sensor_locations.html
+    # 3d locations should be specified in meters, so that's what's tested below
+    # Detector locations are stored in the third three loc values
+    allowed_dist_error = 0.0002
+    locs = [ch['loc'][6:9] for ch in raw.info['chs']]
+    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_allclose(
+        mni_locs[0], [-0.0841, -0.0464, -0.0129], atol=allowed_dist_error)
+
+    assert raw.info['ch_names'][4][3:5] == 'D3'
+    assert_allclose(
+        mni_locs[4], [0.0846, -0.0142, -0.0156], atol=allowed_dist_error)
+
+    assert raw.info['ch_names'][8][3:5] == 'D2'
+    assert_allclose(
+        mni_locs[8], [0.0207, -0.1062, 0.0484], atol=allowed_dist_error)
+
+    assert raw.info['ch_names'][12][3:5] == 'D4'
+    assert_allclose(
+        mni_locs[12], [-0.0196, 0.0821, 0.0275], atol=allowed_dist_error)
+
+    assert raw.info['ch_names'][16][3:5] == 'D5'
+    assert_allclose(
+        mni_locs[16], [-0.0360, 0.0276, 0.0778], atol=allowed_dist_error)
+
+    assert raw.info['ch_names'][19][3:5] == 'D6'
+    assert_allclose(
+        mni_locs[19], [0.0352, 0.0283, 0.0780], atol=allowed_dist_error)
+
+    assert raw.info['ch_names'][21][3:5] == 'D7'
+    assert_allclose(
+        mni_locs[21], [0.0388, -0.0477, 0.0932], atol=allowed_dist_error)
 
 
 @requires_testing_data

mne/viz/_3d.py

@@ -517,7 +517,7 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
                    surfaces='head', coord_frame='head',
                    meg=None, eeg='original', fwd=None,
                    dig=False, ecog=True, src=None, mri_fiducials=False,
-                   bem=None, seeg=True, show_axes=False, fig=None,
+                   bem=None, seeg=True, fnirs=True, show_axes=False, fig=None,
                    interaction='trackball', verbose=None):
     """Plot head, sensor, and source space alignment in 3D.
 
@@ -580,6 +580,10 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
         the subjects bem and bem/flash folders are searched. Defaults to None.
     seeg : bool
         If True (default), show sEEG electrodes.
+    fnirs : bool
+        If True (default), show fNIRS electrodes.
+
+        .. versionadded:: 0.20
     show_axes : bool
         If True (default False), coordinate frame axis indicators will be
         shown:
@@ -701,9 +705,11 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
     ref_meg = 'ref' in meg
     meg_picks = pick_types(info, meg=True, ref_meg=ref_meg)
     eeg_picks = pick_types(info, meg=False, eeg=True, ref_meg=False)
-    ecog_picks = pick_types(info, meg=False, ecog=True, ref_meg=False)
-    seeg_picks = pick_types(info, meg=False, seeg=True, ref_meg=False)
-
+    other_bools = dict(ecog=ecog, seeg=seeg, fnirs=fnirs)
+    del ecog, seeg, fnirs
+    other_keys = sorted(other_bools.keys())
+    other_picks = {key: pick_types(info, meg=False, ref_meg=False,
+                                   **{key: True}) for key in other_keys}
     if trans == 'auto':
         # let's try to do this in MRI coordinates so they're easy to plot
         subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
@@ -942,13 +948,12 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
 
     # determine points
     meg_rrs, meg_tris = list(), list()
-    ecog_loc = list()
-    seeg_loc = list()
     hpi_loc = list()
     ext_loc = list()
     car_loc = list()
     eeg_loc = list()
     eegp_loc = list()
+    other_loc = {key: list() for key in other_keys}
     if len(eeg) > 0:
         eeg_loc = np.array([info['chs'][k]['loc'][:3] for k in eeg_picks])
         if len(eeg_loc) > 0:
@@ -1006,12 +1011,12 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
         if len(car_loc) == len(ext_loc) == len(hpi_loc) == 0:
             warn('Digitization points not found. Cannot plot digitization.')
     del dig
-    if len(ecog_picks) > 0 and ecog:
-        ecog_loc = np.array([info['chs'][pick]['loc'][:3]
-                             for pick in ecog_picks])
-    if len(seeg_picks) > 0 and seeg:
-        seeg_loc = np.array([info['chs'][pick]['loc'][:3]
-                             for pick in seeg_picks])
+    for key, picks in other_picks.items():
+        if other_bools[key] and len(picks):
+            other_loc[key] = np.array([info['chs'][pick]['loc'][:3]
+                                       for pick in picks])
+            logger.info('Plotting %d %s location%s'
+                        % (len(other_loc[key]), key, _pl(other_loc[key])))
 
     # initialize figure
     renderer = _Renderer(fig, bgcolor=(0.5, 0.5, 0.5), size=(800, 800))
@@ -1054,20 +1059,19 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
     defaults = DEFAULTS['coreg']
     datas = [eeg_loc,
              hpi_loc,
-             ext_loc, ecog_loc, seeg_loc]
+             ext_loc] + list(other_loc[key] for key in other_keys)
     colors = [defaults['eeg_color'],
               defaults['hpi_color'],
-              defaults['extra_color'],
-              defaults['ecog_color'],
-              defaults['seeg_color']]
+              defaults['extra_color']
+              ] + [defaults[key + '_color'] for key in other_keys]
     alphas = [0.8,
               0.5,
-              0.25, 0.8, 0.8]
+              0.25] + [0.8] * len(other_keys)
     scales = [defaults['eeg_scale'],
               defaults['hpi_scale'],
-              defaults['extra_scale'],
-              defaults['ecog_scale'],
-              defaults['seeg_scale']]
+              defaults['extra_scale']
+              ] + [defaults[key + '_scale'] for key in other_keys]
+    assert len(datas) == len(colors) == len(alphas) == len(scales)
     for kind, loc in (('dig', car_loc), ('mri', fid_loc)):
         if len(loc) > 0:
             datas.extend(loc[:, np.newaxis])

mne/viz/tests/test_3d.py

@@ -21,7 +21,8 @@
                  setup_volume_source_space, read_forward_solution,
                  VolVectorSourceEstimate, convert_forward_solution,
                  compute_source_morph)
-from mne.io import read_raw_ctf, read_raw_bti, read_raw_kit, read_info
+from mne.io import (read_raw_ctf, read_raw_bti, read_raw_kit, read_info,
+                    read_raw_nirx)
 from mne._digitization._utils import write_dig
 from mne.io.pick import pick_info
 from mne.io.constants import FIFF
@@ -46,6 +47,7 @@
                     'sample-oct-6-src.fif')
 dip_fname = op.join(data_dir, 'MEG', 'sample', 'sample_audvis_trunc_set1.dip')
 ctf_fname = op.join(data_dir, 'CTF', 'testdata_ctf.ds')
+nirx_fname = op.join(data_dir, 'NIRx', 'nirx_15_2_recording_w_short')
 
 io_dir = op.join(op.abspath(op.dirname(__file__)), '..', '..', 'io')
 base_dir = op.join(io_dir, 'tests', 'data')
@@ -328,6 +330,13 @@ def test_plot_alignment(tmpdir, renderer):
                    trans=trans_fname, fwd=fwd,
                    surfaces='white', coord_frame='head')
 
+    # fNIRS
+    info = read_raw_nirx(nirx_fname).info
+    with catch_logging() as log:
+        plot_alignment(info, subject='fsaverage', surfaces=(), verbose=True)
+    log = log.getvalue()
+    assert '26 fnirs locations' in log
+
     renderer._close_all()