eigenstrapping.VolumetricEigenstrapping

class eigenstrapping.VolumetricEigenstrapping(data, volume, label=None, aseg=False, norm_file=None, normalization=None, normalization_factor=None, evals=None, emodes=None, num_modes=200, seed=None, decomp_method='matrix', randomize=False, resample=False, n_jobs=1, permute=False, verbose=True)

Compute the eigenmodes and eigenvalues of the volume in volume and resample the hypersphere bounded by the eigengroups contained within emodes, to reconstruct the data using coefficients conditioned on the original data Based on the degenerate solutions of solving the Laplace-Beltrami operator on the cortex. The power spectrum is perfectly retained (the square of the eigenvalues). If evals and emodes are given (i.e., precomputed) then eigenmodes are not computed. Performs amplitude adjustment by default (see resample)

@author: Nikitas C. Koussis, School of Psychological Sciences,

University of Newcastle & Systems Neuroscience Group Newcastle

Michael Breakspear, School of Psychological Sciences, University of Newcastle & Systems Neuroscience Group Newcastle

Important

Several software packages must be installed apart from python dependencies to use this class:

• FreeSurfer >5

• Gmsh

• fsl >5 (if aseg not passed)

Process

• the orthonormal eigenvectors n within eigengroup l give the surface

  of the hyperellipse of dimensions l-1 NOTE: for eigengroup 0 with the zeroth mode, we ignore it, and for eigengroup 1, this is resampling the surface of the 2-sphere

• the axes of the hyperellipse are given by the sqrt of the eigenvalues

  corresponding to eigenvectors n

• linear transform the eigenvectors N to the hypersphere by dividing

  by the ellipsoid axes

• finds the set of points p on the hypersphere given by the basis

  modes (the eigenmodes) by normalizing them to unit length

• rotate the set of points p by cos(angle) for

  even dimensions and sin(angle) for odd dims (resampling step)

• find the unit vectors of the points p by dividing by the

  Euclidean norm (equivalent to the eigenmodes)

• make the new unit vectors orthonormal using Gram-Schmidt process

• return the new eigenmodes of that eigengroup until all eigengroups

  are computed

• reconstruct the null data by multiplying the original coefficients

  by the new eigenmodes

• resamples the null data by performing rank-ordering to replicate

  the reconstructed data term, then adds the noise term back into the null data to produce a surrogate that replicates the original variance of the empirical data. * Only performed if resample = True.

References

1. Robinson, P. A. et al., (2016), Eigenmodes of brain activity: Neural field theory predictions and comparison with experiment. NeuroImage 142, 79-98. <https://dx.doi.org/10.1016/j.neuroimage.2016.04.050>

2. Jorgensen, M., (2014), Volumes of n-dimensional spheres and ellipsoids. <https://www.whitman.edu/documents/Academics/Mathematics/2014/jorgenmd.pdf>

3. <https://math.stackexchange.com/questions/3316373/ellipsoid-in-n-dimension>

4. Blumenson, L. E. (1960). A derivation of n-dimensional spherical coordinates. The American Mathematical Monthly, 67(1), 63-66. <https://www.jstor.org/stable/2308932>

5. Trefethen, Lloyd N., Bau, David III, (1997). Numerical linear algebra. Philadelphia, PA: Society for Industrial and Applied Mathematics. ISBN 978-0-898713-61-9.

6. https://en.wikipedia.org/wiki/QR_decomposition

7. Chen, Y. C. et al., (2022). The individuality of shape asymmetries of the human cerebral cortex. Elife, 11, e75056. <https://doi.org/10.7554/eLife.75056>

param volume:

Filename of volume ROI with integer intensity values. Must be in the same space (i.e. registered) to the volume file in data. Extensions are those commonly accepted by nibabel (see nibabel file types). Volume is extracted using FreeSurfer or fslmaths if aseg is not passed.

type volume:

str to image

param data:

Filename of empirical data to resample

type data:

str to image or str to text file

param label:

Label of value in volume for ROI extraction. If default None, volume is considered to be a mask - i.e., all values equal to 1 are extracted and transformed to tetrahedral space, and those equal to 0 or otherwise are left out. Default is None

type label:

int or list of ints or array of ints, optional

param aseg:

Specify whether file in volume is FreeSurfer preprocessed “aseg.mgz” file. Extra steps need to be performed for computation, though labels can be combined into structures (e.g., label=[11, 12, 26] with aseg=True would combine the left caudate, left putamen, and left NA to make the left striatum. See the FreeSurferColorLUT.txt for more details about labeling). Default is False

type aseg:

bool, optional

param evals:

Eigenvalues corresponding to the number of eigenmodes n in emodes. Default is None

type evals:

np.ndarray of shape (n,), optional

param emodes:

Eigenmodes that are the solution to the generalized eigenvalue problem or Helmholtz equation in the Laplace-Beltrami operator of the volume (must be in tetrahedral space). Default is None

type emodes:

np.ndarray of shape (n_vertices, n), optional

param num_modes:

If evals and emodes are not given, then they are computed on the surface given in surface. This variable controls how many modes to compute. Cannot exceed the number of vertices after tetrahedral resampling of volume. Default is 200.

type num_modes:

int, optional

param seed:

Specify the seed for random angle generation (or random state instance) default None.

type seed:

None or int or np.random.RandomState instance

param decomp_method:

method of calculation of coefficients: ‘matrix’, ‘matrix_separate’, ‘regression’.

The default is ‘matrix’.

type decomp_method:

str, optional

param resample:

Set whether to resample surrogate map from original map to preserve values, default True

type resample:

bool, optional

param randomize:

Set whether to shuffle coefficients calculated from minimizing least-squares error to reconstruct surrogates. Results in better randomization of surrogates at the cost of poorer replication of empirical spatial autocorrelation.

type randomize:

bool, optional

param n_jobs:

Number of workers to use for parallelization. Default 1.

type n_jobs:

int, optional

param use_cholmod:

Specify whether or not to use sksparse subroutines - requires scikit-sparse and libsuitesparse-dev to be installed. See https://github.com/scikit-sparse/scikit-sparse for more information.

type use_cholmod:

bool, optional

:raises ValueError : Inappropriate inputs:

calc_volume(nifti_volume=None, label=None)

Calculate the physical volume of the ROI in the nifti file

Parameters:

• nifti_volume (str, optional) – Filename of input volume in nifti format

• label (int, optional) – Label value for input ROI in nifti_volume

Returns:

• ROI_number (int) – Total number of non-zero voxels

• ROI_volume (float) – Total volume of non-zero voxels in physical dimensions

calc_volume_eigenmodes(tetra, num_modes=200, use_cholmod=False)

Calculate the eigenvalues and eigenmodes of the ROI volume.

Parameters:

• tetra (lapy compatible object) – Loaded lapy object corresponding to tetrahedral surface of ROI in vtk format

• num_modes (int, optional) – Number of eigenmodes to compute. The default is 200

• use_cholmod (bool, optional) – Specify whether to use scikit-sparse libraries to compute eigenmodes. Much faster but requires libraries to be installed.

Returns:

• evals (np.ndarray of shape (num_modes,)) – Eigenvalues

• emodes (np.ndarray of shape (number of tetra points, num_modes)) – Eigenmodes

calculate_distance_matrix(return_data=False)

Calculates Euclidean distance matrix and index memmap files of ROI, also returns data masked by ROI.

Parameters:

outdir (str, optional) – Output directory of distance and index np.memmap files.

Returns:

• masked_data ((N,))

• dist (np.memmap) – Distance memmap matrix.

• index (str) – Distance index memmap matrix

eigen_spectrum(show=True)

Plot the modal decomposition power spectrum. See geometry.eigen_decomposition() Returns fig and ax handles from matplotlib.pyplot for further use.

Parameters:

show (bool, optional) – Flag whether to show plot, default True

Returns:

• fig (matplotlib.pyplot class) – Figure handle

• ax (matplotlib.pyplot class) – Axes handle

generate(output_modes=False)

Generate eigensphere resampled surrogate.

Parameters:

output_modes (bool, optional) – Output resampled modes for debugging. The default is False

Returns:

(n,) – Surrogate data in tetrahedral space

Return type:

np.ndarray

make_tetra(volume, label=None, aseg=False, verbose=True)

Generate tetrahedral version of the ROI in the nifti file. Can specify label value.

Parameters:

• volume (str) – Filename of input volume where the relevant ROI has voxel values = 1 or label

• label (int or list of ints, optional) – Label value(s) of input volume. Extracts surface from voxels that have this(ese) intensity value(s). If None, defaults to 1

• aseg (bool, optional) – Specify whether input volume and label is from FreeSurfer aseg.nii.gz (post-FS preprocessing).

Returns:

tetra_file – Filename of output tetrahedral vtk file

Return type:

str

Raises:

RuntimeError – Multiple labels given but aseg not passed.

normalize_vtk(tetra, normalization_type=None, normalization_factor=None)

Normalize tetrahedral surface

Parameters:

• tetra (lapy compatible object) – Loaded vtk object corresponding to a surface tetrahedral mesh

• normalization_type (str or None, optional) – Type of normalization. The default is None.

• normalization_factor (float or None, optional) – Factor to be used in “constant” normalization. The default is None.

Returns:

tetra_norm – Loaded vtk object corresponding to normalized tetrahedral surface

Return type:

lapy compatible object

project_to_tetra(tetra, data)

Project data in volumetric space onto tetrahedral surface.

Parameters:

• tetra (lapy compatible object) – Loaded object corresponding to tetrahedral surface.

• data (np.ndarray of shape (xx, yy, zz)) – Data to project to tetrahedral surface with values inside ROI.

• interpolation (str, optional) – Interpolation kind. Default is ‘linear’

Returns:

tetra_data – Resampled data

Return type:

np.ndarray of shape (number of tetra points,)

resample_data(data, method='linear', vector=False)

Resample data on tetrahedral surface to volumetric space.

Parameters:

• data (np.ndarray of shape (n_tetra_points, N)) – Data to resample

• method (str, optional) – Interpolation method for scipy.interpolate.griddata. The default is ‘linear’.

• vector (bool, optional) – Resample data to masked vector (values inside ROI only). Default False

Returns:

resampled_data – Data in volumetric space.

Return type:

np.ndarray of shape (A,B,C,N)

rotate_modes(emodes)

Rotate modes using random rotations.

Parameters:

vectors (np.ndarray of (n_vertices, mu)) – Orthogonal vectors of n_vertices and mu modes

Returns:

rotated_vectors – Array containing rotated vectors for each pair.

Return type:

np.ndarray of (n_vertices, mu)

surfplot(data, view='dorsal', cmap='viridis', vrange=[-4, 4])

Plotting helper function for volumetric data to triangular surfaces.

Parameters:

• data (np.ndarray of shape (n_vertices,)) – Volumetric data projected to tetrahedral mesh. See self.project_to_tetra().

• view (str or pair of float, optional) – Specify which direction to place camera and view from. See nilearn.plotting.plot_surf() for accepted strings. If pair of float, must be a list or tuple of (elevation, azimuthal) of angles in degrees. The default is ‘dorsal’.

• cmap (str, optional) – Colormap for values in data. Must be in matplotlib colormaps. The default is ‘viridis’

• vrange (tuple of two floats, optional) – Range to truncate minimum and maximum values. The default is [-4, 4].

• annotations (bool, optional) – Specify whether to add annotation labels for direction etc

Returns:

• fig (matplotlib.pyplot.Figure handle)

• ax (matplotlib.pyplot.Axes handle)