Aluode/PerceptionLabPortable
0
1"""Freesurfer handling functions."""2 3# Authors: The MNE-Python contributors.4# License: BSD-3-Clause5# Copyright the MNE-Python contributors.6 7import os.path as op8from gzip import GzipFile9from pathlib import Path10 11import numpy as np12 13from ._fiff.constants import FIFF14from ._fiff.meas_info import read_fiducials15from .surface import _read_mri_surface, read_surface16from .transforms import (17 Transform,18 _ensure_trans,19 apply_trans,20 combine_transforms,21 invert_transform,22 read_ras_mni_t,23)24from .utils import (25 _check_fname,26 _check_option,27 _import_nibabel,28 _validate_type,29 get_subjects_dir,30 logger,31 verbose,32)33 34 35def _check_subject_dir(subject, subjects_dir):36 """Check that the Freesurfer subject directory is as expected."""37 subjects_dir = Path(get_subjects_dir(subjects_dir, raise_error=True))38 for img_name in ("T1", "brain", "aseg"):39 if not (subjects_dir / subject / "mri" / f"{img_name}.mgz").is_file():40 raise ValueError(41 "Freesurfer recon-all subject folder "42 "is incorrect or improperly formatted, "43 f"got {subjects_dir / subject}"44 )45 return subjects_dir / subject46 47 48def _get_aseg(aseg, subject, subjects_dir):49 """Check that the anatomical segmentation file exists and load it."""50 nib = _import_nibabel("load aseg")51 subjects_dir = Path(get_subjects_dir(subjects_dir, raise_error=True))52 if aseg == "auto": # use aparc+aseg if auto53 aseg = _check_fname(54 subjects_dir / subject / "mri" / "aparc+aseg.mgz",55 overwrite="read",56 must_exist=False,57 )58 if not aseg: # if doesn't exist use wmparc59 aseg = subjects_dir / subject / "mri" / "wmparc.mgz"60 else:61 aseg = subjects_dir / subject / "mri" / f"{aseg}.mgz"62 _check_fname(aseg, overwrite="read", must_exist=True)63 aseg = nib.load(aseg)64 aseg_data = np.array(aseg.dataobj)65 return aseg, aseg_data66 67 68def _reorient_image(img, axcodes="RAS"):69 """Reorient an image to a given orientation.70 71 Parameters72 ----------73 img : instance of SpatialImage74 The MRI image.75 axcodes : tuple | str76 The axis codes specifying the orientation, e.g. "RAS".77 See :func:`nibabel.orientations.aff2axcodes`.78 79 Returns80 -------81 img_data : ndarray82 The reoriented image data.83 vox_ras_t : ndarray84 The new transform from the new voxels to surface RAS.85 86 Notes87 -----88 .. versionadded:: 0.2489 """90 nib = _import_nibabel("reorient MRI image")91 orig_data = np.array(img.dataobj).astype(np.float32)92 # reorient data to RAS93 ornt = nib.orientations.axcodes2ornt(94 nib.orientations.aff2axcodes(img.affine)95 ).astype(int)96 ras_ornt = nib.orientations.axcodes2ornt(axcodes)97 ornt_trans = nib.orientations.ornt_transform(ornt, ras_ornt)98 img_data = nib.orientations.apply_orientation(orig_data, ornt_trans)99 orig_mgh = nib.MGHImage(orig_data, img.affine)100 aff_trans = nib.orientations.inv_ornt_aff(ornt_trans, img.shape)101 vox_ras_t = np.dot(orig_mgh.header.get_vox2ras_tkr(), aff_trans)102 return img_data, vox_ras_t103 104 105def _mri_orientation(orientation):106 """Get MRI orientation information from an image.107 108 Parameters109 ----------110 orientation : str111 Orientation that you want. Can be "axial", "sagittal", or "coronal".112 113 Returns114 -------115 axis : int116 The dimension of the axis to take slices over when plotting.117 x : int118 The dimension of the x axis.119 y : int120 The dimension of the y axis.121 122 Notes123 -----124 .. versionadded:: 0.21125 .. versionchanged:: 0.24126 """127 _check_option("orientation", orientation, ("coronal", "axial", "sagittal"))128 axis = dict(coronal=1, axial=2, sagittal=0)[orientation]129 x, y = sorted(set([0, 1, 2]).difference(set([axis])))130 return axis, x, y131 132 133def _get_mri_info_data(mri, data):134 # Read the segmentation data using nibabel135 if data:136 _import_nibabel("load MRI atlas data")137 out = dict()138 _, out["vox_mri_t"], out["mri_ras_t"], dims, _, mgz = _read_mri_info(139 mri, return_img=True140 )141 out.update(142 mri_width=dims[0], mri_height=dims[1], mri_depth=dims[1], mri_volume_name=mri143 )144 if data:145 assert mgz is not None146 out["mri_vox_t"] = invert_transform(out["vox_mri_t"])147 out["data"] = np.asarray(mgz.dataobj)148 return out149 150 151def _get_mgz_header(fname):152 """Adapted from nibabel to quickly extract header info."""153 fname = _check_fname(fname, overwrite="read", must_exist=True, name="MRI image")154 if fname.suffix != ".mgz":155 raise OSError("Filename must end with .mgz")156 header_dtd = [157 ("version", ">i4"),158 ("dims", ">i4", (4,)),159 ("type", ">i4"),160 ("dof", ">i4"),161 ("goodRASFlag", ">i2"),162 ("delta", ">f4", (3,)),163 ("Mdc", ">f4", (3, 3)),164 ("Pxyz_c", ">f4", (3,)),165 ]166 header_dtype = np.dtype(header_dtd)167 with GzipFile(fname, "rb") as fid:168 hdr_str = fid.read(header_dtype.itemsize)169 header = np.ndarray(shape=(), dtype=header_dtype, buffer=hdr_str)170 # dims171 dims = header["dims"].astype(int)172 dims = dims[:3] if len(dims) == 4 else dims173 # vox2ras_tkr174 delta = header["delta"]175 ds = np.array(delta, float)176 ns = np.array(dims * ds) / 2.0177 v2rtkr = np.array(178 [179 [-ds[0], 0, 0, ns[0]],180 [0, 0, ds[2], -ns[2]],181 [0, -ds[1], 0, ns[1]],182 [0, 0, 0, 1],183 ],184 dtype=np.float32,185 )186 # ras2vox187 d = np.diag(delta)188 pcrs_c = dims / 2.0189 Mdc = header["Mdc"].T190 pxyz_0 = header["Pxyz_c"] - np.dot(Mdc, np.dot(d, pcrs_c))191 M = np.eye(4, 4)192 M[0:3, 0:3] = np.dot(Mdc, d)193 M[0:3, 3] = pxyz_0.T194 header = dict(dims=dims, vox2ras_tkr=v2rtkr, vox2ras=M, zooms=header["delta"])195 return header196 197 198def _get_atlas_values(vol_info, rr):199 # Transform MRI coordinates (where our surfaces live) to voxels200 rr_vox = apply_trans(vol_info["mri_vox_t"], rr)201 good = (202 (rr_vox >= -0.5) & (rr_vox < np.array(vol_info["data"].shape, int) - 0.5)203 ).all(-1)204 idx = np.round(rr_vox[good].T).astype(np.int64)205 values = np.full(rr.shape[0], np.nan)206 values[good] = vol_info["data"][tuple(idx)]207 return values208 209 210def get_volume_labels_from_aseg(mgz_fname, return_colors=False, atlas_ids=None):211 """Return a list of names and colors of segmented volumes.212 213 Parameters214 ----------215 mgz_fname : path-like216 Filename to read. Typically ``aseg.mgz`` or some variant in the217 freesurfer pipeline.218 return_colors : bool219 If True returns also the labels colors.220 atlas_ids : dict | None221 A lookup table providing a mapping from region names (str) to ID values222 (int). Can be None to use the standard Freesurfer LUT.223 224 .. versionadded:: 0.21.0225 226 Returns227 -------228 label_names : list of str229 The names of segmented volumes included in this mgz file.230 label_colors : list of str231 The RGB colors of the labels included in this mgz file.232 233 See Also234 --------235 read_freesurfer_lut236 237 Notes238 -----239 .. versionchanged:: 0.21.0240 The label names are now sorted in the same order as their corresponding241 values in the MRI file.242 243 .. versionadded:: 0.9.0244 """245 nib = _import_nibabel("load MRI atlas data")246 mgz_fname = _check_fname(247 mgz_fname, overwrite="read", must_exist=True, name="mgz_fname"248 )249 atlas = nib.load(mgz_fname)250 data = np.asarray(atlas.dataobj) # don't need float here251 want = np.unique(data)252 if atlas_ids is None:253 atlas_ids, colors = read_freesurfer_lut()254 elif return_colors:255 raise ValueError("return_colors must be False if atlas_ids are provided")256 # restrict to the ones in the MRI, sorted by label name257 keep = np.isin(list(atlas_ids.values()), want)258 keys = sorted(259 (key for ki, key in enumerate(atlas_ids.keys()) if keep[ki]),260 key=lambda x: atlas_ids[x],261 )262 if return_colors:263 colors = [colors[k] for k in keys]264 out = keys, colors265 else:266 out = keys267 return out268 269 270##############################################################################271# Head to MRI volume conversion272 273 274@verbose275def head_to_mri(276 pos,277 subject,278 mri_head_t,279 subjects_dir=None,280 *,281 kind="mri",282 unscale=False,283 verbose=None,284):285 """Convert pos from head coordinate system to MRI ones.286 287 Parameters288 ----------289 pos : array, shape (n_pos, 3)290 The coordinates (in m) in head coordinate system.291 %(subject)s292 mri_head_t : instance of Transform293 MRI<->Head coordinate transformation.294 %(subjects_dir)s295 kind : str296 The MRI coordinate frame kind, can be ``'mri'`` (default) for297 FreeSurfer surface RAS or ``'ras'`` (default in 1.2) to use MRI RAS298 (scanner RAS).299 300 .. versionadded:: 1.2301 unscale : bool302 For surrogate MRIs (e.g., scaled using ``mne coreg``), if True303 (default False), use the MRI scaling parameters to obtain points in304 the original/surrogate subject's MRI space.305 306 .. versionadded:: 1.2307 %(verbose)s308 309 Returns310 -------311 coordinates : array, shape (n_pos, 3)312 The MRI RAS coordinates (in mm) of pos.313 314 Notes315 -----316 This function requires nibabel.317 """318 from .coreg import read_mri_cfg319 320 _validate_type(kind, str, "kind")321 _check_option("kind", kind, ("ras", "mri"))322 subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)323 t1_fname = subjects_dir / subject / "mri" / "T1.mgz"324 head_mri_t = _ensure_trans(mri_head_t, "head", "mri")325 if kind == "ras":326 _, _, mri_ras_t, _, _ = _read_mri_info(t1_fname)327 head_ras_t = combine_transforms(head_mri_t, mri_ras_t, "head", "ras")328 head_dest_t = head_ras_t329 else:330 assert kind == "mri"331 head_dest_t = head_mri_t332 pos_dest = apply_trans(head_dest_t, pos)333 # unscale if requested334 if unscale:335 params = read_mri_cfg(subject, subjects_dir)336 pos_dest /= params["scale"]337 pos_dest *= 1e3 # mm338 return pos_dest339 340 341##############################################################################342# Surface to MNI conversion343 344 345@verbose346def vertex_to_mni(vertices, hemis, subject, subjects_dir=None, verbose=None):347 """Convert the array of vertices for a hemisphere to MNI coordinates.348 349 Parameters350 ----------351 vertices : int, or list of int352 Vertex number(s) to convert.353 hemis : int, or list of int354 Hemisphere(s) the vertices belong to.355 %(subject)s356 subjects_dir : str, or None357 Path to ``SUBJECTS_DIR`` if it is not set in the environment.358 %(verbose)s359 360 Returns361 -------362 coordinates : array, shape (n_vertices, 3)363 The MNI coordinates (in mm) of the vertices.364 """365 singleton = False366 if not isinstance(vertices, list) and not isinstance(vertices, np.ndarray):367 singleton = True368 vertices = [vertices]369 370 if not isinstance(hemis, list) and not isinstance(hemis, np.ndarray):371 hemis = [hemis] * len(vertices)372 373 if not len(hemis) == len(vertices):374 raise ValueError("hemi and vertices must match in length")375 376 subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)377 378 surfs = [subjects_dir / subject / "surf" / f"{h}.white" for h in ["lh", "rh"]]379 380 # read surface locations in MRI space381 rr = [read_surface(s)[0] for s in surfs]382 383 # take point locations in MRI space and convert to MNI coordinates384 xfm = read_talxfm(subject, subjects_dir)385 xfm["trans"][:3, 3] *= 1000.0 # m->mm386 data = np.array([rr[h][v, :] for h, v in zip(hemis, vertices)])387 if singleton:388 data = data[0]389 return apply_trans(xfm["trans"], data)390 391 392##############################################################################393# Volume to MNI conversion394 395 396@verbose397def head_to_mni(pos, subject, mri_head_t, subjects_dir=None, verbose=None):398 """Convert pos from head coordinate system to MNI ones.399 400 Parameters401 ----------402 pos : array, shape (n_pos, 3)403 The coordinates (in m) in head coordinate system.404 %(subject)s405 mri_head_t : instance of Transform406 MRI<->Head coordinate transformation.407 %(subjects_dir)s408 %(verbose)s409 410 Returns411 -------412 coordinates : array, shape (n_pos, 3)413 The MNI coordinates (in mm) of pos.414 415 Notes416 -----417 This function requires either nibabel.418 """419 subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)420 421 # before we go from head to MRI (surface RAS)422 head_mni_t = combine_transforms(423 _ensure_trans(mri_head_t, "head", "mri"),424 read_talxfm(subject, subjects_dir),425 "head",426 "mni_tal",427 )428 return apply_trans(head_mni_t, pos) * 1000.0429 430 431@verbose432def get_mni_fiducials(subject, subjects_dir=None, verbose=None):433 """Estimate fiducials for a subject.434 435 Parameters436 ----------437 %(subject)s438 %(subjects_dir)s439 %(verbose)s440 441 Returns442 -------443 fids_mri : list444 List of estimated fiducials (each point in a dict), in the order445 LPA, nasion, RPA.446 447 Notes448 -----449 This takes the ``fsaverage-fiducials.fif`` file included with MNE—which450 contain the LPA, nasion, and RPA for the ``fsaverage`` subject—and451 transforms them to the given FreeSurfer subject's MRI space.452 The MRI of ``fsaverage`` is already in MNI Talairach space, so applying453 the inverse of the given subject's MNI Talairach affine transformation454 (``$SUBJECTS_DIR/$SUBJECT/mri/transforms/talairach.xfm``) is used455 to estimate the subject's fiducial locations.456 457 For more details about the coordinate systems and transformations involved,458 see https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems and459 :ref:`tut-source-alignment`.460 """461 # Eventually we might want to allow using the MNI Talairach with-skull462 # transformation rather than the standard brain-based MNI Talaranch463 # transformation, and/or project the points onto the head surface464 # (if available).465 fname_fids_fs = (466 Path(__file__).parent / "data" / "fsaverage" / "fsaverage-fiducials.fif"467 )468 469 # Read fsaverage fiducials file and subject Talairach.470 fids, coord_frame = read_fiducials(fname_fids_fs)471 assert coord_frame == FIFF.FIFFV_COORD_MRI472 if subject == "fsaverage":473 return fids # special short-circuit for fsaverage474 mni_mri_t = invert_transform(read_talxfm(subject, subjects_dir))475 for f in fids:476 f["r"] = apply_trans(mni_mri_t, f["r"])477 return fids478 479 480@verbose481def estimate_head_mri_t(subject, subjects_dir=None, verbose=None):482 """Estimate the head->mri transform from fsaverage fiducials.483 484 A subject's fiducials can be estimated given a Freesurfer ``recon-all``485 by transforming ``fsaverage`` fiducials using the inverse Talairach486 transform, see :func:`mne.coreg.get_mni_fiducials`.487 488 Parameters489 ----------490 %(subject)s491 %(subjects_dir)s492 %(verbose)s493 494 Returns495 -------496 %(trans_not_none)s497 """498 from .channels.montage import compute_native_head_t, make_dig_montage499 500 subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)501 lpa, nasion, rpa = get_mni_fiducials(subject, subjects_dir)502 montage = make_dig_montage(503 lpa=lpa["r"], nasion=nasion["r"], rpa=rpa["r"], coord_frame="mri"504 )505 return invert_transform(compute_native_head_t(montage))506 507 508def _get_affine_from_lta_info(lines):509 """Get the vox2ras affine from lta file info."""510 volume_data = np.loadtxt([line.split("=")[1] for line in lines])511 # get the size of the volume (number of voxels), slice resolution.512 # the matrix of directional cosines and the ras at the center of the bore513 dims, deltas, dir_cos, center_ras = (514 volume_data[0],515 volume_data[1],516 volume_data[2:5],517 volume_data[5],518 )519 dir_cos_delta = dir_cos.T * deltas520 vol_center = (dir_cos_delta @ dims[:3]) / 2521 affine = np.eye(4)522 affine[:3, :3] = dir_cos_delta523 affine[:3, 3] = center_ras - vol_center524 return affine525 526 527@verbose528def read_lta(fname, verbose=None):529 """Read a Freesurfer linear transform array file.530 531 Parameters532 ----------533 fname : path-like534 The transform filename.535 %(verbose)s536 537 Returns538 -------539 affine : ndarray540 The affine transformation described by the lta file.541 """542 _check_fname(fname, "read", must_exist=True)543 with open(fname) as fid:544 lines = fid.readlines()545 # 0 is linear vox2vox, 1 is linear ras2ras546 trans_type = int(lines[0].split("=")[1].strip()[0])547 assert trans_type in (0, 1)548 affine = np.loadtxt(lines[5:9])549 if trans_type == 1:550 return affine551 552 src_affine = _get_affine_from_lta_info(lines[12:18])553 dst_affine = _get_affine_from_lta_info(lines[21:27])554 555 # don't compute if src and dst are already identical556 if np.allclose(src_affine, dst_affine):557 return affine558 559 ras2ras = src_affine @ np.linalg.inv(affine) @ np.linalg.inv(dst_affine)560 affine = np.linalg.inv(np.linalg.inv(src_affine) @ ras2ras @ src_affine)561 return affine562 563 564@verbose565def read_talxfm(subject, subjects_dir=None, verbose=None):566 """Compute MRI-to-MNI transform from FreeSurfer talairach.xfm file.567 568 Parameters569 ----------570 %(subject)s571 %(subjects_dir)s572 %(verbose)s573 574 Returns575 -------576 mri_mni_t : instance of Transform577 The affine transformation from MRI to MNI space for the subject.578 """579 # Adapted from freesurfer m-files. Altered to deal with Norig580 # and Torig correctly581 subjects_dir = get_subjects_dir(subjects_dir)582 # Setup the RAS to MNI transform583 ras_mni_t = read_ras_mni_t(subject, subjects_dir)584 ras_mni_t["trans"][:3, 3] /= 1000.0 # mm->m585 586 # We want to get from Freesurfer surface RAS ('mri') to MNI ('mni_tal').587 # This file only gives us RAS (non-zero origin) ('ras') to MNI ('mni_tal').588 # Se we need to get the ras->mri transform from the MRI headers.589 590 # To do this, we get Norig and Torig591 # (i.e. vox_ras_t and vox_mri_t, respectively)592 path = subjects_dir / subject / "mri" / "orig.mgz"593 if not path.is_file():594 path = subjects_dir / subject / "mri" / "T1.mgz"595 if not path.is_file():596 raise OSError(f"mri not found: {path}")597 _, _, mri_ras_t, _, _ = _read_mri_info(path)598 mri_mni_t = combine_transforms(mri_ras_t, ras_mni_t, "mri", "mni_tal")599 return mri_mni_t600 601 602def _check_mri(mri, subject, subjects_dir) -> str:603 """Check whether an mri exists in the Freesurfer subject directory."""604 _validate_type(mri, "path-like", mri)605 mri = Path(mri)606 if mri.is_file() and mri.name != mri:607 return str(mri)608 elif not mri.is_file():609 if subject is None:610 raise FileNotFoundError(611 f"MRI file {mri!r} not found and no subject provided."612 )613 subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)614 mri = subjects_dir / subject / "mri" / mri615 if not mri.is_file():616 raise FileNotFoundError(617 f"MRI file {mri!r} not found in the subjects directory "618 f"{subjects_dir!r} for subject {subject}."619 )620 if mri.name == mri:621 raise OSError(622 f"Ambiguous filename - found {mri!r} in current folder. "623 "If this is correct prefix name with relative or absolute path."624 )625 return str(mri)626 627 628def _read_mri_info(path, units="m", return_img=False, use_nibabel=False):629 # This is equivalent but 100x slower, so only use nibabel if we need to630 # (later):631 if use_nibabel:632 nib = _import_nibabel()633 hdr = nib.load(path).header634 n_orig = hdr.get_vox2ras()635 t_orig = hdr.get_vox2ras_tkr()636 dims = hdr.get_data_shape()637 zooms = hdr.get_zooms()[:3]638 else:639 hdr = _get_mgz_header(path)640 n_orig = hdr["vox2ras"]641 t_orig = hdr["vox2ras_tkr"]642 dims = hdr["dims"]643 zooms = hdr["zooms"]644 645 # extract the MRI_VOXEL to RAS (non-zero origin) transform646 vox_ras_t = Transform("mri_voxel", "ras", n_orig)647 648 # extract the MRI_VOXEL to MRI transform649 vox_mri_t = Transform("mri_voxel", "mri", t_orig)650 651 # construct the MRI to RAS (non-zero origin) transform652 mri_ras_t = combine_transforms(invert_transform(vox_mri_t), vox_ras_t, "mri", "ras")653 654 assert units in ("m", "mm")655 if units == "m":656 conv = np.array([[1e-3, 1e-3, 1e-3, 1]]).T657 # scaling and translation terms658 vox_ras_t["trans"] *= conv659 vox_mri_t["trans"] *= conv660 # just the translation term661 mri_ras_t["trans"][:, 3:4] *= conv662 663 out = (vox_ras_t, vox_mri_t, mri_ras_t, dims, zooms)664 if return_img:665 nibabel = _import_nibabel()666 out += (nibabel.load(path),)667 return out668 669 670def read_freesurfer_lut(fname=None):671 """Read a Freesurfer-formatted LUT.672 673 Parameters674 ----------675 fname : path-like | None676 The filename. Can be None to read the standard Freesurfer LUT.677 678 Returns679 -------680 atlas_ids : dict681 Mapping from label names to IDs.682 colors : dict683 Mapping from label names to colors.684 """685 lut = _get_lut(fname)686 names, ids = lut["name"], lut["id"]687 colors = np.array([lut["R"], lut["G"], lut["B"], lut["A"]], float).T688 atlas_ids = dict(zip(names, ids))689 colors = dict(zip(names, colors))690 return atlas_ids, colors691 692 693def _get_lut(fname=None):694 """Get a FreeSurfer LUT."""695 if fname is None:696 fname = Path(__file__).parent / "data" / "FreeSurferColorLUT.txt"697 _check_fname(fname, "read", must_exist=True)698 dtype = [699 ("id", "<i8"),700 ("name", "U"),701 ("R", "<i8"),702 ("G", "<i8"),703 ("B", "<i8"),704 ("A", "<i8"),705 ]706 lut = {d[0]: list() for d in dtype}707 with open(fname) as fid:708 for line in fid:709 line = line.strip()710 if line.startswith("#") or not line:711 continue712 line = line.split()713 if len(line) != len(dtype):714 raise RuntimeError(f"LUT is improperly formatted: {fname}")715 for d, part in zip(dtype, line):716 lut[d[0]].append(part)717 lut = {d[0]: np.array(lut[d[0]], dtype=d[1]) for d in dtype}718 assert len(lut["name"]) > 0719 lut["name"] = [str(name) for name in lut["name"]]720 return lut721 722 723@verbose724def _get_head_surface(surf, subject, subjects_dir, bem=None, verbose=None):725 """Get a head surface from the Freesurfer subject directory.726 727 Parameters728 ----------729 surf : str730 The name of the surface 'auto', 'head', 'outer_skin', 'head-dense'731 or 'seghead'.732 %(subject)s733 %(subjects_dir)s734 bem : mne.bem.ConductorModel | None735 The conductor model that stores information about the head surface.736 %(verbose)s737 738 Returns739 -------740 head_surf : dict | None741 A dictionary with keys 'rr', 'tris', 'ntri', 'use_tris', 'np'742 and 'coord_frame' that store information for mesh plotting and other743 useful information about the head surface.744 745 Notes746 -----747 .. versionadded: 0.24748 """749 from .bem import _bem_find_surface, read_bem_surfaces750 751 _check_option("surf", surf, ("auto", "head", "outer_skin", "head-dense", "seghead"))752 if surf in ("auto", "head", "outer_skin"):753 if bem is not None:754 try:755 return _bem_find_surface(bem, "head")756 except RuntimeError:757 logger.info(758 "Could not find the surface for "759 "head in the provided BEM model, "760 "looking in the subject directory."761 )762 if subject is None:763 if surf == "auto":764 return765 raise ValueError(766 "To plot the head surface, the BEM/sphere"767 " model must contain a head surface "768 'or "subject" must be provided (got '769 "None)"770 )771 subject_dir = op.join(get_subjects_dir(subjects_dir, raise_error=True), subject)772 if surf in ("head-dense", "seghead"):773 try_fnames = [774 op.join(subject_dir, "bem", f"{subject}-head-dense.fif"),775 op.join(subject_dir, "surf", "lh.seghead"),776 ]777 else:778 try_fnames = [779 op.join(subject_dir, "bem", "outer_skin.surf"),780 op.join(subject_dir, "bem", "flash", "outer_skin.surf"),781 op.join(subject_dir, "bem", f"{subject}-head-sparse.fif"),782 op.join(subject_dir, "bem", f"{subject}-head.fif"),783 ]784 for fname in try_fnames:785 if op.exists(fname):786 logger.info(f"Using {op.basename(fname)} for head surface.")787 if op.splitext(fname)[-1] == ".fif":788 return read_bem_surfaces(fname, on_defects="warn")[0]789 else:790 return _read_mri_surface(fname)791 raise OSError(792 "No head surface found for subject "793 f"{subject} after trying:\n" + "\n".join(try_fnames)794 )795 796 797@verbose798def _get_skull_surface(surf, subject, subjects_dir, bem=None, verbose=None):799 """Get a skull surface from the Freesurfer subject directory.800 801 Parameters802 ----------803 surf : str804 The name of the surface 'outer' or 'inner'.805 %(subject)s806 %(subjects_dir)s807 bem : mne.bem.ConductorModel | None808 The conductor model that stores information about the skull surface.809 %(verbose)s810 811 Returns812 -------813 skull_surf : dict | None814 A dictionary with keys 'rr', 'tris', 'ntri', 'use_tris', 'np'815 and 'coord_frame' that store information for mesh plotting and other816 useful information about the head surface.817 818 Notes819 -----820 .. versionadded: 0.24821 """822 from .bem import _bem_find_surface823 824 if bem is not None:825 try:826 return _bem_find_surface(bem, surf + "_skull")827 except RuntimeError:828 logger.info(829 "Could not find the surface for "830 "skull in the provided BEM model, "831 "looking in the subject directory."832 )833 subjects_dir = Path(get_subjects_dir(subjects_dir, raise_error=True))834 fname = _check_fname(835 subjects_dir / subject / "bem" / (surf + "_skull.surf"),836 overwrite="read",837 must_exist=True,838 name=f"{surf} skull surface",839 )840 return _read_mri_surface(fname)841 842 843def _estimate_talxfm_rigid(subject, subjects_dir):844 from .coreg import _trans_from_params, fit_matched_points845 846 xfm = read_talxfm(subject, subjects_dir)847 # XYZ+origin + halfway848 pts_tal = np.concatenate([np.eye(4)[:, :3], np.eye(3) * 0.5])849 pts_subj = apply_trans(invert_transform(xfm), pts_tal)850 # we fit with scaling enabled, but then discard it (we just need851 # the rigid-body components)852 params = fit_matched_points(pts_subj, pts_tal, scale=3, out="params")853 rigid = _trans_from_params((True, True, False), params[:6])854 return rigid855