From 372b18e624411f163872a454afe5c6644ad61957 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Mon, 17 Nov 2014 17:12:42 +0200 Subject: [PATCH 01/12] Add more referencing options This adds a base function 'apply_reference' that can apply various different referencing schemes. Two convenience functions are also added that provide a convenient interface to the base function: - set_eeg_reference: apply a simple reference to all EEG channels - set_bipolar_reference: apply a bipolar reference (or more) --- mne/io/__init__.py | 3 +- mne/io/base.py | 266 +++++++++++++++++++++++++++++----- mne/io/fiff/tests/test_raw.py | 82 ++++++++++- 3 files changed, 313 insertions(+), 38 deletions(-) diff --git a/mne/io/__init__.py b/mne/io/__init__.py index 0d74a626a65..5c4f4fb99db 100644 --- a/mne/io/__init__.py +++ b/mne/io/__init__.py @@ -30,4 +30,5 @@ # for backward compatibility from .fiff import RawFIFF from .fiff import RawFIFF as Raw -from .base import concatenate_raws, get_chpi_positions, set_eeg_reference +from .base import (concatenate_raws, get_chpi_positions, apply_reference, + set_eeg_reference, set_bipolar_reference) diff --git a/mne/io/base.py b/mne/io/base.py index 885343dea37..b1d671300de 100644 --- a/mne/io/base.py +++ b/mne/io/base.py @@ -397,7 +397,8 @@ def filter(self, l_freq, h_freq, picks=None, filter_length='10s', if l_freq < h_freq: logger.info('Band-pass filtering from %0.2g - %0.2g Hz' % (l_freq, h_freq)) - self._data = band_pass_filter(self._data, fs, l_freq, h_freq, + self._data = band_pass_filter( + self._data, fs, l_freq, h_freq, filter_length=filter_length, l_trans_bandwidth=l_trans_bandwidth, h_trans_bandwidth=h_trans_bandwidth, @@ -406,7 +407,8 @@ def filter(self, l_freq, h_freq, picks=None, filter_length='10s', else: logger.info('Band-stop filtering from %0.2g - %0.2g Hz' % (h_freq, l_freq)) - self._data = band_stop_filter(self._data, fs, h_freq, l_freq, + self._data = band_stop_filter( + self._data, fs, h_freq, l_freq, filter_length=filter_length, l_trans_bandwidth=h_trans_bandwidth, h_trans_bandwidth=l_trans_bandwidth, method=method, @@ -741,7 +743,7 @@ def save(self, fname, picks=None, tmin=0, tmax=None, buffer_size_sec=10, int=FIFF.FIFFT_INT, single=FIFF.FIFFT_FLOAT, double=FIFF.FIFFT_DOUBLE) - if not format in type_dict.keys(): + if format not in type_dict.keys(): raise ValueError('format must be "short", "int", "single", ' 'or "double"') reset_dict = dict(short=False, int=False, single=True, double=True) @@ -1330,19 +1332,20 @@ def add_events(self, events, stim_channel=None): self._data[pick, idx - self.first_samp] += events[:, 2] -def set_eeg_reference(raw, ref_channels, copy=True): - """Rereference eeg channels to new reference channel(s). - - If multiple reference channels are specified, they will be averaged. +def apply_reference(raw, ref_from, ref_to=None, copy=True): + """ + Calculate a reference signal by taking the mean of a set of channels and + apply the reference to another set of channels. Parameters ---------- raw : instance of Raw - Instance of Raw with eeg channels and reference channel(s). - - ref_channels : list of str - The name(s) of the reference channel(s). - + Instance of Raw with EEG channels and reference channel(s). + ref_from : list of str + The names of the channels to use to construct the reference. + ref_to : list of str | None + The names of the channels to apply the reference to. By default, + all EEG channels are chosen. copy : bool Specifies whether instance of Raw will be copied or modified in place. @@ -1350,40 +1353,230 @@ def set_eeg_reference(raw, ref_channels, copy=True): ------- raw : instance of Raw Instance of Raw with eeg channels rereferenced. - ref_data : array - Array of reference data subtracted from eeg channels. + Array of reference data subtracted from EEG channels. + + Notes + ----- + 1. Do not use this function to apply a common average. By default, an + average reference projection has already been added upon loading raw + data. + + 2. If the reference is applied to any EEG channels, this function removes + any pre-existing average reference projections. + + 3. During source localization, the EEG signal should have a common average + reference. + + 4. The data must be preloaded. + + See also + -------- + set_eeg_reference : Convenience function for creating an EEG reference. + set_bipolar_reference : Convenience function for creating a bipolar + reference. """ # Check to see that raw data is preloaded if not raw.preload: raise RuntimeError('Raw data needs to be preloaded. Use ' 'preload=True (or string) in the constructor.') - # Make sure that reference channels are loaded as list of string - if not isinstance(ref_channels, list): - raise IOError('Reference channel(s) must be a list of string. ' - 'If using a single reference channel, enter as ' - 'a list with one element.') - # Find the indices to the reference electrodes - ref_idx = [raw.ch_names.index(c) for c in ref_channels] - - # Get the reference array - ref_data = raw._data[ref_idx].mean(0) - - # Get the indices to the eeg channels using the pick_types function - eeg_idx = pick_types(raw.info, exclude="bads", eeg=True, meg=False, - ref_meg=False) - - # Copy raw data or modify raw data in place - if copy: # copy data + + eeg_idx = pick_types(raw.info, eeg=True, meg=False, ref_meg=False) + + ref_from = [raw.ch_names.index(ch) for ch in ref_from] + if ref_to is None: + ref_to = eeg_idx + else: + ref_to = [raw.ch_names.index(ch) for ch in ref_to] + + # Compute reference + ref_data = raw._data[ref_from].mean(0) + + if copy: raw = raw.copy() - # Rereference the eeg channels - raw._data[eeg_idx] -= ref_data + raw._data[ref_to] -= ref_data + + # If the reference touches EEG electrodes, remove any pre-existing common + # reference and note in the info that a non-CAR has been applied. + if len(np.intersect1d(ref_from, eeg_idx)) > 0: + raw.info['custom_reference'] = True + + # Remove any existing average reference projections + for i, proj in enumerate(raw.info['projs']): + if proj['desc'] == 'Average EEG reference' or \ + proj['kind'] == FIFF.FIFFV_MNE_PROJ_ITEM_EEG_AVREF: + logger.info('Removing existing average EEG reference ' + 'projection.') + del raw.info['projs'][i] - # Return rereferenced data and reference array return raw, ref_data +def set_eeg_reference(raw, ref_channels=None, copy=True): + """Rereference EEG channels to new reference channel(s). + + If multiple reference channels are specified, they will be averaged. + + Parameters + ---------- + raw : instance of Raw + Instance of Raw with EEG channels and reference channel(s). + ref_channels : list of str | None + The names of the channels to use to construct the reference. If None is + specified here, a common average reference will be applied in the form + of an SSP projector. + copy : bool + Specifies whether instance of Raw will be copied (True) or modified in + place (False). Defaults to True. + + Returns + ------- + raw : instance of Raw + Instance of Raw with eeg channels rereferenced. + ref_data : array + Array of reference data subtracted from EEG channels. + + Notes + ----- + 1. If a reference is requested that is not the common average reference, + this function removes any pre-existing average reference projections. + + 2. During source localization, the EEG signal should have a common average + reference. + + 3. In order to apply a reference other than a common average reference, the + data must be preloaded. + + See also + -------- + apply_reference : Base function for rereferencing. + set_bipolar_reference : Convenience function for creating bipolar + references. + """ + if ref_channels is None: + # CAR requested + if _has_eeg_average_ref_proj(raw.info['projs']): + logger.warning('A common average reference projection was already ' + 'added. The data has been left untouched.') + return raw, None + else: + raw.add_proj(make_eeg_average_ref_proj(raw.info), activate=False) + return raw, None + else: + logger.info('Applying a custom EEG reference.') + return apply_reference(raw, ref_channels, copy=copy) + + +def set_bipolar_reference(raw, anode, cathode, ch_name, ch_info=None, + copy=True): + """Rereference selected channels using a bipolar referencing scheme. + + A bipolar reference takes the difference between two channels (the anode + minus the cathode) and adds it as a new virtual channel. The original + channels will be dropped. + + Multiple anodes and cathodes can be specified, in which case multiple + vitual channels will be created. The 1st anode will be substracted from the + 1st cathode, the 2nd anode from the 2nd cathode, etc. + + By default, the virtual channels will be annotated with channel info of + the anodes and its location set to (0, 0, 0)) and coil type set to + EEG_BIPOLAR. + + Parameters + ---------- + raw : instance of Raw + Instance of Raw containing the unreferenced channels. + anode : str | list of str + The name(s) of the channel(s) to use as anode in the bipolar reference. + cathode : str | list of str + The name(s) of the channel(s) to use as cathode in the bipolar + reference. + ch_name : str | list of str + The channel name(s) for the virtual channel(s) containing the resulting + signal. + ch_info : dict | list of dict | None + This parameter can be used to supply a dictionary (or a dictionary for + each bipolar channel) containing channel information to merge in, + overwriting the default values. + copy : bool + Whether to operate on a copy of the raw data (True) or modify it + in-place (False). Defaults to True. + + Returns + ------- + raw : instance of Raw + Instance of Raw with the specified channels rereferenced. + + See also + -------- + apply_reference : Base function for rereferencing. + set_eeg_reference : Convenience function for creating an EEG reference. + + Notes + ----- + 1. If the anodes contain any EEG channels, this function removes + any pre-existing average reference projections. + + 2. During source localization, the EEG signal should have a common average + reference. + + 3. The data must be preloaded. + """ + if not isinstance(anode, list): + anode = [anode] + + if not isinstance(cathode, list): + cathode = [cathode] + + assert len(anode) == len(cathode), ( + 'Number of anodes must equal the number of cathodes.') + + if not isinstance(ch_name, list): + ch_name = [ch_name] + assert len(ch_name) == len(anode), ( + 'Number of channel names must equal the number of anodes/cathodes.') + + if ch_info is None: + ch_info = [{} for an in anode] + elif not isinstance(ch_info, list): + ch_info = [ch_info] + assert len(ch_info) == len(anode), ( + 'Number of channel info dictionaries must equal the number of ' + 'anodes/cathodes.') + + # Merge specified and anode channel information dictionaries + new_ch_info = [] + for an, ci in zip(anode, ch_info): + new_info = raw.info['chs'][raw.ch_names.index(an)].copy() + + # Set channel location and coil type + if 'eeg_loc' in new_info: + new_info['eeg_loc'] = np.zeros((3, 2)) + new_info['loc'] = np.zeros(12) + new_info['coil_type'] = FIFF.FIFFV_COIL_EEG_BIPOLAR + + new_info.update(ci) + new_ch_info.append(new_info) + + if copy: + raw = raw.copy() + + # Perform bipolar referencing + for an, ca, name, info in zip(anode, cathode, ch_name, new_ch_info): + raw, _ = apply_reference(raw, [ca], [an], copy=False) + an_idx = raw.ch_names.index(an) + raw.info['chs'][an_idx] = info + raw.info['chs'][an_idx]['ch_name'] = name + raw.info['ch_names'][an_idx] = name + + # Drop cathode channels + raw.drop_channels(cathode) + + return raw + + def _allocate_data(data, data_buffer, data_shape, dtype): if data is None: # if not already done, allocate array with right type @@ -1526,7 +1719,8 @@ def _write_raw(fname, raw, info, picks, format, data_type, reset_range, start, # Split files if necessary, leave some space for next file info if pos >= split_size - this_buff_size_bytes - 2 ** 20: - next_fname, next_idx = _write_raw(fname, raw, info, picks, format, + next_fname, next_idx = _write_raw( + fname, raw, info, picks, format, data_type, reset_range, first + buffer_size, stop, buffer_size, projector, inv_comp, drop_small_buffer, split_size, part_idx + 1, use_fname) @@ -1646,7 +1840,7 @@ def _write_raw_buffer(fid, buf, cals, format, inv_comp): if buf.shape[0] != len(cals): raise ValueError('buffer and calibration sizes do not match') - if not format in ['short', 'int', 'single', 'double']: + if format not in ['short', 'int', 'single', 'double']: raise ValueError('format must be "short", "single", or "double"') if np.isrealobj(buf): diff --git a/mne/io/fiff/tests/test_raw.py b/mne/io/fiff/tests/test_raw.py index c67fe043d8d..8d61c39c0ed 100644 --- a/mne/io/fiff/tests/test_raw.py +++ b/mne/io/fiff/tests/test_raw.py @@ -21,7 +21,8 @@ from mne.datasets import testing from mne.io.constants import FIFF from mne.io import (Raw, concatenate_raws, - get_chpi_positions, set_eeg_reference) + get_chpi_positions, apply_reference, set_eeg_reference, + set_bipolar_reference) from mne import concatenate_events, find_events, equalize_channels from mne.utils import (_TempDir, requires_nitime, requires_pandas, requires_mne, run_subprocess, run_tests_if_main) @@ -982,6 +983,40 @@ def compensate_mne(fname, grad): assert_allclose(raw_py._data, raw_c._data, rtol=1e-6, atol=1e-17) +@testing.requires_testing_data +def test_apply_reference(): + """ Test base function for rereferencing""" + raw = Raw(fif_fname, preload=True) + + # Separate EEG channels from other channel types + picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') + picks_other = pick_types(raw.info, meg=True, eeg=False, eog=True, + stim=True, exclude='bads') + + # Rereference raw data by creating a copy of original data + reref, ref_data = apply_reference(raw, ref_from=['EEG 001', 'EEG 002'], + copy=True) + + # Get the raw EEG data and other channel data + raw_eeg_data = raw[picks_eeg][0] + raw_other_data = raw[picks_other][0] + + # Get the rereferenced EEG data and channel other + reref_eeg_data = reref[picks_eeg][0] + unref_eeg_data = reref_eeg_data + ref_data + # Undo rereferencing of EEG channels + reref_other_data = reref[picks_other][0] + + # Check that both EEG data and other data is the same + assert_allclose(raw_eeg_data, unref_eeg_data, 1e-6, atol=1e-15) + assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) + + # Test that data is modified in place when copy=False + reref, ref_data = apply_reference(raw, ['EEG 001', 'EEG 002'], + copy=False) + assert_true(raw is reref) + + @testing.requires_testing_data def test_set_eeg_reference(): """ Test rereference eeg data""" @@ -1015,6 +1050,51 @@ def test_set_eeg_reference(): assert_true(raw is reref) +@testing.requires_testing_data +def test_set_bipolar_reference(): + """ Test bipolar referencing""" + raw = Raw(fif_fname, preload=True) + reref = set_bipolar_reference(raw, 'EEG 001', 'EEG 002', 'bipolar', + {'kind': FIFF.FIFFV_EOG_CH, + 'extra': 'some extra value'}) + + # Compare result to a manual calculation + a = raw.pick_channels(['EEG 001', 'EEG 002'], copy=True) + a = a._data[0, :] - a._data[1, :] + b = reref.pick_channels(['bipolar'], copy=True)._data[0, :] + assert_allclose(a, b) + + # Original channels should be replaced by a virtual one + assert_true('EEG 001' not in reref.ch_names) + assert_true('EEG 002' not in reref.ch_names) + assert_true('bipolar' in reref.ch_names) + + # Check channel information + bp_info = reref.info['chs'][reref.ch_names.index('bipolar')] + an_info = reref.info['chs'][raw.ch_names.index('EEG 001')] + for key in bp_info: + if key == 'loc' or key == 'eeg_loc': + assert_array_equal(bp_info[key], 0) + elif key == 'coil_type': + assert_equal(bp_info[key], FIFF.FIFFV_COIL_EEG_BIPOLAR) + elif key == 'kind': + assert_equal(bp_info[key], FIFF.FIFFV_EOG_CH) + else: + assert_equal(bp_info[key], an_info[key]) + assert_equal(bp_info['extra'], 'some extra value') + + # Test a battery of invalid inputs + assert_raises(AssertionError, set_bipolar_reference, raw, + 'EEG 001', ['EEG 002', 'EEG 003'], 'bipolar') + assert_raises(AssertionError, set_bipolar_reference, raw, + ['EEG 001', 'EEG 002'], 'EEG 003', 'bipolar') + assert_raises(AssertionError, set_bipolar_reference, raw, + 'EEG 001', 'EEG 002', ['bipolar1', 'bipolar2']) + assert_raises(AssertionError, set_bipolar_reference, raw, + 'EEG 001', 'EEG 002', 'bipolar', + [{'foo': 'bar'}, {'foo': 'bar'}]) + + @testing.requires_testing_data def test_drop_channels_mixin(): """Test channels-dropping functionality From 82ebc6e4e6fe6528436f6002179471ccdc78a046 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Tue, 2 Dec 2014 10:26:05 +0200 Subject: [PATCH 02/12] Add flag that a custom reference has been applied Raw and Epoch classes will no longer add an average reference if the `custom_ref_applied` flag was set. SSP projections will still attempt this, but a new assertion in `make_eeg_average_ref_proj` will raise an error if the flag was set. Before a CAR projection is added, the `info['custom_ref_applied']` flag is checked. --- mne/forward/forward.py | 5 +++ mne/io/base.py | 67 ++++++++++++++----------------- mne/io/brainvision/brainvision.py | 1 + mne/io/bti/bti.py | 29 ++++++------- mne/io/constants.py | 1 + mne/io/edf/edf.py | 2 +- mne/io/egi/egi.py | 1 + mne/io/fiff/raw.py | 13 ++++-- mne/io/fiff/tests/test_raw.py | 23 +++++++++-- mne/io/kit/kit.py | 1 + mne/io/meas_info.py | 17 ++++++-- mne/io/proj.py | 39 ++++++++++++------ mne/tests/test_epochs.py | 6 +++ mne/tests/test_proj.py | 22 +++++++++- mne/viz/raw.py | 2 +- 15 files changed, 152 insertions(+), 77 deletions(-) diff --git a/mne/forward/forward.py b/mne/forward/forward.py index 68e6d409d32..81566160f7b 100644 --- a/mne/forward/forward.py +++ b/mne/forward/forward.py @@ -346,6 +346,11 @@ def _read_forward_meas_info(tree, fid): raise ValueError('MEG/head coordinate transformation not found') info['bads'] = read_bad_channels(fid, parent_meg) + + # Check if a custom reference has been applied + tag = find_tag(fid, parent_mri, FIFF.FIFF_CUSTOM_REF) + info['custom_ref_applied'] = bool(tag.data) if tag is not None else False + return info diff --git a/mne/io/base.py b/mne/io/base.py index b1d671300de..7054b3ce858 100644 --- a/mne/io/base.py +++ b/mne/io/base.py @@ -77,15 +77,6 @@ def __hash__(self): raise RuntimeError('Cannot hash raw unless preloaded') return object_hash(dict(info=self.info, data=self._data)) - def _add_eeg_ref(self, add_eeg_ref): - """Helper to add an average EEG reference""" - if add_eeg_ref: - eegs = pick_types(self.info, meg=False, eeg=True, ref_meg=False) - projs = self.info['projs'] - if len(eegs) > 0 and not _has_eeg_average_ref_proj(projs): - eeg_ref = make_eeg_average_ref_proj(self.info, activate=False) - projs.append(eeg_ref) - def _parse_get_set_params(self, item): # make sure item is a tuple if not isinstance(item, tuple): # only channel selection passed @@ -1333,16 +1324,17 @@ def add_events(self, events, stim_channel=None): def apply_reference(raw, ref_from, ref_to=None, copy=True): - """ - Calculate a reference signal by taking the mean of a set of channels and - apply the reference to another set of channels. + """Apply a custom EEG referencing scheme. + + Calculates a reference signal by taking the mean of a set of channels and + applies the reference to another set of channels. Parameters ---------- raw : instance of Raw Instance of Raw with EEG channels and reference channel(s). ref_from : list of str - The names of the channels to use to construct the reference. + The names of the channels to use to construct the reference. ref_to : list of str | None The names of the channels to apply the reference to. By default, all EEG channels are chosen. @@ -1352,20 +1344,20 @@ def apply_reference(raw, ref_from, ref_to=None, copy=True): Returns ------- raw : instance of Raw - Instance of Raw with eeg channels rereferenced. + Instance of Raw with EEG channels rereferenced. ref_data : array Array of reference data subtracted from EEG channels. Notes ----- - 1. Do not use this function to apply a common average. By default, an + 1. Do not use this function to apply an average reference. By default, an average reference projection has already been added upon loading raw data. 2. If the reference is applied to any EEG channels, this function removes any pre-existing average reference projections. - 3. During source localization, the EEG signal should have a common average + 3. During source localization, the EEG signal should have an average reference. 4. The data must be preloaded. @@ -1399,8 +1391,8 @@ def apply_reference(raw, ref_from, ref_to=None, copy=True): # If the reference touches EEG electrodes, remove any pre-existing common # reference and note in the info that a non-CAR has been applied. - if len(np.intersect1d(ref_from, eeg_idx)) > 0: - raw.info['custom_reference'] = True + if len(np.intersect1d(ref_to, eeg_idx)) > 0: + raw.info['custom_ref_applied'] = True # Remove any existing average reference projections for i, proj in enumerate(raw.info['projs']): @@ -1424,8 +1416,8 @@ def set_eeg_reference(raw, ref_channels=None, copy=True): Instance of Raw with EEG channels and reference channel(s). ref_channels : list of str | None The names of the channels to use to construct the reference. If None is - specified here, a common average reference will be applied in the form - of an SSP projector. + specified here, an average reference will be applied in the form of an + SSP projector. copy : bool Specifies whether instance of Raw will be copied (True) or modified in place (False). Defaults to True. @@ -1439,14 +1431,14 @@ def set_eeg_reference(raw, ref_channels=None, copy=True): Notes ----- - 1. If a reference is requested that is not the common average reference, - this function removes any pre-existing average reference projections. + 1. If a reference is requested that is not the average reference, this + function removes any pre-existing average reference projections. - 2. During source localization, the EEG signal should have a common average + 2. During source localization, the EEG signal should have an average reference. - 3. In order to apply a reference other than a common average reference, the - data must be preloaded. + 3. In order to apply a reference other than an average reference, the data + must be preloaded. See also -------- @@ -1457,7 +1449,7 @@ def set_eeg_reference(raw, ref_channels=None, copy=True): if ref_channels is None: # CAR requested if _has_eeg_average_ref_proj(raw.info['projs']): - logger.warning('A common average reference projection was already ' + logger.warning('An average reference projection was already ' 'added. The data has been left untouched.') return raw, None else: @@ -1519,7 +1511,7 @@ def set_bipolar_reference(raw, anode, cathode, ch_name, ch_info=None, 1. If the anodes contain any EEG channels, this function removes any pre-existing average reference projections. - 2. During source localization, the EEG signal should have a common average + 2. During source localization, the EEG signal should have an average reference. 3. The data must be preloaded. @@ -1530,21 +1522,22 @@ def set_bipolar_reference(raw, anode, cathode, ch_name, ch_info=None, if not isinstance(cathode, list): cathode = [cathode] - assert len(anode) == len(cathode), ( - 'Number of anodes must equal the number of cathodes.') + if len(anode) != len(cathode): + raise ValueError('Number of anodes must equal the number of cathodes.') if not isinstance(ch_name, list): ch_name = [ch_name] - assert len(ch_name) == len(anode), ( - 'Number of channel names must equal the number of anodes/cathodes.') + if len(ch_name) != len(anode): + raise ValueError('Number of channel names must equal the number of ' + 'anodes/cathodes.') if ch_info is None: ch_info = [{} for an in anode] elif not isinstance(ch_info, list): ch_info = [ch_info] - assert len(ch_info) == len(anode), ( - 'Number of channel info dictionaries must equal the number of ' - 'anodes/cathodes.') + if len(ch_info) != len(anode): + raise ValueError('Number of channel info dictionaries must equal the ' + 'number of anodes/cathodes.') # Merge specified and anode channel information dictionaries new_ch_info = [] @@ -2058,7 +2051,7 @@ def _check_update_montage(info, montage): # raise error if positions are missing if missing_positions: err = ("The following positions are missing from the montage " - "definitions: %s. If those channels lack positions " - "because they are EOG channels use the eog parameter." - % str(missing_positions)) + "definitions: %s. If those channels lack positions " + "because they are EOG channels use the eog parameter." + % str(missing_positions)) raise KeyError(err) diff --git a/mne/io/brainvision/brainvision.py b/mne/io/brainvision/brainvision.py index b1857e14355..1bf3a870440 100644 --- a/mne/io/brainvision/brainvision.py +++ b/mne/io/brainvision/brainvision.py @@ -436,6 +436,7 @@ def _get_eeg_info(vhdr_fname, reference, eog, misc): info['orig_blocks'] = None info['line_freq'] = None info['subject_info'] = None + info['custom_ref_applied'] = False eeg_info = {} diff --git a/mne/io/bti/bti.py b/mne/io/bti/bti.py index df528939167..7b1db0bf2b6 100644 --- a/mne/io/bti/bti.py +++ b/mne/io/bti/bti.py @@ -781,22 +781,22 @@ def _read_bti_header(pdf_fname, config_fname): fid.seek(1, 1) info.update({'data_format': read_int16(fid), - 'acq_mode': read_int16(fid), - 'total_epochs': read_int32(fid), - 'input_epochs': read_int32(fid), - 'total_events': read_int32(fid), - 'total_fixed_events': read_int32(fid), - 'sample_period': read_float(fid), - 'xaxis_label': read_str(fid, 16), - 'total_processes': read_int32(fid), - 'total_chans': read_int16(fid)}) + 'acq_mode': read_int16(fid), + 'total_epochs': read_int32(fid), + 'input_epochs': read_int32(fid), + 'total_events': read_int32(fid), + 'total_fixed_events': read_int32(fid), + 'sample_period': read_float(fid), + 'xaxis_label': read_str(fid, 16), + 'total_processes': read_int32(fid), + 'total_chans': read_int16(fid)}) fid.seek(2, 1) info.update({'checksum': read_int32(fid), - 'total_ed_classes': read_int32(fid), - 'total_associated_files': read_int16(fid), - 'last_file_index': read_int16(fid), - 'timestamp': read_int32(fid)}) + 'total_ed_classes': read_int32(fid), + 'total_associated_files': read_int16(fid), + 'last_file_index': read_int16(fid), + 'timestamp': read_int32(fid)}) fid.seek(20, 1) _correct_offset(fid) @@ -988,7 +988,7 @@ def __init__(self, pdf_fname, config_fname='config', logger.info('Reading 4D PDF file %s...' % pdf_fname) bti_info = _read_bti_header(pdf_fname, config_fname) - # XXX indx is informed guess. Normally only one transform is stored. + # XXX indx is informed guess. Normally only one transform is stored. dev_ctf_t = bti_info['bti_transform'][0].astype('>f8') bti_to_nm = bti_to_vv_trans(adjust=rotation_x, translation=translation, dtype='>f8') @@ -1022,6 +1022,7 @@ def __init__(self, pdf_fname, config_fname='config', info['lowpass'] = lp info['acq_pars'], info['acq_stim'] = None, None info['filename'] = None + info['custom_ref_applied'] = False chs = [] ch_names = [ch['name'] for ch in bti_info['chs']] diff --git a/mne/io/constants.py b/mne/io/constants.py index caad8b7170c..224df3d9377 100644 --- a/mne/io/constants.py +++ b/mne/io/constants.py @@ -124,6 +124,7 @@ def __init__(self, **kwargs): FIFF.FIFF_DESCRIPTION = FIFF.FIFF_COMMENT # (Textual) Description of an object FIFF.FIFF_DIG_STRING = 234 # String of digitized points FIFF.FIFF_LINE_FREQ = 235 # Line frequency +FIFF.FIFF_CUSTOM_REF = 236 # Whether a custom reference was applied to the data # # HPI fitting program tags # diff --git a/mne/io/edf/edf.py b/mne/io/edf/edf.py index a3e535b691a..6862b64de2a 100644 --- a/mne/io/edf/edf.py +++ b/mne/io/edf/edf.py @@ -13,7 +13,6 @@ import re import warnings from math import ceil, floor -import warnings import numpy as np from scipy.interpolate import interp1d @@ -422,6 +421,7 @@ def _get_edf_info(fname, stim_channel, annot, annotmap, tal_channel, info['experimenter'] = None info['line_freq'] = None info['subject_info'] = None + info['custom_ref_applied'] = False edf_info = dict() edf_info['annot'] = annot diff --git a/mne/io/egi/egi.py b/mne/io/egi/egi.py index 1936b22ea82..de9bbcbe2bb 100644 --- a/mne/io/egi/egi.py +++ b/mne/io/egi/egi.py @@ -282,6 +282,7 @@ def __init__(self, input_fname, montage=None, eog=None, misc=None, info['ch_names'] = ch_names info['bads'] = [] info['comps'] = [] + info['custom_ref_applied'] = False for ii, ch_name in enumerate(ch_names): ch_info = {'cal': cal, 'logno': ii + 1, diff --git a/mne/io/fiff/raw.py b/mne/io/fiff/raw.py index 0ed334860fb..c1479628f08 100644 --- a/mne/io/fiff/raw.py +++ b/mne/io/fiff/raw.py @@ -17,7 +17,8 @@ from ..meas_info import read_meas_info from ..tree import dir_tree_find from ..tag import read_tag -from ..proj import proj_equal +from ..proj import (proj_equal, make_eeg_average_ref_proj, + _needs_eeg_average_ref_proj) from ..compensator import get_current_comp, set_current_comp, make_compensator from ..base import _BaseRaw @@ -127,7 +128,10 @@ def __init__(self, fnames, allow_maxshield=False, preload=False, self.verbose = verbose self.orig_format = raws[0].orig_format self.proj = False - self._add_eeg_ref(add_eeg_ref) + + if add_eeg_ref and _needs_eeg_average_ref_proj(self.info): + eeg_ref = make_eeg_average_ref_proj(self.info, activate=False) + self.add_proj(eeg_ref) if preload: self._preload_data(preload) @@ -305,8 +309,9 @@ def _read_raw_file(self, fname, allow_maxshield, preload, compensation, idx2 = base.rfind('-') if idx2 < 0 and next_num == 1: # this is the first file, which may not be numbered - next_fname = op.join(path, '%s-%d.%s' % (base[:idx], - next_num, base[idx + 1:])) + next_fname = op.join( + path, '%s-%d.%s' % (base[:idx], next_num, + base[idx + 1:])) continue num_str = base[idx2 + 1:idx] if not num_str.isdigit(): diff --git a/mne/io/fiff/tests/test_raw.py b/mne/io/fiff/tests/test_raw.py index 8d61c39c0ed..caafb456551 100644 --- a/mne/io/fiff/tests/test_raw.py +++ b/mne/io/fiff/tests/test_raw.py @@ -23,6 +23,7 @@ from mne.io import (Raw, concatenate_raws, get_chpi_positions, apply_reference, set_eeg_reference, set_bipolar_reference) +from mne.io.proj import _has_eeg_average_ref_proj from mne import concatenate_events, find_events, equalize_channels from mne.utils import (_TempDir, requires_nitime, requires_pandas, requires_mne, run_subprocess, run_tests_if_main) @@ -986,7 +987,7 @@ def compensate_mne(fname, grad): @testing.requires_testing_data def test_apply_reference(): """ Test base function for rereferencing""" - raw = Raw(fif_fname, preload=True) + raw = Raw(fif_fname, preload=True, add_eeg_ref=True) # Separate EEG channels from other channel types picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') @@ -996,6 +997,14 @@ def test_apply_reference(): # Rereference raw data by creating a copy of original data reref, ref_data = apply_reference(raw, ref_from=['EEG 001', 'EEG 002'], copy=True) + assert_true(reref.info['custom_ref_applied']) + + # The CAR reference projection should have been removed by the function + assert_true(not _has_eeg_average_ref_proj(reref.info['projs'])) + + # Check that the ref has been properly computed + ref_ch_idx = [raw.ch_names.index(ch) for ch in ['EEG 001', 'EEG 002']] + assert_array_equal(ref_data, raw[ref_ch_idx, :][0].mean(0)) # Get the raw EEG data and other channel data raw_eeg_data = raw[picks_eeg][0] @@ -1003,8 +1012,8 @@ def test_apply_reference(): # Get the rereferenced EEG data and channel other reref_eeg_data = reref[picks_eeg][0] - unref_eeg_data = reref_eeg_data + ref_data # Undo rereferencing of EEG channels + unref_eeg_data = reref_eeg_data + ref_data reref_other_data = reref[picks_other][0] # Check that both EEG data and other data is the same @@ -1024,6 +1033,7 @@ def test_set_eeg_reference(): # Rereference raw data by creating a copy of original data reref, ref_data = set_eeg_reference(raw, ['EEG 001', 'EEG 002'], copy=True) + assert_true(reref.info['custom_ref_applied']) # Separate EEG channels from other channel types picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') @@ -1036,8 +1046,8 @@ def test_set_eeg_reference(): # Get the rereferenced EEG data and channel other reref_eeg_data = reref[picks_eeg][0] - unref_eeg_data = reref_eeg_data + ref_data # Undo rereferencing of EEG channels + unref_eeg_data = reref_eeg_data + ref_data reref_other_data = reref[picks_other][0] # Check that both EEG data and other data is the same @@ -1057,6 +1067,7 @@ def test_set_bipolar_reference(): reref = set_bipolar_reference(raw, 'EEG 001', 'EEG 002', 'bipolar', {'kind': FIFF.FIFFV_EOG_CH, 'extra': 'some extra value'}) + assert_true(reref.info['custom_ref_applied']) # Compare result to a manual calculation a = raw.pick_channels(['EEG 001', 'EEG 002'], copy=True) @@ -1083,6 +1094,12 @@ def test_set_bipolar_reference(): assert_equal(bp_info[key], an_info[key]) assert_equal(bp_info['extra'], 'some extra value') + # Test creating a bipolar reference that doesn't involve EEG channels: + # it should not set the custom_ref_applied flag + reref = set_bipolar_reference(raw, 'MEG 0111', 'MEG 0112', 'bipolar', + {'kind': FIFF.FIFFV_MEG_CH}) + assert_true(not reref.info['custom_ref_applied']) + # Test a battery of invalid inputs assert_raises(AssertionError, set_bipolar_reference, raw, 'EEG 001', ['EEG 002', 'EEG 003'], 'bipolar') diff --git a/mne/io/kit/kit.py b/mne/io/kit/kit.py index 9adb75d1d99..1cf3f353747 100644 --- a/mne/io/kit/kit.py +++ b/mne/io/kit/kit.py @@ -116,6 +116,7 @@ def __init__(self, input_fname, mrk=None, elp=None, hsp=None, stim='>', self.cals = np.ones(self.info['nchan']) self.orig_format = 'double' self.rawdirs = list() + self.info['custom_ref_applied'] = False if isinstance(mrk, list): mrk = [read_mrk(marker) if isinstance(marker, string_types) diff --git a/mne/io/meas_info.py b/mne/io/meas_info.py index 15853711744..d2010d4855a 100644 --- a/mne/io/meas_info.py +++ b/mne/io/meas_info.py @@ -405,6 +405,7 @@ def read_meas_info(fid, tree, verbose=None): proj_id = None proj_name = None line_freq = None + custom_ref_applied = False p = 0 for k in range(meas_info['nent']): kind = meas_info['directory'][k].kind @@ -452,6 +453,9 @@ def read_meas_info(fid, tree, verbose=None): elif kind == FIFF.FIFF_LINE_FREQ: tag = read_tag(fid, pos) line_freq = float(tag.data) + elif kind == FIFF.FIFF_CUSTOM_REF: + tag = read_tag(fid, pos) + custom_ref_applied = bool(tag.data) # Check that we have everything we need if nchan is None: @@ -627,6 +631,7 @@ def read_meas_info(fid, tree, verbose=None): info['comps'] = comps info['acq_pars'] = acq_pars info['acq_stim'] = acq_stim + info['custom_ref_applied'] = custom_ref_applied return info, meas @@ -745,6 +750,8 @@ def write_meas_info(fid, info, data_type=None, reset_range=True): write_float(fid, FIFF.FIFF_LINE_FREQ, info['line_freq']) if data_type is not None: write_int(fid, FIFF.FIFF_DATA_PACK, data_type) + if info.get('custom_ref_applied'): + write_int(fid, FIFF.FIFF_CUSTOM_REF, info['custom_ref_applied']) # Channel information for k, c in enumerate(info['chs']): @@ -927,10 +934,11 @@ def _merge_info(infos, verbose=None): trans_name) raise ValueError(msg) other_fields = ['acq_pars', 'acq_stim', 'bads', 'buffer_size_sec', - 'comps', 'description', 'dig', 'experimenter', 'file_id', - 'filename', 'highpass', 'line_freq', 'lowpass', - 'meas_date', 'meas_id', 'orig_blocks', 'proj_id', - 'proj_name', 'projs', 'sfreq', 'subject_info', 'sfreq'] + 'comps', 'custom_ref_applied', 'description', 'dig', + 'experimenter', 'file_id', 'filename', 'highpass', + 'line_freq', 'lowpass', 'meas_date', 'meas_id', + 'orig_blocks', 'proj_id', 'proj_name', 'projs', 'sfreq', + 'subject_info', 'sfreq'] for k in other_fields: info[k] = _merge_dict_values(infos, k) @@ -997,4 +1005,5 @@ def create_info(ch_names, sfreq, ch_types=None): info['dev_head_t'] = None info['dev_ctf_t'] = None info['ctf_head_t'] = None + info['custom_ref_applied'] = False return info diff --git a/mne/io/proj.py b/mne/io/proj.py index a13f3a2418d..0eb67a6aef6 100644 --- a/mne/io/proj.py +++ b/mne/io/proj.py @@ -165,10 +165,10 @@ def del_proj(self, idx): self.info['projs'].pop(idx) return self - + def plot_projs_topomap(self, ch_type=None, layout=None): """Plot SSP vector - + Parameters ---------- ch_type : 'mag' | 'grad' | 'planar1' | 'planar2' | 'eeg' | None | List @@ -197,7 +197,7 @@ def plot_projs_topomap(self, ch_type=None, layout=None): if ch_type is None: ch_type = [ch for ch in ['meg', 'eeg'] if ch in self] elif isinstance(ch_type, str): - ch_type = [ch_type] + ch_type = [ch_type] for ch in ch_type: if ch in self: layout.append(find_layout(self.info, ch, exclude=[])) @@ -207,7 +207,7 @@ def plot_projs_topomap(self, ch_type=None, layout=None): fig = plot_projs_topomap(self.info['projs'], layout) else: raise ValueError("Info is missing projs. Nothing to plot.") - + return fig @@ -584,6 +584,11 @@ def make_eeg_average_ref_proj(info, activate=True, verbose=None): eeg_proj: instance of Projection The SSP/PCA projector. """ + if info['custom_ref_applied']: + raise RuntimeError('Cannot add an average EEG reference projection ' + 'since a custom reference has been applied to the ' + 'data earlier.') + logger.info("Adding average EEG reference projection.") eeg_sel = pick_types(info, meg=False, eeg=True, ref_meg=False, exclude='bads') @@ -605,12 +610,25 @@ def make_eeg_average_ref_proj(info, activate=True, verbose=None): def _has_eeg_average_ref_proj(projs): """Determine if a list of projectors has an average EEG ref""" for proj in projs: - if proj['desc'] == 'Average EEG reference' or \ - proj['kind'] == FIFF.FIFFV_MNE_PROJ_ITEM_EEG_AVREF: + if (proj['desc'] == 'Average EEG reference' or + proj['kind'] == FIFF.FIFFV_MNE_PROJ_ITEM_EEG_AVREF): return True return False +def _needs_eeg_average_ref_proj(info): + """Determine if the EEG needs an averge EEG reference + + This returns True if no custom reference has been applied and no average + reference projection is present in the list of projections. + """ + eeg_sel = pick_types(info, meg=False, eeg=True, ref_meg=False, + exclude='bads') + return (len(eeg_sel) > 0 and + not info['custom_ref_applied'] and + not _has_eeg_average_ref_proj(info['projs'])) + + @verbose def setup_proj(info, add_eeg_ref=True, activate=True, verbose=None): @@ -636,14 +654,11 @@ def setup_proj(info, add_eeg_ref=True, activate=True, The modified measurement info (Warning: info is modified inplace). """ # Add EEG ref reference proj if necessary - eeg_sel = pick_types(info, meg=False, eeg=True, ref_meg=False, - exclude='bads') - if len(eeg_sel) > 0 and not _has_eeg_average_ref_proj(info['projs']) \ - and add_eeg_ref is True: + if _needs_eeg_average_ref_proj(info) and add_eeg_ref: eeg_proj = make_eeg_average_ref_proj(info, activate=activate) info['projs'].append(eeg_proj) - # Create the projector + # Create the projector projector, nproj = make_projector_info(info) if nproj == 0: if verbose: @@ -654,7 +669,7 @@ def setup_proj(info, add_eeg_ref=True, activate=True, logger.info('Created an SSP operator (subspace dimension = %d)' % nproj) - # The projection items have been activated + # The projection items have been activated if activate: info['projs'] = activate_proj(info['projs'], copy=False) diff --git a/mne/tests/test_epochs.py b/mne/tests/test_epochs.py index 8115c65bdc8..c85ecd6aebc 100644 --- a/mne/tests/test_epochs.py +++ b/mne/tests/test_epochs.py @@ -276,6 +276,12 @@ def test_epochs_proj(): baseline=(None, 0), proj=True, add_eeg_ref=False) assert_true(not _has_eeg_average_ref_proj(epochs.info['projs'])) + # make sure we don't add avg ref when a custom ref has been applied + raw.info['custom_ref_applied'] = True + epochs = Epochs(raw, events[:4], event_id, tmin, tmax, picks=this_picks, + baseline=(None, 0), proj=True) + assert_true(not _has_eeg_average_ref_proj(epochs.info['projs'])) + def test_evoked_arithmetic(): """Test arithmetic of evoked data diff --git a/mne/tests/test_proj.py b/mne/tests/test_proj.py index b8148624d66..1e735634eac 100644 --- a/mne/tests/test_proj.py +++ b/mne/tests/test_proj.py @@ -1,5 +1,5 @@ import os.path as op -from nose.tools import assert_true +from nose.tools import assert_true, assert_raises import warnings import numpy as np @@ -217,4 +217,24 @@ def test_compute_proj_raw(): assert_array_almost_equal(proj, np.eye(len(raw.ch_names))) +def test_make_eeg_average_ref_proj(): + ''' Test EEG average reference projection ''' + raw = Raw(raw_fname, add_eeg_ref=False, preload=True) + eeg = mne.pick_types(raw.info, meg=False, eeg=True) + + # No average EEG reference + assert_true(not np.all(raw._data[eeg].mean(axis=0) < 1e-19)) + + # Apply average EEG reference + car = make_eeg_average_ref_proj(raw.info) + reref = raw.copy() + reref.add_proj(car) + reref.apply_proj() + assert_array_almost_equal(reref._data[eeg].mean(axis=0), 0, decimal=19) + + # Error when custom reference has already been applied + raw.info['custom_ref_applied'] = True + assert_raises(AssertionError, make_eeg_average_ref_proj, raw.info) + + run_tests_if_main() diff --git a/mne/viz/raw.py b/mne/viz/raw.py index 5b1f9b5210a..b2eb31944ac 100644 --- a/mne/viz/raw.py +++ b/mne/viz/raw.py @@ -401,7 +401,7 @@ def plot_raw(raw, events=None, duration=10.0, start=0.0, n_channels=None, else: if lowpass <= highpass: raise ValueError('lowpass (%s) must be > highpass (%s)' - %(lowpass, highpass)) + % (lowpass, highpass)) ba = butter(filtorder, [highpass / nyq, lowpass / nyq], 'bandpass', analog=False) From ff0b1ebb0a8a98ea5a16636150d4d55486c2c499 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Tue, 16 Dec 2014 10:05:47 +0200 Subject: [PATCH 03/12] Add default channel name for bipolar references - Also check for duplicate channel names - Bilpolar channels can be named after the anode/cathode - Inform user of the chosen bipolar channel name --- mne/io/base.py | 20 +++++++++++++++---- mne/io/fiff/tests/test_raw.py | 17 +++++++++------- mne/tests/test_proj.py | 37 +++++++++++++++++++++++++++++++---- 3 files changed, 59 insertions(+), 15 deletions(-) diff --git a/mne/io/base.py b/mne/io/base.py index 7054b3ce858..1596ccd3a3e 100644 --- a/mne/io/base.py +++ b/mne/io/base.py @@ -1460,7 +1460,7 @@ def set_eeg_reference(raw, ref_channels=None, copy=True): return apply_reference(raw, ref_channels, copy=copy) -def set_bipolar_reference(raw, anode, cathode, ch_name, ch_info=None, +def set_bipolar_reference(raw, anode, cathode, ch_name=None, ch_info=None, copy=True): """Rereference selected channels using a bipolar referencing scheme. @@ -1485,9 +1485,10 @@ def set_bipolar_reference(raw, anode, cathode, ch_name, ch_info=None, cathode : str | list of str The name(s) of the channel(s) to use as cathode in the bipolar reference. - ch_name : str | list of str + ch_name : str | list of str | None The channel name(s) for the virtual channel(s) containing the resulting - signal. + signal. By default, bipolar channels are named after the anode and + cathode, but it is recommended to supply a more meaningful name. ch_info : dict | list of dict | None This parameter can be used to supply a dictionary (or a dictionary for each bipolar channel) containing channel information to merge in, @@ -1525,12 +1526,22 @@ def set_bipolar_reference(raw, anode, cathode, ch_name, ch_info=None, if len(anode) != len(cathode): raise ValueError('Number of anodes must equal the number of cathodes.') - if not isinstance(ch_name, list): + if ch_name is None: + ch_name = ['%s-%s' % ac for ac in zip(anode, cathode)] + elif not isinstance(ch_name, list): ch_name = [ch_name] if len(ch_name) != len(anode): raise ValueError('Number of channel names must equal the number of ' 'anodes/cathodes.') + # Check for duplicate channel names (it is allowed to give the name of the + # anode or cathode channel, as they will be replaced). + for ch, a, c in zip(ch_name, anode, cathode): + if ch not in [a, c] and ch in raw.ch_names: + raise ValueError('There is already a channel named "%s", please ' + 'specify a different name for the bipolar ' + 'channel using the ch_name parameter.' % ch) + if ch_info is None: ch_info = [{} for an in anode] elif not isinstance(ch_info, list): @@ -1563,6 +1574,7 @@ def set_bipolar_reference(raw, anode, cathode, ch_name, ch_info=None, raw.info['chs'][an_idx] = info raw.info['chs'][an_idx]['ch_name'] = name raw.info['ch_names'][an_idx] = name + logger.info('Bipolar channel added as "%s".' % name) # Drop cathode channels raw.drop_channels(cathode) diff --git a/mne/io/fiff/tests/test_raw.py b/mne/io/fiff/tests/test_raw.py index caafb456551..b1fd2f53626 100644 --- a/mne/io/fiff/tests/test_raw.py +++ b/mne/io/fiff/tests/test_raw.py @@ -1096,20 +1096,23 @@ def test_set_bipolar_reference(): # Test creating a bipolar reference that doesn't involve EEG channels: # it should not set the custom_ref_applied flag - reref = set_bipolar_reference(raw, 'MEG 0111', 'MEG 0112', 'bipolar', - {'kind': FIFF.FIFFV_MEG_CH}) + reref = set_bipolar_reference(raw, 'MEG 0111', 'MEG 0112', + ch_info={'kind': FIFF.FIFFV_MEG_CH}) assert_true(not reref.info['custom_ref_applied']) + assert_true('MEG 0111-MEG 0112' in reref.ch_names) # Test a battery of invalid inputs - assert_raises(AssertionError, set_bipolar_reference, raw, + assert_raises(ValueError, set_bipolar_reference, raw, 'EEG 001', ['EEG 002', 'EEG 003'], 'bipolar') - assert_raises(AssertionError, set_bipolar_reference, raw, + assert_raises(ValueError, set_bipolar_reference, raw, ['EEG 001', 'EEG 002'], 'EEG 003', 'bipolar') - assert_raises(AssertionError, set_bipolar_reference, raw, + assert_raises(ValueError, set_bipolar_reference, raw, 'EEG 001', 'EEG 002', ['bipolar1', 'bipolar2']) - assert_raises(AssertionError, set_bipolar_reference, raw, + assert_raises(ValueError, set_bipolar_reference, raw, 'EEG 001', 'EEG 002', 'bipolar', - [{'foo': 'bar'}, {'foo': 'bar'}]) + ch_info=[{'foo': 'bar'}, {'foo': 'bar'}]) + assert_raises(ValueError, set_bipolar_reference, raw, + 'EEG 001', 'EEG 002', ch_name='EEG 003') @testing.requires_testing_data diff --git a/mne/tests/test_proj.py b/mne/tests/test_proj.py index 1e735634eac..1363f180953 100644 --- a/mne/tests/test_proj.py +++ b/mne/tests/test_proj.py @@ -13,8 +13,10 @@ from mne import pick_types from mne.io import Raw from mne import compute_proj_epochs, compute_proj_evoked, compute_proj_raw -from mne.io.proj import make_projector, activate_proj -from mne.proj import read_proj, write_proj, make_eeg_average_ref_proj +from mne.io.proj import (make_projector, activate_proj, + _needs_eeg_average_ref_proj) +from mne.proj import (read_proj, write_proj, make_eeg_average_ref_proj, + _has_eeg_average_ref_proj) from mne import read_events, Epochs, sensitivity_map, read_source_estimate from mne.utils import _TempDir, run_tests_if_main, clean_warning_registry @@ -218,7 +220,7 @@ def test_compute_proj_raw(): def test_make_eeg_average_ref_proj(): - ''' Test EEG average reference projection ''' + """Test EEG average reference projection""" raw = Raw(raw_fname, add_eeg_ref=False, preload=True) eeg = mne.pick_types(raw.info, meg=False, eeg=True) @@ -234,7 +236,34 @@ def test_make_eeg_average_ref_proj(): # Error when custom reference has already been applied raw.info['custom_ref_applied'] = True - assert_raises(AssertionError, make_eeg_average_ref_proj, raw.info) + assert_raises(RuntimeError, make_eeg_average_ref_proj, raw.info) +def test_has_eeg_average_ref_proj(): + """Test checking whether an EEG average reference exists""" + assert_true(not _has_eeg_average_ref_proj([])) + + raw = Raw(raw_fname, add_eeg_ref=True, preload=False) + assert_true(_has_eeg_average_ref_proj(raw.info['projs'])) + + +def test_needs_eeg_average_ref_proj(): + """Test checking whether a recording needs an EEG average reference""" + raw = Raw(raw_fname, add_eeg_ref=False, preload=False) + assert_true(_needs_eeg_average_ref_proj(raw.info)) + + raw = Raw(raw_fname, add_eeg_ref=True, preload=False) + assert_true(not _needs_eeg_average_ref_proj(raw.info)) + + # No EEG channels + raw = Raw(raw_fname, add_eeg_ref=False, preload=True) + eeg = [raw.ch_names[c] for c in pick_types(raw.info, meg=False, eeg=True)] + raw.drop_channels(eeg) + assert_true(not _needs_eeg_average_ref_proj(raw.info)) + + # Custom ref flag set + raw = Raw(raw_fname, add_eeg_ref=False, preload=False) + raw.info['custom_ref_applied'] = True + assert_true(not _needs_eeg_average_ref_proj(raw.info)) + run_tests_if_main() From c0a00a8a50ef8152694a5a2ad8f36198372e8542 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Tue, 16 Dec 2014 11:03:55 +0200 Subject: [PATCH 04/12] Updating ref and what's new --- doc/source/python_reference.rst | 12 ++++++++++++ doc/source/whats_new.rst | 4 +++- 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/doc/source/python_reference.rst b/doc/source/python_reference.rst index 6b4a9fb45dc..dd833656687 100644 --- a/doc/source/python_reference.rst +++ b/doc/source/python_reference.rst @@ -358,6 +358,18 @@ Manipulate channels and set sensors locations for processing and plotting: run_ica infomax +EEG referencing: + +.. currentmodule:: mne.io + +.. autosummary:: + :toctree: generated/ + :template: function.rst + + apply_reference + set_bipolar_reference + set_eeg_reference + :py:mod:`mne.filter`: .. automodule:: mne.filter diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst index fdb00d27e9d..6d95c8c3b62 100644 --- a/doc/source/whats_new.rst +++ b/doc/source/whats_new.rst @@ -48,7 +48,7 @@ Changelog by `Denis Engemann`_ - Add parameter to `whiten_evoked`, `compute_whitener` and `prepare_noise_cov` to set the exact rank by `Martin Luessi`_ and `Denis Engemann`_ - + - Add support for bipolar referencing by `Marijn van Vliet`_ BUG ~~~ @@ -80,6 +80,8 @@ BUG - Fix bug in `mne.viz.plot_topo` if ylim was passed for single sensor layouts by `Denis Engemann`_ + - Average reference projections will no longer by automatically added after applying a custom EEG reference by `Marijn van Vliet`_ + API ~~~ From 51ac222a8a8e70d09868015ed336609ef5d01439 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Tue, 16 Dec 2014 11:29:43 +0200 Subject: [PATCH 05/12] Some improvements to the docs --- examples/plot_rereference_eeg.py | 62 ++++++++++++++++++++++++++++++++ mne/io/base.py | 23 +++++++----- mne/io/fiff/tests/test_raw.py | 3 ++ 3 files changed, 80 insertions(+), 8 deletions(-) create mode 100644 examples/plot_rereference_eeg.py diff --git a/examples/plot_rereference_eeg.py b/examples/plot_rereference_eeg.py new file mode 100644 index 00000000000..69bd0e5791d --- /dev/null +++ b/examples/plot_rereference_eeg.py @@ -0,0 +1,62 @@ +""" +============================= +Re-referencing the EEG signal +============================= + +Load raw data and apply some EEG referencing schemes +""" +# Author: Marijn van Vliet +# +# License: BSD (3-clause) + +import mne +from mne.datasets import sample +from matplotlib import pyplot as plt + +# Setup for reading the raw data +data_path = sample.data_path() +raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif' +event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif' +event_id, tmin, tmax = 1, -0.2, 0.5 + +# Read the raw data +raw = mne.io.Raw(raw_fname, preload=True) +events = mne.read_events(event_fname) + +# Apply different EEG referencing schemes +######################################### + +# No reference. This assumes that the EEG has already been referenced properly. +# This explicitly prevents MNE from adding a default EEG reference. +raw_none = mne.io.set_eeg_reference(raw, [])[0] + +# Average reference. This is normally added by default, but can also be added +# explicitly. +raw_car = mne.io.set_eeg_reference(raw)[0] + +# Use the mean of channels EEG 001 and EEG 002 as a reference +raw_custom = mne.io.set_eeg_reference(raw, ['EEG 001', 'EEG 002'])[0] + +# Plot the evokeds to see the effect of the EEG reference +######################################################### +picks = mne.pick_types(raw.info, meg=False, eeg=True, exclude='bads') +ev_none = mne.Epochs(raw_none, events, event_id, tmin, tmax).average(picks) +ev_car = mne.Epochs(raw_car, events, event_id, tmin, tmax).average(picks) +ev_custom = mne.Epochs(raw_custom, events, event_id, tmin, tmax).average(picks) + +plt.figure(figsize=(6, 6)) + +ax = plt.subplot(3, 1, 1) +ev_none.plot(axes=ax) +plt.title('Original reference') + +ax = plt.subplot(3, 1, 2) +ev_car.plot(axes=ax) +plt.title('Average reference') + +ax = plt.subplot(3, 1, 3) +ev_custom.plot(axes=ax) +plt.title('Custom reference') + +plt.tight_layout() +plt.show() diff --git a/mne/io/base.py b/mne/io/base.py index 1596ccd3a3e..edb3fd6e666 100644 --- a/mne/io/base.py +++ b/mne/io/base.py @@ -1334,7 +1334,9 @@ def apply_reference(raw, ref_from, ref_to=None, copy=True): raw : instance of Raw Instance of Raw with EEG channels and reference channel(s). ref_from : list of str - The names of the channels to use to construct the reference. + The names of the channels to use to construct the reference. If an + empty list is specified, the data is assumed to already have a proper + reference and MNE will not attempt any re-referencing of the data. ref_to : list of str | None The names of the channels to apply the reference to. By default, all EEG channels are chosen. @@ -1381,13 +1383,15 @@ def apply_reference(raw, ref_from, ref_to=None, copy=True): else: ref_to = [raw.ch_names.index(ch) for ch in ref_to] - # Compute reference - ref_data = raw._data[ref_from].mean(0) - if copy: raw = raw.copy() - raw._data[ref_to] -= ref_data + # Compute reference + if len(ref_from) > 0: + ref_data = raw._data[ref_from].mean(0) + raw._data[ref_to] -= ref_data + else: + ref_data = None # If the reference touches EEG electrodes, remove any pre-existing common # reference and note in the info that a non-CAR has been applied. @@ -1408,7 +1412,8 @@ def apply_reference(raw, ref_from, ref_to=None, copy=True): def set_eeg_reference(raw, ref_channels=None, copy=True): """Rereference EEG channels to new reference channel(s). - If multiple reference channels are specified, they will be averaged. + If multiple reference channels are specified, they will be averaged. If + no reference channels are specified, an average reference will be applied. Parameters ---------- @@ -1417,7 +1422,9 @@ def set_eeg_reference(raw, ref_channels=None, copy=True): ref_channels : list of str | None The names of the channels to use to construct the reference. If None is specified here, an average reference will be applied in the form of an - SSP projector. + SSP projector. If an empty list is specified, the data is assumed to + already have a proper reference and MNE will not attempt any + re-referencing of the data. copy : bool Specifies whether instance of Raw will be copied (True) or modified in place (False). Defaults to True. @@ -1473,7 +1480,7 @@ def set_bipolar_reference(raw, anode, cathode, ch_name=None, ch_info=None, 1st cathode, the 2nd anode from the 2nd cathode, etc. By default, the virtual channels will be annotated with channel info of - the anodes and its location set to (0, 0, 0)) and coil type set to + the anodes, their locations set to (0, 0, 0) and coil types set to EEG_BIPOLAR. Parameters diff --git a/mne/io/fiff/tests/test_raw.py b/mne/io/fiff/tests/test_raw.py index b1fd2f53626..3abdff408f8 100644 --- a/mne/io/fiff/tests/test_raw.py +++ b/mne/io/fiff/tests/test_raw.py @@ -1020,6 +1020,9 @@ def test_apply_reference(): assert_allclose(raw_eeg_data, unref_eeg_data, 1e-6, atol=1e-15) assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) + # Test that disabling the reference does not break anything + reref, ref_data = apply_reference(raw, []) + # Test that data is modified in place when copy=False reref, ref_data = apply_reference(raw, ['EEG 001', 'EEG 002'], copy=False) From 21c2227a4bea84c8f1928077c303d6cd6871a2d7 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Tue, 16 Dec 2014 16:12:37 +0200 Subject: [PATCH 06/12] Make apply_reference a private function --- doc/source/python_reference.rst | 1 - mne/io/__init__.py | 4 ++-- mne/io/base.py | 8 +++----- mne/io/fiff/tests/test_raw.py | 14 +++++++------- 4 files changed, 12 insertions(+), 15 deletions(-) diff --git a/doc/source/python_reference.rst b/doc/source/python_reference.rst index dd833656687..3df27b8e3c0 100644 --- a/doc/source/python_reference.rst +++ b/doc/source/python_reference.rst @@ -366,7 +366,6 @@ EEG referencing: :toctree: generated/ :template: function.rst - apply_reference set_bipolar_reference set_eeg_reference diff --git a/mne/io/__init__.py b/mne/io/__init__.py index 5c4f4fb99db..e277b6f22ee 100644 --- a/mne/io/__init__.py +++ b/mne/io/__init__.py @@ -30,5 +30,5 @@ # for backward compatibility from .fiff import RawFIFF from .fiff import RawFIFF as Raw -from .base import (concatenate_raws, get_chpi_positions, apply_reference, - set_eeg_reference, set_bipolar_reference) +from .base import (concatenate_raws, get_chpi_positions, set_eeg_reference, + set_bipolar_reference) diff --git a/mne/io/base.py b/mne/io/base.py index edb3fd6e666..37b1ded1d5e 100644 --- a/mne/io/base.py +++ b/mne/io/base.py @@ -1323,7 +1323,7 @@ def add_events(self, events, stim_channel=None): self._data[pick, idx - self.first_samp] += events[:, 2] -def apply_reference(raw, ref_from, ref_to=None, copy=True): +def _apply_reference(raw, ref_from, ref_to=None, copy=True): """Apply a custom EEG referencing scheme. Calculates a reference signal by taking the mean of a set of channels and @@ -1449,7 +1449,6 @@ def set_eeg_reference(raw, ref_channels=None, copy=True): See also -------- - apply_reference : Base function for rereferencing. set_bipolar_reference : Convenience function for creating bipolar references. """ @@ -1464,7 +1463,7 @@ def set_eeg_reference(raw, ref_channels=None, copy=True): return raw, None else: logger.info('Applying a custom EEG reference.') - return apply_reference(raw, ref_channels, copy=copy) + return _apply_reference(raw, ref_channels, copy=copy) def set_bipolar_reference(raw, anode, cathode, ch_name=None, ch_info=None, @@ -1511,7 +1510,6 @@ def set_bipolar_reference(raw, anode, cathode, ch_name=None, ch_info=None, See also -------- - apply_reference : Base function for rereferencing. set_eeg_reference : Convenience function for creating an EEG reference. Notes @@ -1576,7 +1574,7 @@ def set_bipolar_reference(raw, anode, cathode, ch_name=None, ch_info=None, # Perform bipolar referencing for an, ca, name, info in zip(anode, cathode, ch_name, new_ch_info): - raw, _ = apply_reference(raw, [ca], [an], copy=False) + raw, _ = _apply_reference(raw, [ca], [an], copy=False) an_idx = raw.ch_names.index(an) raw.info['chs'][an_idx] = info raw.info['chs'][an_idx]['ch_name'] = name diff --git a/mne/io/fiff/tests/test_raw.py b/mne/io/fiff/tests/test_raw.py index 3abdff408f8..66dcd189b40 100644 --- a/mne/io/fiff/tests/test_raw.py +++ b/mne/io/fiff/tests/test_raw.py @@ -20,10 +20,10 @@ from mne import pick_types, pick_channels from mne.datasets import testing from mne.io.constants import FIFF -from mne.io import (Raw, concatenate_raws, - get_chpi_positions, apply_reference, set_eeg_reference, - set_bipolar_reference) +from mne.io import (Raw, concatenate_raws, get_chpi_positions, + set_eeg_reference, set_bipolar_reference) from mne.io.proj import _has_eeg_average_ref_proj +from mne.io.base import _apply_reference from mne import concatenate_events, find_events, equalize_channels from mne.utils import (_TempDir, requires_nitime, requires_pandas, requires_mne, run_subprocess, run_tests_if_main) @@ -995,7 +995,7 @@ def test_apply_reference(): stim=True, exclude='bads') # Rereference raw data by creating a copy of original data - reref, ref_data = apply_reference(raw, ref_from=['EEG 001', 'EEG 002'], + reref, ref_data = _apply_reference(raw, ref_from=['EEG 001', 'EEG 002'], copy=True) assert_true(reref.info['custom_ref_applied']) @@ -1021,11 +1021,11 @@ def test_apply_reference(): assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) # Test that disabling the reference does not break anything - reref, ref_data = apply_reference(raw, []) + reref, ref_data = _apply_reference(raw, []) # Test that data is modified in place when copy=False - reref, ref_data = apply_reference(raw, ['EEG 001', 'EEG 002'], - copy=False) + reref, ref_data = _apply_reference(raw, ['EEG 001', 'EEG 002'], + copy=False) assert_true(raw is reref) From cad7e35d12eaac00a0b26b948e6cd89099a0ead6 Mon Sep 17 00:00:00 2001 From: Marijn van Vliet Date: Mon, 22 Dec 2014 10:17:38 +0200 Subject: [PATCH 07/12] Reduced memory consumption of rerefencing example --- examples/plot_rereference_eeg.py | 36 +++++++++++++++++--------------- 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/examples/plot_rereference_eeg.py b/examples/plot_rereference_eeg.py index 69bd0e5791d..60bc898b864 100644 --- a/examples/plot_rereference_eeg.py +++ b/examples/plot_rereference_eeg.py @@ -23,37 +23,39 @@ raw = mne.io.Raw(raw_fname, preload=True) events = mne.read_events(event_fname) -# Apply different EEG referencing schemes -######################################### +# The EEG channels will be plotted to visulualize the difference in referencing +# schemes. +picks = mne.pick_types(raw.info, meg=False, eeg=True, exclude='bads') +plt.figure(figsize=(6, 6)) + +# Apply different EEG referencing schemes and plot the resulting evokeds +######################################################################## # No reference. This assumes that the EEG has already been referenced properly. # This explicitly prevents MNE from adding a default EEG reference. raw_none = mne.io.set_eeg_reference(raw, [])[0] - -# Average reference. This is normally added by default, but can also be added -# explicitly. -raw_car = mne.io.set_eeg_reference(raw)[0] - -# Use the mean of channels EEG 001 and EEG 002 as a reference -raw_custom = mne.io.set_eeg_reference(raw, ['EEG 001', 'EEG 002'])[0] - -# Plot the evokeds to see the effect of the EEG reference -######################################################### -picks = mne.pick_types(raw.info, meg=False, eeg=True, exclude='bads') ev_none = mne.Epochs(raw_none, events, event_id, tmin, tmax).average(picks) -ev_car = mne.Epochs(raw_car, events, event_id, tmin, tmax).average(picks) -ev_custom = mne.Epochs(raw_custom, events, event_id, tmin, tmax).average(picks) - -plt.figure(figsize=(6, 6)) +del raw_none # Free memory ax = plt.subplot(3, 1, 1) ev_none.plot(axes=ax) plt.title('Original reference') +# Average reference. This is normally added by default, but can also be added +# explicitly. +raw_car = mne.io.set_eeg_reference(raw)[0] +ev_car = mne.Epochs(raw_car, events, event_id, tmin, tmax).average(picks) +del raw_car + ax = plt.subplot(3, 1, 2) ev_car.plot(axes=ax) plt.title('Average reference') +# Use the mean of channels EEG 001 and EEG 002 as a reference +raw_custom = mne.io.set_eeg_reference(raw, ['EEG 001', 'EEG 002'])[0] +ev_custom = mne.Epochs(raw_custom, events, event_id, tmin, tmax).average(picks) +del raw_custom + ax = plt.subplot(3, 1, 3) ev_custom.plot(axes=ax) plt.title('Custom reference') From f4f0e6f51c1a4d57dc474adf8712522979a8868e Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Tue, 23 Dec 2014 14:11:24 +0100 Subject: [PATCH 08/12] cosmit in example --- examples/plot_rereference_eeg.py | 39 +++++++++++++++----------------- 1 file changed, 18 insertions(+), 21 deletions(-) diff --git a/examples/plot_rereference_eeg.py b/examples/plot_rereference_eeg.py index 60bc898b864..366076d713a 100644 --- a/examples/plot_rereference_eeg.py +++ b/examples/plot_rereference_eeg.py @@ -25,40 +25,37 @@ # The EEG channels will be plotted to visulualize the difference in referencing # schemes. -picks = mne.pick_types(raw.info, meg=False, eeg=True, exclude='bads') -plt.figure(figsize=(6, 6)) +picks = mne.pick_types(raw.info, meg=False, eeg=True, eog=True, exclude='bads') +############################################################################### # Apply different EEG referencing schemes and plot the resulting evokeds -######################################################################## + +reject = dict(eeg=180e-6, eog=150e-6) +epochs_params = dict(events=events, event_id=event_id, tmin=tmin, tmax=tmax, + picks=picks) # No reference. This assumes that the EEG has already been referenced properly. # This explicitly prevents MNE from adding a default EEG reference. -raw_none = mne.io.set_eeg_reference(raw, [])[0] -ev_none = mne.Epochs(raw_none, events, event_id, tmin, tmax).average(picks) -del raw_none # Free memory +raw_no_ref, _ = mne.io.set_eeg_reference(raw, []) +evoved_no_ref = mne.Epochs(raw_no_ref, **epochs_params).average() +del raw_no_ref # Free memory -ax = plt.subplot(3, 1, 1) -ev_none.plot(axes=ax) -plt.title('Original reference') +fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1, sharex=True) +evoved_no_ref.plot(axes=ax1, titles=dict(eeg='EEG Original reference')) # Average reference. This is normally added by default, but can also be added # explicitly. -raw_car = mne.io.set_eeg_reference(raw)[0] -ev_car = mne.Epochs(raw_car, events, event_id, tmin, tmax).average(picks) +raw_car, _ = mne.io.set_eeg_reference(raw) +evoked_car = mne.Epochs(raw_car, **epochs_params).average() del raw_car -ax = plt.subplot(3, 1, 2) -ev_car.plot(axes=ax) -plt.title('Average reference') +evoked_car.plot(axes=ax2, titles=dict(eeg='EEG Average reference')) # Use the mean of channels EEG 001 and EEG 002 as a reference -raw_custom = mne.io.set_eeg_reference(raw, ['EEG 001', 'EEG 002'])[0] -ev_custom = mne.Epochs(raw_custom, events, event_id, tmin, tmax).average(picks) +raw_custom, _ = mne.io.set_eeg_reference(raw, ['EEG 001', 'EEG 002']) +evoked_custom = mne.Epochs(raw_custom, **epochs_params).average() del raw_custom -ax = plt.subplot(3, 1, 3) -ev_custom.plot(axes=ax) -plt.title('Custom reference') +evoked_custom.plot(axes=ax3, titles=dict(eeg='EEG Custom reference')) -plt.tight_layout() -plt.show() +mne.viz.tight_layout() From 0b008a9fcf2c1440e1a86c420f17eb5979a5d33b Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Tue, 23 Dec 2014 14:26:16 +0100 Subject: [PATCH 09/12] move EEG referencing code to separate file --- doc/source/whats_new.rst | 1 + mne/io/__init__.py | 4 +- mne/io/base.py | 268 +------------------------------- mne/io/fiff/tests/test_raw.py | 139 +---------------- mne/io/reference.py | 275 +++++++++++++++++++++++++++++++++ mne/io/tests/test_reference.py | 157 +++++++++++++++++++ 6 files changed, 438 insertions(+), 406 deletions(-) create mode 100644 mne/io/reference.py create mode 100644 mne/io/tests/test_reference.py diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst index 6d95c8c3b62..862203ee024 100644 --- a/doc/source/whats_new.rst +++ b/doc/source/whats_new.rst @@ -48,6 +48,7 @@ Changelog by `Denis Engemann`_ - Add parameter to `whiten_evoked`, `compute_whitener` and `prepare_noise_cov` to set the exact rank by `Martin Luessi`_ and `Denis Engemann`_ + - Add support for bipolar referencing by `Marijn van Vliet`_ BUG diff --git a/mne/io/__init__.py b/mne/io/__init__.py index e277b6f22ee..0dedf4d78f5 100644 --- a/mne/io/__init__.py +++ b/mne/io/__init__.py @@ -30,5 +30,5 @@ # for backward compatibility from .fiff import RawFIFF from .fiff import RawFIFF as Raw -from .base import (concatenate_raws, get_chpi_positions, set_eeg_reference, - set_bipolar_reference) +from .base import concatenate_raws, get_chpi_positions +from .reference import set_eeg_reference, set_bipolar_reference diff --git a/mne/io/base.py b/mne/io/base.py index 37b1ded1d5e..f5fde7cea14 100644 --- a/mne/io/base.py +++ b/mne/io/base.py @@ -3,6 +3,7 @@ # Martin Luessi # Denis Engemann # Teon Brooks +# Marijn van Vliet # # License: BSD (3-clause) @@ -20,8 +21,7 @@ from .constants import FIFF from .pick import pick_types, channel_type, pick_channels from .meas_info import write_meas_info -from .proj import (setup_proj, activate_proj, proj_equal, ProjMixin, - _has_eeg_average_ref_proj, make_eeg_average_ref_proj) +from .proj import setup_proj, activate_proj, proj_equal, ProjMixin from ..channels.channels import (ContainsMixin, PickDropChannelsMixin, SetChannelsMixin, InterpolationMixin) from ..channels.layout import read_montage, apply_montage, Montage @@ -1323,270 +1323,6 @@ def add_events(self, events, stim_channel=None): self._data[pick, idx - self.first_samp] += events[:, 2] -def _apply_reference(raw, ref_from, ref_to=None, copy=True): - """Apply a custom EEG referencing scheme. - - Calculates a reference signal by taking the mean of a set of channels and - applies the reference to another set of channels. - - Parameters - ---------- - raw : instance of Raw - Instance of Raw with EEG channels and reference channel(s). - ref_from : list of str - The names of the channels to use to construct the reference. If an - empty list is specified, the data is assumed to already have a proper - reference and MNE will not attempt any re-referencing of the data. - ref_to : list of str | None - The names of the channels to apply the reference to. By default, - all EEG channels are chosen. - copy : bool - Specifies whether instance of Raw will be copied or modified in place. - - Returns - ------- - raw : instance of Raw - Instance of Raw with EEG channels rereferenced. - ref_data : array - Array of reference data subtracted from EEG channels. - - Notes - ----- - 1. Do not use this function to apply an average reference. By default, an - average reference projection has already been added upon loading raw - data. - - 2. If the reference is applied to any EEG channels, this function removes - any pre-existing average reference projections. - - 3. During source localization, the EEG signal should have an average - reference. - - 4. The data must be preloaded. - - See also - -------- - set_eeg_reference : Convenience function for creating an EEG reference. - set_bipolar_reference : Convenience function for creating a bipolar - reference. - """ - # Check to see that raw data is preloaded - if not raw.preload: - raise RuntimeError('Raw data needs to be preloaded. Use ' - 'preload=True (or string) in the constructor.') - - eeg_idx = pick_types(raw.info, eeg=True, meg=False, ref_meg=False) - - ref_from = [raw.ch_names.index(ch) for ch in ref_from] - if ref_to is None: - ref_to = eeg_idx - else: - ref_to = [raw.ch_names.index(ch) for ch in ref_to] - - if copy: - raw = raw.copy() - - # Compute reference - if len(ref_from) > 0: - ref_data = raw._data[ref_from].mean(0) - raw._data[ref_to] -= ref_data - else: - ref_data = None - - # If the reference touches EEG electrodes, remove any pre-existing common - # reference and note in the info that a non-CAR has been applied. - if len(np.intersect1d(ref_to, eeg_idx)) > 0: - raw.info['custom_ref_applied'] = True - - # Remove any existing average reference projections - for i, proj in enumerate(raw.info['projs']): - if proj['desc'] == 'Average EEG reference' or \ - proj['kind'] == FIFF.FIFFV_MNE_PROJ_ITEM_EEG_AVREF: - logger.info('Removing existing average EEG reference ' - 'projection.') - del raw.info['projs'][i] - - return raw, ref_data - - -def set_eeg_reference(raw, ref_channels=None, copy=True): - """Rereference EEG channels to new reference channel(s). - - If multiple reference channels are specified, they will be averaged. If - no reference channels are specified, an average reference will be applied. - - Parameters - ---------- - raw : instance of Raw - Instance of Raw with EEG channels and reference channel(s). - ref_channels : list of str | None - The names of the channels to use to construct the reference. If None is - specified here, an average reference will be applied in the form of an - SSP projector. If an empty list is specified, the data is assumed to - already have a proper reference and MNE will not attempt any - re-referencing of the data. - copy : bool - Specifies whether instance of Raw will be copied (True) or modified in - place (False). Defaults to True. - - Returns - ------- - raw : instance of Raw - Instance of Raw with eeg channels rereferenced. - ref_data : array - Array of reference data subtracted from EEG channels. - - Notes - ----- - 1. If a reference is requested that is not the average reference, this - function removes any pre-existing average reference projections. - - 2. During source localization, the EEG signal should have an average - reference. - - 3. In order to apply a reference other than an average reference, the data - must be preloaded. - - See also - -------- - set_bipolar_reference : Convenience function for creating bipolar - references. - """ - if ref_channels is None: - # CAR requested - if _has_eeg_average_ref_proj(raw.info['projs']): - logger.warning('An average reference projection was already ' - 'added. The data has been left untouched.') - return raw, None - else: - raw.add_proj(make_eeg_average_ref_proj(raw.info), activate=False) - return raw, None - else: - logger.info('Applying a custom EEG reference.') - return _apply_reference(raw, ref_channels, copy=copy) - - -def set_bipolar_reference(raw, anode, cathode, ch_name=None, ch_info=None, - copy=True): - """Rereference selected channels using a bipolar referencing scheme. - - A bipolar reference takes the difference between two channels (the anode - minus the cathode) and adds it as a new virtual channel. The original - channels will be dropped. - - Multiple anodes and cathodes can be specified, in which case multiple - vitual channels will be created. The 1st anode will be substracted from the - 1st cathode, the 2nd anode from the 2nd cathode, etc. - - By default, the virtual channels will be annotated with channel info of - the anodes, their locations set to (0, 0, 0) and coil types set to - EEG_BIPOLAR. - - Parameters - ---------- - raw : instance of Raw - Instance of Raw containing the unreferenced channels. - anode : str | list of str - The name(s) of the channel(s) to use as anode in the bipolar reference. - cathode : str | list of str - The name(s) of the channel(s) to use as cathode in the bipolar - reference. - ch_name : str | list of str | None - The channel name(s) for the virtual channel(s) containing the resulting - signal. By default, bipolar channels are named after the anode and - cathode, but it is recommended to supply a more meaningful name. - ch_info : dict | list of dict | None - This parameter can be used to supply a dictionary (or a dictionary for - each bipolar channel) containing channel information to merge in, - overwriting the default values. - copy : bool - Whether to operate on a copy of the raw data (True) or modify it - in-place (False). Defaults to True. - - Returns - ------- - raw : instance of Raw - Instance of Raw with the specified channels rereferenced. - - See also - -------- - set_eeg_reference : Convenience function for creating an EEG reference. - - Notes - ----- - 1. If the anodes contain any EEG channels, this function removes - any pre-existing average reference projections. - - 2. During source localization, the EEG signal should have an average - reference. - - 3. The data must be preloaded. - """ - if not isinstance(anode, list): - anode = [anode] - - if not isinstance(cathode, list): - cathode = [cathode] - - if len(anode) != len(cathode): - raise ValueError('Number of anodes must equal the number of cathodes.') - - if ch_name is None: - ch_name = ['%s-%s' % ac for ac in zip(anode, cathode)] - elif not isinstance(ch_name, list): - ch_name = [ch_name] - if len(ch_name) != len(anode): - raise ValueError('Number of channel names must equal the number of ' - 'anodes/cathodes.') - - # Check for duplicate channel names (it is allowed to give the name of the - # anode or cathode channel, as they will be replaced). - for ch, a, c in zip(ch_name, anode, cathode): - if ch not in [a, c] and ch in raw.ch_names: - raise ValueError('There is already a channel named "%s", please ' - 'specify a different name for the bipolar ' - 'channel using the ch_name parameter.' % ch) - - if ch_info is None: - ch_info = [{} for an in anode] - elif not isinstance(ch_info, list): - ch_info = [ch_info] - if len(ch_info) != len(anode): - raise ValueError('Number of channel info dictionaries must equal the ' - 'number of anodes/cathodes.') - - # Merge specified and anode channel information dictionaries - new_ch_info = [] - for an, ci in zip(anode, ch_info): - new_info = raw.info['chs'][raw.ch_names.index(an)].copy() - - # Set channel location and coil type - if 'eeg_loc' in new_info: - new_info['eeg_loc'] = np.zeros((3, 2)) - new_info['loc'] = np.zeros(12) - new_info['coil_type'] = FIFF.FIFFV_COIL_EEG_BIPOLAR - - new_info.update(ci) - new_ch_info.append(new_info) - - if copy: - raw = raw.copy() - - # Perform bipolar referencing - for an, ca, name, info in zip(anode, cathode, ch_name, new_ch_info): - raw, _ = _apply_reference(raw, [ca], [an], copy=False) - an_idx = raw.ch_names.index(an) - raw.info['chs'][an_idx] = info - raw.info['chs'][an_idx]['ch_name'] = name - raw.info['ch_names'][an_idx] = name - logger.info('Bipolar channel added as "%s".' % name) - - # Drop cathode channels - raw.drop_channels(cathode) - - return raw - - def _allocate_data(data, data_buffer, data_shape, dtype): if data is None: # if not already done, allocate array with right type diff --git a/mne/io/fiff/tests/test_raw.py b/mne/io/fiff/tests/test_raw.py index 66dcd189b40..408a7a7ddc3 100644 --- a/mne/io/fiff/tests/test_raw.py +++ b/mne/io/fiff/tests/test_raw.py @@ -20,10 +20,7 @@ from mne import pick_types, pick_channels from mne.datasets import testing from mne.io.constants import FIFF -from mne.io import (Raw, concatenate_raws, get_chpi_positions, - set_eeg_reference, set_bipolar_reference) -from mne.io.proj import _has_eeg_average_ref_proj -from mne.io.base import _apply_reference +from mne.io import Raw, concatenate_raws, get_chpi_positions from mne import concatenate_events, find_events, equalize_channels from mne.utils import (_TempDir, requires_nitime, requires_pandas, requires_mne, run_subprocess, run_tests_if_main) @@ -984,140 +981,6 @@ def compensate_mne(fname, grad): assert_allclose(raw_py._data, raw_c._data, rtol=1e-6, atol=1e-17) -@testing.requires_testing_data -def test_apply_reference(): - """ Test base function for rereferencing""" - raw = Raw(fif_fname, preload=True, add_eeg_ref=True) - - # Separate EEG channels from other channel types - picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') - picks_other = pick_types(raw.info, meg=True, eeg=False, eog=True, - stim=True, exclude='bads') - - # Rereference raw data by creating a copy of original data - reref, ref_data = _apply_reference(raw, ref_from=['EEG 001', 'EEG 002'], - copy=True) - assert_true(reref.info['custom_ref_applied']) - - # The CAR reference projection should have been removed by the function - assert_true(not _has_eeg_average_ref_proj(reref.info['projs'])) - - # Check that the ref has been properly computed - ref_ch_idx = [raw.ch_names.index(ch) for ch in ['EEG 001', 'EEG 002']] - assert_array_equal(ref_data, raw[ref_ch_idx, :][0].mean(0)) - - # Get the raw EEG data and other channel data - raw_eeg_data = raw[picks_eeg][0] - raw_other_data = raw[picks_other][0] - - # Get the rereferenced EEG data and channel other - reref_eeg_data = reref[picks_eeg][0] - # Undo rereferencing of EEG channels - unref_eeg_data = reref_eeg_data + ref_data - reref_other_data = reref[picks_other][0] - - # Check that both EEG data and other data is the same - assert_allclose(raw_eeg_data, unref_eeg_data, 1e-6, atol=1e-15) - assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) - - # Test that disabling the reference does not break anything - reref, ref_data = _apply_reference(raw, []) - - # Test that data is modified in place when copy=False - reref, ref_data = _apply_reference(raw, ['EEG 001', 'EEG 002'], - copy=False) - assert_true(raw is reref) - - -@testing.requires_testing_data -def test_set_eeg_reference(): - """ Test rereference eeg data""" - raw = Raw(fif_fname, preload=True) - - # Rereference raw data by creating a copy of original data - reref, ref_data = set_eeg_reference(raw, ['EEG 001', 'EEG 002'], copy=True) - assert_true(reref.info['custom_ref_applied']) - - # Separate EEG channels from other channel types - picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') - picks_other = pick_types(raw.info, meg=True, eeg=False, eog=True, - stim=True, exclude='bads') - - # Get the raw EEG data and other channel data - raw_eeg_data = raw[picks_eeg][0] - raw_other_data = raw[picks_other][0] - - # Get the rereferenced EEG data and channel other - reref_eeg_data = reref[picks_eeg][0] - # Undo rereferencing of EEG channels - unref_eeg_data = reref_eeg_data + ref_data - reref_other_data = reref[picks_other][0] - - # Check that both EEG data and other data is the same - assert_allclose(raw_eeg_data, unref_eeg_data, 1e-6, atol=1e-15) - assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) - - # Test that data is modified in place when copy=False - reref, ref_data = set_eeg_reference(raw, ['EEG 001', 'EEG 002'], - copy=False) - assert_true(raw is reref) - - -@testing.requires_testing_data -def test_set_bipolar_reference(): - """ Test bipolar referencing""" - raw = Raw(fif_fname, preload=True) - reref = set_bipolar_reference(raw, 'EEG 001', 'EEG 002', 'bipolar', - {'kind': FIFF.FIFFV_EOG_CH, - 'extra': 'some extra value'}) - assert_true(reref.info['custom_ref_applied']) - - # Compare result to a manual calculation - a = raw.pick_channels(['EEG 001', 'EEG 002'], copy=True) - a = a._data[0, :] - a._data[1, :] - b = reref.pick_channels(['bipolar'], copy=True)._data[0, :] - assert_allclose(a, b) - - # Original channels should be replaced by a virtual one - assert_true('EEG 001' not in reref.ch_names) - assert_true('EEG 002' not in reref.ch_names) - assert_true('bipolar' in reref.ch_names) - - # Check channel information - bp_info = reref.info['chs'][reref.ch_names.index('bipolar')] - an_info = reref.info['chs'][raw.ch_names.index('EEG 001')] - for key in bp_info: - if key == 'loc' or key == 'eeg_loc': - assert_array_equal(bp_info[key], 0) - elif key == 'coil_type': - assert_equal(bp_info[key], FIFF.FIFFV_COIL_EEG_BIPOLAR) - elif key == 'kind': - assert_equal(bp_info[key], FIFF.FIFFV_EOG_CH) - else: - assert_equal(bp_info[key], an_info[key]) - assert_equal(bp_info['extra'], 'some extra value') - - # Test creating a bipolar reference that doesn't involve EEG channels: - # it should not set the custom_ref_applied flag - reref = set_bipolar_reference(raw, 'MEG 0111', 'MEG 0112', - ch_info={'kind': FIFF.FIFFV_MEG_CH}) - assert_true(not reref.info['custom_ref_applied']) - assert_true('MEG 0111-MEG 0112' in reref.ch_names) - - # Test a battery of invalid inputs - assert_raises(ValueError, set_bipolar_reference, raw, - 'EEG 001', ['EEG 002', 'EEG 003'], 'bipolar') - assert_raises(ValueError, set_bipolar_reference, raw, - ['EEG 001', 'EEG 002'], 'EEG 003', 'bipolar') - assert_raises(ValueError, set_bipolar_reference, raw, - 'EEG 001', 'EEG 002', ['bipolar1', 'bipolar2']) - assert_raises(ValueError, set_bipolar_reference, raw, - 'EEG 001', 'EEG 002', 'bipolar', - ch_info=[{'foo': 'bar'}, {'foo': 'bar'}]) - assert_raises(ValueError, set_bipolar_reference, raw, - 'EEG 001', 'EEG 002', ch_name='EEG 003') - - @testing.requires_testing_data def test_drop_channels_mixin(): """Test channels-dropping functionality diff --git a/mne/io/reference.py b/mne/io/reference.py new file mode 100644 index 00000000000..78f0085e8d8 --- /dev/null +++ b/mne/io/reference.py @@ -0,0 +1,275 @@ +# Authors: Marijn van Vliet +# Alexandre Gramfort +# +# License: BSD (3-clause) + +import numpy as np + +from .constants import FIFF +from .proj import _has_eeg_average_ref_proj, make_eeg_average_ref_proj +from .pick import pick_types +from ..utils import logger + + +def _apply_reference(raw, ref_from, ref_to=None, copy=True): + """Apply a custom EEG referencing scheme. + + Calculates a reference signal by taking the mean of a set of channels and + applies the reference to another set of channels. + + Parameters + ---------- + raw : instance of Raw + Instance of Raw with EEG channels and reference channel(s). + ref_from : list of str + The names of the channels to use to construct the reference. If an + empty list is specified, the data is assumed to already have a proper + reference and MNE will not attempt any re-referencing of the data. + ref_to : list of str | None + The names of the channels to apply the reference to. By default, + all EEG channels are chosen. + copy : bool + Specifies whether instance of Raw will be copied or modified in place. + + Returns + ------- + raw : instance of Raw + Instance of Raw with EEG channels rereferenced. + ref_data : array + Array of reference data subtracted from EEG channels. + + Notes + ----- + 1. Do not use this function to apply an average reference. By default, an + average reference projection has already been added upon loading raw + data. + + 2. If the reference is applied to any EEG channels, this function removes + any pre-existing average reference projections. + + 3. During source localization, the EEG signal should have an average + reference. + + 4. The data must be preloaded. + + See also + -------- + set_eeg_reference : Convenience function for creating an EEG reference. + set_bipolar_reference : Convenience function for creating a bipolar + reference. + """ + # Check to see that raw data is preloaded + if not raw.preload: + raise RuntimeError('Raw data needs to be preloaded. Use ' + 'preload=True (or string) in the constructor.') + + eeg_idx = pick_types(raw.info, eeg=True, meg=False, ref_meg=False) + + ref_from = [raw.ch_names.index(ch) for ch in ref_from] + if ref_to is None: + ref_to = eeg_idx + else: + ref_to = [raw.ch_names.index(ch) for ch in ref_to] + + if copy: + raw = raw.copy() + + # Compute reference + if len(ref_from) > 0: + ref_data = raw._data[ref_from].mean(0) + raw._data[ref_to] -= ref_data + else: + ref_data = None + + # If the reference touches EEG electrodes, remove any pre-existing common + # reference and note in the info that a non-CAR has been applied. + if len(np.intersect1d(ref_to, eeg_idx)) > 0: + raw.info['custom_ref_applied'] = True + + # Remove any existing average reference projections + for i, proj in enumerate(raw.info['projs']): + if proj['desc'] == 'Average EEG reference' or \ + proj['kind'] == FIFF.FIFFV_MNE_PROJ_ITEM_EEG_AVREF: + logger.info('Removing existing average EEG reference ' + 'projection.') + del raw.info['projs'][i] + + return raw, ref_data + + +def set_eeg_reference(raw, ref_channels=None, copy=True): + """Rereference EEG channels to new reference channel(s). + + If multiple reference channels are specified, they will be averaged. If + no reference channels are specified, an average reference will be applied. + + Parameters + ---------- + raw : instance of Raw + Instance of Raw with EEG channels and reference channel(s). + ref_channels : list of str | None + The names of the channels to use to construct the reference. If None is + specified here, an average reference will be applied in the form of an + SSP projector. If an empty list is specified, the data is assumed to + already have a proper reference and MNE will not attempt any + re-referencing of the data. + copy : bool + Specifies whether instance of Raw will be copied (True) or modified in + place (False). Defaults to True. + + Returns + ------- + raw : instance of Raw + Instance of Raw with eeg channels rereferenced. + ref_data : array + Array of reference data subtracted from EEG channels. + + Notes + ----- + 1. If a reference is requested that is not the average reference, this + function removes any pre-existing average reference projections. + + 2. During source localization, the EEG signal should have an average + reference. + + 3. In order to apply a reference other than an average reference, the data + must be preloaded. + + See also + -------- + set_bipolar_reference : Convenience function for creating bipolar + references. + """ + if ref_channels is None: + # CAR requested + if _has_eeg_average_ref_proj(raw.info['projs']): + logger.warning('An average reference projection was already ' + 'added. The data has been left untouched.') + return raw, None + else: + raw.add_proj(make_eeg_average_ref_proj(raw.info), activate=False) + return raw, None + else: + logger.info('Applying a custom EEG reference.') + return _apply_reference(raw, ref_channels, copy=copy) + + +def set_bipolar_reference(raw, anode, cathode, ch_name=None, ch_info=None, + copy=True): + """Rereference selected channels using a bipolar referencing scheme. + + A bipolar reference takes the difference between two channels (the anode + minus the cathode) and adds it as a new virtual channel. The original + channels will be dropped. + + Multiple anodes and cathodes can be specified, in which case multiple + vitual channels will be created. The 1st anode will be substracted from the + 1st cathode, the 2nd anode from the 2nd cathode, etc. + + By default, the virtual channels will be annotated with channel info of + the anodes, their locations set to (0, 0, 0) and coil types set to + EEG_BIPOLAR. + + Parameters + ---------- + raw : instance of Raw + Instance of Raw containing the unreferenced channels. + anode : str | list of str + The name(s) of the channel(s) to use as anode in the bipolar reference. + cathode : str | list of str + The name(s) of the channel(s) to use as cathode in the bipolar + reference. + ch_name : str | list of str | None + The channel name(s) for the virtual channel(s) containing the resulting + signal. By default, bipolar channels are named after the anode and + cathode, but it is recommended to supply a more meaningful name. + ch_info : dict | list of dict | None + This parameter can be used to supply a dictionary (or a dictionary for + each bipolar channel) containing channel information to merge in, + overwriting the default values. + copy : bool + Whether to operate on a copy of the raw data (True) or modify it + in-place (False). Defaults to True. + + Returns + ------- + raw : instance of Raw + Instance of Raw with the specified channels rereferenced. + + See also + -------- + set_eeg_reference : Convenience function for creating an EEG reference. + + Notes + ----- + 1. If the anodes contain any EEG channels, this function removes + any pre-existing average reference projections. + + 2. During source localization, the EEG signal should have an average + reference. + + 3. The data must be preloaded. + """ + if not isinstance(anode, list): + anode = [anode] + + if not isinstance(cathode, list): + cathode = [cathode] + + if len(anode) != len(cathode): + raise ValueError('Number of anodes must equal the number of cathodes.') + + if ch_name is None: + ch_name = ['%s-%s' % ac for ac in zip(anode, cathode)] + elif not isinstance(ch_name, list): + ch_name = [ch_name] + if len(ch_name) != len(anode): + raise ValueError('Number of channel names must equal the number of ' + 'anodes/cathodes.') + + # Check for duplicate channel names (it is allowed to give the name of the + # anode or cathode channel, as they will be replaced). + for ch, a, c in zip(ch_name, anode, cathode): + if ch not in [a, c] and ch in raw.ch_names: + raise ValueError('There is already a channel named "%s", please ' + 'specify a different name for the bipolar ' + 'channel using the ch_name parameter.' % ch) + + if ch_info is None: + ch_info = [{} for an in anode] + elif not isinstance(ch_info, list): + ch_info = [ch_info] + if len(ch_info) != len(anode): + raise ValueError('Number of channel info dictionaries must equal the ' + 'number of anodes/cathodes.') + + # Merge specified and anode channel information dictionaries + new_ch_info = [] + for an, ci in zip(anode, ch_info): + new_info = raw.info['chs'][raw.ch_names.index(an)].copy() + + # Set channel location and coil type + if 'eeg_loc' in new_info: + new_info['eeg_loc'] = np.zeros((3, 2)) + new_info['loc'] = np.zeros(12) + new_info['coil_type'] = FIFF.FIFFV_COIL_EEG_BIPOLAR + + new_info.update(ci) + new_ch_info.append(new_info) + + if copy: + raw = raw.copy() + + # Perform bipolar referencing + for an, ca, name, info in zip(anode, cathode, ch_name, new_ch_info): + raw, _ = _apply_reference(raw, [ca], [an], copy=False) + an_idx = raw.ch_names.index(an) + raw.info['chs'][an_idx] = info + raw.info['chs'][an_idx]['ch_name'] = name + raw.info['ch_names'][an_idx] = name + logger.info('Bipolar channel added as "%s".' % name) + + # Drop cathode channels + raw.drop_channels(cathode) + + return raw diff --git a/mne/io/tests/test_reference.py b/mne/io/tests/test_reference.py new file mode 100644 index 00000000000..ad4b0db8b09 --- /dev/null +++ b/mne/io/tests/test_reference.py @@ -0,0 +1,157 @@ +# Authors: Marijn van Vliet +# Alexandre Gramfort +# +# License: BSD (3-clause) + +import os.path as op +import warnings + +from nose.tools import assert_true, assert_equal, assert_raises +from numpy.testing import assert_array_equal, assert_allclose + +from mne.io.constants import FIFF +from mne.io import set_eeg_reference, set_bipolar_reference +from mne.io.proj import _has_eeg_average_ref_proj +from mne.io.reference import _apply_reference +from mne.datasets import testing +from mne.io import Raw +from mne import pick_types + +warnings.simplefilter('always') # enable b/c these tests throw warnings + +data_dir = op.join(testing.data_path(download=False), 'MEG', 'sample') +fif_fname = op.join(data_dir, 'sample_audvis_trunc_raw.fif') + + +@testing.requires_testing_data +def test_apply_reference(): + """ Test base function for rereferencing""" + raw = Raw(fif_fname, preload=True, add_eeg_ref=True) + + # Separate EEG channels from other channel types + picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') + picks_other = pick_types(raw.info, meg=True, eeg=False, eog=True, + stim=True, exclude='bads') + + # Rereference raw data by creating a copy of original data + reref, ref_data = _apply_reference(raw, ref_from=['EEG 001', 'EEG 002'], + copy=True) + assert_true(reref.info['custom_ref_applied']) + + # The CAR reference projection should have been removed by the function + assert_true(not _has_eeg_average_ref_proj(reref.info['projs'])) + + # Check that the ref has been properly computed + ref_ch_idx = [raw.ch_names.index(ch) for ch in ['EEG 001', 'EEG 002']] + assert_array_equal(ref_data, raw[ref_ch_idx, :][0].mean(0)) + + # Get the raw EEG data and other channel data + raw_eeg_data = raw[picks_eeg][0] + raw_other_data = raw[picks_other][0] + + # Get the rereferenced EEG data and channel other + reref_eeg_data = reref[picks_eeg][0] + # Undo rereferencing of EEG channels + unref_eeg_data = reref_eeg_data + ref_data + reref_other_data = reref[picks_other][0] + + # Check that both EEG data and other data is the same + assert_allclose(raw_eeg_data, unref_eeg_data, 1e-6, atol=1e-15) + assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) + + # Test that disabling the reference does not break anything + reref, ref_data = _apply_reference(raw, []) + + # Test that data is modified in place when copy=False + reref, ref_data = _apply_reference(raw, ['EEG 001', 'EEG 002'], + copy=False) + assert_true(raw is reref) + + +@testing.requires_testing_data +def test_set_eeg_reference(): + """ Test rereference eeg data""" + raw = Raw(fif_fname, preload=True) + + # Rereference raw data by creating a copy of original data + reref, ref_data = set_eeg_reference(raw, ['EEG 001', 'EEG 002'], copy=True) + assert_true(reref.info['custom_ref_applied']) + + # Separate EEG channels from other channel types + picks_eeg = pick_types(raw.info, meg=False, eeg=True, exclude='bads') + picks_other = pick_types(raw.info, meg=True, eeg=False, eog=True, + stim=True, exclude='bads') + + # Get the raw EEG data and other channel data + raw_eeg_data = raw[picks_eeg][0] + raw_other_data = raw[picks_other][0] + + # Get the rereferenced EEG data and channel other + reref_eeg_data = reref[picks_eeg][0] + # Undo rereferencing of EEG channels + unref_eeg_data = reref_eeg_data + ref_data + reref_other_data = reref[picks_other][0] + + # Check that both EEG data and other data is the same + assert_allclose(raw_eeg_data, unref_eeg_data, 1e-6, atol=1e-15) + assert_allclose(raw_other_data, reref_other_data, 1e-6, atol=1e-15) + + # Test that data is modified in place when copy=False + reref, ref_data = set_eeg_reference(raw, ['EEG 001', 'EEG 002'], + copy=False) + assert_true(raw is reref) + + +@testing.requires_testing_data +def test_set_bipolar_reference(): + """ Test bipolar referencing""" + raw = Raw(fif_fname, preload=True) + reref = set_bipolar_reference(raw, 'EEG 001', 'EEG 002', 'bipolar', + {'kind': FIFF.FIFFV_EOG_CH, + 'extra': 'some extra value'}) + assert_true(reref.info['custom_ref_applied']) + + # Compare result to a manual calculation + a = raw.pick_channels(['EEG 001', 'EEG 002'], copy=True) + a = a._data[0, :] - a._data[1, :] + b = reref.pick_channels(['bipolar'], copy=True)._data[0, :] + assert_allclose(a, b) + + # Original channels should be replaced by a virtual one + assert_true('EEG 001' not in reref.ch_names) + assert_true('EEG 002' not in reref.ch_names) + assert_true('bipolar' in reref.ch_names) + + # Check channel information + bp_info = reref.info['chs'][reref.ch_names.index('bipolar')] + an_info = reref.info['chs'][raw.ch_names.index('EEG 001')] + for key in bp_info: + if key == 'loc' or key == 'eeg_loc': + assert_array_equal(bp_info[key], 0) + elif key == 'coil_type': + assert_equal(bp_info[key], FIFF.FIFFV_COIL_EEG_BIPOLAR) + elif key == 'kind': + assert_equal(bp_info[key], FIFF.FIFFV_EOG_CH) + else: + assert_equal(bp_info[key], an_info[key]) + assert_equal(bp_info['extra'], 'some extra value') + + # Test creating a bipolar reference that doesn't involve EEG channels: + # it should not set the custom_ref_applied flag + reref = set_bipolar_reference(raw, 'MEG 0111', 'MEG 0112', + ch_info={'kind': FIFF.FIFFV_MEG_CH}) + assert_true(not reref.info['custom_ref_applied']) + assert_true('MEG 0111-MEG 0112' in reref.ch_names) + + # Test a battery of invalid inputs + assert_raises(ValueError, set_bipolar_reference, raw, + 'EEG 001', ['EEG 002', 'EEG 003'], 'bipolar') + assert_raises(ValueError, set_bipolar_reference, raw, + ['EEG 001', 'EEG 002'], 'EEG 003', 'bipolar') + assert_raises(ValueError, set_bipolar_reference, raw, + 'EEG 001', 'EEG 002', ['bipolar1', 'bipolar2']) + assert_raises(ValueError, set_bipolar_reference, raw, + 'EEG 001', 'EEG 002', 'bipolar', + ch_info=[{'foo': 'bar'}, {'foo': 'bar'}]) + assert_raises(ValueError, set_bipolar_reference, raw, + 'EEG 001', 'EEG 002', ch_name='EEG 003') From 8fa8dc9a3748c088108ddd9d220383fb8d66bfa2 Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Tue, 23 Dec 2014 14:33:04 +0100 Subject: [PATCH 10/12] typos --- examples/plot_rereference_eeg.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/plot_rereference_eeg.py b/examples/plot_rereference_eeg.py index 366076d713a..b9e18d3dd6d 100644 --- a/examples/plot_rereference_eeg.py +++ b/examples/plot_rereference_eeg.py @@ -3,7 +3,7 @@ Re-referencing the EEG signal ============================= -Load raw data and apply some EEG referencing schemes +Load raw data and apply some EEG referencing schemes. """ # Author: Marijn van Vliet # @@ -23,7 +23,7 @@ raw = mne.io.Raw(raw_fname, preload=True) events = mne.read_events(event_fname) -# The EEG channels will be plotted to visulualize the difference in referencing +# The EEG channels will be plotted to visualize the difference in referencing # schemes. picks = mne.pick_types(raw.info, meg=False, eeg=True, eog=True, exclude='bads') From 8e150c203de1a9ed73f425d0e7336068131480cb Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Tue, 23 Dec 2014 14:49:01 +0100 Subject: [PATCH 11/12] prevent inverse modeling with bad or no EEG reference --- mne/beamformer/_dics.py | 5 ++++- mne/beamformer/_lcmv.py | 6 +++++- mne/inverse_sparse/_gamma_map.py | 4 +++- mne/inverse_sparse/mxne_inverse.py | 5 +++++ mne/minimum_norm/inverse.py | 13 +++++++++++++ mne/minimum_norm/tests/test_inverse.py | 7 +++++++ 6 files changed, 37 insertions(+), 3 deletions(-) diff --git a/mne/beamformer/_dics.py b/mne/beamformer/_dics.py index 84aabff85d4..9fd7a726e82 100644 --- a/mne/beamformer/_dics.py +++ b/mne/beamformer/_dics.py @@ -14,7 +14,7 @@ from ..utils import logger, verbose from ..io.pick import pick_types from ..forward import _subject_from_forward -from ..minimum_norm.inverse import combine_xyz +from ..minimum_norm.inverse import combine_xyz, _check_reference from ..source_estimate import SourceEstimate from ..time_frequency import CrossSpectralDensity, compute_epochs_csd from ._lcmv import _prepare_beamformer_input @@ -186,6 +186,7 @@ def dics(evoked, forward, noise_csd, data_csd, reg=0.01, label=None, Gross et al. Dynamic imaging of coherent sources: Studying neural interactions in the human brain. PNAS (2001) vol. 98 (2) pp. 694-699 """ + _check_reference(evoked) info = evoked.info data = evoked.data tmin = evoked.times[0] @@ -244,6 +245,7 @@ def dics_epochs(epochs, forward, noise_csd, data_csd, reg=0.01, label=None, Gross et al. Dynamic imaging of coherent sources: Studying neural interactions in the human brain. PNAS (2001) vol. 98 (2) pp. 694-699 """ + _check_reference(epochs) info = epochs.info tmin = epochs.times[0] @@ -481,6 +483,7 @@ def tf_dics(epochs, forward, noise_csds, tmin, tmax, tstep, win_lengths, NOTE : Dalal et al. used a synthetic aperture magnetometry beamformer (SAM) in each time-frequency window instead of DICS. """ + _check_reference(epochs) if pick_ori not in [None, 'normal']: raise ValueError('Unrecognized orientation option in pick_ori, ' diff --git a/mne/beamformer/_lcmv.py b/mne/beamformer/_lcmv.py index 878645af6d2..c9bf0e9a716 100644 --- a/mne/beamformer/_lcmv.py +++ b/mne/beamformer/_lcmv.py @@ -15,7 +15,7 @@ from ..io.proj import make_projector from ..io.pick import pick_types, pick_channels_forward, pick_channels_cov from ..forward import _subject_from_forward -from ..minimum_norm.inverse import _get_vertno, combine_xyz +from ..minimum_norm.inverse import _get_vertno, combine_xyz, _check_reference from ..cov import compute_whitener, compute_covariance from ..source_estimate import _make_stc, SourceEstimate from ..source_space import label_src_vertno_sel @@ -285,6 +285,7 @@ def lcmv(evoked, forward, noise_cov, data_cov, reg=0.01, label=None, beamformers for neuromagnetic source reconstruction. Biomedical Engineering (2004) vol. 51 (10) pp. 1726--34 """ + _check_reference(evoked) info = evoked.info data = evoked.data @@ -348,6 +349,7 @@ def lcmv_epochs(epochs, forward, noise_cov, data_cov, reg=0.01, label=None, beamformers for neuromagnetic source reconstruction. Biomedical Engineering (2004) vol. 51 (10) pp. 1726--34 """ + _check_reference(epochs) info = epochs.info tmin = epochs.times[0] @@ -422,6 +424,7 @@ def lcmv_raw(raw, forward, noise_cov, data_cov, reg=0.01, label=None, beamformers for neuromagnetic source reconstruction. Biomedical Engineering (2004) vol. 51 (10) pp. 1726--34 """ + _check_reference(raw) info = raw.info @@ -610,6 +613,7 @@ def tf_lcmv(epochs, forward, noise_covs, tmin, tmax, tstep, win_lengths, time-frequency dynamics of cortical activity. NeuroImage (2008) vol. 40 (4) pp. 1686-1700 """ + _check_reference(epochs) if pick_ori not in [None, 'normal']: raise ValueError('Unrecognized orientation option in pick_ori, ' diff --git a/mne/inverse_sparse/_gamma_map.py b/mne/inverse_sparse/_gamma_map.py index 01d11e58e7f..eab2e220615 100644 --- a/mne/inverse_sparse/_gamma_map.py +++ b/mne/inverse_sparse/_gamma_map.py @@ -8,7 +8,7 @@ from ..forward import is_fixed_orient, _to_fixed_ori from ..io.pick import pick_channels_evoked -from ..minimum_norm.inverse import _prepare_forward +from ..minimum_norm.inverse import _prepare_forward, _check_reference from ..utils import logger, verbose from .mxne_inverse import _make_sparse_stc, _prepare_gain @@ -231,6 +231,8 @@ def gamma_map(evoked, forward, noise_cov, alpha, loose=0.2, depth=0.8, Wipf et al. A unified Bayesian framework for MEG/EEG source imaging, NeuroImage, vol. 44, no. 3, pp. 947-66, Mar. 2009. """ + _check_reference(evoked) + # make forward solution in fixed orientation if necessary if loose is None and not is_fixed_orient(forward): forward = deepcopy(forward) diff --git a/mne/inverse_sparse/mxne_inverse.py b/mne/inverse_sparse/mxne_inverse.py index 7223096c161..bb9dc7a79e9 100644 --- a/mne/inverse_sparse/mxne_inverse.py +++ b/mne/inverse_sparse/mxne_inverse.py @@ -8,6 +8,7 @@ from ..source_estimate import SourceEstimate from ..minimum_norm.inverse import combine_xyz, _prepare_forward +from ..minimum_norm.inverse import _check_reference from ..forward import compute_orient_prior, is_fixed_orient, _to_fixed_ori from ..io.pick import pick_channels_evoked from .mxne_optim import mixed_norm_solver, norm_l2inf, tf_mixed_norm_solver @@ -159,6 +160,8 @@ def mixed_norm(evoked, forward, noise_cov, alpha, loose=0.2, depth=0.8, if not isinstance(evoked, list): evoked = [evoked] + _check_reference(evoked[0]) + all_ch_names = evoked[0].ch_names if not all(all_ch_names == evoked[i].ch_names for i in range(1, len(evoked))): @@ -361,6 +364,8 @@ def tf_mixed_norm(evoked, forward, noise_cov, alpha_space, alpha_time, The residual a.k.a. data not explained by the sources. Only returned if return_residual is True. """ + _check_reference(evoked) + all_ch_names = evoked.ch_names info = evoked.info diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index e3e3b863a67..609811e1251 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -15,6 +15,7 @@ from ..io.matrix import (_read_named_matrix, _transpose_named_matrix, write_named_matrix) from ..io.proj import _read_proj, make_projector, _write_proj +from ..io.proj import _has_eeg_average_ref_proj from ..io.tree import dir_tree_find from ..io.write import (write_int, write_float_matrix, start_file, start_block, end_block, end_file, write_float, @@ -716,6 +717,15 @@ def _check_ori(pick_ori, pick_normal): return pick_ori +def _check_reference(inst): + if "eeg" in inst and not _has_eeg_average_ref_proj(inst.info['projs']): + raise ValueError('EEG average reference is mandatory for inverse ' + 'modeling.') + if inst.info['custom_ref_applied']: + raise ValueError('Custom EEG reference is not allowed for inverse ' + 'modeling.') + + def _subject_from_inverse(inverse_operator): """Get subject id from inverse operator""" return inverse_operator['src'][0].get('subject_his_id', None) @@ -752,6 +762,7 @@ def apply_inverse(evoked, inverse_operator, lambda2, method="dSPM", stc : SourceEstimate | VolSourceEstimate The source estimates """ + _check_reference(evoked) method = _check_method(method) pick_ori = _check_ori(pick_ori, pick_normal) # @@ -850,6 +861,7 @@ def apply_inverse_raw(raw, inverse_operator, lambda2, method="dSPM", stc : SourceEstimate | VolSourceEstimate The source estimates. """ + _check_reference(raw) method = _check_method(method) pick_ori = _check_ori(pick_ori, pick_normal) @@ -1017,6 +1029,7 @@ def apply_inverse_epochs(epochs, inverse_operator, lambda2, method="dSPM", stc : list of SourceEstimate or VolSourceEstimate The source estimates for all epochs. """ + _check_reference(epochs) stcs = _apply_inverse_epochs_gen(epochs, inverse_operator, lambda2, method=method, label=label, nave=nave, pick_ori=pick_ori, verbose=verbose, diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index 5cc3e9953d8..755963c7a75 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -232,6 +232,13 @@ def test_apply_inverse_operator(): assert_true(stc.data.max() < 35) assert_true(stc.data.mean() > 0.1) + # Test we get errors when using custom ref or no average proj is present + evoked.info['custom_ref_applied'] = True + assert_raises(ValueError, apply_inverse, evoked, inv_op, lambda2, "MNE") + evoked.info['custom_ref_applied'] = False + evoked.info['projs'] = [] # remove EEG proj + assert_raises(ValueError, apply_inverse, evoked, inv_op, lambda2, "MNE") + @testing.requires_testing_data def test_make_inverse_operator_fixed(): From d9406b3d8d3694bdde28ef631b250dec6ad83803 Mon Sep 17 00:00:00 2001 From: Alexandre Gramfort Date: Mon, 29 Dec 2014 12:43:50 +0100 Subject: [PATCH 12/12] address tiny commnents --- mne/io/proj.py | 3 ++- mne/io/reference.py | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/mne/io/proj.py b/mne/io/proj.py index 0eb67a6aef6..b66c35a093a 100644 --- a/mne/io/proj.py +++ b/mne/io/proj.py @@ -17,6 +17,7 @@ from .constants import FIFF from .pick import pick_types from ..utils import logger, verbose +from ..externals.six import string_type class Projection(dict): @@ -196,7 +197,7 @@ def plot_projs_topomap(self, ch_type=None, layout=None): layout = [] if ch_type is None: ch_type = [ch for ch in ['meg', 'eeg'] if ch in self] - elif isinstance(ch_type, str): + elif isinstance(ch_type, string_type): ch_type = [ch_type] for ch in ch_type: if ch in self: diff --git a/mne/io/reference.py b/mne/io/reference.py index 78f0085e8d8..651aeec5255 100644 --- a/mne/io/reference.py +++ b/mne/io/reference.py @@ -30,12 +30,13 @@ def _apply_reference(raw, ref_from, ref_to=None, copy=True): all EEG channels are chosen. copy : bool Specifies whether instance of Raw will be copied or modified in place. + Defaults to True. Returns ------- raw : instance of Raw Instance of Raw with EEG channels rereferenced. - ref_data : array + ref_data : array, shape (n_times,) Array of reference data subtracted from EEG channels. Notes