From 99217ffba2931d8563e0cdbfc1f03399bdb214e9 Mon Sep 17 00:00:00 2001 From: Xan McPherson Date: Thu, 11 Sep 2025 10:49:16 -0700 Subject: [PATCH 01/11] add new preprocessing function fit_spheres_to_mri, begin mSSS implementation in maxwell --- mne/preprocessing/fit_spheres_to_mri.py | 101 ++++++++++++++++++++++++ mne/preprocessing/maxwell.py | 19 ++++- 2 files changed, 119 insertions(+), 1 deletion(-) create mode 100644 mne/preprocessing/fit_spheres_to_mri.py diff --git a/mne/preprocessing/fit_spheres_to_mri.py b/mne/preprocessing/fit_spheres_to_mri.py new file mode 100644 index 00000000000..4464b93c81d --- /dev/null +++ b/mne/preprocessing/fit_spheres_to_mri.py @@ -0,0 +1,101 @@ +from scipy.special import KDTree + +import numpy as np +import vedo +import nibabel as nib + +from .._fiff.constants import FIFF +from ..surface import _CheckInside +from ..transforms import ( + invert_transform, + apply_trans, + read_trans, +) + +def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): + mindist = 2e-3 + assert bem[0]['id'] == FIFF.FIFFV_BEM_SURF_ID_HEAD + assert bem[2]['id'] == FIFF.FIFFV_BEM_SURF_ID_BRAIN + scalp, _, inner_skull = bem + inside_scalp = _CheckInside(scalp, mode='pyvista') + inside_skull = _CheckInside(inner_skull, mode='pyvista') + m3_to_cc = 100 ** 3 + assert inside_scalp(inner_skull['rr']).all() + assert not inside_skull(scalp['rr']).any() + b = vedo.Mesh([inner_skull['rr'], inner_skull['tris']]) + s = vedo.Mesh([scalp['rr'], scalp['tris']]) + s_tree = KDTree(scalp['rr']) + brain_volume = b.volume() + print(f'Brain vedo: {brain_volume * m3_to_cc:8.2f} cc') + brain_vol = nib.load(subjects_dir / subject / 'mri' / 'brainmask.mgz') + brain_rr = np.array(np.where(brain_vol.get_fdata())).T + brain_rr = apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000. #apply a transformation matrix + del brain_vol + brain_rr = brain_rr[inside_skull(brain_rr)] + vox_to_m3 = 1e-9 + brain_volume_vox = len(brain_rr) * vox_to_m3 + + def _print_q(title, got, want): + title = f'{title}:'.ljust(15) + print(f'{title} {got * m3_to_cc:8.2f} cc ({(want - got) / want * 100:6.2f} %)') + + _print_q('Brain vox', brain_volume_vox, brain_volume_vox) + + # 1. Compute a naive sphere using the center of mass of brain surf verts + naive_c = np.mean(inner_skull['rr'], axis=0) + + # 2. Define optimization functions + from scipy.optimize import fmin_cobyla + def _cost(c): + cs = c.reshape(-1, 3) + rs = np.maximum(s_tree.query(cs)[0] - mindist, 0.) + resid = brain_volume + mask = None + for c, r in zip(cs, rs): + if not (r and s.contains(c)): #was is_inside + continue + m = np.linalg.norm(brain_rr - c, axis=1) <= r + if mask is None: + mask = m + else: + mask |= m + resid = brain_volume_vox + if mask is not None: + resid = resid - np.sum(mask) * vox_to_m3 + return resid + + def _cons(c): + cs = c.reshape(-1, 3) + sign = np.array([2 * s.contains(c) - 1 for c in cs], float) #was "is_inside" + cons = sign * s_tree.query(cs)[0] - mindist + return cons + + # 3. Now optimize spheres and find centers + if n_spheres ==1: + x = naive_c + c_opt = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) + + elif n_spheres ==2: + c_opt_1 = fmin_cobyla(_cost, naive_c, _cons, rhobeg=1e-2, rhoend=1e-4) + x = np.concatenate([c_opt_1, naive_c]) + c_opt = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) + + elif n_spheres ==3: + print("WARNING: mSSS method has been optimized with two origins") + c_opt_1 = fmin_cobyla(_cost, naive_c, _cons, rhobeg=1e-2, rhoend=1e-4) + x = np.concatenate([c_opt_1, naive_c]) + c_opt_2 = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) + x = np.concatenate([c_opt_2, naive_c]) + c_opt = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) + else: + raise ValueError(f"Implementation is for 3 or less origins") + + + #4. transform centers for return using "trans" matrix + mri_head_t = invert_transform(read_trans(trans)) + assert mri_head_t['from'] == FIFF.FIFFV_COORD_MRI, mri_head_t['from'] + centers = apply_trans(mri_head_t, c_opt.reshape(-1, 3)) + return centers + + + \ No newline at end of file diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py index 8c270252bb2..40975755949 100644 --- a/mne/preprocessing/maxwell.py +++ b/mne/preprocessing/maxwell.py @@ -12,6 +12,7 @@ from scipy import linalg from scipy.special import lpmv + from .. import __version__ from .._fiff.compensator import make_compensator from .._fiff.constants import FIFF, FWD @@ -60,6 +61,7 @@ warn, ) + # Note: MF uses single precision and some algorithms might use # truncated versions of constants (e.g., μ0), which could lead to small # differences between algorithms @@ -1865,6 +1867,22 @@ def _sss_basis(exp, all_coils): ) return S_tot +def _combine_sss_basis(S_in1,S_in2): + """mSSS calculations using optimized multi-centers + TODO: Add some "if" statement where two different S_in basis are + calculated if "origin = more than 1D array" based on centers calculated with "prprocessing.fit_spheres_to_mri" + + """ + thresh = 0.005 + U, s, Vh = linalg.svd([S_in1,S_in2]) + S_tot = []; + #apply threshold to limit dimensions of resulting basis + for i in range(0, np.shape(s)[0]): + ratio = s[i]/s[1] + if ratio >= thresh: + S_tot = np.append(S_tot,U[:,i]) + return S_tot + def _integrate_points( cos_az, sin_az, cos_pol, sin_pol, b_r, b_az, b_pol, cosmags, bins, n_coils @@ -2919,7 +2937,6 @@ def _read_cross_talk(cross_talk, ch_names): sss_ctc["decoupler"] = sss_ctc["decoupler"].T.tocsc() return ctc, sss_ctc - @verbose def compute_maxwell_basis( info, From 9a0d348de30a7ac0095f7504cc576ff6e084f707 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 11 Sep 2025 18:09:38 +0000 Subject: [PATCH 02/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- mne/preprocessing/fit_spheres_to_mri.py | 79 ++++++++++++------------- mne/preprocessing/maxwell.py | 20 +++---- 2 files changed, 49 insertions(+), 50 deletions(-) diff --git a/mne/preprocessing/fit_spheres_to_mri.py b/mne/preprocessing/fit_spheres_to_mri.py index 4464b93c81d..cfeb04fb269 100644 --- a/mne/preprocessing/fit_spheres_to_mri.py +++ b/mne/preprocessing/fit_spheres_to_mri.py @@ -1,58 +1,61 @@ -from scipy.special import KDTree - +import nibabel as nib import numpy as np import vedo -import nibabel as nib +from scipy.special import KDTree from .._fiff.constants import FIFF from ..surface import _CheckInside from ..transforms import ( - invert_transform, apply_trans, + invert_transform, read_trans, ) + def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): mindist = 2e-3 - assert bem[0]['id'] == FIFF.FIFFV_BEM_SURF_ID_HEAD - assert bem[2]['id'] == FIFF.FIFFV_BEM_SURF_ID_BRAIN + assert bem[0]["id"] == FIFF.FIFFV_BEM_SURF_ID_HEAD + assert bem[2]["id"] == FIFF.FIFFV_BEM_SURF_ID_BRAIN scalp, _, inner_skull = bem - inside_scalp = _CheckInside(scalp, mode='pyvista') - inside_skull = _CheckInside(inner_skull, mode='pyvista') - m3_to_cc = 100 ** 3 - assert inside_scalp(inner_skull['rr']).all() - assert not inside_skull(scalp['rr']).any() - b = vedo.Mesh([inner_skull['rr'], inner_skull['tris']]) - s = vedo.Mesh([scalp['rr'], scalp['tris']]) - s_tree = KDTree(scalp['rr']) + inside_scalp = _CheckInside(scalp, mode="pyvista") + inside_skull = _CheckInside(inner_skull, mode="pyvista") + m3_to_cc = 100**3 + assert inside_scalp(inner_skull["rr"]).all() + assert not inside_skull(scalp["rr"]).any() + b = vedo.Mesh([inner_skull["rr"], inner_skull["tris"]]) + s = vedo.Mesh([scalp["rr"], scalp["tris"]]) + s_tree = KDTree(scalp["rr"]) brain_volume = b.volume() - print(f'Brain vedo: {brain_volume * m3_to_cc:8.2f} cc') - brain_vol = nib.load(subjects_dir / subject / 'mri' / 'brainmask.mgz') + print(f"Brain vedo: {brain_volume * m3_to_cc:8.2f} cc") + brain_vol = nib.load(subjects_dir / subject / "mri" / "brainmask.mgz") brain_rr = np.array(np.where(brain_vol.get_fdata())).T - brain_rr = apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000. #apply a transformation matrix + brain_rr = ( + apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000.0 + ) # apply a transformation matrix del brain_vol brain_rr = brain_rr[inside_skull(brain_rr)] vox_to_m3 = 1e-9 brain_volume_vox = len(brain_rr) * vox_to_m3 def _print_q(title, got, want): - title = f'{title}:'.ljust(15) - print(f'{title} {got * m3_to_cc:8.2f} cc ({(want - got) / want * 100:6.2f} %)') + title = f"{title}:".ljust(15) + print(f"{title} {got * m3_to_cc:8.2f} cc ({(want - got) / want * 100:6.2f} %)") - _print_q('Brain vox', brain_volume_vox, brain_volume_vox) + _print_q("Brain vox", brain_volume_vox, brain_volume_vox) # 1. Compute a naive sphere using the center of mass of brain surf verts - naive_c = np.mean(inner_skull['rr'], axis=0) - + naive_c = np.mean(inner_skull["rr"], axis=0) + # 2. Define optimization functions - from scipy.optimize import fmin_cobyla + from scipy.optimize import fmin_cobyla + def _cost(c): cs = c.reshape(-1, 3) - rs = np.maximum(s_tree.query(cs)[0] - mindist, 0.) + rs = np.maximum(s_tree.query(cs)[0] - mindist, 0.0) resid = brain_volume mask = None for c, r in zip(cs, rs): - if not (r and s.contains(c)): #was is_inside + if not (r and s.contains(c)): # was is_inside continue m = np.linalg.norm(brain_rr - c, axis=1) <= r if mask is None: @@ -66,21 +69,21 @@ def _cost(c): def _cons(c): cs = c.reshape(-1, 3) - sign = np.array([2 * s.contains(c) - 1 for c in cs], float) #was "is_inside" + sign = np.array([2 * s.contains(c) - 1 for c in cs], float) # was "is_inside" cons = sign * s_tree.query(cs)[0] - mindist return cons - + # 3. Now optimize spheres and find centers - if n_spheres ==1: + if n_spheres == 1: x = naive_c c_opt = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) - elif n_spheres ==2: + elif n_spheres == 2: c_opt_1 = fmin_cobyla(_cost, naive_c, _cons, rhobeg=1e-2, rhoend=1e-4) x = np.concatenate([c_opt_1, naive_c]) c_opt = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) - - elif n_spheres ==3: + + elif n_spheres == 3: print("WARNING: mSSS method has been optimized with two origins") c_opt_1 = fmin_cobyla(_cost, naive_c, _cons, rhobeg=1e-2, rhoend=1e-4) x = np.concatenate([c_opt_1, naive_c]) @@ -88,14 +91,10 @@ def _cons(c): x = np.concatenate([c_opt_2, naive_c]) c_opt = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) else: - raise ValueError(f"Implementation is for 3 or less origins") - - - #4. transform centers for return using "trans" matrix - mri_head_t = invert_transform(read_trans(trans)) - assert mri_head_t['from'] == FIFF.FIFFV_COORD_MRI, mri_head_t['from'] + raise ValueError("Implementation is for 3 or less origins") + + # 4. transform centers for return using "trans" matrix + mri_head_t = invert_transform(read_trans(trans)) + assert mri_head_t["from"] == FIFF.FIFFV_COORD_MRI, mri_head_t["from"] centers = apply_trans(mri_head_t, c_opt.reshape(-1, 3)) return centers - - - \ No newline at end of file diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py index 40975755949..3baf5535521 100644 --- a/mne/preprocessing/maxwell.py +++ b/mne/preprocessing/maxwell.py @@ -12,7 +12,6 @@ from scipy import linalg from scipy.special import lpmv - from .. import __version__ from .._fiff.compensator import make_compensator from .._fiff.constants import FIFF, FWD @@ -61,7 +60,6 @@ warn, ) - # Note: MF uses single precision and some algorithms might use # truncated versions of constants (e.g., μ0), which could lead to small # differences between algorithms @@ -1867,20 +1865,21 @@ def _sss_basis(exp, all_coils): ) return S_tot -def _combine_sss_basis(S_in1,S_in2): - """mSSS calculations using optimized multi-centers + +def _combine_sss_basis(S_in1, S_in2): + """MSSS calculations using optimized multi-centers TODO: Add some "if" statement where two different S_in basis are calculated if "origin = more than 1D array" based on centers calculated with "prprocessing.fit_spheres_to_mri" - + """ thresh = 0.005 - U, s, Vh = linalg.svd([S_in1,S_in2]) - S_tot = []; - #apply threshold to limit dimensions of resulting basis + U, s, Vh = linalg.svd([S_in1, S_in2]) + S_tot = [] + # apply threshold to limit dimensions of resulting basis for i in range(0, np.shape(s)[0]): - ratio = s[i]/s[1] + ratio = s[i] / s[1] if ratio >= thresh: - S_tot = np.append(S_tot,U[:,i]) + S_tot = np.append(S_tot, U[:, i]) return S_tot @@ -2937,6 +2936,7 @@ def _read_cross_talk(cross_talk, ch_names): sss_ctc["decoupler"] = sss_ctc["decoupler"].T.tocsc() return ctc, sss_ctc + @verbose def compute_maxwell_basis( info, From ee8244c4a19660b3f2cd621a8f81c15be7438186 Mon Sep 17 00:00:00 2001 From: Xan McPherson <123395570+xannnimal@users.noreply.github.com> Date: Tue, 11 Nov 2025 13:59:04 -0800 Subject: [PATCH 03/11] Update "_combine_sss_basis()" updated threshold from Matlab to Python --- mne/preprocessing/maxwell.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py index 3baf5535521..07dbdc3b07e 100644 --- a/mne/preprocessing/maxwell.py +++ b/mne/preprocessing/maxwell.py @@ -1872,15 +1872,15 @@ def _combine_sss_basis(S_in1, S_in2): calculated if "origin = more than 1D array" based on centers calculated with "prprocessing.fit_spheres_to_mri" """ - thresh = 0.005 - U, s, Vh = linalg.svd([S_in1, S_in2]) - S_tot = [] - # apply threshold to limit dimensions of resulting basis + S_tot = []; + thresh = 5e-7 #0.005 in Matlab + U, s, Vh = np.linalg.svd(np.concatenate((S_in1,S_in2),axis=1)) + #apply threshold to limit dimensions of resulting basis for i in range(0, np.shape(s)[0]): - ratio = s[i] / s[1] + ratio = s[i]/s[0] if ratio >= thresh: - S_tot = np.append(S_tot, U[:, i]) - return S_tot + S_tot.append(U[:,i]) + return np.transpose(np.array(S_tot)) def _integrate_points( From c57d9d06d3ae09ae4f11838a43a0dd44e08c5e4a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 11 Nov 2025 21:59:29 +0000 Subject: [PATCH 04/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- mne/preprocessing/maxwell.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py index 07dbdc3b07e..d4b85b6fafc 100644 --- a/mne/preprocessing/maxwell.py +++ b/mne/preprocessing/maxwell.py @@ -1872,14 +1872,14 @@ def _combine_sss_basis(S_in1, S_in2): calculated if "origin = more than 1D array" based on centers calculated with "prprocessing.fit_spheres_to_mri" """ - S_tot = []; - thresh = 5e-7 #0.005 in Matlab - U, s, Vh = np.linalg.svd(np.concatenate((S_in1,S_in2),axis=1)) - #apply threshold to limit dimensions of resulting basis + S_tot = [] + thresh = 5e-7 # 0.005 in Matlab + U, s, Vh = np.linalg.svd(np.concatenate((S_in1, S_in2), axis=1)) + # apply threshold to limit dimensions of resulting basis for i in range(0, np.shape(s)[0]): - ratio = s[i]/s[0] + ratio = s[i] / s[0] if ratio >= thresh: - S_tot.append(U[:,i]) + S_tot.append(U[:, i]) return np.transpose(np.array(S_tot)) From ac5f4a224fa8a159a88c8eda0fa677c6bd52c788 Mon Sep 17 00:00:00 2001 From: Shubhankar Seth Date: Wed, 25 Mar 2026 11:17:04 -0700 Subject: [PATCH 05/11] Implement two-origin mSSS: attempt 1 --- mne/preprocessing/fit_spheres_to_mri.py | 6 +- mne/preprocessing/maxwell.py | 101 ++++++++++++++++-------- 2 files changed, 71 insertions(+), 36 deletions(-) diff --git a/mne/preprocessing/fit_spheres_to_mri.py b/mne/preprocessing/fit_spheres_to_mri.py index cfeb04fb269..da3bf1d7a61 100644 --- a/mne/preprocessing/fit_spheres_to_mri.py +++ b/mne/preprocessing/fit_spheres_to_mri.py @@ -1,7 +1,7 @@ -import nibabel as nib -import numpy as np import vedo -from scipy.special import KDTree +import numpy as np +import nibabel as nib +from scipy.spatial import KDTree from .._fiff.constants import FIFF from ..surface import _CheckInside diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py index d4b85b6fafc..f619e2de22b 100644 --- a/mne/preprocessing/maxwell.py +++ b/mne/preprocessing/maxwell.py @@ -403,34 +403,41 @@ def maxwell_filter( .. footbibliography:: """ # noqa: E501 logger.info("Maxwell filtering raw data") - params = _prep_maxwell_filter( - raw=raw, - origin=origin, - int_order=int_order, - ext_order=ext_order, - calibration=calibration, - cross_talk=cross_talk, - st_duration=st_duration, - st_correlation=st_correlation, - coord_frame=coord_frame, - destination=destination, - regularize=regularize, - ignore_ref=ignore_ref, - bad_condition=bad_condition, - head_pos=head_pos, - st_fixed=st_fixed, - st_only=st_only, - mag_scale=mag_scale, - skip_by_annotation=skip_by_annotation, - extended_proj=extended_proj, - st_overlap=st_overlap, - mc_interp=mc_interp, - ) - raw_sss = _run_maxwell_filter(raw, **params) - # Update info - _update_sss_info(raw_sss, **params["update_kwargs"]) - logger.info("[done]") - return raw_sss + # TODO: fix + if isinstance(origin, np.ndarray): + raw_mSSS = _run_mSSS(raw, origin) + # Update info _update_sss_info(raw_sss, **params["update_kwargs"]) ?? + return raw_mSSS + + else: + params = _prep_maxwell_filter( + raw=raw, + origin=origin, + int_order=int_order, + ext_order=ext_order, + calibration=calibration, + cross_talk=cross_talk, + st_duration=st_duration, + st_correlation=st_correlation, + coord_frame=coord_frame, + destination=destination, + regularize=regularize, + ignore_ref=ignore_ref, + bad_condition=bad_condition, + head_pos=head_pos, + st_fixed=st_fixed, + st_only=st_only, + mag_scale=mag_scale, + skip_by_annotation=skip_by_annotation, + extended_proj=extended_proj, + st_overlap=st_overlap, + mc_interp=mc_interp, + ) + raw_sss = _run_maxwell_filter(raw, **params) + # Update info + _update_sss_info(raw_sss, **params["update_kwargs"]) + logger.info("[done]") + return raw_sss @verbose @@ -834,6 +841,37 @@ def _run_maxwell_filter( return raw_sss +def _run_mSSS(raw, origin): + if len(origin) > 2: + # TODO: fix error msg + raise ValueError("n > 2") + + S_in = [] + for i in range(len(origin)): + [S_i, _, _, moments_i] = compute_maxwell_basis( + raw.info, origin=origin[i], regularize=None + ) + S_in.append(S_i[:, :moments_i]) + + S_tot = _combine_sss_basis(S_in[0], S_in[1]) + + raw_mSSS = maxwell_filter( + raw, origin="auto", regularize=None, ignore_ref=True, bad_condition="ignore" + ) + + meg_picks = pick_types(raw_mSSS.info, meg=True) + + # reconstruct + phi_0 = raw.get_data(picks="meg") + pS = np.linalg.pinv(S_tot) + XN = pS @ phi_0 + mSSS_data = np.real(S_tot @ XN) + + raw_mSSS._data[meg_picks] = mSSS_data + + return raw_mSSS + + class _MoveComp: """Perform movement compensation.""" @@ -1867,11 +1905,8 @@ def _sss_basis(exp, all_coils): def _combine_sss_basis(S_in1, S_in2): - """MSSS calculations using optimized multi-centers - TODO: Add some "if" statement where two different S_in basis are - calculated if "origin = more than 1D array" based on centers calculated with "prprocessing.fit_spheres_to_mri" - - """ + """mSSS calculations using optimized multi-centers""" + # TODO: n > 2 S_tot = [] thresh = 5e-7 # 0.005 in Matlab U, s, Vh = np.linalg.svd(np.concatenate((S_in1, S_in2), axis=1)) From a21b90b71803363bd701e4c9be7402aa4728d748 Mon Sep 17 00:00:00 2001 From: xannnimal Date: Thu, 9 Apr 2026 17:11:51 -0700 Subject: [PATCH 06/11] change bad condition error to ignore in mSSS calc --- mne/preprocessing/fit_spheres_to_mri.py | 28 +++++++++++++++++++++++++ mne/preprocessing/maxwell.py | 2 +- 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/mne/preprocessing/fit_spheres_to_mri.py b/mne/preprocessing/fit_spheres_to_mri.py index da3bf1d7a61..ba88186095c 100644 --- a/mne/preprocessing/fit_spheres_to_mri.py +++ b/mne/preprocessing/fit_spheres_to_mri.py @@ -13,6 +13,34 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): + """Fits two spheres to MRI using BEM, such that spheres fit while brain but + do not encroach on sensors. For use with Milti-SSS Maxwell Filtering + + Parameters + ---------- + subjects dir: str + director to Freesurfer subjects + subject: str + Subject ID + bem: bem.ConductorModel + istance of bem.ConductorModel, must be three shell conductivity profiles + trans: str + path to trans file, mri_dev_t information + n_spheres: int + number of spheres to fit, recommended 2 + + Returns + ------- + centers: np.ndarray + 2D array containing the two centers in cartesian coordinates + + Notes + ----- + + * Must have vedo and nibabel installed + * Must have run mne watershed BEM using freesurfer segmentation + """ + mindist = 2e-3 assert bem[0]["id"] == FIFF.FIFFV_BEM_SURF_ID_HEAD assert bem[2]["id"] == FIFF.FIFFV_BEM_SURF_ID_BRAIN diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py index f619e2de22b..dc5757d4f4d 100644 --- a/mne/preprocessing/maxwell.py +++ b/mne/preprocessing/maxwell.py @@ -849,7 +849,7 @@ def _run_mSSS(raw, origin): S_in = [] for i in range(len(origin)): [S_i, _, _, moments_i] = compute_maxwell_basis( - raw.info, origin=origin[i], regularize=None + raw.info, origin=origin[i], regularize=None, bad_condition="ignore" ) S_in.append(S_i[:, :moments_i]) From 8b8ebc5b4f42fcda8a5edabbc7996aee5a9ae9b1 Mon Sep 17 00:00:00 2001 From: xannnimal Date: Wed, 15 Apr 2026 13:11:57 -0700 Subject: [PATCH 07/11] update fit_spheres for MRI_HEAD_T inversion check --- mne/preprocessing/fit_spheres_to_mri.py | 56 +++++++++++++++---------- 1 file changed, 34 insertions(+), 22 deletions(-) diff --git a/mne/preprocessing/fit_spheres_to_mri.py b/mne/preprocessing/fit_spheres_to_mri.py index ba88186095c..95cd53b912c 100644 --- a/mne/preprocessing/fit_spheres_to_mri.py +++ b/mne/preprocessing/fit_spheres_to_mri.py @@ -1,17 +1,3 @@ -import vedo -import numpy as np -import nibabel as nib -from scipy.spatial import KDTree - -from .._fiff.constants import FIFF -from ..surface import _CheckInside -from ..transforms import ( - apply_trans, - invert_transform, - read_trans, -) - - def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): """Fits two spheres to MRI using BEM, such that spheres fit while brain but do not encroach on sensors. For use with Milti-SSS Maxwell Filtering @@ -22,8 +8,8 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): director to Freesurfer subjects subject: str Subject ID - bem: bem.ConductorModel - istance of bem.ConductorModel, must be three shell conductivity profiles + bem: list + output of mne.make_bem_model(), must be three shell conductivity profiles trans: str path to trans file, mri_dev_t information n_spheres: int @@ -32,15 +18,35 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): Returns ------- centers: np.ndarray - 2D array containing the two centers in cartesian coordinates + 2x3 array containing the two centers in HEAD coordinate space + can be directly fed into: + raw_msss = mne.preprocessing.maxwell_filter(raw, origin=origins, ...) + for multi-SSS preprocessing + Notes ----- - * Must have vedo and nibabel installed * Must have run mne watershed BEM using freesurfer segmentation """ + ## --- required imports + import vedo + import os + import numpy as np + import nibabel as nib + from scipy.spatial import KDTree + + from .._fiff.constants import FIFF + from ..surface import _CheckInside + from ..transforms import ( + apply_trans, + invert_transform, + read_trans, + ) + + + ## --- begin mindist = 2e-3 assert bem[0]["id"] == FIFF.FIFFV_BEM_SURF_ID_HEAD assert bem[2]["id"] == FIFF.FIFFV_BEM_SURF_ID_BRAIN @@ -55,7 +61,7 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): s_tree = KDTree(scalp["rr"]) brain_volume = b.volume() print(f"Brain vedo: {brain_volume * m3_to_cc:8.2f} cc") - brain_vol = nib.load(subjects_dir / subject / "mri" / "brainmask.mgz") + brain_vol = nib.load(os.path.join(subjects_dir, subject, 'mri','brainmask.mgz')) brain_rr = np.array(np.where(brain_vol.get_fdata())).T brain_rr = ( apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000.0 @@ -121,8 +127,14 @@ def _cons(c): else: raise ValueError("Implementation is for 3 or less origins") - # 4. transform centers for return using "trans" matrix - mri_head_t = invert_transform(read_trans(trans)) - assert mri_head_t["from"] == FIFF.FIFFV_COORD_MRI, mri_head_t["from"] + # 4. transform centers for return using "mri_head_t" matrix + # Output 2-sphere case + mri_head_t = invert_transform(read_trans(trans)) #trans is the raw_fif_trans file + if mri_head_t["from"] == FIFF.FIFFV_COORD_HEAD: + mri_head_t = invert_transform(mri_head_t) + assert mri_head_t['from'] == FIFF.FIFFV_COORD_MRI, mri_head_t['from'] centers = apply_trans(mri_head_t, c_opt.reshape(-1, 3)) return centers + + + From e5ef98372b2740643f20310f765342b73468e5fc Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 16 Apr 2026 18:09:41 +0000 Subject: [PATCH 08/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- mne/preprocessing/fit_spheres_to_mri.py | 20 ++++++++------------ mne/preprocessing/maxwell.py | 2 +- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/mne/preprocessing/fit_spheres_to_mri.py b/mne/preprocessing/fit_spheres_to_mri.py index 95cd53b912c..8e38425cf06 100644 --- a/mne/preprocessing/fit_spheres_to_mri.py +++ b/mne/preprocessing/fit_spheres_to_mri.py @@ -1,5 +1,5 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): - """Fits two spheres to MRI using BEM, such that spheres fit while brain but + """Fits two spheres to MRI using BEM, such that spheres fit while brain but do not encroach on sensors. For use with Milti-SSS Maxwell Filtering Parameters @@ -28,13 +28,13 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): ----- * Must have vedo and nibabel installed * Must have run mne watershed BEM using freesurfer segmentation - """ - + """ ## --- required imports - import vedo import os - import numpy as np + import nibabel as nib + import numpy as np + import vedo from scipy.spatial import KDTree from .._fiff.constants import FIFF @@ -45,7 +45,6 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): read_trans, ) - ## --- begin mindist = 2e-3 assert bem[0]["id"] == FIFF.FIFFV_BEM_SURF_ID_HEAD @@ -61,7 +60,7 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): s_tree = KDTree(scalp["rr"]) brain_volume = b.volume() print(f"Brain vedo: {brain_volume * m3_to_cc:8.2f} cc") - brain_vol = nib.load(os.path.join(subjects_dir, subject, 'mri','brainmask.mgz')) + brain_vol = nib.load(os.path.join(subjects_dir, subject, "mri", "brainmask.mgz")) brain_rr = np.array(np.where(brain_vol.get_fdata())).T brain_rr = ( apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000.0 @@ -129,12 +128,9 @@ def _cons(c): # 4. transform centers for return using "mri_head_t" matrix # Output 2-sphere case - mri_head_t = invert_transform(read_trans(trans)) #trans is the raw_fif_trans file + mri_head_t = invert_transform(read_trans(trans)) # trans is the raw_fif_trans file if mri_head_t["from"] == FIFF.FIFFV_COORD_HEAD: mri_head_t = invert_transform(mri_head_t) - assert mri_head_t['from'] == FIFF.FIFFV_COORD_MRI, mri_head_t['from'] + assert mri_head_t["from"] == FIFF.FIFFV_COORD_MRI, mri_head_t["from"] centers = apply_trans(mri_head_t, c_opt.reshape(-1, 3)) return centers - - - diff --git a/mne/preprocessing/maxwell.py b/mne/preprocessing/maxwell.py index dc5757d4f4d..7f835400eb8 100644 --- a/mne/preprocessing/maxwell.py +++ b/mne/preprocessing/maxwell.py @@ -1905,7 +1905,7 @@ def _sss_basis(exp, all_coils): def _combine_sss_basis(S_in1, S_in2): - """mSSS calculations using optimized multi-centers""" + """MSSS calculations using optimized multi-centers""" # TODO: n > 2 S_tot = [] thresh = 5e-7 # 0.005 in Matlab From ef1e45718f5ce83fed310b1bb4a2050b3585cf05 Mon Sep 17 00:00:00 2001 From: xannnimal Date: Mon, 4 May 2026 12:59:38 -0700 Subject: [PATCH 09/11] update fit_spheres to include option for vizualization of results --- mne/preprocessing/fit_spheres_to_mri.py | 155 +++++++++++++++--------- 1 file changed, 99 insertions(+), 56 deletions(-) diff --git a/mne/preprocessing/fit_spheres_to_mri.py b/mne/preprocessing/fit_spheres_to_mri.py index 8e38425cf06..643ef04c20a 100644 --- a/mne/preprocessing/fit_spheres_to_mri.py +++ b/mne/preprocessing/fit_spheres_to_mri.py @@ -1,4 +1,4 @@ -def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): +def fit_spheres_to_mri(subjects_dir, subject, bem_surf, trans, n_spheres,show_spheres=False): """Fits two spheres to MRI using BEM, such that spheres fit while brain but do not encroach on sensors. For use with Milti-SSS Maxwell Filtering @@ -8,12 +8,14 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): director to Freesurfer subjects subject: str Subject ID - bem: list + bem_surf: list output of mne.make_bem_model(), must be three shell conductivity profiles trans: str path to trans file, mri_dev_t information n_spheres: int number of spheres to fit, recommended 2 + show_spheres: bool + show pyvista plot of the origins and optimized spheres overlayed with the head Returns ------- @@ -30,7 +32,6 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): * Must have run mne watershed BEM using freesurfer segmentation """ ## --- required imports - import os import nibabel as nib import numpy as np @@ -47,48 +48,58 @@ def fit_spheres_to_mri(subjects_dir, subject, bem, trans, n_spheres): ## --- begin mindist = 2e-3 - assert bem[0]["id"] == FIFF.FIFFV_BEM_SURF_ID_HEAD - assert bem[2]["id"] == FIFF.FIFFV_BEM_SURF_ID_BRAIN - scalp, _, inner_skull = bem - inside_scalp = _CheckInside(scalp, mode="pyvista") - inside_skull = _CheckInside(inner_skull, mode="pyvista") - m3_to_cc = 100**3 - assert inside_scalp(inner_skull["rr"]).all() - assert not inside_skull(scalp["rr"]).any() - b = vedo.Mesh([inner_skull["rr"], inner_skull["tris"]]) - s = vedo.Mesh([scalp["rr"], scalp["tris"]]) - s_tree = KDTree(scalp["rr"]) + assert bem_surf[0]['id'] == FIFF.FIFFV_BEM_SURF_ID_HEAD + assert bem_surf[2]['id'] == FIFF.FIFFV_BEM_SURF_ID_BRAIN + scalp, _, inner_skull = bem_surf + inside_scalp = _CheckInside(scalp, mode='pyvista') + inside_skull = _CheckInside(inner_skull, mode='pyvista') + m3_to_cc = 100 ** 3 + assert inside_scalp(inner_skull['rr']).all() + assert not inside_skull(scalp['rr']).any() + b = vedo.Mesh([inner_skull['rr'], inner_skull['tris']]) + s = vedo.Mesh([scalp['rr'], scalp['tris']]) + s_tree = KDTree(scalp['rr']) brain_volume = b.volume() - print(f"Brain vedo: {brain_volume * m3_to_cc:8.2f} cc") - brain_vol = nib.load(os.path.join(subjects_dir, subject, "mri", "brainmask.mgz")) + print(f'Brain vedo: {brain_volume * m3_to_cc:8.2f} cc') + brain_vol = nib.load(subjects_dir / subject / 'mri' / 'brainmask.mgz') brain_rr = np.array(np.where(brain_vol.get_fdata())).T - brain_rr = ( - apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000.0 - ) # apply a transformation matrix - del brain_vol + brain_rr = apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000. #apply a transformation matrix + del brain_vol #delete brain volume brain_rr = brain_rr[inside_skull(brain_rr)] vox_to_m3 = 1e-9 brain_volume_vox = len(brain_rr) * vox_to_m3 def _print_q(title, got, want): - title = f"{title}:".ljust(15) - print(f"{title} {got * m3_to_cc:8.2f} cc ({(want - got) / want * 100:6.2f} %)") + title = f'{title}:'.ljust(15) + print(f'{title} {got * m3_to_cc:8.2f} cc ({(want - got) / want * 100:6.2f} %)') - _print_q("Brain vox", brain_volume_vox, brain_volume_vox) + _print_q('Brain vox', brain_volume_vox, brain_volume_vox) # 1. Compute a naive sphere using the center of mass of brain surf verts - naive_c = np.mean(inner_skull["rr"], axis=0) - - # 2. Define optimization functions - from scipy.optimize import fmin_cobyla + naive_c = np.mean(inner_skull['rr'], axis=0) + naive_r = np.min(np.linalg.norm(inner_skull['rr'] - naive_c, axis=1)) + naive_v = 4 / 3 * np.pi * naive_r ** 3 + _print_q('Naive sphere', naive_v, brain_volume) + s1 = vedo.Sphere(naive_c, naive_r, res=100) + _print_q('Naive vedo', s1.volume(), brain_volume) + + # 2. Now use the larger radius (to head) plus mesh arithmetic + better_r = s_tree.query(naive_c)[0] - mindist + s1 = vedo.Sphere(naive_c, better_r, res=24) + _print_q('Better vedo', s1.boolean("intersect", b).volume(), brain_volume) + v = np.sum(np.linalg.norm(brain_rr - naive_c, axis=1) <= better_r) * vox_to_m3 + _print_q('Better vox', v, brain_volume_vox) + + # 3. Now optimize one sphere + from scipy.optimize import fmin_cobyla #constrained optimization by linear approximation def _cost(c): cs = c.reshape(-1, 3) - rs = np.maximum(s_tree.query(cs)[0] - mindist, 0.0) + rs = np.maximum(s_tree.query(cs)[0] - mindist, 0.) resid = brain_volume mask = None for c, r in zip(cs, rs): - if not (r and s.contains(c)): # was is_inside + if not (r and s.contains(c)): #was is_inside continue m = np.linalg.norm(brain_rr - c, axis=1) <= r if mask is None: @@ -102,35 +113,67 @@ def _cost(c): def _cons(c): cs = c.reshape(-1, 3) - sign = np.array([2 * s.contains(c) - 1 for c in cs], float) # was "is_inside" + sign = np.array([2 * s.contains(c) - 1 for c in cs], float) #was "is_inside" cons = sign * s_tree.query(cs)[0] - mindist return cons - # 3. Now optimize spheres and find centers - if n_spheres == 1: - x = naive_c - c_opt = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) - - elif n_spheres == 2: - c_opt_1 = fmin_cobyla(_cost, naive_c, _cons, rhobeg=1e-2, rhoend=1e-4) - x = np.concatenate([c_opt_1, naive_c]) - c_opt = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) - - elif n_spheres == 3: - print("WARNING: mSSS method has been optimized with two origins") - c_opt_1 = fmin_cobyla(_cost, naive_c, _cons, rhobeg=1e-2, rhoend=1e-4) - x = np.concatenate([c_opt_1, naive_c]) - c_opt_2 = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) - x = np.concatenate([c_opt_2, naive_c]) - c_opt = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) - else: - raise ValueError("Implementation is for 3 or less origins") - - # 4. transform centers for return using "mri_head_t" matrix - # Output 2-sphere case - mri_head_t = invert_transform(read_trans(trans)) # trans is the raw_fif_trans file + x = naive_c + c_opt_1 = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) + v_opt_1 = brain_volume_vox - _cost(c_opt_1) + _print_q('COBYLA 1', v_opt_1, brain_volume_vox) + + # 4. Now optimize two spheres + x = np.concatenate([c_opt_1, naive_c]) + c_opt_2 = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) + v_opt_2 = brain_volume_vox - _cost(c_opt_2) + _print_q('COBYLA 2', v_opt_2, brain_volume_vox) + + # 4. Finally, three spheres (not perfect, not global opt) + x = np.concatenate([c_opt_2, naive_c]) + c_opt_3 = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) + v_opt_3 = brain_volume_vox - _cost(c_opt_3) + _print_q('COBYLA 3', v_opt_3, brain_volume_vox) + + if show_spheres: + import pyvista as pv + import pyvistaqt + import matplotlib + plotter = pyvistaqt.BackgroundPlotter( + shape=(1, 2), window_size=(1200, 300), + editor=False, menu_bar=False, toolbar=False) + plotter.background_color = 'w' + brain_mesh = pv.make_tri_mesh(inner_skull['rr'], inner_skull['tris']) + scalp_mesh = pv.make_tri_mesh(scalp['rr'], scalp['tris']) + colors = matplotlib.rcParams['axes.prop_cycle'].by_key()['color'] + mesh_kwargs = dict(render=False, reset_camera=False, smooth_shading=True) + for ci, cs in enumerate((c_opt_1, c_opt_2, c_opt_3)): + plotter.subplot(0, ci) + plotter.camera.position = (0., -0.5, 0) + plotter.camera.focal_point = (0., 0., 0.) + plotter.camera.azimuth = 90 + plotter.camera.elevation = 0 + plotter.camera.up = (0., 0., 1.) + plotter.add_mesh(brain_mesh, opacity=0.2, color='k', **mesh_kwargs) + plotter.add_mesh(scalp_mesh, opacity=0.1, color='tan', **mesh_kwargs) + for c, color in zip(cs.reshape(-1, 3), colors): + sphere = pv.Sphere(s_tree.query(c)[0] - mindist, c) + plotter.add_mesh(sphere, opacity=0.5, color=color, **mesh_kwargs) + plotter.show() + + # Ready centers to output, transform into device space + mri_head_t = invert_transform(read_trans(trans)) if mri_head_t["from"] == FIFF.FIFFV_COORD_HEAD: mri_head_t = invert_transform(mri_head_t) - assert mri_head_t["from"] == FIFF.FIFFV_COORD_MRI, mri_head_t["from"] - centers = apply_trans(mri_head_t, c_opt.reshape(-1, 3)) - return centers + assert mri_head_t['from'] == FIFF.FIFFV_COORD_MRI, mri_head_t['from'] + centers=[] + for use in (c_opt_1,c_opt_2,c_opt_3): + centers.append(apply_trans(mri_head_t, use.reshape(-1, 3))) + if n_spheres==1: + return centers[0] + if n_spheres==2: + return centers[1] + if n_spheres==3: + print("Warning: use of mSSS with three origins and expansions is not tested or recommended") + return centers[2] + + From bc2ce19d3431a67af02b6f9bc5f589f37d38c96e Mon Sep 17 00:00:00 2001 From: xannnimal Date: Wed, 6 May 2026 14:33:29 -0700 Subject: [PATCH 10/11] fix path os parsing --- mne/preprocessing/fit_spheres_to_mri.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mne/preprocessing/fit_spheres_to_mri.py b/mne/preprocessing/fit_spheres_to_mri.py index 643ef04c20a..2bf24465017 100644 --- a/mne/preprocessing/fit_spheres_to_mri.py +++ b/mne/preprocessing/fit_spheres_to_mri.py @@ -34,6 +34,7 @@ def fit_spheres_to_mri(subjects_dir, subject, bem_surf, trans, n_spheres,show_sp ## --- required imports import nibabel as nib + import os import numpy as np import vedo from scipy.spatial import KDTree @@ -61,7 +62,7 @@ def fit_spheres_to_mri(subjects_dir, subject, bem_surf, trans, n_spheres,show_sp s_tree = KDTree(scalp['rr']) brain_volume = b.volume() print(f'Brain vedo: {brain_volume * m3_to_cc:8.2f} cc') - brain_vol = nib.load(subjects_dir / subject / 'mri' / 'brainmask.mgz') + brain_vol = nib.load(os.path.join(subjects_dir,subject,'mri','brainmask.mgz')) brain_rr = np.array(np.where(brain_vol.get_fdata())).T brain_rr = apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000. #apply a transformation matrix del brain_vol #delete brain volume From 9501204154cad568bf218714999555b3c081645a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 6 May 2026 21:34:57 +0000 Subject: [PATCH 11/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- mne/preprocessing/fit_spheres_to_mri.py | 122 +++++++++++++----------- 1 file changed, 67 insertions(+), 55 deletions(-) diff --git a/mne/preprocessing/fit_spheres_to_mri.py b/mne/preprocessing/fit_spheres_to_mri.py index 2bf24465017..396954d8e3d 100644 --- a/mne/preprocessing/fit_spheres_to_mri.py +++ b/mne/preprocessing/fit_spheres_to_mri.py @@ -1,4 +1,6 @@ -def fit_spheres_to_mri(subjects_dir, subject, bem_surf, trans, n_spheres,show_spheres=False): +def fit_spheres_to_mri( + subjects_dir, subject, bem_surf, trans, n_spheres, show_spheres=False +): """Fits two spheres to MRI using BEM, such that spheres fit while brain but do not encroach on sensors. For use with Milti-SSS Maxwell Filtering @@ -33,8 +35,9 @@ def fit_spheres_to_mri(subjects_dir, subject, bem_surf, trans, n_spheres,show_sp """ ## --- required imports - import nibabel as nib import os + + import nibabel as nib import numpy as np import vedo from scipy.spatial import KDTree @@ -49,58 +52,62 @@ def fit_spheres_to_mri(subjects_dir, subject, bem_surf, trans, n_spheres,show_sp ## --- begin mindist = 2e-3 - assert bem_surf[0]['id'] == FIFF.FIFFV_BEM_SURF_ID_HEAD - assert bem_surf[2]['id'] == FIFF.FIFFV_BEM_SURF_ID_BRAIN + assert bem_surf[0]["id"] == FIFF.FIFFV_BEM_SURF_ID_HEAD + assert bem_surf[2]["id"] == FIFF.FIFFV_BEM_SURF_ID_BRAIN scalp, _, inner_skull = bem_surf - inside_scalp = _CheckInside(scalp, mode='pyvista') - inside_skull = _CheckInside(inner_skull, mode='pyvista') - m3_to_cc = 100 ** 3 - assert inside_scalp(inner_skull['rr']).all() - assert not inside_skull(scalp['rr']).any() - b = vedo.Mesh([inner_skull['rr'], inner_skull['tris']]) - s = vedo.Mesh([scalp['rr'], scalp['tris']]) - s_tree = KDTree(scalp['rr']) + inside_scalp = _CheckInside(scalp, mode="pyvista") + inside_skull = _CheckInside(inner_skull, mode="pyvista") + m3_to_cc = 100**3 + assert inside_scalp(inner_skull["rr"]).all() + assert not inside_skull(scalp["rr"]).any() + b = vedo.Mesh([inner_skull["rr"], inner_skull["tris"]]) + s = vedo.Mesh([scalp["rr"], scalp["tris"]]) + s_tree = KDTree(scalp["rr"]) brain_volume = b.volume() - print(f'Brain vedo: {brain_volume * m3_to_cc:8.2f} cc') - brain_vol = nib.load(os.path.join(subjects_dir,subject,'mri','brainmask.mgz')) + print(f"Brain vedo: {brain_volume * m3_to_cc:8.2f} cc") + brain_vol = nib.load(os.path.join(subjects_dir, subject, "mri", "brainmask.mgz")) brain_rr = np.array(np.where(brain_vol.get_fdata())).T - brain_rr = apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000. #apply a transformation matrix - del brain_vol #delete brain volume + brain_rr = ( + apply_trans(brain_vol.header.get_vox2ras_tkr(), brain_rr) / 1000.0 + ) # apply a transformation matrix + del brain_vol # delete brain volume brain_rr = brain_rr[inside_skull(brain_rr)] vox_to_m3 = 1e-9 brain_volume_vox = len(brain_rr) * vox_to_m3 def _print_q(title, got, want): - title = f'{title}:'.ljust(15) - print(f'{title} {got * m3_to_cc:8.2f} cc ({(want - got) / want * 100:6.2f} %)') + title = f"{title}:".ljust(15) + print(f"{title} {got * m3_to_cc:8.2f} cc ({(want - got) / want * 100:6.2f} %)") - _print_q('Brain vox', brain_volume_vox, brain_volume_vox) + _print_q("Brain vox", brain_volume_vox, brain_volume_vox) # 1. Compute a naive sphere using the center of mass of brain surf verts - naive_c = np.mean(inner_skull['rr'], axis=0) - naive_r = np.min(np.linalg.norm(inner_skull['rr'] - naive_c, axis=1)) - naive_v = 4 / 3 * np.pi * naive_r ** 3 - _print_q('Naive sphere', naive_v, brain_volume) + naive_c = np.mean(inner_skull["rr"], axis=0) + naive_r = np.min(np.linalg.norm(inner_skull["rr"] - naive_c, axis=1)) + naive_v = 4 / 3 * np.pi * naive_r**3 + _print_q("Naive sphere", naive_v, brain_volume) s1 = vedo.Sphere(naive_c, naive_r, res=100) - _print_q('Naive vedo', s1.volume(), brain_volume) + _print_q("Naive vedo", s1.volume(), brain_volume) # 2. Now use the larger radius (to head) plus mesh arithmetic better_r = s_tree.query(naive_c)[0] - mindist s1 = vedo.Sphere(naive_c, better_r, res=24) - _print_q('Better vedo', s1.boolean("intersect", b).volume(), brain_volume) + _print_q("Better vedo", s1.boolean("intersect", b).volume(), brain_volume) v = np.sum(np.linalg.norm(brain_rr - naive_c, axis=1) <= better_r) * vox_to_m3 - _print_q('Better vox', v, brain_volume_vox) + _print_q("Better vox", v, brain_volume_vox) # 3. Now optimize one sphere - from scipy.optimize import fmin_cobyla #constrained optimization by linear approximation + from scipy.optimize import ( + fmin_cobyla, # constrained optimization by linear approximation + ) def _cost(c): cs = c.reshape(-1, 3) - rs = np.maximum(s_tree.query(cs)[0] - mindist, 0.) + rs = np.maximum(s_tree.query(cs)[0] - mindist, 0.0) resid = brain_volume mask = None for c, r in zip(cs, rs): - if not (r and s.contains(c)): #was is_inside + if not (r and s.contains(c)): # was is_inside continue m = np.linalg.norm(brain_rr - c, axis=1) <= r if mask is None: @@ -114,67 +121,72 @@ def _cost(c): def _cons(c): cs = c.reshape(-1, 3) - sign = np.array([2 * s.contains(c) - 1 for c in cs], float) #was "is_inside" + sign = np.array([2 * s.contains(c) - 1 for c in cs], float) # was "is_inside" cons = sign * s_tree.query(cs)[0] - mindist return cons x = naive_c c_opt_1 = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) v_opt_1 = brain_volume_vox - _cost(c_opt_1) - _print_q('COBYLA 1', v_opt_1, brain_volume_vox) + _print_q("COBYLA 1", v_opt_1, brain_volume_vox) # 4. Now optimize two spheres x = np.concatenate([c_opt_1, naive_c]) c_opt_2 = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) v_opt_2 = brain_volume_vox - _cost(c_opt_2) - _print_q('COBYLA 2', v_opt_2, brain_volume_vox) + _print_q("COBYLA 2", v_opt_2, brain_volume_vox) # 4. Finally, three spheres (not perfect, not global opt) x = np.concatenate([c_opt_2, naive_c]) c_opt_3 = fmin_cobyla(_cost, x, _cons, rhobeg=1e-2, rhoend=1e-4) v_opt_3 = brain_volume_vox - _cost(c_opt_3) - _print_q('COBYLA 3', v_opt_3, brain_volume_vox) + _print_q("COBYLA 3", v_opt_3, brain_volume_vox) if show_spheres: + import matplotlib import pyvista as pv import pyvistaqt - import matplotlib + plotter = pyvistaqt.BackgroundPlotter( - shape=(1, 2), window_size=(1200, 300), - editor=False, menu_bar=False, toolbar=False) - plotter.background_color = 'w' - brain_mesh = pv.make_tri_mesh(inner_skull['rr'], inner_skull['tris']) - scalp_mesh = pv.make_tri_mesh(scalp['rr'], scalp['tris']) - colors = matplotlib.rcParams['axes.prop_cycle'].by_key()['color'] + shape=(1, 2), + window_size=(1200, 300), + editor=False, + menu_bar=False, + toolbar=False, + ) + plotter.background_color = "w" + brain_mesh = pv.make_tri_mesh(inner_skull["rr"], inner_skull["tris"]) + scalp_mesh = pv.make_tri_mesh(scalp["rr"], scalp["tris"]) + colors = matplotlib.rcParams["axes.prop_cycle"].by_key()["color"] mesh_kwargs = dict(render=False, reset_camera=False, smooth_shading=True) for ci, cs in enumerate((c_opt_1, c_opt_2, c_opt_3)): plotter.subplot(0, ci) - plotter.camera.position = (0., -0.5, 0) - plotter.camera.focal_point = (0., 0., 0.) + plotter.camera.position = (0.0, -0.5, 0) + plotter.camera.focal_point = (0.0, 0.0, 0.0) plotter.camera.azimuth = 90 plotter.camera.elevation = 0 - plotter.camera.up = (0., 0., 1.) - plotter.add_mesh(brain_mesh, opacity=0.2, color='k', **mesh_kwargs) - plotter.add_mesh(scalp_mesh, opacity=0.1, color='tan', **mesh_kwargs) + plotter.camera.up = (0.0, 0.0, 1.0) + plotter.add_mesh(brain_mesh, opacity=0.2, color="k", **mesh_kwargs) + plotter.add_mesh(scalp_mesh, opacity=0.1, color="tan", **mesh_kwargs) for c, color in zip(cs.reshape(-1, 3), colors): sphere = pv.Sphere(s_tree.query(c)[0] - mindist, c) plotter.add_mesh(sphere, opacity=0.5, color=color, **mesh_kwargs) plotter.show() # Ready centers to output, transform into device space - mri_head_t = invert_transform(read_trans(trans)) + mri_head_t = invert_transform(read_trans(trans)) if mri_head_t["from"] == FIFF.FIFFV_COORD_HEAD: mri_head_t = invert_transform(mri_head_t) - assert mri_head_t['from'] == FIFF.FIFFV_COORD_MRI, mri_head_t['from'] - centers=[] - for use in (c_opt_1,c_opt_2,c_opt_3): + assert mri_head_t["from"] == FIFF.FIFFV_COORD_MRI, mri_head_t["from"] + centers = [] + for use in (c_opt_1, c_opt_2, c_opt_3): centers.append(apply_trans(mri_head_t, use.reshape(-1, 3))) - if n_spheres==1: + if n_spheres == 1: return centers[0] - if n_spheres==2: + if n_spheres == 2: return centers[1] - if n_spheres==3: - print("Warning: use of mSSS with three origins and expansions is not tested or recommended") + if n_spheres == 3: + print( + "Warning: use of mSSS with three origins and expansions is not tested or recommended" + ) return centers[2] - -