- MNE version: 1.8.0
- operating system: Windows 11
I’m trying to simulate raw data with self defined OPM sensor configuration. Following this post and this post I created a coil_def.dat file containing 64 sensors like this:
# coil_def.dat
1 9999 2 1 1e-02 0.000e+00 "OPM sensor. Amplitude of Bg: 5.000e-05"
1.000 -1.196e-02 -1.096e-02 1.290e-01 0.000e+00 0.000e+00 1.000e+00
1 9999 2 1 1e-02 0.000e+00 "OPM sensor. Amplitude of Bg: 5.000e-05"
1.000 2.446e-03 2.787e-02 1.270e-01 0.000e+00 0.000e+00 1.000e+00
...
The sensor class, id and accuracy are set as 1, 9999, 2 which represent
magnetomter, OPM and accurate representation, respectively (according to mne doc).
After reading the file and extracting sensor_number and sensor_locs, I created my info,
info = mne.create_info(ch_names=sensor_number,
sfreq=100,
ch_types="mag"
)
# info['dev_head_t'] = mne.transforms.Transform('meg', 'head', np.eye(4))
trans = mne.transforms.Transform("head","mri")
for i in range(sensor_number):
info["chs"][i]["loc"] = sensor_locs[i] # sensor position and oritation
info["chs"][i]["coil_type"] = 9999 # It seems this type represent a custom sensor and can be used for OPM.
To verify I checked info['chs'], it’s a list of length 64 and every item has a structure as below
{'loc': array([-0.01195882, -0.01095526, 0.12898438, 1. , 0. ,
0. , 0. , 1. , 0. , 0. ,
0. , 1. ]),
'unit_mul': 0 (FIFF_UNITM_NONE),
'range': 1.0,
'cal': 1.0,
'kind': 1 (FIFFV_MEG_CH),
'coil_type': 9999,
'unit': 112 (FIFF_UNIT_T),
'coord_frame': 1 (FIFFV_COORD_DEVICE),
'ch_name': '0',
'scanno': 1,
'logno': 1}
A sample anatomy from mne.datasets was chosen for simplicity.
data_path = mne.datasets.opm.data_path()
subject = "OPM_sample"
subjects_dir = data_path / "subjects"
raw_fname = data_path / "MEG" / "OPM" / "OPM_SEF_raw.fif"
bem_fname = subjects_dir / subject / "bem" / f"{subject}-5120-5120-5120-bem-sol.fif"
src = mne.setup_source_space(subject,add_dist=False,subjects_dir=subjects_dir)
bem = mne.read_bem_solution(bem_fname)
Then I tried to compute forward model with coil_def.dat and plot alignment,
with mne.use_coil_def(opm_coil_def_fname):
fwd = mne.make_forward_solution(info, trans, src, bem, eeg=False)
fwd = mne.convert_forward_solution(fwd, force_fixed=True, surf_ori=True, use_cps=True)
fig = mne.viz.plot_alignment(
info,
trans=trans,
subject=subject,
subjects_dir=subjects_dir,
# dig=True,
coord_frame="mri",
surfaces=("head", "white"),
bem=bem,
)
# kwargs = dict(azimuth=0, elevation=90, distance=0.6, focalpoint=(0.0, 0.0, 0.0))
# mne.viz.set_3d_view(figure=fig, **kwargs)
And I got a QhullError
---------------------------------------------------------------------------
QhullError Traceback (most recent call last)
File d:\Anaconda\envs\mne250306\lib\site-packages\mne\viz\_3d.py:1738, in _sensor_shape(coil)
1737 try:
-> 1738 tris = _reorder_ccw(rrs, ConvexHull(rrs).simplices)
1739 except QhullError: # 2D geometry likely
File _qhull.pyx:2425, in scipy.spatial._qhull.ConvexHull.__init__()
File _qhull.pyx:343, in scipy.spatial._qhull._Qhull.__init__()
QhullError: QH6214 qhull input error: not enough points(1) to construct initial simplex (need 4)
While executing: | qhull i Qt
Options selected for Qhull 2019.1.r 2019/06/21:
run-id 1078245579 incidence Qtriangulate _pre-merge _zero-centrum
_maxoutside 0
During handling of the above exception, another exception occurred:
IndexError Traceback (most recent call last)
Cell In[8], line 4
2 fwd = mne.make_forward_solution(info, trans, src, bem, eeg=False)
3 fwd = mne.convert_forward_solution(fwd, force_fixed=True, surf_ori=True, use_cps=True)
----> 4 fig = mne.viz.plot_alignment(
5 info,
6 trans=trans,
7 subject=subject,
8 subjects_dir=subjects_dir,
9 # dig=True,
10 coord_frame="mri",
11 surfaces=("head", "white"),
12 bem=bem,
13 )
14 # kwargs = dict(azimuth=0, elevation=90, distance=0.6, focalpoint=(0.0, 0.0, 0.0))
15 # mne.viz.set_3d_view(figure=fig, **kwargs)
File <decorator-gen-185>:12, in plot_alignment(info, trans, subject, subjects_dir, surfaces, coord_frame, meg, eeg, fwd, dig, ecog, src, mri_fiducials, bem, seeg, fnirs, show_axes, dbs, fig, interaction, sensor_colors, verbose)
File d:\Anaconda\envs\mne250306\lib\site-packages\mne\viz\_3d.py:896, in plot_alignment(info, trans, subject, subjects_dir, surfaces, coord_frame, meg, eeg, fwd, dig, ecog, src, mri_fiducials, bem, seeg, fnirs, show_axes, dbs, fig, interaction, sensor_colors, verbose)
893 # plot sensors (NB snapshot_brain_montage relies on the last thing being
894 # plotted being the sensors, so we need to do this after the surfaces)
895 if picks.size > 0:
--> 896 _plot_sensors_3d(
897 renderer,
898 info,
899 to_cf_t,
900 picks,
901 meg,
902 eeg,
903 fnirs,
904 warn_meg,
905 head_surf,
906 "m",
907 sensor_alpha=sensor_alpha,
908 sensor_colors=sensor_colors,
909 )
911 if src is not None:
912 atlas_ids, colors = read_freesurfer_lut()
File d:\Anaconda\envs\mne250306\lib\site-packages\mne\viz\_3d.py:1488, in _plot_sensors_3d(renderer, info, to_cf_t, picks, meg, eeg, fnirs, warn_meg, head_surf, units, sensor_alpha, orient_glyphs, scale_by_distance, project_points, surf, check_inside, nearest, sensor_colors)
1485 from matplotlib.colors import to_rgba_array
1487 defaults = DEFAULTS["coreg"]
-> 1488 ch_pos, sources, detectors = _ch_pos_in_coord_frame(
1489 pick_info(info, picks), to_cf_t=to_cf_t, warn_meg=warn_meg
1490 )
1492 actors = defaultdict(lambda: list())
1493 locs = defaultdict(lambda: list())
File <decorator-gen-186>:12, in _ch_pos_in_coord_frame(info, to_cf_t, warn_meg, verbose)
File d:\Anaconda\envs\mne250306\lib\site-packages\mne\viz\_3d.py:1066, in _ch_pos_in_coord_frame(info, to_cf_t, warn_meg, verbose)
1062 coil = _create_meg_coils(this_coil, acc="normal", coilset=coilset)[
1063 0
1064 ]
1065 # store verts as ch_coord
-> 1066 ch_coord, triangles = _sensor_shape(coil)
1067 ch_coord = apply_trans(coil_trans, ch_coord)
1068 if len(ch_coord) == 0 and warn_meg:
File d:\Anaconda\envs\mne250306\lib\site-packages\mne\viz\_3d.py:1742, in _sensor_shape(coil)
1740 logger.debug("Falling back to planar geometry")
1741 u, _, _ = np.linalg.svd(rrs.T, full_matrices=False)
-> 1742 u[:, 2] = 0
1743 rr_rot = rrs @ u
1744 tris = Delaunay(rr_rot[:, :2]).simplices
IndexError: index 2 is out of bounds for axis 1 with size 1
I navigated to _3d.py but couldn’t settle the problem.
I’ll be appreciated if anybody could tell me what caused the bug or how to fix it. I’m not sure whether I defined my coil_def.dat correctly (for example the class, id, accuracy, weight, etc.)
Thanks for any advices.