Aluode/PerceptionLabPortable
0
1# Authors: The MNE-Python contributors.2# License: BSD-3-Clause3# Copyright the MNE-Python contributors.4 5import copy6import os.path as op7import warnings8 9import numpy as np10from scipy import sparse11 12from .fixes import _eye_array, _get_img_fdata13from .morph_map import read_morph_map14from .parallel import parallel_func15from .source_estimate import (16 _BaseSourceEstimate,17 _BaseSurfaceSourceEstimate,18 _BaseVolSourceEstimate,19 _get_ico_tris,20)21from .source_space._source_space import SourceSpaces, _ensure_src, _grid_interp22from .surface import _compute_nearest, mesh_edges, read_surface23from .utils import (24 BunchConst,25 ProgressBar,26 _check_fname,27 _check_option,28 _custom_lru_cache,29 _ensure_int,30 _import_h5io_funcs,31 _import_nibabel,32 _validate_type,33 check_version,34 fill_doc,35 get_subjects_dir,36 logger,37 use_log_level,38 verbose,39 warn,40)41from .utils import (42 warn as warn_,43)44 45 46@verbose47def compute_source_morph(48 src,49 subject_from=None,50 subject_to="fsaverage",51 subjects_dir=None,52 zooms="auto",53 niter_affine=(100, 100, 10),54 niter_sdr=(5, 5, 3),55 spacing=5,56 smooth=None,57 warn=True,58 xhemi=False,59 sparse=False,60 src_to=None,61 precompute=False,62 verbose=None,63):64 """Create a SourceMorph from one subject to another.65 66 Method is based on spherical morphing by FreeSurfer for surface67 cortical estimates :footcite:`GreveEtAl2013` and68 Symmetric Diffeomorphic Registration for volumic data69 :footcite:`AvantsEtAl2008`.70 71 Parameters72 ----------73 src : instance of SourceSpaces | instance of SourceEstimate74 The SourceSpaces of subject_from (can be a75 SourceEstimate if only using a surface source space).76 subject_from : str | None77 Name of the original subject as named in the SUBJECTS_DIR.78 If None (default), then ``src[0]['subject_his_id]'`` will be used.79 subject_to : str | None80 Name of the subject to which to morph as named in the SUBJECTS_DIR.81 Default is ``'fsaverage'``. If None, ``src_to[0]['subject_his_id']``82 will be used.83 84 .. versionchanged:: 0.2085 Support for subject_to=None.86 %(subjects_dir)s87 zooms : float | tuple | str | None88 The voxel size of volume for each spatial dimension in mm.89 If spacing is None, MRIs won't be resliced, and both volumes90 must have the same number of spatial dimensions.91 Can also be ``'auto'`` to use ``5.`` if ``src_to is None`` and92 the zooms from ``src_to`` otherwise.93 94 .. versionchanged:: 0.2095 Support for 'auto' mode.96 niter_affine : tuple of int97 Number of levels (``len(niter_affine)``) and number of98 iterations per level - for each successive stage of iterative99 refinement - to perform the affine transform.100 Default is niter_affine=(100, 100, 10).101 niter_sdr : tuple of int102 Number of levels (``len(niter_sdr)``) and number of103 iterations per level - for each successive stage of iterative104 refinement - to perform the Symmetric Diffeomorphic Registration (sdr)105 transform. Default is niter_sdr=(5, 5, 3).106 spacing : int | list | None107 The resolution of the icosahedral mesh (typically 5).108 If None, all vertices will be used (potentially filling the109 surface). If a list, then values will be morphed to the set of110 vertices specified in in ``spacing[0]`` and ``spacing[1]``.111 This will be ignored if ``src_to`` is supplied.112 113 .. versionchanged:: 0.21114 src_to, if provided, takes precedence.115 smooth : int | str | None116 Number of iterations for the smoothing of the surface data.117 If None, smooth is automatically defined to fill the surface118 with non-zero values. Can also be ``'nearest'`` to use the nearest119 vertices on the surface.120 121 .. versionchanged:: 0.20122 Added support for 'nearest'.123 warn : bool124 If True, warn if not all vertices were used. The default is warn=True.125 xhemi : bool126 Morph across hemisphere. Currently only implemented for127 ``subject_to == subject_from``. See notes below.128 The default is xhemi=False.129 sparse : bool130 Morph as a sparse source estimate. Works only with (Vector)131 SourceEstimate. If True the only parameters used are subject_to and132 subject_from, and spacing has to be None. Default is sparse=False.133 src_to : instance of SourceSpaces | None134 The destination source space.135 136 - For surface-based morphing, this is the preferred over ``spacing``137 for providing the vertices.138 - For volumetric morphing, this should be passed so that 1) the139 resultingmorph volume is properly constrained to the brain volume,140 and 2) STCs from multiple subjects morphed to the same destination141 subject/source space have the vertices.142 - For mixed (surface + volume) morphing, this is required.143 144 .. versionadded:: 0.20145 precompute : bool146 If True (default False), compute the sparse matrix representation of147 the volumetric morph (if present). This takes a long time to148 compute, but can make morphs faster when thousands of points are used.149 See :meth:`mne.SourceMorph.compute_vol_morph_mat` (which can be called150 later if desired) for more information.151 152 .. versionadded:: 0.22153 %(verbose)s154 155 Returns156 -------157 morph : instance of SourceMorph158 The :class:`mne.SourceMorph` object.159 160 Notes161 -----162 This function can be used to morph surface data between hemispheres by163 setting ``xhemi=True``. The full cross-hemisphere morph matrix maps left164 to right and right to left. A matrix for cross-mapping only one hemisphere165 can be constructed by specifying the appropriate vertices, for example, to166 map the right hemisphere to the left::167 168 vertices_from=[[], vert_rh], vertices_to=[vert_lh, []]169 170 Cross-hemisphere mapping requires appropriate ``sphere.left_right``171 morph-maps in the subject's directory. These morph maps are included172 with the ``fsaverage_sym`` FreeSurfer subject, and can be created for other173 subjects with the ``mris_left_right_register`` FreeSurfer command. The174 ``fsaverage_sym`` subject is included with FreeSurfer > 5.1 and can be175 obtained as described `here176 <https://surfer.nmr.mgh.harvard.edu/fswiki/Xhemi>`_. For statistical177 comparisons between hemispheres, use of the symmetric ``fsaverage_sym``178 model is recommended to minimize bias :footcite:`GreveEtAl2013`.179 180 .. versionadded:: 0.17.0181 182 .. versionadded:: 0.21.0183 Support for morphing mixed source estimates.184 185 References186 ----------187 .. footbibliography::188 """189 src_data, kind, src_subject = _get_src_data(src)190 subject_from = _check_subject_src(subject_from, src_subject, warn_none=True)191 del src192 _validate_type(src_to, (SourceSpaces, None), "src_to")193 _validate_type(subject_to, (str, None), "subject_to")194 if src_to is None and subject_to is None:195 raise ValueError("subject_to cannot be None when src_to is None")196 subject_to = _check_subject_src(subject_to, src_to, "subject_to")197 198 # Params199 warn = False if sparse else warn200 201 if kind not in "surface" and xhemi:202 raise ValueError(203 "Inter-hemispheric morphing can only be used with surface source estimates."204 )205 if sparse and kind != "surface":206 raise ValueError("Only surface source estimates can compute a sparse morph.")207 208 subjects_dir = str(get_subjects_dir(subjects_dir, raise_error=True))209 shape = affine = pre_affine = sdr_morph = morph_mat = None210 vertices_to_surf, vertices_to_vol = list(), list()211 212 if kind in ("volume", "mixed"):213 _check_dep(nibabel="2.1.0", dipy="0.10.1")214 nib = _import_nibabel("work with a volume source space")215 216 logger.info("Volume source space(s) present...")217 218 # load moving MRI219 mri_subpath = op.join("mri", "brain.mgz")220 mri_path_from = op.join(subjects_dir, subject_from, mri_subpath)221 222 logger.info(f' Loading {mri_path_from} as "from" volume')223 with warnings.catch_warnings():224 mri_from = nib.load(mri_path_from)225 226 # eventually we could let this be some other volume, but for now227 # let's KISS and use `brain.mgz`, too228 mri_path_to = op.join(subjects_dir, subject_to, mri_subpath)229 if not op.isfile(mri_path_to):230 raise OSError(f"cannot read file: {mri_path_to}")231 logger.info(f' Loading {mri_path_to} as "to" volume')232 with warnings.catch_warnings():233 mri_to = nib.load(mri_path_to)234 235 # deal with `src_to` subsampling236 zooms_src_to = None237 if src_to is None:238 if kind == "mixed":239 raise ValueError(240 "src_to must be provided when using a mixed source space"241 )242 else:243 surf_offset = 2 if src_to.kind == "mixed" else 0244 # All of our computations are in RAS (like img.affine), so we need245 # to get the transformation from RAS to the source space246 # subsampling of vox (src), not MRI (FreeSurfer surface RAS) to src247 src_ras_t = np.dot(248 src_to[-1]["mri_ras_t"]["trans"], src_to[-1]["src_mri_t"]["trans"]249 )250 src_ras_t[:3] *= 1e3251 src_data["to_vox_map"] = (src_to[-1]["shape"], src_ras_t)252 vertices_to_vol = [s["vertno"] for s in src_to[surf_offset:]]253 zooms_src_to = np.diag(src_to[-1]["src_mri_t"]["trans"])[:3] * 1000254 zooms_src_to = tuple(zooms_src_to)255 256 # pre-compute non-linear morph257 zooms = _check_zooms(mri_from, zooms, zooms_src_to)258 shape, zooms, affine, pre_affine, sdr_morph = _compute_morph_sdr(259 mri_from, mri_to, niter_affine, niter_sdr, zooms260 )261 262 if kind in ("surface", "mixed"):263 logger.info("surface source space present ...")264 vertices_from = src_data["vertices_from"]265 if sparse:266 if spacing is not None:267 raise ValueError("spacing must be set to None if sparse=True.")268 if xhemi:269 raise ValueError("xhemi=True can only be used with sparse=False")270 vertices_to_surf, morph_mat = _compute_sparse_morph(271 vertices_from, subject_from, subject_to, subjects_dir272 )273 else:274 if src_to is not None:275 assert src_to.kind in ("surface", "mixed")276 vertices_to_surf = [s["vertno"].copy() for s in src_to[:2]]277 else:278 vertices_to_surf = grade_to_vertices(279 subject_to, spacing, subjects_dir, 1280 )281 morph_mat = _compute_morph_matrix(282 subject_from=subject_from,283 subject_to=subject_to,284 vertices_from=vertices_from,285 vertices_to=vertices_to_surf,286 subjects_dir=subjects_dir,287 smooth=smooth,288 warn=warn,289 xhemi=xhemi,290 )291 n_verts = sum(len(v) for v in vertices_to_surf)292 assert morph_mat.shape[0] == n_verts293 294 vertices_to = vertices_to_surf + vertices_to_vol295 if src_to is not None:296 assert len(vertices_to) == len(src_to)297 morph = SourceMorph(298 subject_from,299 subject_to,300 kind,301 zooms,302 niter_affine,303 niter_sdr,304 spacing,305 smooth,306 xhemi,307 morph_mat,308 vertices_to,309 shape,310 affine,311 pre_affine,312 sdr_morph,313 src_data,314 None,315 )316 if precompute:317 morph.compute_vol_morph_mat()318 logger.info("[done]")319 return morph320 321 322def _compute_sparse_morph(vertices_from, subject_from, subject_to, subjects_dir=None):323 """Get nearest vertices from one subject to another."""324 from scipy import sparse325 326 maps = read_morph_map(subject_to, subject_from, subjects_dir)327 cnt = 0328 vertices = list()329 cols = list()330 for verts, map_hemi in zip(vertices_from, maps):331 vertno_h = _sparse_argmax_nnz_row(map_hemi[verts])332 order = np.argsort(vertno_h)333 cols.append(cnt + order)334 vertices.append(vertno_h[order])335 cnt += len(vertno_h)336 cols = np.concatenate(cols)337 rows = np.arange(len(cols))338 data = np.ones(len(cols))339 morph_mat = sparse.coo_array(340 (data, (rows, cols)), shape=(len(cols), len(cols))341 ).tocsr()342 return vertices, morph_mat343 344 345_SOURCE_MORPH_ATTRIBUTES = [ # used in writing346 "subject_from",347 "subject_to",348 "kind",349 "zooms",350 "niter_affine",351 "niter_sdr",352 "spacing",353 "smooth",354 "xhemi",355 "morph_mat",356 "vertices_to",357 "shape",358 "affine",359 "pre_affine",360 "sdr_morph",361 "src_data",362 "vol_morph_mat",363]364 365 366@fill_doc367class SourceMorph:368 """Morph source space data from one subject to another.369 370 .. note::371 This class should not be instantiated directly via372 ``mne.SourceMorph(...)``. Instead, use one of the functions373 listed in the See Also section below.374 375 Parameters376 ----------377 subject_from : str | None378 Name of the subject from which to morph as named in the SUBJECTS_DIR.379 subject_to : str | array | list of array380 Name of the subject on which to morph as named in the SUBJECTS_DIR.381 The default is 'fsaverage'. If morphing a volume source space,382 subject_to can be the path to a MRI volume. Can also be a list of383 two arrays if morphing to hemisphere surfaces.384 kind : str | None385 Kind of source estimate. E.g. ``'volume'`` or ``'surface'``.386 zooms : float | tuple387 See :func:`mne.compute_source_morph`.388 niter_affine : tuple of int389 Number of levels (``len(niter_affine)``) and number of390 iterations per level - for each successive stage of iterative391 refinement - to perform the affine transform.392 niter_sdr : tuple of int393 Number of levels (``len(niter_sdr)``) and number of394 iterations per level - for each successive stage of iterative395 refinement - to perform the Symmetric Diffeomorphic Registration (sdr)396 transform :footcite:`AvantsEtAl2008`.397 spacing : int | list | None398 See :func:`mne.compute_source_morph`.399 smooth : int | str | None400 See :func:`mne.compute_source_morph`.401 xhemi : bool402 Morph across hemisphere.403 morph_mat : scipy.sparse.csr_array404 The sparse surface morphing matrix for spherical surface405 based morphing :footcite:`GreveEtAl2013`.406 vertices_to : list of ndarray407 The destination surface vertices.408 shape : tuple409 The volume MRI shape.410 affine : ndarray411 The volume MRI affine.412 pre_affine : instance of dipy.align.AffineMap413 The transformation that is applied before the before ``sdr_morph``.414 sdr_morph : instance of dipy.align.DiffeomorphicMap415 The class that applies the the symmetric diffeomorphic registration416 (SDR) morph.417 src_data : dict418 Additional source data necessary to perform morphing.419 vol_morph_mat : scipy.sparse.csr_array | None420 The volumetric morph matrix, if :meth:`compute_vol_morph_mat`421 was used.422 %(verbose)s423 424 See Also425 --------426 compute_source_morph427 read_source_morph428 429 Notes430 -----431 .. versionadded:: 0.17432 433 References434 ----------435 .. footbibliography::436 """437 438 @verbose439 def __init__(440 self,441 subject_from,442 subject_to,443 kind,444 zooms,445 niter_affine,446 niter_sdr,447 spacing,448 smooth,449 xhemi,450 morph_mat,451 vertices_to,452 shape,453 affine,454 pre_affine,455 sdr_morph,456 src_data,457 vol_morph_mat,458 *,459 verbose=None,460 ):461 # universal462 self.subject_from = subject_from463 self.subject_to = subject_to464 self.kind = kind465 # vol input466 self.zooms = zooms467 self.niter_affine = niter_affine468 self.niter_sdr = niter_sdr469 # surf input470 self.spacing = spacing471 self.smooth = smooth472 self.xhemi = xhemi473 # surf computed474 self.morph_mat = morph_mat475 # vol computed476 self.shape = shape477 self.affine = affine478 self.sdr_morph = sdr_morph479 self.pre_affine = pre_affine480 # used by both481 self.src_data = src_data482 self.vol_morph_mat = vol_morph_mat483 # compute vertices_to here (partly for backward compat and no src484 # provided)485 if vertices_to is None or len(vertices_to) == 0 and kind == "volume":486 assert src_data["to_vox_map"] is None487 vertices_to = self._get_vol_vertices_to_nz()488 self.vertices_to = vertices_to489 490 @property491 def _vol_vertices_from(self):492 assert isinstance(self.src_data["inuse"], list)493 vertices_from = [np.where(in_)[0] for in_ in self.src_data["inuse"]]494 return vertices_from495 496 @property497 def _vol_vertices_to(self):498 return self.vertices_to[0 if self.kind == "volume" else 2 :]499 500 def _get_vol_vertices_to_nz(self):501 logger.info("Computing nonzero vertices after morph ...")502 n_vertices = sum(len(v) for v in self._vol_vertices_from)503 ones = np.ones((n_vertices, 1))504 with use_log_level(False):505 return [np.where(self._morph_vols(ones, "", subselect=False))[0]]506 507 @verbose508 def apply(509 self, stc_from, output="stc", mri_resolution=False, mri_space=None, verbose=None510 ):511 """Morph source space data.512 513 Parameters514 ----------515 stc_from : VolSourceEstimate | VolVectorSourceEstimate | SourceEstimate | VectorSourceEstimate516 The source estimate to morph.517 output : str518 Can be ``'stc'`` (default) or possibly ``'nifti1'``, or519 ``'nifti2'`` when working with a volume source space defined on a520 regular grid.521 mri_resolution : bool | tuple | int | float522 If True the image is saved in MRI resolution. Default False.523 524 .. warning: If you have many time points the file produced can be525 huge. The default is ``mri_resolution=False``.526 mri_space : bool | None527 Whether the image to world registration should be in mri space. The528 default (None) is mri_space=mri_resolution.529 %(verbose)s530 531 Returns532 -------533 stc_to : VolSourceEstimate | SourceEstimate | VectorSourceEstimate | Nifti1Image | Nifti2Image534 The morphed source estimates.535 """ # noqa: E501536 _validate_type(output, str, "output")537 _validate_type(stc_from, _BaseSourceEstimate, "stc_from", "source estimate")538 if isinstance(stc_from, _BaseSurfaceSourceEstimate):539 allowed_kinds = ("stc",)540 extra = "when stc is a surface source estimate"541 else:542 allowed_kinds = ("stc", "nifti1", "nifti2")543 extra = ""544 _check_option("output", output, allowed_kinds, extra)545 stc = copy.deepcopy(stc_from)546 547 mri_space = mri_resolution if mri_space is None else mri_space548 if stc.subject is None:549 stc.subject = self.subject_from550 if self.subject_from is None:551 self.subject_from = stc.subject552 if stc.subject != self.subject_from:553 raise ValueError(554 "stc_from.subject and "555 "morph.subject_from "556 f"must match. ({stc.subject} != {self.subject_from})"557 )558 out = _apply_morph_data(self, stc)559 if output != "stc": # convert to volume560 out = _morphed_stc_as_volume(561 self,562 out,563 mri_resolution=mri_resolution,564 mri_space=mri_space,565 output=output,566 )567 return out568 569 @verbose570 def compute_vol_morph_mat(self, *, verbose=None):571 """Compute the sparse matrix representation of the volumetric morph.572 573 Parameters574 ----------575 %(verbose)s576 577 Returns578 -------579 morph : instance of SourceMorph580 The instance (modified in-place).581 582 Notes583 -----584 For a volumetric morph, this will compute the morph for an identity585 source volume, i.e., with one source vertex active at a time, and store586 the result as a :class:`sparse <scipy.sparse.csr_array>`587 morphing matrix. This takes a long time (minutes) to compute initially,588 but drastically speeds up :meth:`apply` for STCs, so it can be589 beneficial when many time points or many morphs (i.e., greater than590 the number of volumetric ``src_from`` vertices) will be performed.591 592 When calling :meth:`save`, this sparse morphing matrix is saved with593 the instance, so this only needs to be called once. This function does594 nothing if the morph matrix has already been computed, or if there is595 no volume morphing necessary.596 597 .. versionadded:: 0.22598 """599 if self.affine is None or self.vol_morph_mat is not None:600 return601 logger.info("Computing sparse volumetric morph matrix (will take some time...)")602 self.vol_morph_mat = self._morph_vols(None, "Vertex")603 return self604 605 def _morph_vols(self, vols, mesg, subselect=True):606 from dipy.align.reslice import reslice607 608 interp = self.src_data["interpolator"].tocsc()[609 :, np.concatenate(self._vol_vertices_from)610 ]611 n_vols = interp.shape[1] if vols is None else vols.shape[1]612 attrs = ("real", "imag") if np.iscomplexobj(vols) else ("real",)613 dtype = np.complex128 if len(attrs) == 2 else np.float64614 if vols is None: # sparse -> sparse mode615 img_to = (list(), list(), [0]) # data, indices, indptr616 assert subselect617 else: # dense -> dense mode618 img_to = None619 if subselect:620 vol_verts = np.concatenate(self._vol_vertices_to)621 else:622 vol_verts = slice(None)623 # morph data624 from_affine = np.dot(625 self.src_data["src_affine_ras"], # mri_ras_t626 self.src_data["src_affine_vox"],627 ) # vox_mri_t628 from_affine[:3] *= 1000.0629 # equivalent of:630 # _resample_from_to(img_real, from_affine,631 # (self.pre_affine.codomain_shape,632 # (self.pre_affine.codomain_grid2world))633 src_shape = self.src_data["src_shape_full"][::-1]634 resamp_0 = _grid_interp(635 src_shape,636 self.pre_affine.codomain_shape,637 np.linalg.inv(from_affine) @ self.pre_affine.codomain_grid2world,638 )639 # reslice to match what was used during the morph640 # (brain.mgz and whatever was used to create the source space641 # will not necessarily have the same domain/zooms)642 # equivalent of:643 # pre_affine.transform(img_real)644 resamp_1 = _grid_interp(645 self.pre_affine.codomain_shape,646 self.pre_affine.domain_shape,647 np.linalg.inv(self.pre_affine.codomain_grid2world)648 @ self.pre_affine.affine649 @ self.pre_affine.domain_grid2world,650 )651 resamp_0_1 = resamp_1 @ resamp_0652 resamp_2 = None653 for ii in ProgressBar(list(range(n_vols)), mesg=mesg):654 for attr in attrs:655 # transform from source space to mri_from resolution/space656 if vols is None:657 img_real = interp[:, [ii]]658 else:659 img_real = interp @ getattr(vols[:, ii], attr)660 _debug_img(img_real, from_affine, "From", src_shape)661 662 img_real = resamp_0_1 @ img_real663 if sparse.issparse(img_real):664 img_real = img_real.toarray()665 img_real = img_real.reshape(self.pre_affine.domain_shape, order="F")666 if self.sdr_morph is not None:667 img_real = self.sdr_morph.transform(img_real)668 _debug_img(img_real, self.affine, "From-reslice-transform")669 670 # subselect the correct cube if src_to is provided671 if self.src_data["to_vox_map"] is not None:672 affine = self.affine673 to_zooms = np.diag(self.src_data["to_vox_map"][1])[:3]674 # There might be some sparse equivalent to this but675 # not sure...676 if not np.allclose(self.zooms, to_zooms, atol=1e-3):677 img_real, affine = reslice(678 img_real, self.affine, self.zooms, to_zooms679 )680 _debug_img(img_real, affine, "From-reslice-transform-src")681 if resamp_2 is None:682 resamp_2 = _grid_interp(683 img_real.shape,684 self.src_data["to_vox_map"][0],685 np.linalg.inv(affine) @ self.src_data["to_vox_map"][1],686 )687 # Equivalent to:688 # _resample_from_to(689 # img_real, affine, self.src_data['to_vox_map'])690 img_real = resamp_2 @ img_real.ravel(order="F")691 _debug_img(692 img_real,693 self.src_data["to_vox_map"][1],694 "From-reslice-transform-src-subselect",695 self.src_data["to_vox_map"][0],696 )697 698 # This can be used to help debug, but it really should just699 # show the brain filling the volume:700 # img_want = np.zeros(np.prod(img_real.shape))701 # img_want[np.concatenate(self._vol_vertices_to)] = 1.702 # img_want = np.reshape(703 # img_want, self.src_data['src_shape'][::-1], order='F')704 # _debug_img(img_want, self.src_data['to_vox_map'][1],705 # 'To mask')706 # raise RuntimeError('Check')707 708 # combine real and complex parts709 img_real = img_real.ravel(order="F")[vol_verts]710 711 # initialize output712 if img_to is None and vols is not None:713 img_to = np.zeros((img_real.size, n_vols), dtype=dtype)714 715 if vols is None:716 idx = np.where(img_real)[0]717 img_to[0].extend(img_real[idx])718 img_to[1].extend(idx)719 img_to[2].append(img_to[2][-1] + len(idx))720 else:721 if attr == "real":722 img_to[:, ii] = img_to[:, ii] + img_real723 else:724 img_to[:, ii] = img_to[:, ii] + 1j * img_real725 726 if vols is None:727 img_to = sparse.csc_array(img_to, shape=(len(vol_verts), n_vols)).tocsr()728 729 return img_to730 731 def __repr__(self): # noqa: D105732 s = f"{self.kind}"733 s += f", {self.subject_from} -> {self.subject_to}"734 if self.kind == "volume":735 s += f", zooms : {self.zooms}"736 s += f", niter_affine : {self.niter_affine}"737 s += f", niter_sdr : {self.niter_sdr}"738 elif self.kind in ("surface", "vector"):739 s += f", spacing : {self.spacing}"740 s += f", smooth : {self.smooth}"741 s += ", xhemi" if self.xhemi else ""742 743 return f"<SourceMorph | {s}>"744 745 @verbose746 def save(self, fname, overwrite=False, verbose=None):747 """Save the morph for source estimates to a file.748 749 Parameters750 ----------751 fname : path-like752 The path to the file. ``'-morph.h5'`` will be added if fname does753 not end with ``'.h5'``.754 %(overwrite)s755 %(verbose)s756 """757 _, write_hdf5 = _import_h5io_funcs()758 fname = _check_fname(fname, overwrite=overwrite, must_exist=False)759 if fname.suffix != ".h5":760 fname = fname.with_name(f"{fname.name}-morph.h5")761 762 out_dict = {k: getattr(self, k) for k in _SOURCE_MORPH_ATTRIBUTES}763 for key in ("pre_affine", "sdr_morph"): # classes764 if out_dict[key] is not None:765 out_dict[key] = out_dict[key].__dict__766 write_hdf5(fname, out_dict, overwrite=overwrite)767 768 769_slicers = list()770 771 772def _debug_img(data, affine, title, shape=None):773 # Uncomment these lines for debugging help with volume morph:774 #775 # import nibabel as nib776 # if sparse.issparse(data):777 # data = data.toarray()778 # data = np.asarray(data)779 # if shape is not None:780 # data = np.reshape(data, shape, order='F')781 # _slicers.append(nib.viewers.OrthoSlicer3D(782 # data, affine, axes=None, title=title))783 # _slicers[-1].figs[0].suptitle(title, color='r')784 return785 786 787def _check_zooms(mri_from, zooms, zooms_src_to):788 # use voxel size of mri_from789 if isinstance(zooms, str) and zooms == "auto":790 zooms = zooms_src_to if zooms_src_to is not None else 5.0791 if zooms is None:792 zooms = mri_from.header.get_zooms()[:3]793 zooms = np.atleast_1d(zooms).astype(float)794 if zooms.shape == (1,):795 zooms = np.repeat(zooms, 3)796 if zooms.shape != (3,):797 raise ValueError(798 "zooms must be None, a singleton, or have shape (3,),"799 f" got shape {zooms.shape}"800 )801 zooms = tuple(zooms)802 return zooms803 804 805# def _resample_from_to(img, affine, to_vox_map):806# # Wrap to dipy for speed, equivalent to:807# # from nibabel.processing import resample_from_to808# # from nibabel.spatialimages import SpatialImage809# # return _get_img_fdata(810# # resample_from_to(SpatialImage(img, affine), to_vox_map, order=1))811# import dipy.align.imaffine812#813# return dipy.align.imaffine.AffineMap(814# None, to_vox_map[0], to_vox_map[1], img.shape, affine815# ).transform(img, resample_only=True)816 817 818###############################################################################819# I/O820def _check_subject_src(821 subject, src, name="subject_from", src_name="src", *, warn_none=False822):823 if isinstance(src, str):824 subject_check = src825 elif src is None: # assume it's correct although dangerous but unlikely826 subject_check = subject827 else:828 subject_check = src._subject829 warn_none = True830 if subject_check is None and warn_none:831 warn(832 "The source space does not contain the subject name, we "833 "recommend regenerating the source space (and forward / "834 "inverse if applicable) for better code reliability"835 )836 if subject is None:837 subject = subject_check838 elif subject_check is not None and subject != subject_check:839 raise ValueError(840 f"{name} does not match {src_name} subject ({subject} != {subject_check})"841 )842 if subject is None:843 raise ValueError(844 f"{name} could not be inferred from {src_name}, it must be specified"845 )846 return subject847 848 849def read_source_morph(fname):850 """Load the morph for source estimates from a file.851 852 Parameters853 ----------854 fname : path-like855 Path to the file containing the morph source estimates.856 857 Returns858 -------859 source_morph : instance of SourceMorph860 The loaded morph.861 """862 read_hdf5, _ = _import_h5io_funcs()863 vals = read_hdf5(fname)864 if vals["pre_affine"] is not None: # reconstruct865 from dipy.align.imaffine import AffineMap866 867 affine = vals["pre_affine"]868 vals["pre_affine"] = AffineMap(None)869 vals["pre_affine"].__dict__ = affine870 if vals["sdr_morph"] is not None:871 from dipy.align.imwarp import DiffeomorphicMap872 873 morph = vals["sdr_morph"]874 vals["sdr_morph"] = DiffeomorphicMap(None, [])875 vals["sdr_morph"].__dict__ = morph876 # Backward compat with when it used to be a list877 if isinstance(vals["vertices_to"], np.ndarray):878 vals["vertices_to"] = [vals["vertices_to"]]879 # Backward compat with when it used to be a single array880 if isinstance(vals["src_data"].get("inuse", None), np.ndarray):881 vals["src_data"]["inuse"] = [vals["src_data"]["inuse"]]882 # added with compute_vol_morph_mat in 0.22:883 vals["vol_morph_mat"] = vals.get("vol_morph_mat", None)884 return SourceMorph(**vals)885 886 887###############################################################################888# Helper functions for SourceMorph methods889def _check_dep(nibabel="2.1.0", dipy="0.10.1"):890 """Check dependencies."""891 for lib, ver in zip(["nibabel", "dipy"], [nibabel, dipy]):892 passed = True if not ver else check_version(lib, ver)893 894 if not passed:895 raise ImportError(896 f"{lib} {ver} or higher must be correctly "897 "installed and accessible from Python"898 )899 900 901def _morphed_stc_as_volume(morph, stc, mri_resolution, mri_space, output):902 """Return volume source space as Nifti1Image and/or save to disk."""903 assert isinstance(stc, _BaseVolSourceEstimate) # should be guaranteed904 if stc._data_ndim == 3:905 stc = stc.magnitude()906 _check_dep(nibabel="2.1.0", dipy=False)907 908 NiftiImage, NiftiHeader = _triage_output(output)909 910 # if MRI resolution is set manually as a single value, convert to tuple911 if isinstance(mri_resolution, int | float):912 # use iso voxel size913 new_zooms = (float(mri_resolution),) * 3914 elif isinstance(mri_resolution, tuple):915 new_zooms = mri_resolution916 # if full MRI resolution, compute zooms from shape and MRI zooms917 if isinstance(mri_resolution, bool):918 new_zooms = _get_zooms_orig(morph) if mri_resolution else None919 920 # create header921 hdr = NiftiHeader()922 hdr.set_xyzt_units("mm", "msec")923 hdr["pixdim"][4] = 1e3 * stc.tstep924 925 # setup empty volume926 if morph.src_data["to_vox_map"] is not None:927 shape = morph.src_data["to_vox_map"][0]928 affine = morph.src_data["to_vox_map"][1]929 else:930 shape = morph.shape931 affine = morph.affine932 assert stc.data.ndim == 2933 n_times = stc.data.shape[1]934 img = np.zeros((np.prod(shape), n_times))935 img[stc.vertices[0], :] = stc.data936 img = img.reshape(shape + (n_times,), order="F") # match order='F' above937 del shape938 939 # make nifti from data940 with warnings.catch_warnings(): # nibabel<->numpy warning941 img = NiftiImage(img, affine, header=hdr)942 943 # reslice in case of manually defined voxel size944 zooms = morph.zooms[:3]945 if new_zooms is not None:946 from dipy.align.reslice import reslice947 948 new_zooms = new_zooms[:3]949 img, affine = reslice(950 _get_img_fdata(img),951 img.affine, # MRI to world registration952 zooms, # old voxel size in mm953 new_zooms,954 ) # new voxel size in mm955 with warnings.catch_warnings(): # nibabel<->numpy warning956 img = NiftiImage(img, affine)957 zooms = new_zooms958 959 # set zooms in header960 img.header.set_zooms(tuple(zooms) + (1,))961 return img962 963 964def _get_src_data(src, mri_resolution=True):965 # copy data to avoid conflicts966 _validate_type(967 src,968 (_BaseSurfaceSourceEstimate, "path-like", SourceSpaces),969 "src",970 "source space or surface source estimate",971 )972 if isinstance(src, _BaseSurfaceSourceEstimate):973 src_t = [dict(vertno=src.vertices[0]), dict(vertno=src.vertices[1])]974 src_kind = "surface"975 src_subject = src.subject976 else:977 src_t = _ensure_src(src).copy()978 src_kind = src_t.kind979 src_subject = src_t._subject980 del src981 _check_option("src kind", src_kind, ("surface", "volume", "mixed"))982 983 # extract all relevant data for volume operations984 src_data = dict()985 if src_kind in ("volume", "mixed"):986 use_src = src_t[-1]987 shape = use_src["shape"]988 start = 0 if src_kind == "volume" else 2989 for si, s in enumerate(src_t[start:], start):990 if s.get("interpolator", None) is None:991 if mri_resolution:992 raise RuntimeError(993 f"MRI interpolator not present in src[{si}], "994 "cannot use mri_resolution=True"995 )996 interpolator = None997 break998 else:999 interpolator = sum((s["interpolator"] for s in src_t[start:]), 0.0)1000 inuses = [s["inuse"] for s in src_t[start:]]1001 src_data.update(1002 {1003 "src_shape": (shape[2], shape[1], shape[0]), # SAR1004 "src_affine_vox": use_src["vox_mri_t"]["trans"],1005 "src_affine_src": use_src["src_mri_t"]["trans"],1006 "src_affine_ras": use_src["mri_ras_t"]["trans"],1007 "src_shape_full": ( # SAR1008 use_src["mri_height"],1009 use_src["mri_depth"],1010 use_src["mri_width"],1011 ),1012 "interpolator": interpolator,1013 "inuse": inuses,1014 "to_vox_map": None,1015 }1016 )1017 if src_kind in ("surface", "mixed"):1018 src_data.update(vertices_from=[s["vertno"].copy() for s in src_t[:2]])1019 1020 # delete copy1021 return src_data, src_kind, src_subject1022 1023 1024def _triage_output(output):1025 _check_option("output", output, ["nifti", "nifti1", "nifti2"])1026 if output in ("nifti", "nifti1"):1027 from nibabel import Nifti1Header as NiftiHeader1028 from nibabel import Nifti1Image as NiftiImage1029 else:1030 assert output == "nifti2"1031 from nibabel import Nifti2Header as NiftiHeader1032 from nibabel import Nifti2Image as NiftiImage1033 return NiftiImage, NiftiHeader1034 1035 1036def _interpolate_data(stc, morph, mri_resolution, mri_space, output):1037 """Interpolate source estimate data to MRI."""1038 _check_dep(nibabel="2.1.0", dipy=False)1039 NiftiImage, NiftiHeader = _triage_output(output)1040 _validate_type(stc, _BaseVolSourceEstimate, "stc", "volume source estimate")1041 assert morph.kind in ("volume", "mixed")1042 1043 voxel_size_defined = False1044 1045 if isinstance(mri_resolution, int | float) and not isinstance(mri_resolution, bool):1046 # use iso voxel size1047 mri_resolution = (float(mri_resolution),) * 31048 1049 if isinstance(mri_resolution, tuple):1050 _check_dep(nibabel=False, dipy="0.10.1") # nibabel was already checked1051 from dipy.align.reslice import reslice1052 1053 voxel_size = mri_resolution1054 voxel_size_defined = True1055 mri_resolution = True1056 1057 # if data wasn't morphed yet - necessary for call of1058 # stc_unmorphed.as_volume. Since only the shape of src is known, it cannot1059 # be resliced to a given voxel size without knowing the original.1060 if isinstance(morph, SourceSpaces):1061 assert morph.kind in ("volume", "mixed")1062 offset = 2 if morph.kind == "mixed" else 01063 if voxel_size_defined:1064 raise ValueError(1065 "Cannot infer original voxel size for reslicing... "1066 "set mri_resolution to boolean value or apply morph first."1067 )1068 # Now deal with the fact that we may have multiple sub-volumes1069 inuse = [s["inuse"] for s in morph[offset:]]1070 src_shape = [s["shape"] for s in morph[offset:]]1071 assert len(set(map(tuple, src_shape))) == 11072 src_subject = morph._subject1073 morph = BunchConst(src_data=_get_src_data(morph, mri_resolution)[0])1074 else:1075 # Make a list as we may have many inuse when using multiple sub-volumes1076 inuse = morph.src_data["inuse"]1077 src_subject = morph.subject_from1078 assert isinstance(inuse, list)1079 if stc.subject is not None:1080 _check_subject_src(stc.subject, src_subject, "stc.subject")1081 1082 n_times = stc.data.shape[1]1083 shape = morph.src_data["src_shape"][::-1] + (n_times,) # SAR->RAST1084 dtype = np.complex128 if np.iscomplexobj(stc.data) else np.float641085 # order='F' so that F-order flattening is faster1086 vols = np.zeros((np.prod(shape[:3]), shape[3]), dtype=dtype, order="F")1087 n_vertices_seen = 01088 for this_inuse in inuse:1089 this_inuse = this_inuse.astype(bool)1090 n_vertices = np.sum(this_inuse)1091 stc_slice = slice(n_vertices_seen, n_vertices_seen + n_vertices)1092 vols[this_inuse] = stc.data[stc_slice]1093 n_vertices_seen += n_vertices1094 1095 # use mri resolution as represented in src1096 if mri_resolution:1097 if morph.src_data["interpolator"] is None:1098 raise RuntimeError(1099 "Cannot morph with mri_resolution when add_interpolator=False "1100 "was used with setup_volume_source_space"1101 )1102 shape = morph.src_data["src_shape_full"][::-1] + (n_times,)1103 vols = morph.src_data["interpolator"] @ vols1104 1105 # reshape back to proper shape1106 vols = np.reshape(vols, shape, order="F")1107 1108 # set correct space1109 if mri_resolution:1110 affine = morph.src_data["src_affine_vox"]1111 else:1112 affine = morph.src_data["src_affine_src"]1113 1114 if mri_space:1115 affine = np.dot(morph.src_data["src_affine_ras"], affine)1116 1117 affine[:3] *= 1e31118 1119 # pre-define header1120 header = NiftiHeader()1121 header.set_xyzt_units("mm", "msec")1122 header["pixdim"][4] = 1e3 * stc.tstep1123 1124 # if a specific voxel size was targeted (only possible after morphing)1125 if voxel_size_defined:1126 # reslice mri1127 vols, affine = reslice(vols, affine, _get_zooms_orig(morph), voxel_size)1128 1129 with warnings.catch_warnings(): # nibabel<->numpy warning1130 vols = NiftiImage(vols, affine, header=header)1131 1132 return vols1133 1134 1135###############################################################################1136# Morph for VolSourceEstimate1137 1138 1139def _compute_morph_sdr(mri_from, mri_to, niter_affine, niter_sdr, zooms):1140 """Get a matrix that morphs data from one subject to another."""1141 from dipy.align.imaffine import AffineMap1142 1143 from .transforms import _compute_volume_registration1144 1145 pipeline = "all" if niter_sdr else "affines"1146 niter = dict(1147 translation=niter_affine,1148 rigid=niter_affine,1149 affine=niter_affine,1150 sdr=niter_sdr if niter_sdr else (1,),1151 )1152 (1153 pre_affine,1154 sdr_morph,1155 to_shape,1156 to_affine,1157 from_shape,1158 from_affine,1159 ) = _compute_volume_registration(1160 mri_from, mri_to, zooms=zooms, niter=niter, pipeline=pipeline1161 )1162 pre_affine = AffineMap(1163 pre_affine,1164 domain_grid_shape=to_shape,1165 domain_grid2world=to_affine,1166 codomain_grid_shape=from_shape,1167 codomain_grid2world=from_affine,1168 )1169 return to_shape, zooms, to_affine, pre_affine, sdr_morph1170 1171 1172def _compute_morph_matrix(1173 subject_from,1174 subject_to,1175 vertices_from,1176 vertices_to,1177 smooth=None,1178 subjects_dir=None,1179 warn=True,1180 xhemi=False,1181):1182 """Compute morph matrix."""1183 logger.info("Computing morph matrix...")1184 subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)1185 1186 tris = _get_subject_sphere_tris(subject_from, subjects_dir)1187 maps = read_morph_map(subject_from, subject_to, subjects_dir, xhemi)1188 1189 # morph the data1190 1191 morpher = []1192 for hemi_to in range(2): # iterate over to / block-rows of CSR matrix1193 hemi_from = (1 - hemi_to) if xhemi else hemi_to1194 morpher.append(1195 _hemi_morph(1196 tris[hemi_from],1197 vertices_to[hemi_to],1198 vertices_from[hemi_from],1199 smooth,1200 maps[hemi_from],