OSCR

Multilayer MEG source modelling enables depth-resolved laminar inference in humans.

Code ↔ Paper

16 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 16 matches
  1. [1] § Methods › Simulations ↔ lameg/surf.py, lines 681 to 806 · score 0.93 · cortical patch statistics, surface normal vectors, original surface normal, downsampled surface normal, mesh normals, pial surface · Report a mismatch
  2. [2] § Methods › Empirical data preprocessing ↔ lameg/surf.py, lines 858 to 915 · score 0.83 · intermediate surfaces, FreeSurfer, removes vertices, pial surface, link vectors, corresponding vertices · Report a mismatch
  3. [3] § Results › Cortical column orientation and laminar inference ↔ lameg/surf.py, lines 681 to 806 · score 0.81 · cortical patch statistics, white matter surfaces, link vectors, surface normals, corresponding vertices, dipole orientation · Report a mismatch
  4. [4] § Methods › Empirical data analysis ↔ lameg/invert.py, lines 474 to 554 · score 0.79 · multiple sparse priors, Hann windowing, sliding window, temporal mode, overlapping, free energy · Report a mismatch
  5. [5] § Methods › Simulations ↔ lameg/util.py, lines 895 to 943 · score 0.79 · left preauricular, right preauricular, fiducial coordinates, nasion, scanner, LPA · Report a mismatch
  6. [6] § Methods › Empirical data analysis ↔ run_motor_epoch_01_mf.py, lines 294 to 379 · score 0.75 · composite anatomical suitability, lead field RMS, cortical thickness, motor, score, ranked · Report a mismatch
  7. [7] § Methods › Empirical data analysis ↔ run_motor_epoch_02_mefI.py, lines 293 to 378 · score 0.75 · composite anatomical suitability, lead field RMS, cortical thickness, motor, score, ranked · Report a mismatch
  8. [8] § Methods › Simulation analysis ↔ lameg/invert.py, lines 474 to 554 · score 0.73 · multiple sparse priors, Hann windowed, temporal modes, free energy, source reconstruction, MSP · Report a mismatch
  9. [9] § Methods › Simulation analysis ↔ 06_run_multilaminar_VEF.ipynb, lines 148 to 198 · score 0.72 · dot product, scalp normal, scalp surface, cortical column, trees, cross · Report a mismatch
  10. [10] § Methods › Simulation analysis ↔ pipeline/run_motor_epoch_01_mf.py, lines 112 to 162 · score 0.72 · dot product, scalp normal, scalp surface, cortical column, trees, cross · Report a mismatch
  11. [11] § Methods › Empirical data preprocessing ↔ lameg/surf.py, lines 917 to 957 · score 0.68 · intermediate surface, deep vertices, layer mesh, mris, hemispheres · Report a mismatch
  12. [12] § Methods › Simulations ↔ lameg/simulate.py, lines 34 to 151 · score 0.65 · noise ratio, dipole moment, White noise, forward model, density, patch · Report a mismatch
  13. [13] § Methods › Empirical data analysis ↔ lameg/invert.py, lines 288 to 359 · score 0.63 · multiple sparse priors, temporal modes, beamformer, MSP, channels, patch · Report a mismatch
  14. [14] § Methods › Simulations ↔ config/spm_cfg_opm_create.m, lines 108 to 170 · score 0.61 · inner skull, outer skull, scalp, sensor, cortical · Report a mismatch
  15. [15] § Results › SNR and co-registration constraints on laminar inference ↔ lameg/laminar.py, lines 25 to 120 · score 0.53 · empirical Bayesian beamformer, metric, EBB, free energy, algorithm, MSP · Report a mismatch
  16. [16] § Results › Anatomical features determine laminar inference accuracy ↔ figures/plot_figure_08_anatomy.ipynb, lines 589 to 639 · score 0.50 · Optimal thresholds, zero error, marginal, joint, Figure 8, mm · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,216 lines · 46 KB · GPL-3.0 · 4 matches

  1. """
  2. This module provides a set of tools for handling and manipulating surface mesh data. The
  3. functionalities include computing mesh normals, interpolating mesh data, handling non-manifold
  4. edges, creating and manipulating GIFTI surface files, and downsampling meshes using the VTK
  5. library.
  7. Key functionalities include:
  8. - Normalization of vectors to unit length.
  9. - Calculation of normals for mesh surfaces using both Delaunay triangulation and custom methods.
  10. - Creation of GIFTI images from mesh data.
  11. - Removal of specified vertices from a mesh and updating the mesh topology accordingly.
  12. - Identification and handling of non-manifold edges to ensure mesh manifoldness.
  13. - Interpolation of data from a downsampled mesh back to its original high-resolution mesh.
  14. - Downsampling of meshes using VTK's decimation algorithms.
  15. - Combination and adjustment of multiple surface meshes into a single mesh.
  16. """
  18. # pylint: disable=C0302
  19. import os
  20. import copy
  21. import subprocess
  22. from collections import defaultdict
  24. from joblib import Parallel, delayed
  26. import nibabel as nib
  27. import numpy as np
  28. from scipy.spatial import cKDTree # pylint: disable=E0611
  29. from scipy.sparse import csr_matrix
  31. # pylint: disable=E0611
  32. from vtkmodules.vtkCommonCore import vtkPoints
  33. from vtkmodules.vtkCommonDataModel import vtkPolyData, vtkCellArray
  34. from vtkmodules.vtkFiltersCore import vtkDecimatePro
  35. from vtkmodules.util.numpy_support import vtk_to_numpy
  37. def _normit(vectors):
  38. """
  39. Normalize a numpy array of vectors.
  41. This function normalizes each row in the array of vectors to have unit length. If the length of
  42. a vector is below a certain threshold (machine epsilon), it is set to 1 to avoid division by
  43. zero.
  45. Parameters
  46. ----------
  47. vectors : ndarray
  48. Array of vectors to be normalized. Each row represents a vector.
  50. Returns
  51. -------
  52. ndarray
  53. Normalized array of vectors where each row has unit length.
  54. """
  56. norm_n = np.sqrt(np.sum(vectors ** 2, axis=1))
  57. norm_n[norm_n < np.finfo(float).eps] = 1
  58. return vectors / norm_n[:, np.newaxis]
  61. def mesh_normals(vertices, faces, unit=False):
  62. """
  63. Normalize a numpy array of vectors.
  65. This function normalizes each row in the array of vectors to have unit length. If the length of
  66. a vector is below a certain threshold (machine epsilon), it is set to 1 to avoid division by
  67. zero.
  69. Parameters
  70. ----------
  71. vectors : ndarray
  72. Array of vectors to be normalized. Each row represents a vector.
  74. Returns
  75. -------
  76. ndarray
  77. Normalized array of vectors where each row has unit length.
  78. """
  80. face_normal = np.cross(
  81. vertices[faces[:, 1], :] - vertices[faces[:, 0], :],
  82. vertices[faces[:, 2], :] - vertices[faces[:, 0], :]
  83. )
  84. face_normal = _normit(face_normal)
  86. vertex_normal = np.zeros_like(vertices)
  87. for i in range(len(faces)):
  88. for j in range(3):
  89. vertex_normal[faces[i, j], :] += face_normal[i, :]
  91. centered_vertices = vertices - np.mean(vertices, axis=0)
  92. if np.count_nonzero(np.sign(np.sum(centered_vertices * vertex_normal, axis=1))) > \
  93. len(centered_vertices) / 2:
  94. vertex_normal = -vertex_normal
  95. face_normal = -face_normal
  97. if unit:
  98. vertex_normal = _normit(vertex_normal)
  99. face_normal = _normit(face_normal)
  101. return vertex_normal, face_normal
  104. def create_surf_gifti(vertices, faces, normals=None):
  105. """
  106. Create a Gifti image object from surface mesh data.
  108. This function creates a GiftiImage object from the provided vertices, faces, and optional
  109. normals. The vertices and faces are required, while normals are optional. If normals are
  110. provided, they are added to the Gifti image. The function returns the GiftiImage object.
  112. Parameters
  113. ----------
  114. vertices : numpy.ndarray
  115. Array of vertices. Each row represents a vertex with its x, y, z coordinates.
  116. faces : numpy.ndarray
  117. Array of faces. Each row represents a face with three integers corresponding to vertex
  118. indices.
  119. normals : numpy.ndarray, optional
  120. Array of vertex normals. Each row represents a normal vector corresponding to a vertex.
  122. Returns
  123. -------
  124. new_gifti : nibabel.gifti.GiftiImage
  125. The GiftiImage object created from the provided mesh data.
  127. Notes
  128. -----
  129. - Vertex, face, and normal arrays should be NumPy arrays.
  130. - Vertices and normals should be in float32 format, and faces should be in int32 format.
  132. Example
  133. -------
  134. >>> import numpy as np
  135. >>> import nibabel as nib
  136. >>> vertices = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]])
  137. >>> faces = np.array([[0, 1, 2], [0, 2, 3]])
  138. >>> normals = np.array([[0, 0, 1], [0, 0, 1], [0, 0, 1], [0, 0, 1]])
  139. >>> gifti_img = create_surf_gifti(vertices, faces, normals)
  140. """
  142. # Create new gifti object
  143. new_gifti = nib.gifti.GiftiImage()
  145. # Cast vertices and faces to the appropriate data types
  146. vertices = vertices.astype(np.float32)
  147. faces = faces.astype(np.int32)
  149. # Add the vertices and faces to the gifti object
  150. new_gifti.add_gifti_data_array(
  151. nib.gifti.GiftiDataArray(
  152. data=vertices,
  153. intent=nib.nifti1.intent_codes['NIFTI_INTENT_POINTSET']
  154. )
  155. )
  156. new_gifti.add_gifti_data_array(
  157. nib.gifti.GiftiDataArray(
  158. data=faces,
  159. intent=nib.nifti1.intent_codes['NIFTI_INTENT_TRIANGLE']
  160. )
  161. )
  163. # If normals are provided and not empty, cast them to float32 and add them to the Gifti image
  164. if normals is not None:
  165. normals = np.array(normals).astype(np.float32)
  166. new_gifti.add_gifti_data_array(
  167. nib.gifti.GiftiDataArray(
  168. data=normals,
  169. intent=nib.nifti1.intent_codes['NIFTI_INTENT_VECTOR']
  170. ))
  172. return new_gifti
  175. def remove_unconnected_vertices(gifti_surf):
  176. """
  177. Remove vertices that are not connected to any faces from a Gifti surface object.
  179. Parameters
  180. ----------
  181. gifti_surf : nibabel.gifti.GiftiImage
  182. The Gifti surface object to be processed.
  184. Returns
  185. -------
  186. cleaned_gifti_surf : nibabel.gifti.GiftiImage
  187. A new GiftiImage object with unconnected vertices removed.
  188. """
  190. # Get the pointset (vertices) and triangle array (faces) from the Gifti surface
  191. vertices = gifti_surf.darrays[0].data
  192. faces = gifti_surf.darrays[1].data
  194. # Find all unique vertex indices that are used in faces
  195. connected_vertices = np.unique(faces.flatten())
  197. # Determine which vertices are not connected to any faces
  198. all_vertices = np.arange(vertices.shape[0])
  199. unconnected_vertices = np.setdiff1d(all_vertices, connected_vertices)
  201. # Remove unconnected vertices using the provided remove_vertices function
  202. cleaned_gifti_surf = remove_vertices(gifti_surf, unconnected_vertices)
  204. return cleaned_gifti_surf
  207. def remove_vertices(gifti_surf, vertices_to_remove):
  208. """
  209. Remove specified vertices from a Gifti surface and update the faces accordingly.
  211. This function modifies a Gifti surface by removing the specified vertices. It also updates
  212. the faces of the surface so that they only reference the remaining vertices. If normals
  213. are present in the surface, they are also updated to correspond to the new set of vertices.
  215. Parameters
  216. ----------
  217. gifti_surf : nibabel.gifti.GiftiImage
  218. The Gifti surface object from which vertices will be removed.
  219. vertices_to_remove : array_like
  220. An array of vertex indices to be removed from the surface.
  222. Returns
  223. -------
  224. new_gifti : nibabel.gifti.GiftiImage
  225. A new GiftiImage object with the specified vertices removed and faces updated.
  227. Notes
  228. -----
  229. - The function assumes that the GiftiImage object contains at least two data arrays: one for
  230. vertices and one for faces. If normals are present, they are also updated.
  231. - Vertex indices in `vertices_to_remove` should be zero-based (following Python's indexing
  232. convention).
  233. - The returned GiftiImage object is a new object; the original `gifti_surf` object is not
  234. modified in place.
  236. Example
  237. -------
  238. >>> import nibabel as nib
  239. >>> gifti_surf = nib.load('path_to_gifti_file.gii')
  240. >>> vertices_to_remove = np.array([0, 2, 5]) # Indices of vertices to remove
  241. >>> new_gifti_surf = remove_vertices(gifti_surf, vertices_to_remove)
  242. """
  244. # Extract vertices and faces from the gifti object
  245. vertices_data = [da for da in gifti_surf.darrays
  246. if da.intent == nib.nifti1.intent_codes['NIFTI_INTENT_POINTSET']][0]
  247. faces_data = [da for da in gifti_surf.darrays
  248. if da.intent == nib.nifti1.intent_codes['NIFTI_INTENT_TRIANGLE']][0]
  250. vertices = vertices_data.data
  251. faces = faces_data.data
  253. # Determine vertices to keep
  254. vertices_to_keep = np.setdiff1d(np.arange(vertices.shape[0]), vertices_to_remove)
  256. # Create new array of vertices
  257. new_vertices = vertices[vertices_to_keep, :]
  259. # Find which faces to keep - ones that point to kept vertices
  260. if faces.shape[0]>0:
  261. face_x = np.isin(faces[:, 0], vertices_to_keep)
  262. face_y = np.isin(faces[:, 1], vertices_to_keep)
  263. face_z = np.isin(faces[:, 2], vertices_to_keep)
  264. faces_to_keep = np.where(face_x & face_y & face_z)[0]
  266. # Re-index faces
  267. x_faces = faces[faces_to_keep, :].reshape(-1)
  268. idxs = np.searchsorted(vertices_to_keep, x_faces)
  269. new_faces = idxs.reshape(-1, 3)
  270. else:
  271. new_faces = faces
  273. # Create new gifti object
  274. normals = None
  275. normals_data = [da for da in gifti_surf.darrays
  276. if da.intent == nib.nifti1.intent_codes['NIFTI_INTENT_VECTOR']]
  277. if normals_data:
  278. normals = normals_data[0].data[vertices_to_keep, :]
  279. new_gifti = create_surf_gifti(new_vertices, new_faces, normals=normals)
  281. return new_gifti
  284. def find_non_manifold_edges(faces):
  285. """
  286. Identify non-manifold edges in a given mesh represented by its faces.
  288. A non-manifold edge is defined as an edge that is shared by more than two faces. This function
  289. processes an array of faces, each face represented by a tuple of vertex indices, and identifies
  290. edges that meet the non-manifold criteria.
  292. Parameters
  293. ----------
  294. faces : np.ndarray
  295. An array where each row represents a face as a tuple of three vertex indices.
  297. Returns
  298. -------
  299. non_manifold_edges : dict
  300. A dictionary where keys are tuples representing non-manifold edges (vertex indices are
  301. sorted), and values are lists of face indices that share the edge.
  303. Notes
  304. -----
  305. The function uses a `defaultdict` to collect face indices for each edge encountered in the mesh.
  306. It then filters out edges that are associated with more than two faces, identifying them as
  307. non-manifold.
  308. """
  310. edge_faces = defaultdict(list)
  312. for i, (vertex_1, vertex_2, vertex_3) in enumerate(faces):
  313. for edge in [(vertex_1, vertex_2), (vertex_2, vertex_3), (vertex_3, vertex_1)]:
  314. edge_faces[tuple(sorted(edge))].append(i)
  316. non_manifold_edges = {edge: fcs for edge, fcs in edge_faces.items() if len(fcs) > 2}
  317. return non_manifold_edges
  320. def fix_non_manifold_edges(vertices, faces):
  321. """
  322. Remove faces associated with non-manifold edges from a mesh defined by vertices and faces.
  324. Non-manifold edges are edges that are shared by more than two faces, which can cause issues
  325. in various mesh processing tasks such as mesh simplification, smoothing, or 3D printing. This
  326. function identifies such edges and removes all faces associated with them to ensure
  327. manifoldness of the mesh.
  329. Parameters
  330. ----------
  331. vertices : np.ndarray
  332. An array of vertices, where each row represents a vertex as [x, y, z] coordinates.
  333. faces : np.ndarray
  334. An array of faces, where each row represents a face as indices into the vertices array.
  336. Returns
  337. -------
  338. vertices : np.ndarray
  339. The unchanged array of vertices.
  340. new_faces : np.ndarray
  341. The modified array of faces, with faces associated with non-manifold edges removed.
  343. Examples
  344. --------
  345. >>> vertices = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]])
  346. >>> faces = np.array([[0, 1, 2], [0, 2, 3], [1, 2, 3]])
  347. >>> new_vertices, new_faces = fix_non_manifold_edges(vertices, faces)
  348. >>> new_faces
  349. array([[0, 1, 2], [0, 2, 3]]) # Assuming face [1, 2, 3] was associated with a non-manifold
  350. edge.
  351. """
  353. non_manifold_edges = find_non_manifold_edges(faces)
  354. conflicting_faces = set()
  355. for faces_list in non_manifold_edges.values():
  356. conflicting_faces.update(faces_list)
  358. # Create a new face list excluding the conflicting faces
  359. new_faces = np.array(
  360. [face for i, face in enumerate(faces) if i not in conflicting_faces],
  361. dtype=np.int32
  362. )
  363. return vertices, new_faces
  366. def downsample_single_surface(gifti_surf, ds_factor=0.1):
  367. """
  368. Downsample a Gifti surface using the VTK library.
  370. This function takes a Gifti surface defined by its vertices and faces, and downsamples it using
  371. VTK's `vtkDecimatePro` algorithm. The reduction ratio determines the degree of downsampling.
  372. The function returns the downsampled Gifti surface.
  374. Parameters
  375. ----------
  376. gifti_surf : nibabel.gifti.GiftiImage
  377. The Gifti surface object to be downsampled.
  378. reduction_ratio : float
  379. The proportion of the mesh to remove. For example, a reduction ratio of 0.1 retains 90% of
  380. the original mesh.
  382. Returns
  383. -------
  384. new_gifti_surf : nibabel.gifti.GiftiImage
  385. A new GiftiImage object with the downsampled surface.
  387. Notes
  388. -----
  389. - The input faces array should be triangulated, i.e., each face should consist of exactly three
  390. vertex indices.
  391. - The VTK library is used for mesh decimation, which must be installed and properly configured.
  392. - The returned GiftiImage object is a new object; the original `gifti_surf` object is not
  393. modified in place.
  395. Example
  396. -------
  397. >>> import numpy as np
  398. >>> gifti_surf = nib.load('path_to_gifti_file.gii')
  399. >>> new_gifti_surf = downsample_single_surface(gifti_surf, 0.1)
  400. """
  402. vertices = gifti_surf.darrays[0].data
  403. faces = gifti_surf.darrays[1].data
  405. # Convert vertices and faces to a VTK PolyData object
  406. points = vtkPoints()
  407. for point in vertices:
  408. points.InsertNextPoint(point)
  410. cells = vtkCellArray()
  411. for face in faces:
  412. cells.InsertNextCell(len(face))
  413. for vertex in face:
  414. cells.InsertCellPoint(vertex)
  416. polydata = vtkPolyData()
  417. polydata.SetPoints(points)
  418. polydata.SetPolys(cells)
  420. # Apply vtkDecimatePro for decimation
  421. decimate = vtkDecimatePro()
  422. decimate.SetInputData(polydata)
  423. decimate.SetTargetReduction(1 - ds_factor)
  424. decimate.Update()
  426. # Extract the decimated mesh
  427. decimated_polydata = decimate.GetOutput()
  429. # Convert back to numpy arrays
  430. reduced_vertices = vtk_to_numpy(decimated_polydata.GetPoints().GetData())
  432. # Extract and reshape the face data
  433. face_data = vtk_to_numpy(decimated_polydata.GetPolys().GetData())
  434. # Assuming the mesh is triangulated, every fourth item is the size (3), followed by three
  435. # vertex indices
  436. reduced_faces = face_data.reshape(-1, 4)[:, 1:4]
  438. # Find the original vertices closest to the downsampled vertices
  439. kdtree = cKDTree(gifti_surf.darrays[0].data)
  440. _, orig_vert_idx = kdtree.query(reduced_vertices, k=1)
  441. orig_vert_idx = np.squeeze(orig_vert_idx)
  443. reduced_normals = None
  444. if len(gifti_surf.darrays) > 2 and \
  445. gifti_surf.darrays[2].intent == nib.nifti1.intent_codes['NIFTI_INTENT_VECTOR'] and \
  446. gifti_surf.darrays[2].data.shape[0]==gifti_surf.darrays[0].data.shape[0]:
  447. reduced_normals = gifti_surf.darrays[2].data[orig_vert_idx]
  449. new_gifti_surf = create_surf_gifti(reduced_vertices, reduced_faces, normals=reduced_normals)
  451. return new_gifti_surf
  454. def iterative_downsample_single_surface(gifti_surf, ds_factor=0.1):
  455. """
  456. Iteratively downsample a single surface mesh to a target number of vertices.
  458. This function reduces the number of vertices in a surface mesh (in GIFTI format) to a specified
  459. fraction of its original size. Downsampling is performed iteratively until the target number of
  460. vertices is reached or closely approximated.
  462. Parameters
  463. ----------
  464. gifti_surf : nibabel.gifti.GiftiImage
  465. The surface mesh to be downsampled, provided as a GIFTI image object.
  466. ds_factor : float, optional
  467. The downsampling factor representing the target fraction of the original number of vertices.
  468. Default is 0.1.
  470. Returns
  471. -------
  472. current_surf : nibabel.gifti.GiftiImage
  473. The downsampled surface mesh as a GIFTI image object.
  475. Notes
  476. -----
  477. - The downsampling process is iterative. In each iteration, the mesh is downsampled by a factor
  478. calculated to approach the target number of vertices.
  479. - If the calculated downsampling factor in an iteration equals or exceeds 1, the process is
  480. terminated to prevent upsampling or infinite loops.
  481. """
  483. current_surf = gifti_surf
  484. current_vertices = gifti_surf.darrays[0].data.shape[0]
  485. target_vertices = int(current_vertices * ds_factor)
  486. current_ds_factor = target_vertices / current_vertices
  488. while current_vertices > target_vertices:
  489. # Downsample the mesh
  490. current_surf = downsample_single_surface(current_surf, ds_factor=current_ds_factor)
  492. # Update the current vertices
  493. current_vertices = current_surf.darrays[0].data.shape[0]
  495. current_ds_factor = (target_vertices / current_vertices) * 1.25
  496. if current_ds_factor >= 1:
  497. break
  499. # Remove non-manifold edges
  500. ds_vertices = current_surf.darrays[0].data
  501. ds_faces = current_surf.darrays[1].data
  502. nonmani_vertices, nonmani_faces = fix_non_manifold_edges(ds_vertices, ds_faces)
  504. normals = None
  505. if len(current_surf.darrays) > 2 and \
  506. current_surf.darrays[2].intent == nib.nifti1.intent_codes['NIFTI_INTENT_VECTOR'] and \
  507. current_surf.darrays[2].data.shape[0]==current_surf.darrays[0].data.shape[0]:
  508. normals = current_surf.darrays[2].data
  510. current_surf = create_surf_gifti(nonmani_vertices, nonmani_faces, normals=normals)
  512. # Remove unconnected vertices
  513. current_surf = remove_unconnected_vertices(current_surf)
  515. return current_surf
  518. def downsample_multiple_surfaces(in_surfs, ds_factor):
  519. """
  520. Downsample multiple surface meshes using the VTK decimation algorithm.
  522. This function takes a list of input surface meshes (in Gifti format) and applies a downsampling
  523. process to each surface. The downsampling is performed using VTK's `vtkDecimatePro` algorithm.
  524. The first surface in the list is downsampled, and its vertex mapping is then applied to all
  525. other surfaces in the list. The function returns a list of downsampled surface meshes.
  527. Parameters
  528. ----------
  529. in_surfs : list of nibabel.gifti.GiftiImage
  530. Input Gifti surface meshes to be downsampled.
  531. ratio : float
  532. The reduction ratio for the downsampling process. For example, a ratio of 0.1
  533. implies that the mesh will be reduced to 90% of its original size.
  535. Returns
  536. -------
  537. out_surfs : list of nibabel.gifti.GiftiImage
  538. List of downsampled Gifti surface meshes.
  540. Notes
  541. -----
  542. - The function prints the percentage of vertices retained in the first surface after
  543. downsampling.
  544. - If normals are present in the input surfaces, they are also downsampled and mapped to the
  545. new surfaces.
  546. - The resulting surfaces maintain the original topology and are suitable for visualization and
  547. further processing.
  549. Example
  550. -------
  551. >>> import nibabel as nib
  552. >>> in_surfs = [nib.load('path/to/input_surf1.gii'), nib.load('path/to/input_surf2.gii')]
  553. >>> ratio = 0.1
  554. >>> out_surfs = downsample_multiple_surfaces(in_surfs, ratio)
  555. >>> for i, ds_surf in enumerate(out_surfs):
  556. ... nib.save(ds_surf, f'path/to/output_surf{i+1}.gii')
  557. """
  559. out_surfs = []
  561. primary_surf = in_surfs[0]
  562. ds_primary_surf = iterative_downsample_single_surface(primary_surf, ds_factor=ds_factor)
  563. reduced_vertices = ds_primary_surf.darrays[0].data
  564. reduced_faces = ds_primary_surf.darrays[1].data
  566. # Find the original vertices closest to the downsampled vertices
  567. kdtree = cKDTree(primary_surf.darrays[0].data)
  568. # Calculate the percentage of vertices retained
  569. decim_orig_dist, orig_vert_idx = kdtree.query(reduced_vertices, k=1)
  570. orig_vert_idx = np.squeeze(orig_vert_idx)
  572. print(
  573. f"{(1 - np.mean(decim_orig_dist > 0)) * 100}% of the vertices in the decimated surface "
  574. f"belong to the original surface."
  575. )
  577. # Save the downsampled primary surface with normals
  578. out_surfs.append(ds_primary_surf)
  580. # Process other surfaces
  581. for i in range(1, len(in_surfs)):
  582. surf = in_surfs[i]
  584. reduced_normals = None
  585. if len(surf.darrays) > 2 and \
  586. surf.darrays[2].intent == nib.nifti1.intent_codes['NIFTI_INTENT_VECTOR'] and \
  587. surf.darrays[2].data.shape[0] == surf.darrays[0].data.shape[0]:
  588. reduced_normals = surf.darrays[2].data[orig_vert_idx]
  590. surf_verts=surf.darrays[0].data[orig_vert_idx, :]
  592. nonmani_vertices, nonmani_faces = fix_non_manifold_edges(surf_verts, reduced_faces)
  594. ds_surf = create_surf_gifti(nonmani_vertices, nonmani_faces, normals=reduced_normals)
  596. out_surfs.append(ds_surf)
  597. return out_surfs
  600. def combine_surfaces(surfaces):
  601. """
  602. Combine multiple surface meshes into a single surface mesh.
  604. This function takes a list of Gifti surface meshes and combines them into a single surface
  605. mesh. It concatenates the vertices, faces, and normals (if present) from each surface. The
  606. faces are re-indexed appropriately to maintain the correct references to the combined vertex
  607. array.
  609. Parameters
  610. ----------
  611. surfaces : list of nibabel.gifti.GiftiImage
  612. List of Gifti surface meshes to be combined.
  614. Returns
  615. -------
  616. combined_surf : nibabel.gifti.GiftiImage
  617. A single combined Gifti surface mesh.
  619. Notes
  620. -----
  621. - The vertices, faces, and normals (if present) from each surface are concatenated.
  622. - The faces are re-indexed to reference the correct vertices in the combined vertex array.
  623. - If normals are present in any of the input surfaces, they are also combined.
  625. Raises
  626. ------
  627. ValueError
  628. If the vertex or face arrays do not have the expected dimensions.
  630. Example
  631. -------
  632. >>> import nibabel as nib
  633. >>> surfaces = [nib.load('path/to/surface1.gii'), nib.load('path/to/surface2.gii')]
  634. >>> combined_surf = combine_surfaces(surfaces)
  635. >>> nib.save(combined_surf, 'path/to/combined_surface.gii')
  636. """
  638. combined_vertices = []
  639. combined_faces = []
  640. combined_normals = []
  642. face_offset = 0
  644. for mesh in surfaces:
  646. # Extract vertices and faces
  647. vertices = np.concatenate(
  648. [da.data for da in mesh.darrays
  649. if da.intent == nib.nifti1.intent_codes['NIFTI_INTENT_POINTSET']])
  650. faces = np.concatenate(
  651. [da.data for da in mesh.darrays
  652. if da.intent == nib.nifti1.intent_codes['NIFTI_INTENT_TRIANGLE']])
  654. # Check for normals
  655. normal_arrays = [da.data for da in mesh.darrays
  656. if da.intent == nib.nifti1.intent_codes['NIFTI_INTENT_VECTOR']]
  657. normals = np.concatenate(normal_arrays) if normal_arrays else np.array([])
  659. combined_vertices.append(vertices)
  660. combined_faces.append(faces + face_offset)
  661. if normals.size:
  662. combined_normals.append(normals)
  664. face_offset += vertices.shape[0]
  666. # Combine the arrays
  667. combined_vertices = np.vstack(combined_vertices).astype(np.float32)
  668. combined_faces = np.vstack(combined_faces).astype(np.int32)
  669. if combined_normals:
  670. combined_normals = np.vstack(combined_normals).astype(np.float32)
  672. combined_surf = create_surf_gifti(
  673. combined_vertices,
  674. combined_faces,
  675. normals=combined_normals
  676. )
  678. return combined_surf
  681. # pylint: disable=R0912,R0915
  682. def compute_dipole_orientations(method, layer_names, surf_dir, fixed=True):
  683. """
  684. Compute dipole orientations for cortical layers using different methods.
  686. Parameters
  687. ----------
  688. method : str
  689. Method for computing dipole orientations ('link_vector', 'ds_surf_norm', 'orig_surf_norm',
  690. or 'cps').
  691. link_vector: Vectors connecting pial vertices to corresponding white matter vertices.
  692. ds_surf_norm: Surface normal vectors computed from the downsampled surface.
  693. orig_surf_norm: Surface normal vectors computed from the original (non-downsampled) surface.
  694. cps: Cortical patch statistics - mean surface normal vectors from connected vertices in the
  695. original (non-downsampled) surface.
  696. layer_names : list
  697. Names of the cortical layers.
  698. surf_dir : str
  699. Directory where the surface files are stored.
  700. fixed : bool, optional
  701. Flag to ensure that orientation of corresponding vertices across layers is the same (True
  702. by default). If True, for 'ds_surf_norm', 'orig_surf_norm', and 'cps', orientations
  703. computed from the pial surface are used for all layers.
  705. Returns
  706. -------
  707. orientations : np.ndarray
  708. An array of dipole orientations for each vertex in each layer.
  710. Raises
  711. ------
  712. ValueError
  713. If the number of vertices in pial and white matter surfaces do not match.
  714. """
  716. orientations = None
  718. if method == 'link_vector':
  719. # Method: Use link vectors between pial and white surfaces as dipole orientations
  720. # Load downsampled pial and white surfaces
  721. pial_surf = nib.load(os.path.join(surf_dir, 'pial.ds.gii'))
  722. white_surf = nib.load(os.path.join(surf_dir, 'white.ds.gii'))
  724. # Extract vertices
  725. pial_vertices = pial_surf.darrays[0].data
  726. white_vertices = white_surf.darrays[0].data
  728. # Ensure same number of vertices in pial and white surfaces
  729. if pial_vertices.shape[0] != white_vertices.shape[0]:
  730. raise ValueError("Pial and white surfaces must have the same number of vertices")
  732. # Compute link vectors
  733. link_vectors = white_vertices - pial_vertices
  734. link_vectors = _normit(link_vectors)
  736. # Replicate link vectors for each layer
  737. orientations = np.tile(link_vectors, (len(layer_names), 1, 1))
  739. elif method == 'ds_surf_norm':
  740. # Method: Use normals of the downsampled surfaces
  741. orientations = []
  742. for l_idx, layer_name in enumerate(layer_names):
  743. if l_idx == 0 or not fixed:
  744. in_surf_path = os.path.join(surf_dir, f'{layer_name}.ds.gii')
  745. surf = nib.load(in_surf_path)
  746. vtx_norms, _ = mesh_normals(
  747. surf.darrays[0].data,
  748. surf.darrays[1].data,
  749. unit=True
  750. )
  751. orientations.append(vtx_norms)
  752. orientations = np.array(orientations)
  754. elif method == 'orig_surf_norm':
  755. # Method: Use normals of the original surfaces, mapped to downsampled surfaces
  756. orientations = []
  757. for l_idx, layer_name in enumerate(layer_names):
  758. if l_idx == 0 or not fixed:
  759. in_surf_path = os.path.join(surf_dir, f'{layer_name}.gii')
  760. orig_surf = nib.load(in_surf_path)
  761. ds_surf_path = os.path.join(surf_dir, f'{layer_name}.ds.gii')
  762. ds_surf = nib.load(ds_surf_path)
  763. kdtree = cKDTree(orig_surf.darrays[0].data)
  764. _, orig_vert_idx = kdtree.query(ds_surf.darrays[0].data, k=1)
  765. orig_vert_idx = np.squeeze(orig_vert_idx)
  767. vtx_norms, _ = mesh_normals(
  768. orig_surf.darrays[0].data,
  769. orig_surf.darrays[1].data,
  770. unit=True
  771. )
  772. orientations.append(vtx_norms[orig_vert_idx, :])
  773. orientations = np.array(orientations)
  775. elif method == 'cps':
  776. # Method: Use cortical patch statistics for normals
  777. orientations = []
  778. for l_idx, layer_name in enumerate(layer_names):
  779. if l_idx == 0 or not fixed:
  780. in_surf_path = os.path.join(surf_dir, f'{layer_name}.gii')
  781. orig_surf = nib.load(in_surf_path)
  782. ds_surf_path = os.path.join(surf_dir, f'{layer_name}.ds.gii')
  783. ds_surf = nib.load(ds_surf_path)
  784. kdtree = cKDTree(ds_surf.darrays[0].data)
  785. _, ds_vert_idx = kdtree.query(orig_surf.darrays[0].data, k=1)
  786. ds_vert_idx = np.squeeze(ds_vert_idx)
  788. orig_vtx_norms, _ = mesh_normals(
  789. orig_surf.darrays[0].data,
  790. orig_surf.darrays[1].data,
  791. unit=True
  792. )
  793. vtx_norms, _ = mesh_normals(
  794. ds_surf.darrays[0].data,
  795. ds_surf.darrays[1].data,
  796. unit=True
  797. )
  798. for v_idx in range(vtx_norms.shape[0]):
  799. orig_idxs = np.where(ds_vert_idx == v_idx)[0]
  800. if len(orig_idxs):
  801. vtx_norms[v_idx, :] = np.mean(orig_vtx_norms[orig_idxs, :], axis=0)
  802. vtx_norms = _normit(vtx_norms)
  803. orientations.append(vtx_norms)
  804. orientations = np.array(orientations)
  806. return orientations
  809. def create_layer_mesh(layer, hemispheres, fs_subject_dir) -> None:
  810. """
  811. Create or retrieve a specified cortical layer mesh file name or path based on the provided
  812. layer proportional thickness or identifier.
  814. Parameters
  815. ----------
  816. layer : float or int
  817. Specifies the cortical layer. The value 1 corresponds to the 'pial' surface, values between
  818. 0 and 1 (exclusive) correspond to intermediate layers (specified as a decimal), and the
  819. value 0 corresponds to the 'white matter' surface.
  820. hemispheres : list of str
  821. A list of hemisphere identifiers (e.g., ['lh', 'rh']) for which meshes should be created or
  822. retrieved.
  823. fs_subject_dir : str
  824. Path to the subject directory within the FreeSurfer environment. This directory should
  825. include a 'surf' directory where mesh files are stored.
  827. Returns
  828. -------
  829. layer_name : str or None
  830. Returns a string representing the mesh layer ('pial', 'white', or a specific intermediate
  831. layer as a formatted string). Returns None if the input layer does not match any recognized
  832. pattern (e.g., a negative number or a number greater than 1).
  834. Notes
  835. -----
  836. For intermediate layers (0 < layer < 1), the function will check for the existence of the mesh
  837. file corresponding to each hemisphere. If it does not exist, it uses 'mris_expand' to generate
  838. it using the white matter surface file. If the layer exactly matches 0 or 1, it returns the
  839. corresponding standard FreeSurfer mesh identifier ('white' or 'pial').
  840. """
  842. if layer == 1:
  843. return 'pial'
  844. if 0 < layer < 1:
  845. layer_name = f'{layer:.3f}'
  846. for hemi in hemispheres:
  847. wm_file = os.path.join(fs_subject_dir, 'surf', f'{hemi}.white')
  848. out_file = os.path.join(fs_subject_dir, 'surf', f'{hemi}.{layer_name}')
  849. if not os.path.exists(out_file):
  850. cmd = ['mris_expand', '-thickness', wm_file, f'{layer}', out_file]
  851. subprocess.run(cmd, check=True)
  852. return layer_name
  853. if layer == 0:
  854. return 'white'
  855. return None
  858. # pylint: disable=R0912,R0915
  859. def postprocess_freesurfer_surfaces(subj_id,
  860. out_dir,
  861. out_fname,
  862. n_surfaces=11,
  863. ds_factor=0.1,
  864. orientation='link_vector',
  865. fix_orientation=True,
  866. remove_deep=True,
  867. n_jobs=-1):
  868. """
  869. Process and combine FreeSurfer surface meshes for a subject.
  871. This function processes FreeSurfer surface meshes for a given subject by creating intermediate
  872. surfaces, adjusting for RAS offset, removing deep vertices, combining hemispheres,
  873. downsampling, and computing link vectors. The resulting surfaces are combined and saved to a
  874. specified output file.
  876. Parameters
  877. ----------
  878. subj_id : str
  879. Subject ID corresponding to the FreeSurfer subject directory.
  880. out_dir : str
  881. Output directory where the processed files will be saved.
  882. out_fname : str
  883. Filename for the final combined surface mesh.
  884. n_surfaces : int, optional
  885. Number of intermediate surfaces to create between white and pial surfaces.
  886. ds_factor : float, optional
  887. Downsampling factor for surface decimation.
  888. orientation : str, optional
  889. Method to compute orientation vectors ('link_vector' for pial-white link, 'ds_surf_norm'
  890. for downsampled surface normals, 'orig_surf_norm' for original surface normals, and 'cps'
  891. for cortical patch statistics).
  892. fix_orientation : bool, optional
  893. Flag to ensure that orientation of corresponding vertices across layers is the same (True
  894. by default).
  895. remove_deep : bool, optional
  896. Flag to remove vertices located in deep regions (labeled as 'unknown').
  897. n_jobs : int, optional
  898. Number of parallel processes to run. -1 for all available cores (default is -1).
  900. Notes
  901. -----
  902. - This function assumes the FreeSurfer 'SUBJECTS_DIR' environment variable is set.
  903. - Surfaces are processed in Gifti format and combined into a single surface mesh.
  905. Example
  906. -------
  907. >>> postprocess_freesurfer_surfaces('subject1', '/path/to/output', 'combined_surface.gii')
  908. """
  910. hemispheres = ['lh', 'rh']
  911. fs_subjects_dir = os.getenv('SUBJECTS_DIR')
  913. fs_subject_dir = os.path.join(fs_subjects_dir, subj_id)
  915. layers = np.linspace(1, 0, n_surfaces)
  917. ## Create intermediate surfaces if needed
  918. layer_names = Parallel(n_jobs=n_jobs, backend='loky')(
  919. delayed(create_layer_mesh)(layer, hemispheres, fs_subject_dir) for layer in layers
  920. )
  922. orig_mgz = os.path.join(fs_subject_dir, 'mri', 'orig.mgz')
  924. def _mat(flag):
  925. out = subprocess.check_output(['mri_info', flag, orig_mgz]).decode().split()
  926. return np.array([float(x) for x in out]).reshape(4, 4)
  928. t_orig = _mat('--vox2ras-tkr') # vox -> tkRAS (mm)
  929. n_orig = _mat('--vox2ras') # vox -> scanner RAS (mm)
  930. inv_t_orig = np.linalg.inv(t_orig)
  932. def tkras_to_scanner_ras(coords_mm):
  933. # coords_mm: (N,3) tkRAS
  934. xyz1 = np.c_[coords_mm, np.ones((coords_mm.shape[0], 1))]
  935. out = (n_orig @ (inv_t_orig @ xyz1.T)).T[:, :3]
  936. return out
  938. ## Convert to gifti, adjust for RAS offset, and remove deep vertices
  939. surfaces_to_process = copy.copy(layer_names)
  940. surfaces_to_process.append('inflated')
  942. # pylint: disable=R1702
  943. for surface_name in surfaces_to_process:
  944. for hemi in hemispheres:
  945. # Construct the original and new file names
  946. orig_name = os.path.join(fs_subject_dir, 'surf', f'{hemi}.{surface_name}')
  947. new_name = os.path.join(out_dir, f'{hemi}.{surface_name}.gii')
  948. rm_deep_name = os.path.join(out_dir, f'{hemi}.{surface_name}.nodeep.gii')
  950. # Convert the surface file to Gifti format
  951. subprocess.run(['mris_convert', orig_name, new_name], check=True)
  953. # Load the Gifti file
  954. surf_g = nib.load(new_name)
  956. # Set transformation matrix to identity
  957. surf_g.affine = np.eye(4)
  959. # Adjust for RAS offset
  960. n_vertices = 0
  961. for data_array in surf_g.darrays:
  962. if data_array.intent == nib.nifti1.intent_codes['NIFTI_INTENT_POINTSET']:
  963. data_array.data = tkras_to_scanner_ras(data_array.data)
  964. n_vertices = data_array.data.shape[0]
  965. nib.save(surf_g, new_name)
  967. annotation = os.path.join(fs_subject_dir, 'label', f'{hemi}.aparc.annot')
  968. label, _, names = nib.freesurfer.read_annot(annotation)
  970. # Remove vertices created by cutting the hemispheres
  971. if remove_deep:
  972. vertices_to_remove = []
  973. for vtx in range(n_vertices):
  974. region_name = names[label[vtx]]
  975. if isinstance(region_name, bytes):
  976. region_name = region_name.decode('utf-8', errors='ignore')
  977. if label[vtx] <= 0 or region_name.lower() == 'unknown':
  978. vertices_to_remove.append(vtx)
  979. surf_g = remove_vertices(surf_g, np.array(vertices_to_remove))
  981. # Save the modified Gifti file
  982. nib.save(surf_g, rm_deep_name)
  984. ## Combine hemispheres
  985. for surface_name in surfaces_to_process:
  986. # Load left and right hemisphere surfaces
  987. lh_fname = os.path.join(out_dir, f'lh.{surface_name}.nodeep.gii')
  988. l_hemi = nib.load(lh_fname)
  989. rh_fname = os.path.join(out_dir, f'rh.{surface_name}.nodeep.gii')
  990. r_hemi = nib.load(rh_fname)
  991. if surface_name == 'inflated':
  992. lh_width = np.max(l_hemi.darrays[0].data[:, 0])-np.min(l_hemi.darrays[0].data[:, 0])
  993. r_hemi.darrays[0].data[:, 0] = r_hemi.darrays[0].data[:, 0] + \
  994. np.max(l_hemi.darrays[0].data[:, 0]) + \
  995. (.5*lh_width)
  997. # Combine the surfaces
  998. combined = combine_surfaces([l_hemi, r_hemi])
  999. combined_fname = os.path.join(out_dir, f'{surface_name}.gii')
  1000. nib.save(combined, combined_fname)
  1002. ## Downsample surfaces at the same time
  1003. # Get list of surfaces
  1004. in_surfs = []
  1005. for surface_name in surfaces_to_process:
  1006. in_surf_fname = os.path.join(out_dir, f'{surface_name}.gii')
  1007. in_surf = nib.load(in_surf_fname)
  1008. in_surfs.append(in_surf)
  1010. # Downsample multiple surfaces
  1011. out_surfs = downsample_multiple_surfaces(in_surfs, ds_factor)
  1012. for surface_name, out_surf in zip(surfaces_to_process, out_surfs):
  1013. out_surf_path = os.path.join(out_dir, f'{surface_name}.ds.gii')
  1014. nib.save(out_surf, out_surf_path)
  1016. ## Compute dipole orientations
  1017. orientations = compute_dipole_orientations(
  1018. orientation,
  1019. layer_names,
  1020. out_dir,
  1021. fixed=fix_orientation
  1022. )
  1024. base_fname = f'ds.{orientation}.not_fixed'
  1025. if fix_orientation:
  1026. base_fname = f'ds.{orientation}.fixed'
  1027. for l_idx, layer_name in enumerate(layer_names):
  1028. in_surf_path = os.path.join(out_dir, f'{layer_name}.ds.gii')
  1029. surf = nib.load(in_surf_path)
  1031. # Set these link vectors as the normals for the downsampled surface
  1032. ori_array=nib.gifti.GiftiDataArray(data=orientations[l_idx, :, :],
  1033. intent=nib.nifti1.intent_codes['NIFTI_INTENT_VECTOR'])
  1034. surf.add_gifti_data_array(ori_array)
  1036. # Save the modified downsampled surface with link vectors as normals
  1037. out_surf_path = os.path.join(out_dir, f'{layer_name}.{base_fname}.gii')
  1038. nib.save(surf, out_surf_path)
  1040. ## Combine layers
  1041. all_surfs = []
  1042. for layer_name in layer_names:
  1043. surf_path = os.path.join(out_dir, f'{layer_name}.{base_fname}.gii')
  1044. surf = nib.load(surf_path)
  1045. all_surfs.append(surf)
  1047. combined = combine_surfaces(all_surfs)
  1048. nib.save(combined, os.path.join(out_dir, out_fname))
  1051. def mesh_adjacency(faces):
  1052. """
  1053. Compute the adjacency matrix of a triangle mesh.
  1055. Parameters
  1056. ----------
  1057. faces : np.array, shape (f, 3)
  1058. The mesh faces, where `f` is the number of faces, and each face is represented by a tuple of
  1059. three vertex indices.
  1061. Returns
  1062. -------
  1063. adjacency : np.array, shape (v, v)
  1064. The adjacency matrix of the mesh, where `v` is the number of vertices. Each entry (i, j)
  1065. indicates whether vertices i and j are connected by an edge.
  1066. """
  1068. faces = np.asarray(faces, dtype=int)
  1069. n_vertices = np.max(faces)+1 # Assuming max vertex index represents the number of vertices
  1071. # Flatten the indices to create row and column indices for the adjacency matrix
  1072. row_indices = np.hstack([faces[:, 0], faces[:, 0], faces[:, 1],
  1073. faces[:, 1], faces[:, 2], faces[:, 2]])
  1074. col_indices = np.hstack([faces[:, 1], faces[:, 2], faces[:, 0],
  1075. faces[:, 2], faces[:, 0], faces[:, 1]])
  1077. # Create a sparse matrix from row and column indices
  1078. adjacency = csr_matrix(
  1079. (np.ones_like(row_indices), (row_indices, col_indices)),
  1080. shape=(n_vertices, n_vertices)
  1081. )
  1083. # Ensure the adjacency matrix is binary
  1084. adjacency = (adjacency > 0).astype(int)
  1086. return adjacency
  1089. def interpolate_data(original_mesh, downsampled_mesh, downsampled_data, adjacency_matrix=None,
  1090. k_neighbors=5):
  1091. """
  1092. Interpolate vertex data from a downsampled mesh back to the original mesh using weighted
  1093. averaging of the k-nearest downsampled vertices. After initial interpolation, vertices that
  1094. are in the downsampled mesh are refined using a weighted sum of their neighbors.
  1096. Parameters
  1097. ----------
  1098. original_mesh : nibabel.gifti.GiftiImage
  1099. The original high-resolution mesh as a nibabel Gifti object.
  1100. downsampled_mesh : nibabel.gifti.GiftiImage
  1101. The downsampled version of the original mesh.
  1102. downsampled_data : np.ndarray
  1103. Data values at each vertex in 'downsampled_mesh'.
  1104. adjacency_matrix : scipy.sparse matrix, optional
  1105. The adjacency matrix of the original mesh.
  1106. If None, it will be computed.
  1107. k_neighbors : int, optional
  1108. Number of nearest downsampled neighbors to use for interpolation.
  1110. Returns
  1111. -------
  1112. vertex_data : np.ndarray
  1113. An array of interpolated data for each vertex in 'original_mesh'.
  1114. """
  1115. original_vertices = original_mesh.darrays[0].data
  1116. downsampled_vertices = downsampled_mesh.darrays[0].data
  1118. # Build a KD-tree for the downsampled mesh
  1119. tree = cKDTree(downsampled_vertices)
  1121. # Find the k nearest downsampled vertices for each original vertex
  1122. distances, indices = tree.query(original_vertices, k=k_neighbors)
  1123. indices = np.squeeze(indices)
  1125. # Initialize interpolated data array
  1126. vertex_data = np.full(len(original_vertices), np.nan)
  1128. # Compute weighted interpolation from k-nearest downsampled vertices
  1129. valid_mask = indices[:, 0] < len(downsampled_data) # Ensure indices are within bounds
  1130. for i in np.where(valid_mask)[0]:
  1131. valid_indices = indices[i, :]
  1132. valid_distances = distances[i, :]
  1134. # Avoid division by zero (replace zero distances with a small number)
  1135. valid_distances[valid_distances == 0] = 1e-6
  1136. weights = 1 / valid_distances
  1137. weights /= weights.sum() # Normalize weights
  1139. vertex_data[i] = np.dot(weights, downsampled_data[valid_indices])
  1141. # Second pass: Refine values for vertices that were in the downsampled mesh
  1142. if adjacency_matrix is None:
  1143. adjacency_matrix = mesh_adjacency(original_mesh.darrays[1].data)
  1145. downsampled_vertex_mask = np.isin(np.arange(len(original_vertices)), indices[:, 0])
  1147. for i in np.where(downsampled_vertex_mask)[0]:
  1148. neighbors = adjacency_matrix[i].nonzero()[1]
  1149. if len(neighbors) > 0:
  1150. distances = np.linalg.norm(original_vertices[neighbors] - original_vertices[i], axis=1)
  1151. distances[distances == 0] = 1e-6 # Avoid division by zero
  1152. weights = 1 / distances
  1153. weights /= weights.sum() # Normalize weights
  1155. vertex_data[i] = np.dot(weights, vertex_data[neighbors])
  1157. return vertex_data
  1160. def split_fv(faces, vertices):
  1161. """
  1162. Split faces and vertices into connected pieces based on the connectivity of the faces.
  1164. Parameters
  1165. ----------
  1166. faces : np.ndarray
  1167. A 2D numpy array of faces, where each row represents a face and each element is an index to
  1168. a vertex in `vertices`.
  1169. vertices : np.ndarray
  1170. A 2D numpy array of vertices, where each row represents a vertex.
  1172. Returns
  1173. -------
  1174. fv_out : list of dict
  1175. A list where each element is a dictionary with keys 'faces' and 'vertices'. Each dictionary
  1176. represents a separately connected patch of the mesh.
  1178. Examples
  1179. --------
  1180. >>> faces = np.array([[1, 2, 3], [1, 3, 4], [5, 6, 1], [7, 8, 9], [11, 10, 4]])
  1181. >>> vertices = np.array([[2, 4], [2, 8], [8, 4], [8, 0], [0, 4], [2, 6], [2, 2], [4, 2],
  1182. >>> [4, 0], [5, 2], [5, 0]])
  1183. >>> split_patches = split_fv(faces, vertices)
  1185. Notes
  1186. -----
  1187. Faces and vertices should be defined such that faces sharing a vertex reference the same vertex
  1188. number. This function does not explicitly test for duplicate vertices at the same location.
  1189. """
  1191. num_faces = faces.shape[0]
  1192. f_sets = np.zeros(num_faces, dtype=np.uint32)
  1193. current_set = 0
  1195. while np.any(f_sets == 0):
  1196. current_set += 1
  1197. next_avail_face = np.where(f_sets == 0)[0][0]
  1198. open_vertices = faces[next_avail_face]
  1200. while open_vertices.size:
  1201. avail_face_inds = np.where(f_sets == 0)[0]
  1202. is_member = np.isin(faces[avail_face_inds], open_vertices)
  1203. avail_face_sub = np.where(np.any(is_member, axis=1))[0]
  1204. f_sets[avail_face_inds[avail_face_sub]] = current_set
  1205. open_vertices = np.unique(faces[avail_face_inds[avail_face_sub]])
  1207. fv_out = []
  1208. for set_num in range(1, current_set + 1):
  1209. set_f = faces[f_sets == set_num]
  1210. unique_vertices, new_vertex_indices = np.unique(set_f, return_inverse=True)
  1211. fv_out.append({
  1212. 'faces': new_vertex_indices.reshape(set_f.shape),
  1213. 'vertices': vertices[unique_vertices]
  1214. })
  1216. return fv_out

surf.py, under GPL-3.0 · at the source

Overview

Authors: Maciej J. Szul1,2,3, Matteo G. Maspoli1,2,4,5, Danila Shelepenkov2,4,5, Ishita Agarwal1,2,6, Quentin Moreau1,2, Solène Gailhard1,2, Maxime Ferez2,4,5, Carolina Fernandez Pujol7,8, Yunkai Zhu8,9, Bassem Hiba1,2, Sebastien Daligault10, Franck Lamberton10,11, Denis Schwartz4,5,10, Jérémie Mattout2,4,5, Mathilde Bonnefond2,4,5, Andrew R. Dykstra8,9, Sven Bestmann12,13, Gareth R. Barnes12, James J. Bonaiuto1,2
ORCID iDs: Maciej J. Szul, Matteo G. Maspoli, Maxime Ferez, Sebastien Daligault, Denis Schwartz, Jérémie Mattout, Mathilde Bonnefond, Andrew R. Dykstra, Gareth R. Barnes, James J. Bonaiuto
13 affiliations
  1. Decision, Action, and Neural Computation (DANC) team, Institut des Sciences Cognitives Marc Jeannerod,CNRS UMR5229 Lyon, France
  2. Université Claude Bernard Lyon 1, Université de Lyon,Lyon, France
  3. Department of Psychiatry and Psychotherapy, Faculty of Medicine, University of Tübingen,Tübingen, Germany
  4. Lyon Neuroscience Research Center, CRNL, INSERM U1028, CNRS UMR5292,Lyon, France
  5. INRIA, Computation, Cognition and Neurophysiology (COPHY) team,Lyon, France
  6. Department of Psychological Sciences, Purdue University,West Lafayette, IN USA
  7. Department of Neuroscience, Carney Institute of Brain Science, Brown University,Providence, RI USA
  8. Department of Biomedical Engineering, University of Miami,Coral Gables, FL USA
  9. School of Communication Sciences and Disorders, University of Central Florida,Orlando, FL USA
  10. CERMEP-Imagerie du Vivant, Lyon, France
  11. Lyon-East Health Research Federation, CNRS UAR3453, INSERM US7, Université Claude Bernard Lyon 1,Lyon, France
  12. Department of Imaging Neuroscience, UCL Queen Square Institute of Neurology, University College London,London, UK
  13. Department of Clinical and Movement Neurosciences, UCL Queen Square Institute of Neurology, University College London,London, UK
Institutions: Institut des Sciences Cognitives Marc Jeannerod (France); Université de Lyon (France); University of Tübingen (Germany); Centre de Recherche en Neurosciences de Lyon (France); Calcul, Cognition et Neurophysiologie (France); Purdue University West Lafayette (United States); Brown University (United States); University of Miami (United States); University of Central Florida (United States); Centre d'Exploration et de Recherche Médicale par Emission de Positons (France); Lyon 1 Université (France); University College London (United Kingdom)
Journal: Nature communications, volume 17, issue 1, article 10128
Dates: received 20 March 2026; accepted 31 August 2026; published online 30 September 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41467-026-77782-9 · PMID 42816533 · PMCID PMC13627721 · OpenAlex W7214994117
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: MEG (modality), human (organism), systems (subfield)
Methods: Connectivity, Statistics, Machine learning, Preprocessing, Spectral & time-frequency, fMRI & imaging
Keywords: Neural circuits, Neuroscience
MeSH: Magnetoencephalography*, Models, Neurological*, Adult, Brain, Brain Mapping, Cerebral Cortex, Computer Simulation, Female, Humans, Male, Signal-To-Noise Ratio (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: European Research Council (864550, 716862); Agence Nationale de la Recherche (ANR-20-CE17-0023, ANR-25-CE37-0908-01); Fondation pour l’Audition (RD-2022-03); Wellcome Trust (Wellcome) (226793/Z/22/Z)
Citations: not cited yet (Europe PMC); 69 references in the paper

Abstract

Neural dynamics at the laminar level are critical for cortical computation. However, in humans, non-invasive methods to probe such dynamics have been limited to coarse distinctions between deep and superficial layers. Here, we present a multilayer magnetoencephalography source reconstruction framework and evaluate the conditions under which depth-resolved laminar inference may be feasible. Using simulations, we systematically assess the limits of magnetoencephalography depth resolution, showing that laminar discrimination depends on sufficiently high signal-to-noise ratio, precise co-registration, and accurate specification of cortical column orientation. We demonstrate that regional variations in cortical anatomy influence reconstruction fidelity, with lead-field separability emerging as a key determinant. We then apply this framework to empirical data from three independent datasets and find laminar activation patterns that align with canonical feedforward and feedback motifs in visual and sensorimotor circuits, supporting the plausibility of laminar inference under favorable conditions and offering opportunities to bridge invasive electrophysiology and human neuroimaging.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 16 matches between paragraphs and lines of code.

danclab/spm

License: GPL-2.0
State: the link answers, verified on 3 October 2026
Evidence: files inventoried
Commit: 2a3ea29c41014b0230fb78ff34401ffbe68fb1ab, 29 September 2026
Languages: MATLAB (4597), C (126), C/C++ (34), CUDA (13), C++ (5), JavaScript (2), Java (1)
Size: 5,433 files, 4,778 scripts
Software Heritage: archived
Found in: the text, “Simulations”
Holds: README, license file, continuous integration
Not found: CITATION.cff, environment file, tests, documentation
Tools: FieldTrip (1,004 files), SPM (190 files), Signal Processing Toolbox (20 files), GIfTI library for MATLAB (13 files), Statistics and Machine Learning Toolbox (11 files), FreeSurfer (6 files), Image Processing Toolbox (6 files), Optimization Toolbox (3 files), Neurodata Without Borders (PyNWB, MatNWB) (3 files), Brain Connectivity Toolbox (1 file), Brainstorm (1 file), EEGLAB (1 file)
Availability: 1 check, the latest on 3 October 2026: the link answers
  • 3 October 2026: the link answers
2,000 files
At the source: github.com/danclab/spm

danclab/multilaminar_sim

License: none: the authors keep all their rights
State: the link answers, verified on 3 October 2026
Evidence: files inventoried
Commit: cdc682481d6afc522ac3c8de534f8b5eb394261b, 4 August 2026
Languages: Jupyter (31), Python (23), Shell (15), MATLAB (5)
Size: 78 files, 74 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, 31 notebooks
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (9 files), Matplotlib (7 files), NiBabel (7 files), pandas (7 files), SciPy (7 files), Pillow (3 files), h5py (1 file), scikit-learn (1 file), statsmodels (1 file)
Availability: 1 check, the latest on 3 October 2026: the link answers
  • 3 October 2026: the link answers
10 files
At the source: github.com/danclab/multilaminar_sim/

jbonaiuto/cued_action_selection

License: GPL-3.0
State: the link answers, verified on 3 October 2026
Evidence: files inventoried
Commit: 58ea009b36a03992609bf68bf493158972ac8102, 14 September 2022
Languages: MATLAB (15)
Size: 17 files, 15 scripts
Software Heritage: archived
Found in: “Code availability”
Holds: README, license file
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Psychtoolbox (1 file)
Availability: 1 check, the latest on 3 October 2026: the link answers
  • 3 October 2026: the link answers
17 files
At the source: github.com/jbonaiuto/cued_action_selection

danclab/auditory_laminar

License: none: the authors keep all their rights
State: the link answers, verified on 3 October 2026
Evidence: files inventoried
Commit: ecf322b04d855c143ffb2d86b160d11230e3cd3c, 15 June 2026
Languages: Python (24), Jupyter (2), MATLAB (2)
Size: 385 files, 28 scripts
Software Heritage: archived
Found in: “Code availability”
Holds: 2 notebooks
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (14 files), MNE-Python (12 files), Matplotlib (11 files), pandas (7 files), SciPy (7 files), NiBabel (5 files), PsychoPy (5 files), Psychtoolbox (5 files), autoreject (2 files), h5py (2 files), FreeSurfer (1 file), FSL (1 file), MEEGkit (1 file)
Availability: 1 check, the latest on 3 October 2026: the link answers
  • 3 October 2026: the link answers
28 files
At the source: github.com/danclab/auditory_laminar

TommyClausner/meglyon-experiment-runner

License: GPL-3.0
State: the link answers, verified on 3 October 2026
Evidence: files inventoried
Commit: 3142259d00c028187537b91c54fe6151f49f31fc, 25 January 2024
Languages: Python (17)
Size: 2,089 files, 17 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (requirements.txt, setup.cfg, setup.py)
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: NumPy (13 files), PsychoPy (7 files), Matplotlib (2 files)
Availability: 1 check, the latest on 3 October 2026: the link answers
  • 3 October 2026: the link answers
19 files
At the source: gitlab.com/TommyClausner/meglyon-experiment-runner

danclab/laminar_motor_erf

License: GPL-3.0
State: the link answers, verified on 3 October 2026
Evidence: files inventoried
Commit: c4a2f09d40615349d710ed32e73ea22668ef26a1, 15 June 2026
Languages: Python (16), Jupyter (2), MATLAB (1)
Size: 22 files, 19 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: license file, 2 notebooks
Not found: README, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (12 files), MNE-Python (10 files), Matplotlib (7 files), SciPy (7 files), NiBabel (5 files), pandas (4 files), h5py (2 files), autoreject (1 file), MEEGkit (1 file)
Availability: 1 check, the latest on 3 October 2026: the link answers
  • 3 October 2026: the link answers
20 files
At the source: github.com/danclab/laminar_motor_erf/

cophyteam/laminar_visual_erf

License: GPL-3.0
State: the link answers, verified on 3 October 2026
Evidence: files inventoried
Commit: 7361a1c56cbf0ac1edf8f5fb5ffa2f2f0d328fc0, 5 August 2026
Languages: Jupyter (9), Python (2), MATLAB (2)
Size: 17 files, 13 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 9 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (5 files), NumPy (5 files), SciPy (4 files), MNE-Python (3 files), NiBabel (3 files), pandas (3 files), autoreject (1 file), FreeSurfer (1 file), h5py (1 file)
Availability: 1 check, the latest on 3 October 2026: the link answers
  • 3 October 2026: the link answers
8 files
At the source: github.com/cophyteam/laminar_visual_erf

Zenodo 21638150

License: GPL-3.0
State: the link answers, verified on 3 October 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (20 files), NiBabel (13 files), Matplotlib (10 files), SciPy (7 files), h5py (4 files), MNE-Python (3 files), Elephant (2 files), FreeSurfer (2 files), Neo (2 files)
Availability: 1 check, the latest on 3 October 2026: the link answers (HTTP 200)
  • 3 October 2026: the link answers (HTTP 200)
27 files
At the source:

Code availability

The laMEG software package used for multilayer source reconstruction and laminar inference has been deposited in Zenodo65. Experimental stimulus presentation and behavioral response collection were performed using task-specific software available at: https://github.com/jbonaiuto/cued_action_selection, https://github.com/danclab/auditory_laminar, and https://gitlab.com/TommyClausner/meglyon-experiment-runner. All custom data analysis code used in this study is available at https://github.com/danclab/multilaminar_sim/, https://github.com/danclab/laminar_motor_erf/, https://github.com/danclab/auditory_laminar, and https://github.com/cophyteam/laminar_visual_erf.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 8 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 2,119 scripts, each with its path and the digest of its content;
  • 16 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data availability

The preprocessed MEG data, structural MRI scans, cortical surface reconstructions, and fiducial locations underlying the empirical analyses have been deposited in the Zenodo repository (https://doi.org/10.5281/zenodo.20721185). Raw MEG and MRI data are not publicly available because they contain identifiable human participant data and are subject to institutional ethical and data protection requirements. Access to the raw data may be granted to qualified researchers for research purposes, subject to approval by the relevant ethics committee and completion of an appropriate data-sharing agreement. Requests should be directed to the corresponding author, who will respond within a reasonable timeframe. Access will be provided for the duration of the approved research project. Source data are provided with this paper.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 3 October 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 19 authors, 2 keywords, 11 MeSH terms, 4 funders, 66 references.

Cite

This paper

Szul, M. J., Maspoli, M. G., Shelepenkov, D., Agarwal, I., Moreau, Q., Gailhard, S., Ferez, M., Fernandez Pujol, C., Zhu, Y., Hiba, B., Daligault, S., Lamberton, F., Schwartz, D., Mattout, J., Bonnefond, M., Dykstra, A. R., Bestmann, S., Barnes, G. R., & Bonaiuto, J. J. (2026). Multilayer MEG source modelling enables depth-resolved laminar inference in humans. Nature communications, 17(1), 10128. https://doi.org/10.1038/s41467-026-77782-9

BibTeX

@article{szul2026multilayer,
author = {Szul, Maciej J. and Maspoli, Matteo G. and Shelepenkov, Danila and Agarwal, Ishita and Moreau, Quentin and Gailhard, Solène and Ferez, Maxime and Fernandez Pujol, Carolina and Zhu, Yunkai and Hiba, Bassem and Daligault, Sebastien and Lamberton, Franck and Schwartz, Denis and Mattout, Jérémie and Bonnefond, Mathilde and Dykstra, Andrew R. and Bestmann, Sven and Barnes, Gareth R. and Bonaiuto, James J.},
title = {{Multilayer MEG source modelling enables depth-resolved laminar inference in humans}},
journal = {Nature communications},
year = {2026},
month = sep,
volume = {17},
number = {1},
pages = {10128},
publisher = {Nature Publishing Group},
issn = {2041-1723},
doi = {10.1038/s41467-026-77782-9},
url = {https://doi.org/10.1038/s41467-026-77782-9},
pmid = {42816533},
pmcid = {PMC13627721}
}

RIS

TY - JOUR
AU - Szul, Maciej J.
AU - Maspoli, Matteo G.
AU - Shelepenkov, Danila
AU - Agarwal, Ishita
AU - Moreau, Quentin
AU - Gailhard, Solène
AU - Ferez, Maxime
AU - Fernandez Pujol, Carolina
AU - Zhu, Yunkai
AU - Hiba, Bassem
AU - Daligault, Sebastien
AU - Lamberton, Franck
AU - Schwartz, Denis
AU - Mattout, Jérémie
AU - Bonnefond, Mathilde
AU - Dykstra, Andrew R.
AU - Bestmann, Sven
AU - Barnes, Gareth R.
AU - Bonaiuto, James J.
TI - Multilayer MEG source modelling enables depth-resolved laminar inference in humans
T2 - Nature communications
J2 - Nat Commun
PY - 2026
DA - 2026/09/30
VL - 17
IS - 1
SP - 10128
SN - 2041-1723
PB - Nature Publishing Group
DO - 10.1038/s41467-026-77782-9
UR - https://doi.org/10.1038/s41467-026-77782-9
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41467-026-77782-9",
"type": "article-journal",
"title": "Multilayer MEG source modelling enables depth-resolved laminar inference in humans",
"container-title": "Nature communications",
"author": [
{
"family": "Szul",
"given": "Maciej J."
},
{
"family": "Maspoli",
"given": "Matteo G."
},
{
"family": "Shelepenkov",
"given": "Danila"
},
{
"family": "Agarwal",
"given": "Ishita"
},
{
"family": "Moreau",
"given": "Quentin"
},
{
"family": "Gailhard",
"given": "Solène"
},
{
"family": "Ferez",
"given": "Maxime"
},
{
"family": "Fernandez Pujol",
"given": "Carolina"
},
{
"family": "Zhu",
"given": "Yunkai"
},
{
"family": "Hiba",
"given": "Bassem"
},
{
"family": "Daligault",
"given": "Sebastien"
},
{
"family": "Lamberton",
"given": "Franck"
},
{
"family": "Schwartz",
"given": "Denis"
},
{
"family": "Mattout",
"given": "Jérémie"
},
{
"family": "Bonnefond",
"given": "Mathilde"
},
{
"family": "Dykstra",
"given": "Andrew R."
},
{
"family": "Bestmann",
"given": "Sven"
},
{
"family": "Barnes",
"given": "Gareth R."
},
{
"family": "Bonaiuto",
"given": "James J."
}
],
"container-title-short": "Nat Commun",
"volume": "17",
"issue": "1",
"page": "10128",
"DOI": "10.1038/s41467-026-77782-9",
"PMID": "42816533",
"PMCID": "PMC13627721",
"ISSN": "2041-1723",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41467-026-77782-9",
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
30
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1101/2025.05.28.656642 [code]
Multilayer MEG source modelling enables depth-resolved laminar inference in humans
Journal: bioRxiv (preprint)
In common: MEEGkit, autoreject, Brainstorm, 24 other tools, MEG, systems, 38 references, 6 authors
[2] doi:10.1162/imag_a_00391 [code]
Surfing beta burst waveforms to improve motor imagery-based BCI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: MEEGkit, autoreject, MNE-Python, 5 other tools, 2 references, 3 authors
[3] doi:10.64898/2026.03.06.710026 [code]
Distinct beta burst motifs exhibit opposing error relationships during motor adaptation
Journal: bioRxiv (preprint)
In common: MNE-Python, FieldTrip, statsmodels, 7 other tools, 2 references, 3 authors
[4] doi:10.1523/jneurosci.0886-23.2023 [code]
Bursting with Potential: How Sensorimotor Beta Bursts Develop from Infancy to Adulthood.
Journal: The Journal of neuroscience : the official journal of the Society for Neuroscience
In common: MEEGkit, EEGLAB, MNE-Python, 9 other tools, 2 references, 2 authors
[5] doi:10.7554/elife.108408 [code]
Frequency and laminar profile of feature-specific visual activity revealed by interleaved EEG-fMRI.
Journal: eLife
In common: FieldTrip, FreeSurfer, FSL, 11 other tools, systems, 5 references, author Mathilde Bonnefond
[6] doi:10.1162/imag.a.1189 [code]
Distorting anatomy to test MEG models and metrics.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: GIfTI library for MATLAB, SPM, Statistics and Machine Learning Toolbox, MEG, 10 references, author Gareth R. Barnes
[7] doi:10.1162/imag.a.1040 [code]
Novel 4 He-OPMs support waveform-specific beta burst analysis comparable to SQUID-MEG
Journal: n/a
In common: MEEGkit, autoreject, MNE-Python, 8 other tools, MEG, 2 references, author Denis Schwartz
[8] doi:10.1038/s41467-023-41689-6 [code]
Neurophysiological signatures of cortical micro-architecture.
Journal: Nature communications
In common: Brainstorm, GIfTI library for MATLAB, FieldTrip, 13 other tools, MEG, 2 references
[9] doi:10.1523/eneuro.0288-25.2025 [code]
Visual Speech Reduces Cognitive Effort as Measured by EEG Theta Power and Pupil Dilation.
Journal: eNeuro
In common: Brainstorm, Neurodata Without Borders (PyNWB, MatNWB), Psychtoolbox, 12 other tools
[10] doi:10.1162/imag_a_00410 [code]
Inferring laminar origins of MEG signals with optically pumped magnetometers (OPMs): A simulation study.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: GIfTI library for MATLAB, FieldTrip, SPM, 2 other tools, MEG, systems, 9 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1038/s41467-026-77782-9.svg)](https://openscicode.org/badge/doi/10.1038/s41467-026-77782-9)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1038/s41467-026-77782-9.svg
   :target: https://openscicode.org/badge/doi/10.1038/s41467-026-77782-9
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1038/s41467-026-77782-9"><img src="https://openscicode.org/badge/doi/10.1038/s41467-026-77782-9.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.