Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Combine gradiometer channels #11082

Draft
wants to merge 3 commits into
base: main
Choose a base branch
from

Conversation

AJQuinn
Copy link
Contributor

@AJQuinn AJQuinn commented Aug 24, 2022

This is a work-in-progress start on general functions to combine meg gradiometer channels.

Idea is to allow users to combine gradiometer channels in the standard data instances, and eventually to replace the channel combining code used in mne/viz with this more general solution.

Would be great to get some feedback etc - any/all opinions welcome! @drammock @larsoner @agramfort

Some points for next steps...

  • modify coil_type and kind to something indicating a combined gradiometer
  • modify channel units, depending on combination method
  • currently _combine_meg_grads assumes inst only contains 'grads'
  • outsource return instance creation in _combine_meg_grads to specific class methods
  • need to handle noise_cov? projections?
  • addtional combination methods? eg PCA/SVD

Reference issue

Relevant: #1280

What does this implement/fix?

I've implemented a private function in mne.channels.layout and linked it up to epochs.as_type as a first use-case

Copy link
Member

@larsoner larsoner left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some preliminary comments/ideas!

def _combine_meg_grads(inst, method='rms'):
"""Combine pairs of gradiometer channels."""
# TODO:
# - modify coil_type and kind to a combined gradiometer
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

kind should still be MEG I think. For coil_type I think we'd maybe want to add a 3016 here:

https://github.com/mne-tools/mne-python/blob/main/mne/io/constants.py#L929

We'll need to make a PR to fif-constants for this. We can do this later, but to make mne/io/tests/test_constants.py happy you can for fif-constants, make a new branch with a commit that adds it, push to your fork, and tell test_constants to use your commit like here:

https://github.com/mne-tools/mne-python/pull/11064/files#diff-d545e84fd98e16f320b83c21e8c35a67ed64d378ded8fb35d3734e9f2196b53eR23-R25

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks - I think this is done in 5fa25d1 and 1e87305 referencing this fork of fiff-constants https://github.com/AJQuinn/fiff-constants/tree/grad_combined

Not sure how critical the naming choice is?

"""Combine pairs of gradiometer channels."""
# TODO:
# - modify coil_type and kind to a combined gradiometer
# - modify units depending on method
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think the units should change, right?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Makes sense for mean and rms - I'll ignore this for now.

# TODO:
# - modify coil_type and kind to a combined gradiometer
# - modify units depending on method
# - need to handle noise_cov? projections?
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes -- both should raise an error as you shouldn't need these. Same if you try to make_forward_solution... but we can tackle these separately because there will be a lot of these "don't do this" paths

# - modify coil_type and kind to a combined gradiometer
# - modify units depending on method
# - need to handle noise_cov? projections?
# - currently assumes inst only contains 'grads'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To start you could do inst.copy().pick(...), but in theory this should already live in as_type in some reasonable way

# - modify units depending on method
# - need to handle noise_cov? projections?
# - currently assumes inst only contains 'grads'
# - outsource return instance creation to specific methods
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should also already live in as_type. It would be good to reuse the code, deduplicating / making shared private functions as necessary

@AJQuinn
Copy link
Contributor Author

AJQuinn commented Aug 24, 2022

Ok - have some more updates which can be explored with this code.

import mne

sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = (sample_data_folder / 'MEG' / 'sample' /
                        'sample_audvis_filt-0-40_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file)
events = mne.find_events(raw, stim_channel='STI 014')

event_dict = {'auditory/left': 1, 'auditory/right': 2, 'visual/left': 3,
              'visual/right': 4, 'smiley': 5, 'buttonpress': 32}
reject_criteria = dict(grad=4000e-13)

epochs = mne.Epochs(raw, events, event_id=event_dict, tmin=-0.2, tmax=0.5,
                    reject=reject_criteria, preload=True, picks='grad')
cmb_epochs = epochs.as_type(ch_type='combined_grad', mode='rms')

epochs['visual'].average().plot_joint()
cmb_epochs['visual'].average().plot_joint()

Producing this standard ERF....
image

and this combined grads ERF
image

@AJQuinn
Copy link
Contributor Author

AJQuinn commented Aug 24, 2022

and a bit cleaner if you remember to baseline correct....

cmb_epochs['visual'].average().apply_baseline().plot_joint()

image

Copy link
Member

@agramfort agramfort left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

just a quick suggestion

picks = np.array(picks).reshape(-1, 2) # now [npairs x 2]
new_chs = list()
new_bads = list()
for idx in range(picks.shape[0]):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you can do:

Suggested change
for idx in range(picks.shape[0]):
for pair in picks:

then do below

Comment on lines +1017 to +1018
ch1 = inst.info['chs'][picks[idx, 0]]
ch2 = inst.info['chs'][picks[idx, 1]]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
ch1 = inst.info['chs'][picks[idx, 0]]
ch2 = inst.info['chs'][picks[idx, 1]]
ch1 = inst.info['chs'][pair[0]]
ch2 = inst.info['chs'][pair[1]]

@drammock
Copy link
Member

and a bit cleaner if you remember to baseline correct....

something is off about the plot. Since RMS makes the signal all positive, it should be plotting with a one-sided colormap. In our plotting code this will be handled by a variable called norm usually.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants