Aluode/PerceptionLabPortable

1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5# The computations in this code were primarily derived from Matti Hämäläinen's6# C code.7 8import inspect9from copy import deepcopy10 11import numpy as np12from scipy.interpolate import interp1d13 14from .._fiff.constants import FIFF15from .._fiff.meas_info import _simplify_info16from .._fiff.pick import pick_info, pick_types17from .._fiff.proj import _has_eeg_average_ref_proj, make_projector18from ..bem import _check_origin19from ..cov import make_ad_hoc_cov20from ..epochs import BaseEpochs, EpochsArray21from ..evoked import Evoked, EvokedArray22from ..fixes import _safe_svd23from ..surface import get_head_surf, get_meg_helmet_surf24from ..transforms import _find_trans, transform_surface_to25from ..utils import _check_fname, _check_option, _pl, _reg_pinv, logger, verbose, warn26from ._lead_dots import (27    _do_cross_dots,28    _do_self_dots,29    _do_surface_dots,30    _get_legen_table,31)32from ._make_forward import _create_eeg_els, _create_meg_coils, _read_coil_defs33 34 35def _setup_dots(mode, info, coils, ch_type):36    """Set up dot products."""37    int_rad = 0.0638    noise = make_ad_hoc_cov(info, dict(mag=20e-15, grad=5e-13, eeg=1e-6))39    n_coeff, interp = (50, "nearest") if mode == "fast" else (100, "linear")40    lut, n_fact = _get_legen_table(ch_type, False, n_coeff, verbose=False)41    lut_fun = interp1d(np.linspace(-1, 1, lut.shape[0]), lut, interp, axis=0)42    return int_rad, noise, lut_fun, n_fact43 44 45def _compute_mapping_matrix(fmd, info):46    """Do the hairy computations."""47    logger.info("    Preparing the mapping matrix...")48    # assemble a projector and apply it to the data49    ch_names = fmd["ch_names"]50    projs = info.get("projs", list())51    proj_op = make_projector(projs, ch_names)[0]52    proj_dots = np.dot(proj_op.T, np.dot(fmd["self_dots"], proj_op))53 54    noise_cov = fmd["noise"]55    # Whiten56    if not noise_cov["diag"]:57        raise NotImplementedError  # this shouldn't happen58    whitener = np.diag(1.0 / np.sqrt(noise_cov["data"].ravel()))59    whitened_dots = np.dot(whitener.T, np.dot(proj_dots, whitener))60 61    # SVD is numerically better than the eigenvalue composition even if62    # mat is supposed to be symmetric and positive definite63    if fmd.get("pinv_method", "tsvd") == "tsvd":64        inv, fmd["nest"] = _pinv_trunc(whitened_dots, fmd["miss"])65    else:66        assert fmd["pinv_method"] == "tikhonov", fmd["pinv_method"]67        inv, fmd["nest"] = _pinv_tikhonov(whitened_dots, fmd["miss"])68 69    # Sandwich with the whitener70    inv_whitened = np.dot(whitener.T, np.dot(inv, whitener))71 72    # Take into account that the lead fields used to compute73    # d->surface_dots were unprojected74    inv_whitened_proj = proj_op.T @ inv_whitened75 76    # Finally sandwich in the selection matrix77    # This one picks up the correct lead field projection78    mapping_mat = np.dot(fmd["surface_dots"], inv_whitened_proj)79 80    # Optionally apply the average electrode reference to the final field map81    if fmd["kind"] == "eeg" and _has_eeg_average_ref_proj(info):82        logger.info(83            "    The map has an average electrode reference "84            f"({mapping_mat.shape[0]} channels)"85        )86        mapping_mat -= np.mean(mapping_mat, axis=0)87    return mapping_mat88 89 90def _pinv_trunc(x, miss):91    """Compute pseudoinverse, truncating at most "miss" fraction of varexp."""92    u, s, v = _safe_svd(x, full_matrices=False)93 94    # Eigenvalue truncation95    varexp = np.cumsum(s)96    varexp /= varexp[-1]97    n = np.where(varexp >= (1.0 - miss))[0][0] + 198    logger.info(99        "    Truncating at %d/%d components to omit less than %g (%0.2g)",100        n,101        len(s),102        miss,103        1.0 - varexp[n - 1],104    )105    s = 1.0 / s[:n]106    inv = ((u[:, :n] * s) @ v[:n]).T107    return inv, n108 109 110def _pinv_tikhonov(x, reg):111    # _reg_pinv requires square Hermitian, which we have here112    inv, _, n = _reg_pinv(x, reg=reg, rank=None)113    logger.info(114        f"    Truncating at {n}/{len(x)} components and regularizing with α={reg:0.1e}"115    )116    return inv, n117 118 119def _map_meg_or_eeg_channels(info_from, info_to, mode, *, origin, miss=None):120    """Find mapping from one set of channels to another.121 122    Parameters123    ----------124    info_from : instance of Info125        The measurement data to interpolate from.126    info_to : instance of Info127        The measurement info to interpolate to.128    mode : str129        Either `'accurate'` or `'fast'`, determines the quality of the130        Legendre polynomial expansion used. `'fast'` should be sufficient131        for most applications.132    origin : array-like, shape (3,) | str133        Origin of the sphere in the head coordinate frame and in meters.134        Can be ``'auto'``, which means a head-digitization-based origin135        fit.136 137    Returns138    -------139    mapping : array, shape (n_to, n_from)140        A mapping matrix.141    """142    assert origin is not None  # should be assured elsewhere143 144    # no need to apply trans because both from and to coils are in device145    # coordinates146    info_kinds = set(ch["kind"] for ch in info_to["chs"])147    info_kinds |= set(ch["kind"] for ch in info_from["chs"])148    if FIFF.FIFFV_REF_MEG_CH in info_kinds:  # refs same as MEG149        info_kinds |= set([FIFF.FIFFV_MEG_CH])150        info_kinds -= set([FIFF.FIFFV_REF_MEG_CH])151    info_kinds = sorted(info_kinds)152    # This should be guaranteed by the callers153    assert len(info_kinds) == 1 and info_kinds[0] in (154        FIFF.FIFFV_MEG_CH,155        FIFF.FIFFV_EEG_CH,156    )157    kind = "eeg" if info_kinds[0] == FIFF.FIFFV_EEG_CH else "meg"158 159    #160    # Step 1. Prepare the coil definitions161    #162    if kind == "meg":163        templates = _read_coil_defs(verbose=False)164        coils_from = _create_meg_coils(165            info_from["chs"], "normal", info_from["dev_head_t"], templates166        )167        coils_to = _create_meg_coils(168            info_to["chs"], "normal", info_to["dev_head_t"], templates169        )170        pinv_method = "tsvd"171        miss = 1e-4172    else:173        coils_from = _create_eeg_els(info_from["chs"])174        coils_to = _create_eeg_els(info_to["chs"])175        pinv_method = "tikhonov"176        miss = 1e-1177        if _has_eeg_average_ref_proj(info_from) and not _has_eeg_average_ref_proj(178            info_to179        ):180            raise RuntimeError(181                "info_to must have an average EEG reference projector if "182                "info_from has one"183            )184    origin = _check_origin(origin, info_from)185    #186    # Step 2. Calculate the dot products187    #188    int_rad, noise, lut_fun, n_fact = _setup_dots(mode, info_from, coils_from, kind)189    logger.info(190        f"    Computing dot products for {len(coils_from)} "191        f"{kind.upper()} channel{_pl(coils_from)}..."192    )193    self_dots = _do_self_dots(194        int_rad, False, coils_from, origin, kind, lut_fun, n_fact, n_jobs=None195    )196    logger.info(197        f"    Computing cross products for {len(coils_from)} → "198        f"{len(coils_to)} {kind.upper()} channel{_pl(coils_to)}..."199    )200    cross_dots = _do_cross_dots(201        int_rad, False, coils_from, coils_to, origin, kind, lut_fun, n_fact202    ).T203 204    ch_names = [c["ch_name"] for c in info_from["chs"]]205    fmd = dict(206        kind=kind,207        ch_names=ch_names,208        origin=origin,209        noise=noise,210        self_dots=self_dots,211        surface_dots=cross_dots,212        int_rad=int_rad,213        miss=miss,214        pinv_method=pinv_method,215    )216 217    #218    # Step 3. Compute the mapping matrix219    #220    mapping = _compute_mapping_matrix(fmd, info_from)221    return mapping222 223 224def _as_meg_type_inst(inst, ch_type="grad", mode="fast"):225    """Compute virtual evoked using interpolated fields in mag/grad channels.226 227    Parameters228    ----------229    inst : instance of mne.Evoked or mne.Epochs230        The evoked or epochs object.231    ch_type : str232        The destination channel type. It can be 'mag' or 'grad'.233    mode : str234        Either `'accurate'` or `'fast'`, determines the quality of the235        Legendre polynomial expansion used. `'fast'` should be sufficient236        for most applications.237 238    Returns239    -------240    inst : instance of mne.EvokedArray or mne.EpochsArray241        The transformed evoked object containing only virtual channels.242    """243    _check_option("ch_type", ch_type, ["mag", "grad"])244 245    # pick the original and destination channels246    pick_from = pick_types(inst.info, meg=True, eeg=False, ref_meg=False)247    pick_to = pick_types(inst.info, meg=ch_type, eeg=False, ref_meg=False)248 249    if len(pick_to) == 0:250        raise ValueError(251            "No channels matching the destination channel type"252            " found in info. Please pass an evoked containing"253            "both the original and destination channels. Only the"254            " locations of the destination channels will be used"255            " for interpolation."256        )257 258    info_from = pick_info(inst.info, pick_from)259    info_to = pick_info(inst.info, pick_to)260    # XXX someday we should probably expose the origin261    mapping = _map_meg_or_eeg_channels(262        info_from, info_to, origin=(0.0, 0.0, 0.04), mode=mode263    )264 265    # compute data by multiplying by the 'gain matrix' from266    # original sensors to virtual sensors267    if hasattr(inst, "get_data"):268        kwargs = dict()269        if "copy" in inspect.getfullargspec(inst.get_data).kwonlyargs:270            kwargs["copy"] = False271        data = inst.get_data(**kwargs)272    else:273        data = inst.data274 275    ndim = data.ndim276    if ndim == 2:277        data = data[np.newaxis, :, :]278 279    data_ = np.empty((data.shape[0], len(mapping), data.shape[2]), dtype=data.dtype)280    for d, d_ in zip(data, data_):281        d_[:] = np.dot(mapping, d[pick_from])282 283    # keep only the destination channel types284    info = pick_info(inst.info, sel=pick_to, copy=True)285 286    # change channel names to emphasize they contain interpolated data287    for ch in info["chs"]:288        ch["ch_name"] += "_v"289    info._update_redundant()290    info._check_consistency()291    if isinstance(inst, Evoked):292        assert ndim == 2293        data_ = data_[0]  # undo new axis294        inst_ = EvokedArray(295            data_, info, tmin=inst.times[0], comment=inst.comment, nave=inst.nave296        )297    else:298        assert isinstance(inst, BaseEpochs)299        inst_ = EpochsArray(300            data_,301            info,302            tmin=inst.tmin,303            events=inst.events,304            event_id=inst.event_id,305            metadata=inst.metadata,306        )307 308    return inst_309 310 311@verbose312def _make_surface_mapping(313    info,314    surf,315    ch_type="meg",316    trans=None,317    mode="fast",318    n_jobs=None,319    *,320    origin,321    verbose=None,322):323    """Re-map M/EEG data to a surface.324 325    Parameters326    ----------327    %(info_not_none)s328    surf : dict329        The surface to map the data to. The required fields are `'rr'`,330        `'nn'`, and `'coord_frame'`. Must be in head coordinates.331    ch_type : str332        Must be either `'meg'` or `'eeg'`, determines the type of field.333    trans : None | dict334        If None, no transformation applied. Should be a Head<->MRI335        transformation.336    mode : str337        Either `'accurate'` or `'fast'`, determines the quality of the338        Legendre polynomial expansion used. `'fast'` should be sufficient339        for most applications.340    %(n_jobs)s341    origin : array-like, shape (3,) | str342        Origin of the sphere in the head coordinate frame and in meters.343    %(verbose)s344 345    Returns346    -------347    mapping : array348        A n_vertices x n_sensors array that remaps the MEG or EEG data,349        as `new_data = np.dot(mapping, data)`.350    """351    assert origin is not None  # should be assured elsewhere352 353    if not all(key in surf for key in ["rr", "nn"]):354        raise KeyError('surf must have both "rr" and "nn"')355    if "coord_frame" not in surf:356        raise KeyError(357            'The surface coordinate frame must be specified in surf["coord_frame"]'358        )359    _check_option("mode", mode, ["accurate", "fast"])360 361    # deal with coordinate frames here -- always go to "head" (easiest)362    orig_surf = surf363    surf = transform_surface_to(deepcopy(surf), "head", trans)364    origin = _check_origin(origin, info)365 366    #367    # Step 1. Prepare the coil definitions368    # Do the dot products, assume surf in head coords369    #370    _check_option("ch_type", ch_type, ["meg", "eeg"])371    if ch_type == "meg":372        picks = pick_types(info, meg=True, eeg=False, ref_meg=False)373        logger.info("Prepare MEG mapping...")374    else:375        picks = pick_types(info, meg=False, eeg=True, ref_meg=False)376        logger.info("Prepare EEG mapping...")377    if len(picks) == 0:378        raise RuntimeError("cannot map, no channels found")379    # XXX this code does not do any checking for compensation channels,380    # but it seems like this must be intentional from the ref_meg=False381    # (presumably from the C code)382    dev_head_t = info["dev_head_t"]383    info = pick_info(_simplify_info(info), picks)384    info["dev_head_t"] = dev_head_t385 386    # create coil defs in head coordinates387    if ch_type == "meg":388        # Put them in head coordinates389        coils = _create_meg_coils(info["chs"], "normal", info["dev_head_t"])390        type_str = "coils"391        miss = 1e-4  # Smoothing criterion for MEG392    else:  # EEG393        coils = _create_eeg_els(info["chs"])394        type_str = "electrodes"395        miss = 1e-3  # Smoothing criterion for EEG396 397    #398    # Step 2. Calculate the dot products399    #400    int_rad, noise, lut_fun, n_fact = _setup_dots(mode, info, coils, ch_type)401    logger.info("Computing dot products for %i %s...", len(coils), type_str)402    self_dots = _do_self_dots(403        int_rad, False, coils, origin, ch_type, lut_fun, n_fact, n_jobs404    )405    sel = np.arange(len(surf["rr"]))  # eventually we should do sub-selection406    logger.info("Computing dot products for %i surface locations...", len(sel))407    surface_dots = _do_surface_dots(408        int_rad, False, coils, surf, sel, origin, ch_type, lut_fun, n_fact, n_jobs409    )410 411    #412    # Step 4. Return the result413    #414    fmd = dict(415        kind=ch_type,416        surf=surf,417        ch_names=info["ch_names"],418        coils=coils,419        origin=origin,420        noise=noise,421        self_dots=self_dots,422        surface_dots=surface_dots,423        int_rad=int_rad,424        miss=miss,425    )426    logger.info("Field mapping data ready")427 428    fmd["data"] = _compute_mapping_matrix(fmd, info)429    # bring the original back, whatever coord frame it was in430    fmd["surf"] = orig_surf431 432    # Remove some unnecessary fields433    del fmd["self_dots"]434    del fmd["surface_dots"]435    del fmd["int_rad"]436    del fmd["miss"]437    return fmd438 439 440@verbose441def make_field_map(442    evoked,443    trans="auto",444    subject=None,445    subjects_dir=None,446    ch_type=None,447    mode="fast",448    meg_surf="helmet",449    origin=None,450    n_jobs=None,451    *,452    upsampling=1,453    head_source=("bem", "head"),454    verbose=None,455):456    """Compute surface maps used for field display in 3D.457 458    Parameters459    ----------460    evoked : Evoked | Epochs | Raw461        The measurement file. Need to have info attribute.462    %(trans)s ``"auto"`` (default) will load trans from the FreeSurfer463        directory specified by ``subject`` and ``subjects_dir`` parameters.464 465        .. versionchanged:: 0.19466            Support for ``'fsaverage'`` argument.467    subject : str | None468        The subject name corresponding to FreeSurfer environment469        variable SUBJECT. If None, map for EEG data will not be available.470    subjects_dir : path-like471        The path to the freesurfer subjects reconstructions.472        It corresponds to Freesurfer environment variable SUBJECTS_DIR.473    ch_type : None | ``'eeg'`` | ``'meg'``474        If None, a map for each available channel type will be returned.475        Else only the specified type will be used.476    mode : ``'accurate'`` | ``'fast'``477        Either ``'accurate'`` or ``'fast'``, determines the quality of the478        Legendre polynomial expansion used. ``'fast'`` should be sufficient479        for most applications.480    meg_surf : 'helmet' | 'head'481        Should be ``'helmet'`` or ``'head'`` to specify in which surface482        to compute the MEG field map. The default value is ``'helmet'``.483    origin : array-like, shape (3,) | 'auto'484        Origin of the sphere in the head coordinate frame and in meters.485        Can be ``'auto'``, which means a head-digitization-based origin486        fit. Default is ``(0., 0., 0.04)``.487 488        .. versionadded:: 0.11489        .. versionchanged:: 1.12490           In 1.12 the default value is "auto".491           In 1.11 and prior versions, it is ``(0., 0., 0.04)``.492    %(n_jobs)s493    %(helmet_upsampling)s494 495        .. versionadded:: 1.10496    %(head_source)s497 498        .. versionadded:: 1.1499    %(verbose)s500 501    Returns502    -------503    surf_maps : list504        The surface maps to be used for field plots. The list contains505        separate ones for MEG and EEG (if both MEG and EEG are present).506    """507    if origin is None:508        warn_message = (509            'Default value for origin is "(0.0, 0.0, 0.04)" in version 1.11 '510            'but will be changed to "auto" in 1.12. Set the origin parameter '511            "explicitly to avoid this warning."512        )513        warn(warn_message, FutureWarning)514        origin = (0.0, 0.0, 0.04)515 516    info = evoked.info517 518    if ch_type is None:519        types = [t for t in ["eeg", "meg"] if t in evoked]520    else:521        _check_option("ch_type", ch_type, ["eeg", "meg"])522        types = [ch_type]523 524    if subjects_dir is not None:525        subjects_dir = _check_fname(526            subjects_dir,527            overwrite="read",528            must_exist=True,529            name="subjects_dir",530            need_dir=True,531        )532 533    trans, trans_type = _find_trans(534        trans=trans,535        subject=subject,536        subjects_dir=subjects_dir,537    )538 539    if "eeg" in types and trans_type == "identity":540        logger.info("No trans file available. EEG data ignored.")541        types.remove("eeg")542 543    if len(types) == 0:544        raise RuntimeError("No data available for mapping.")545 546    _check_option("meg_surf", meg_surf, ["helmet", "head"])547 548    surfs = []549    for this_type in types:550        if this_type == "meg" and meg_surf == "helmet":551            surf = get_meg_helmet_surf(info, trans, upsampling=upsampling)552        else:553            surf = get_head_surf(subject, source=head_source, subjects_dir=subjects_dir)554        surfs.append(surf)555 556    surf_maps = list()557 558    for this_type, this_surf in zip(types, surfs):559        this_map = _make_surface_mapping(560            evoked.info,561            this_surf,562            this_type,563            trans,564            n_jobs=n_jobs,565            origin=origin,566            mode=mode,567        )568        surf_maps.append(this_map)569 570    return surf_maps571