# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
from inspect import signature
import numpy as np
import pytest
from numpy.testing import assert_allclose, assert_array_equal, assert_array_less
from scipy.sparse import csr_array
from scipy.sparse import eye as speye
from scipy.spatial.distance import cdist
import mne
from mne import (
SourceEstimate,
SourceMorph,
VectorSourceEstimate,
VolSourceEstimate,
VolVectorSourceEstimate,
compute_source_morph,
get_volume_labels_from_aseg,
grade_to_vertices,
make_ad_hoc_cov,
make_forward_solution,
make_sphere_model,
read_evokeds,
read_forward_solution,
read_freesurfer_lut,
read_source_estimate,
read_source_morph,
setup_volume_source_space,
)
from mne._freesurfer import _get_atlas_values, _get_mri_info_data
from mne.datasets import testing
from mne.fixes import _get_img_fdata
from mne.minimum_norm import apply_inverse, make_inverse_operator, read_inverse_operator
from mne.source_space._source_space import _add_interpolator, _grid_interp
from mne.transforms import quat_to_rot
from mne.utils import _record_warnings, catch_logging
# Setup paths
data_path = testing.data_path(download=False)
sample_dir = data_path / "MEG" / "sample"
subjects_dir = data_path / "subjects"
fname_evoked = sample_dir / "sample_audvis-ave.fif"
fname_trans = sample_dir / "sample_audvis_trunc-trans.fif"
fname_inv_vol = sample_dir / "sample_audvis_trunc-meg-vol-7-meg-inv.fif"
fname_fwd_vol = sample_dir / "sample_audvis_trunc-meg-vol-7-fwd.fif"
fname_vol_w = sample_dir / "sample_audvis_trunc-grad-vol-7-fwd-sensmap-vol.w"
fname_inv_surf = sample_dir / "sample_audvis_trunc-meg-eeg-oct-6-meg-inv.fif"
fname_aseg = subjects_dir / "sample" / "mri" / "aseg.mgz"
fname_fmorph = data_path / "MEG" / "sample" / "fsaverage_audvis_trunc-meg"
fname_smorph = sample_dir / "sample_audvis_trunc-meg"
fname_t1 = subjects_dir / "sample" / "mri" / "T1.mgz"
fname_vol = subjects_dir / "sample" / "bem" / "sample-volume-7mm-src.fif"
fname_aseg = subjects_dir / "sample" / "mri" / "aseg.mgz"
fname_fs_vol = subjects_dir / "fsaverage" / "bem" / "fsaverage-vol7-nointerp-src.fif.gz"
fname_aseg_fs = subjects_dir / "fsaverage" / "mri" / "aseg.mgz"
fname_stc = sample_dir / "fsaverage_audvis_trunc-meg"
pytest.importorskip("nibabel")
def _real_vec_stc():
inv = read_inverse_operator(fname_inv_surf)
evoked = read_evokeds(fname_evoked, baseline=(None, 0))[0].crop(0, 0.01)
return apply_inverse(evoked, inv, pick_ori="vector")
def test_sourcemorph_consistency():
"""Test SourceMorph class consistency."""
assert (
list(signature(SourceMorph.__init__).parameters)[1:-1]
== mne.morph._SOURCE_MORPH_ATTRIBUTES
)
@testing.requires_testing_data
def test_sparse_morph():
"""Test sparse morphing."""
rng = np.random.RandomState(0)
vertices_fs = [
np.sort(rng.permutation(np.arange(10242))[:4]),
np.sort(rng.permutation(np.arange(10242))[:6]),
]
data = rng.randn(10, 1)
stc_fs = SourceEstimate(data, vertices_fs, 1, 1, "fsaverage")
spheres_fs = [
mne.read_surface(subjects_dir / "fsaverage" / "surf" / f"{hemi}.sphere.reg")[0]
for hemi in ("lh", "rh")
]
spheres_sample = [
mne.read_surface(subjects_dir / "sample" / "surf" / f"{hemi}.sphere.reg")[0]
for hemi in ("lh", "rh")
]
morph_fs_sample = compute_source_morph(
stc_fs,
"fsaverage",
"sample",
sparse=True,
spacing=None,
subjects_dir=subjects_dir,
)
stc_sample = morph_fs_sample.apply(stc_fs)
offset = 0
orders = list()
for v1, s1, v2, s2 in zip(
stc_fs.vertices, spheres_fs, stc_sample.vertices, spheres_sample
):
dists = cdist(s1[v1], s2[v2])
order = np.argmin(dists, axis=-1)
assert_array_less(dists[np.arange(len(order)), order], 1.5) # mm
orders.append(order + offset)
offset += len(order)
assert_allclose(stc_fs.data, stc_sample.data[np.concatenate(orders)])
# Return
morph_sample_fs = compute_source_morph(
stc_sample,
"sample",
"fsaverage",
sparse=True,
spacing=None,
subjects_dir=subjects_dir,
)
stc_fs_return = morph_sample_fs.apply(stc_sample)
offset = 0
orders = list()
for v1, s, v2 in zip(stc_fs.vertices, spheres_fs, stc_fs_return.vertices):
dists = cdist(s[v1], s[v2])
order = np.argmin(dists, axis=-1)
assert_array_less(dists[np.arange(len(order)), order], 1.5) # mm
orders.append(order + offset)
offset += len(order)
assert_allclose(stc_fs.data, stc_fs_return.data[np.concatenate(orders)])
@testing.requires_testing_data
def test_xhemi_morph():
"""Test cross-hemisphere morphing."""
stc = read_source_estimate(fname_stc, subject="sample")
# smooth 1 for speed where possible
smooth = 4
spacing = 4
n_grade_verts = 2562
stc = compute_source_morph(
stc,
"sample",
"fsaverage_sym",
smooth=smooth,
warn=False,
spacing=spacing,
subjects_dir=subjects_dir,
).apply(stc)
morph = compute_source_morph(
stc,
"fsaverage_sym",
"fsaverage_sym",
smooth=1,
xhemi=True,
warn=False,
spacing=[stc.vertices[0], []],
subjects_dir=subjects_dir,
)
stc_xhemi = morph.apply(stc)
assert stc_xhemi.data.shape[0] == n_grade_verts
assert stc_xhemi.rh_data.shape[0] == 0
assert len(stc_xhemi.vertices[1]) == 0
assert stc_xhemi.lh_data.shape[0] == n_grade_verts
assert len(stc_xhemi.vertices[0]) == n_grade_verts
# complete reversal mapping
morph = compute_source_morph(
stc,
"fsaverage_sym",
"fsaverage_sym",
smooth=smooth,
xhemi=True,
warn=False,
spacing=stc.vertices,
subjects_dir=subjects_dir,
)
mm = morph.morph_mat
assert mm.shape == (n_grade_verts * 2,) * 2
assert mm.size > n_grade_verts * 2
assert mm[:n_grade_verts, :n_grade_verts].size == 0 # L to L
assert mm[n_grade_verts:, n_grade_verts:].size == 0 # R to L
assert mm[n_grade_verts:, :n_grade_verts].size > n_grade_verts # L to R
assert mm[:n_grade_verts, n_grade_verts:].size > n_grade_verts # R to L
# more complicated reversal mapping
vertices_use = [stc.vertices[0], np.arange(10242)]
n_src_verts = len(vertices_use[1])
assert vertices_use[0].shape == (n_grade_verts,)
assert vertices_use[1].shape == (n_src_verts,)
# ensure it's sufficiently diffirent to manifest round-trip errors
assert np.isin(vertices_use[1], stc.vertices[1]).mean() < 0.3
morph = compute_source_morph(
stc,
"fsaverage_sym",
"fsaverage_sym",
smooth=smooth,
xhemi=True,
warn=False,
spacing=vertices_use,
subjects_dir=subjects_dir,
)
mm = morph.morph_mat
assert mm.shape == (n_grade_verts + n_src_verts, n_grade_verts * 2)
assert mm[:n_grade_verts, :n_grade_verts].size == 0
assert mm[n_grade_verts:, n_grade_verts:].size == 0
assert mm[:n_grade_verts, n_grade_verts:].size > n_grade_verts
assert mm[n_grade_verts:, :n_grade_verts].size > n_src_verts
# morph forward then back
stc_xhemi = morph.apply(stc)
morph = compute_source_morph(
stc_xhemi,
"fsaverage_sym",
"fsaverage_sym",
smooth=smooth,
xhemi=True,
warn=False,
spacing=stc.vertices,
subjects_dir=subjects_dir,
)
stc_return = morph.apply(stc_xhemi)
for hi in range(2):
assert_array_equal(stc_return.vertices[hi], stc.vertices[hi])
correlation = np.corrcoef(stc.data.ravel(), stc_return.data.ravel())[0, 1]
assert correlation > 0.9 # not great b/c of sparse grade + small smooth
@testing.requires_testing_data
@pytest.mark.parametrize(
"smooth, lower, upper, n_warn, dtype",
[
(None, 0.959, 0.963, 0, float),
(3, 0.968, 0.971, 2, complex),
("nearest", 0.98, 0.99, 0, float),
],
)
def test_surface_source_morph_round_trip(smooth, lower, upper, n_warn, dtype):
"""Test round-trip morphing yields similar STCs."""
kwargs = dict(smooth=smooth, warn=True, subjects_dir=subjects_dir)
stc = mne.read_source_estimate(fname_smorph)
if dtype is complex:
stc.data = 1j * stc.data
assert_array_equal(stc.data.real, 0.0)
with _record_warnings() as w:
morph = compute_source_morph(stc, "sample", "fsaverage", **kwargs)
w = [ww for ww in w if "vertices not included" in str(ww.message)]
assert len(w) == n_warn
assert morph.morph_mat.shape == (20484, len(stc.data))
stc_fs = morph.apply(stc)
morph_back = compute_source_morph(
stc_fs, "fsaverage", "sample", spacing=stc.vertices, **kwargs
)
assert morph_back.morph_mat.shape == (len(stc.data), 20484)
stc_back = morph_back.apply(stc_fs)
corr = np.corrcoef(stc.data.ravel(), stc_back.data.ravel())[0, 1]
assert lower <= corr <= upper
# check the round-trip power
assert_power_preserved(stc, stc_back)
@testing.requires_testing_data
def test_surface_source_morph_shortcut():
"""Test that our shortcut for smooth=0 works."""
stc = mne.read_source_estimate(fname_smorph)
morph_identity = compute_source_morph(
stc,
"sample",
"sample",
spacing=stc.vertices,
smooth=0,
subjects_dir=subjects_dir,
)
stc_back = morph_identity.apply(stc)
assert_allclose(stc_back.data, stc.data, rtol=1e-4)
abs_sum = morph_identity.morph_mat - speye(len(stc.data), format="csc")
abs_sum = np.abs(abs_sum.data).sum()
assert abs_sum < 1e-4
def assert_power_preserved(orig, new, limits=(1.0, 1.05)):
"""Assert that the power is preserved during a round-trip morph."""
__tracebackhide__ = True
for kind in ("real", "imag"):
numer = np.linalg.norm(getattr(orig.data, kind))
denom = np.linalg.norm(getattr(new.data, kind))
if numer == denom == 0.0: # no data of this type
continue
power_ratio = numer / denom
min_, max_ = limits
assert min_ < power_ratio < max_, f"Power ratio {kind} = {power_ratio}"
@testing.requires_testing_data
def test_surface_vector_source_morph(tmp_path):
"""Test surface and vector source estimate morph."""
pytest.importorskip("h5io")
inverse_operator_surf = read_inverse_operator(fname_inv_surf)
stc_surf = read_source_estimate(fname_smorph, subject="sample")
stc_surf.crop(0.09, 0.1) # for faster computation
stc_vec = _real_vec_stc()
source_morph_surf = compute_source_morph(
inverse_operator_surf["src"], subjects_dir=subjects_dir, smooth=1, warn=False
) # smooth 1 for speed
assert source_morph_surf.subject_from == "sample"
assert source_morph_surf.subject_to == "fsaverage"
assert source_morph_surf.kind == "surface"
assert isinstance(source_morph_surf.src_data, dict)
assert isinstance(source_morph_surf.src_data["vertices_from"], list)
assert isinstance(source_morph_surf, SourceMorph)
stc_surf_morphed = source_morph_surf.apply(stc_surf)
assert isinstance(stc_surf_morphed, SourceEstimate)
stc_vec_morphed = source_morph_surf.apply(stc_vec)
with pytest.raises(ValueError, match="Invalid value for the 'output'"):
source_morph_surf.apply(stc_surf, output="nifti1")
# check if correct class after morphing
assert isinstance(stc_surf_morphed, SourceEstimate)
assert isinstance(stc_vec_morphed, VectorSourceEstimate)
# check __repr__
assert "surface" in repr(source_morph_surf)
# check loading and saving for surf
source_morph_surf.save(tmp_path / "42.h5")
source_morph_surf_r = read_source_morph(tmp_path / "42.h5")
assert all(
[
read == saved
for read, saved in zip(
sorted(source_morph_surf_r.__dict__), sorted(source_morph_surf.__dict__)
)
]
)
# check wrong subject correction
stc_surf.subject = None
assert isinstance(source_morph_surf.apply(stc_surf), SourceEstimate)
# degenerate
stc_vol = read_source_estimate(fname_vol_w, "sample")
with pytest.raises(TypeError, match="stc_from must be an instance"):
source_morph_surf.apply(stc_vol)
@pytest.mark.slowtest
@testing.requires_testing_data
def test_volume_source_morph_basic(tmp_path):
"""Test volume source estimate morph, special cases and exceptions."""
nib = pytest.importorskip("nibabel")
pytest.importorskip("h5io")
pytest.importorskip("dipy")
inverse_operator_vol = read_inverse_operator(fname_inv_vol)
stc_vol = read_source_estimate(fname_vol_w, "sample")
# check for invalid input type
with pytest.raises(TypeError, match="src must be"):
compute_source_morph(src=42)
# check for raising an error if neither
# inverse_operator_vol['src'][0]['subject_his_id'] nor subject_from is set,
# but attempting to perform a volume morph
src = inverse_operator_vol["src"]
assert src._subject is None # already None on disk (old!)
with pytest.raises(ValueError, match="subject_from could not be inferred"):
with pytest.warns(RuntimeWarning, match="recommend regenerating"):
compute_source_morph(src=src, subjects_dir=subjects_dir)
# check infer subject_from from src[0]['subject_his_id']
src[0]["subject_his_id"] = "sample"
with pytest.raises(ValueError, match="Inter-hemispheric morphing"):
compute_source_morph(src=src, subjects_dir=subjects_dir, xhemi=True)
with pytest.raises(ValueError, match="Only surface.*sparse morph"):
compute_source_morph(src=src, sparse=True, subjects_dir=subjects_dir)
# terrible quality but fast
zooms = 20
kwargs = dict(zooms=zooms, niter_sdr=(1,), niter_affine=(1,))
source_morph_vol = compute_source_morph(
subjects_dir=subjects_dir, src=fname_inv_vol, subject_from="sample", **kwargs
)
shape = (13,) * 3 # for the given zooms
assert source_morph_vol.subject_from == "sample"
# the brain used in sample data has shape (255, 255, 255)
assert tuple(source_morph_vol.sdr_morph.domain_shape) == shape
assert tuple(source_morph_vol.pre_affine.domain_shape) == shape
# proofs the above
assert_array_equal(source_morph_vol.zooms, (zooms,) * 3)
# assure proper src shape
mri_size = (src[0]["mri_height"], src[0]["mri_depth"], src[0]["mri_width"])
assert source_morph_vol.src_data["src_shape_full"] == mri_size
fwd = read_forward_solution(fname_fwd_vol)
fwd["src"][0]["subject_his_id"] = "sample" # avoid further warnings
source_morph_vol = compute_source_morph(
fwd["src"], "sample", "sample", subjects_dir=subjects_dir, **kwargs
)
# check wrong subject_to
with pytest.raises(OSError, match="cannot read file"):
compute_source_morph(fwd["src"], "sample", "42", subjects_dir=subjects_dir)
# two different ways of saving
source_morph_vol.save(tmp_path / "vol")
# check loading
source_morph_vol_r = read_source_morph(tmp_path / "vol-morph.h5")
# check for invalid file name handling ()
with pytest.raises(OSError, match="not found"):
read_source_morph(tmp_path / "42")
# check morph
stc_vol_morphed = source_morph_vol.apply(stc_vol)
# old way, verts do not match
assert not np.array_equal(stc_vol_morphed.vertices[0], stc_vol.vertices[0])
# vector
stc_vol_vec = VolVectorSourceEstimate(
np.tile(stc_vol.data[:, np.newaxis], (1, 3, 1)), stc_vol.vertices, 0, 1
)
stc_vol_vec_morphed = source_morph_vol.apply(stc_vol_vec)
assert isinstance(stc_vol_vec_morphed, VolVectorSourceEstimate)
for ii in range(3):
assert_allclose(stc_vol_vec_morphed.data[:, ii], stc_vol_morphed.data)
# check output as NIfTI
assert isinstance(
source_morph_vol.apply(stc_vol_vec, output="nifti2"), nib.Nifti2Image
)
# check for subject_from mismatch
source_morph_vol_r.subject_from = "42"
with pytest.raises(ValueError, match="subject_from must match"):
source_morph_vol_r.apply(stc_vol_morphed)
# check if nifti is in grid morph space with voxel_size == spacing
img_morph_res = source_morph_vol.apply(stc_vol, output="nifti1")
# assure morph spacing
assert isinstance(img_morph_res, nib.Nifti1Image)
assert img_morph_res.header.get_zooms()[:3] == (zooms,) * 3
# assure src shape
img_mri_res = source_morph_vol.apply(stc_vol, output="nifti1", mri_resolution=True)
assert isinstance(img_mri_res, nib.Nifti1Image)
assert img_mri_res.shape == (
src[0]["mri_height"],
src[0]["mri_depth"],
src[0]["mri_width"],
) + (img_mri_res.shape[3],)
# check if nifti is defined resolution with voxel_size == (5., 5., 5.)
img_any_res = source_morph_vol.apply(
stc_vol, output="nifti1", mri_resolution=(5.0, 5.0, 5.0)
)
assert isinstance(img_any_res, nib.Nifti1Image)
assert img_any_res.header.get_zooms()[:3] == (5.0, 5.0, 5.0)
# check if morph outputs correct data
assert isinstance(stc_vol_morphed, VolSourceEstimate)
# check if loaded and saved objects contain the same
assert all(
[
read == saved
for read, saved in zip(
sorted(source_morph_vol_r.__dict__), sorted(source_morph_vol.__dict__)
)
]
)
# check __repr__
assert "volume" in repr(source_morph_vol)
# check Nifti2Image
assert isinstance(
source_morph_vol.apply(
stc_vol, mri_resolution=True, mri_space=True, output="nifti2"
),
nib.Nifti2Image,
)
# Degenerate conditions
with pytest.raises(TypeError, match="output must be"):
source_morph_vol.apply(stc_vol, output=1)
with pytest.raises(ValueError, match="subject_from does not match"):
compute_source_morph(src=src, subject_from="42")
with pytest.raises(ValueError, match="output"):
source_morph_vol.apply(stc_vol, output="42")
with pytest.raises(ValueError, match="subject_to cannot be None"):
compute_source_morph(src, "sample", None, subjects_dir=subjects_dir)
# Check if not morphed, but voxel size not boolean, raise ValueError.
# Note that this check requires dipy to not raise the dipy ImportError
# before checking if the actual voxel size error will raise.
with pytest.raises(ValueError, match="Cannot infer original voxel size"):
stc_vol.as_volume(inverse_operator_vol["src"], mri_resolution=4)
stc_surf = read_source_estimate(fname_stc, "sample")
with pytest.raises(TypeError, match="stc_from must be an instance"):
source_morph_vol.apply(stc_surf)
# src_to
source_morph_vol = compute_source_morph(
fwd["src"],
subject_from="sample",
src_to=fwd["src"],
subject_to="sample",
subjects_dir=subjects_dir,
**kwargs,
)
stc_vol_2 = source_morph_vol.apply(stc_vol)
# new way, verts match
assert_array_equal(stc_vol.vertices[0], stc_vol_2.vertices[0])
stc_vol_bad = VolSourceEstimate(
stc_vol.data[:-1], [stc_vol.vertices[0][:-1]], stc_vol.tmin, stc_vol.tstep
)
match = (
"vertices do not match between morph \\(4157\\) and stc \\(4156\\).*"
"\n.*\n.*\n.*Vertices were likely excluded during forward computatio.*"
)
with pytest.raises(ValueError, match=match):
source_morph_vol.apply(stc_vol_bad)
# nifti outputs and stc equiv
img_vol = source_morph_vol.apply(stc_vol, output="nifti1")
img_vol_2 = stc_vol_2.as_volume(src=fwd["src"], mri_resolution=False)
assert_allclose(img_vol.affine, img_vol_2.affine)
img_vol = img_vol.get_fdata()
img_vol_2 = img_vol_2.get_fdata()
assert img_vol.shape == img_vol_2.shape
assert_allclose(img_vol, img_vol_2)
@pytest.mark.slowtest
@testing.requires_testing_data
@pytest.mark.parametrize(
"subject_from, subject_to, lower, upper, dtype, morph_mat",
[
("sample", "fsaverage", 5.9, 6.1, float, False),
("fsaverage", "fsaverage", 0.0, 0.1, float, False),
("sample", "sample", 0.0, 0.1, complex, False),
("sample", "sample", 0.0, 0.1, float, True), # morph_mat
("sample", "fsaverage", 10, 12, float, True), # morph_mat
],
)
def test_volume_source_morph_round_trip(
tmp_path, subject_from, subject_to, lower, upper, dtype, morph_mat, monkeypatch
):
"""Test volume source estimate morph round-trips well."""
nib = pytest.importorskip("nibabel")
pytest.importorskip("h5io")
pytest.importorskip("dipy")
from nibabel.processing import resample_from_to
src = dict()
if morph_mat:
# ~1.5 minutes with pos=7. (4157 morphs!) for sample, so only test
# morph_mat computation mode with a few labels
label_names = sorted(get_volume_labels_from_aseg(fname_aseg))[1:2]
if "sample" in (subject_from, subject_to):
src["sample"] = setup_volume_source_space(
"sample",
subjects_dir=subjects_dir,
volume_label=label_names,
mri=fname_aseg,
)
assert sum(s["nuse"] for s in src["sample"]) == 12
if "fsaverage" in (subject_from, subject_to):
src["fsaverage"] = setup_volume_source_space(
"fsaverage",
subjects_dir=subjects_dir,
volume_label=label_names[:3],
mri=fname_aseg_fs,
)
assert sum(s["nuse"] for s in src["fsaverage"]) == 16
else:
assert not morph_mat
if "sample" in (subject_from, subject_to):
src["sample"] = mne.read_source_spaces(fname_vol)
src["sample"][0]["subject_his_id"] = "sample"
assert src["sample"][0]["nuse"] == 4157
if "fsaverage" in (subject_from, subject_to):
# Created to save space with:
#
# bem = op.join(op.dirname(mne.__file__), 'data', 'fsaverage',
# 'fsaverage-inner_skull-bem.fif')
# src_fsaverage = mne.setup_volume_source_space(
# 'fsaverage', pos=7., bem=bem, mindist=0,
# subjects_dir=subjects_dir, add_interpolator=False)
# mne.write_source_spaces(fname_fs_vol, src_fsaverage,
# overwrite=True)
#
# For speed we do it without the interpolator because it's huge.
src["fsaverage"] = mne.read_source_spaces(fname_fs_vol)
src["fsaverage"][0].update(
vol_dims=np.array([23, 29, 25]), seg_name="brain"
)
_add_interpolator(src["fsaverage"])
assert src["fsaverage"][0]["nuse"] == 6379
src_to, src_from = src[subject_to], src[subject_from]
del src
# No SDR just for speed once everything works
kwargs = dict(
niter_sdr=(), niter_affine=(1,), subjects_dir=subjects_dir, verbose=True
)
morph_from_to = compute_source_morph(
src=src_from, src_to=src_to, subject_to=subject_to, **kwargs
)
morph_to_from = compute_source_morph(
src=src_to, src_to=src_from, subject_to=subject_from, **kwargs
)
nuse = sum(s["nuse"] for s in src_from)
assert nuse > 10
use = np.linspace(0, nuse - 1, 10).round().astype(int)
data = np.eye(nuse)[:, use]
if dtype is complex:
data = data * 1j
vertices = [s["vertno"] for s in src_from]
stc_from = VolSourceEstimate(data, vertices, 0, 1)
with catch_logging() as log:
stc_from_rt = morph_to_from.apply(
morph_from_to.apply(stc_from, verbose="debug")
)
log = log.getvalue()
assert "individual volume morph" in log
maxs = np.argmax(stc_from_rt.data, axis=0)
src_rr = np.concatenate([s["rr"][s["vertno"]] for s in src_from])
dists = 1000 * np.linalg.norm(src_rr[use] - src_rr[maxs], axis=1)
mu = np.mean(dists)
# fsaverage=5.99; 7.97 without additional src_ras_t fix
# fsaverage=7.97; 25.4 without src_ras_t fix
assert lower <= mu < upper, f"round-trip distance {mu}"
# check that pre_affine is close to identity when subject_to==subject_from
if subject_to == subject_from:
for morph in (morph_to_from, morph_from_to):
assert_allclose(morph.pre_affine.affine, np.eye(4), atol=1e-2)
# check that power is more or less preserved (labelizing messes with this)
if morph_mat:
if subject_to == "fsaverage":
limits = (18, 18.5)
else:
limits = (7, 7.5)
else:
limits = (1, 1.2)
stc_from_unit = stc_from.copy().crop(0, 0)
stc_from_unit._data.fill(1.0)
stc_from_unit_rt = morph_to_from.apply(morph_from_to.apply(stc_from_unit))
assert_power_preserved(stc_from_unit, stc_from_unit_rt, limits=limits)
if morph_mat:
fname = tmp_path / "temp-morph.h5"
morph_to_from.save(fname)
morph_to_from = read_source_morph(fname)
assert morph_to_from.vol_morph_mat is None
morph_to_from.compute_vol_morph_mat(verbose=True)
morph_to_from.save(fname, overwrite=True)
morph_to_from = read_source_morph(fname)
assert isinstance(morph_to_from.vol_morph_mat, csr_array), "csr"
# equivalence (plus automatic calling)
assert morph_from_to.vol_morph_mat is None
monkeypatch.setattr(mne.morph, "_VOL_MAT_CHECK_RATIO", 0.0)
with catch_logging() as log:
with pytest.warns(RuntimeWarning, match=r"calling morph\.compute"):
stc_from_rt_lin = morph_to_from.apply(
morph_from_to.apply(stc_from, verbose="debug")
)
assert isinstance(morph_from_to.vol_morph_mat, csr_array), "csr"
log = log.getvalue()
assert "sparse volume morph matrix" in log
assert_allclose(stc_from_rt.data, stc_from_rt_lin.data)
del stc_from_rt_lin
stc_from_unit_rt_lin = morph_to_from.apply(morph_from_to.apply(stc_from_unit))
assert_allclose(stc_from_unit_rt.data, stc_from_unit_rt_lin.data)
del stc_from_unit_rt_lin
del stc_from, stc_from_rt
# before and after morph, check the proportion of vertices
# that are inside and outside the brainmask.mgz
brain = nib.load(subjects_dir / subject_from / "mri" / "brain.mgz")
mask = _get_img_fdata(brain) > 0
if subject_from == subject_to == "sample":
for stc in [stc_from_unit, stc_from_unit_rt]:
img = stc.as_volume(src_from, mri_resolution=True)
img = nib.Nifti1Image( # abs to convert complex
np.abs(_get_img_fdata(img)[:, :, :, 0]), img.affine
)
img = _get_img_fdata(resample_from_to(img, brain, order=1))
assert img.shape == mask.shape
in_ = img[mask].astype(bool).mean()
out = img[~mask].astype(bool).mean()
if morph_mat:
out_max = 0.001
in_min, in_max = 0.005, 0.007
else:
out_max = 0.02
in_min, in_max = 0.97, 0.98
assert out < out_max, f"proportion out of volume {out}"
assert in_min < in_ < in_max, f"proportion inside volume {in_}"
@pytest.mark.slowtest
@testing.requires_testing_data
def test_morph_stc_dense():
"""Test morphing stc."""
subject_from = "sample"
subject_to = "fsaverage"
stc_from = read_source_estimate(fname_smorph, subject="sample")
stc_to = read_source_estimate(fname_fmorph)
# make sure we can specify grade
stc_from.crop(0.09, 0.1) # for faster computation
stc_to.crop(0.09, 0.1) # for faster computation
assert_array_equal(
stc_to.time_as_index([0.09, 0.1], use_rounding=True), [0, len(stc_to.times) - 1]
)
# After dep change this to:
morph = compute_source_morph(
subject_to=subject_to,
spacing=3,
smooth=12,
src=stc_from,
subjects_dir=subjects_dir,
precompute=True,
)
assert morph.vol_morph_mat is None # a no-op for surface
stc_to1 = morph.apply(stc_from)
assert_allclose(stc_to.data, stc_to1.data, atol=1e-5)
mean_from = stc_from.data.mean(axis=0)
mean_to = stc_to1.data.mean(axis=0)
assert np.corrcoef(mean_to, mean_from).min() > 0.999
vertices_to = grade_to_vertices(subject_to, grade=3, subjects_dir=subjects_dir)
# make sure we can fill by morphing
with pytest.warns(RuntimeWarning, match="consider increasing"):
morph = compute_source_morph(
stc_from,
subject_from,
subject_to,
spacing=None,
smooth=1,
subjects_dir=subjects_dir,
)
stc_to5 = morph.apply(stc_from)
assert stc_to5.data.shape[0] == 163842 + 163842
# Morph vector data
stc_vec = _real_vec_stc()
stc_vec_to1 = compute_source_morph(
stc_vec,
subject_from,
subject_to,
subjects_dir=subjects_dir,
spacing=vertices_to,
smooth=1,
warn=False,
).apply(stc_vec)
assert stc_vec_to1.subject == subject_to
assert stc_vec_to1.tmin == stc_vec.tmin
assert stc_vec_to1.tstep == stc_vec.tstep
assert len(stc_vec_to1.lh_vertno) == 642
assert len(stc_vec_to1.rh_vertno) == 642
# Degenerate conditions
# Morphing to a density that is too high should raise an informative error
# (here we need to push to grade=6, but for some subjects even grade=5
# will break)
with pytest.raises(ValueError, match="Cannot use icosahedral grade 6 "):
compute_source_morph(
stc_to1,
subject_from=subject_to,
subject_to=subject_from,
spacing=6,
subjects_dir=subjects_dir,
)
del stc_to1
with pytest.raises(ValueError, match="smooth.* has to be at least 0"):
compute_source_morph(
stc_from,
subject_from,
subject_to,
spacing=5,
smooth=-1,
subjects_dir=subjects_dir,
)
# subject from mismatch
with pytest.raises(ValueError, match="subject_from does not match"):
compute_source_morph(stc_from, subject_from="foo", subjects_dir=subjects_dir)
# only one set of vertices
with pytest.raises(ValueError, match="grade.*list must have two elements"):
compute_source_morph(
stc_from,
subject_from=subject_from,
spacing=[vertices_to[0]],
subjects_dir=subjects_dir,
)
@testing.requires_testing_data
def test_morph_stc_sparse():
"""Test morphing stc with sparse=True."""
subject_from = "sample"
subject_to = "fsaverage"
# Morph sparse data
# Make a sparse stc
stc_from = read_source_estimate(fname_smorph, subject="sample")
stc_from.vertices[0] = stc_from.vertices[0][[100, 500]]
stc_from.vertices[1] = stc_from.vertices[1][[200]]
stc_from._data = stc_from._data[:3]
stc_to_sparse = compute_source_morph(
stc_from,
subject_from=subject_from,
subject_to=subject_to,
spacing=None,
sparse=True,
subjects_dir=subjects_dir,
).apply(stc_from)
assert_allclose(
np.sort(stc_from.data.sum(axis=1)), np.sort(stc_to_sparse.data.sum(axis=1))
)
assert len(stc_from.rh_vertno) == len(stc_to_sparse.rh_vertno)
assert len(stc_from.lh_vertno) == len(stc_to_sparse.lh_vertno)
assert stc_to_sparse.subject == subject_to
assert stc_from.tmin == stc_from.tmin
assert stc_from.tstep == stc_from.tstep
stc_from.vertices[0] = np.array([], dtype=np.int64)
stc_from._data = stc_from._data[:1]
stc_to_sparse = compute_source_morph(
stc_from,
subject_from,
subject_to,
spacing=None,
sparse=True,
subjects_dir=subjects_dir,
).apply(stc_from)
assert_allclose(
np.sort(stc_from.data.sum(axis=1)), np.sort(stc_to_sparse.data.sum(axis=1))
)
assert len(stc_from.rh_vertno) == len(stc_to_sparse.rh_vertno)
assert len(stc_from.lh_vertno) == len(stc_to_sparse.lh_vertno)
assert stc_to_sparse.subject == subject_to
assert stc_from.tmin == stc_from.tmin
assert stc_from.tstep == stc_from.tstep
# Degenerate cases
with pytest.raises(ValueError, match="spacing must be set to None"):
compute_source_morph(
stc_from,
subject_from=subject_from,
subject_to=subject_to,
spacing=5,
sparse=True,
subjects_dir=subjects_dir,
)
with pytest.raises(ValueError, match="xhemi=True can only be used with"):
compute_source_morph(
stc_from,
subject_from=subject_from,
subject_to=subject_to,
spacing=None,
sparse=True,
xhemi=True,
subjects_dir=subjects_dir,
)
@testing.requires_testing_data
@pytest.mark.parametrize(
"sl, n_real, n_mri, n_orig",
[
# First and last should add up, middle can have overlap should be <= sum
(slice(0, 1), 37, 138, 8),
(slice(1, 2), 51, 204, 12),
(slice(0, 2), 88, 324, 20),
],
)
def test_volume_labels_morph(tmp_path, sl, n_real, n_mri, n_orig):
"""Test generating a source space from volume label."""
nib = pytest.importorskip("nibabel")
n_use = (sl.stop - sl.start) // (sl.step or 1)
# see gh-5224
evoked = mne.read_evokeds(fname_evoked)[0].crop(0, 0)
evoked.pick(evoked.ch_names[:306:8])
evoked.info.normalize_proj()
n_ch = len(evoked.ch_names)
lut, _ = read_freesurfer_lut()
label_names = sorted(get_volume_labels_from_aseg(fname_aseg))
use_label_names = label_names[sl]
src = setup_volume_source_space(
"sample",
subjects_dir=subjects_dir,
volume_label=use_label_names,
mri=fname_aseg,
)
assert len(src) == n_use
assert src.kind == "volume"
n_src = sum(s["nuse"] for s in src)
sphere = make_sphere_model("auto", "auto", evoked.info)
fwd = make_forward_solution(evoked.info, fname_trans, src, sphere)
assert fwd["sol"]["data"].shape == (n_ch, n_src * 3)
inv = make_inverse_operator(
evoked.info, fwd, make_ad_hoc_cov(evoked.info), loose=1.0
)
stc = apply_inverse(evoked, inv)
assert stc.data.shape == (n_src, 1)
img = stc.as_volume(src, mri_resolution=True)
assert img.shape == (86, 86, 86, 1)
n_on = np.array(img.dataobj).astype(bool).sum()
aseg_img = _get_img_fdata(nib.load(fname_aseg))
n_got_real = np.isin(
aseg_img.ravel(), [lut[name] for name in use_label_names]
).sum()
assert n_got_real == n_real
# - This was 291 on `main` before gh-5590
# - Refactoring transforms it became 279 with a < 1e-8 change in vox_mri_t
# - Dropped to 123 once nearest-voxel was used in gh-7653
# - Jumped back up to 330 with morphing fixes actually correctly
# interpolating across all volumes
assert aseg_img.shape == img.shape[:3]
assert n_on == n_mri
for ii in range(2):
# should work with (ii=0) or without (ii=1) the interpolator
if ii:
src[0]["interpolator"] = None
img = stc.as_volume(src, mri_resolution=False)
n_on = np.array(img.dataobj).astype(bool).sum()
# was 20 on `main` before gh-5590
# then 44 before gh-7653, which took it back to 20
assert n_on == n_orig
# without the interpolator, this should fail
assert src[0]["interpolator"] is None
with pytest.raises(RuntimeError, match=r".*src\[0\], .* mri_resolution"):
stc.as_volume(src, mri_resolution=True)
@pytest.fixture(scope="session", params=[testing._pytest_param()])
def _mixed_morph_srcs():
pytest.importorskip("nibabel")
pytest.importorskip("dipy")
# create a mixed source space
labels_vol = ["Left-Cerebellum-Cortex", "Right-Cerebellum-Cortex"]
src = mne.setup_source_space(
"sample", spacing="oct3", add_dist=False, subjects_dir=subjects_dir
)
src += mne.setup_volume_source_space(
"sample",
mri=fname_aseg,
pos=10.0,
volume_label=labels_vol,
subjects_dir=subjects_dir,
add_interpolator=True,
verbose=True,
)
# create the destination space
src_fs = mne.read_source_spaces(
subjects_dir / "fsaverage" / "bem" / "fsaverage-ico-5-src.fif"
)
src_fs += mne.setup_volume_source_space(
"fsaverage",
pos=7.0,
volume_label=labels_vol,
subjects_dir=subjects_dir,
add_interpolator=False,
verbose=True,
)
del labels_vol
with pytest.raises(ValueError, match="src_to must be provided .* mixed"):
mne.compute_source_morph(
src=src,
subject_from="sample",
subject_to="fsaverage",
subjects_dir=subjects_dir,
)
with pytest.warns(RuntimeWarning, match="not included in smoothing"):
morph = mne.compute_source_morph(
src=src,
subject_from="sample",
subject_to="fsaverage",
subjects_dir=subjects_dir,
niter_affine=[1, 0, 0],
niter_sdr=[1, 0, 0],
src_to=src_fs,
smooth=5,
verbose=True,
)
return morph, src, src_fs
@pytest.mark.slowtest
@pytest.mark.parametrize("vector", (False, True))
def test_mixed_source_morph(_mixed_morph_srcs, vector):
"""Test mixed source space morphing."""
nib = pytest.importorskip("nibabel")
pytest.importorskip("dipy")
morph, src, src_fs = _mixed_morph_srcs
# Test some basic properties in the subject's own space
lut, _ = read_freesurfer_lut()
ids = [lut[s["seg_name"]] for s in src[2:]]
del lut
vertices = [s["vertno"] for s in src]
n_vertices = sum(len(v) for v in vertices)
data = np.zeros((n_vertices, 3, 1))
data[:, 1] = 1.0
klass = mne.MixedVectorSourceEstimate
if not vector:
data = data[:, 1]
klass = klass._scalar_class
stc = klass(data, vertices, 0, 1, "sample")
vol_info = _get_mri_info_data(fname_aseg, data=True)
rrs = np.concatenate([src[2]["rr"][sp["vertno"]] for sp in src[2:]])
n_want = np.isin(_get_atlas_values(vol_info, rrs), ids).sum()
img = _get_img_fdata(stc.volume().as_volume(src, mri_resolution=False))
assert img.astype(bool).sum() == n_want
img_res = nib.load(fname_aseg)
n_want = np.isin(_get_img_fdata(img_res), ids).sum()
img = _get_img_fdata(stc.volume().as_volume(src, mri_resolution=True))
assert img.astype(bool).sum() > n_want # way more get interpolated into
with pytest.raises(TypeError, match="stc_from must be an instance"):
morph.apply(1.0)
# Now actually morph
stc_fs = morph.apply(stc)
vol_info = _get_mri_info_data(fname_aseg_fs, data=True)
rrs = np.concatenate([src_fs[2]["rr"][sp["vertno"]] for sp in src_fs[2:]])
n_want = np.isin(_get_atlas_values(vol_info, rrs), ids).sum()
with pytest.raises(ValueError, match=r"stc\.subject does not match src s"):
stc_fs.volume().as_volume(src, mri_resolution=False)
img = _get_img_fdata(stc_fs.volume().as_volume(src_fs, mri_resolution=False))
assert img.astype(bool).sum() == n_want # correct number of voxels
# Morph separate parts and compare to morphing the entire one
stc_fs_surf = morph.apply(stc.surface())
stc_fs_vol = morph.apply(stc.volume())
stc_fs_2 = stc_fs.__class__(
np.concatenate([stc_fs_surf.data, stc_fs_vol.data]),
stc_fs_surf.vertices + stc_fs_vol.vertices,
stc_fs.tmin,
stc_fs.tstep,
stc_fs.subject,
)
assert_allclose(stc_fs.data, stc_fs_2.data)
def _rand_affine(rng):
quat = rng.randn(3)
quat /= 5 * np.linalg.norm(quat)
affine = np.eye(4)
affine[:3, 3] = rng.randn(3) / 5.0
affine[:3, :3] = quat_to_rot(quat)
return affine
_shapes = (
(10, 10, 10),
(20, 5, 10),
(5, 10, 20),
)
_affines = (
[[2, 0, 0, 1], [0, 0, 1, -1], [0, -1, 0, 2], [0, 0, 0, 1]],
np.eye(4),
np.eye(4)[[0, 2, 1, 3]],
"rand",
)
@pytest.mark.parametrize("from_shape", _shapes)
@pytest.mark.parametrize("from_affine", _affines)
@pytest.mark.parametrize("to_shape", _shapes)
@pytest.mark.parametrize("to_affine", _affines)
@pytest.mark.parametrize("order", [0, 1])
@pytest.mark.parametrize("seed", [0, 1])
def test_resample_equiv(from_shape, from_affine, to_shape, to_affine, order, seed):
"""Test resampling equivalences."""
pytest.importorskip("nibabel")
pytest.importorskip("dipy")
rng = np.random.RandomState(seed)
from_data = rng.randn(*from_shape)
is_rand = False
if isinstance(to_affine, str):
assert to_affine == "rand"
to_affine = _rand_affine(rng)
is_rand = True
if isinstance(from_affine, str):
assert from_affine == "rand"
from_affine = _rand_affine(rng)
is_rand = True
to_affine = np.array(to_affine, float)
assert to_affine.shape == (4, 4)
from_affine = np.array(from_affine, float)
assert from_affine.shape == (4, 4)
#
# 1. nibabel.processing.resample_from_to
#
# for a 1mm iso / 256 -> 5mm / 51 one sample takes ~486 ms
from nibabel.processing import resample_from_to
from nibabel.spatialimages import SpatialImage
start = np.linalg.norm(from_data)
got_nibabel = resample_from_to(
SpatialImage(from_data, from_affine), (to_shape, to_affine), order=order
).get_fdata()
end = np.linalg.norm(got_nibabel)
assert end > 0.05 * start # not too much power lost
#
# 2. dipy.align.imaffine
#
# ~366 ms
import dipy.align.imaffine
interp = "linear" if order == 1 else "nearest"
got_dipy = dipy.align.imaffine.AffineMap(
None,
domain_grid_shape=to_shape,
domain_grid2world=to_affine,
codomain_grid_shape=from_shape,
codomain_grid2world=from_affine,
).transform(from_data, interpolation=interp, resample_only=True)
# XXX possibly some error in dipy or nibabel (/SciPy), or some boundary
# condition?
nib_different = (is_rand and order == 1) or (
from_affine[0, 0] == 2.0 and not np.allclose(from_affine, to_affine)
)
nib_different = nib_different and not (
is_rand and from_affine[0, 0] == 2 and order == 0
)
if nib_different:
assert not np.allclose(got_dipy, got_nibabel), "nibabel fixed"
else:
assert_allclose(got_dipy, got_nibabel, err_msg="dipy<->nibabel")
#
# 3. mne.source_space._grid_interp
#
# ~339 ms
trans = np.linalg.inv(from_affine) @ to_affine # to -> from
interp = _grid_interp(from_shape, to_shape, trans, order=order)
got_mne = np.asarray(interp @ from_data.ravel(order="F")).reshape(
to_shape, order="F"
)
if order == 1:
assert_allclose(got_mne, got_dipy, err_msg="MNE<->dipy")
else:
perc = 100 * np.isclose(got_mne, got_dipy).mean()
assert 83 < perc <= 100
Datasets
Models
