SciPy

pynibs.mesh package

Submodules

pynibs.mesh.mesh_struct module

The mesh_struct.py module provides classes and methods for handling and manipulating mesh structures in the context of neuroimaging and brain stimulation studies. It includes classes for handling tetrahedral meshes and regions of interest (ROIs).

The module includes the following classes:

  • TetrahedraLinear: This class represents a mesh consisting of linear tetrahedra. It provides methods for calculating various quantities of interest (QOIs) in the mesh, interpolating data between nodes and elements, calculating the electric field and current density, and more.

  • Mesh: This class is a general mesh class to initialize default attributes and provides methods to fill default values based on the approach, write the mesh data to an HDF5 file, and print the mesh information.

  • ROI: This class represents a region of interest (ROI) in the mesh. It provides methods to write the ROI data to an HDF5 file and print the ROI information.

Each class and method in this module is documented with docstrings providing more detailed information about its purpose, parameters, and return values.

This module is primarily used for handling and visualizing data related to neuroimaging and brain stimulation studies.

class pynibs.mesh.mesh_struct.Mesh(mesh_name, subject_id, subject_folder)

Bases: object

Mesh class to initialize default attributes.

fill_defaults(approach)

Initializes attributes for a headreco mesh.

Parameters:

approach (str) – ‘headreco’ ‘mri2mesh’ ‘charm’

print()

Print self information.

write_to_hdf5(fn_hdf5, check_file_exist=False, verbose=False)

Write this mesh’ attributes to .hdf5 file.

Parameters:

  • fn_hdf5 (str) –

  • check_file_exist (bool) – Check if provided filenames exist, warn if not.

  • verbose (bool) – Print self information.

class pynibs.mesh.mesh_struct.ROI(subject_id, roi_name, mesh_name)

Bases: object

Region of interest class to initialize default attributes.

print()

Print self information.

write_to_hdf5(fn_hdf5, check_file_exist=False, verbose=False)

Write this mesh’ attributes to .hdf5 file.

Parameters:

  • fn_hdf5 (str) –

  • check_file_exist (bool) – Check if provided filenames exist, warn if not.

  • verbose (bool) – Print self information.

class pynibs.mesh.mesh_struct.TetrahedraLinear(points, triangles, triangles_regions, tetrahedra, tetrahedra_regions)

Bases: object

Mesh, consisting of linear tetrahedra.

Parameters:

  • points (array of float [N_points x 3]) – Vertices of FE mesh.

  • triangles (np.ndarray of int [N_tri x 3]) – Connectivity of points forming triangles.

  • triangles_regions (np.ndarray of int [N_tri x 1]) – Region identifiers of triangles.

  • tetrahedra (np.ndarray of int [N_tet x 4]) – Connectivity of points forming tetrahedra.

  • tetrahedra_regions (np.ndarray of int [N_tet x 1]) – Region identifiers of tetrahedra.

N_points

Number of vertices.

Type:

int

N_tet

Number of tetrahedra.

Type:

int

N_tri

Number of triangles.

Type:

int

N_region

Number of regions.

Type:

int

region

Region labels.

Type:

np.ndarray of int

tetrahedra_volume

Volumes of tetrahedra.

Type:

np.ndarray of float [N_tet x 1]

tetrahedra_center

Center of tetrahedra.

Type:

np.ndarray of float [N_tet x 1]

triangles_center

Center of triangles.

Type:

np.ndarray of float [N_tri x 1]

triangles_normal

Normal components of triangles pointing outwards.

Type:

np.ndarray of float [N_tri x 3]

calc_E(grad_phi, omegaA)

Calculate electric field with gradient of electric potential and omega-scaled magnetic vector potential A.

\mathbf{E}=-\nabla\varphi-\omega\mathbf{A}

Parameters:

  • grad_phi (np.ndarray of float) – (N_tet, 3) Gradient of Scalar DOF in tetrahedra center.

  • omegaA (np.ndarray of float) – (N_tet, 3) Magnetic vector potential in tetrahedra center (scaled with angular frequency omega).

Returns:

E – (N_tet, 3) Electric field in tetrahedra center.

Return type:

np.ndarray of float

calc_E_normal_tangential_surface(E, fname)

Calculate normal and tangential component of electric field on given surfaces of mesh instance.

Parameters:

  • E (np.ndarray of float [N_tri x 3]) – Electric field data on surfaces.

  • fname (str) – Filename of the .txt file containing the tetrahedra indices, which are adjacent to the surface triangles generated by the method “calc_surface_adjacent_tetrahedra_idx_list(self, fname)”.

Returns:

  • En_pos (np.ndarray of float [N_tri x 3]) – Normal component of electric field of top side (outside) of surface.

  • En_neg (np.ndarray of float [N_tri x 3]) – Normal component of electric field of bottom side (inside) of surface.

  • n (np.ndarray of float [N_tri x 3]) – Normal vector

  • Et (np.ndarray of float [N_tri x 3]) – Tangential component of electric field lying in surface.

  • t (np.ndarray of float [N_tri x 3]) – Tangential vector.

calc_E_on_GM_WM_surface(E, roi)

Determines the normal and tangential component of the induced electric field on a GM-WM surface using nearest neighbour principle.

Parameters:

  • E (np.ndarray of float [N_tri x 3]) – Induced electric field given in the tetrahedra centre of the mesh instance.

  • roi (pynibs.roi.RegionOfInterestSurface) – RegionOfInterestSurface object class instance.

Returns:

  • E_normal (np.ndarray of float [N_points x 3]) – Normal vector of electric field on GM-WM surface.

  • E_tangential (np.ndarray of float [N_points x 3]) – Tangential vector of electric field on GM-WM surface.

calc_E_on_GM_WM_surface3(phi, dAdt, roi, verbose=True, mode='components')

Determines the normal and tangential component of the induced electric field on a GM-WM surface by recalculating phi and dA/dt in an epsilon environment around the GM/WM surface (upper and lower GM-WM surface).

Parameters:

  • phi (np.ndarray of float) – (N_nodes, 1) Scalar electric potential given in the nodes of the mesh.

  • dAdt (np.ndarray of float) – (N_nodes, 3) Magnetic vector potential given in the nodes of the mesh.

  • roi (pynibs.mesh.mesh_struct.ORI) – RegionOfInterestSurface object class instance.

  • verbose (bool) – Print information to stdout.

  • mode (str) – Select mode of output: - “components” : return x, y, and z component of tangential and normal components - “magnitude” : return magnitude of tangential and normal component (normal with sign for direction)

Returns:

  • E_normal (np.ndarray of float) – (N_nodes, 3) Normal vector of electric field on GM-WM surface.

  • E_tangential (np.ndarray of float) – (N_nodes, 3) Tangential vector of electric field on GM-WM surface.

calc_E_on_GM_WM_surface_simnibs(phi, dAdt, roi, subject, verbose=False, mesh_idx=0)

Determines the normal and tangential component of the induced electric field on a GM-WM surface by recalculating phi and dA/dt in an epsilon environment around the GM/WM surface (upper and lower GM-WM surface) or by using the Simnibs interpolation function.

Parameters:

  • phi (np.ndarray of float) – (N_nodes, 1) Scalar electric potential given in the nodes of the mesh.

  • dAdt (np.ndarray of float) – (N_nodes, 3) Magnetic vector potential given in the nodes of the mesh.

  • roi (pynibs.mesh.mesh_struct.ROI) – RegionOfInterestSurface object class instance.

  • subject (pynibs.subject.Subject) – Subject object loaded from .hdf5 file.

  • verbose (bool) – Print information to stdout.

  • mesh_idx (int) – Mesh index.

Returns:

  • E_normal (np.ndarray of float) – (N_points, 3) Normal vector of electric field on GM-WM surface.

  • E_tangential (np.ndarray of float) – (N_points, 3) Tangential vector of electric field on GM-WM surface.

calc_E_on_GM_WM_surface_simnibs_KW(phi, dAdt, roi, subject, verbose=False, mesh_idx=0)

Determines the normal and tangential component of the induced electric field on a GM-WM surface by recalculating phi and dA/dt in an epsilon environment around the GM/WM surface (upper and lower GM-WM surface) or by using the Simnibs interpolation function.

Parameters:

  • phi (np.ndarray of float) – (N_nodes, 1) Scalar electric potential given in the nodes of the mesh.

  • dAdt (np.ndarray of float) – (N_nodes, 1) Magnetic vector potential given in the nodes of the mesh.

  • roi (pynibs.mesh.mesh_struct.ROI) – RegionOfInterestSurface object class instance.

  • subject (pynibs.subject.Subject) – Subject object loaded from .hdf5 file.

  • verbose (bool) – Print information to stdout.

  • mesh_idx (int) – Mesh index.

Returns:

  • E_normal (np.ndarray of float) – (N_points, 3) Normal vector of electric field on GM-WM surface.

  • E_tangential (np.ndarray of float) – (N_points, 3) Tangential vector of electric field on GM-WM surface.

calc_J(E, sigma)

Calculate current density J. The conductivity sigma is a list of np.arrays containing conductivities of regions (scalar and/or tensor).

\mathbf{J} = [\sigma]\mathbf{E}

Parameters:

  • E (np.ndarray of float) – (N_tet, 3) Electric field in tetrahedra center.

  • sigma (list of np.ndarray of float) – [N_regions](3, 3) Conductivities of regions (scalar and/or tensor).

Returns:

E – (N_tet, 3) Electric field in tetrahedra center.

Return type:

np.ndarray of float

calc_QOI_in_points(qoi, points_out)

Calculate QOI_out in points_out using the mesh instance and the quantity of interest (QOI).

Parameters:

  • qoi (np.ndarray of float) – Quantity of interest in nodes of tetrahedra mesh instance.

  • points_out (np.ndarray of float) – Point coordinates (x, y, z) where the qoi is going to be interpolated by linear basis functions.

Returns:

qoi_out – Quantity of interest in points_out.

Return type:

np.ndarray of float

calc_QOI_in_points_tet_idx(qoi, points_out, tet_idx)

Calculate QOI_out in points_out sitting in tet_idx using the mesh instance and the quantity of interest (QOI).

Parameters:

  • qoi (np.ndarray of float) – Quantity of interest in nodes of tetrahedra mesh instance.

  • points_out (np.ndarray of float) – Point coordinates (x, y, z) where the qoi is going to be interpolated by linear basis functions.

  • tet_idx (np.ndarray of int) – Element indices where the points_out are sitting.

Returns:

qoi_out – Quantity of interest in points_out

Return type:

np.ndarray of float

calc_gradient(phi)

Calculate gradient of scalar DOF in tetrahedra center.

Parameters:

phi (np.ndarray of float [N_nodes]) – Scalar DOF the gradient is calculated for.

Returns:

grad_phi – Gradient of Scalar DOF in tetrahedra center.

Return type:

np.ndarray of float [N_tet x 3]

calc_surface_adjacent_tetrahedra_idx_list(fname)

Determine the indices of the tetrahedra touching the surfaces and save the indices into a .txt file specified with fname.

Parameters:

fname (str) – Filename of output .txt file.

Returns:

<File> – Element indices of the tetrahedra touching the surfaces (outer-most elements).

Return type:

.txt file

data_elements2nodes(data)

Transforms data in tetrahedra into the nodes after Zienkiewicz et al. (1992) [1]. Can only transform volume data, i.e. needs the data in the surrounding tetrahedra to average it to the nodes. Will not work well for discontinuous fields (like E, if several tissues are used).

Parameters:

data (np.ndarray [N_elements x N_data]) – Data in tetrahedra.

Returns:

data_nodes – Data in nodes.

Return type:

np.ndarray [N_nodes x N_data]

Notes

[1]

Zienkiewicz, Olgierd Cecil, and Jian Zhong Zhu. “The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique.” International Journal for Numerical Methods in Engineering 33.7 (1992): 1331-1364.

data_nodes2elements(data)

Interpolate data given in the nodes to the tetrahedra center.

Parameters:

data (np.ndarray [N_nodes x N_data]) – Data in nodes.

Returns:

data_elements – Data in elements.

Return type:

np.ndarray [N_elements x N_data]

get_faces(tetrahedra_indexes=None)

Creates a list of nodes in each face and a list of faces in each tetrahedron.

Parameters:

tetrahedra_indexes (np.ndarray) – Indices of the tetrehedra where the faces are to be determined (default: all tetrahedra).

Returns:

  • faces (np.ndarray) – List of nodes in faces, in arbitrary order

  • th_faces (np.ndarray.) – List of faces in each tetrahedron, starts at 0, order=((0, 2, 1), (0, 1, 3), (0, 3, 2), (1, 2, 3)).

  • face_adjacency_list (np.ndarray) – List of tetrahedron adjacent to each face, filled with -1 if a face is in a single tetrahedron. Not in the normal element ordering, but only in the order the tetrahedra are presented.

get_outside_faces(tetrahedra_indices=None)

Creates a list of nodes in each face that are in the outer volume.

Parameters:

tetrahedra_indices (np.ndarray) – Indices of the tetrahedra where the outer volume is to be determined (default: all tetrahedra).

Returns:

faces – List of nodes in faces in arbitrary order.

Return type:

np.ndarray

pynibs.mesh.transformations module

pynibs.mesh.transformations.cell_data_to_point_data(tris, data_tris, nodes, method='nearest')

A wrapper for scipy.interpolate.griddata to interpolate cell data to node data.

Parameters:

  • tris (np.ndarray) – (n_tri, 3) element number list.

  • data_tris (np.ndarray) – (n_tri x 3) data in tris.

  • nodes (np.ndarray) – (n_nodes, 3) nodes coordinates.

  • method (str, default: 'nearest') – Which method to use for interpolation. Default uses NearestNDInterpolator.

Returns:

data_nodes – Data in nodes.

Return type:

np.ndarray

pynibs.mesh.transformations.cell_data_to_point_data_vtk(mesh=None, nodes=None, con=None, cell_data=None, fn=None)

Convert cell data to point data in a VTK unstructured grid and save the result to a file.

Parameters:

  • mesh (meshio.Mesh) – The mesh object containing points and cells.

  • nodes (np.ndarray of float, optional) – (N_points, 3) Coordinates of the nodes.

  • con (np.ndarray of int, optional) – (N_elements, 3) Connectivity index list forming the elements.

  • cell_data (dict, optional) – Cell data to be transformed to point data. keys: str, values: np.ndarray

  • fn (str, optional) – If provided, vtk file is written out to this file.

Returns:

dict – All data sets from mesh transformed to point data.

Return type:

point_data

pynibs.mesh.transformations.data_elements2nodes(data, con, precise=False)

Transforms data in elements (triangles or tetrahedra) to nodes. Data can be list of multiple data arrays.

Parameters:

  • data (np.ndarray of float or list of np.ndarray) – (N_elements, N_data) Data given in the elements (multiple datasets who fit to con may be passed in a list).

  • con (np.ndarray of int) – triangles: (N_elements. 3). tetrahedra: (N_elements, 4). Connectivity index list forming the elements.

  • precise (bool, default: False) – Compute data transformation precisely but slow. Better for near-0 values.

Returns:

out – (N_nodes, N_data) Data in nodes.

Return type:

np.ndarray of float or list of np.ndarray

pynibs.mesh.transformations.data_nodes2elements(data, con)

Transforms data in nodes to elements (triangles or tetrahedra).

Parameters:

  • data (np.ndarray of float) – (N_nodes, N_data) Data given in the nodes.

  • con (np.ndarray of int) – triangles: (N_elements, 3). tetrahedra: (N_elements, 4). Connectivity index list forming the elements.

Returns:

out – (N_elements, N_data) Data given in the element centers.

Return type:

np.ndarray of float

pynibs.mesh.transformations.map_data_to_surface(datasets, points_datasets, con_datasets, fname_fsl_gm, fname_fsl_wm, fname_midlayer=None, delta=0.5, input_data_in_center=True, return_data_in_center=True, data_substitute=-1)

Maps data from ROI of fsl surface (wm, gm, or midlayer) to given Freesurfer brain surface (wm, gm, inflated).

Parameters:

  • datasets (np.ndarray of float [N_points x N_data] or list of np.ndarray) – Data in nodes or center of triangles in ROI (specify this in “data_in_center”).

  • points_datasets (np.ndarray of float [N_points x 3] or list of np.ndarray) – Point coordinates (x,y,z) of ROI where data in datasets list is given, the points have to be a subset of the GM/WM surface (has to be provided for each dataset).

  • con_datasets (np.ndarray of int [N_tri x 3] or list of np.ndarray) – Connectivity matrix of dataset points (has to be provided for each dataset).

  • fname_fsl_gm (str or list of str or list of None) – Filename of pial surface fsl file(s) (one or two hemispheres) e.g. in mri2msh: …/fs_ID/surf/lh.pial

  • fname_fsl_wm (str or list of str or list of None) – Filename of wm surface fsl file(s) (one or two hemispheres) e.g. in mri2msh: …/fs_ID/surf/lh.white

  • fname_midlayer (str or list of str) – Filename of midlayer surface fsl file(s) (one or two hemispheres) e.g. in headreco: …/fs_ID/surf/lh.central

  • delta (float) – Distance parameter where gm-wm surface was generated 0…1 (default: 0.5) 0 -> WM surface 1 -> GM surface

  • input_data_in_center (bool) – Flag if data in datasets in given in triangle centers or in points (Default: True).

  • return_data_in_center (bool) – Flag if data should be returned in nodes or in elements (Default: True).

  • data_substitute (float) – Data substitute with this number for all points in the inflated brain, which do not belong to the given data set.

Returns:

data_mapped – Mapped data to target brain surface. In points or elements.

Return type:

np.ndarray of float [N_points_inf x N_data]

pynibs.mesh.transformations.map_nii2surf(tri_node_coords, nii, datatype=<class 'float'>, nonzeroneighbors=False)

Map nifti data to mesh surface.

For each surface triangle,

  • the center coordinate is determined

  • and the corresponding voxel value from the nifti file is extracted.

Raw, unsmoothed data is returned. Nifti file and mesh needs to be in the same space.

Nifti volume on mesh surface

Volume data from nifti file is mapped to mesh surface.

fn_roi = "geo.hdf5"
fn_nii = "some_T_map.nii"

nii = nib.load(fn_nii)

with h5py.File(fn_roi, 'r') as f:
    triangle_number_list = f['/mesh/elm/triangle_number_list'][:]
    node_coord = f['/mesh/nodes/node_coord'][:]

nii_data_on_surf = pynibs.map_nii2surf(node_coord[triangle_number_list], nii)
output_fn = "some_T_map_on_surf.vtk"
point_data = {"T": [nii_data_on_surf]}
meshio.Mesh(np.squeeze(node_coord), [('triangle', triangle_number_list)], cell_data=point_data).write(output_fn)
Parameters:

  • tri_node_coords (np.ndarray of int) – (n,3,3) list of node coordinates for n trianlges (mesh surface).

  • nii (nibabel.Nifti1Image or str) – Nifti image or filename of nifti.

  • datatype (type, default: float) – Data type of the output array. For masks, use int.

  • nonzeroneighbors (boolean, default=False) – If True, zero datapoints are replaced by largest absolute neihbors voxel data. This is useful for masked niftis, where data outside of the mask is zero and you want to pull data only from inside the mask.

Returns:

(n,) array with voxel values from nifti file for each triangle center.

Return type:

np.ndarray of float

pynibs.mesh.transformations.midlayer_2_surf(midlayer_data, coords_target, coords_midlayer, midlayer_con=None, midlayer_data_in_nodes=False, max_dist=5, outside_roi_val=0, precise_map=True)

Convert midlayer data to whole-brain surface data, e.g. grey matter. Output is returned as data in nodes.

Parameters:

  • midlayer_data (np.ndarray of float) – (n_elm_midlayer,) or (n_nodes_midlayer,), the data in the midlayer.

  • coords_target (np.ndarray of float) – (n_nodes_target, 3) Coordinates of the nodes of the target surface.

  • coords_midlayer (np.ndarray of float) – (n_nodes_midlayer, 3) Coordinates of the nodes of the midlayer surface.

  • midlayer_con (np.ndarray of int, optional) – (n_elm_midlayer, 3) Connectivity of the midlayer elements. Provide if data_in_points == True.

  • midlayer_data_in_nodes (bool, default=False) – If midlayer data is provided in nodes, set to True and provide midlayer_con.

  • max_dist (float, default=5) – Maximum distance between target and midlayer nodes to pull data from midlayer_data for.

  • outside_roi_val (float, default=0) – Areas outside of max_dist are filled with outside_roi_val.

  • precise_map (bool, default=True) – If elements to nodes mapping is done, perform this precise and slow or not.

Returns:

data_target – (n_nodes_target, 1) The data in nodes of the target surface.

Return type:

np.ndarray

pynibs.mesh.transformations.point_data_to_cell_data_vtk(mesh=None, nodes=None, con=None, point_data=None, fn=None)

Convert point data to cell data in a VTK unstructured grid and save the result to a file.

Parameters:

  • mesh (meshio.Mesh, optional) – The mesh object containing points and cells.

  • nodes (np.ndarray of float, optional) – (N_points, 3) Coordinates of the nodes.

  • con (np.ndarray of int, optional) – (N_elements, 3) Connectivity index list forming the elements.

  • point_data (dict, optional) – Point data to be transformed to cell data.

  • fn (str, optional) – If provided, vtk file is written out to this file.

Returns:

dict – All data sets from mesh transformed to cell_data.

Return type:

cell_data

pynibs.mesh.transformations.project_on_scalp(coords, mesh, scalp_tag=1005)

Find the node in the scalp closest to each coordinate.

Parameters:

  • coords (nx3 np.ndarray) – Vectors to be transformed.

  • mesh (pynibs.TetrahedraLinear or simnibs.msh.mesh_io.Msh) – Mesh structure in simnibs or pynibs format.

  • scalp_tag (int, default: 1005) – Tag in the mesh where the scalp is to be set.

Returns:

points_closest – (n, 3) coordinates projected scalp (closest skin points).

Return type:

np.ndarry

pynibs.mesh.transformations.project_on_scalp_hdf5(coords, mesh, scalp_tag=1005)

Find the node in the scalp closest to each coordinate.

Parameters:

  • coords (np.ndarray) – (n, 3) Vectors to be transformed.

  • mesh (str or pynibs.TetrahedraLinear) – Filename of mesh in .hdf5 format or Mesh structure.

  • scalp_tag (int, default: 1005) – Tag in the mesh where the scalp is to be set.

Returns:

points_closest – (n, 3) coordinates projected scalp (closest skin points).

Return type:

np.ndarray

pynibs.mesh.transformations.refine_surface(fn_surf, fn_surf_refined, center, radius, repair=True, remesh=True, verbose=True)

Refines surface (.stl) in spherical ROI and saves as .stl file.

Parameters:

  • fn_surf (str) – Input filename (.stl).

  • fn_surf_refined (str) – Output filename (.stl).

  • center (np.ndarray of float) –

    1. Center of spherical ROI (x,y,z).

  • radius (float) – Radius of ROI.

  • repair (bool, default: True) – Repair surface mesh to ensure that it is watertight and forms a volume.

  • remesh (bool, default: False) – Perform remeshing with meshfix (also removes possibly overlapping facets and intersections).

  • verbose (bool, default: True) – Print output messages.

Returns:

<file>

Return type:

.stl file

pynibs.mesh.transformations.write_vtu(fn, vtk_grid)

pynibs.mesh.utils module

pynibs.mesh.utils.calc_distances(coords, mesh_fn, tissues=None)

Calculates the distances between coords and tissue types.

Parameters:

  • coords (list of list or list or np.ndarray) – Coordinates (X, Y, Z) to compute depths for.

  • mesh_fn (str) – pynibs.Mesh hdf5 filename.

  • tissues (list of int, optional) – Which tissue types to compute depths for. If none, distances to all tissue types are computed.

Returns:

distances – colunms: coorrd, tissue_type, distance.

Return type:

pd.Dataframe()

pynibs.mesh.utils.calc_gradient_surface(phi, points, triangles)

Calculate gradient of potential phi on surface (i.e. tangential component) given in vertices of a triangular mesh forming a 2D surface.

Parameters:

  • phi (np.ndarray of float [N_points x 1]) – Potential in nodes.

  • points (np.ndarray of float [N_points x 3]) – Coordinates of nodes (x,y,z).

  • triangles (np.ndarray of int32 [N_tri x 3]) – Connectivity of triangular mesh.

Returns:

grad_phi – Gradient of potential phi on surface.

Return type:

np.ndarray of float [N_tri x 3]

pynibs.mesh.utils.calc_tet_volume(points, abs=True)

Calculate tetrahedra volumes.

Parameters:

  • points (np.ndarray) –

    shape: (n_tets,4,3)

    [[[Ax, Ay, Az],
  [Bx, By, Bz],
  [Cx, Cy, Cz],
  [Dx, Dy, Dy]],

 [[Ax, Ay, Az],
  [Bx, By, Bz],
  [Cx, Cy, Cz],
  [Dx, Dy, Dy]],
 ...
]

  • abs (bool, default: true) – Return magnitude-

Returns:

volume – shape: (n_tets)-

Return type:

np.ndarray

pynibs.mesh.utils.calc_tetrahedra_volume_cross(P1, P2, P3, P4)

Calculates volume of tetrahedra specified by the 4 points P1…P4 multiple tetrahedra can be defined by P1…P4 as 2-D np.ndarrays using the cross and vector dot product.

P1=\begin{bmatrix} x_{{tet}_1} & y_{{tet}_1} & z_{{tet}_1} \\ x_{{tet}_2} & y_{{tet}_2} & z_{{tet}_2} \\ ... & ... & ... \\ x_{{tet}_N} & y_{{tet}_N} & z_{{tet}_N} \\ \end{bmatrix}

Parameters:

  • P1 (np.ndarray of float [N_tet x 3]) – Coordinates of first point of tetrahedra.

  • P2 (np.ndarray of float [N_tet x 3]) – Coordinates of second point of tetrahedra.

  • P3 (np.ndarray of float [N_tet x 3]) – Coordinates of third point of tetrahedra.

  • P4 (np.ndarray of float [N_tet x 3]) – Coordinates of fourth point of tetrahedra.

Returns:

tetrahedra_volume – Volumes of tetrahedra.

Return type:

np.ndarray of float [N_tet x 1]

pynibs.mesh.utils.calc_tetrahedra_volume_det(P1, P2, P3, P4)

Calculate volume of tetrahedron specified by 4 points P1…P4 multiple tetrahedra can be defined by P1…P4 as 2-D np.ndarray using the determinant.

P1=\begin{bmatrix} x_{{tet}_1} & y_{{tet}_1} & z_{{tet}_1} \\ x_{{tet}_2} & y_{{tet}_2} & z_{{tet}_2} \\ ... & ... & ... \\ x_{{tet}_N} & y_{{tet}_N} & z_{{tet}_N} \\ \end{bmatrix}

Parameters:

  • P1 (np.ndarray of float [N_tet x 3]) – Coordinates of first point of tetrahedra.

  • P2 (np.ndarray of float [N_tet x 3]) – Coordinates of second point of tetrahedra.

  • P3 (np.ndarray of float [N_tet x 3]) – Coordinates of third point of tetrahedra.

  • P4 (np.ndarray of float [N_tet x 3]) – Coordinates of fourth point of tetrahedra.

Returns:

tetrahedra_volume – Volumes of tetrahedra.

Return type:

np.ndarray of float [N_tet x 1]

pynibs.mesh.utils.calc_tri_surface(points)

Calculate triangle surface areas.

Parameters:

points (np.ndarray) – (n_triangles,3,3)

Returns:

triangle_area

Return type:

np.ndarray

pynibs.mesh.utils.check_islands_for_single_elm(source_elm, connectivity=None, adjacency=None, island_crit=1)

This identifies islands in a mesh for a given element. An island is a set of elements, that is only connect via a single node to another set of elements. These islands usually crash the FEM solver and should be removed.

  1. Find all elements connect to source_elm via one node (1-node-neighbor)

  2. Start with source_elm and visit all 2-node-neighbors (‘shared-edge)

  3. Continue recursively with all 2-node-neighbors and visit their 2-node-neighbors

  4. See if any 1-node-neighbors have not been visited with this strategy. If so, an island has been found

Parameters:

  • source_elm (int) – The source element to check

  • connectivity (np.ndarray, optional) – Connectivity (‘node_number_list’) starting with 0. Can be triangles or tetrahedra (n_elms, 3) or (n_elms_4).

  • adjacency (np.ndparray, optional) – Adjenceny matrix (n_elm, n_elm). Weights are supposed to be number of shared nodes. Computed from neighbors if not provided.

  • island_crit (int, default: 'any') – How many nodes to define islands? ‘any’ -> Elements connected via a single node or single edge are defined as an island. ‘node’ -> Elements connected via a single _node_ are defined as an island. ‘edge’ -> Elements connected via a single _edge_ are defined as an island.

Returns:

  • n_visited (int)

  • n_not_visited (int)

  • neighbors_visited (dict, which neighbors have been visited and which have not)

pynibs.mesh.utils.cortical_depth(mesh_fn, geo_fn=None, skin_surface_id=1005, mesh_fn_out=None, tissue_type=None, only_tri=False, verbose=False)

Compute skin-cortex-distance (SCD) for surface and volume data in mesh_fn.

Visualized cortical depth.

Cortical depth computed against skin surface.

Parameters:

  • mesh_fn (str) – TetrahedraLinear mesh file.

  • geo_fn (str, optional) – TetrahedraLinear mesh file with geometric data. If provided, geometric information is read from here.

  • skin_surface_id (int, default: 1005) – Which tissue type nr to compute distance against.

  • mesh_fn_out (str, optional) – If provided, SCD information is written to this file.

  • tissue_type (list of int, optional) – For which tissue types to compute SCD. If None, all tissue types are computed.

  • only_tri (bool, default: False) – If True, only triangles are computed.

  • verbose (bool, default: False) – Print some verbosity information.

Returns:

  • distances_tri (np.ndarray) – Distances of triangles to skin surface.

  • distances_tets (np.ndarray) – Distances of tetrahedra to skin surface.

  • <file> (.hdf5) – mesh_fn or geo_fn with SCD information in /data/tris/Cortex_dist and /data/tets/Cortex_dist.

  • <file> (.xdmf) – Only if write_xdmf == True.

pynibs.mesh.utils.determine_e_midlayer(fn_e_results, fn_mesh_hdf5, subject, mesh_idx, roi_idx, n_cpu=4, midlayer_fun='simnibs', phi_scaling=1.0, verbose=False)

Parallel version to determine the midlayer e-fields from a list of .hdf5 results files

Parameters:

  • fn_e_results (list of str) – List of results filenames (.hdf5 format).

  • fn_mesh_hdf5 (str) – Filename of corresponding mesh file.

  • subject (pynibs.Subject) – Subject object.

  • mesh_idx (int) – Mesh index.

  • roi_idx (int) – ROI index.

  • n_cpu (int, default: 4) – Number of parallel computations.

  • midlayer_fun (str, default: "simnibs") – Method to determine the midlayer e-fields (“pynibs” or “simnibs”).

  • phi_scaling (float, default: 1.0) – Scaling factor of scalar potential to change between “m” and “mm”.

Returns:

Adds midlayer e-field results to ROI.

Return type:

<File> .hdf5 file

pynibs.mesh.utils.determine_e_midlayer_workhorse(fn_e_results, subject, mesh_idx, midlayer_fun, fn_mesh_hdf5, roi_idx, phi_scaling=1.0, verbose=False)
phi_scaling: float

simnibs < 3.0 : 1000. simnibs >= 3.0 : 1. (Default)

pynibs.mesh.utils.find_element_idx_by_points(nodes, con, points)

Finds the tetrahedral element index of an arbitrary point in the FEM mesh.

Parameters:

  • nodes (np.ndarray [N_nodes x 3]) – Coordinates (x, y, z) of the nodes.

  • con (np.ndarray [N_tet x 4]) – Connectivity matrix.

  • points (np.ndarray [N_points x 3]) – Points for which the element indices are found.

Returns:

ele_idx – Element indices of tetrahedra where corresponding ‘points’ are lying in.

Return type:

np.ndarray [N_points]

pynibs.mesh.utils.find_island_elms(connectivity=None, adjacency=None, verbose=False, island_crit='edge', decision='cumulative')

Searches for islands in a mesh and returns element indices of the smallest island. Island is defind as a set of elements, which share a single node and/or single edge with the rest of the mesh.

Parameters:

  • connectivity (np.ndarray, optional) – Connectivity (‘node_number_list’) starting with 0. Can be triangles or tetrahedra (n_elms, 3) or (n_elms_4).

  • adjacency (np.ndparray, optional) – Adjenceny matrix (n_elm, n_elm). Weights are supposed to be number of shared nodes. Computed from neighbors if not provided.

  • island_crit (int, default: 'edge') – How many nodes to define islands? ‘node’ -> Elements connected via a single _node_ are defined as an island. ‘edge’ -> Elements connected via a single _edge_ are defined as an island.

  • decision (str, default: cumulative) – ‘cumulative’ -> Return all element indices that are not visited any times ‘smallest’ -> Return smallest island.

  • verbose (bool, optional) – Print some verbosity information. Default: False

Returns:

island

Return type:

list of island-elms.

pynibs.mesh.utils.find_islands(connectivity=None, adjacency=None, island_crit='any', verbose=False, largest=False)

This identifies islands in a mesh. An island is a set of elements, that is only connect via a single node to another set of elements. These islands usually crash the FEM solver and should be removed.

For each element:

  1. Find all elements connect to source_elm via one node (1-node-neighbor)

  2. Start with source_elm and visit all 2-node-neighbors (‘shared-edge)

  3. Continue recursively with all 2-node-neighbors and visit their 2-node-neighbors

  4. See if any 1-node-neighbors have not been visited with this strategy. If so, an island has been found

Island detection

Islands are groups of elements that are only connected via a single node/edge to another group.

Parameters:

  • connectivity (np.ndarray, optional) – Connectivity (‘node_number_list’) starting with 0. Can be triangles or tetrahedra (n_elms, 3) or (n_elms_4).

  • adjacency (np.ndparray, optional) – Adjenceny matrix (n_elm, n_elm). Weights are supposed to be number of shared nodes. Computed from neighbors if not provided.

  • island_crit (int or str, default: 'any') – How many nodes to define islands? ‘any’ -> Elements connected via a single node or single edge are defined as an island. ‘node’ -> Elements connected via a single _node_ are defined as an island. ‘edge’ -> Elements connected via a single _edge_ are defined as an island.

  • largest (book, default: False) – Only return largest island, speeds up computation quite a bit if only one large, and many small islands exist.

  • verbose (bool, optional) – Print some verbosity information. Default: False.

Returns:

  • elms_with_island (list) – Elements with neighboring islands.

  • counter_visited (np.ndarray) – shape = (n_elms). How often as each element been visited.

  • counter_not_visited (np.ndarray) – shape = (n_elms). How often as each element not been visited.

pynibs.mesh.utils.find_nearest(array, value)

Given an “array”, and given a “value” , returns an index j such that “value” is between array[j] and array[j+1]. “array” must be monotonic increasing. j=-1 or j=len(array) is returned to indicate that “value” is out of range below and above respectively.

Parameters:

  • array (np.ndarray of float) – Monotonic increasing array.

  • value (float) – Target value the nearest neighbor index in array is computed for.

Returns:

idx – Index j such that “value” is between array[j] and array[j+1].

Return type:

int

pynibs.mesh.utils.get_indices_discontinuous_data(data, con, neighbor=False, deviation_factor=2, min_val=None, not_fitted_elms=None, crit='median', neigh_style='point')

Get element indices (and the best neighbor index), where the data is discontinuous.

Parameters:

  • data (np.ndarray of float [n_data]) – Data array to analyze given in the element center.

  • con (np.ndarray of float [n_data, 3 or 4]) – Connectivity matrix.

  • neighbor (bool, default: False) – Return also the element index of the “best” neighbor (w.r.t. median of data).

  • deviation_factor (float) – Allows data deviation from 1/deviation_factor < data[i]/median < deviation_factor.

  • min_val (float, optional) – If given, only return elements which have a neighbor with data higher than min_val.

  • not_fitted_elms (np.ndarray) – If given, these elements are not used as neighbors.

  • crit (str, default: median) – Criterium for best neighbor. Either median or max value.

  • neigh_style (str, default: 'point') – Should neighbors share point or ‘edge’.

Returns:

  • idx_disc (list of int [n_disc]) – Index list containing the indices of the discontinuous elements.

  • idx_neighbor (list of int [n_disc]) – Index list containing the indices of the “best” neighbors of the discontinuous elements.

pynibs.mesh.utils.get_sphere(mesh=None, mesh_fn=None, target=None, radius=None, roi_idx=None, roi=None, elmtype='tris', domain=None)

Return element idx of elements within a certain distance to provided target. Element indices are 0-based (tris and tets start at 0, ‘pynibs’ style) Elements might be ‘tris’ (default) or ‘tets’

If roi object / idx and mesh fn is provided, the roi is expected to have midlayer information and the roi geometry is used.

Parameters:

  • mesh (pynibs.mesh.mesh_struct.TetrahedraLinear, optional) –

  • mesh_fn (str, optional) – Filename to SimNIBS .msh or pyNIBS .hdf5 mesh file.

  • target (np.ndarray of float or list of float) – (3,) X, Y, Z coordinates of target.

  • radius (float) – Sphere radius im mm.

  • roi_idx (str or int, optional) – ROI name.

  • elmtype (str, default: 'tris') – Return triangles or tetrahedra in sphere around target. One of (‘tris’, ‘tets’).

Returns:

elms_in_sphere – (n_elements): Indices of elements found in ROI-

Return type:

np.ndarray

pynibs.mesh.utils.head_circumference(fn_hdf5, z_center=None, result_folder=None)

Computes the head circumference of a mesh.

Head circumference

For each point on the circumference, the closest point on the skin surface is determined. Distances from point to point are summed up to get the head circumference.

Parameters:

  • fn_hdf5 (str) – Filename of the mesh in hdf5 format.

  • z_center (float, optional) – Z-coordinate of the center of the head. If not provided, it is computed from the mesh.

  • result_folder (str, optional) – Folder to save the results. If not provided, the results are not saved.

Returns:

circ – Head circumference in mm.

Return type:

float

pynibs.mesh.utils.in_hull(points, hull)

Test if points in points are in hull. points should be a [N x K] coordinates of N points in K dimensions. hull is either a scipy.spatial.Delaunay object or the [M x K] array of the coordinates of M points in Kdimensions for which Delaunay triangulation will be computed.

Parameters:

  • points (np.ndarray) – (N_points x 3) Set of floating point data to test whether they are lying inside the hull or not.

  • hull (scipy.spatial.Delaunay or np.ndarray) – (M x K) Surface data.

Returns:

inside – TRUE: point inside the hull. FALSE: point outside the hull.

Return type:

np.ndarray of bool

pynibs.mesh.utils.sample_sphere(n_points, r)

Creates n_points evenly spread in a sphere of radius r.

Parameters:

  • n_points (int) – Number of points to be spread, must be odd.

  • r (float) – Radius of sphere.

Returns:

points – (N x 3), Evenly spread points in a unit sphere.

Return type:

np.ndarray of float

pynibs.mesh.utils.tets_in_sphere(mesh, target, radius, roi, domain=None)

Worker function for get_sphere()

Returns element idx of elements within a certain distance to provided target. If roi object / idx and mesh fn is provided, the roi is expected to have midlayer information and the roi geometry is used.

If radius is None or 0, the nearest element is returned.

Parameters:

  • mesh (pynibs.TetrahedraLinear, optional) –

  • target (np.ndarray of float, optional) – (3,) X, Y, Z coordinates of target-

  • radius (float, optional) – Sphere radius im mm-

  • roi (pynibs.mesh.ROI, optional) – Region of interest-

Returns:

tets_in_sphere – (n_tets): Indices of elements found in ROI-

Return type:

np.ndarray

pynibs.mesh.utils.tris_in_sphere(mesh, target, radius, roi)

Worker function for get_sphere().

Returns triangle idx of elements within a certain distance to provided target. If roi object / idx and mesh fn is provided, the roi is expected to have midlayer information and the roi geometry is used.

If radius is None or 0, the nearest element is returned.

Parameters:
Returns:

tris_in_sphere – (n_triangles): Indices of elements found in sphere.

Return type:

np.ndarray

Table of Contents

Previous topic

pynibs.hdf5_io package

This Page

Edit page
  • © Copyright 2026, Konstantin Weise, Ole Numssen.
  • Created using Sphinx 6.2.1.