Source code for medusa.meeg.meeg

# Built-in imports
import math

# External imports
import warnings

import scipy.io
import numpy as np

# Medusa imports
from medusa import components
from medusa import meeg


[docs]class MEEGChannel(components.SerializableComponent): """This class implements a M/EEG channel. """ # TODO: This class is still not being used for compatibility reasons, but it # should be introduced in the next major update of medusa kernel to remove # channel dictionaries and simplify the management of MEEG signals
[docs] def __init__(self, label, coordinates=None, reference=None): """Constructor for class MEEGChannel. Parameters ---------- label: str Label of the channel. coordinates: dict, optional Dict with the coordinates of the electrode. It is strongly recommended to set this optional parameter in order to use advanced features of MEDUSA (e.g., topographic plots). reference: MEEGChannel, optional Use only for bipolar montages to define the reference electrode. """ # Check errors if reference is not None: if not isinstance(reference, MEEGChannel): raise ValueError( 'Parameter reference must be of type MEEGChannel or None') if coordinates is None: warnings.warn('Channel coordinates are used by some advanced ' 'features of MEDUSA (e.g., topographic plots). ' 'Please consider setting the coordinates or use ' 'function MEEGChannel.from_standard_set to load ' 'standard coordinates.') # Set attributes self.label = label self.coordinates = coordinates self.reference = reference
[docs] def to_serializable_obj(self): pass
[docs] @classmethod def from_standard_set(cls, label, standard='10-5'): pass
[docs] @classmethod def from_serializable_obj(cls, data): pass
[docs]class EEGChannelSet(components.SerializableComponent): """Class to represent an EEG montage with ordered channels in specific coordinates. It also provides functionality to load channels from EEG standards 10-20, 10-10 and 10-5 directly from the labels. """
[docs] def __init__(self, reference_method='common', dim='2D', coord_system='spherical'): """Constructor of class EEGChannelSet Parameters ---------- reference_method: str {'common'|'average'|'bipolar'} Reference method. Recordings with common reference are referenced to the same channel (e.g., ear lobe, mastoid). Recordings with average reference are referenced to the average of all or several channels. Finally, bipolar reference is the subtraction of 2 channels. dim: str {'2D'|'3D'} Dimensions of the coordinates plane. coord_system: str {'cartesian'|'spherical'} Coordinates system. Take into account that, if dim = '2D' spherical refers to polar coordinates """ # Check errors if reference_method not in ('common', 'average', 'bipolar'): raise ValueError('Unknown reference method %s' % reference_method) if dim not in ('2D', '3D'): raise ValueError('Unknown number of dimensions %s' % dim) if coord_system not in ('cartesian', 'spherical'): raise ValueError('Unknown coordinates system %s' % coord_system) # Set attributes self.reference_method = reference_method self.dim = dim self.coord_system = coord_system self.channels = None self.n_cha = None self.l_cha = None self.montage = None self.ground = None self.allow_unlocated_channels = False
[docs] def set_ground(self, ground): """Sets the ground of the montage Parameters ---------- ground: dict Dict with the ground data. The easiest way to calculate it is using function get_standard_channel_data_from_label. Keys: - label: channel label - coordinates: depends on the coordinate system and dimensions: - Dim = 2D: - Cartesian coordinates. Ex: {'x': 0.5, 'y': 0} - Spherical coordinates Ex: {'r': 0.5, 'theta': np.pi/2} - Dim = 3D: - Cartesian coordinates. Ex: {'x': 0.5, 'y': 0, 'z': 0.8} - Spherical coordinates. Ex: {'r': 0.5, 'theta': np.pi/2, 'phi': np.pi/4} """ self.ground = ground
[docs] def add_channel(self, channel, reference): """Function to add a channel to the end of the current montage. Take into account that the order of the channels is important! Parameters ---------- channel: dict Dict with the channel data. The easiest way to calculate it is using function get_standard_channel_data_from_label. Keys: - label: channel label - coordinates: depends on the coordinate system and dimensions: - Dim = 2D: - Cartesian coordinates. Ex: {'x': 0.5, 'y': 0} - Spherical coordinates Ex: {'r': 0.5, 'theta': np.pi/2} - Dim = 3D: - Cartesian coordinates. Ex: {'x': 0.5, 'y': 0, 'z': 0.8} - Spherical coordinates. Ex: {'r': 0.5, 'theta': np.pi/2, 'phi': np.pi/4} reference: dict, list of dicts or str, optional For common and bipolar reference modes, reference must be a dict with the label and coordinates of the reference. For average reference, it can be 'all' to indicate a common average reference, or a list of dicts (with label and coordinates) of the averaged reference for each channel See Also -------- get_standard_channel_data_from_label: returns channel data given the channel label and the standard. It can be used to get the reference """ # TODO: check input channel['reference'] = reference channels = list() if self.channels is None else self.channels channels.append(channel) # Check channels self.__check_channels(channels) # Store attributes self.channels = channels self.n_cha = len(self.channels) self.l_cha = [cha['label'] for cha in self.channels]
[docs] def set_montage(self, channels, ground=None, allow_unlocated_channels=False): """Sets a custom montage, overwriting the previous one. Add single channels more easily using function add_custom_channel and add_standard_channel. Parameters ---------- channels : list List of dicts, each of them representing a channel. The dict must contain the label, coordinates according to parameters dim an coord_system, and reference. For common reference mode, the reference must be a single channel with label and coordinates. For average reference mode, reference can be a list of dicts with the channels (label and coordinates) or "all", to specify common average reference. For bipolar reference mode, the reference must be a single channel with label and coordinates. In all cases you can set the reference to None, but this is not recommended. The definition of the coordinates may be skipped depending on the value of allow_unlocated_channels. ground : dict Dict containing the label and coordinates of the ground electrode allow_unlocated_channels: bool If False, the coordinates of all channels must be defined within channels dict. If True, the channels may not have known coordinates, This may be convenient if the localization of a channel is not known or it has a non-standard label, but the behaviour in functions that need coordinates (e.g., topographic_plots) is unpredictable. See Also -------- set_standard_montage: preferred choice in most cases """ # Check errors self.allow_unlocated_channels = allow_unlocated_channels self.__check_channels(channels, ground) # Set attributes self.channels = channels self.ground = ground self.n_cha = len(self.channels) self.l_cha = [cha['label'] for cha in self.channels]
[docs] def set_standard_montage(self, l_cha=None, l_reference=None, l_ground=None, montage='10-05', drop_landmarks=True, allow_unlocated_channels=False): """Set standard EEG channels with common reference. In 3 dimensions, the equator is taken a Nz-T10-Iz-T9. Parameters ---------- l_cha : list, optional List of channels labels. The data will be returned keeping the same order. If None, the channels will be returned in the same order as they appear in the corresponding standard in medusa.meeg l_reference : str, optional Label of the reference. Usual choices are left preauricular point (LPA) and right preauricular point (RPA). Leave to None if you do not want to specify the reference (not recommended). Use only if reference_method is 'common'. l_ground : str, optional Label of the ground. Usual choices are AFz or FPz. montage : str {'10-20'|'10-10'|'10-05'} or dict EEG standard. If its a string, the corresponding labels and locations of the standard channels will be loaded using the files in medusa.meeg. To load a different montage, use function as meeg.read_montage_file and pass the returned dict here. drop_landmarks : bool Drop landmarks: nasion (NAS), left preauricular point (LPA) and right preauricular point (RPA) or mastoids (M1, M2). These are usually used as reference or ground points, so they are usually removed from the channel set for data analysis. Only used if l_cha is None. allow_unlocated_channels: bool If False, all the labels in parameter l_cha must be contained in the montage, which contains the corresponding coordinates. If True, the channels may not be in the standard, and they will be saved with no coordinates. This allows to save labels that are not defined in the montage, but the behaviour in functions that need locations (e.g., topographic_plots) is unpredictable. """ # Check errors if self.reference_method == 'bipolar': raise ValueError('Function set_standard_channels is not available ' 'for recordings with bipolar reference. For ' 'custom montages use set_custom_channels') assert self.dim == '2D' or self.dim == '3D', \ 'Incorrect input on dim parameter' assert self.coord_system == 'cartesian' or \ self.coord_system == 'spherical', \ 'Incorrect input on coord_system parameter' # Get montage if isinstance(montage, str): # Load standard montage self.montage = montage montage = meeg.get_standard_montage(standard=montage, dim=self.dim, coord_system=self.coord_system) else: # Set custom montage montage = montage.copy() self.montage = montage # Reference if l_reference is not None: if l_reference in montage: reference = montage[l_reference] else: if allow_unlocated_channels: reference = dict() warnings.warn('Reference not defined in montage') else: raise meeg.ChannelNotFound(l_reference) reference['label'] = l_reference else: reference = None # Get list of labels to get labels = montage.keys() if l_cha is None \ else [l.upper().strip() for l in l_cha] # Get channels channels = list() for label in labels: # Drop landmarks if l_cha is None: if drop_landmarks: if label in ('NAS', 'LPA', 'RPA', 'M1', 'M2'): continue # Append info if label in montage: channel_data = montage[label] else: if allow_unlocated_channels: channel_data = dict() warnings.warn('Channel %s not defined in montage' % label) else: raise meeg.ChannelNotFound(label) channel_data['label'] = label channel_data['reference'] = reference channels.append(channel_data) # Ground if l_ground is not None: if l_ground in montage: ground = montage[l_ground] else: if allow_unlocated_channels: ground = dict() warnings.warn('Ground not defined in montage') else: raise meeg.ChannelNotFound(l_ground) ground['label'] = ground else: ground = None # Check channels self.set_montage(channels, ground=ground, allow_unlocated_channels=allow_unlocated_channels)
def __check_channels(self, channels, ground=None): # Get mandatory and coordinates keys for each dim and coord_system cha_keys = ['label', 'reference'] ref_keys = ['label'] gnd_keys = ['label'] if self.dim == '2D': if self.coord_system == 'cartesian': coord_keys = ['x', 'y'] else: coord_keys = ['r', 'theta'] else: if self.coord_system == 'cartesian': coord_keys = ['x', 'y', 'z'] else: coord_keys = ['r', 'theta', 'phi'] if not self.allow_unlocated_channels: cha_keys += coord_keys ref_keys += ref_keys gnd_keys += gnd_keys # Check keys references = list() for cha in channels: if not all(k in cha for k in cha_keys): raise ValueError('Malformed channel %s. Dict keys must be %s' % (str(cha), str(cha_keys))) # Check reference reference = cha['reference'] references.append(reference) if reference is not None: if self.reference_method == 'common': # Single reference assert isinstance(reference, dict), \ 'Reference must be None or dict' if not all(k in reference for k in ref_keys): raise ValueError('Malformed reference in channel %s. ' 'Reference keys must be %s' % (str(cha), str(ref_keys))) # All references must be identical if not all(ref == references[0] for ref in references): raise ValueError('All references must be identical ' 'in common reference mode') elif self.reference_method == 'average': # Average reference assert isinstance(reference, list), \ 'Reference must be None or list of dicts' for r in reference: if not all(k in r for k in ref_keys): raise ValueError('Malformed reference in ' 'channel %s. Reference keys ' 'must be %s' % (str(cha), str(ref_keys))) elif self.reference_method == 'bipolar': # Single reference assert isinstance(reference, dict), \ 'Reference must be of type dict' if not all(k in reference for k in ref_keys): raise ValueError('Malformed reference in channel %s. ' 'Reference keys must be %s' % (str(cha), str(ref_keys))) # Check ground if ground is not None: if not all(k in ground for k in gnd_keys): raise ValueError( 'Malformed ground. Dict keys must be %s' % (str(gnd_keys)))
[docs] def get_cha_idx_from_labels(self, labels): """Returns the position of the channels given the labels Parameters ---------- labels : list Labels to check. The order matters Returns ------- indexes : np.ndarray Indexes of the channels in the set """ return [self.l_cha.index(l) for l in labels]
[docs] def check_channels_labels(self, labels, strict=False): """Checks the order and labels of the channels Parameters ---------- labels : list Labels to check. The order matters strict : bool If True, comparison is strict. The function will check that the channel set contains the channels given by parameter labels and in the same order. If false, the function checks that the channels are contained in the channel set, but they could be in different order and the set could contain more channels Returns ------- check : bool True if the labels and order are the same. False otherwise """ if strict: check = True for i in range(len(labels)): if self.l_cha[i] != labels[i]: check = False else: check = True for l in labels: if l not in self.l_cha: check = False return check
[docs] def subset(self, cha_idx): """Selects the channels given the indexes, creating a subset. The order of the channels will be updated Parameters ---------- cha_idx : np.ndarray Indexes of the channels to select. The order matters """ self.channels = [self.channels[idx] for idx in cha_idx] self.n_cha = len(self.channels) self.l_cha = [cha['label'] for cha in self.channels]
[docs] def sort_nearest_channels(self): """ Sorts the nearest channels for each possible channel Return ---------- sorted_dist_ch: dict() Dictionary that includes, for each channel (as a key), the rest of channels sorted by its closeness to that channel. """ # TODO: Añadir la opción de hacer un sort de una coordenada específica, # o en otra función if not self.channels: raise Exception( 'Cannot compute the nearest channels if channel set ' 'is not initialized!') # Compute distance matrix dist_matrix = self.compute_dist_matrix() # Create the dictionary and format it sorted_dist_ch = dict() for i, label in enumerate(self.l_cha): d_sorted = np.sort(dist_matrix[i, :]) idx_sorted = np.argsort(dist_matrix[i, :]) if d_sorted[0] != 0: raise Exception('[meeg/EEGChannelSet/sort_nearest_channels] ' 'Something ocurred, distance between the same ' 'channel is not zero!') sorted_dist_ch[label] = [] for j in range(1, len(d_sorted)): sorted_dist_ch[label].append({ "dist": d_sorted[j], "channel": self.channels[idx_sorted[j]] }) return sorted_dist_ch
[docs] def compute_dist_matrix(self): """This function computes the distances between all channels in the channel set and stores them into a matrix. Returns ------------- dist_matrix: ndarray of dimensions [ncha x ncha] Distances between all the channels. """ if not self.channels: raise Exception( 'Cannot compute the distance matrix if channel set ' 'is not initialized!') # Instantiate matrix of distances dist_matrix = np.empty((self.n_cha, self.n_cha)) # For 2D coordinates if self.dim == '2D': if self.coord_system == 'cartesian': for i, cha in enumerate(self.channels): # Find location of channel cha_pos = np.array([cha['x'], cha['y']]) # Find the location of the rest of the channels for j, temp_cha in enumerate(self.channels): temp_cha_pos = np.array( [temp_cha['x'], temp_cha['y']] ) d = np.sqrt(np.sum(np.power(temp_cha_pos - cha_pos, 2))) dist_matrix[i, j] = d elif self.coord_system == 'spherical': for i, cha in enumerate(self.channels): # Find location of channel r_cha, theta_cha = cha['r'], cha['theta'] # Find the location of the rest of the channels for j, temp_cha in enumerate(self.channels): r_temp_cha, theta_temp_cha = \ temp_cha['r'], temp_cha['theta'] d = np.abs(np.sqrt(r_cha ** 2 + r_temp_cha ** 2 - 2 * r_cha * r_temp_cha * np.cos(theta_temp_cha - theta_cha))) dist_matrix[i, j] = d # For 3D coordinates elif self.dim == '3D': if self.coord_system == 'cartesian': for i, cha in enumerate(self.channels): # Find location of channel cha_pos = np.array([cha['x'], cha['y'], cha['z']]) # Find the location of the rest of the channels for j, temp_cha in enumerate(self.channels): temp_cha_pos = np.array( [temp_cha['x'], temp_cha['y'], temp_cha['z']] ) d = np.sqrt(np.sum(np.power(temp_cha_pos - cha_pos, 2))) dist_matrix[i, j] = d elif self.coord_system == 'spherical': for i, cha in enumerate(self.channels): # Find location of channel r_cha, theta_cha, phi_cha = \ cha['r'], cha['theta'], cha['phi'] # Find the location of the rest of the channels for j, temp_cha in enumerate(self.channels): r_temp_cha, theta_temp_cha, phi_temp_cha = \ temp_cha['r'], temp_cha['theta'], temp_cha['phi'] d = np.abs( np.sqrt( r_cha ** 2 + r_temp_cha ** 2 - 2 * r_cha * r_temp_cha * ( np.sin(theta_cha) * np.sin(theta_temp_cha) * np.cos(phi_cha - phi_temp_cha) + np.cos(theta_temp_cha) * np.cos(theta_cha)) ) ) dist_matrix[i, j] = d return dist_matrix
[docs] def to_serializable_obj(self): return self.__dict__
[docs] @classmethod def from_serializable_obj(cls, dict_data): inst = cls() inst.__dict__.update(dict_data) return inst
[docs]class MEGChannelSet(components.SerializableComponent): # TODO
[docs] def __init__(self): self.channels = None self.n_cha = None self.l_cha = None
[docs] def to_serializable_obj(self): return self.__dict__
[docs] @classmethod def from_serializable_obj(cls, dict_data): inst = cls() inst.__dict__.update(dict_data) return inst
[docs]class EEG(components.BiosignalData): """Electroencephalography (EEG) biosignal """
[docs] def __init__(self, times, signal, fs, channel_set, **kwargs): """EEG constructor Parameters ---------- times : list or numpy.ndarray 1D numpy array [n_samples]. Timestamps of each sample. If they are not available, generate them artificially. Nevertheless, all signals and events must have the same temporal origin signal : list or numpy.ndarray 2D numpy array [n_samples x n_channels]. EEG samples (the units should be defined using kwargs) fs : int or float Sample rate of the recording. channel_set : meeg_standards.EEGChannelSet EEG channel set kwargs: kwargs Any other parameter provided will be saved in the class (e.g., equipment description) """ # To numpy arrays times = np.array(times) signal = np.array(signal) # Check errors if signal.shape[1] != channel_set.n_cha: raise Exception("Signal with shape [samples x channels] does not " "match with the number of channels") if times.shape[0] != signal.shape[0]: raise Exception("Parameters times (shape: %s) and signal (shape: " "%s) must have the same length" % ( str(times.shape), str(signal.shape)) ) # Standard attributes self.times = times self.signal = signal self.fs = fs self.channel_set = channel_set # Optional attributes for key, value in kwargs.items(): setattr(self, key, value)
[docs] def change_channel_set(self, channel_set): """Smart change of channel set, updating the signal and all related attributes Parameters ---------- channel_set : meeg_standards.EEGChannelSet EEG channel set """ # Get the index of the channels cha_idx = self.channel_set.get_cha_idx_from_labels(channel_set.l_cha) # Select and reorganize channels channels self.channel_set.subset(cha_idx) # Reorganize signal self.signal = self.signal[:, cha_idx]
[docs] def to_serializable_obj(self): rec_dict = self.__dict__ for key in rec_dict.keys(): if type(rec_dict[key]) == np.ndarray: rec_dict[key] = rec_dict[key].tolist() if type(rec_dict[key]) == EEGChannelSet: rec_dict[key] = rec_dict[key].to_serializable_obj() return rec_dict
[docs] @classmethod def from_serializable_obj(cls, dict_data): # Load channel set dict dict_data['channel_set'] = EEGChannelSet.from_serializable_obj( dict_data['channel_set'] ) return cls(**dict_data)
[docs]class MEG(components.BiosignalData): # TODO check everything """Magnetoencephalography (MEG) biosignal """
[docs] def __init__(self, times, signal, fs, channel_set, **kwargs): """MEG constructor Parameters ---------- times : list or numpy.ndarray 1D numpy array [n_samples]. Timestamps of each sample. If they are not available, generate them artificially. Nevertheless, all signals and events must have the same temporal origin signal : list or numpy.ndarray 2D numpy array [n_samples x n_channels]. MEG samples (the units should be defined using kwargs) fs : int or float Sample rate of the recording. channel_set : list or meeg_standards.MEGChannelSet MEG channel set. kwargs: kwargs Any other parameter provided will be saved in the class (e.g., equipment description) """ # To numpy arrays times = np.array(times) signal = np.array(signal) # Check errors if signal.shape[1] != channel_set.n_cha: raise Exception("Signal with shape [samples x channels] does not " "match with the number of channels") if times.shape[0] != signal.shape[0]: raise Exception("Parameters times and signal must have the same " "length") # Standard attributes self.times = times self.signal = signal self.fs = fs self.channel_set = channel_set # Optional attributes for key, value in kwargs.items(): setattr(self, key, value)
[docs] def change_channel_set(self, channel_set): """Smart change of channel set, updating the signal and all related attributes Parameters ---------- channel_set : meeg_standards.MEGChannelSet MEG channel set """ raise NotImplementedError
[docs] def to_serializable_obj(self): rec_dict = self.__dict__ for key in rec_dict.keys(): if type(rec_dict[key]) == np.ndarray: rec_dict[key] = rec_dict[key].tolist() return rec_dict
[docs] @classmethod def from_serializable_obj(cls, dict_data): return cls(**dict_data)
[docs] @staticmethod def load_meg_signal_from_spm_file(path, fs): # TODO check everything """Function to load a MEG recording from a spm file Parameters ---------- path : str Path of the file fs : int or float Sample rate of the recording """ data = scipy.io.loadmat(path) info = data['D'][0, 0] datatype = info['datatype'][0] num_chan = np.size(info['channels']) num_samples = info['Nsamples'][0, 0] raw = np.fromfile(path[0:len(path) - 3] + 'dat', datatype[0:len(datatype) - 3]) signal = np.reshape(raw, [num_chan, num_samples], order='F') times = np.linspace(0, signal[0] / fs, signal[0]) channels = None return MEG(times, signal, fs, channels)
[docs]class Connecitivity: # TODO: check everything """Customizable class with connectivity info from EEG/MEG recordings """
[docs] def __init__(self, data, trial_len, parameter, filt_mode, **kwargs): """Class constructor Parameters ---------- data : bla bla Bla bla trial_len : bla bla Bla bla parameter : bla bla Bla bla filt_mode : bla bla Bla bla kwargs Optional information of the EEG recording (e.g. subject, amplifier, etc) """ # Params data = np.array(data) if not (filt_mode == 'all' or filt_mode == 'bands' or filt_mode == 'win'): raise ValueError("Unknown filtering mode") self.data = data self.trial_len = trial_len self.parameter = parameter self.filt_mode = filt_mode # Optional attributes for key, value in kwargs.items(): setattr(self, key, value)
[docs]class UnknownStandardChannel(Exception):
[docs] def __init__(self, msg=None): """Class constructor Parameters ---------- msg: string or None Custom message """ if msg is None: msg = 'Unknown standard channel' super().__init__(msg)