import json
import time
from collections import OrderedDict
from copy import deepcopy
from functools import lru_cache, partial
from glob import glob
from os import path as op
from pathlib import Path
import numpy as np
from ._fiff.constants import FIFF
from ._fiff.pick import pick_types
from .parallel import parallel_func
from .transforms import (
Transform,
_angle_between_quats,
_cart_to_sph,
_get_trans,
_MatchedDisplacementFieldInterpolator,
_pol_to_cart,
apply_trans,
transform_surface_to,
)
from .utils import (
_check_fname,
_check_freesurfer_home,
_check_option,
_ensure_int,
_hashable_ndarray,
_import_nibabel,
_pl,
_soft_import,
_TempDir,
_validate_type,
_verbose_control,
fill_doc_static,
get_subjects_dir,
logger,
run_subprocess,
verbose_static,
warn,
)
_helmet_path = Path(__file__).parent / "data" / "helmets"
@verbose_static("on_defects")
def get_head_surf(
subject, source=("bem", "head"), subjects_dir=None, on_defects="raise", verbose=None
):
"""Load the subject head surface.
Parameters
----------
subject : str
Subject name.
source : str | list of str
Type to load. Common choices would be ``'bem'`` or ``'head'``. We first
try loading ``'$SUBJECTS_DIR/$SUBJECT/bem/$SUBJECT-$SOURCE.fif'``, and
then look for ``'$SUBJECT*$SOURCE.fif'`` in the same directory by going
through all files matching the pattern. The head surface will be read
from the first file containing a head surface. Can also be a list
to try multiple strings.
subjects_dir : path-like | None
Path to the ``SUBJECTS_DIR``. If None, the path is obtained by using
the environment variable ``SUBJECTS_DIR``.
on_defects : 'raise' | 'warn' | 'ignore'
What to do if the surface is found to have topological defects.
Can be ``'raise'`` (default) to raise an error, ``'warn'`` to emit a
warning, or ``'ignore'`` to ignore when one or more defects are found.
Note that a lot of computations in MNE-Python assume the surfaces to be
topologically correct, topological defects may still make other
computations (e.g., `mne.make_bem_model` and `mne.make_bem_solution`)
fail irrespective of this parameter.
.. versionadded:: 1.0
verbose : bool | str | int | None
Control verbosity of the logging output. If ``None``, use the default
verbosity level. See the :ref:`logging documentation <tut-logging>` and
:func:`mne.verbose` for details. Should only be passed as a keyword
argument.
Returns
-------
surf : dict
The head surface.
"""
return _get_head_surface(
subject=subject, source=source, subjects_dir=subjects_dir, on_defects=on_defects
)
def _get_head_surface(subject, source, subjects_dir, on_defects, raise_error=True):
"""Load the subject head surface."""
from .bem import read_bem_surfaces
subjects_dir = str(get_subjects_dir(subjects_dir, raise_error=True))
_validate_type(subject, str, "subject")
if isinstance(source, str):
source = [source]
surf = None
for this_source in source:
this_head = op.realpath(
op.join(subjects_dir, subject, "bem", f"{subject}-{this_source}.fif")
)
if op.exists(this_head):
surf = read_bem_surfaces(
this_head,
True,
FIFF.FIFFV_BEM_SURF_ID_HEAD,
on_defects=on_defects,
verbose=False,
)
else:
path = op.join(subjects_dir, subject, "bem")
if not op.isdir(path):
raise OSError(f'Subject bem directory "{path}" does not exist.')
files = sorted(glob(op.join(path, f"{subject}*{this_source}.fif")))
for this_head in files:
try:
surf = read_bem_surfaces(
this_head,
True,
FIFF.FIFFV_BEM_SURF_ID_HEAD,
on_defects=on_defects,
verbose=False,
)
except ValueError:
pass
else:
break
if surf is not None:
break
if surf is None:
if raise_error:
raise OSError(
f'No file matching "{subject}*{this_source}" and containing a head '
"surface found."
)
else:
return surf
logger.info(f"Using surface from {this_head}.")
return surf
@verbose_static("info_not_none", "helmet_upsampling")
def get_meg_helmet_surf(info, trans=None, *, upsampling=1, verbose=None):
"""Load the MEG helmet associated with the MEG sensors.
Parameters
----------
info : mne.Info
The :class:`mne.Info` object with information about the
sensors and methods of measurement.
trans : dict
The head<->MRI transformation, usually obtained using
read_trans(). Can be None, in which case the surface will
be in head coordinates instead of MRI coordinates.
upsampling : int
The upsampling factor to use for the helmet mesh. The default (1) does no
upsampling. Larger integers lead to more densely sampled helmet surfaces, and
the number of vertices increases as a factor of ``4**(upsampling-1)``.
verbose : bool | str | int | None
Control verbosity of the logging output. If ``None``, use the default
verbosity level. See the :ref:`logging documentation <tut-logging>` and
:func:`mne.verbose` for details. Should only be passed as a keyword
argument.
Returns
-------
surf : dict
The MEG helmet as a surface.
Notes
-----
A built-in helmet is loaded if possible. If not, a helmet surface
will be approximated based on the sensor locations.
"""
from scipy.spatial import ConvexHull, Delaunay
from .bem import _fit_sphere, read_bem_surfaces
from .channels.channels import _get_meg_system
_validate_type(upsampling, "int", "upsampling")
system, have_helmet = _get_meg_system(info)
incomplete = False
if have_helmet:
logger.info(f"Getting helmet for system {system}")
fname = _helmet_path / f"{system}.fif.gz"
surf = read_bem_surfaces(
fname, False, FIFF.FIFFV_MNE_SURF_MEG_HELMET, verbose=False
)
surf = _scale_helmet_to_sensors(system, surf, info)
else:
rr = np.array(
[
info["chs"][pick]["loc"][:3]
for pick in pick_types(info, meg=True, ref_meg=False, exclude=())
]
)
logger.info(
"Getting helmet for system %s (derived from %d MEG channel locations)",
system,
len(rr),
)
hull = ConvexHull(rr)
rr = rr[np.unique(hull.simplices)]
R, center = _fit_sphere(rr)
sph = _cart_to_sph(rr - center)[:, 1:]
front_sph = [[np.pi / 2.0, sph[:, 1].max()]]
sph = np.concatenate((sph, front_sph))
xy = _pol_to_cart(sph[:, ::-1])
tris = Delaunay(xy).simplices
tris = tris[(tris != len(sph) - 1).all(-1)]
tris = _reorder_ccw(rr, tris)
surf = dict(rr=rr, tris=tris)
incomplete = True
if upsampling > 1:
pv = _soft_import("pyvista", "upsample a mesh")
factor = 4 ** (upsampling - 1)
rr, tris = surf["rr"], surf["tris"]
logger.info(
f"Upsampling from {len(rr)} to {len(rr) * factor} vertices ({upsampling=})"
)
tris = np.c_[np.full(len(tris), 3), tris]
mesh = pv.PolyData(rr, tris)
mesh = mesh.subdivide(upsampling - 1, subfilter="linear")
rr, tris = mesh.points, mesh.faces.reshape(-1, 4)[:, 1:]
tris = _reorder_ccw(rr, tris)
surf = dict(rr=rr, tris=tris)
incomplete = True
if incomplete:
complete_surface_info(surf, copy=False, verbose=False)
surf["coord_frame"] = FIFF.FIFFV_COORD_DEVICE
dev_head_t = info["dev_head_t"]
if dev_head_t is None:
dev_head_t = Transform("meg", "head")
transform_surface_to(surf, "head", dev_head_t)
if trans is not None:
transform_surface_to(surf, "mri", trans)
return surf
def _scale_helmet_to_sensors(system, surf, info):
fname = _helmet_path / f"{system}_ch_pos.txt"
if not fname.is_file():
return surf
with open(fname) as fid:
ch_pos_from = json.load(fid)
fro, to = list(), list()
for key, f_ in ch_pos_from.items():
t_ = [ch["loc"][:3] for ch in info["chs"] if ch["ch_name"].startswith(key)]
if not len(t_):
continue
fro.append(f_)
to.append(np.mean(t_, axis=0))
if len(fro) < 4:
logger.info(
"Using untransformed helmet, not enough sensors found to deform to match "
f"acquisition based on sensor positions (got {len(fro)}, need at least 4)"
)
return surf
fro = np.array(fro, float)
to = np.array(to, float)
delta = np.ptp(surf["rr"], axis=0) * 0.1
extrema = np.array([surf["rr"].min(0) - delta, surf["rr"].max(0) + delta])
interp = _MatchedDisplacementFieldInterpolator(fro, to, extrema=extrema)
new_rr = interp(surf["rr"])
try:
from ._transforms_numba import _fit_matched_points
quat, sc = _fit_matched_points(surf["rr"], new_rr)
except np.linalg.LinAlgError as exc:
logger.info(
f"Using untransformed helmet, deformation using {len(fro)} points "
f"failed ({exc})"
)
return surf
rot = np.rad2deg(_angle_between_quats(quat[:3]))
tr = 1000 * np.linalg.norm(quat[3:])
logger.info(
f" Deforming CAD helmet to match {len(fro)} acquisition sensor positions:"
)
logger.info(f" 1. Affine: {rot:0.1f}°, {tr:0.1f} mm, {sc:0.2f}× scale")
deltas = interp._last_deltas * 1000
mu, mx = np.mean(deltas), np.max(deltas)
logger.info(f" 2. Nonlinear displacement: mean={mu:0.1f}, max={mx:0.1f} mm")
surf["rr"] = new_rr
complete_surface_info(surf, copy=False, verbose=False)
return surf
def _reorder_ccw(rrs, tris):
"""Reorder tris of a convex hull to be wound counter-clockwise."""
com = np.mean(rrs, axis=0)
rr_tris = rrs[tris]
dirs = (
np.cross(rr_tris[:, 1] - rr_tris[:, 0], rr_tris[:, 2] - rr_tris[:, 0])
* (rr_tris[:, 0] - com)
).sum(-1)
return np.where(dirs[:, np.newaxis] < 0, tris[:, ::-1], tris)
def fast_cross_3d(x, y):
"""Compute cross product between list of 3D vectors.
Equivalent to :func:`numpy.cross`.
Parameters
----------
x : array
Input array 1, shape (..., 3).
y : array
Input array 2, shape (..., 3).
Returns
-------
z : array, shape (..., 3)
Cross product of x and y along the last dimension.
Notes
-----
x and y must broadcast against each other.
"""
assert x.ndim >= 1
assert y.ndim >= 1
assert x.shape[-1] == 3
assert y.shape[-1] == 3
return np.cross(x, y)
def _triangle_neighbors(tris, npts):
"""Efficiently compute vertex neighboring triangles."""
from scipy.sparse import coo_array
rows = tris.ravel()
cols = np.repeat(np.arange(len(tris)), 3)
data = np.ones(len(cols))
csr = coo_array((data, (rows, cols)), shape=(npts, len(tris))).tocsr()
neighbor_tri = [
csr.indices[start:stop] for start, stop in zip(csr.indptr[:-1], csr.indptr[1:])
]
assert len(neighbor_tri) == npts
return neighbor_tri
def _project_onto_surface(
rrs, surf, project_rrs=False, return_nn=False, method="accurate"
):
"""Project points onto (scalp) surface."""
from ._surface_numba import _accumulate_normals, _find_nearest_tri_pts
if method == "accurate":
surf_geom = _get_tri_supp_geom(surf)
pt_tris = np.empty((0,), int)
pt_lens = np.zeros(len(rrs) + 1, int)
out = _find_nearest_tri_pts(rrs, pt_tris, pt_lens, reproject=True, **surf_geom)
if project_rrs:
out += (np.einsum("ij,ijk->ik", out[0], surf["rr"][surf["tris"][out[1]]]),)
if return_nn:
out += (surf_geom["nn"][out[1]],)
else:
assert project_rrs
idx = _compute_nearest(surf["rr"], rrs)
out = (None, None, surf["rr"][idx])
if return_nn:
surf_geom = _get_tri_supp_geom(surf)
nn = _accumulate_normals(
surf["tris"].astype(int), surf_geom["nn"], len(surf["rr"])
)
nn = nn[idx]
_normalize_vectors(nn)
out += (nn,)
return out
def _normal_orth(nn):
"""Compute orthogonal basis given a normal."""
assert nn.shape[-1:] == (3,)
prod = np.einsum("...i,...j->...ij", nn, nn)
_, u = np.linalg.eigh(np.eye(3) - prod)
u = u[..., ::-1]
signs = np.sign(np.matmul(nn[..., np.newaxis, :], u[..., -1:]))
signs[signs == 0] = 1
u *= signs
return u.swapaxes(-1, -2)
@verbose_static()
def complete_surface_info(
surf, do_neighbor_vert=False, copy=True, do_neighbor_tri=True, *, verbose=None
):
"""Complete surface information.
Parameters
----------
surf : dict
The surface.
do_neighbor_vert : bool
If True (default False), add neighbor vertex information.
copy : bool
If True (default), make a copy. If False, operate in-place.
do_neighbor_tri : bool
If True (default), compute triangle neighbors.
verbose : bool | str | int | None
Control verbosity of the logging output. If ``None``, use the default
verbosity level. See the :ref:`logging documentation <tut-logging>` and
:func:`mne.verbose` for details. Should only be passed as a keyword
argument.
Returns
-------
surf : dict
The transformed surface.
"""
if copy:
surf = deepcopy(surf)
surf["ntri"] = surf.get("ntri", len(surf["tris"]))
surf["np"] = surf.get("np", len(surf["rr"]))
surf["tri_area"] = np.zeros(surf["ntri"])
r1 = surf["rr"][surf["tris"][:, 0], :]
r2 = surf["rr"][surf["tris"][:, 1], :]
r3 = surf["rr"][surf["tris"][:, 2], :]
surf["tri_cent"] = (r1 + r2 + r3) / 3.0
surf["tri_nn"] = fast_cross_3d((r2 - r1), (r3 - r1))
surf["tri_area"] = _normalize_vectors(surf["tri_nn"]) / 2.0
zidx = np.where(surf["tri_area"] == 0)[0]
if len(zidx) > 0:
logger.info(f" Warning: zero size triangles: {zidx}")
from ._surface_numba import _accumulate_normals
logger.info(" Triangle neighbors and vertex normals...")
surf["nn"] = _accumulate_normals(
surf["tris"].astype(int), surf["tri_nn"], surf["np"]
)
_normalize_vectors(surf["nn"])
if do_neighbor_tri:
surf["neighbor_tri"] = _triangle_neighbors(surf["tris"], surf["np"])
zero, fewer = list(), list()
for ni, n in enumerate(surf["neighbor_tri"]):
if len(n) < 3:
if len(n) == 0:
zero.append(ni)
else:
fewer.append(ni)
surf["neighbor_tri"][ni] = np.array([], int)
if len(zero) > 0:
logger.info(
" Vertices do not have any neighboring triangles: "
f"[{', '.join(str(z) for z in zero)}]"
)
if len(fewer) > 0:
fewer = ", ".join(str(f) for f in fewer)
logger.info(
" Vertices have fewer than three neighboring triangles, removing "
f"neighbors: [{fewer}]"
)
if do_neighbor_vert is True:
logger.info(" Vertex neighbors...")
edges = mesh_edges(surf["tris"])
n_vert = np.diff(edges.indptr)
n_tri = np.array([len(n) for n in surf["neighbor_tri"]])
bad = np.where(n_vert > n_tri)[0]
if len(bad):
raise RuntimeError(f"Too many neighbors for vertex {bad[0]}.")
for k in np.where(n_vert != n_tri)[0]:
logger.info(
" Incorrect number of distinct neighbors for vertex"
" %d (%d instead of %d) [fixed].",
k,
n_vert[k],
n_tri[k],
)
indices = edges.indices.copy()
surf["neighbor_vert"] = [
indices[start:stop]
for start, stop in zip(edges.indptr[:-1], edges.indptr[1:])
]
return surf
def _get_surf_neighbors(surf, k):
"""Calculate the surface neighbors based on triangulation."""
verts = set()
for v in surf["tris"][surf["neighbor_tri"][k]].flat:
verts.add(v)
verts.remove(k)
verts = np.array(sorted(verts))
assert np.all(verts < surf["np"])
nneighbors = len(verts)
nneigh_max = len(surf["neighbor_tri"][k])
if nneighbors > nneigh_max:
raise RuntimeError(f"Too many neighbors for vertex {k}.")
elif nneighbors != nneigh_max:
logger.info(
" Incorrect number of distinct neighbors for vertex"
" %d (%d instead of %d) [fixed].",
k,
nneighbors,
nneigh_max,
)
return verts
def _normalize_vectors(rr):
"""Normalize surface vertices."""
size = np.linalg.norm(rr, axis=1)
np.divide(rr, size[:, np.newaxis], out=rr, where=size[:, np.newaxis] > 0)
return size
class _CDist:
"""Wrapper for cdist that uses a Tree-like pattern."""
def __init__(self, xhs):
self._xhs = xhs
def query(self, rr):
from scipy.spatial.distance import cdist
nearest = list()
dists = list()
for r in rr:
d = cdist(r[np.newaxis, :], self._xhs)
idx = np.argmin(d)
nearest.append(idx)
dists.append(d[0, idx])
return np.array(dists), np.array(nearest)
def _compute_nearest(xhs, rr, method="BallTree", return_dists=False):
"""Find nearest neighbors.
Parameters
----------
xhs : array, shape=(n_samples, n_dim)
Points of data set.
rr : array, shape=(n_query, n_dim)
Points to find nearest neighbors for.
method : str
The query method. If scikit-learn and scipy<1.0 are installed,
it will fall back to the slow brute-force search.
return_dists : bool
If True, return associated distances.
Returns
-------
nearest : array, shape=(n_query,)
Index of nearest neighbor in xhs for every point in rr.
distances : array, shape=(n_query,)
The distances. Only returned if return_dists is True.
"""
if xhs.size == 0 or rr.size == 0:
if return_dists:
return np.array([], int), np.array([])
return np.array([], int)
tree = _DistanceQuery(xhs, method=method)
out = tree.query(rr)
return out[::-1] if return_dists else out[1]
def _safe_query(rr, func, reduce=False, **kwargs):
if len(rr) == 0:
return np.array([]), np.array([], int)
out = func(rr)
out = [out[0][:, 0], out[1][:, 0]] if reduce else out
return out
class _DistanceQuery:
"""Wrapper for fast distance queries."""
def __init__(self, xhs, method="BallTree"):
assert method in ("BallTree", "KDTree", "cdist")
if method == "BallTree":
try:
from sklearn.neighbors import BallTree
except ImportError:
logger.info(
"Nearest-neighbor searches will be significantly "
"faster if scikit-learn is installed."
)
method = "KDTree"
else:
self.query = partial(
_safe_query,
func=BallTree(xhs).query,
reduce=True,
return_distance=True,
)
if method == "KDTree":
from scipy.spatial import KDTree
self.query = KDTree(xhs).query
if method == "cdist":
self.query = _CDist(xhs).query
self.data = xhs
@_verbose_control
def _points_outside_surface(rr, surf, n_jobs=None, verbose=None):
"""Check whether points are outside a surface.
Parameters
----------
rr : ndarray
Nx3 array of points to check.
surf : dict
Surface with entries "rr" and "tris".
Returns
-------
outside : ndarray
1D logical array of size N for which points are outside the surface.
"""
from ._surface_numba import _get_solids
rr = np.atleast_2d(rr)
assert rr.shape[1] == 3
parallel, p_fun, n_jobs = parallel_func(_get_solids, n_jobs)
tot_angles = parallel(
p_fun(surf["rr"][tris], rr) for tris in np.array_split(surf["tris"], n_jobs)
)
return np.abs(np.sum(tot_angles, axis=0) / (2 * np.pi) - 1.0) > 1e-5
def _surface_to_polydata(surf):
import pyvista as pv
vertices = np.array(surf["rr"])
if "tris" not in surf:
return pv.PolyData(vertices)
else:
triangles = np.array(surf["tris"])
triangles = np.c_[np.full(len(triangles), 3), triangles]
return pv.PolyData(vertices, triangles)
def _polydata_to_surface(pd, normals=True):
from pyvista import PolyData
if not isinstance(pd, PolyData):
pd = PolyData(pd)
out = dict(rr=pd.points, tris=pd.faces.reshape(-1, 4)[:, 1:])
if normals:
out["nn"] = pd.point_normals
return out
class _CheckInside:
"""Efficiently check if points are inside a surface."""
@_verbose_control
def __init__(self, surf, *, mode="old", verbose=None):
assert mode in ("pyvista", "old")
self.mode = mode
t0 = time.time()
self.surf = surf
if self.mode == "pyvista":
self._init_pyvista()
else:
self._init_old()
logger.debug(
f"Setting up {mode} interior check for {len(self.surf['rr'])} "
f"points took {(time.time() - t0) * 1000:0.1f} ms"
)
def _init_old(self):
from scipy.spatial import Delaunay
self.inner_r = None
self.center = self.surf["rr"].mean(0)
self.del_tri = Delaunay(self.surf["rr"])
if self.del_tri.find_simplex(self.center) >= 0:
dists = np.linalg.norm(self.surf["rr"] - self.center, axis=-1)
self.inner_r = dists.min()
self.outer_r = dists.max()
def _init_pyvista(self):
if not isinstance(self.surf, dict):
self.pdata = self.surf
self.surf = _polydata_to_surface(self.pdata)
else:
self.pdata = _surface_to_polydata(self.surf).clean()
@_verbose_control
def __call__(self, rr, *, n_jobs=None, verbose=None):
n_orig = len(rr)
logger.info(
f"Checking surface interior status for {n_orig} point{_pl(n_orig, ' ')}..."
)
t0 = time.time()
if self.mode == "pyvista":
inside = self._call_pyvista(rr)
else:
inside = self._call_old(rr, n_jobs)
n = inside.sum()
logger.info(f" Total {n}/{n_orig} point{_pl(n, ' ')} inside the surface")
logger.info(f"Interior check completed in {(time.time() - t0) * 1000:0.1f} ms")
return inside
def query(self, rr):
"""Get the distance to the nearest point."""
if not hasattr(self, "_tree"):
from scipy.spatial import KDTree
self._tree = KDTree(self.surf["rr"])
return self._tree.query(rr)
def _call_pyvista(self, rr):
pdata = _surface_to_polydata(dict(rr=rr))
if hasattr(pdata, "select_interior_points"):
meth = pdata.select_interior_points
key = "selected_points"
else:
meth = pdata.select_enclosed_points
key = "SelectedPoints"
out = meth(self.pdata, check_surface=False)
return out[key].astype(bool)
def _call_old(self, rr, n_jobs):
n_orig = len(rr)
prec = int(np.ceil(np.log10(max(n_orig, 10))))
inside = np.ones(n_orig, bool)
idx = np.arange(n_orig)
if self.inner_r is not None:
dists = np.linalg.norm(rr - self.center, axis=-1)
in_mask = dists < self.inner_r
n = (in_mask).sum()
n_pad = str(n).rjust(prec)
logger.info(
f" Found {n_pad}/{n_orig} point{_pl(n, ' ')} "
f"inside an interior sphere of radius "
f"{1000 * self.inner_r:6.1f} mm"
)
out_mask = dists > self.outer_r
inside[out_mask] = False
n = (out_mask).sum()
n_pad = str(n).rjust(prec)
logger.info(
f" Found {n_pad}/{n_orig} point{_pl(n, ' ')} "
f"outside an exterior sphere of radius "
f"{1000 * self.outer_r:6.1f} mm"
)
mask = (~in_mask) & (~out_mask)
idx = idx[mask]
rr = rr[mask]
del_outside = self.del_tri.find_simplex(rr) < 0
n = del_outside.sum()
inside[idx[del_outside]] = False
idx = idx[~del_outside]
rr = rr[~del_outside]
n_pad = str(n).rjust(prec)
check_pad = str(len(del_outside)).rjust(prec)
logger.info(
f" Found {n_pad}/{check_pad} point{_pl(n, ' ')} outside using "
"surface Qhull"
)
solid_outside = _points_outside_surface(rr, self.surf, n_jobs)
n = np.sum(solid_outside)
n_pad = str(n).rjust(prec)
check_pad = str(len(solid_outside)).rjust(prec)
logger.info(
f" Found {n_pad}/{check_pad} point{_pl(n, ' ')} outside using "
"solid angles"
)
inside[idx[solid_outside]] = False
return inside
class _CheckInsideSphere:
def __init__(self, sphere, *, check="inner"):
from .bem import ConductorModel
assert isinstance(sphere, ConductorModel) and sphere["is_sphere"]
self.center = sphere["r0"]
assert isinstance(check, str) and check in ("inner", "outer"), check
self.check = check
if len(sphere["layers"]):
self.inner_r = sphere["layers"][0]["rad"]
self.outer_r = sphere["layers"][-1]["rad"]
else:
self.inner_r = self.outer_r = None
def __call__(self, rr, *, n_jobs=None, verbose=None):
assert isinstance(rr, np.ndarray), type(rr)
assert rr.ndim == 2 and rr.shape[1] == 3
if self.inner_r is None:
return np.ones(rr.shape[0], bool)
else:
return np.linalg.norm(rr - self.center, axis=-1) <= self._check_r
@property
def _check_r(self):
return self.inner_r if self.check == "inner" else self.outer_r
def query(self, rr):
"""Return the distance to the sphere surface for each point."""
assert isinstance(rr, np.ndarray), type(rr)
assert rr.ndim == 2 and rr.shape[1] == 3, rr.shape
idx = np.zeros(rr.shape[0], int)
if self.inner_r is None:
dists = np.full(rr.shape[0], np.inf)
else:
dists = np.abs(np.linalg.norm(rr - self.center, axis=-1) - self._check_r)
return dists, idx
def _fread3(fobj):
"""Read 3 bytes and adjust."""
b1, b2, b3 = np.fromfile(fobj, ">u1", 3).astype(np.int64)
return (b1 << 16) + (b2 << 8) + b3
def read_curvature(filepath, binary=True):
"""Load in curvature values from the ?h.curv file.
Parameters
----------
filepath : path-like
Input path to the ``.curv`` file.
binary : bool
Specify if the output array is to hold binary values. Defaults to True.
Returns
-------
curv : array of shape (n_vertices,)
The curvature values loaded from the user given file.
"""
filepath = _check_fname(filepath, "read", must_exist=True, name="Curvature file")
with open(filepath, "rb") as fobj:
magic = _fread3(fobj)
if magic == 16777215:
vnum = np.fromfile(fobj, ">i4", 3)[0]
curv = np.fromfile(fobj, ">f4", vnum)
else:
vnum = magic
_fread3(fobj)
curv = np.fromfile(fobj, ">i2", vnum) / 100
if binary:
return 1 - np.array(curv != 0, np.int64)
else:
return curv
@verbose_static()
def read_surface(
fname, read_metadata=False, return_dict=False, file_format="auto", verbose=None
):
"""Load a FreeSurfer surface mesh in triangular format.
Parameters
----------
fname : path-like
The name of the file containing the surface.
read_metadata : bool
Read metadata as key-value pairs. Only works when reading a FreeSurfer
surface file. For .obj files this dictionary will be empty.
Valid keys:
* 'head' : array of int
* 'valid' : str
* 'filename' : str
* 'volume' : array of int, shape (3,)
* 'voxelsize' : array of float, shape (3,)
* 'xras' : array of float, shape (3,)
* 'yras' : array of float, shape (3,)
* 'zras' : array of float, shape (3,)
* 'cras' : array of float, shape (3,)
.. versionadded:: 0.13.0
return_dict : bool
If True, a dictionary with surface parameters is returned.
file_format : 'auto' | 'freesurfer' | 'obj'
File format to use. Can be 'freesurfer' to read a FreeSurfer surface
file, or 'obj' to read a Wavefront .obj file (common format for
importing in other software), or 'auto' to attempt to infer from the
file name. Defaults to 'auto'.
.. versionadded:: 0.21.0
verbose : bool | str | int | None
Control verbosity of the logging output. If ``None``, use the default
verbosity level. See the :ref:`logging documentation <tut-logging>` and
:func:`mne.verbose` for details. Should only be passed as a keyword
argument.
Returns
-------
rr : array, shape=(n_vertices, 3)
Coordinate points.
tris : int array, shape=(n_faces, 3)
Triangulation (each line contains indices for three points which
together form a face).
volume_info : dict-like
If read_metadata is true, key-value pairs found in the geometry file.
surf : dict
The surface parameters. Only returned if ``return_dict`` is True.
See Also
--------
write_surface
read_tri
"""
fname = _check_fname(fname, "read", True)
_check_option("file_format", file_format, ["auto", "freesurfer", "obj"])
if file_format == "auto":
if fname.suffix.lower() == ".obj":
file_format = "obj"
else:
file_format = "freesurfer"
if file_format == "freesurfer":
_import_nibabel("read surface geometry")
from nibabel.freesurfer import read_geometry
ret = read_geometry(fname, read_metadata=read_metadata)
elif file_format == "obj":
ret = _read_wavefront_obj(fname)
if read_metadata:
ret += (dict(),)
if return_dict:
ret += (_rr_tris_dict(ret[0], ret[1]),)
return ret
def _rr_tris_dict(rr, tris):
return dict(rr=rr, tris=tris, ntri=len(tris), use_tris=tris, np=len(rr))
def _read_mri_surface(fname):
surf = read_surface(fname, return_dict=True)[2]
surf["rr"] /= 1000.0
surf.update(coord_frame=FIFF.FIFFV_COORD_MRI)
return surf
def _read_wavefront_obj(fname):
"""Read a surface form a Wavefront .obj file.
Parameters
----------
fname : str
Name of the .obj file to read.
Returns
-------
coords : ndarray, shape (n_points, 3)
The XYZ coordinates of each vertex.
faces : ndarray, shape (n_faces, 3)
For each face of the mesh, the integer indices of the vertices that
make up the face.
"""
coords = []
faces = []
with open(fname) as f:
for line in f:
line = line.strip()
if len(line) == 0 or line[0] == "#":
continue
split = line.split()
if split[0] == "v":
coords.append([float(item) for item in split[1:]])
elif split[0] == "f":
dat = [int(item.split("/")[0]) for item in split[1:]]
if len(dat) != 3:
raise RuntimeError("Only triangle faces allowed.")
faces.append([d - 1 for d in dat])
return np.array(coords), np.array(faces)
def _read_patch(fname):
"""Load a FreeSurfer binary patch file.
Parameters
----------
fname : str
The filename.
Returns
-------
rrs : ndarray, shape (n_vertices, 3)
The points.
tris : ndarray, shape (n_tris, 3)
The patches. Not all vertices will be present.
"""
patch = dict()
with open(fname) as fid:
ver = np.fromfile(fid, dtype=">i4", count=1).item()
if ver != -1:
raise RuntimeError(f"incorrect version # {ver} (not -1) found")
npts = np.fromfile(fid, dtype=">i4", count=1).item()
dtype = np.dtype([("vertno", ">i4"), ("x", ">f"), ("y", ">f"), ("z", ">f")])
recs = np.fromfile(fid, dtype=dtype, count=npts)
patch = {key: recs[key] for key in dtype.fields.keys()}
patch["vertno"] -= 1
rrs, tris = read_surface(
op.join(op.dirname(fname), op.basename(fname)[:3] + "sphere")
)
orig_tris = tris
is_vert = patch["vertno"] > 0
verts = patch["vertno"][is_vert]
mask = np.zeros((len(rrs),), dtype=bool)
mask[verts] = True
rrs[~mask] = 0.0
tris = tris[mask[tris].all(1)]
for ii, key in enumerate(["x", "y", "z"]):
rrs[verts, ii] = patch[key][is_vert]
return rrs, tris, orig_tris
def _get_ico_surface(grade, patch_stats=False):
"""Return an icosahedral surface of the desired grade."""
from .bem import read_bem_surfaces
ico_file_name = op.join(op.dirname(__file__), "data", "icos.fif.gz")
ico = read_bem_surfaces(
ico_file_name, patch_stats, s_id=9000 + grade, verbose=False
)
return ico
def _tessellate_sphere_surf(level, rad=1.0):
"""Return a surface structure instead of the details."""
rr, tris = _tessellate_sphere(level)
npt = len(rr)
ntri = len(tris)
nn = rr.copy()
rr *= rad
s = dict(
rr=rr,
np=npt,
tris=tris,
use_tris=tris,
ntri=ntri,
nuse=npt,
nn=nn,
inuse=np.ones(npt, int),
)
return s
def _norm_midpt(ai, bi, rr):
"""Get normalized midpoint."""
c = rr[ai]
c += rr[bi]
_normalize_vectors(c)
return c
def _tessellate_sphere(mylevel):
"""Create a tessellation of a unit sphere."""
rr = np.array(
[
[1, 0, 0],
[-1, 0, 0],
[0, 1, 0],
[0, -1, 0],
[0, 0, 1],
[0, 0, -1],
],
float,
)
tris = np.array(
[
[0, 4, 2],
[2, 4, 1],
[1, 4, 3],
[3, 4, 0],
[0, 2, 5],
[2, 1, 5],
[1, 3, 5],
[3, 0, 5],
],
int,
)
if mylevel < 1:
raise ValueError("oct subdivision must be >= 1")
tris = tris[:, [2, 1, 0]]
for _ in range(1, mylevel):
r"""
Subdivide each triangle in the old approximation and normalize
the new points thus generated to lie on the surface of the unit
sphere.
Each input triangle with vertices labelled [0,1,2] as shown
below will be turned into four new triangles:
Make new points
a = (0+2)/2
b = (0+1)/2
c = (1+2)/2
1
/\ Normalize a, b, c
/ \
b/____\c Construct new triangles
/\ /\ [0,b,a]
/ \ / \ [b,1,c]
/____\/____\ [a,b,c]
0 a 2 [a,c,2]
"""
a = _norm_midpt(tris[:, 0], tris[:, 2], rr)
b = _norm_midpt(tris[:, 0], tris[:, 1], rr)
c = _norm_midpt(tris[:, 1], tris[:, 2], rr)
lims = np.cumsum([len(rr), len(a), len(b), len(c)])
aidx = np.arange(lims[0], lims[1])
bidx = np.arange(lims[1], lims[2])
cidx = np.arange(lims[2], lims[3])
rr = np.concatenate((rr, a, b, c))
tris = np.array(
(
np.c_[tris[:, 0], bidx, aidx],
np.c_[bidx, tris[:, 1], cidx],
np.c_[aidx, bidx, cidx],
np.c_[aidx, cidx, tris[:, 2]],
),
int,
).swapaxes(0, 1)
tris = np.reshape(tris, (np.prod(tris.shape[:2]), 3))
flat = tris.ravel()
_, first, inverse = np.unique(
rr[flat], axis=0, return_index=True, return_inverse=True
)
order = np.argsort(first)
tris = np.argsort(order)[inverse].reshape(tris.shape)
rr = rr[flat[first[order]]]
return rr, tris
def _decimate_surface_ico_oct(subject, subjects_dir, hemi, surf, spacing):
from .source_space._source_space import _check_spacing
stype, _, ico_surf, _ = _check_spacing(spacing, verbose=False)
subjects_dir = Path(subjects_dir)
surf_fname = subjects_dir / subject / "surf" / f"{hemi}.{surf}"
dec = _create_surf_spacing(surf_fname, hemi, subject, stype, ico_surf, subjects_dir)
vertno, use_tris = dec["vertno"], dec["use_tris"]
lut = np.zeros(dec["np"], int)
lut[vertno] = np.arange(len(vertno))
tris = lut[use_tris]
return vertno, tris
def _create_surf_spacing(surf, hemi, subject, stype, ico_surf, subjects_dir):
"""Load a surf and use the subdivided icosahedron to get points."""
surf = read_surface(surf, return_dict=True)[-1]
do_neighbor_vert = stype == "spacing"
complete_surface_info(surf, do_neighbor_vert, copy=False)
if stype == "all":
surf["inuse"] = np.ones(surf["np"], int)
surf["use_tris"] = None
elif stype == "spacing":
_decimate_surface_spacing(surf, ico_surf)
surf["use_tris"] = None
del surf["neighbor_vert"]
else:
surf_name = subjects_dir / subject / "surf" / f"{hemi}.sphere"
logger.info(f"Loading geometry from {surf_name}...")
from_surf = read_surface(surf_name, return_dict=True)[-1]
_normalize_vectors(from_surf["rr"])
if from_surf["np"] != surf["np"]:
raise RuntimeError(
"Mismatch between number of surface vertices, "
"possible parcellation error?"
)
_normalize_vectors(ico_surf["rr"])
mmap = _compute_nearest(from_surf["rr"], ico_surf["rr"])
nmap = len(mmap)
surf["inuse"] = np.zeros(surf["np"], int)
for k in range(nmap):
if surf["inuse"][mmap[k]]:
neigh = _get_surf_neighbors(surf, mmap[k])
was = mmap[k]
inds = np.where(np.logical_not(surf["inuse"][neigh]))[0]
if len(inds) == 0:
raise RuntimeError(
f"Could not find neighbor for vertex {k} / {nmap}."
)
else:
mmap[k] = neigh[inds[-1]]
logger.info(
" Source space vertex moved from %d to %d "
"because of double occupation",
was,
mmap[k],
)
elif mmap[k] < 0 or mmap[k] > surf["np"]:
raise RuntimeError(
f"Map number out of range ({mmap[k]}), this is probably due to "
"inconsistent surfaces. Parts of the FreeSurfer reconstruction "
"need to be redone."
)
surf["inuse"][mmap[k]] = True
logger.info("Setting up the triangulation for the decimated surface...")
surf["use_tris"] = mmap[ico_surf["tris"]].astype(np.int32)
if surf["use_tris"] is not None:
surf["nuse_tri"] = len(surf["use_tris"])
else:
surf["nuse_tri"] = 0
surf["nuse"] = np.sum(surf["inuse"])
surf["vertno"] = np.where(surf["inuse"])[0]
sizes = _normalize_vectors(surf["nn"])
surf["inuse"][sizes <= 0] = False
surf["nuse"] = np.sum(surf["inuse"])
surf["subject_his_id"] = subject
return surf
def _decimate_surface_spacing(surf, spacing):
assert isinstance(spacing, int)
assert spacing > 0
logger.info(" Decimating...")
from ._surface_numba import _decimate_spacing
edges = mesh_edges(surf["tris"])
inuse = _decimate_spacing(edges.indptr, edges.indices, spacing)
surf["inuse"] = inuse.astype(int)
return surf
@verbose_static("overwrite")
def write_surface(
fname,
coords,
faces,
create_stamp="",
volume_info=None,
file_format="auto",
overwrite=False,
*,
verbose=None,
):
"""Write a triangular FreeSurfer surface mesh.
Accepts the same data format as is returned by read_surface().
Parameters
----------
fname : path-like
File to write.
coords : array, shape=(n_vertices, 3)
Coordinate points.
faces : int array, shape=(n_faces, 3)
Triangulation (each line contains indices for three points which
together form a face).
create_stamp : str
Comment that is written to the beginning of the file. Can not contain
line breaks.
volume_info : dict-like or None
Key-value pairs to encode at the end of the file.
Valid keys:
* 'head' : array of int
* 'valid' : str
* 'filename' : str
* 'volume' : array of int, shape (3,)
* 'voxelsize' : array of float, shape (3,)
* 'xras' : array of float, shape (3,)
* 'yras' : array of float, shape (3,)
* 'zras' : array of float, shape (3,)
* 'cras' : array of float, shape (3,)
.. versionadded:: 0.13.0
file_format : 'auto' | 'freesurfer' | 'obj'
File format to use. Can be 'freesurfer' to write a FreeSurfer surface
file, or 'obj' to write a Wavefront .obj file (common format for
importing in other software), or 'auto' to attempt to infer from the
file name. Defaults to 'auto'.
.. versionadded:: 0.21.0
overwrite : bool
If True (default False), overwrite the destination file if it
exists.
verbose : bool | str | int | None
Control verbosity of the logging output. If ``None``, use the default
verbosity level. See the :ref:`logging documentation <tut-logging>` and
:func:`mne.verbose` for details. Should only be passed as a keyword
argument.
See Also
--------
read_surface
read_tri
"""
fname = _check_fname(fname, overwrite=overwrite)
_check_option("file_format", file_format, ["auto", "freesurfer", "obj"])
if file_format == "auto":
if fname.suffix.lower() == ".obj":
file_format = "obj"
else:
file_format = "freesurfer"
if file_format == "freesurfer":
_import_nibabel("write surface geometry")
from nibabel.freesurfer import write_geometry
write_geometry(
fname, coords, faces, create_stamp=create_stamp, volume_info=volume_info
)
else:
assert file_format == "obj"
with open(fname, "w") as fid:
for line in create_stamp.splitlines():
fid.write(f"# {line}\n")
for v in coords:
fid.write(f"v {v[0]} {v[1]} {v[2]}\n")
for f in faces:
fid.write(f"f {f[0] + 1} {f[1] + 1} {f[2] + 1}\n")
def _decimate_surface_vtk(points, triangles, n_triangles):
"""Aux function."""
try:
from vtkmodules.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray
from vtkmodules.vtkCommonCore import vtkPoints
from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkPolyData
from vtkmodules.vtkFiltersCore import vtkQuadricDecimation
except ImportError:
raise ValueError("This function requires the VTK package to be installed")
if triangles.max() > len(points) - 1:
raise ValueError(
"The triangles refer to undefined points. Please check your mesh."
)
src = vtkPolyData()
vtkpoints = vtkPoints()
vtkpoints.SetData(numpy_to_vtk(points.astype(np.float64)))
src.SetPoints(vtkpoints)
vtkcells = vtkCellArray()
triangles_ = np.pad(triangles, ((0, 0), (1, 0)), "constant", constant_values=3)
idarr = numpy_to_vtkIdTypeArray(triangles_.ravel().astype(np.int64))
vtkcells.ImportLegacyFormat(idarr)
src.SetPolys(vtkcells)
decimate = vtkQuadricDecimation()
decimate.VolumePreservationOn()
decimate.SetInputData(src)
reduction = 1 - (float(n_triangles) / len(triangles))
decimate.SetTargetReduction(reduction)
decimate.Update()
out = _polydata_to_surface(decimate.GetOutput(), normals=False)
return out["rr"], out["tris"]
def _decimate_surface_sphere(rr, tris, n_triangles):
_check_freesurfer_home()
map_ = {}
ico_levels = [20, 80, 320, 1280, 5120, 20480]
map_.update({n_tri: ("ico", ii) for ii, n_tri in enumerate(ico_levels)})
oct_levels = 2 ** (2 * np.arange(7) + 3)
map_.update({n_tri: ("oct", ii) for ii, n_tri in enumerate(oct_levels, 1)})
_check_option(
"n_triangles", n_triangles, sorted(map_), extra=' when method="sphere"'
)
func_map = dict(ico=_get_ico_surface, oct=_tessellate_sphere_surf)
kind, level = map_[n_triangles]
logger.info(f"Decimating using FreeSurfer spherical {kind}{level} downsampling")
ico_surf = func_map[kind](level)
assert len(ico_surf["tris"]) == n_triangles
tempdir = _TempDir()
orig = op.join(tempdir, "lh.temp")
write_surface(orig, rr, tris)
logger.info(" Extracting main mesh component ...")
run_subprocess(["mris_extract_main_component", orig, orig], verbose="error")
logger.info(" Smoothing ...")
smooth = orig + ".smooth"
run_subprocess(["mris_smooth", "-nw", orig, smooth], verbose="error")
logger.info(" Inflating ...")
inflated = orig + ".inflated"
run_subprocess(["mris_inflate", "-no-save-sulc", smooth, inflated], verbose="error")
logger.info(" Sphere ...")
qsphere = orig + ".qsphere"
run_subprocess(["mris_sphere", "-q", inflated, qsphere], verbose="error")
sphere_rr, _ = read_surface(qsphere)
norms = np.linalg.norm(sphere_rr, axis=1, keepdims=True)
sphere_rr /= norms
idx = _compute_nearest(sphere_rr, ico_surf["rr"], method="KDTree")
n_dup = len(idx) - len(np.unique(idx))
if n_dup:
raise RuntimeError(
f"Could not reduce to {n_triangles} triangles using ico, "
f"{n_dup}/{len(idx)} vertices were duplicates."
)
logger.info("[done]")
return rr[idx], ico_surf["tris"]
@verbose_static()
def decimate_surface(points, triangles, n_triangles, method="quadric", *, verbose=None):
"""Decimate surface data.
Parameters
----------
points : ndarray
The surface to be decimated, a 3 x number of points array.
triangles : ndarray
The surface to be decimated, a 3 x number of triangles array.
n_triangles : int
The desired number of triangles.
method : str
Can be "quadric" or "sphere". "sphere" will inflate the surface to a
sphere using FreeSurfer and downsample to an icosahedral or
octahedral mesh.
.. versionadded:: 0.20
verbose : bool | str | int | None
Control verbosity of the logging output. If ``None``, use the default
verbosity level. See the :ref:`logging documentation <tut-logging>` and
:func:`mne.verbose` for details. Should only be passed as a keyword
argument.
Returns
-------
points : ndarray
The decimated points.
triangles : ndarray
The decimated triangles.
Notes
-----
**"quadric" mode**
This requires VTK. If an odd target number was requested,
the ``'decimation'`` algorithm used results in the
next even number of triangles. For example a reduction request
to 30001 triangles may result in 30000 triangles.
**"sphere" mode**
This requires FreeSurfer to be installed and available in the
environment. The destination number of triangles must be one of
``[20, 80, 320, 1280, 5120, 20480]`` for ico (0-5) downsampling or one of
``[8, 32, 128, 512, 2048, 8192, 32768]`` for oct (1-7) downsampling.
This mode is slower, but could be more suitable for decimating meshes for
BEM creation (recommended ``n_triangles=5120``) due to better topological
property preservation.
"""
n_triangles = _ensure_int(n_triangles)
method_map = dict(quadric=_decimate_surface_vtk, sphere=_decimate_surface_sphere)
_check_option("method", method, sorted(method_map))
if n_triangles > len(triangles):
raise ValueError(
f"Requested n_triangles ({n_triangles}) exceeds number of "
f"original triangles ({len(triangles)})"
)
return method_map[method](points, triangles, n_triangles)
def _get_tri_supp_geom(surf):
"""Create supplementary geometry information using tris and rrs."""
r1 = surf["rr"][surf["tris"][:, 0], :]
r12 = surf["rr"][surf["tris"][:, 1], :] - r1
r13 = surf["rr"][surf["tris"][:, 2], :] - r1
r1213 = np.ascontiguousarray(np.array([r12, r13]).swapaxes(0, 1))
a = np.einsum("ij,ij->i", r12, r12)
b = np.einsum("ij,ij->i", r13, r13)
c = np.einsum("ij,ij->i", r12, r13)
mat = np.ascontiguousarray(np.rollaxis(np.array([[b, -c], [-c, a]]), 2))
norm = a * b - c * c
norm[norm == 0] = 1.0
mat /= norm[:, np.newaxis, np.newaxis]
nn = fast_cross_3d(r12, r13)
_normalize_vectors(nn)
return dict(r1=r1, r12=r12, r13=r13, r1213=r1213, a=a, b=b, c=c, mat=mat, nn=nn)
def mesh_edges(tris):
"""Return sparse matrix with edges as an adjacency matrix.
Parameters
----------
tris : array of shape [n_triangles x 3]
The triangles.
Returns
-------
edges : scipy.sparse.sparray
The adjacency matrix.
"""
tris = _hashable_ndarray(tris)
return _mesh_edges(tris=tris)
@lru_cache(maxsize=10)
def _mesh_edges(tris=None):
from scipy.sparse import coo_array
if np.max(tris) > len(np.unique(tris)):
raise ValueError("Cannot compute adjacency on a selection of triangles.")
npoints = np.max(tris) + 1
ones_ntris = np.ones(3 * len(tris))
a, b, c = tris.T
x = np.concatenate((a, b, c))
y = np.concatenate((b, c, a))
edges = coo_array((ones_ntris, (x, y)), shape=(npoints, npoints))
edges = edges.tocsr()
edges = edges + edges.T
return edges
def mesh_dist(tris, vert):
"""Compute adjacency matrix weighted by distances.
It generates an adjacency matrix where the entries are the distances
between neighboring vertices.
Parameters
----------
tris : array (n_tris x 3)
Mesh triangulation.
vert : array (n_vert x 3)
Vertex locations.
Returns
-------
dist_matrix : scipy.sparse.csr_array
Sparse matrix with distances between adjacent vertices.
"""
from scipy.sparse import csr_array
edges = mesh_edges(tris).tocoo()
dist = np.linalg.norm(vert[edges.row, :] - vert[edges.col, :], axis=1)
dist_matrix = csr_array((dist, (edges.row, edges.col)), shape=edges.shape)
return dist_matrix
@verbose_static()
def read_tri(fname_in, swap=False, verbose=None):
"""Read triangle definitions from an ascii file.
Parameters
----------
fname_in : path-like
Path to surface ASCII file (ending with ``'.tri'``).
swap : bool
Assume the ASCII file vertex ordering is clockwise instead of
counterclockwise.
verbose : bool | str | int | None
Control verbosity of the logging output. If ``None``, use the default
verbosity level. See the :ref:`logging documentation <tut-logging>` and
:func:`mne.verbose` for details. Should only be passed as a keyword
argument.
Returns
-------
rr : array, shape=(n_vertices, 3)
Coordinate points.
tris : int array, shape=(n_faces, 3)
Triangulation (each line contains indices for three points which
together form a face).
See Also
--------
read_surface
write_surface
Notes
-----
.. versionadded:: 0.13.0
"""
with open(fname_in) as fid:
lines = fid.readlines()
n_nodes = int(lines[0])
n_tris = int(lines[n_nodes + 1])
n_items = len(lines[1].split())
if n_items in [3, 6, 14, 17]:
inds = range(3)
elif n_items in [4, 7]:
inds = range(1, 4)
else:
raise OSError("Unrecognized format of data.")
rr = np.array(
[
np.array([float(v) for v in line.split()])[inds]
for line in lines[1 : n_nodes + 1]
]
)
tris = np.array(
[
np.array([int(v) for v in line.split()])[inds]
for line in lines[n_nodes + 2 : n_nodes + 2 + n_tris]
]
)
if swap:
tris[:, [2, 1]] = tris[:, [1, 2]]
tris -= 1
logger.info(
f"Loaded surface from {fname_in} with {n_nodes} nodes and {n_tris} triangles."
)
if n_items in [3, 4]:
logger.info("Node normals were not included in the source file.")
else:
warn("Node normals were not read.")
return (rr, tris)
def _get_solids_numpy(tri_rrs, fros):
"""Compute _sum_solids_div total angle in chunks (NumPy fallback)."""
tot_angle = np.zeros(len(fros))
for ti in range(len(tri_rrs)):
tri_rr = tri_rrs[ti]
v1 = fros - tri_rr[0]
v2 = fros - tri_rr[1]
v3 = fros - tri_rr[2]
triple = np.sum(np.cross(v1, v2) * v3, axis=1)
l1 = np.sqrt(np.sum(v1 * v1, axis=1))
l2 = np.sqrt(np.sum(v2 * v2, axis=1))
l3 = np.sqrt(np.sum(v3 * v3, axis=1))
s = (
l1 * l2 * l3
+ np.sum(v1 * v2, axis=1) * l3
+ np.sum(v1 * v3, axis=1) * l2
+ np.sum(v2 * v3, axis=1) * l1
)
tot_angle -= np.arctan2(triple, s)
return tot_angle
def _complete_sphere_surf(sphere, idx, level, complete=True):
"""Convert sphere conductor model to surface."""
rad = sphere["layers"][idx]["rad"]
r0 = sphere["r0"]
surf = _tessellate_sphere_surf(level, rad=rad)
surf["rr"] += r0
if complete:
complete_surface_info(surf, copy=False)
surf["coord_frame"] = sphere["coord_frame"]
return surf
@verbose_static("info_not_none", "dig_kinds", "exclude_frontal", "on_defects")
def dig_mri_distances(
info,
trans,
subject,
subjects_dir=None,
dig_kinds="auto",
exclude_frontal=False,
on_defects="raise",
verbose=None,
):
"""Compute distances between head shape points and the scalp surface.
This function is useful to check that coregistration is correct.
Unless outliers are present in the head shape points,
one can assume an average distance around 2-3 mm.
Parameters
----------
info : mne.Info
The :class:`mne.Info` object with information about the
sensors and methods of measurement.
Must contain the head shape points in ``info['dig']``.
trans : str | instance of Transform
The head<->MRI transform. If str is passed it is the
path to file on disk.
subject : str
The name of the subject.
subjects_dir : str | None
Directory containing subjects data. If None use
the FreeSurfer SUBJECTS_DIR environment variable.
dig_kinds : list of str | str
Kind of digitization points to use in the fitting. These can be any
combination of ('cardinal', 'hpi', 'eeg', 'extra'). Can also
be 'auto' (default), which will use only the 'extra' points if
enough (more than 4) are available, and if not, uses 'extra' and
'eeg' points.
exclude_frontal : bool
If True, exclude points that have both negative Z values
(below the nasion) and positive Y values (in front of the LPA/RPA).
Default is False.
on_defects : 'raise' | 'warn' | 'ignore'
What to do if the surface is found to have topological defects.
Can be ``'raise'`` (default) to raise an error, ``'warn'`` to emit a
warning, or ``'ignore'`` to ignore when one or more defects are found.
Note that a lot of computations in MNE-Python assume the surfaces to be
topologically correct, topological defects may still make other
computations (e.g., `mne.make_bem_model` and `mne.make_bem_solution`)
fail irrespective of this parameter.
.. versionadded:: 1.0
verbose : bool | str | int | None
Control verbosity of the logging output. If ``None``, use the default
verbosity level. See the :ref:`logging documentation <tut-logging>` and
:func:`mne.verbose` for details. Should only be passed as a keyword
argument.
Returns
-------
dists : array, shape (n_points,)
The distances.
See Also
--------
mne.bem.get_fitting_dig
Notes
-----
.. versionadded:: 0.19
"""
from .bem import get_fitting_dig
pts = get_head_surf(
subject,
("head-dense", "head", "bem"),
subjects_dir=subjects_dir,
on_defects=on_defects,
)["rr"]
trans = _get_trans(trans, fro="mri", to="head")[0]
pts = apply_trans(trans, pts)
info_dig = get_fitting_dig(info, dig_kinds, exclude_frontal=exclude_frontal)
dists = _compute_nearest(pts, info_dig, return_dists=True)[1]
return dists
def _mesh_borders(tris, mask):
assert isinstance(mask, np.ndarray) and mask.ndim == 1
edges = mesh_edges(tris)
edges = edges.tocoo()
border_edges = mask[edges.row] != mask[edges.col]
return np.unique(edges.row[border_edges])
def _marching_cubes(image, level, smooth=0, fill_hole_size=None, use_flying_edges=True):
"""Compute marching cubes on a 3D image."""
from scipy.ndimage import binary_dilation
from vtkmodules.util.numpy_support import numpy_to_vtk, vtk_to_numpy
from vtkmodules.vtkCommonDataModel import vtkDataSetAttributes, vtkImageData
from vtkmodules.vtkFiltersCore import vtkThreshold
from vtkmodules.vtkFiltersGeneral import (
vtkDiscreteFlyingEdges3D,
vtkDiscreteMarchingCubes,
)
from vtkmodules.vtkFiltersGeometry import vtkGeometryFilter
if image.ndim != 3:
raise ValueError(f"3D data must be supplied, got {image.shape}")
level = np.array(level)
if level.ndim != 1 or level.size == 0 or level.dtype.kind not in "ui":
raise TypeError(
"level must be non-empty numeric or 1D array-like of int, "
f"got {level.ndim}D array-like of {level.dtype} with "
f"{level.size} elements"
)
image = np.array(image, dtype=float, order="F")
image_shape = image.shape
if fill_hole_size is not None:
for val in level:
bin_image = image == val
mask = image == 0
bin_image = binary_dilation(bin_image, iterations=fill_hole_size, mask=mask)
image[bin_image] = val
data_vtk = numpy_to_vtk(image.ravel(order="F"), deep=False)
mc = vtkDiscreteFlyingEdges3D() if use_flying_edges else vtkDiscreteMarchingCubes()
imdata = vtkImageData()
imdata.SetDimensions(image_shape)
imdata.SetSpacing([1, 1, 1])
imdata.SetOrigin([0, 0, 0])
imdata.GetPointData().SetScalars(data_vtk)
mc.SetNumberOfContours(len(level))
for li, lev in enumerate(level):
mc.SetValue(li, lev)
mc.SetInputData(imdata)
mc.Update()
mc = _vtk_smooth(mc.GetOutput(), smooth)
selector = vtkThreshold()
selector.SetInputData(mc)
dsa = vtkDataSetAttributes()
selector.SetInputArrayToProcess(
0,
0,
0,
imdata.FIELD_ASSOCIATION_POINTS
if use_flying_edges
else imdata.FIELD_ASSOCIATION_CELLS,
dsa.SCALARS,
)
geometry = vtkGeometryFilter()
geometry.SetInputConnection(selector.GetOutputPort())
out = list()
for val in level:
try:
selector.SetLowerThreshold
except AttributeError:
selector.ThresholdBetween(val, val)
else:
selector.SetLowerThreshold(val)
selector.SetUpperThreshold(val)
geometry.Update()
polydata = geometry.GetOutput()
rr = vtk_to_numpy(polydata.GetPoints().GetData())
tris = vtk_to_numpy(polydata.GetPolys().GetConnectivityArray()).reshape(-1, 3)
rr = np.ascontiguousarray(rr)
tris = np.ascontiguousarray(tris)
out.append((rr, tris))
return out
def _keep_largest_component(rr, tris):
"""Keep only the largest connected component of a triangulated mesh."""
from scipy.sparse.csgraph import connected_components
if len(tris) == 0:
return rr, tris
n_comp, labels = connected_components(mesh_edges(tris), directed=False)
if n_comp == 1:
return rr, tris
largest = np.argmax(np.bincount(labels))
keep = labels == largest
new_index = np.full(len(rr), -1, int)
new_index[keep] = np.arange(keep.sum())
tris = new_index[tris[keep[tris].all(axis=1)]]
return rr[keep], tris
@_verbose_control
def _vtk_smooth(pd, smooth, *, verbose=None):
_validate_type(smooth, "numeric", smooth)
smooth = float(smooth)
if not 0 <= smooth < 1:
raise ValueError(
"smoothing factor must be between 0 (inclusive) and "
f"1 (exclusive), got {smooth}"
)
if smooth == 0:
return pd
from vtkmodules.vtkFiltersCore import vtkWindowedSincPolyDataFilter
logger.info(f" Smoothing by a factor of {smooth}")
return_ndarray = False
if isinstance(pd, dict):
pd = _surface_to_polydata(pd)
return_ndarray = True
smoother = vtkWindowedSincPolyDataFilter()
smoother.SetInputData(pd)
smoother.SetNumberOfIterations(100)
smoother.BoundarySmoothingOff()
smoother.FeatureEdgeSmoothingOff()
smoother.SetFeatureAngle(120.0)
smoother.SetPassBand(1 - smooth)
smoother.NonManifoldSmoothingOn()
smoother.NormalizeCoordinatesOff()
smoother.Update()
out = smoother.GetOutput()
if return_ndarray:
out = _polydata_to_surface(out, normals=False)
return out
_VOXELS_MAX = 1000
@fill_doc_static("montage", "subject", "subjects_dir", "aseg")
def get_montage_volume_labels(montage, subject, subjects_dir=None, aseg="auto", dist=2):
"""Get regions of interest near channels from a FreeSurfer parcellation.
.. note:: This is applicable for channels inside the brain
(intracranial electrodes).
Parameters
----------
montage : None | str | DigMontage
A montage containing channel positions. If a string or
:class:`~mne.channels.DigMontage` is
specified, the existing channel information will be updated with the
channel positions from the montage. Valid strings are the names of the
built-in montages that ship with MNE-Python; you can list those via
:func:`mne.channels.get_builtin_montages`.
If ``None`` (default), the channel positions will be removed from the
:class:`~mne.Info`.
subject : str
The FreeSurfer subject name.
subjects_dir : path-like | None
The path to the directory containing the FreeSurfer subjects
reconstructions. If ``None``, defaults to the ``SUBJECTS_DIR`` environment
variable.
aseg : str
The anatomical segmentation file. Default ``auto`` uses ``aparc+aseg``
if available and ``wmparc`` if not. This may be any anatomical
segmentation file in the mri subdirectory of the FreeSurfer subject
directory.
.. versionchanged:: 1.8
Added support for the new default ``'auto'``.
dist : float
The distance in mm to use for identifying regions of interest.
Returns
-------
labels : dict
The regions of interest labels within ``dist`` of each channel.
colors : dict
The FreeSurfer lookup table colors for the labels.
"""
from ._freesurfer import _get_aseg, read_freesurfer_lut
from .channels import DigMontage
_validate_type(montage, DigMontage, "montage")
_validate_type(dist, (int, float), "dist")
if dist < 0 or dist > 10:
raise ValueError("`dist` must be between 0 and 10")
aseg, aseg_data = _get_aseg(aseg, subject, subjects_dir)
lut, fs_colors = read_freesurfer_lut()
label_lut = {v: k for k, v in lut.items()}
assert all([idx in label_lut for idx in np.unique(aseg_data)])
vox2ras_tkr = aseg.header.get_vox2ras_tkr()
ch_dict = montage.get_positions()
if ch_dict["coord_frame"] != "mri":
raise RuntimeError(
"Coordinate frame not supported, expected "
'"mri", got ' + str(ch_dict["coord_frame"])
)
ch_coords = np.array(list(ch_dict["ch_pos"].values()))
ch_coords = apply_trans(
np.linalg.inv(aseg.header.get_vox2ras_tkr()), ch_coords * 1000
)
labels = OrderedDict()
for ch_name, ch_coord in zip(montage.ch_names, ch_coords):
if np.isnan(ch_coord).any():
labels[ch_name] = list()
else:
voxels = _voxel_neighbors(
ch_coord,
aseg_data,
dist=dist,
vox2ras_tkr=vox2ras_tkr,
voxels_max=_VOXELS_MAX,
)
label_idxs = set([aseg_data[tuple(voxel)].astype(int) for voxel in voxels])
labels[ch_name] = [label_lut[idx] for idx in label_idxs]
all_labels = set([label for val in labels.values() for label in val])
colors = {label: tuple(fs_colors[label][:3] / 255) + (1.0,) for label in all_labels}
return labels, colors
def _get_neighbors(loc, image, voxels, thresh, dist_params):
"""Find all the neighbors above a threshold near a voxel."""
neighbors = set()
for axis in range(len(loc)):
for i in (-1, 1):
next_loc = np.array(loc)
next_loc[axis] += i
if thresh is not None:
assert dist_params is None
next_loc = tuple(next_loc)
if (
image[next_loc] > thresh
and image[next_loc] <= image[loc]
and next_loc not in voxels
):
neighbors.add(next_loc)
else:
assert thresh is None
dist, seed_fs_ras, vox2ras_tkr = dist_params
next_loc_fs_ras = apply_trans(vox2ras_tkr, next_loc + 0.5)
if np.linalg.norm(seed_fs_ras - next_loc_fs_ras) <= dist:
neighbors.add(tuple(next_loc))
return neighbors
def _voxel_neighbors(
seed,
image,
thresh=None,
max_peak_dist=1,
use_relative=True,
dist=None,
vox2ras_tkr=None,
voxels_max=100,
):
"""Find voxels above a threshold contiguous with a seed location.
Parameters
----------
seed : tuple | ndarray
The location in image coordinated to seed the algorithm.
image : ndarray
The image to search.
thresh : float
The threshold to use as a cutoff for what qualifies as a neighbor.
Will be relative to the peak if ``use_relative`` or absolute if not.
max_peak_dist : int
The maximum number of voxels to search for the peak near
the seed location.
use_relative : bool
If ``True``, the threshold will be relative to the peak, if
``False``, the threshold will be absolute.
dist : float
The distance in mm to include surrounding voxels.
vox2ras_tkr : ndarray
The voxel to surface RAS affine. Must not be None if ``dist``
if not None.
voxels_max : int
The maximum size of the output ``voxels``.
Returns
-------
voxels : set
The set of locations including the ``seed`` voxel and
surrounding that meet the criteria.
.. note:: Either ``dist`` or ``thesh`` may be used but not both.
When ``thresh`` is used, first a peak nearby the seed
location is found and then voxels are only included if they
decrease monotonically from the peak. When ``dist`` is used,
only voxels within ``dist`` mm of the seed are included.
"""
seed = np.array(seed).round().astype(int)
assert ((dist is not None) + (thresh is not None)) == 1
if thresh is not None:
dist_params = None
check_grid = image[
tuple([slice(idx - max_peak_dist, idx + max_peak_dist + 1) for idx in seed])
]
peak = (
np.array(np.unravel_index(np.argmax(check_grid), check_grid.shape))
- max_peak_dist
+ seed
)
voxels = neighbors = set([tuple(peak)])
if use_relative:
thresh *= image[tuple(peak)]
else:
assert vox2ras_tkr is not None
seed_fs_ras = apply_trans(vox2ras_tkr, seed + 0.5)
dist_params = (dist, seed_fs_ras, vox2ras_tkr)
voxels = neighbors = set([tuple(seed)])
while neighbors and len(voxels) <= voxels_max:
next_neighbors = set()
for next_loc in neighbors:
voxel_neighbors = _get_neighbors(
next_loc, image, voxels, thresh, dist_params
)
voxel_neighbors = voxel_neighbors.difference(voxels)
next_neighbors = next_neighbors.union(voxel_neighbors)
voxels = voxels.union(voxel_neighbors)
if len(voxels) > voxels_max:
break
neighbors = next_neighbors
return voxels
def __getattr__(name):
if name not in (
"_accumulate_normals",
"_fast_cross_nd_sum",
"_find_nearest_tri_pts",
"_get_solids",
"_get_tri_dist",
"_nearest_tri_edge",
"_triangle_coords",
):
raise AttributeError(f"module {__name__!r} has no attribute {name!r}")
import importlib
return getattr(importlib.import_module("._surface_numba", __package__), name)