Aluode/PerceptionLabPortable

1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5# Many of the computations in this code were derived from Matti Hämäläinen's6# C code.7 8import os9 10import numpy as np11from scipy.sparse import csr_array12 13from ._fiff.constants import FIFF14from ._fiff.open import fiff_open15from ._fiff.tag import find_tag16from ._fiff.tree import dir_tree_find17from ._fiff.write import (18    end_block,19    start_and_end_file,20    start_block,21    write_float_sparse_rcs,22    write_int,23    write_string,24)25from .fixes import _eye_array26from .surface import (27    _compute_nearest,28    _find_nearest_tri_pts,29    _get_tri_supp_geom,30    _normalize_vectors,31    _triangle_neighbors,32    read_surface,33)34from .utils import get_subjects_dir, logger, verbose, warn35 36 37@verbose38def read_morph_map(39    subject_from, subject_to, subjects_dir=None, xhemi=False, verbose=None40):41    """Read morph map.42 43    Morph maps can be generated with mne_make_morph_maps. If one isn't44    available, it will be generated automatically and saved to the45    ``subjects_dir/morph_maps`` directory.46 47    Parameters48    ----------49    subject_from : str50        Name of the original subject as named in the ``SUBJECTS_DIR``.51    subject_to : str52        Name of the subject on which to morph as named in the ``SUBJECTS_DIR``.53    subjects_dir : path-like54        Path to ``SUBJECTS_DIR`` is not set in the environment.55    xhemi : bool56        Morph across hemisphere. Currently only implemented for57        ``subject_to == subject_from``. See notes of58        :func:`mne.compute_source_morph`.59    %(verbose)s60 61    Returns62    -------63    left_map, right_map : ~scipy.sparse.csr_array64        The morph maps for the 2 hemispheres.65    """66    subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)67 68    # First check for morph-map dir existence69    mmap_dir = subjects_dir / "morph-maps"70    if not mmap_dir.is_dir():71        try:72            os.mkdir(mmap_dir)73        except Exception:74            warn(f'Could not find or make morph map directory "{mmap_dir}"')75 76    # filename components77    if xhemi:78        if subject_to != subject_from:79            raise NotImplementedError(80                "Morph-maps between hemispheres are currently only "81                "implemented for subject_to == subject_from"82            )83        map_name_temp = "%s-%s-xhemi"84        log_msg = "Creating morph map %s -> %s xhemi"85    else:86        map_name_temp = "%s-%s"87        log_msg = "Creating morph map %s -> %s"88 89    map_names = [90        map_name_temp % (subject_from, subject_to),91        map_name_temp % (subject_to, subject_from),92    ]93 94    # find existing file95    fname = None96    for map_name in map_names:97        fname = mmap_dir / f"{map_name}-morph.fif"98        if fname.exists():99            return _read_morph_map(fname, subject_from, subject_to)100    # if file does not exist, make it101    logger.info(102        f'Morph map "{fname}" does not exist, creating it and saving it to disk'103    )104    logger.info(log_msg % (subject_from, subject_to))105    mmap_1 = _make_morph_map(subject_from, subject_to, subjects_dir, xhemi)106    if subject_to == subject_from:107        mmap_2 = None108    else:109        logger.info(log_msg % (subject_to, subject_from))110        mmap_2 = _make_morph_map(subject_to, subject_from, subjects_dir, xhemi)111    _write_morph_map(fname, subject_from, subject_to, mmap_1, mmap_2)112    return mmap_1113 114 115def _read_morph_map(fname, subject_from, subject_to):116    """Read a morph map from disk."""117    f, tree, _ = fiff_open(fname)118    with f as fid:119        # Locate all maps120        maps = dir_tree_find(tree, FIFF.FIFFB_MNE_MORPH_MAP)121        if len(maps) == 0:122            raise ValueError("Morphing map data not found")123 124        # Find the correct ones125        left_map = None126        right_map = None127        for m in maps:128            tag = find_tag(fid, m, FIFF.FIFF_MNE_MORPH_MAP_FROM)129            if tag.data == subject_from:130                tag = find_tag(fid, m, FIFF.FIFF_MNE_MORPH_MAP_TO)131                if tag.data == subject_to:132                    #  Names match: which hemishere is this?133                    tag = find_tag(fid, m, FIFF.FIFF_MNE_HEMI)134                    if tag.data == FIFF.FIFFV_MNE_SURF_LEFT_HEMI:135                        tag = find_tag(fid, m, FIFF.FIFF_MNE_MORPH_MAP)136                        left_map = tag.data137                        logger.info("    Left-hemisphere map read.")138                    elif tag.data == FIFF.FIFFV_MNE_SURF_RIGHT_HEMI:139                        tag = find_tag(fid, m, FIFF.FIFF_MNE_MORPH_MAP)140                        right_map = tag.data141                        logger.info("    Right-hemisphere map read.")142 143    if left_map is None or right_map is None:144        raise ValueError(f"Could not find both hemispheres in {fname}")145 146    return left_map, right_map147 148 149def _write_morph_map(fname, subject_from, subject_to, mmap_1, mmap_2):150    """Write a morph map to disk."""151    try:152        with start_and_end_file(fname) as fid:153            _write_morph_map_(fid, subject_from, subject_to, mmap_1, mmap_2)154    except Exception as exp:155        warn(f'Could not write morph-map file "{fname}" (error: {exp})')156 157 158def _write_morph_map_(fid, subject_from, subject_to, mmap_1, mmap_2):159    assert len(mmap_1) == 2160    hemis = [FIFF.FIFFV_MNE_SURF_LEFT_HEMI, FIFF.FIFFV_MNE_SURF_RIGHT_HEMI]161    for m, hemi in zip(mmap_1, hemis):162        start_block(fid, FIFF.FIFFB_MNE_MORPH_MAP)163        write_string(fid, FIFF.FIFF_MNE_MORPH_MAP_FROM, subject_from)164        write_string(fid, FIFF.FIFF_MNE_MORPH_MAP_TO, subject_to)165        write_int(fid, FIFF.FIFF_MNE_HEMI, hemi)166        write_float_sparse_rcs(fid, FIFF.FIFF_MNE_MORPH_MAP, m)167        end_block(fid, FIFF.FIFFB_MNE_MORPH_MAP)168    # don't write mmap_2 if it is identical (subject_to == subject_from)169    if mmap_2 is not None:170        assert len(mmap_2) == 2171        for m, hemi in zip(mmap_2, hemis):172            start_block(fid, FIFF.FIFFB_MNE_MORPH_MAP)173            write_string(fid, FIFF.FIFF_MNE_MORPH_MAP_FROM, subject_to)174            write_string(fid, FIFF.FIFF_MNE_MORPH_MAP_TO, subject_from)175            write_int(fid, FIFF.FIFF_MNE_HEMI, hemi)176            write_float_sparse_rcs(fid, FIFF.FIFF_MNE_MORPH_MAP, m)177            end_block(fid, FIFF.FIFFB_MNE_MORPH_MAP)178 179 180def _make_morph_map(subject_from, subject_to, subjects_dir, xhemi):181    """Construct morph map from one subject to another.182 183    Note that this is close, but not exactly like the C version.184    For example, parts are more accurate due to double precision,185    so expect some small morph-map differences!186 187    Note: This seems easily parallelizable, but the overhead188    of pickling all the data structures makes it less efficient189    than just running on a single core :(190    """191    subjects_dir = get_subjects_dir(subjects_dir)192    if xhemi:193        reg = "%s.sphere.left_right"194        hemis = (("lh", "rh"), ("rh", "lh"))195    else:196        reg = "%s.sphere.reg"197        hemis = (("lh", "lh"), ("rh", "rh"))198 199    return [200        _make_morph_map_hemi(201            subject_from, subject_to, subjects_dir, reg % hemi_from, reg % hemi_to202        )203        for hemi_from, hemi_to in hemis204    ]205 206 207def _make_morph_map_hemi(subject_from, subject_to, subjects_dir, reg_from, reg_to):208    """Construct morph map for one hemisphere."""209    # add speedy short-circuit for self-maps210    if subject_from == subject_to and reg_from == reg_to:211        fname = subjects_dir / subject_from / "surf" / reg_from212        n_pts = len(read_surface(fname, verbose=False)[0])213        return _eye_array(n_pts, format="csr")214 215    # load surfaces and normalize points to be on unit sphere216    fname = subjects_dir / subject_from / "surf" / reg_from217    from_rr, from_tri = read_surface(fname, verbose=False)218    fname = subjects_dir / subject_to / "surf" / reg_to219    to_rr = read_surface(fname, verbose=False)[0]220    _normalize_vectors(from_rr)221    _normalize_vectors(to_rr)222 223    # from surface: get nearest neighbors, find triangles for each vertex224    nn_pts_idx = _compute_nearest(from_rr, to_rr, method="KDTree")225    from_pt_tris = _triangle_neighbors(from_tri, len(from_rr))226    from_pt_tris = [from_pt_tris[pt_idx].astype(int) for pt_idx in nn_pts_idx]227    from_pt_lens = np.cumsum([0] + [len(x) for x in from_pt_tris])228    from_pt_tris = np.concatenate(from_pt_tris)229    assert from_pt_tris.ndim == 1230    assert from_pt_lens[-1] == len(from_pt_tris)231 232    # find triangle in which point lies and assoc. weights233    tri_inds = []234    weights = []235    tri_geom = _get_tri_supp_geom(dict(rr=from_rr, tris=from_tri))236    weights, tri_inds = _find_nearest_tri_pts(237        to_rr, from_pt_tris, from_pt_lens, run_all=False, reproject=False, **tri_geom238    )239 240    nn_idx = from_tri[tri_inds]241    weights = np.array(weights)242 243    row_ind = np.repeat(np.arange(len(to_rr)), 3)244    this_map = csr_array(245        (weights.ravel(), (row_ind, nn_idx.ravel())), shape=(len(to_rr), len(from_rr))246    )247    return this_map248