Calling maxwell_filter on Yokogawa recordings from Brainstorm

💬 Support & Discussions
1 
Platform:         Windows-10-10.0.22621-SP0
Python:           3.10.10 (tags/v3.10.10:aad5f6a, Feb  7 2023, 17:20:36) [MSC v.1929 64 bit (AMD64)]
mne:              1.3.1

Problem initially posted on the Brainstorm forum:

I could reproduce the error on my end with an example dataset provided by the user.
I confirmed that the problem is not coming too obviously from Brainstorm in the creation of the RawArray object to be passed to maxwell_filter, as I can all other MNE methods on it.

Example dataset to reproduce the error (zipped pickle - first 10s of the Yokogawa .con file shared by the user):

Code:

import pickle 
import mne
filehandler = open('pyRaw10s.pickle', 'rb') 
raw = pickle.load(filehandler)
filehandler.close()
mne.preprocessing.maxwell_filter(raw, int_order=8, ext_order=3, origin='auto', coord_frame='head', regularize='in', ignore_ref=False, st_duration=10, st_correlation=0.98, st_fixed=True, st_only=False, mag_scale=100, skip_by_annotation=['edge','bad_acq_skip'])

Output (same error as from Brainstorm):

>>> mne.preprocessing.maxwell_filter(raw, int_order=8, ext_order=3, origin='auto', coord_frame='head', regularize='in', ignore_ref=False, st_duration=10, st_correlation=0.98, st_fixed=True, st_only=False, mag_scale=100, skip_by_annotation=['edge','bad_acq_skip'])
Maxwell filtering raw data
    No bad MEG channels
    Processing 0 gradiometers and 128 magnetometers (of which 125 are actually KIT gradiometers)
<stdin>:1: RuntimeWarning: Only 5 head digitization points of the specified kind ("extra",), fitting may be inaccurate
<stdin>:1: RuntimeWarning: (X, Y) fit (-7.9, 42.7) more than 20 mm from head frame origin
    Automatic origin fit: head of radius 73.8 mm
    Using origin -7.9, 42.7, 5.3 mm in the head frame
    Processing data using tSSS with st_duration=10.0
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<decorator-gen-457>", line 12, in maxwell_filter
  File "C:\Users\franc\AppData\Local\Programs\Python\Python310\lib\site-packages\mne\preprocessing\maxwell.py", line 353, in maxwell_filter
    raw_sss = _run_maxwell_filter(raw, **params)
  File "C:\Users\franc\AppData\Local\Programs\Python\Python310\lib\site-packages\mne\preprocessing\maxwell.py", line 555, in _run_maxwell_filter
    _get_this_decomp_trans(info['dev_head_t'], t=0.)
  File "C:\Users\franc\AppData\Local\Programs\Python\Python310\lib\site-packages\mne\preprocessing\maxwell.py", line 1012, in _get_decomp
    S_decomp, reg_moments, n_use_in = _regularize(
  File "<decorator-gen-460>", line 12, in _regularize
  File "C:\Users\franc\AppData\Local\Programs\Python\Python310\lib\site-packages\mne\preprocessing\maxwell.py", line 1066, in _regularize
    in_removes, out_removes = _regularize_in(
  File "C:\Users\franc\AppData\Local\Programs\Python\Python310\lib\site-packages\mne\preprocessing\maxwell.py", line 1947, in _regularize_in
    u, s, v = _safe_svd(this_S, full_matrices=False, **check_disable)
  File "C:\Users\franc\AppData\Local\Programs\Python\Python310\lib\site-packages\mne\fixes.py", line 84, in _safe_svd
    return linalg.svd(A, **kwargs)
  File "C:\Users\franc\AppData\Local\Programs\Python\Python310\lib\site-packages\scipy\linalg\_decomp_svd.py", line 133, in svd
    raise ValueError('illegal value in %dth argument of internal gesdd'
ValueError: illegal value in 4th argument of internal gesdd

Is the problem in the creation of the RawArray object in Brainstorm, or in MNE?

2

Hi Francois -

I don’t know if this is the correct way to do it, but it does run. If you drop the reference sensors (In [2] below) - the maxfilter runs. With Elekta/MEGIN not having external reference channels - I suspect it wasn’t designed to handle these – thats just a guess though.

In [1]: import pickle
   ...: import mne
   ...: filehandler = open('pyRaw10s.pickle', 'rb')
   ...: raw = pickle.load(filehandler)
   ...: filehandler.close()

In [2]: raw.pick_types(meg=True, ref_meg=False)
Out[2]: <RawArray | 125 x 50001 (10.0 s), ~47.8 MB, data loaded>

In [3]: mne.preprocessing.maxwell_filter(raw, int_order=8, ext_order=3, origin='auto', coord_frame='head', regularize='in', ignore_ref=False, st_dura
   ...: tion=10, st_correlation=0.98, st_fixed=True, st_only=False, mag_scale=100, skip_by_annotation=['edge','bad_acq_skip'])
Maxwell filtering raw data
    No bad MEG channels
    Processing 0 gradiometers and 125 magnetometers (of which 125 are actually KIT gradiometers)
<ipython-input-3-283f85684eea>:1: RuntimeWarning: Only 5 head digitization points of the specified kind ("extra",), fitting may be inaccurate
  mne.preprocessing.maxwell_filter(raw, int_order=8, ext_order=3, origin='auto', coord_frame='head', regularize='in', ignore_ref=False, st_duration=10, st_correlation=0.98, st_fixed=True, st_only=False, mag_scale=100, skip_by_annotation=['edge','bad_acq_skip'])
<ipython-input-3-283f85684eea>:1: RuntimeWarning: (X, Y) fit (-7.9, 42.7) more than 20 mm from head frame origin
  mne.preprocessing.maxwell_filter(raw, int_order=8, ext_order=3, origin='auto', coord_frame='head', regularize='in', ignore_ref=False, st_duration=10, st_correlation=0.98, st_fixed=True, st_only=False, mag_scale=100, skip_by_annotation=['edge','bad_acq_skip'])
    Automatic origin fit: head of radius 73.8 mm
    Using origin -7.9, 42.7, 5.3 mm in the head frame
    Processing data using tSSS with st_duration=10.0
        Using 53/95 harmonic components for    0.000  (41/80 in, 12/15 out)
    Using loaded raw data
    Processing 1 data chunk
        Projecting  3 intersecting tSSS components for    0.000 -   10.000 sec (#1/1)
[done]
Out[3]: <RawArray | 125 x 50001 (10.0 s), ~47.8 MB, data loaded>

Also - I would verify the origin fit - since it only found 5 points (I am assuming those are localizer coils). You can set those to (0, 0, 0.04) - versus fitting to the headshape.

–Jeff

3

Hi Jeff,

Thank you for this information.
However, it is not a nice thing to drop the reference sensors, as they are used in the estimation of the forward model, later in the analysis.

Note that we can run tSSS on CTF recordings, which also include reference sensors.

Any other suggestion?

4

Should I report this as a github issue?

5

Francois,

I am not too familiar with the backend of that algorithm, but I was able to reproduce your error when keeping the ref sensors.

@larsoner @agramfort @richard @drammock – thoughts on the issue?

6

@Francois yes feel free to open a GitHub issue for this. It would be good to include a link to the original KIT data. That way we can try it on read_raw_kit-ed data as well as the RawArray that BrainStorm creates, and see if they match. I’m guessing it’s some problem with the encoding of the reference channel positions, but we’ll see !

7

Done: Calling maxwell_filter on Yokogawa recordings (from Brainstorm) · Issue #11548 · mne-tools/mne-python · GitHub