Aluode/PerceptionLabPortable
0
1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5import logging6from collections import defaultdict7from copy import deepcopy8from itertools import combinations9from pathlib import Path10 11import numpy as np12from scipy.spatial.distance import pdist, squareform13 14from .._fiff.constants import FIFF15from .._fiff.meas_info import Info16from .._fiff.pick import _FNIRS_CH_TYPES_SPLIT, _picks_to_idx, pick_types17from ..transforms import _cart_to_sph, _pol_to_cart18from ..utils import (19 _check_ch_locs,20 _check_fname,21 _check_option,22 _check_sphere,23 _clean_names,24 _ensure_int,25 fill_doc,26 logger,27 verbose,28 warn,29)30from ..viz.topomap import plot_layout31from .channels import _get_ch_info32 33 34class Layout:35 """Sensor layouts.36 37 Layouts are typically loaded from a file using38 :func:`~mne.channels.read_layout`. Only use this class directly if you're39 constructing a new layout.40 41 Parameters42 ----------43 box : tuple of length 444 The box dimension (x_min, x_max, y_min, y_max).45 pos : array, shape=(n_channels, 4)46 The unit-normalized positions of the channels in 2d47 (x, y, width, height).48 names : list of str49 The channel names.50 ids : array-like of int51 The channel ids.52 kind : str53 The type of Layout (e.g. 'Vectorview-all').54 """55 56 def __init__(self, box, pos, names, ids, kind):57 self.box = box58 self.pos = pos59 self.names = names60 self.ids = np.array(ids)61 if self.ids.ndim != 1:62 raise ValueError("The channel indices should be a 1D array-like.")63 self.kind = kind64 65 def copy(self):66 """Return a copy of the layout.67 68 Returns69 -------70 layout : instance of Layout71 A deepcopy of the layout.72 73 Notes74 -----75 .. versionadded:: 1.776 """77 return deepcopy(self)78 79 def save(self, fname, overwrite=False):80 """Save Layout to disk.81 82 Parameters83 ----------84 fname : path-like85 The file name (must end with either ``.lout`` or ``.lay``).86 overwrite : bool87 If True, overwrites the destination file if it exists.88 89 See Also90 --------91 read_layout92 """93 x = self.pos[:, 0]94 y = self.pos[:, 1]95 width = self.pos[:, 2]96 height = self.pos[:, 3]97 fname = _check_fname(fname, overwrite=overwrite, name=fname)98 if fname.suffix == ".lout":99 out_str = "{:8.2f} {:8.2f} {:8.2f} {:8.2f}\n".format(*self.box)100 elif fname.suffix == ".lay":101 out_str = ""102 else:103 raise ValueError("Unknown layout type. Should be of type .lout or .lay.")104 105 for ii in range(x.shape[0]):106 out_str += (107 f"{self.ids[ii]:03d} {x[ii]:8.2f} {y[ii]:8.2f} "108 f"{width[ii]:8.2f} {height[ii]:8.2f} {self.names[ii]}\n"109 )110 111 f = open(fname, "w")112 f.write(out_str)113 f.close()114 115 def __repr__(self):116 """Return the string representation."""117 return "<Layout | {} - Channels: {} ...>".format(118 self.kind,119 ", ".join(self.names[:3]),120 )121 122 @fill_doc123 def plot(self, picks=None, show_axes=False, show=True):124 """Plot the sensor positions.125 126 Parameters127 ----------128 %(picks_nostr)s129 show_axes : bool130 Show layout axes if True. Defaults to False.131 show : bool132 Show figure if True. Defaults to True.133 134 Returns135 -------136 fig : instance of matplotlib.figure.Figure137 Figure containing the sensor topography.138 139 Notes140 -----141 .. versionadded:: 0.12.0142 """143 return plot_layout(self, picks=picks, show_axes=show_axes, show=show)144 145 @verbose146 def pick(self, picks=None, exclude=(), *, verbose=None):147 """Pick a subset of channels.148 149 Parameters150 ----------151 %(picks_layout)s152 exclude : str | int | array-like of str or int153 Set of channels to exclude, only used when ``picks`` is set to ``'all'`` or154 ``None``. Exclude will not drop channels explicitly provided in ``picks``.155 %(verbose)s156 157 Returns158 -------159 layout : instance of Layout160 The modified layout.161 162 Notes163 -----164 .. versionadded:: 1.7165 """166 # TODO: all the picking functions operates on an 'info' object which is missing167 # for a layout, thus we have to do the extra work here. The logic below can be168 # replaced when https://github.com/mne-tools/mne-python/issues/11913 is solved.169 if (isinstance(picks, str) and picks == "all") or (picks is None):170 picks = deepcopy(self.names)171 apply_exclude = True172 elif isinstance(picks, str):173 picks = [picks]174 apply_exclude = False175 elif isinstance(picks, slice):176 try:177 picks = np.arange(len(self.names))[picks]178 except TypeError:179 raise TypeError(180 "If a slice is provided, it must be a slice of integers."181 )182 apply_exclude = False183 else:184 try:185 picks = [_ensure_int(picks)]186 except TypeError:187 picks = (188 list(picks) if isinstance(picks, tuple | set) else deepcopy(picks)189 )190 apply_exclude = False191 if apply_exclude:192 if isinstance(exclude, str):193 exclude = [exclude]194 else:195 try:196 exclude = [_ensure_int(exclude)]197 except TypeError:198 exclude = (199 list(exclude)200 if isinstance(exclude, tuple | set)201 else deepcopy(exclude)202 )203 for var, var_name in ((picks, "picks"), (exclude, "exclude")):204 if var_name == "exclude" and not apply_exclude:205 continue206 if not isinstance(var, list | tuple | set | np.ndarray):207 raise TypeError(208 f"'{var_name}' must be a list, tuple, set or ndarray. "209 f"Got {type(var)} instead."210 )211 if isinstance(var, np.ndarray) and var.ndim != 1:212 raise ValueError(213 f"'{var_name}' must be a 1D array-like. Got {var.ndim}D instead."214 )215 for k, elt in enumerate(var):216 if isinstance(elt, str) and elt in self.names:217 var[k] = self.names.index(elt)218 continue219 elif isinstance(elt, str):220 raise ValueError(221 f"The channel name {elt} provided in {var_name} does not match "222 "any channels from the layout."223 )224 try:225 var[k] = _ensure_int(elt)226 except TypeError:227 raise TypeError(228 f"All elements in '{var_name}' must be integers or strings."229 )230 if not (0 <= var[k] < len(self.names)):231 raise ValueError(232 f"The value {elt} provided in {var_name} does not match any "233 f"channels from the layout. The layout has {len(self.names)} "234 "channels."235 )236 if len(var) != len(set(var)):237 warn(238 f"The provided '{var_name}' has duplicates which will be ignored.",239 RuntimeWarning,240 )241 picks = picks.astype(int) if isinstance(picks, np.ndarray) else picks242 exclude = exclude.astype(int) if isinstance(exclude, np.ndarray) else exclude243 if apply_exclude:244 picks = np.array(list(set(picks) - set(exclude)), dtype=int)245 if len(picks) == 0:246 raise RuntimeError(247 "The channel selection yielded no remaining channels. Please edit "248 "the arguments 'picks' and 'exclude' to include at least one "249 "channel."250 )251 else:252 picks = np.array(list(set(picks)), dtype=int)253 self.pos = self.pos[picks]254 self.ids = self.ids[picks]255 self.names = [self.names[k] for k in picks]256 return self257 258 259def _read_lout(fname):260 """Aux function."""261 with open(fname) as f:262 box_line = f.readline() # first line contains box dimension263 box = tuple(map(float, box_line.split()))264 names, pos, ids = [], [], []265 for line in f:266 splits = line.split()267 if len(splits) == 7:268 cid, x, y, dx, dy, chkind, nb = splits269 name = chkind + " " + nb270 else:271 cid, x, y, dx, dy, name = splits272 pos.append(np.array([x, y, dx, dy], dtype=np.float64))273 names.append(name)274 ids.append(int(cid))275 276 pos = np.array(pos)277 278 return box, pos, names, ids279 280 281def _read_lay(fname):282 """Aux function."""283 with open(fname) as f:284 box = None285 names, pos, ids = [], [], []286 for line in f:287 splits = line.split()288 if len(splits) == 7:289 cid, x, y, dx, dy, chkind, nb = splits290 name = chkind + " " + nb291 else:292 cid, x, y, dx, dy, name = splits293 pos.append(np.array([x, y, dx, dy], dtype=np.float64))294 names.append(name)295 ids.append(int(cid))296 297 pos = np.array(pos)298 299 return box, pos, names, ids300 301 302def read_layout(fname=None, *, scale=True):303 """Read layout from a file.304 305 Parameters306 ----------307 fname : path-like | str308 Either the path to a ``.lout`` or ``.lay`` file or the name of a309 built-in layout. See Notes for a list of the available built-in310 layouts.311 scale : bool312 Apply useful scaling for out the box plotting using ``layout.pos``.313 Defaults to True.314 315 Returns316 -------317 layout : instance of Layout318 The layout.319 320 See Also321 --------322 Layout.save323 324 Notes325 -----326 Valid ``fname`` arguments are:327 328 .. table::329 :widths: auto330 331 +----------------------+332 | Kind |333 +======================+334 | biosemi |335 +----------------------+336 | CTF151 |337 +----------------------+338 | CTF275 |339 +----------------------+340 | CTF-275 |341 +----------------------+342 | EEG1005 |343 +----------------------+344 | EGI256 |345 +----------------------+346 | GeodesicHeadWeb-130 |347 +----------------------+348 | GeodesicHeadWeb-280 |349 +----------------------+350 | KIT-125 |351 +----------------------+352 | KIT-157 |353 +----------------------+354 | KIT-160 |355 +----------------------+356 | KIT-AD |357 +----------------------+358 | KIT-AS-2008 |359 +----------------------+360 | KIT-UMD-3 |361 +----------------------+362 | magnesWH3600 |363 +----------------------+364 | Neuromag_122 |365 +----------------------+366 | Vectorview-all |367 +----------------------+368 | Vectorview-grad |369 +----------------------+370 | Vectorview-grad_norm |371 +----------------------+372 | Vectorview-mag |373 +----------------------+374 """375 readers = {".lout": _read_lout, ".lay": _read_lay}376 377 if isinstance(fname, str):378 # is it a built-in layout?379 directory = Path(__file__).parent / "data" / "layouts"380 for suffix in ("", ".lout", ".lay"):381 _fname = (directory / fname).with_suffix(suffix)382 if _fname.exists():383 fname = _fname384 break385 # if not, it must be a valid path provided as str or Path386 fname = _check_fname(fname, "read", must_exist=True, name="layout")387 # and it must have a valid extension388 _check_option("fname extension", fname.suffix, readers)389 kind = fname.stem390 box, pos, names, ids = readers[fname.suffix](fname)391 392 if scale:393 pos[:, 0] -= np.min(pos[:, 0])394 pos[:, 1] -= np.min(pos[:, 1])395 scaling = max(np.max(pos[:, 0]), np.max(pos[:, 1])) + pos[0, 2]396 pos /= scaling397 pos[:, :2] += 0.03398 pos[:, :2] *= 0.97 / 1.03399 pos[:, 2:] *= 0.94400 401 return Layout(box=box, pos=pos, names=names, kind=kind, ids=ids)402 403 404@fill_doc405def make_eeg_layout(406 info, radius=0.5, width=None, height=None, exclude="bads", csd=False407):408 """Make a Layout object based on EEG electrode digitization.409 410 Parameters411 ----------412 %(info_not_none)s413 radius : float414 Viewport radius as a fraction of main figure height. Defaults to 0.5.415 width : float | None416 Width of sensor axes as a fraction of main figure height. By default,417 this will be the maximum width possible without axes overlapping.418 height : float | None419 Height of sensor axes as a fraction of main figure height. By default,420 this will be the maximum height possible without axes overlapping.421 exclude : list of str | str422 List of channels to exclude. If empty do not exclude any.423 If 'bads', exclude channels in info['bads'] (default).424 csd : bool425 Whether the channels contain current-source-density-transformed data.426 427 Returns428 -------429 layout : Layout430 The generated Layout.431 432 See Also433 --------434 make_grid_layout, generate_2d_layout435 """436 if not (0 <= radius <= 0.5):437 raise ValueError("The radius parameter should be between 0 and 0.5.")438 if width is not None and not (0 <= width <= 1.0):439 raise ValueError("The width parameter should be between 0 and 1.")440 if height is not None and not (0 <= height <= 1.0):441 raise ValueError("The height parameter should be between 0 and 1.")442 443 pick_kwargs = dict(meg=False, eeg=True, ref_meg=False, exclude=exclude)444 if csd:445 pick_kwargs.update(csd=True, eeg=False)446 picks = pick_types(info, **pick_kwargs)447 loc2d = _find_topomap_coords(info, picks)448 names = [info["chs"][i]["ch_name"] for i in picks]449 450 # Scale [x, y] to be in the range [-0.5, 0.5]451 # Don't mess with the origin or aspect ratio452 scale = np.maximum(-np.min(loc2d, axis=0), np.max(loc2d, axis=0)).max() * 2453 loc2d /= scale454 455 # If no width or height specified, calculate the maximum value possible456 # without axes overlapping.457 if width is None or height is None:458 width, height = _box_size(loc2d, width, height, padding=0.1)459 460 # Scale to viewport radius461 loc2d *= 2 * radius462 463 # Some subplot centers will be at the figure edge. Shrink everything so it464 # fits in the figure.465 scaling = min(1 / (1.0 + width), 1 / (1.0 + height))466 loc2d *= scaling467 width *= scaling468 height *= scaling469 470 # Shift to center471 loc2d += 0.5472 473 n_channels = loc2d.shape[0]474 pos = np.c_[475 loc2d[:, 0] - 0.5 * width,476 loc2d[:, 1] - 0.5 * height,477 width * np.ones(n_channels),478 height * np.ones(n_channels),479 ]480 481 box = (0, 1, 0, 1)482 ids = 1 + np.arange(n_channels)483 layout = Layout(box=box, pos=pos, names=names, kind="EEG", ids=ids)484 return layout485 486 487@fill_doc488def make_grid_layout(info, picks=None, n_col=None):489 """Make a grid Layout object.490 491 This can be helpful to plot custom data such as ICA sources.492 493 Parameters494 ----------495 %(info_not_none)s496 %(picks_base)s all good misc channels.497 n_col : int | None498 Number of columns to generate. If None, a square grid will be produced.499 500 Returns501 -------502 layout : Layout503 The generated layout.504 505 See Also506 --------507 make_eeg_layout, generate_2d_layout508 """509 picks = _picks_to_idx(info, picks, "misc")510 511 names = [info["chs"][k]["ch_name"] for k in picks]512 513 if not names:514 raise ValueError("No misc data channels found.")515 516 ids = list(range(len(picks)))517 size = len(picks)518 519 if n_col is None:520 # prepare square-like layout521 n_row = n_col = np.sqrt(size) # try square522 if n_col % 1:523 # try n * (n-1) rectangle524 n_col, n_row = int(n_col + 1), int(n_row)525 526 if n_col * n_row < size: # jump to the next full square527 n_row += 1528 else:529 n_row = int(np.ceil(size / float(n_col)))530 531 # setup position grid532 x, y = np.meshgrid(np.linspace(-0.5, 0.5, n_col), np.linspace(-0.5, 0.5, n_row))533 x, y = x.ravel()[:size], y.ravel()[:size]534 width, height = _box_size(np.c_[x, y], padding=0.1)535 536 # Some axes will be at the figure edge. Shrink everything so it fits in the537 # figure. Add 0.01 border around everything538 border_x, border_y = (0.01, 0.01)539 x_scaling = 1 / (1.0 + width + border_x)540 y_scaling = 1 / (1.0 + height + border_y)541 x = x * x_scaling542 y = y * y_scaling543 width *= x_scaling544 height *= y_scaling545 546 # Shift to center547 x += 0.5548 y += 0.5549 550 # calculate pos551 pos = np.c_[552 x - 0.5 * width, y - 0.5 * height, width * np.ones(size), height * np.ones(size)553 ]554 box = (0, 1, 0, 1)555 556 layout = Layout(box=box, pos=pos, names=names, kind="grid-misc", ids=ids)557 return layout558 559 560@fill_doc561def find_layout(info, ch_type=None, exclude="bads"):562 """Choose a layout based on the channels in the info 'chs' field.563 564 Parameters565 ----------566 %(info_not_none)s567 ch_type : {'mag', 'grad', 'meg', 'eeg'} | None568 The channel type for selecting single channel layouts.569 Defaults to None. Note, this argument will only be considered for570 VectorView type layout. Use ``'meg'`` to force using the full layout571 in situations where the info does only contain one sensor type.572 exclude : list of str | str573 List of channels to exclude. If empty do not exclude any.574 If 'bads', exclude channels in info['bads'] (default).575 576 Returns577 -------578 layout : Layout instance | None579 None if layout not found.580 """581 _check_option("ch_type", ch_type, [None, "mag", "grad", "meg", "eeg", "csd"])582 583 (584 has_vv_mag,585 has_vv_grad,586 is_old_vv,587 has_4D_mag,588 ctf_other_types,589 has_CTF_grad,590 n_kit_grads,591 has_any_meg,592 has_eeg_coils,593 has_eeg_coils_and_meg,594 has_eeg_coils_only,595 has_neuromag_122_grad,596 has_csd_coils,597 ) = _get_ch_info(info)598 has_vv_meg = has_vv_mag and has_vv_grad599 has_vv_only_mag = has_vv_mag and not has_vv_grad600 has_vv_only_grad = has_vv_grad and not has_vv_mag601 if ch_type == "meg" and not has_any_meg:602 raise RuntimeError("No MEG channels present. Cannot find MEG layout.")603 604 if ch_type == "eeg" and not has_eeg_coils:605 raise RuntimeError("No EEG channels present. Cannot find EEG layout.")606 607 layout_name = None608 if (has_vv_meg and ch_type is None) or (609 any([has_vv_mag, has_vv_grad]) and ch_type == "meg"610 ):611 layout_name = "Vectorview-all"612 elif has_vv_only_mag or (has_vv_meg and ch_type == "mag"):613 layout_name = "Vectorview-mag"614 elif has_vv_only_grad or (has_vv_meg and ch_type == "grad"):615 if info["ch_names"][0].endswith("X"):616 layout_name = "Vectorview-grad_norm"617 else:618 layout_name = "Vectorview-grad"619 elif has_neuromag_122_grad:620 layout_name = "Neuromag_122"621 elif (has_eeg_coils_only and ch_type in [None, "eeg"]) or (622 has_eeg_coils_and_meg and ch_type == "eeg"623 ):624 if not isinstance(info, dict | Info):625 raise RuntimeError(626 "Cannot make EEG layout, no measurement info "627 "was passed to `find_layout`"628 )629 return make_eeg_layout(info, exclude=exclude)630 elif has_csd_coils and ch_type in [None, "csd"]:631 return make_eeg_layout(info, exclude=exclude, csd=True)632 elif has_4D_mag:633 layout_name = "magnesWH3600"634 elif has_CTF_grad:635 layout_name = "CTF-275"636 elif n_kit_grads > 0:637 layout_name = _find_kit_layout(info, n_kit_grads)638 639 # If no known layout is found, fall back on automatic layout640 if layout_name is None:641 picks = _picks_to_idx(info, "data", exclude=(), with_ref_meg=False)642 ch_names = [info["ch_names"][pick] for pick in picks]643 xy = _find_topomap_coords(info, picks=picks, ignore_overlap=True)644 return generate_2d_layout(xy, ch_names=ch_names, name="custom", normalize=True)645 646 layout = read_layout(fname=layout_name)647 if not is_old_vv:648 layout.names = _clean_names(layout.names, remove_whitespace=True)649 if has_CTF_grad:650 layout.names = _clean_names(layout.names, before_dash=True)651 652 # Apply mask for excluded channels.653 if exclude == "bads":654 exclude = info["bads"]655 idx = [ii for ii, name in enumerate(layout.names) if name not in exclude]656 layout.names = [layout.names[ii] for ii in idx]657 layout.pos = layout.pos[idx]658 layout.ids = layout.ids[idx]659 660 return layout661 662 663@fill_doc664def _find_kit_layout(info, n_grads):665 """Determine the KIT layout.666 667 Parameters668 ----------669 %(info_not_none)s670 n_grads : int671 Number of KIT-gradiometers in the info.672 673 Returns674 -------675 kit_layout : str | None676 String naming the detected KIT layout or ``None`` if layout is missing.677 """678 from ..io.kit.constants import KIT_LAYOUT679 680 if info["kit_system_id"] is not None:681 # avoid circular import682 return KIT_LAYOUT.get(info["kit_system_id"])683 elif n_grads == 160:684 return "KIT-160"685 elif n_grads == 125:686 return "KIT-125"687 elif n_grads > 157:688 return "KIT-AD"689 690 # channels which are on the left hemisphere for NY and right for UMD691 test_chs = (692 "MEG 13",693 "MEG 14",694 "MEG 15",695 "MEG 16",696 "MEG 25",697 "MEG 26",698 "MEG 27",699 "MEG 28",700 "MEG 29",701 "MEG 30",702 "MEG 31",703 "MEG 32",704 "MEG 57",705 "MEG 60",706 "MEG 61",707 "MEG 62",708 "MEG 63",709 "MEG 64",710 "MEG 73",711 "MEG 90",712 "MEG 93",713 "MEG 95",714 "MEG 96",715 "MEG 105",716 "MEG 112",717 "MEG 120",718 "MEG 121",719 "MEG 122",720 "MEG 123",721 "MEG 124",722 "MEG 125",723 "MEG 126",724 "MEG 142",725 "MEG 144",726 "MEG 153",727 "MEG 154",728 "MEG 155",729 "MEG 156",730 )731 x = [ch["loc"][0] < 0 for ch in info["chs"] if ch["ch_name"] in test_chs]732 if np.all(x):733 return "KIT-157" # KIT-NY734 elif np.all(np.invert(x)):735 raise NotImplementedError(736 "Guessing sensor layout for legacy UMD "737 "files is not implemented. Please convert "738 "your files using MNE-Python 0.13 or "739 "higher."740 )741 else:742 raise RuntimeError("KIT system could not be determined for data")743 744 745def _box_size(points, width=None, height=None, padding=0.0):746 """Given a series of points, calculate an appropriate box size.747 748 Parameters749 ----------750 points : array, shape (n_points, 2)751 The centers of the axes as a list of (x, y) coordinate pairs. Normally752 these are points in the range [0, 1] centered at 0.5.753 width : float | None754 An optional box width to enforce. When set, only the box height will be755 calculated by the function.756 height : float | None757 An optional box height to enforce. When set, only the box width will be758 calculated by the function.759 padding : float760 Portion of the box to reserve for padding. The value can range between761 0.0 (boxes will touch, default) to 1.0 (boxes consist of only padding).762 763 Returns764 -------765 width : float766 Width of the box767 height : float768 Height of the box769 """770 771 def xdiff(a, b):772 return np.abs(a[0] - b[0])773 774 def ydiff(a, b):775 return np.abs(a[1] - b[1])776 777 points = np.asarray(points)778 all_combinations = list(combinations(points, 2))779 780 if width is None and height is None:781 if len(points) <= 1:782 # Trivial case first783 width = 1.0784 height = 1.0785 else:786 # Find the closest two points A and B.787 a, b = all_combinations[np.argmin(pdist(points))]788 789 # The closest points define either the max width or max height.790 w, h = xdiff(a, b), ydiff(a, b)791 if w > h:792 width = w793 else:794 height = h795 796 # At this point, either width or height is known, or both are known.797 if height is None:798 # Find all axes that could potentially overlap horizontally.799 hdist = pdist(points, xdiff)800 candidates = [all_combinations[i] for i, d in enumerate(hdist) if d < width]801 802 if len(candidates) == 0:803 # No axes overlap, take all the height you want.804 height = 1.0805 else:806 # Find an appropriate height so all none of the found axes will807 # overlap.808 height = np.min([ydiff(*c) for c in candidates])809 810 elif width is None:811 # Find all axes that could potentially overlap vertically.812 vdist = pdist(points, ydiff)813 candidates = [all_combinations[i] for i, d in enumerate(vdist) if d < height]814 815 if len(candidates) == 0:816 # No axes overlap, take all the width you want.817 width = 1.0818 else:819 # Find an appropriate width so all none of the found axes will820 # overlap.821 width = np.min([xdiff(*c) for c in candidates])822 823 # Add a bit of padding between boxes824 width *= 1 - padding825 height *= 1 - padding826 827 return width, height828 829 830@fill_doc831def _find_topomap_coords(832 info, picks, layout=None, ignore_overlap=False, to_sphere=True, sphere=None833):834 """Guess the E/MEG layout and return appropriate topomap coordinates.835 836 Parameters837 ----------838 %(info_not_none)s839 picks : str | list | slice | None840 None will choose all channels.841 layout : None | instance of Layout842 Enforce using a specific layout. With None, a new map is generated843 and a layout is chosen based on the channels in the picks844 parameter.845 sphere : array-like | str846 Definition of the head sphere.847 848 Returns849 -------850 coords : array, shape = (n_chs, 2)851 2 dimensional coordinates for each sensor for a topomap plot.852 """853 picks = _picks_to_idx(info, picks, "all", exclude=(), allow_empty=False)854 855 if layout is not None:856 chs = [info["chs"][i] for i in picks]857 pos = [layout.pos[layout.names.index(ch["ch_name"])] for ch in chs]858 pos = np.asarray(pos)859 else:860 pos = _auto_topomap_coords(861 info,862 picks,863 ignore_overlap=ignore_overlap,864 to_sphere=to_sphere,865 sphere=sphere,866 )867 868 return pos869 870 871@fill_doc872def _auto_topomap_coords(info, picks, ignore_overlap, to_sphere, sphere):873 """Make a 2 dimensional sensor map from sensor positions in an info dict.874 875 The default is to use the electrode locations. The fallback option is to876 attempt using digitization points of kind FIFFV_POINT_EEG. This only works877 with EEG and requires an equal number of digitization points and sensors.878 879 Parameters880 ----------881 %(info_not_none)s882 picks : list | str | slice | None883 None will pick all channels.884 ignore_overlap : bool885 Whether to ignore overlapping positions in the layout. If False and886 positions overlap, an error is thrown.887 to_sphere : bool888 If True, the radial distance of spherical coordinates is ignored, in889 effect fitting the xyz-coordinates to a sphere.890 sphere : array-like | str891 The head sphere definition.892 893 Returns894 -------895 locs : array, shape = (n_sensors, 2)896 An array of positions of the 2 dimensional map.897 """898 sphere = _check_sphere(sphere, info)899 logger.debug(f"Generating coords using: {sphere}")900 901 picks = _picks_to_idx(info, picks, "all", exclude=(), allow_empty=False)902 chs = [info["chs"][i] for i in picks]903 904 # Use channel locations if available905 locs3d = np.array([ch["loc"][:3] for ch in chs])906 907 # If electrode locations are not available, use digitization points908 if not _check_ch_locs(info=info, picks=picks):909 logging.warning(910 "Did not find any electrode locations (in the info "911 "object), will attempt to use digitization points "912 "instead. However, if digitization points do not "913 "correspond to the EEG electrodes, this will lead to "914 "bad results. Please verify that the sensor locations "915 "in the plot are accurate."916 )917 918 # MEG/EOG/ECG sensors don't have digitization points; all requested919 # channels must be EEG920 for ch in chs:921 if ch["kind"] != FIFF.FIFFV_EEG_CH:922 raise ValueError(923 "Cannot determine location of MEG/EOG/ECG "924 "channels using digitization points."925 )926 927 eeg_ch_names = [928 ch["ch_name"] for ch in info["chs"] if ch["kind"] == FIFF.FIFFV_EEG_CH929 ]930 931 # Get EEG digitization points932 if info["dig"] is None or len(info["dig"]) == 0:933 raise RuntimeError("No digitization points found.")934 935 locs3d = np.array(936 [937 point["r"]938 for point in info["dig"]939 if point["kind"] == FIFF.FIFFV_POINT_EEG940 ]941 )942 943 if len(locs3d) == 0:944 raise RuntimeError(945 "Did not find any digitization points of "946 f"kind {FIFF.FIFFV_POINT_EEG} in the info."947 )948 949 if len(locs3d) != len(eeg_ch_names):950 raise ValueError(951 f"Number of EEG digitization points ({len(locs3d)}) doesn't match the "952 f"number of EEG channels ({len(eeg_ch_names)})"953 )954 955 # We no longer center digitization points on head origin, as we work956 # in head coordinates always957 958 # Match the digitization points with the requested959 # channels.960 eeg_ch_locs = dict(zip(eeg_ch_names, locs3d))961 locs3d = np.array([eeg_ch_locs[ch["ch_name"]] for ch in chs])962 963 # Sometimes we can get nans964 locs3d[~np.isfinite(locs3d)] = 0.0965 966 # Duplicate points cause all kinds of trouble during visualization967 dist = pdist(locs3d)968 if len(locs3d) > 1 and np.min(dist) < 1e-10 and not ignore_overlap:969 problematic_electrodes = [970 chs[elec_i]["ch_name"]971 for elec_i in squareform(dist < 1e-10).any(axis=0).nonzero()[0]972 ]973 974 raise ValueError(975 "The following electrodes have overlapping positions,"976 " which causes problems during visualization:\n"977 + ", ".join(problematic_electrodes)978 )979 980 if to_sphere:981 # translate to sphere origin, transform/flatten Z, translate back982 locs3d -= sphere[:3]983 # use spherical (theta, pol) as (r, theta) for polar->cartesian984 cart_coords = _cart_to_sph(locs3d)985 out = _pol_to_cart(cart_coords[:, 1:][:, ::-1])986 # scale from radians to mm987 out *= cart_coords[:, [0]] / (np.pi / 2.0)988 out += sphere[:2]989 else:990 out = _pol_to_cart(_cart_to_sph(locs3d))991 return out992 993 994def _topo_to_sphere(pos, eegs):995 """Transform xy-coordinates to sphere.996 997 Parameters998 ----------999 pos : array-like, shape (n_channels, 2)1000 xy-oordinates to transform.1001 eegs : list of int1002 Indices of EEG channels that are included when calculating the sphere.1003 1004 Returns1005 -------1006 coords : array, shape (n_channels, 3)1007 xyz-coordinates.1008 """1009 xs, ys = np.array(pos).T1010 1011 sqs = np.max(np.sqrt((xs[eegs] ** 2) + (ys[eegs] ** 2)))1012 xs /= sqs # Shape to a sphere and normalize1013 ys /= sqs1014 1015 xs += 0.5 - np.mean(xs[eegs]) # Center the points1016 ys += 0.5 - np.mean(ys[eegs])1017 1018 xs = xs * 2.0 - 1.0 # Values ranging from -1 to 11019 ys = ys * 2.0 - 1.01020 1021 rs = np.clip(np.sqrt(xs**2 + ys**2), 0.0, 1.0)1022 alphas = np.arccos(rs)1023 zs = np.sin(alphas)1024 return np.column_stack([xs, ys, zs])1025 1026 1027@fill_doc1028def _pair_grad_sensors(1029 info, layout=None, topomap_coords=True, exclude="bads", raise_error=True1030):1031 """Find the picks for pairing grad channels.1032 1033 Parameters1034 ----------1035 %(info_not_none)s1036 layout : Layout | None1037 The layout if available. Defaults to None.1038 topomap_coords : bool1039 Return the coordinates for a topomap plot along with the picks. If1040 False, only picks are returned. Defaults to True.1041 exclude : list of str | str1042 List of channels to exclude. If empty, do not exclude any.1043 If 'bads', exclude channels in info['bads']. Defaults to 'bads'.1044 raise_error : bool1045 Whether to raise an error when no pairs are found. If False, raises a1046 warning.1047 1048 Returns1049 -------1050 picks : array of int1051 Picks for the grad channels, ordered in pairs.1052 coords : array, shape = (n_grad_channels, 3)1053 Coordinates for a topomap plot (optional, only returned if1054 topomap_coords == True).1055 """1056 # find all complete pairs of grad channels1057 pairs = defaultdict(list)1058 grad_picks = pick_types(info, meg="grad", ref_meg=False, exclude=exclude)1059 1060 _, has_vv_grad, *_, has_neuromag_122_grad, _ = _get_ch_info(info)1061 1062 for i in grad_picks:1063 ch = info["chs"][i]1064 name = ch["ch_name"]1065 if has_vv_grad and name.startswith("MEG"):1066 if name.endswith(("2", "3")):1067 key = name[-4:-1]1068 pairs[key].append(ch)1069 if has_neuromag_122_grad and name.startswith("MEG"):1070 key = (int(name[-3:]) - 1) // 21071 pairs[key].append(ch)1072 1073 pairs = [p for p in pairs.values() if len(p) == 2]1074 if len(pairs) == 0:1075 if raise_error:1076 raise ValueError("No 'grad' channel pairs found.")1077 else:1078 warn("No 'grad' channel pairs found.")1079 return list()1080 1081 # find the picks corresponding to the grad channels1082 grad_chs = sum(pairs, [])1083 ch_names = info["ch_names"]1084 picks = [ch_names.index(c["ch_name"]) for c in grad_chs]1085 1086 if topomap_coords:1087 shape = (len(pairs), 2, -1)1088 coords = _find_topomap_coords(info, picks, layout).reshape(shape).mean(axis=1)1089 return picks, coords1090 else:1091 return picks1092 1093 1094def _merge_ch_data(data, ch_type, names, method="rms", *, modality="opm"):1095 """Merge data from channel pairs.1096 1097 Parameters1098 ----------1099 data : array, shape = (n_channels, ..., n_times)1100 Data for channels, ordered in pairs.1101 ch_type : str1102 Channel type.1103 names : list1104 List of channel names.1105 method : str1106 Can be 'rms' or 'mean'.1107 modality : str1108 The modality of the data, either 'grad', 'fnirs', or 'opm'1109 1110 Returns1111 -------1112 data : array, shape = (n_channels / 2, ..., n_times)1113 The root mean square or mean for each pair.1114 names : list1115 List of channel names.1116 """1117 if ch_type == "grad":1118 data = _merge_grad_data(data, method)1119 elif modality == "fnirs" or ch_type in _FNIRS_CH_TYPES_SPLIT:1120 data, names = _merge_nirs_data(data, names)1121 elif modality == "opm" and ch_type == "mag":1122 data, names = _merge_opm_data(data, names)1123 else:1124 raise ValueError(f"Unknown modality {modality} for channel type {ch_type}")1125 1126 return data, names1127 1128 1129def _merge_grad_data(data, method="rms"):1130 """Merge data from channel pairs using the RMS or mean.1131 1132 Parameters1133 ----------1134 data : array, shape = (n_channels, ..., n_times)1135 Data for channels, ordered in pairs.1136 method : str1137 Can be 'rms' or 'mean'.1138 1139 Returns1140 -------1141 data : array, shape = (n_channels / 2, ..., n_times)1142 The root mean square or mean for each pair.1143 """1144 data, orig_shape = data.reshape((len(data) // 2, 2, -1)), data.shape1145 if method == "mean":1146 data = np.mean(data, axis=1)1147 elif method == "rms":1148 data = np.sqrt(np.sum(data**2, axis=1) / 2)1149 else:1150 raise ValueError(f'method must be "rms" or "mean", got {method}.')1151 return data.reshape(data.shape[:1] + orig_shape[1:])1152 1153 1154def _merge_nirs_data(data, merged_names):1155 """Merge data from multiple nirs channel using the mean.1156 1157 Channel names that have an x in them will be merged. The first channel in1158 the name is replaced with the mean of all listed channels. The other1159 channels are removed.1160 1161 Parameters1162 ----------1163 data : array, shape = (n_channels, ..., n_times)1164 Data for channels.1165 merged_names : list1166 List of strings containing the channel names. Channels that are to be1167 merged contain an x between them.1168 1169 Returns1170 -------1171 data : array1172 Data for channels with requested channels merged. Channels used in the1173 merge are removed from the array.1174 """1175 to_remove = np.empty(0, dtype=np.int32)1176 for idx, ch in enumerate(merged_names):1177 if "x" in ch:1178 indices = np.empty(0, dtype=np.int32)1179 channels = ch.split("x")1180 for sub_ch in channels[1:]:1181 indices = np.append(indices, merged_names.index(sub_ch))1182 data[idx] = np.mean(data[np.append(idx, indices)], axis=0)1183 to_remove = np.append(to_remove, indices)1184 to_remove = np.unique(to_remove)1185 for rem in sorted(to_remove, reverse=True):1186 del merged_names[rem]1187 data = np.delete(data, rem, 0)1188 return data, merged_names1189 1190 1191def _merge_opm_data(data, merged_names):1192 """Merge data from multiple opm channel by just using the radial component.1193 1194 Channel names that end in "MERGE_REMOVE" (ie non-radial channels) will be1195 removed. Only the the radial channel is kept.1196 1197 Parameters1198 ----------1199 data : array, shape = (n_channels, ..., n_times)1200 Data for channels.