Hi!
I am currently processing MEG empty room data for an experiment. For some of the participants, the digitization points were kept in the acquisition software while doing the empty room, which leads to issues later in the pipeline.
A solution we found was to use the method raw.set_montage(None), which as I understood would remove the digitization points from the data.
However, when running the following:
auto_noisy_channels_er, auto_flat_channels_er, auto_scores_er = mne.preprocessing.find_bad_channels_maxwell(
raw_bids_empty_room, # Raw data to process
cross_talk = crosstalk_file, # Specific to the MEG machine
calibration = fine_cal_file, # Specific to the MEG machine
return_scores = True,
coord_frame = "meg", # Necessary for empty room recordings
duration = 30, # Duration of each window (in seconds)
min_count = 10, # Number of window appearances to be counted as bad channel
verbose = False,
)
I get the following error:
ValueError Traceback (most recent call last)
/1_Preprocessing_batch.ipynb Cell 15 line 2
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=251'>252</a> print(\"Running the automatic detection of bad channels...\")
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=252'>253</a> t_step = dt.now()
--> <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=253'>254</a> auto_noisy_channels_er, auto_flat_channels_er, auto_scores_er = mne.preprocessing.find_bad_channels_maxwell(
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=254'>255</a> raw_bids_empty_room, # Raw data to process
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=255'>256</a> cross_talk = subject.crosstalk_file,
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=256'>257</a> calibration = subject.fine_cal_file,
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=257'>258</a> return_scores = True,
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=258'>259</a> coord_frame = \"meg\", # Necessary for empty room recordings
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=259'>260</a> duration = 30, # Duration of each window (in seconds)
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=260'>261</a> min_count = 10, # Number of window appearances to be counted as bad channel
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=261'>262</a> verbose = False,
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=262'>263</a> )
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=263'>264</a> if show_timings:
<a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=264'>265</a> print(f\"\\tPerformed in {dt.now() - t_step}\")
File <decorator-gen-433>:10, in find_bad_channels_maxwell(raw, limit, duration, min_count, return_scores, origin, int_order, ext_order, calibration, cross_talk, coord_frame, regularize, ignore_ref, bad_condition, head_pos, mag_scale, skip_by_annotation, h_freq, extended_proj, verbose)
File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:2723, in find_bad_channels_maxwell(***failed resolving arguments***)
2721 delta = chunk_raw.get_data(these_picks)
2722 with use_log_level(False):
-> 2723 _run_maxwell_filter(chunk_raw, reconstruct=\"orig\", copy=False, **params)
2725 if n_iter == 1 and len(chunk_flats):
2726 logger.info(
2727 \" Flat (%2d): %s\",
2728 len(chunk_flats),
2729 \" \".join(chunk_flats),
2730 )
File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:685, in _run_maxwell_filter(raw, skip_by_annotation, st_duration, st_correlation, st_only, st_when, ctc, coil_scale, this_pos_quat, meg_picks, good_mask, grad_picks, head_pos, info, _get_this_decomp_trans, S_recon, update_kwargs, ignore_ref, reconstruct, copy)
657 def _run_maxwell_filter(
658 raw,
659 skip_by_annotation,
(...)
683 # The time it takes to recompute S and pS themselves is roughly on par
684 # with the np.dot with the data, so not a huge gain to be made there.
--> 685 S_decomp, S_decomp_full, pS_decomp, reg_moments, n_use_in = _get_this_decomp_trans(
686 info[\"dev_head_t\"], t=0.0
687 )
688 update_kwargs.update(reg_moments=reg_moments.copy())
689 if ctc is not None:
File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:1240, in _get_decomp(trans, all_coils, cal, regularize, exp, ignore_ref, coil_scale, grad_picks, mag_picks, good_mask, mag_or_fine, bad_condition, t, mag_scale, mult)
1235 del extended_proj
1237 #
1238 # Regularization
1239 #
-> 1240 S_decomp, reg_moments, n_use_in = _regularize(
1241 regularize, exp, S_decomp, mag_or_fine, extended_remove, t=t
1242 )
1243 S_decomp_full = S_decomp_full.take(reg_moments, axis=1)
1245 #
1246 # Pseudo-inverse of total multipolar moment basis set (Part of Eq. 37)
1247 #
File <decorator-gen-431>:12, in _regularize(regularize, exp, S_decomp, mag_or_fine, extended_remove, t, verbose)
File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:1297, in _regularize(regularize, exp, S_decomp, mag_or_fine, extended_remove, t, verbose)
1295 t_str = f\"{t:8.3f}\"
1296 if regularize is not None: # regularize='in'
-> 1297 in_removes, out_removes = _regularize_in(
1298 int_order, ext_order, S_decomp, mag_or_fine, extended_remove
1299 )
1300 else:
1301 in_removes = []
File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:2294, in _regularize_in(int_order, ext_order, S_decomp, mag_or_fine, extended_remove)
2292 for ii in range(n_in):
2293 this_S = S_decomp.take(in_keepers + out_keepers, axis=1)
-> 2294 u, s, v = _safe_svd(this_S, full_matrices=False, **check_disable)
2295 del this_S
2296 eigs[ii] = s[[0, -1]]
File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/fixes.py:92, in _safe_svd(A, **kwargs)
90 raise ValueError(\"Cannot set overwrite_a=True with this function\")
91 try:
---> 92 return linalg.svd(A, **kwargs)
93 except np.linalg.LinAlgError as exp:
94 from .utils import warn
File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/scipy/linalg/_decomp_svd.py:166, in svd(a, full_matrices, compute_uv, overwrite_a, check_finite, lapack_driver)
164 raise LinAlgError(\"SVD did not converge\")
165 if info < 0:
--> 166 raise ValueError('illegal value in %dth argument of internal gesdd'
167 % -info)
168 if compute_uv:
169 return u, s, v
ValueError: illegal value in 4th argument of internal gesdd
Now, I checked and it seems like using raw_bids_empty_room.set_montage(None) turned the gradiometers and magnetometers positions to nans. As an example, an excerpt of the output of raw_bids_empty_room.info["chs"] is:
[{'scanno': 1, 'logno': 113, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 3.250000046861601e-09, 'coil_type': 3012 (FIFFV_COIL_VV_PLANAR_T1), 'loc': array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]), 'unit': 201 (FIFF_UNIT_T_M), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0113', 'coord_frame': 1 (FIFFV_COORD_DEVICE)},
{'scanno': 2, 'logno': 112, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 3.250000046861601e-09, 'coil_type': 3012 (FIFFV_COIL_VV_PLANAR_T1), 'loc': array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]), 'unit': 201 (FIFF_UNIT_T_M), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0112', 'coord_frame': 1 (FIFFV_COORD_DEVICE)},
{'scanno': 3, 'logno': 111, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 4.1400001482561066e-11, 'coil_type': 3022 (FIFFV_COIL_VV_MAG_T1), 'loc': array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]), 'unit': 112 (FIFF_UNIT_T), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0111', 'coord_frame': 1 (FIFFV_COORD_DEVICE)},
I checked, and printing raw_bids_empty_room.info["chs"] before calling raw_bids_empty_room.set_montage(None) proves that the loc data was there to begin with:
[{'scanno': 1, 'logno': 113, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 3.250000046861601e-09, 'coil_type': 3012 (FIFFV_COIL_VV_PLANAR_T1), 'loc': array([-0.1066 , 0.0464 , -0.0604 , -0.0127 , 0.0057 ,
-0.99990302, -0.186801 , -0.98240298, -0.0033 , -0.98232698,
0.18674099, 0.013541 ]), 'unit': 201 (FIFF_UNIT_T_M), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0113', 'coord_frame': 1 (FIFFV_COORD_DEVICE)},
{'scanno': 2, 'logno': 112, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 3.250000046861601e-09, 'coil_type': 3012 (FIFFV_COIL_VV_PLANAR_T1), 'loc': array([-0.1066 , 0.0464 , -0.0604 , -0.186801 , -0.98240298,
-0.0033 , 0.0127 , -0.0057 , 0.99990302, -0.98232698,
0.18674099, 0.013541 ]), 'unit': 201 (FIFF_UNIT_T_M), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0112', 'coord_frame': 1 (FIFFV_COORD_DEVICE)},
{'scanno': 3, 'logno': 111, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 4.1400001482561066e-11, 'coil_type': 3022 (FIFFV_COIL_VV_MAG_T1), 'loc': array([-0.1066 , 0.0464 , -0.0604 , -0.0127 , 0.0057 ,
-0.99990302, -0.186801 , -0.98240298, -0.0033 , -0.98232698,
0.18674099, 0.013541 ]), 'unit': 112 (FIFF_UNIT_T), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0111', 'coord_frame': 1 (FIFFV_COORD_DEVICE)},
So, I tried to do something else to remove the digitization data:
montage = empty_room_fif.get_montage() # We get the current montage
montage.remove_fiducials() # We remove the fiducials
montage.dig = [] # We remove the digitization points
# At this point, the digitization is made of 0 points
empty_room_fif.set_montage(montage) # We set the new montage
This worked, and allowed to remove the digitization points and coil data without removing the MEG channels positions.
I would like to know if this solution has a risk of breaking anything down the line - considering I am directly modifying an attribute instead of using a function, there might be some things I do not know about that I overlooked.
I compared the info with that of an empty room that didn’t have any digitization points to begin with. There are a few things that vary, notably dig is None in this one, instead of [] in the one we manually set (and I cannot set dig on None directly, as set_montage() then returns an error).
I am using Python 3.13.1 and MNE 1.9.0, working in a notebook on a Linux system.
Thank you!