Skip to content

BUG: Fix bug with scaling and fids for standard montages #10324

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 4 commits into from
Feb 12, 2022
Merged

BUG: Fix bug with scaling and fids for standard montages #10324

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d9345d3
BUG: Fix bug with scaling and fids for standard montages
larsoner Feb 11, 2022
0f1f05f
FIX: Dont do that vscode
larsoner Feb 11, 2022
4224595
FIX: Test
larsoner Feb 11, 2022
a96ada5
FIX: Test
larsoner Feb 11, 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
2 changes: 2 additions & 0 deletions doc/changes/latest.inc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
35 changes: 25 additions & 10 deletions mne/channels/_standard_montage_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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',
Expand Down
17 changes: 14 additions & 3 deletions mne/channels/montage.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
-------
Expand Down Expand Up @@ -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)


Expand Down
24 changes: 20 additions & 4 deletions mne/channels/tests/test_montage.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
28 changes: 14 additions & 14 deletions mne/channels/tests/test_standard_montage.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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():
Expand Down Expand Up @@ -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
Comment on lines -145 to +146
Copy link
Copy Markdown
Member Author

@larsoner larsoner Feb 11, 2022

Choose a reason for hiding this comment

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

This is a sign that things work better (and that they were broken before by this amount)

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():
Expand All @@ -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])
Copy link
Copy Markdown
Member Author

@larsoner larsoner Feb 11, 2022

Choose a reason for hiding this comment

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

@HanBnrd @rob-luke are the artinis montages added in #9141 in MNI / fsaverage space? If so, these new positions should be more accurate. Could you check if plot_alignment looks more accurate, or at least as accurate as before? The differences might be subtle (things will move slightly more anterior)...

Copy link
Copy Markdown
Contributor

@HanBnrd HanBnrd Feb 11, 2022

Choose a reason for hiding this comment

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

Hi @larsoner, yes artinis-brite23.elc and artinis-octamon.elc are using MNI coordinates. Sorry, what would you like to check exactly? Would you like us to run plot_alignment with fNIRS data for the Artinis montages and compare before and after your PR?

Copy link
Copy Markdown
Member Author

@larsoner larsoner Feb 11, 2022

Choose a reason for hiding this comment

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

Would you like us to run plot_alignment with fNIRS data for the Artinis montages and compare before and after your PR?

Yeah that would be great. The changes suggest that they will shift anterior (and a little inferior) by ~1.5 cm but it should actually be more correct

Copy link
Copy Markdown
Contributor

@HanBnrd HanBnrd Feb 13, 2022

Choose a reason for hiding this comment

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

Hi @larsoner, this is before (left) and after (right) the PR with the Artinis Octamon montage
before_after
before_after_2

It seems that after the PR, the sensors look like they're inside the head in the 3D vis but this is a bit confusing because the 2D viz shows the sensors more anterior indeed. This is the code I used:

raw.set_montage('artinis-octamon')

raw.plot_sensors()

subjects_dir = mne.datasets.sample.data_path() + '/subjects'
mne.datasets.fetch_fsaverage(subjects_dir=subjects_dir, verbose=True)

fig = mne.viz.create_3d_figure(size=(800, 600), bgcolor='white')
fig = mne.viz.plot_alignment(raw.info, show_axes=True,
                             subject='fsaverage', coord_frame='mri',
                             trans='fsaverage', surfaces=['brain', 'head'],
                             fnirs=['channels', 'pairs',
                                    'sources', 'detectors'],
                             dig=True, mri_fiducials=True,
                             subjects_dir=subjects_dir, fig=fig)

mne.viz.set_3d_view(figure=fig, azimuth=70, elevation=100, distance=0.4,
                    focalpoint=(0., -0.01, 0.02))

Copy link
Copy Markdown
Member Author

@larsoner larsoner Feb 14, 2022

Choose a reason for hiding this comment

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

It seems that after the PR, the sensors look like they're inside the head in the 3D vis

It could be that the head surface that the Artinis folks used to create the montage was slightly different. This seems to be the case for our EEG standard montages at least. It's tough to know which one is "more correct" I guess...

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],
Expand All @@ -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],
Expand Down
3 changes: 2 additions & 1 deletion tutorials/forward/35_eeg_no_mri.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down