diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index d025b4649ca..2959af4d429 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -187,6 +187,8 @@ Bugs - :func:`mne.find_events` now uses ``first_samp`` and not ``0`` for initial event when using ``initial_value`` (:gh:`10289`, by `Alex Gramfort`_) +- Fix bug with :func:`mne.channels.make_standard_montage` for ``'standard*'``, ``'mgh*'``, and ``'artinis*'`` montages where the points were incorrectly scaled and fiducials incorrectly set away from the correct values for use with the ``fsaverage`` subject (:gh:`10324` by `Eric Larson`_) + - Fix plotting bug in :ref:`ex-electrode-pos-2d` and make view look more natural in :ref:`ex-movement-detect` (:gh:`10313`, by `Alex Rockhill`_) API changes diff --git a/mne/channels/_standard_montage_utils.py b/mne/channels/_standard_montage_utils.py index bd5d3164f21..6683f393589 100644 --- a/mne/channels/_standard_montage_utils.py +++ b/mne/channels/_standard_montage_utils.py @@ -12,6 +12,7 @@ import xml.etree.ElementTree as ElementTree from .montage import make_dig_montage +from .._freesurfer import get_mni_fiducials from ..transforms import _sph_to_cart from ..utils import warn, _pl from . import __file__ as _CHANNELS_INIT_FILE @@ -93,12 +94,20 @@ def _mgh_or_standard(basename, head_size, coord_frame='unknown'): break ch_names_.append(line.strip(' ').strip('\n')) - pos = np.array(pos) + pos = np.array(pos) / 1000. ch_pos = _check_dupes_odict(ch_names_, pos) nasion, lpa, rpa = [ch_pos.pop(n) for n in fid_names] - scale = head_size / np.median(np.linalg.norm(pos, axis=1)) + if head_size is None: + scale = 1. + else: + scale = head_size / np.median(np.linalg.norm(pos, axis=1)) for value in ch_pos.values(): value *= scale + # if we are in MRI/MNI coordinates, we need to replace nasion, LPA, and RPA + # with those of fsaverage for ``trans='fsaverage'`` to work + if coord_frame == 'mri': + lpa, nasion, rpa = [ + x['r'].copy() for x in get_mni_fiducials('fsaverage')] nasion *= scale lpa *= scale rpa *= scale @@ -130,20 +139,26 @@ def _mgh_or_standard(basename, head_size, coord_frame='unknown'): 'biosemi32': partial(_biosemi, basename='biosemi32.txt'), 'biosemi64': partial(_biosemi, basename='biosemi64.txt'), - 'mgh60': partial(_mgh_or_standard, basename='mgh60.elc'), - 'mgh70': partial(_mgh_or_standard, basename='mgh70.elc'), + 'mgh60': partial(_mgh_or_standard, basename='mgh60.elc', + coord_frame='mri'), + 'mgh70': partial(_mgh_or_standard, basename='mgh70.elc', + coord_frame='mri'), 'standard_1005': partial(_mgh_or_standard, - basename='standard_1005.elc'), + basename='standard_1005.elc', coord_frame='mri'), 'standard_1020': partial(_mgh_or_standard, - basename='standard_1020.elc'), + basename='standard_1020.elc', coord_frame='mri'), 'standard_alphabetic': partial(_mgh_or_standard, - basename='standard_alphabetic.elc'), + basename='standard_alphabetic.elc', + coord_frame='mri'), 'standard_postfixed': partial(_mgh_or_standard, - basename='standard_postfixed.elc'), + basename='standard_postfixed.elc', + coord_frame='mri'), 'standard_prefixed': partial(_mgh_or_standard, - basename='standard_prefixed.elc'), + basename='standard_prefixed.elc', + coord_frame='mri'), 'standard_primed': partial(_mgh_or_standard, - basename='standard_primed.elc'), + basename='standard_primed.elc', + coord_frame='mri'), 'artinis-octamon': partial(_mgh_or_standard, coord_frame='mri', basename='artinis-octamon.elc'), 'artinis-brite23': partial(_mgh_or_standard, coord_frame='mri', diff --git a/mne/channels/montage.py b/mne/channels/montage.py index cd367beb59b..db1815211c5 100644 --- a/mne/channels/montage.py +++ b/mne/channels/montage.py @@ -1483,16 +1483,19 @@ def compute_native_head_t(montage): return Transform(coord_frame, 'head', native_head_t) -def make_standard_montage(kind, head_size=HEAD_SIZE_DEFAULT): +def make_standard_montage(kind, head_size='auto'): """Read a generic (built-in) montage. Parameters ---------- kind : str The name of the montage to use. See notes for valid kinds. - head_size : float + head_size : float | None | str The head size (radius, in meters) to use for spherical montages. - Defaults to 95mm. + Can be None to not scale the read sizes. ``'auto'`` (default) will + use 95mm for all montages except the ``'standard*'``, ``'mgh*'``, and + ``'artinis*'``, which are already in fsaverage's MRI coordinates + (same as MNI). Returns ------- @@ -1566,7 +1569,15 @@ def make_standard_montage(kind, head_size=HEAD_SIZE_DEFAULT): .. versionadded:: 0.19.0 """ from ._standard_montage_utils import standard_montage_look_up_table + _validate_type(kind, str, 'kind') _check_option('kind', kind, _BUILT_IN_MONTAGES) + _validate_type(head_size, ('numeric', str, None), 'head_size') + if isinstance(head_size, str): + _check_option('head_size', head_size, ('auto',), extra='when str') + if kind.startswith(('standard', 'mgh', 'artinis')): + head_size = None + else: + head_size = HEAD_SIZE_DEFAULT return standard_montage_look_up_table[kind](head_size=head_size) diff --git a/mne/channels/tests/test_montage.py b/mne/channels/tests/test_montage.py index d0b37b00603..fd6dad3a695 100644 --- a/mne/channels/tests/test_montage.py +++ b/mne/channels/tests/test_montage.py @@ -35,7 +35,8 @@ from mne.io._digitization import (_format_dig_points, _get_fid_coords, _get_dig_eeg, _count_points_by_type) -from mne.transforms import _ensure_trans, apply_trans, invert_transform +from mne.transforms import (_ensure_trans, apply_trans, invert_transform, + _get_trans) from mne.viz._3d import _fiducial_coords from mne.io.kit import read_mrk @@ -1051,6 +1052,7 @@ def test_set_montage_mgh(rename): orig_pos = np.array([raw.info['chs'][pick]['loc'][:3] for pick in eeg_picks]) atol = 1e-6 + mon = None if rename == 'raw': raw.rename_channels(lambda x: x.replace('EEG ', 'EEG')) raw.set_montage('mgh60') # test loading with string argument @@ -1075,15 +1077,29 @@ def renamer(x): mon.rename_channels(renamer) raw.set_montage(mon) + if mon is not None: + # first two are 'Fz' and 'F2', take them from standard_1020.elc -- + # they should not be changed on load! + want_pos = [[0.3122, 58.5120, 66.4620], [29.5142, 57.6019, 59.5400]] + got_pos = [mon.get_positions()['ch_pos'][f'EEG {x:03d}'] * 1000 + for x in range(1, 3)] + assert_allclose(want_pos, got_pos) + assert mon.dig[0]['coord_frame'] == FIFF.FIFFV_COORD_MRI + trans = compute_native_head_t(mon) + trans_2 = _get_trans('fsaverage', 'mri', 'head')[0] + assert trans['to'] == trans_2['to'] + assert trans['from'] == trans_2['from'] + assert_allclose(trans['trans'], trans_2['trans'], atol=1e-6) + new_pos = np.array([ch['loc'][:3] for ch in raw.info['chs'] if ch['ch_name'].startswith('EEG')]) assert ((orig_pos != new_pos).all()) r0 = _fit_sphere(new_pos)[1] - assert_allclose(r0, [0.000775, 0.006881, 0.047398], atol=1e-3) + assert_allclose(r0, [-0.001021, 0.014554, 0.041404], atol=1e-4) # spot check - assert_allclose(new_pos[:2], [[0.000273, 0.084920, 0.105838], - [0.028822, 0.083529, 0.099164]], atol=atol) + assert_allclose(new_pos[:2], [[-0.001229, 0.093274, 0.102639], + [0.027968, 0.09187, 0.09578]], atol=atol) # XXX: this does not check ch_names + it cannot work because of write_dig diff --git a/mne/channels/tests/test_standard_montage.py b/mne/channels/tests/test_standard_montage.py index aaed9f1962b..9c85874e0bb 100644 --- a/mne/channels/tests/test_standard_montage.py +++ b/mne/channels/tests/test_standard_montage.py @@ -29,10 +29,11 @@ def test_standard_montages_have_fids(kind): for k, v in fids.items(): assert v is not None, k for d in montage.dig: - if kind == 'artinis-octamon' or kind == 'artinis-brite23': - assert d['coord_frame'] == FIFF.FIFFV_COORD_MRI + if kind.startswith(('artinis', 'standard', 'mgh')): + want = FIFF.FIFFV_COORD_MRI else: - assert d['coord_frame'] == FIFF.FIFFV_COORD_UNKNOWN + want = FIFF.FIFFV_COORD_UNKNOWN + assert d['coord_frame'] == want def test_standard_montage_errors(): @@ -142,12 +143,11 @@ def test_set_montage_artinis_fsaverage(kind): assert trans['from'] == trans_fs['from'] translation = 1000 * np.linalg.norm(trans['trans'][:3, 3] - trans_fs['trans'][:3, 3]) - # TODO: This is actually quite big... - assert 15 < translation < 18 # mm + assert 0 < translation < 1 # mm rotation = np.rad2deg( _angle_between_quats(rot_to_quat(trans['trans'][:3, :3]), rot_to_quat(trans_fs['trans'][:3, :3]))) - assert 3 < rotation < 7 # degrees + assert 0 < rotation < 1 # degrees def test_set_montage_artinis_basic(): @@ -172,15 +172,15 @@ def test_set_montage_artinis_basic(): # Check a known location assert_array_almost_equal(raw.info['chs'][0]['loc'][:3], - [0.0616, 0.075398, 0.07347]) + [0.054243, 0.081884, 0.054544]) assert_array_almost_equal(raw.info['chs'][8]['loc'][:3], - [-0.033875, 0.101276, 0.077291]) + [-0.03013, 0.105097, 0.055894]) assert_array_almost_equal(raw.info['chs'][12]['loc'][:3], - [-0.062749, 0.080417, 0.074884]) + [-0.055681, 0.086566, 0.055858]) assert_array_almost_equal(raw_od.info['chs'][12]['loc'][:3], - [-0.062749, 0.080417, 0.074884]) + [-0.055681, 0.086566, 0.055858]) assert_array_almost_equal(raw_hb.info['chs'][12]['loc'][:3], - [-0.062749, 0.080417, 0.074884]) + [-0.055681, 0.086566, 0.055858]) # Check that locations are identical for a pair of channels (all elements # except the 10th which is the wavelength if not hbo and hbr type) assert_array_almost_equal(raw.info['chs'][0]['loc'][:9], @@ -200,11 +200,11 @@ def test_set_montage_artinis_basic(): raw.info['chs'][0]['loc'][:9]) # Check a known location assert_array_almost_equal(raw.info['chs'][0]['loc'][:3], - [0.085583, 0.036275, 0.089426]) + [0.068931, 0.046201, 0.072055]) assert_array_almost_equal(raw.info['chs'][8]['loc'][:3], - [0.069555, 0.078579, 0.069305]) + [0.055196, 0.082757, 0.052165]) assert_array_almost_equal(raw.info['chs'][12]['loc'][:3], - [0.044861, 0.100952, 0.065175]) + [0.033592, 0.102607, 0.047423]) # Check that locations are identical for a pair of channels (all elements # except the 10th which is the wavelength if not hbo and hbr type) assert_array_almost_equal(raw.info['chs'][0]['loc'][:9], diff --git a/tutorials/forward/35_eeg_no_mri.py b/tutorials/forward/35_eeg_no_mri.py index 6ed765be630..e5036c9e907 100644 --- a/tutorials/forward/35_eeg_no_mri.py +++ b/tutorials/forward/35_eeg_no_mri.py @@ -58,7 +58,8 @@ for ch_name in raw.ch_names) raw.rename_channels(new_names) -# Read and set the EEG electrode locations: +# Read and set the EEG electrode locations, which are already in fsaverage's +# space (MNI space) for standard_1020: montage = mne.channels.make_standard_montage('standard_1005') raw.set_montage(montage) raw.set_eeg_reference(projection=True) # needed for inverse modeling