How to realign the CTF MEG runs to a common head position using maxwell_filter?

💬 Support & Discussions
,
1

Hi, MNE experts

Recently I would like to concatenate 12 different MEG runs. It reports an error that I could not concatenate different runs because of different head positions. I could use mne.transforms.Transform to change the head position (epoch.info[‘dev_head_t’] = Transform(‘meg’, ‘head’, np.identity(4))) but the following source localization would not work because of the wrong head position. I notice that the ‘Movement compensation’ function could annotate movement artifacts (Signal-space separation (SSS) and Maxwell filtering — MNE 1.0.2 documentation), but CTF MEG data do not have an additional ‘.pos’ file,
So I would like to ask how to realign the CTF MEG runs to a common head position?

Environment

  • MNE version: e.g. 0.24.0
  • operating system: Ubuntu 18.04
  • data: CTF MEG data.

Relavant code:

for file in os.listdir(savePath):
    if 'filterEpochICAMR_' in file:
        fifpath = pj(savePath, file)
        epoch = mne.read_epochs(fifpath, preload=False, verbose=True)
        epochs_list.append(epoch)
        del epoch
epochs_all = mne.concatenate_epochs(epochs_list)

ERROR:
Exception has occurred: ValueError
epochs[1].info['dev_head_t'] differs. The instances probably come from different runs, and are therefore associated with different head positions. Manually change info['dev_head_t'] to avoid this message but beware that this means the MEG sensors will not be properly spatially aligned. See mne.preprocessing.maxwell_filter to realign the runs to a common head position.
  File "evoked.py", line xx, in <module>
    epochs_all = mne.concatenate_epochs(epochs_list)

Best!
Jinhua

2

Hello @Clemens and welcome to the forum!

Personally I’ve never really worked with CTF data, but generally the approach should be the following:

  • Decide which of your acquisition runs should be the “reference” run. The goal is to transform the head position of all runs to that of the reference run.
  • Extract the dev_head_t from the reference run, which in this example was assigned to the variable raw_ref:
    dev_head_t_ref = raw_ref.info['dev_head_t']
  • Pass this transformation matrix, together with an origin and coord_frame='head' parameter, to maxwell_filter(), e.g.:
    raw_sss = mne.preprocessing.maxwell_filter(
    raw,
    origin='auo',
    coord_frame='head',
    destination=dev_head_t_ref
)

Do this for all your runs; after that, the head positions should be “equalized” across runs. Don’t forget to Maxwell filter the reference run too :slight_smile: You can use the exact same parameters there.

Best wishes,
Richard

2 Likes
3

Hi Mr. Richard,

Thank you so much for your reply, it’s really helpful!

However, an error was reported when I use this code.

ERROR:
Exception has occurred: RuntimeError
Maxwell filter cannot be done on compensated channels, but data have been compensated with grade 3.

In the ‘Brainstorm CTF phantom dataset tutorial’ (Brainstorm CTF phantom dataset tutorial — MNE 1.0.2 documentation), it says that must un-do software compensation first (raw.apply_gradient_compensation(0)). but after that, another error was reported.

Code:
raw.apply_gradient_compensation(0)
raw = mne.preprocessing.maxwell_filter(raw, origin=‘auto’,coord_frame=‘head’, estination=dev_head_t_ref)

ERROR:
Exception has occurred: ValueError
Only 0 head digitization points of the specified kinds (“eeg”, “extra”,), at least 4 required
File “xx.py”, line xx, in
raw = mne.preprocessing.maxwell_filter(raw, rigin=‘auto’,coord_frame=‘head’,destination=dev_head_t_ref)

And when I changed the code as raw.apply_gradient_compensation(3) OR aw.apply_gradient_compensation(0, 0,0) OR raw.apply_gradient_compensation(0, 0, 0, 0) , it still not works.

ERROR:
Maxwell filter cannot be done on compensated channels, but data have been compensated with grade 3.

Another question is which run could I use as the ‘reference run’? Could I just use the first run or average all 12 runs’ [‘dev_head_t’] information together to compute an averaged head position?

Best!
Jinhua

4

hi,

to me both are fine but you can see experimentally how big of a difference it makes

Alex

1 Like
5

In the MNE-BIDS-Pipeline, we take the first run by default; I’ve also worked with researchers who would usually take the “middle” run, e.g. the fourth out of an 8-runs experiment, as the reference.

1 Like
6

@Clemens - Does your dataset have a headshape?

I think the error is related to an auxilliary function - fit_sphere_to_headshape - and if there is no headshape, you get the “Only 0 head digitization points 
”

Here is the error I get with a CTF dataset and no headshape:

mne.preprocessing.maxwell_filter(raw)
Maxwell filtering raw data
    No bad MEG channels
    Processing 0 gradiometers and 299 magnetometers (of which 290 are actually KIT gradiometers)
Traceback (most recent call last):

  File "<ipython-input-35-980e5af9fd01>", line 1, in <module>
    mne.preprocessing.maxwell_filter(raw)

  File "<decorator-gen-410>", line 24, in maxwell_filter

  File "/home/jstout/miniconda3/envs/test5/lib/python3.7/site-packages/mne/preprocessing/maxwell.py", line 214, in maxwell_filter
    skip_by_annotation=skip_by_annotation, extended_proj=extended_proj)

  File "<decorator-gen-411>", line 24, in _prep_maxwell_filter

  File "/home/jstout/miniconda3/envs/test5/lib/python3.7/site-packages/mne/preprocessing/maxwell.py", line 312, in _prep_maxwell_filter
    origin = _check_origin(origin, info, coord_frame, disp=True)

  File "/home/jstout/miniconda3/envs/test5/lib/python3.7/site-packages/mne/bem.py", line 996, in _check_origin
    units='m')[:2]

  File "<decorator-gen-60>", line 22, in fit_sphere_to_headshape

  File "/home/jstout/miniconda3/envs/test5/lib/python3.7/site-packages/mne/bem.py", line 850, in fit_sphere_to_headshape
    info, dig_kinds)

  File "<decorator-gen-62>", line 24, in _fit_sphere_to_headshape

  File "/home/jstout/miniconda3/envs/test5/lib/python3.7/site-packages/mne/bem.py", line 933, in _fit_sphere_to_headshape
    hsp = get_fitting_dig(info, dig_kinds)

  File "<decorator-gen-61>", line 24, in get_fitting_dig

  File "/home/jstout/miniconda3/envs/test5/lib/python3.7/site-packages/mne/bem.py", line 896, in get_fitting_dig
    return get_fitting_dig(info, ('extra', 'eeg'))

  File "<decorator-gen-61>", line 24, in get_fitting_dig

  File "/home/jstout/miniconda3/envs/test5/lib/python3.7/site-packages/mne/bem.py", line 924, in get_fitting_dig
    raise ValueError(msg + ', at least 4 required')

ValueError: Only 2 head digitization points of the specified kinds ("eeg", "extra",), at least 4 required

This is in the documentation:

origin : array-like, shape (3,) | str
    Origin of internal and external multipolar moment space in meters.
    The default is ``'auto'``, which means ``(0., 0., 0.)`` when
    ``coord_frame='meg'``, and a head-digitization-based
    origin fit using :func:`~mne.bem.fit_sphere_to_headshape`
    when ``coord_frame='head'``. If automatic fitting fails (e.g., due
    to having too few digitization points),
    consider separately calling the fitting function with different
    options or specifying the origin manually.

If you manually assign the origin - it works - I don’t know what an appropriate origin would be, but 0,0,0 sounds like it should be good.:

raw_sss = mne.preprocessing.maxwell_filter(
    raw,
    origin=(0,0,0),
    coord_frame='head',
    destination=dev_head_t_ref
)
Maxwell filtering raw data
    No bad MEG channels
    Processing 0 gradiometers and 299 magnetometers (of which 290 are actually KIT gradiometers)
    Using origin 0.0, 0.0, 0.0 mm in the head frame
        Using 63/95 harmonic components for    0.000  (48/80 in, 15/15 out)
    Using loaded raw data
    Processing 45 data chunks
[done]
1 Like
7

Dear MNE experts,

Thank you so much for Mr. Höchenberger and Mr. Gramfort’s reply!

When I use the code as follows, it still does not work:

Environment

  • MNE version: e.g. 0.24.0
  • operating system: Ubuntu 18.04
  • data: CTF MEG data.

Relavant code:
raw = read_raw_ctf(filepath, preload=True)
dev_head_t_ref = raw.info[‘dev_head_t’]
raw.apply_gradient_compensation(0) # un-do software compensation
raw = mne.preprocessing.maxwell_filter(raw, origin=‘auto’,coord_frame=‘head’, destination=dev_head_t_ref)

ERROR:
Exception has occurred: ValueError
Only 0 head digitization points of the specified kinds (“eeg”, “extra”,), at least 4 required
File “xx.py”, line xx, in
raw = mne.preprocessing.maxwell_filter(raw, origin=‘auto’,coord_frame=‘head’, destination=dev_head_t_ref)

When I set the maxwell_filtering parameter “ignore_ref=True”(If True, do not include reference channels in compensation.), it still reports the same error.

One data I used in data processing was provided in google drive: CTFRawData.zip - Google Drive
Thank you for your attention! Hope that experts who are familiar with CTF data could help solve this problem:)

Best wishes!
Jinhua

8

Dear Mr. Stout,

Thank you for your reply! I tried your code and it works!
I used it as follows:

raw.apply_gradient_compensation(0)
raw = mne.preprocessing.maxwell_filter(raw, origin=(0., 0., 0.), coord_frame=‘head’, destination=dev_head_t_ref)

Do I need to undo the software compensation before the maxwell filter?

Best wishes!
Jinhua Tian

1 Like
9

@Clemens - I believe there is a post somewhere about SSS needing the raw data without 3rd order gradient applied - I can’t find it though.

@larsoner , @richard or @agramfort may be able to comment.

1 Like
10

Dear Mr. Stout,

Thank you for your reply!

I’m also confused about if it is appropriate to set the ‘origin’ as “(0., 0., 0.)”. And if I want to set the origin using “mne.bem.fit_sphere_to_headshape(raw.info)” to get the original head center position, it reports the same error. Hope that someone could answer it:)

Relavant code:
raw = read_raw_ctf(filepath, preload=True)
,head_origin, = mne.bem.fit_sphere_to_headshape(raw.info)
dev_head_t_ref = raw.info[‘dev_head_t’]
raw.apply_gradient_compensation(0)
raw = mne.preprocessing.maxwell_filter(raw, origin=head_origin, coord_frame=‘head’, destination=dev_head_t_ref)

Error:
Exception has occurred: ValueError
Only 0 head digitization points of the specified kinds (“eeg”, “extra”,), at least 4 required
File “xx.py”, line xx, in
,head_origin, = mne.bem.fit_sphere_to_headshape(raw.info)

I’d try origin as ‘(0.,0.,0.)’ to see if it works well:)

Best!
Jinhua Tian

11

I think you’re right @jstout211 – you could always try it and see if you get an error.

I think in principle it should be possible to allow Maxwell filtering with gradient compensation applied since it’s just a linear operator – we should be able to apply it to the SSS basis functions and have everything work I think. But again this would need to be implemented and tested


I’m also confused about if it is appropriate to set the ‘origin’ as “(0., 0., 0.)”. And if I want to set the origin using “mne.bem.fit_sphere_to_headshape(raw.info)” to get the original head center position, it reports the same error. Hope that someone could answer it:)

You need digitization points available in order to use the sphere-fitting (‘auto’) mode. (0., 0., 0.) is probably is not an ideal origin for a real human head. You could go with something like (0., 0., 0.04), which IIRC is the default in MEGIN’s maxfilter software and generally reasonable for an adult head. A smaller Z offset is usually good for children.

1 Like
12

@Clemens here is a quick hack using MNE internal functions that could help you

# %%
import mne

raw = mne.io.read_raw_ctf('S01_G15BNU_20220405_01.ds')

# %%

info_from = raw.copy().pick('meg').info
info_to = raw.copy().pick('meg').info
map = mne.forward._map_meg_or_eeg_channels(
    info_from, info_to, mode='fast', origin=(0., 0., 0.04))

map.shape

# %%
data = raw.get_data(picks='meg')
data_new_pos = map @ data

big disclaimer it’s not tested


but it could be a start of a PR to make such a feature part of the public API without relying on maxfilter

Alex

3 Likes
13

Hi Mr. Gramfort,

Thank you for your suggestion! I will try this method:)

Another naive question is, do these two methods (maxwell filter and dot product) have a theoretical difference?

Jinhua.

14

I think it would be great to have this as part of the main API.

Additionally, could this also emulate continuous head tracking? Would it make sense to apply the mapping to each data point individually, given the current head position?

15

to do continuous head position correction you would do this on time segments as if it gets
applied for each sample you’ll mess up the time courses I fear. But I never tried


Alex