OSCR

Siibra: a software tool suite for realizing a Multilevel Human Brain Atlas from complex data resources.

Code ↔ Paper

25 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 25 matches
  1. [1] § Results › Multimodal comparison of brain areas ↔ examples/tutorials/2025-paper-fig5.py, lines 17–38 · score 0.98 · high resolution scans, underlying cloud resource, chosen position, oriented cortical, motor area, full resolution image
  2. [2] § Results › Anatomical evaluation of subcortical maps ↔ examples/tutorials/2025-paper-fig6.py, lines 16–49 · score 0.84 · ventral intermediate nucleus, subcortical maps, structural connectivity, thalamus, diffusion, Precentral
  3. [3] § Results › Multimodal comparison of brain areas ↔ siibra/livequeries/bigbrain.py, lines 253–394 · score 0.79 · intersected layer IV, cortical image patch, cortical layer surface, BigBrain, closest, oriented
  4. [4] § Methods › Memory-efficient representation of sparse regional maps ↔ siibra/volumes/parcellationmap.py, lines 1142–1283 · score 0.77 · Statistical maps, region maps, image volume, parcellation maps, coordinate space, numerical
  5. [5] § Results › Anatomically guided reproducible extraction of microscopy data ↔ examples/tutorials/2025-paper-fig5.py, lines 17–38 · score 0.74 · high resolution image, siibra allows, full resolution, BigBrain, precomputed, microscopy
  6. [6] § Results › Design and architecture of siibra ↔ backend/app/bkwdcompat.py, lines 104–173 · score 0.74 · DiFuMo, MNI Colin, cytoarchitectonic maps, cortical layer, Julich Brain, bundle
  7. [7] § Methods › Support of FAIR principles ↔ siibra/volumes/sparsemap.py, lines 489–545 · score 0.72 · NumPy, Nifti1Image, nibabel, nilearn, fetching, siibra
  8. [8] § Methods › Modeling spatial entities in different reference coordinate systems ↔ siibra/livequeries/bigbrain.py, lines 253–394 · score 0.69 · cortical image patch, point clouds, BigBrain, vertices, probability, locations
  9. [9] § Methods › Memory-efficient representation of sparse regional maps ↔ siibra/volumes/sparsemap.py, lines 338–468 · score 0.67 · bounding boxes, Statistical maps, image volume, sparse, weights, probabilistic
  10. [10] § Methods › Separation of content from code › Live queries ↔ siibra/livequeries/ebrains.py, lines 34–145 · score 0.64 · EBRAINS Knowledge Graph, anatomically anchored, live queries, siibra, Atlas
  11. [11] § Results › Anatomical characterization and multimodal profiling for regions of interest ↔ siibra/volumes/parcellationmap.py, lines 1142–1283 · score 0.63 · statistical maps, image volume, brain regions, reference space, correlation, assignment
  12. [12] § Methods › Modular software architecture ↔ siibra/core/space.py, lines 26–143 · score 0.61 · inflated surface, FreeSurfer, pial, volumetric, space, siibra
  13. [13] § Methods › Modeling spatial entities in different reference coordinate systems ↔ hbp_spatial_backend/__init__.py, lines 150–224 · score 0.61 · folding patterns, transform coordinates, alignment, diffeomorphisms, Brain
  14. [14] § Results › Visually guided exploration from full brain networks to cells ↔ backend/app/bkwdcompat.py, lines 104–173 · score 0.61 · MNI Colin, cytoarchitectonic maps, cortical layer, Julich Brain, parcellation, space
  15. [15] § Results › Expansion with new contents ↔ backend/app/bkwdcompat.py, lines 69–102 · score 0.61 · Allen Mouse Brain, Waxholm Space, Rat Brain, Brain Atlas, parcellation
  16. [16] § Results › Anatomically guided reproducible extraction of microscopy data ↔ siibra/features/image/sections.py, lines 38–107 · score 0.59 · cortical patch, bounding box, BigBrain, fetch, resolutions, volume
  17. [17] § Results › Anatomical characterization and multimodal profiling for regions of interest ↔ siibra/volumes/parcellationmap.py, lines 1285–1383 · score 0.59 · Gaussian blob, uncertain coordinates, kernel, split, assignment, volume
  18. [18] § Results › Anatomical characterization and multimodal profiling for regions of interest ↔ siibra/livequeries/allen.py, lines 69–120 · score 0.59 · Allen Human Brain, gene expression, Human Brain Atlas, microarray, tissue, connectivity
  19. [19] § Results › Anatomical characterization and multimodal profiling for regions of interest ↔ siibra/features/tabular/gene_expression.py, lines 30–170 · score 0.58 · gene expression, molecular, Human Brain Atlas, tabular, microarray, Allen
  20. [20] § Methods › Linking data features to atlas elements ↔ siibra/features/tabular/tabular.py, lines 29–157 · score 0.57 · neurotransmitter receptor, tissue samples, properties, locations, anatomical, brain
  21. [21] § Methods › Separation of content from code › Live queries ↔ siibra/livequeries/allen.py, lines 69–120 · score 0.56 · Allen Human Brain, live queries, microarray, API, interfaces, connecting
  22. [22] § Methods › Linking data features to atlas elements ↔ examples/tutorials/2025-paper-fig6.py, lines 16–49 · score 0.55 · inter subject variability, reliable, histological, anatomical, maps
  23. [23] § Results › Anatomical characterization and multimodal profiling for regions of interest ↔ siibra/features/tabular/gene_expression.py, lines 30–170 · score 0.53 · gene expressions, Human Brain Atlas, tabular, Allen, co, retrieved
  24. [24] § Results › Multimodal comparison of brain areas ↔ siibra/livequeries/allen.py, lines 142–200 · score 0.52 · microarray probe, tissue samples, donors, Allen, genes, MNI
  25. [25] § Methods › Modular software architecture ↔ src/viewerModule/nehuba/nehubaViewerGlue/nehubaViewerGlue.component.spec.ts, lines 1–66 · score 0.51 · user interactions, user interface, touch, reactive, Angular, viewer

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,859 lines · 73 KB · Apache-2.0 · 3 matches

  1. # Copyright 2018-2026
  2. # Institute of Neuroscience and Medicine (INM-1), Forschungszentrum Jülich GmbH
  3. # Licensed under the Apache License, Version 2.0 (the "License");
  4. # you may not use this file except in compliance with the License.
  5. # You may obtain a copy of the License at
  6. # http://www.apache.org/licenses/LICENSE-2.0
  7. # Unless required by applicable law or agreed to in writing, software
  8. # distributed under the License is distributed on an "AS IS" BASIS,
  9. # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  10. # See the License for the specific language governing permissions and
  11. # limitations under the License.
  12. """Provides spatial representations for parcellations and regions."""
  13. from collections import defaultdict
  14. from dataclasses import dataclass, asdict
  15. from typing import Union, Dict, List, TYPE_CHECKING, Iterable, Tuple, Literal, NamedTuple
  16. import numpy as np
  17. import pandas as pd
  18. from scipy.ndimage import distance_transform_edt
  19. from . import volume as _volume
  20. from .providers import provider
  21. from .. import exceptions
  22. from ..commons import (
  23. MapIndex,
  24. MapType,
  25. compare_arrays,
  26. resample_img_to_img,
  27. connected_components,
  28. clear_name,
  29. create_key,
  30. create_gaussian_kernel,
  31. siibra_tqdm,
  32. Species,
  33. CompareMapsResult,
  34. generate_uuid,
  35. logger,
  36. QUIET,
  37. )
  38. from ..core import concept, space, parcellation, region as _region
  39. from ..locations import location, point, pointcloud
  40. if TYPE_CHECKING:
  41. from ..core.region import Region
  42. from nilearn.maskers import NiftiLabelsMasker, SurfaceLabelsMasker
  43. import json
  44. from itertools import groupby
  45. from os import path, rename
  46. from ..retrieval.cache import CACHE
  47. @dataclass
  48. class MapAssignment:
  49. input_structure: int
  50. centroid: Union[Tuple[np.ndarray], point.Point]
  51. volume: int
  52. fragment: str
  53. map_value: np.ndarray
  54. time: Union[int, float, None]
  55. @dataclass
  56. class AssignImageResult(CompareMapsResult, MapAssignment):
  57. pass
  58. class _CompressedMapSpec(NamedTuple):
  59. """What a compression produces, before `compress()` wraps it as a Map."""
  60. volumes: List[_volume.Volume]
  61. indices: Dict[str, List[Dict]]
  62. class Map(concept.AtlasConcept, configuration_folder="maps"):
  63. def __init__(
  64. self,
  65. identifier: str,
  66. name: str,
  67. space_spec: dict,
  68. parcellation_spec: dict,
  69. indices: Dict[str, List[Dict]],
  70. volumes: list = [],
  71. shortname: str = "",
  72. description: str = "",
  73. modality: str = None,
  74. publications: list = [],
  75. datasets: list = [],
  76. prerelease: bool = False,
  77. ):
  78. """
  79. Constructs a new parcellation object.
  80. Parameters
  81. ----------
  82. identifier: str
  83. Unique identifier of the parcellation
  84. name: str
  85. Human-readable name of the parcellation
  86. space_spec: dict
  87. Specification of the space (use @id or name fields)
  88. parcellation_spec: str
  89. Specification of the parcellation (use @id or name fields)
  90. indices: dict
  91. Dictionary of indices for the brain regions.
  92. Keys are exact region names.
  93. Per region name, a list of dictionaries with fields "volume" and "label" is expected,
  94. where "volume" points to the index of the Volume object where this region is mapped,
  95. and optional "label" is the voxel label for that region.
  96. For continuous / probability maps, the "label" can be null or omitted.
  97. For single-volume labelled maps, the "volume" can be null or omitted.
  98. volumes: list[Volume]
  99. parcellation volumes
  100. shortname: str, optional
  101. Shortform of human-readable name
  102. description: str, optional
  103. Textual description of the parcellation
  104. modality: str, default: None
  105. Specification of the modality used for creating the parcellation
  106. publications: list
  107. List of associated publications, each a dictionary with "doi" and/or "citation" fields
  108. datasets : list
  109. datasets associated with this concept
  110. """
  111. concept.AtlasConcept.__init__(
  112. self,
  113. identifier=identifier,
  114. name=name,
  115. species=None, # inherits species from space
  116. shortname=shortname,
  117. description=description,
  118. publications=publications,
  119. datasets=datasets,
  120. modality=modality,
  121. prerelease=prerelease,
  122. )
  123. self._space_spec = space_spec
  124. self._parcellation_spec = parcellation_spec
  125. # Since the volumes might include 4D arrays, where the actual
  126. # volume index points to a z coordinate, we create subvolume
  127. # indexers from the given volume provider if 'z' is specified.
  128. self._indices: Dict[str, List[MapIndex]] = {}
  129. self.volumes: List[_volume.Volume] = []
  130. remap_volumes = {}
  131. # TODO: This assumes knowledge of the preconfigruation specs wrt. z.
  132. # z to subvolume conversion should probably go to the factory.
  133. for regionname, indexlist in indices.items():
  134. k = clear_name(regionname)
  135. self._indices[k] = []
  136. for index in indexlist:
  137. vol = index.get('volume', 0)
  138. assert vol in range(len(volumes))
  139. z = index.get('z')
  140. if (vol, z) not in remap_volumes:
  141. if z is None:
  142. self.volumes.append(volumes[vol])
  143. else:
  144. self.volumes.append(_volume.Subvolume(volumes[vol], z))
  145. remap_volumes[vol, z] = len(self.volumes) - 1
  146. self._indices[k].append(
  147. MapIndex(volume=remap_volumes[vol, z], label=index.get('label'), fragment=index.get('fragment'))
  148. )
  149. # make sure the indices are unique - each map/label pair should appear at most once
  150. all_indices = sum(self._indices.values(), [])
  151. seen = set()
  152. duplicates = {x for x in all_indices if x in seen or seen.add(x)}
  153. self._nonunique_indices = duplicates
  154. self._affine_cached = None
  155. self._compressed_cached: Dict[tuple, "Map"] = {}
  156. @property
  157. def key(self):
  158. _id = self.id
  159. return create_key(_id[len("siibra-map-v0.0.1"):])
  160. @property
  161. def species(self) -> Species:
  162. # lazy implementation
  163. if self._species_cached is None:
  164. self._species_cached = self.space.species
  165. return self.space._species_cached
  166. def get_index(self, region: Union[str, "Region"]):
  167. """
  168. Returns the unique index corresponding to the specified region.
  169. Tip
  170. ----
  171. Use find_indices() method for a less strict search returning all matches.
  172. Parameters
  173. ----------
  174. region: str or Region
  175. Returns
  176. -------
  177. MapIndex
  178. Raises
  179. ------
  180. NonUniqueIndexError
  181. If not unique or not defined in this parcellation map.
  182. """
  183. matches = self.find_indices(region)
  184. if len(matches) > 1:
  185. # if there is an exact match, we still use it. If not, we cannot proceed.
  186. regionname = region.name if isinstance(region, _region.Region) \
  187. else region
  188. for index, matched_name in matches.items():
  189. if matched_name == regionname:
  190. return index
  191. raise exceptions.NonUniqueIndexError(
  192. f"The specification '{region}' matches multiple mapped "
  193. f"structures in {str(self)}: {list(matches.values())}"
  194. )
  195. elif len(matches) == 0:
  196. raise exceptions.NonUniqueIndexError(
  197. f"The specification '{region}' does not match to any structure mapped in {self}"
  198. )
  199. else:
  200. return next(iter(matches))
  201. def find_indices(self, region: Union[str, "Region"]):
  202. """
  203. Returns the volume/label indices in this map which match the given
  204. region specification.
  205. Parameters
  206. ----------
  207. region: str or Region
  208. Returns
  209. -------
  210. dict
  211. - keys: MapIndex
  212. - values: region name
  213. """
  214. if region in self._indices:
  215. return {
  216. idx: region
  217. for idx in self._indices[region]
  218. }
  219. regionname = region.name if isinstance(region, _region.Region) else region
  220. matched_region_names = set(_.name for _ in (self.parcellation.find(regionname)))
  221. matches = matched_region_names & self._indices.keys()
  222. if len(matches) == 0:
  223. logger.warning(f"Region {regionname} not defined in {self}")
  224. return {
  225. idx: regionname
  226. for regionname in matches
  227. for idx in self._indices[regionname]
  228. }
  229. def get_region(self, label: int = None, volume: int = 0, index: MapIndex = None):
  230. """
  231. Returns the region mapped by the given index, if any.
  232. Tip
  233. ----
  234. Use get_index() or find_indices() methods to obtain the MapIndex.
  235. Parameters
  236. ----------
  237. label: int, default: None
  238. volume: int, default: 0
  239. index: MapIndex, default: None
  240. Returns
  241. -------
  242. Region
  243. A region object defined in the parcellation map.
  244. """
  245. if isinstance(label, MapIndex) and index is None:
  246. raise TypeError("Specify MapIndex with 'index' keyword.")
  247. if index is None:
  248. index = MapIndex(volume, label)
  249. matches = [
  250. regionname
  251. for regionname, indexlist in self._indices.items()
  252. if index in indexlist
  253. ]
  254. if len(matches) == 0:
  255. logger.warning(f"Index {index} not defined in {self}")
  256. return None
  257. elif len(matches) == 1:
  258. return self.parcellation.get_region(matches[0])
  259. else:
  260. # this should not happen, already tested in constructor
  261. raise RuntimeError(f"Index {index} is not unique in {self}")
  262. @property
  263. def space(self):
  264. for key in ["@id", "name"]:
  265. if key in self._space_spec:
  266. return space.Space.get_instance(self._space_spec[key])
  267. return space.Space(None, "Unspecified space", species=Species.UNSPECIFIED_SPECIES)
  268. @property
  269. def parcellation(self):
  270. for key in ["@id", "name"]:
  271. if key in self._parcellation_spec:
  272. return parcellation.Parcellation.get_instance(self._parcellation_spec[key])
  273. logger.warning(
  274. f"Cannot determine parcellation of {self.__class__.__name__} "
  275. f"{self.name} from {self._parcellation_spec}"
  276. )
  277. return None
  278. @property
  279. def labels(self):
  280. """
  281. The set of all label indices defined in this map, including "None" if
  282. not defined for one or more regions.
  283. """
  284. return {d.label for v in self._indices.values() for d in v}
  285. @property
  286. def has_unique_labels(self) -> bool:
  287. """
  288. True if every mapped region carries a distinct label across all volumes and
  289. fragments. Surface and fragmented volumetric maps are commonly labelled per
  290. hemisphere, so the same label may denote two different regions.
  291. """
  292. labels = [ix.label for ixs in self._indices.values() for ix in ixs]
  293. return None not in labels and len(labels) == len(set(labels))
  294. @property
  295. def maptype(self) -> MapType:
  296. if all(isinstance(_, int) for _ in self.labels):
  297. return MapType.LABELLED
  298. elif self.labels == {None}:
  299. return MapType.STATISTICAL
  300. else:
  301. raise RuntimeError(
  302. f"Inconsistent label indices encountered in {self}"
  303. )
  304. def __len__(self):
  305. return len(self.volumes)
  306. @property
  307. def regions(self):
  308. return list(self._indices)
  309. def get_volume(
  310. self,
  311. region: Union[str, "Region"] = None,
  312. *,
  313. index: MapIndex = None,
  314. **kwargs,
  315. ) -> Union[_volume.Volume, _volume.FilteredVolume, _volume.Subvolume]:
  316. try:
  317. length = len([arg for arg in [region, index] if arg is not None])
  318. assert length == 1
  319. except AssertionError:
  320. if length > 1:
  321. raise exceptions.ExcessiveArgumentException(
  322. "One and only one of region or index can be defined for `get_volume`."
  323. )
  324. mapindex = None
  325. if region is not None:
  326. try:
  327. assert isinstance(region, (str, _region.Region))
  328. except AssertionError:
  329. raise TypeError(f"Please provide a region name or region instance, not a {type(region)}")
  330. mapindex = self.get_index(region)
  331. if index is not None:
  332. assert isinstance(index, MapIndex)
  333. mapindex = index
  334. if mapindex is None:
  335. if len(self) == 1:
  336. mapindex = MapIndex(volume=0, label=None)
  337. elif len(self) > 1:
  338. assert self.maptype == MapType.LABELLED, f"Cannot merge multiple volumes of map type {self.maptype}. Please specify a region or index."
  339. logger.info(
  340. "Map provides multiple volumes and no region specification is"
  341. " provided. Reducing them to a single volume resampled to space template."
  342. )
  343. labels = list(range(1, len(self.volumes) + 1)) if self.labels == {1} else None
  344. merged_volume = _volume.ReducedVolume(self.volumes, labels)
  345. return merged_volume
  346. else:
  347. raise exceptions.NoVolumeFound("Map provides no volumes.")
  348. kwargs_fragment = kwargs.pop("fragment", None)
  349. if kwargs_fragment is not None:
  350. if (mapindex.fragment is not None) and (kwargs_fragment != mapindex.fragment):
  351. raise exceptions.ConflictingArgumentException(
  352. f"Conflicting specifications for fetching volume fragment{f' for region {region}'}: "
  353. f"supplied: {kwargs_fragment}, preconfigured: {mapindex.fragment}"
  354. )
  355. if mapindex.volume is None:
  356. mapindex.volume = 0
  357. if mapindex.volume >= len(self.volumes):
  358. raise IndexError(
  359. f"{self} provides {len(self)} mapped volumes, but #{mapindex.volume} was requested."
  360. )
  361. if mapindex.label is None and mapindex.fragment is None:
  362. return self.volumes[mapindex.volume]
  363. return _volume.FilteredVolume(
  364. parent_volume=self.volumes[mapindex.volume],
  365. label=mapindex.label,
  366. fragment=kwargs_fragment or mapindex.fragment,
  367. )
  368. def fetch(
  369. self,
  370. region: Union[str, "Region"] = None,
  371. *,
  372. index: MapIndex = None,
  373. **fetch_kwargs
  374. ):
  375. """
  376. Fetches one particular volume of this parcellation map.
  377. If there's only one volume, this is the default, otherwise further
  378. specification is requested:
  379. - the volume index,
  380. - the MapIndex (which results in a regional map being returned)
  381. You might also consider fetch_iter() to iterate the volumes, or
  382. compress() to produce a single-volume parcellation map.
  383. Parameters
  384. ----------
  385. region: str, Region
  386. Specification of a region name, resulting in a regional map
  387. (mask or statistical map) to be returned.
  388. index: MapIndex
  389. Explicit specification of the map index, typically resulting
  390. in a regional map (mask or statistical map) to be returned.
  391. Note that supplying 'region' will result in retrieving the map index of that region
  392. automatically.
  393. **fetch_kwargs
  394. - resolution_mm: resolution in millimeters as float or a tuple of floats
  395. - format: the format of the volume, like "mesh" or "nii"
  396. - voi: a BoundingBox of interest
  397. Note
  398. ----
  399. Not all keyword arguments are supported for volume formats. Format
  400. is restricted by available formats (check formats property).
  401. Returns
  402. -------
  403. An image or mesh
  404. """
  405. vol = self.get_volume(region=region, index=index, **fetch_kwargs)
  406. return vol.fetch(**fetch_kwargs)
  407. def fetch_iter(self, **kwargs):
  408. """
  409. Returns an iterator to fetch all mapped volumes sequentially.
  410. All arguments are passed on to function Map.fetch(). By default, it
  411. will go through all fragments as well.
  412. """
  413. fragments = {kwargs.pop('fragment', None)} or self.fragments or {None}
  414. return (
  415. self.fetch(
  416. index=MapIndex(volume=i, label=None, fragment=frag), **kwargs
  417. )
  418. for frag in fragments
  419. for i in range(len(self))
  420. )
  421. @property
  422. def provides_image(self):
  423. return any(v.provides_image for v in self.volumes)
  424. @property
  425. def fragments(self):
  426. return sorted({
  427. index.fragment
  428. for indices in self._indices.values()
  429. for index in indices
  430. if index.fragment is not None
  431. })
  432. @property
  433. def provides_mesh(self):
  434. return any(v.provides_mesh for v in self.volumes)
  435. @property
  436. def formats(self):
  437. return {f for v in self.volumes for f in v.formats}
  438. @property
  439. def is_labelled(self):
  440. return self.maptype == MapType.LABELLED
  441. @property
  442. def affine(self):
  443. if self._affine_cached is None:
  444. # we compute the affine from a volumetric volume provider
  445. for fmt in _volume.Volume.SUPPORTED_FORMATS:
  446. if fmt not in _volume.Volume.MESH_FORMATS:
  447. if fmt not in self.formats:
  448. continue
  449. try:
  450. self._affine_cached = self.fetch(index=MapIndex(volume=0), format=fmt).affine
  451. break
  452. except Exception:
  453. logger.debug("Caught exceptions:\n", exc_info=1)
  454. continue
  455. else:
  456. raise RuntimeError(f"No volumetric provider in {self} to derive the affine matrix.")
  457. if not isinstance(self._affine_cached, np.ndarray):
  458. logger.error("invalid affine:", self._affine_cached)
  459. return self._affine_cached
  460. def __iter__(self):
  461. return self.fetch_iter()
  462. def compress(self, **kwargs) -> "Map":
  463. """
  464. Convert this map into an equivalent map whose labels are unique across the
  465. whole image or surface, re-labelling regions sequentially from 1.
  466. Volumetric maps are merged into a single labelled volume on the grid of the
  467. space template. Surface maps keep their fragments, since a vertex belongs to
  468. exactly one of them, and are relabelled so that a label identifies one region
  469. rather than one region per hemisphere.
  470. The compressed image is built and kept as a memory-mapped array on disk, so
  471. neither building nor using it holds the full volume in RAM. Results are also
  472. cached per instance, keyed by the fetch arguments.
  473. Note
  474. ----
  475. Labels that no region claims are dropped to background. Surface maps
  476. sometimes carry such labels for technical reasons (e.g. a brainstem label
  477. so the map loads in freesurfer) and they have no name to report.
  478. Parameters
  479. ----------
  480. **kwargs
  481. Fetch arguments applied both to the space template, which defines the
  482. output grid, and to the mapped volumes, so that sources are read at the
  483. same resolution and resampling is usually unnecessary. `variant` is
  484. passed to the template only.
  485. Returns
  486. -------
  487. parcellationmap.Map
  488. Raises
  489. ------
  490. ValueError
  491. If this map is not labelled.
  492. RuntimeError
  493. If there is nothing to merge.
  494. """
  495. if not self.is_labelled:
  496. raise ValueError(f"Compression is not possible for {self.maptype} maps.")
  497. if len(self.volumes) == 1 and (not self.fragments or self.has_unique_labels):
  498. raise RuntimeError(
  499. "The map cannot be compressed: it is already a single, uniquely labelled volume."
  500. )
  501. key = tuple(sorted(kwargs.items()))
  502. try:
  503. hash(key)
  504. except TypeError: # an unhashable fetch argument, e.g. target_affine
  505. logger.debug(f"Cannot cache compression of {self} for {kwargs}.")
  506. key = None
  507. if key is not None and key in self._compressed_cached:
  508. return self._compressed_cached[key]
  509. entries = sorted(
  510. (index.volume, index.fragment or "", index.label, regionname)
  511. for regionname, indices in self._indices.items()
  512. for index in indices
  513. )
  514. relabelling_plan = [
  515. (MapIndex(volume=vol, fragment=frag or None, label=label), regionname, newlabel)
  516. for newlabel, (vol, frag, label, regionname) in enumerate(entries, start=1)
  517. ]
  518. # everything that shapes the result, so a changed configuration or changed
  519. # fetch arguments never reuse a stale artifact
  520. signature = json.dumps(
  521. {
  522. "map": self.id,
  523. "space": self.space.id,
  524. "kwargs": {k: str(v) for k, v in sorted(kwargs.items())},
  525. "indices": [
  526. [index.volume, index.fragment, index.label, regionname]
  527. for index, regionname, _ in relabelling_plan
  528. ],
  529. },
  530. sort_keys=True, # ensure_ascii=True by default: CACHE encodes as ascii
  531. )
  532. if self.provides_image:
  533. spec = self._compress_image_map(relabelling_plan, signature, **kwargs)
  534. elif self.provides_mesh:
  535. if kwargs:
  536. logger.info(f"Fetch arguments {list(kwargs)} are ignored when compressing a surface map.")
  537. spec = self._compress_surface_map(relabelling_plan, signature)
  538. else:
  539. raise NotImplementedError(f"{self} provides neither image nor mesh data to compress.")
  540. compressed = Map(
  541. identifier=f"{create_key(self.name)}_compressed",
  542. name=f"{self.name} compressed",
  543. space_spec=self._space_spec,
  544. parcellation_spec=self._parcellation_spec,
  545. indices=spec.indices,
  546. volumes=spec.volumes,
  547. )
  548. if key is not None:
  549. self._compressed_cached[key] = compressed
  550. return compressed
  551. def _compress_image_map(
  552. self, plan: List[Tuple[MapIndex, str, int]], signature: str, **kwargs
  553. ) -> "_CompressedMapSpec":
  554. """
  555. Merge the volumes and fragments of a volumetric map into a single labelled
  556. volume on the template grid. See `compress()`.
  557. Source volumes are streamed one at a time into a memory-mapped output array,
  558. so peak memory is one source volume rather than the whole map. Where regions
  559. overlap, the later entry of the relabelling plan wins; the number of
  560. contested voxels is reported.
  561. """
  562. cachefile = CACHE.build_filename(signature, suffix=".npy")
  563. metafile = f"{cachefile}.json"
  564. if path.isfile(cachefile) and path.isfile(metafile):
  565. try:
  566. with open(metafile) as f:
  567. meta = json.load(f)
  568. data = np.lib.format.open_memmap(cachefile, mode="r")
  569. logger.debug(f"Reusing the compressed {self} from {cachefile}")
  570. return _CompressedMapSpec(
  571. volumes=[_volume.from_array(
  572. data, np.array(meta["affine"]), self.space.id,
  573. name=self.name + " compressed", cache=False,
  574. )],
  575. indices=meta["indices"],
  576. )
  577. except (ValueError, OSError, KeyError):
  578. logger.debug(f"Discarding unreadable compression cache {cachefile}")
  579. # only the grid of the template is needed - never load its data, which for
  580. # large templates (e.g. BigBrain) dominates both time and memory
  581. variant = kwargs.pop("variant", None)
  582. template_img = self.space.get_template(variant=variant).fetch(**kwargs)
  583. shape, affine = tuple(template_img.shape[:3]), template_img.affine
  584. dtype = np.min_scalar_type(len(plan)) # sized by region count, not by the template
  585. nbytes = int(np.prod(shape)) * np.dtype(dtype).itemsize
  586. gib = nbytes / 1024**3
  587. logger.info(
  588. f"Compressing {self} into a {gib:.2f} GiB labelled volume "
  589. f"(shape {shape}, {np.dtype(dtype).name}) from {len(self.volumes)} volume(s)."
  590. )
  591. units = [(unit, list(entries)) for unit, entries in groupby(
  592. plan, key=lambda entry: (entry[0].volume, entry[0].fragment)
  593. )]
  594. tempfile = f"{cachefile}_temp"
  595. result = np.lib.format.open_memmap(tempfile, mode="w+", dtype=dtype, shape=shape)
  596. region_indices = defaultdict(list)
  597. contested, unmapped = 0, []
  598. try:
  599. with provider.SubvolumeProvider.UseCaching():
  600. for (volume, fragment), entries in siibra_tqdm(
  601. units, total=len(units), unit=" maps",
  602. desc=f"Compressing {len(self.volumes)} volume(s) and "
  603. f"{len(self.fragments) or 1} fragment(s) of {self.name}",
  604. disable=len(units) == 1,
  605. ):
  606. img = self.fetch(index=MapIndex(volume=volume, fragment=fragment), **kwargs)
  607. if tuple(img.shape[:3]) == shape and np.allclose(img.affine, affine):
  608. img_data = np.asanyarray(img.dataobj)
  609. else:
  610. logger.debug(f"Compression requires resampling volume {volume} (nearest)")
  611. img_data = np.asanyarray(resample_img_to_img(img, template_img).dataobj)
  612. observed = set(np.unique(img_data)) - {0}
  613. for index, regionname, newlabel in entries:
  614. update_voxels = img_data == index.label
  615. if not update_voxels.any():
  616. unmapped.append(regionname)
  617. contested += int(np.count_nonzero(result[update_voxels]))
  618. result[update_voxels] = newlabel
  619. region_indices[regionname].append({"volume": 0, "label": newlabel})
  620. observed.discard(index.label)
  621. if observed:
  622. logger.warning(
  623. f"Labels {sorted(observed)} are observed in volume {volume} "
  624. f"(fragment {fragment}) of {self}, but no region is defined for them."
  625. )
  626. del img_data
  627. result.flush()
  628. finally:
  629. del result # close the memmap before renaming, required on Windows
  630. rename(tempfile, cachefile)
  631. with open(metafile, "w") as f:
  632. json.dump({"affine": affine.tolist(), "indices": region_indices}, f)
  633. if contested:
  634. logger.info(
  635. f"{contested} voxel(s) are mapped by more than one region in {self}; "
  636. "the last entry of the relabelling order was kept for each."
  637. )
  638. if unmapped:
  639. logger.warning(f"{len(unmapped)} region(s) have no voxels after compression:\n{unmapped}")
  640. return _CompressedMapSpec(
  641. volumes=[_volume.from_array(
  642. np.lib.format.open_memmap(cachefile, mode="r"), affine, self.space.id,
  643. name=self.name + " compressed", cache=False,
  644. )],
  645. indices=region_indices,
  646. )
  647. def _compress_surface_map(
  648. self, plan: List[Tuple[MapIndex, str, int]], signature: str
  649. ) -> "_CompressedMapSpec":
  650. """
  651. Relabel the fragments of a surface map so that every region has a globally
  652. unique label. Fragments are preserved, since a vertex belongs to exactly one
  653. of them. See `compress()`.
  654. The label arrays are one per-vertex value per fragment, well under a megabyte
  655. even for the densest fsaverage mesh, so they are held in memory rather than
  656. memory-mapped, and written as GIFTI label files.
  657. Results are cached per fragment as GIFTI label files; a cache hit skips the
  658. relabelling and therefore also its warnings about unnamed labels and empty
  659. regions.
  660. """
  661. from nibabel import gifti
  662. from os import replace
  663. if len(self.volumes) > 1:
  664. raise NotImplementedError(
  665. f"{self} provides {len(self.volumes)} surface volumes; compression expects one."
  666. )
  667. prov = self.volumes[0]._providers["gii-label"]
  668. filemap, region_indices = {}, defaultdict(list)
  669. units = [(unit, list(entries)) for unit, entries in groupby(
  670. plan, key=lambda entry: (entry[0].volume, entry[0].fragment)
  671. )]
  672. for (_, fragment), entries in siibra_tqdm(
  673. units, total=len(units), unit=" fragments",
  674. desc=f"Relabelling {len(units)} surface fragment(s) of {self.name}",
  675. ):
  676. filename = CACHE.build_filename(f"{signature}-{fragment}", suffix=".label.gii")
  677. if not path.isfile(filename):
  678. # read from the provider: Volume.fetch() would also pull the template mesh
  679. data = prov.fetch(fragment=fragment)["labels"]
  680. # GIFTI data arrays support uint8, int32 and float32 only
  681. relabelled = np.zeros_like(data, dtype="int32")
  682. observed = set(np.unique(data)) - {0}
  683. unmapped = []
  684. for index, regionname, newlabel in entries:
  685. selection = data == index.label
  686. if not selection.any():
  687. unmapped.append(regionname)
  688. relabelled[selection] = newlabel
  689. observed.discard(index.label)
  690. if observed:
  691. logger.warning(
  692. f"Labels {sorted(observed)} are observed in fragment '{fragment}' of "
  693. f"{self}, but no region is defined for them."
  694. )
  695. if unmapped:
  696. logger.warning(
  697. f"{len(unmapped)} region(s) have no vertices in fragment "
  698. f"'{fragment}' of {self}:\n{unmapped}"
  699. )
  700. tempfile = CACHE.build_filename(f"{signature}-{fragment}-temp", suffix=".label.gii")
  701. gifti.GiftiImage(darrays=[
  702. gifti.GiftiDataArray(relabelled, intent="NIFTI_INTENT_LABEL")
  703. ]).to_filename(tempfile)
  704. replace(tempfile, filename)
  705. for index, regionname, newlabel in entries:
  706. region_indices[regionname].append(
  707. {"volume": 0, "fragment": fragment, "label": newlabel}
  708. )
  709. filemap[fragment] = filename
  710. return _CompressedMapSpec(
  711. volumes=[_volume.from_file(
  712. filemap, space=self.space.id, name=self.name + " compressed", format="gii-label",
  713. )],
  714. indices=region_indices,
  715. )
  716. def compute_centroids(self, split_components: bool = True, **fetch_kwargs) -> Dict[str, pointcloud.PointCloud]:
  717. """
  718. Compute a dictionary of all regions in this map to their centroids.
  719. By default, the regional masks will be split to connected components
  720. and each point in the PointCloud corresponds to a region component.
  721. Parameters
  722. ----------
  723. split_components: bool, default: True
  724. If True, finds the spatial properties for each connected component
  725. found by skimage.measure.label.
  726. Returns
  727. -------
  728. Dict[str, point.Point]
  729. Region names as keys and computed centroids as items.
  730. """
  731. assert self.provides_image, "Centroid computation for meshes is not supported yet."
  732. centroids = dict()
  733. for regionname, indexlist in siibra_tqdm(
  734. self._indices.items(), unit="regions", desc="Computing centroids"
  735. ):
  736. assert regionname not in centroids
  737. # get the mask of the region in this map
  738. with QUIET:
  739. if len(indexlist) >= 1:
  740. merged_volume = _volume.merge(
  741. [
  742. _volume.from_nifti(
  743. self.fetch(index=index, **fetch_kwargs),
  744. self.space,
  745. f"{self.name} - {index}"
  746. )
  747. for index in indexlist
  748. ],
  749. labels=[1] * len(indexlist)
  750. )
  751. mapimg = merged_volume.fetch()
  752. elif len(indexlist) == 1:
  753. index = indexlist[0]
  754. mapimg = self.fetch(index=index, **fetch_kwargs) # returns a mask of the region
  755. props = _volume.ComponentSpatialProperties.compute_from_image(
  756. img=mapimg,
  757. space=self.space,
  758. split_components=split_components,
  759. )
  760. try:
  761. centroids[regionname] = pointcloud.from_points([c.centroid for c in props])
  762. except exceptions.EmptyPointCloudError:
  763. centroids[regionname] = None
  764. return centroids
  765. def get_resampled_template(self, **fetch_kwargs) -> _volume.Volume:
  766. """
  767. Resample the reference space template to fetched map image. Uses
  768. nilearn.image.resample_to_img to resample the template.
  769. Parameters
  770. ----------
  771. **fetch_kwargs: takes the arguments of Map.fetch()
  772. Returns
  773. -------
  774. Volume
  775. """
  776. from nilearn.image import resample_to_img
  777. source_template = self.space.get_template().fetch()
  778. map_image = self.fetch(**fetch_kwargs)
  779. img = resample_to_img(source_template, map_image, interpolation='continuous')
  780. return _volume.from_array(
  781. data=img.dataobj,
  782. affine=img.affine,
  783. space=self.space,
  784. name=f"{source_template} resampled to coordinate system of {self}"
  785. )
  786. def colorize(
  787. self,
  788. values: Union[dict, "pd.Series", "pd.DataFrame"],
  789. background_label: Union[int, float] = 0,
  790. **masker_kwargs
  791. ):
  792. """Colorize the map with the provided regional values.
  793. Parameters
  794. ----------
  795. values : dict
  796. Dictionary mapping regions to values
  797. Return
  798. ------
  799. Nifti1Image
  800. """
  801. if not self.is_labelled:
  802. raise NotImplementedError("Since statistical maps can overlap, this is not yet implemented.")
  803. if isinstance(values, dict):
  804. resolved = {}
  805. for spec, value in values.items():
  806. matched = set(self.find_indices(spec).values()) # {MapIndex: regionname}
  807. if not matched:
  808. logger.warning(f"'{spec}' is not mapped in {self} - skipped in colorization.")
  809. continue
  810. for regionname in matched:
  811. resolved[regionname] = value
  812. if not resolved:
  813. raise ValueError(f"None of the {len(values)} provided keys are mapped in {self}.")
  814. values = pd.Series({r: resolved.get(r, background_label) for r in self.regions})
  815. masker_kwargs.setdefault("background_label", background_label)
  816. masker = self.as_nilearn_masker(**masker_kwargs)
  817. masker.fit()
  818. # ensure the order of columns follow the bids table used for masker
  819. values = values[masker.lut["name"]]
  820. return masker.inverse_transform(values)
  821. def get_colormap(self, region_specs: Iterable = None, *, fill_uncolored: bool = False):
  822. """
  823. Generate a matplotlib colormap from known rgb values of label indices.
  824. Parameters
  825. ----------
  826. region_specs: iterable(regions), optional
  827. Optional parameter to only color the desired regions.
  828. fill_uncolored: bool , optional
  829. If a region has no preconfigured color, a color will be randomly (reproducible) created.
  830. Returns
  831. -------
  832. ListedColormap
  833. """
  834. try:
  835. from matplotlib.colors import ListedColormap
  836. except ImportError as e:
  837. logger.error(
  838. "matplotlib not available. Please install matplotlib to create a matplotlib colormap."
  839. )
  840. raise e
  841. if fill_uncolored:
  842. seed = len(self.regions)
  843. np.random.seed(seed)
  844. logger.info(f"Random colors are allowed for regions without preconfgirued colors. Random seed: {seed}.")
  845. colors = {}
  846. if region_specs is not None:
  847. include_region_names = {
  848. self.parcellation.get_region(region_spec).name for region_spec in region_specs
  849. }
  850. else:
  851. include_region_names = None
  852. use_volindx = False
  853. if len(self.volumes) > 1 and self.labels == {1}:
  854. use_volindx = True
  855. logger.info("Using sequential relabling determined by the order of volumes.")
  856. no_predefined_color = []
  857. for regionname, indices in self._indices.items():
  858. for index in indices:
  859. if index.label is None:
  860. continue
  861. if (include_region_names is not None) and (regionname not in include_region_names):
  862. continue
  863. else:
  864. region = self.get_region(index=index)
  865. if region.rgb is not None:
  866. if use_volindx:
  867. colors[index.volume + 1] = region.rgb
  868. else:
  869. colors[index.label] = region.rgb
  870. elif fill_uncolored:
  871. random_clr = [np.random.randint(0, 255) for r in range(3)]
  872. while random_clr in list(colors.values()):
  873. random_clr = [np.random.randint(0, 255) for r in range(3)]
  874. colors[index.label] = random_clr
  875. else:
  876. no_predefined_color.append(region.name)
  877. if len(colors) == 0:
  878. raise exceptions.NoPredifinedColormapException(
  879. f"There is no predefined/preconfigured colormap for '{self}'."
  880. "Set `fill_uncolored=True` to get a reproducible colormap."
  881. )
  882. if no_predefined_color:
  883. logger.info(
  884. f"No preconfigured color found for the following regions."
  885. "Use `fill_uncolored=True` to display with a non-background color.\n"
  886. f"{no_predefined_color}"
  887. )
  888. max_label_index = max(colors.keys())
  889. palette = np.array(
  890. [
  891. list(colors[i]) + [1] if i in colors else [0, 0, 0, 0]
  892. for i in range(max_label_index + 1)
  893. ]
  894. ) / [255, 255, 255, 1]
  895. return ListedColormap(palette)
  896. def sample_locations(self, regionspec, numpoints: int):
  897. """ Sample 3D locations inside a given region.
  898. The probability distribution is approximated from the region mask based
  899. on the squared distance transform.
  900. Parameters
  901. ----------
  902. regionspec: Region or str
  903. Region to be used
  904. numpoints: int
  905. Number of samples to draw
  906. Returns
  907. -------
  908. PointCloud
  909. Sample points in physical coordinates corresponding to this
  910. parcellationmap
  911. """
  912. index = self.get_index(regionspec)
  913. mask = self.fetch(index=index)
  914. arr = np.asanyarray(mask.dataobj)
  915. if arr.dtype.char in np.typecodes['AllInteger']:
  916. # a binary mask - use distance transform to get sampling weights
  917. W = distance_transform_edt(np.asanyarray(mask.dataobj))**2
  918. else:
  919. # a statistical map - interpret directly as weights
  920. W = arr
  921. p = (W / W.sum()).ravel()
  922. XYZ_ = np.array(
  923. np.unravel_index(np.random.choice(len(p), numpoints, p=p), W.shape)
  924. ).T
  925. XYZ = np.dot(mask.affine, np.c_[XYZ_, np.ones(numpoints)].T)[:3, :].T
  926. return pointcloud.PointCloud(XYZ, space=self.space)
  927. def to_sparse(self):
  928. """
  929. Creates a SparseMap object from this parcellation map object.
  930. Returns
  931. -------
  932. SparseMap
  933. """
  934. from .sparsemap import SparseMap
  935. indices = {
  936. regionname: [
  937. {'volume': idx.volume, 'label': idx.label, 'fragment': idx.fragment}
  938. for idx in indexlist
  939. ]
  940. for regionname, indexlist in self._indices.items()
  941. }
  942. return SparseMap(
  943. identifier=self.id,
  944. name=self.name,
  945. space_spec={'@id': self.space.id},
  946. parcellation_spec={'@id': self.parcellation.id},
  947. indices=indices,
  948. volumes=self.volumes,
  949. shortname=self.shortname,
  950. description=self.description,
  951. modality=self.modality,
  952. publications=self.publications,
  953. datasets=self.datasets
  954. )
  955. def _read_voxel(
  956. self,
  957. x: Union[int, np.ndarray, List],
  958. y: Union[int, np.ndarray, List],
  959. z: Union[int, np.ndarray, List]
  960. ):
  961. def _read_voxels_from_volume(xyz, volimg):
  962. xyz = np.stack(xyz, axis=1)
  963. valid_points_mask = np.all([(0 <= di) & (di < vol_size) for vol_size, di in zip(volimg.shape, xyz.T)], axis=0)
  964. x, y, z = xyz[valid_points_mask].T
  965. valid_points_indices, *_ = np.where(valid_points_mask)
  966. valid_data_points = np.asanyarray(volimg.dataobj)[x, y, z]
  967. return zip(valid_points_indices, valid_data_points)
  968. # integers are just single-element arrays, cast to avoid an extra code branch for integers
  969. x, y, z = [np.array(di) for di in (x, y, z)]
  970. fragments = self.fragments or {None}
  971. return [
  972. (pointindex, volume, fragment, data_point)
  973. for fragment in fragments
  974. for volume, volimg in enumerate(self.fetch_iter(fragment=fragment))
  975. # transformations or user input might produce points outside the volume, filter these out.
  976. for (pointindex, data_point) in _read_voxels_from_volume((x, y, z), volimg)
  977. ]
  978. def _assign(
  979. self,
  980. item: location.Location,
  981. minsize_voxel=1,
  982. lower_threshold=0.0,
  983. **kwargs
  984. ) -> List[Union[MapAssignment, AssignImageResult]]:
  985. """
  986. For internal use only. Returns a dataclass, which provides better static type checking.
  987. """
  988. if isinstance(item, point.Point):
  989. return self._assign_points(
  990. pointcloud.PointCloud([item], item.space, sigma_mm=item.sigma),
  991. lower_threshold
  992. )
  993. if isinstance(item, pointcloud.PointCloud):
  994. return self._assign_points(item, lower_threshold)
  995. if isinstance(item, _volume.Volume):
  996. if isinstance(item, _volume.TimeSeriesVolume):
  997. return self._assign_timeseries_volume(
  998. queryvolume=item,
  999. lower_threshold=lower_threshold,
  1000. minsize_voxel=minsize_voxel,
  1001. **kwargs
  1002. )
  1003. else:
  1004. return self._assign_volume(
  1005. queryvolume=item,
  1006. lower_threshold=lower_threshold,
  1007. minsize_voxel=minsize_voxel,
  1008. **kwargs
  1009. )
  1010. raise RuntimeError(
  1011. f"Items of type {item.__class__.__name__} cannot be used for region assignment."
  1012. )
  1013. def assign(
  1014. self,
  1015. item: location.Location,
  1016. minsize_voxel=1,
  1017. lower_threshold=0.0,
  1018. **kwargs
  1019. ) -> "pd.DataFrame":
  1020. """Assign an input Location to brain regions.
  1021. The input is assumed to be defined in the same coordinate space
  1022. as this parcellation map.
  1023. Parameters
  1024. ----------
  1025. item: Location
  1026. A spatial object defined in the same physical reference space as
  1027. this parcellation map, which could be a point, set of points, or
  1028. image volume. If it is an image, it will be resampled to the same voxel
  1029. space if its affine transformation differs from that of the
  1030. parcellation map. Resampling will use linear interpolation for float
  1031. image types, otherwise nearest neighbor.
  1032. minsize_voxel: int, default: 1
  1033. Minimum voxel size of image components to be taken into account.
  1034. lower_threshold: float, default: 0
  1035. Lower threshold on values in the statistical map. Values smaller
  1036. than this threshold will be excluded from the assignment computation.
  1037. Returns
  1038. -------
  1039. pandas.DataFrame
  1040. A table of associated regions and their scores per component found
  1041. in the input image, or per coordinate provided. The scores are:
  1042. - Value: Maximum value of the voxels in the map covered by an
  1043. input coordinate or input image signal component.
  1044. - Pearson correlation coefficient between the brain region map
  1045. and an input image signal component (NaN for exact coordinates)
  1046. - Contains: Percentage of the brain region map contained in an
  1047. input image signal component, measured from their binarized
  1048. masks as the ratio between the volume of their intersection
  1049. and the volume of the brain region (NaN for exact coordinates)
  1050. - Contained: Percentage of an input image signal component
  1051. contained in the brain region map, measured from their binary
  1052. masks as the ratio between the volume of their intersection and
  1053. the volume of the input image signal component (NaN for exact
  1054. coordinates)
  1055. """
  1056. assignments = self._assign(item, minsize_voxel, lower_threshold, **kwargs)
  1057. # format assignments as pandas dataframe
  1058. columns = [
  1059. "input structure",
  1060. "time",
  1061. "centroid",
  1062. "volume",
  1063. "fragment",
  1064. "region",
  1065. "correlation",
  1066. "intersection over union",
  1067. "map value",
  1068. "map weighted mean",
  1069. "map containedness",
  1070. "input weighted mean",
  1071. "input containedness"
  1072. ]
  1073. if len(assignments) == 0:
  1074. return pd.DataFrame(columns=columns).dropna(axis='columns', how='all')
  1075. # determine the unique set of observed indices in order to do region lookups
  1076. # only once for each map index occurring in the point list
  1077. labelled = self.is_labelled # avoid calling this in a loop
  1078. observed_indices = { # unique set of observed map indices. NOTE: len(observed_indices) << len(assignments)
  1079. (
  1080. a.volume,
  1081. a.fragment,
  1082. a.map_value if labelled else None
  1083. )
  1084. for a in assignments
  1085. }
  1086. region_lut = { # lookup table of observed region objects
  1087. (v, f, l): self.get_region(
  1088. index=MapIndex(
  1089. volume=int(v),
  1090. label=l if l is None else int(l),
  1091. fragment=f
  1092. )
  1093. )
  1094. for v, f, l in observed_indices
  1095. }
  1096. dataframe_list = []
  1097. for a in assignments:
  1098. item_to_append = {
  1099. "input structure": a.input_structure,
  1100. "time": a.time,
  1101. "centroid": a.centroid,
  1102. "volume": a.volume,
  1103. "fragment": a.fragment,
  1104. "region": region_lut[
  1105. a.volume,
  1106. a.fragment,
  1107. a.map_value if labelled else None
  1108. ],
  1109. }
  1110. # because AssignImageResult is a subclass of Assignment
  1111. # need to check for isinstance AssignImageResult first
  1112. if isinstance(a, AssignImageResult):
  1113. item_to_append = {
  1114. **item_to_append,
  1115. **{
  1116. "correlation": a.correlation,
  1117. "intersection over union": a.intersection_over_union,
  1118. "map value": a.map_value,
  1119. "map weighted mean": a.weighted_mean_of_first,
  1120. "map containedness": a.intersection_over_first,
  1121. "input weighted mean": a.weighted_mean_of_second,
  1122. "input containedness": a.intersection_over_second,
  1123. }
  1124. }
  1125. elif isinstance(a, MapAssignment):
  1126. item_to_append = {
  1127. **item_to_append,
  1128. **{
  1129. "correlation": None,
  1130. "intersection over union": None,
  1131. "map value": a.map_value,
  1132. "map weighted mean": None,
  1133. "map containedness": None,
  1134. "input weighted mean": None,
  1135. "input containedness": None,
  1136. }
  1137. }
  1138. else:
  1139. raise RuntimeError("assignments must be of type Assignment or AssignImageResult!")
  1140. dataframe_list.append(item_to_append)
  1141. return (
  1142. pd.DataFrame(dataframe_list)
  1143. .convert_dtypes() # convert will guess numeric column types
  1144. .reindex(columns=columns)
  1145. .dropna(axis='columns', how='all')
  1146. )
  1147. def _assign_points(self, points: pointcloud.PointCloud, lower_threshold: float) -> List[MapAssignment]:
  1148. """
  1149. assign a PointCloud to this parcellation map.
  1150. Parameters
  1151. -----------
  1152. lower_threshold: float, default: 0
  1153. Lower threshold on values in the statistical map. Values smaller than
  1154. this threshold will be excluded from the assignment computation.
  1155. """
  1156. assignments = []
  1157. if points.space != self.space:
  1158. logger.info(
  1159. f"Coordinates will be converted from {points.space.name} "
  1160. f"to {self.space.name} space for assignment."
  1161. )
  1162. # convert sigma to voxel coordinates
  1163. scaling = np.array(
  1164. [np.linalg.norm(self.affine[:, i]) for i in range(3)]
  1165. ).mean()
  1166. phys2vox = np.linalg.inv(self.affine)
  1167. # if all points have the same sigma, and lead to a standard deviation
  1168. # below 3 voxels, we are much faster with a multi-coordinate readout.
  1169. if points.has_constant_sigma:
  1170. sigma_vox = points.sigma[0] / scaling
  1171. if sigma_vox < 3:
  1172. pts_warped = points.warp(self.space.id)
  1173. X, Y, Z = (np.dot(phys2vox, pts_warped.homogeneous.T) + 0.5).astype("int")[:3]
  1174. for pointindex, vol, frag, value in self._read_voxel(X, Y, Z):
  1175. if value > lower_threshold:
  1176. position = pts_warped[pointindex].coordinate
  1177. assignments.append(
  1178. MapAssignment(
  1179. input_structure=pointindex,
  1180. centroid=tuple(position),
  1181. volume=vol,
  1182. fragment=frag,
  1183. map_value=value,
  1184. time=None,
  1185. )
  1186. )
  1187. return assignments
  1188. # if we get here, we need to handle each point independently.
  1189. # This is much slower but more precise in dealing with the uncertainties
  1190. # of the coordinates.
  1191. for pointindex, pt in siibra_tqdm(
  1192. enumerate(points.warp(self.space.id)),
  1193. total=len(points), desc="Assigning points",
  1194. ):
  1195. sigma_vox = pt.sigma / scaling
  1196. if sigma_vox < 3:
  1197. # voxel-precise - just read out the value in the maps
  1198. N = len(self)
  1199. logger.debug(f"Assigning coordinate {tuple(pt)} to {N} maps")
  1200. x, y, z = (np.dot(phys2vox, pt.homogeneous) + 0.5).astype("int")[:3]
  1201. values = self._read_voxel(x, y, z)
  1202. for _, vol, frag, value in values:
  1203. if value > lower_threshold:
  1204. assignments.append(
  1205. MapAssignment(
  1206. input_structure=pointindex,
  1207. centroid=tuple(pt),
  1208. volume=vol,
  1209. fragment=frag,
  1210. map_value=value,
  1211. time=None,
  1212. )
  1213. )
  1214. else:
  1215. logger.debug(
  1216. f"Assigning uncertain coordinate {tuple(pt)} to {len(self)} maps."
  1217. )
  1218. kernel = create_gaussian_kernel(sigma_vox, 3)
  1219. r = int(kernel.shape[0] / 2) # effective radius
  1220. assert pt.homogeneous.shape[0] == 1
  1221. xyz_vox = (np.dot(phys2vox, pt.homogeneous.T) + 0.5).astype("int")
  1222. shift = np.identity(4)
  1223. shift[:3, -1] = xyz_vox[:3, 0] - r
  1224. # build niftiimage with the Gaussian blob,
  1225. # then recurse into this method with the image input
  1226. gaussian_kernel = _volume.from_array(
  1227. data=kernel,
  1228. affine=np.dot(self.affine, shift),
  1229. space=self.space,
  1230. name=f"Gaussian kernel of {pt}",
  1231. cache=False,
  1232. )
  1233. for entry in self._assign(
  1234. item=gaussian_kernel,
  1235. lower_threshold=lower_threshold,
  1236. split_components=False
  1237. ):
  1238. entry.input_structure = pointindex
  1239. entry.centroid = tuple(pt)
  1240. assignments.append(entry)
  1241. return assignments
  1242. def _assign_volume(
  1243. self,
  1244. queryvolume: "_volume.Volume",
  1245. lower_threshold: float,
  1246. split_components: bool = True,
  1247. time: int = None,
  1248. **kwargs,
  1249. ) -> List[AssignImageResult]:
  1250. """
  1251. Assign an image volume to this parcellation map.
  1252. Parameters
  1253. -----------
  1254. queryvolume: Volume
  1255. the volume to be compared with maps
  1256. minsize_voxel: int, default: 1
  1257. Minimum voxel size of image components to be taken into account.
  1258. lower_threshold: float, default: 0
  1259. Lower threshold on values in the statistical map. Values smaller than
  1260. this threshold will be excluded from the assignment computation.
  1261. split_components: bool, default: True
  1262. Whether to split the query volume into disjoint components.
  1263. """
  1264. # TODO: split_components is not known to `assign`
  1265. # TODO: `minsize_voxel` is not used here. Consider the implementation of `assign` again.
  1266. if kwargs:
  1267. logger.info(f"The keywords {[k for k in kwargs]} are not passed on during volume assignment.")
  1268. if queryvolume.space != self.space:
  1269. raise ValueError("Assigned volume must be in the same space as the map.")
  1270. if split_components:
  1271. iter_components = lambda arr: connected_components(arr)
  1272. else:
  1273. iter_components = lambda arr: [(0, arr)]
  1274. queryimg = queryvolume.fetch()
  1275. assignments = []
  1276. all_indices = [
  1277. index
  1278. for regionindices in self._indices.values()
  1279. for index in regionindices
  1280. ]
  1281. with QUIET and provider.SubvolumeProvider.UseCaching():
  1282. for index in siibra_tqdm(
  1283. all_indices,
  1284. desc=f"Assigning {queryvolume} to {self}",
  1285. disable=len(all_indices) < 5,
  1286. unit="map",
  1287. leave=False
  1288. ):
  1289. region_map = self.fetch(index=index)
  1290. region_map_arr = np.asanyarray(region_map.dataobj)
  1291. # the shape and affine are checked by `nilearn.image.resample_to_img()`
  1292. # and returns the original data if resampling is not necessary.
  1293. queryimgarr_res = np.asanyarray(
  1294. resample_img_to_img(queryimg, region_map).dataobj
  1295. )
  1296. for compmode, voxelmask in iter_components(queryimgarr_res):
  1297. scores = compare_arrays(
  1298. voxelmask,
  1299. region_map.affine, # after resampling, both should have the same affine
  1300. region_map_arr,
  1301. region_map.affine
  1302. )
  1303. component_position = np.array(np.where(voxelmask)).T.mean(0)
  1304. if scores.intersection_over_union > lower_threshold:
  1305. assignments.append(
  1306. AssignImageResult(
  1307. input_structure=compmode,
  1308. centroid=tuple(component_position.round(2)),
  1309. volume=index.volume,
  1310. fragment=index.fragment,
  1311. map_value=index.label,
  1312. time=time,
  1313. **asdict(scores)
  1314. )
  1315. )
  1316. return assignments
  1317. def _assign_timeseries_volume(
  1318. self,
  1319. queryvolume: "_volume.TimeSeriesVolume",
  1320. lower_threshold: float,
  1321. split_components: bool = True,
  1322. **kwargs
  1323. ) -> List[AssignImageResult]:
  1324. assignments = []
  1325. for v_t in siibra_tqdm(queryvolume, unit='time point'):
  1326. assignments_t = self._assign_volume(
  1327. v_t,
  1328. lower_threshold=lower_threshold,
  1329. split_components=split_components,
  1330. time=v_t.timepoint,
  1331. **kwargs
  1332. )
  1333. assignments.extend(assignments_t)
  1334. return assignments
  1335. def to_BIDS_lookup_table(self) -> pd.DataFrame:
  1336. """
  1337. Generate a BIDS-compatible lookup table for the labelled map.
  1338. The lookup table associates voxel labels with region names and optional
  1339. RGB colors derived from the corresponding parcellation regions.
  1340. If the map consists of multiple fragments, the fragments are first
  1341. compressed into a single labelled image and relabelled to ensure BIDS
  1342. compatibility.
  1343. Parameters
  1344. ----------
  1345. filepath : str, optional
  1346. Path to a ``.tsv`` file where the lookup table should be written.
  1347. If provided, the table is saved in tab-separated format.
  1348. Returns
  1349. -------
  1350. pandas.DataFrame
  1351. A lookup table with the following columns:
  1352. - ``index``: Integer label value in the image.
  1353. - ``name``: Name of the corresponding brain region.
  1354. - ``color``: Hexadecimal RGB color code associated with the region,
  1355. or ``None`` if no color is defined.
  1356. Raises
  1357. ------
  1358. Exception
  1359. If the map contains more than one volume.
  1360. AssertionError
  1361. If ``filepath`` does not end with ``.tsv``.
  1362. Notes
  1363. -----
  1364. BIDS lookup tables require a single labelled volume. Maps distributed
  1365. across multiple fragments are therefore compressed and relabelled before
  1366. generating the table.
  1367. """
  1368. if not self.is_labelled:
  1369. raise NotImplementedError("Currently, there is not LUT standard defined by BIDS for statistical maps.")
  1370. if len(self.volumes) > 1 or not self.has_unique_labels:
  1371. logger.info(
  1372. f"{self} has {len(self.volumes)} volume(s)/{len(self.fragments)} fragment(s); "
  1373. "siibra will compress and reindex it for BIDS compatibility."
  1374. )
  1375. mp = self.compress()
  1376. else:
  1377. mp = self
  1378. def to_record(regionname: str, index: MapIndex) -> List[Dict]:
  1379. rgb = mp.parcellation.get_region(regionname).rgb
  1380. color = "#{:02x}{:02x}{:02x}".format(*rgb) if rgb else None
  1381. return {
  1382. "index": index.label,
  1383. "name": regionname,
  1384. "color": color,
  1385. }
  1386. table = pd.DataFrame(
  1387. [
  1388. to_record(r, indices[0])
  1389. for r, indices in mp._indices.items()
  1390. ]
  1391. )
  1392. if "gii-label" in mp.formats:
  1393. # read from the provider: mp.fetch() would also pull the template mesh
  1394. prov = mp.volumes[0]._providers["gii-label"]
  1395. observed = set()
  1396. for fragment in (mp.fragments or [None]):
  1397. observed |= set(np.unique(prov.fetch(fragment=fragment)["labels"]))
  1398. for xl in sorted(observed - set(table["index"]) - {0}): # 0 is background
  1399. table.loc[len(table)] = {"name": f"{xl} (unnamed)", "index": xl, "color": None}
  1400. if not table["index"].is_unique:
  1401. duplicated = sorted(table.loc[table["index"].duplicated(), "index"])
  1402. raise RuntimeError(
  1403. f"Labels {duplicated} are assigned to more than one region in {mp}. "
  1404. "A BIDS lookup table requires unique indices."
  1405. )
  1406. return table
  1407. def _as_surfaceimage(self, variant: str = None):
  1408. from nilearn.surface import SurfaceImage, PolyData
  1409. if "gii-label" not in self.formats:
  1410. raise ValueError("`SurfaceImage` representation is only possible for 'gii-label' maps.")
  1411. if len(self.volumes) > 1:
  1412. raise ValueError("`SurfaceImage` representation is only possible for maps with single and hemisphere fragemented maps.")
  1413. giilabel_filemap = {}
  1414. for frag in self.fragments:
  1415. loader = self.volumes[0]._providers["gii-label"]._loaders[frag]
  1416. loader._retrieve()
  1417. giilabel_filemap[frag.replace(' hemisphere', "")] = loader.cachefile
  1418. return SurfaceImage(mesh=self.space._as_polymesh(variant=variant), data=PolyData(**giilabel_filemap))
  1419. def as_nilearn_masker(
  1420. self,
  1421. strategy: Literal[
  1422. "mean",
  1423. "median",
  1424. "sum",
  1425. "minimum",
  1426. "maximum",
  1427. "standard_deviation",
  1428. "variance",
  1429. ] = "mean",
  1430. surface_variant: str = None,
  1431. **masker_kwargs,
  1432. ) -> Union["NiftiLabelsMasker", "SurfaceLabelsMasker"]:
  1433. from nilearn import maskers
  1434. try:
  1435. from ..retrieval.cache import jobmemory_path
  1436. masker_kwargs.setdefault("memory", jobmemory_path)
  1437. except ImportError:
  1438. ...
  1439. if not self.is_labelled:
  1440. raise NotImplementedError(
  1441. f"Nilearn maskers for {self.maptype} maps are provided by SparseMap, "
  1442. "which projects the data onto the maps instead of summarizing labelled "
  1443. "regions. Convert this map with `to_sparse()` first."
  1444. )
  1445. mp = self.compress() if (len(self.volumes) > 1 or self.fragments) else self
  1446. if "lut" in masker_kwargs:
  1447. raise ValueError("siibra handles `lut` parameter based on the map.")
  1448. masker_kwargs.setdefault("verbose", 1)
  1449. masker_kwargs.setdefault("strategy", strategy)
  1450. masker_kwargs["lut"] = mp.to_BIDS_lookup_table()
  1451. if self.provides_image:
  1452. masker = maskers.NiftiLabelsMasker(mp.fetch(), **masker_kwargs)
  1453. else:
  1454. masker = maskers.SurfaceLabelsMasker(
  1455. mp._as_surfaceimage(variant=surface_variant),
  1456. **masker_kwargs,
  1457. )
  1458. return masker
  1459. def extract_signals_with_nilearn(
  1460. self,
  1461. volume: _volume.Volume,
  1462. strategy: Literal[
  1463. "mean",
  1464. "median",
  1465. "sum",
  1466. "minimum",
  1467. "maximum",
  1468. "standard_deviation",
  1469. "variance",
  1470. ] = "mean",
  1471. confounds: np.ndarray = None,
  1472. sample_mask: np.ndarray = None,
  1473. surface_variant: str = None,
  1474. **masker_kwargs,
  1475. ) -> pd.DataFrame:
  1476. """
  1477. Extract region-wise signals from a 3D/4D volume using `nilearn`.
  1478. The regions defined in this map are used to summarize the input volume
  1479. with :class:`nilearn.maskers.NiftiLabelsMasker` (labelled volumetric maps),
  1480. :class:`nilearn.maskers.SurfaceLabelsMasker` (labelled surface maps), or
  1481. :class:`nilearn.maskers.NiftiMapsMasker` (statistical maps). Maps with
  1482. several volumes or fragments are compressed into a single labelled map first.
  1483. Parameters
  1484. ----------
  1485. volume: Volume
  1486. Input 3D or 4D volume from which signals should be extracted, typically
  1487. an fMRI image. Both NIfTI and timeseries GIFTI sources are supported.
  1488. strategy: str, default: "mean"
  1489. How voxels or vertices are summarized within each region. Only applies
  1490. to labelled maps; it is ignored for statistical maps, since
  1491. `NiftiMapsMasker` projects the data onto the (overlapping, continuous)
  1492. maps by least squares instead of summarizing discrete regions.
  1493. confounds: array-like, optional
  1494. Confounds to regress out during extraction.
  1495. See `nilearn.maskers.BaseMasker.transform`.
  1496. sample_mask: array-like, optional
  1497. Mask of samples to include when extracting signals.
  1498. See `nilearn.maskers.BaseMasker.transform`.
  1499. surface_variant: str, optional
  1500. Template surface variant to use for surface maps, e.g. "inflated".
  1501. **masker_kwargs
  1502. Passed on to the nilearn masker constructor.
  1503. Returns
  1504. -------
  1505. pandas.DataFrame
  1506. One column per region mapped in this map, one row per sample. The index
  1507. is the time axis of `volume` when it has one and no `sample_mask` was
  1508. given, otherwise a range index.
  1509. Notes
  1510. -----
  1511. Region names are taken from the BIDS lookup table generated by
  1512. `to_BIDS_lookup_table()`. Regions that disappear when the map is resampled
  1513. onto the input volume are reported and filled with zeros; nilearn's
  1514. `keep_masked_labels` does not prevent them from being dropped, and it is
  1515. deprecated since nilearn 0.14.
  1516. """
  1517. masker = self.as_nilearn_masker(
  1518. strategy=strategy if self.is_labelled else None,
  1519. surface_variant=surface_variant,
  1520. **masker_kwargs
  1521. )
  1522. if self.provides_image and volume.provides_image:
  1523. source = volume.fetch()
  1524. elif "gii-label" in self.formats and "gii-timeseries" in volume.formats:
  1525. source = volume._as_surfaceimage(variant=surface_variant)
  1526. else:
  1527. raise ValueError(
  1528. f"Cannot extract signals from {volume} with {self}: no common representation. "
  1529. f"The map provides {sorted(self.formats)}, the input provides {sorted(volume.formats)}."
  1530. )
  1531. # np.asarray normalizes plain and pandas output alike. (set_output(transform="pandas")
  1532. # raises NotImplementedError before nilearn 0.13, and the column names it
  1533. # produces are the ones we assign below anyway.)
  1534. signals = np.atleast_2d(np.asarray(
  1535. masker.fit_transform(source, confounds=confounds, sample_mask=sample_mask)
  1536. ))
  1537. if self.is_labelled:
  1538. # region_names_ maps output column index -> region name. It is available on
  1539. # both Nifti and Surface labels maskers since nilearn 0.10.4.
  1540. # (masker.labels_ cannot be used here: it includes the background label.
  1541. # get_feature_names_out() is equivalent but only exists from nilearn 0.13.)
  1542. extracted = [name for _, name in sorted(masker.region_names_.items())]
  1543. all_regions = masker.lut["name"].tolist()
  1544. else:
  1545. # NiftiMapsMasker names its columns positionally, so column i corresponds
  1546. # to volume i of the array stacked by _stack_maps().
  1547. extracted = all_regions = [
  1548. regionname for regionname, indices
  1549. in sorted(self._indices.items(), key=lambda kv: kv[1][0].volume)
  1550. ]
  1551. if signals.shape[1] != len(extracted):
  1552. raise RuntimeError(
  1553. f"nilearn returned {signals.shape[1]} signals but {len(extracted)} "
  1554. f"regions were expected for {self}."
  1555. )
  1556. result = pd.DataFrame(signals, columns=extracted)
  1557. missing = [r for r in all_regions if r not in set(extracted)]
  1558. if missing:
  1559. logger.info(
  1560. f"{len(missing)} region(s) were removed when resampling {self} to the "
  1561. f"input volume and are filled with zeros:\n{missing}"
  1562. )
  1563. result[missing] = 0
  1564. result = result[all_regions]
  1565. result.columns.name = "region"
  1566. # sample_mask drops samples, so the time axis would no longer align
  1567. time = getattr(volume, "time", None)
  1568. if time is not None and sample_mask is None and len(time) == len(result):
  1569. result.index = pd.Index(time, name="time")
  1570. return result
  1571. def from_volume(
  1572. name: str,
  1573. volume: Union[_volume.Volume, List[_volume.Volume]],
  1574. regionnames: List[str],
  1575. regionlabels: List[int],
  1576. parcellation_spec: Union[str, "parcellation.Parcellation"] = None
  1577. ) -> 'Map':
  1578. """
  1579. Add a custom labelled parcellation map to siibra from a labelled NIfTI file.
  1580. Parameters
  1581. ------------
  1582. name: str
  1583. Human-readable name of the parcellation.
  1584. volume: Volume, or a list of Volumes.
  1585. space_spec: str, Space
  1586. Specification of the reference space (space object, name, keyword, or id - e.g. 'mni152').
  1587. regionnames: list[str]
  1588. List of human-readable names of the mapped regions.
  1589. regionlabels: list[int]
  1590. List of integer labels in the nifti file corresponding to the list of regions.
  1591. parcellation: str or Parcellation. Optional.
  1592. If the related parcellation already defined or preconfigured in siibra.
  1593. """
  1594. # providers and map indices
  1595. providers = []
  1596. volumes = volume if isinstance(volume, list) else [volume]
  1597. map_space = volumes[0].space
  1598. assert all(v.space == map_space for v in volumes), "Volumes have to be in the same space"
  1599. for vol_idx, vol in enumerate(volumes):
  1600. image = vol.fetch()
  1601. arr = np.asanyarray(image.dataobj)
  1602. labels_in_volume = np.unique(arr)[1:].astype('int')
  1603. # populate region indices from given name/label lists
  1604. indices = dict()
  1605. for label, regionname in zip(regionlabels, regionnames):
  1606. if label not in labels_in_volume:
  1607. logger.warning(
  1608. f"Label {label} not mapped in the provided NIfTI volume -> "
  1609. f"region '{regionname} will not be in the map."
  1610. )
  1611. elif label in [v[0]['label'] for v in indices.values() if v[0]['volume'] == vol_idx]:
  1612. logger.warning(f"Label {label} already defined in the same volume; will not map it to '{regionname}'.")
  1613. else:
  1614. assert regionname not in indices, f"'{regionname}' must be unique in `regionnames`."
  1615. indices[regionname] = [{'volume': vol_idx, 'label': label}]
  1616. # check for any remaining labels in the NIfTI volume
  1617. unnamed_labels = list(set(labels_in_volume) - set(regionlabels))
  1618. if unnamed_labels:
  1619. logger.warning(
  1620. f"The following labels appear in the NIfTI volume {vol_idx}, but not in "
  1621. f"the specified regions: {', '.join(str(lb) for lb in unnamed_labels)}. "
  1622. "They will be removed from the nifti volume."
  1623. )
  1624. for label in unnamed_labels:
  1625. arr[arr == label] = 0
  1626. providers.extend(vol._providers.values())
  1627. # parcellation
  1628. if parcellation_spec is None:
  1629. parcellation_spec = name
  1630. try:
  1631. parcobj = parcellation.Parcellation.registry().get(parcellation_spec)
  1632. logger.info(f"Using '{parcellation_spec}', siibra decoded the parcellation as '{parcobj}'")
  1633. except Exception:
  1634. logger.info(
  1635. f"Using '{parcellation_spec}', siibra could not decode the "
  1636. " parcellation. Building a new parcellation."
  1637. )
  1638. # build a new parcellation
  1639. parcobj = parcellation.Parcellation(
  1640. identifier=generate_uuid(','.join(regionnames)),
  1641. name=name,
  1642. species=vol.space.species,
  1643. regions=list(map(_region.Region, regionnames)),
  1644. )
  1645. if parcobj.key not in list(parcellation.Parcellation.registry()):
  1646. parcellation.Parcellation.registry().add(parcobj.key, parcobj)
  1647. for region in siibra_tqdm(
  1648. indices.keys(),
  1649. desc="Checking if provided regions are defined in the parcellation."
  1650. ):
  1651. try:
  1652. _ = parcobj.get_region(region)
  1653. except Exception:
  1654. logger.warning(f"'{region}' is missing in the parcellation.")
  1655. # build the parcellation map object
  1656. parcmap = Map(
  1657. identifier=generate_uuid(name),
  1658. name=f"{name} map in {map_space.name}",
  1659. space_spec={"@id": map_space.id},
  1660. parcellation_spec={'name': parcobj.name},
  1661. indices=indices,
  1662. volumes=volumes
  1663. )
  1664. # add it to siibra's registry
  1665. Map.registry().add(parcmap.key, parcmap)
  1666. # return the map - note that it has a pointer to the parcellation
  1667. return parcmap

parcellationmap.py at commit 0a6a1e8, under Apache-2.0 · at the source

Overview

Authors: Timo Dickscheid1,2,3, Xiaoyun Gui1, Ahmet N Simsek1, Christian Schiffer1,3, Jean-Francois Mangin4, Yann Leprince4, Viktor Jirsa5, Jan G Bjaalie6, Trygve B Leergaard6, Sebastian Bludau1, Katrin Amunts1,7
  1. Institute of Neuroscience and Medicine (INM-1), Research Centre Jülich, Jülich, Germany
  2. Institute of Computational Visualistics, University of Koblenz, Koblenz, Germany
  3. Helmholtz AI, Research Centre Jülich, Jülich, Germany
  4. NeuroSpin, CEA, Université Paris-Saclay, Gif-sur-Yvette, France
  5. Institut de Neurosciences des Systèmes, INSERM, Aix-Marseille University, Marseille, France
  6. Institute of Basic Medical Sciences, University of Oslo, Oslo, Norway
  7. Medical Faculty & University Hospital Düsseldorf - Cécile & Oskar Vogt Institute of Brain Research, Heinrich Heine University, Düsseldorf, Germany
Journal: Nature methods, volume 23, issue 8, pages 1647-1658
Dates: received 18 March 2025; accepted 15 June 2026; published online 20 July 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41592-026-03159-x · PMID 42477446 · PMCID PMC13441889 · OpenAlex W4410759675
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), methods / tools (subfield)
Methods: Spectral & time-frequency, fMRI & imaging
Keywords: Computational neuroscience, Data integration, Software, Neuroscience
MeSH: Atlases as Topic*, Brain*, Brain Mapping*, Software*, Humans, Image Processing, Computer-Assisted, Internet, Magnetic Resonance Imaging, User-Computer Interface (* major topic)
Topic: Functional Brain Connectivity Studies (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: EC | Horizon 2020 Framework Programme (EU Framework Programme for Research and Innovation H2020) (945539, 101058516); Helmholtz Association (InterLabs-0015, Portfolio Theme SMHB, ZT-I-PF-4-061); Deutsche Forschungsgemeinschaft (German Research Foundation) (501864659, SPP2041)
Citations: cited by 1 paper (Europe PMC); 76 references in the paper
Research resources: RRID:SCR_001362, RRID:SCR_002498, RRID:SCR_002577, Matplotlib RRID:SCR_008624, as NumPy RRID:SCR_008633, stored in the EBRAINS Knowledge Graph RRID:SCR_017612, fetched tabular data as pandas RRID:SCR_018214, napari RRID:SCR_022765, RRID:SCR_023173, RRID:SCR_023498

Abstract

Computational technology opens new possibilities toward understanding the complexity of the human brain, but it requires integrating measurements from different modalities and scales in an anatomical context and exposing them in an interoperable, actionable form. Especially with growing big data resources, accessing information from different scales and modalities coherently for visual exploration, reproducible analysis and application development remains challenging. Here we present siibra, a tool suite that connects diverse data from cloud resources to reference atlases and coordinate spaces. It supports different use cases by making contents accessible through a web viewer, Python library and HTTP application programming interface. Using siibra we implemented a Multilevel Human Brain Atlas linking macro-anatomical concepts and their inter-subject variability with measurements of the microstructural composition and intrinsic variance of brain regions, building on cytoarchitecture as a reference and supporting MRI-based and microscopic templates. The atlas is integrated with the EBRAINS research infrastructure. All software and content are openly accessible.

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 25 matches between paragraphs and lines of code.

siibra-python.readthedocs.io

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: the link answers
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)

HumanBrainProject/hbp-spatial-backend

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 0acac3bc230fb8a65e6cf1a3f9196648606da65f, 1 April 2025
Languages: Python (13), Shell (2)
Size: 54 files, 15 scripts
Software Heritage: not archived
Found in: the text, “Modeling spatial entities in different reference”
Holds: README, license file, environment (Dockerfile, Dockerfile.server, setup.cfg, setup.py, docker-aims/Dockerfile), tests, continuous integration
Not found: CITATION.cff, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
17 files

FZJ-INM1-BDA/siibra-tutorials

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 88580425bf6815235725df94c5246e7f05583787, 22 June 2026
Languages: Jupyter (12)
Size: 20 files, 12 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (requirements.txt, runtime.txt), continuous integration, 12 notebooks
Not found: CITATION.cff, tests, documentation
Tools: Nilearn (12 files), Matplotlib (10 files), NumPy (2 files), pandas (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
14 files

FZJ-INM1-BDA/siibra-python

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 0a6a1e81fb007ed580816a7b963d8d9f89bb46b7, 26 September 2026
Languages: Python (168)
Size: 267 files, 168 scripts
Software Heritage: archived
Found in: “Code availability”
Holds: README, license file, CITATION.cff, environment (Dockerfile, requirements-test.txt, requirements.txt, setup.py, config_schema/requirements.txt, docs/requirements-doc.txt), tests, continuous integration, documentation
Tools: NumPy (45 files), Nilearn (29 files), pandas (20 files), NiBabel (15 files), Matplotlib (14 files), Plotly (6 files), scikit-image (6 files), h5py (2 files), SciPy (2 files), seaborn (2 files), scikit-learn (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
170 files

FZJ-INM1-BDA/siibra-explorer

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: a1b93e091f2457ff8d683dab27001b89332e29f9, 26 July 2026
Languages: TypeScript (574), JavaScript (23), Python (22), Shell (2)
Size: 1,304 files, 621 scripts
Software Heritage: archived
Found in: “Code availability”
Holds: README, license file, environment (Dockerfile, backend/requirements.txt, docs/Dockerfile, docs/requirements.txt, docs/develop/docker-compose.yml), tests, continuous integration, documentation
Not found: CITATION.cff
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
623 files

FZJ-INM1-BDA/siibra-api

License: Apache-2.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 3af86fe2c1f9b7532203d05f47eb4f682b8108b7, 23 September 2026
Languages: Python (489), Shell (2)
Size: 549 files, 491 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (docker-compose.yml, Dockerfile, server.dockerfile, worker-v4.dockerfile, worker.dockerfile, .docker/docker-compose.yaml, requirements/all.txt, requirements/docs.txt, requirements/partial-all.txt, requirements/partial-worker.txt, requirements/server.txt, requirements/siibra.txt), tests, continuous integration, documentation
Not found: CITATION.cff
Tools: NumPy (7 files), NiBabel (2 files), pandas (2 files), Plotly (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
493 files

apache.org/licenses

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: the link answers
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
At the source: apache.org/licenses/

Zenodo 21127648

License: CC-BY-4.0
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file, 0 scripts
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)

Code availability

Source codes for siibra are available on GitHub (https://github.com/FZJ-INM1-BDA/siibra-python; https://github.com/FZJ-INM1-BDA/siibra-explorer and https://github.com/FZJ-INM1-BDA/siibra-api) and provided under the Apache 2.0 license (http://www.apache.org/licenses/). We have integrated all components of siibra with Zenodo, which generates unique DOIs for each new release of the tool suite. An installation of the interactive viewer is publicly accessible at https://atlases.ebrains.eu/viewer as a core component of the brain atlas services in the research infrastructure EBRAINS. The web API can be accessed at https://siibra-api.apps.ebrains.eu/v3_0/docs. Researchers can install siibra-python through the Python Package Index from https://pypi.org/project/siibra. Comprehensive online documentation, including executable codes to reproduce all figures can be accessed via https://siibra-python.readthedocs.io/en/v1. More detailed tutorials notebooks are provided at https://github.com/FZJ-INM1-BDA/siibra-tutorials.

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;
  • 1,307 scripts, each with its path and the digest of its content;
  • 25 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 Statement

All data used in this work are publicly available via the URLs and references provided with the respective figures. Most contents are published as curated datasets with permanent identifiers on the EBRAINS platform (https://ebrains.eu). Gene expression levels were retrieved from the Allen Human Brain microarray data18,19. The current foundational content linked with the atlas is defined and maintained at https://github.com/FZJ-INM1-BDA/siibra-configurations/tree/v1. An overview is provided as Extended Data Table 2.

Source codes for siibra are available on GitHub (https://github.com/FZJ-INM1-BDA/siibra-python; https://github.com/FZJ-INM1-BDA/siibra-explorer and https://github.com/FZJ-INM1-BDA/siibra-api) and provided under the Apache 2.0 license (http://www.apache.org/licenses/). We have integrated all components of siibra with Zenodo, which generates unique DOIs for each new release of the tool suite. An installation of the interactive viewer is publicly accessible at https://atlases.ebrains.eu/viewer as a core component of the brain atlas services in the research infrastructure EBRAINS. The web API can be accessed at https://siibra-api.apps.ebrains.eu/v3_0/docs. Researchers can install siibra-python through the Python Package Index from https://pypi.org/project/siibra. Comprehensive online documentation, including executable codes to reproduce all figures can be accessed via https://siibra-python.readthedocs.io/en/v1. More detailed tutorials notebooks are provided at https://github.com/FZJ-INM1-BDA/siibra-tutorials.

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 2, 28 September 2026

  • Publisher: n/a → Nature Portfolio

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 11 authors, 4 keywords, 9 MeSH terms, 3 funders, 60 references, 10 RRIDs.

Cite

This paper

Dickscheid, T., Gui, X., Simsek, A. N., Schiffer, C., Mangin, J.-F., Leprince, Y., Jirsa, V., Bjaalie, J. G., Leergaard, T. B., Bludau, S., & Amunts, K. (2026). Siibra: a software tool suite for realizing a Multilevel Human Brain Atlas from complex data resources. Nature methods, 23(8), 1647-1658. https://doi.org/10.1038/s41592-026-03159-x

BibTeX

@article{dickscheid2026siibra,
author = {Dickscheid, Timo and Gui, Xiaoyun and Simsek, Ahmet N and Schiffer, Christian and Mangin, Jean-Francois and Leprince, Yann and Jirsa, Viktor and Bjaalie, Jan G and Leergaard, Trygve B and Bludau, Sebastian and Amunts, Katrin},
title = {{Siibra: a software tool suite for realizing a Multilevel Human Brain Atlas from complex data resources}},
journal = {Nature methods},
year = {2026},
month = jul,
volume = {23},
number = {8},
pages = {1647--1658},
publisher = {Nature Portfolio},
issn = {1548-7091},
doi = {10.1038/s41592-026-03159-x},
url = {https://doi.org/10.1038/s41592-026-03159-x},
pmid = {42477446},
pmcid = {PMC13441889}
}

RIS

TY - JOUR
AU - Dickscheid, Timo
AU - Gui, Xiaoyun
AU - Simsek, Ahmet N
AU - Schiffer, Christian
AU - Mangin, Jean-Francois
AU - Leprince, Yann
AU - Jirsa, Viktor
AU - Bjaalie, Jan G
AU - Leergaard, Trygve B
AU - Bludau, Sebastian
AU - Amunts, Katrin
TI - Siibra: a software tool suite for realizing a Multilevel Human Brain Atlas from complex data resources
T2 - Nature methods
J2 - Nat Methods
PY - 2026
DA - 2026/07/20
VL - 23
IS - 8
SP - 1647
EP - 1658
SN - 1548-7091
PB - Nature Portfolio
DO - 10.1038/s41592-026-03159-x
UR - https://doi.org/10.1038/s41592-026-03159-x
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41592-026-03159-x",
"type": "article-journal",
"title": "Siibra: a software tool suite for realizing a Multilevel Human Brain Atlas from complex data resources",
"container-title": "Nature methods",
"author": [
{
"family": "Dickscheid",
"given": "Timo"
},
{
"family": "Gui",
"given": "Xiaoyun"
},
{
"family": "Simsek",
"given": "Ahmet N"
},
{
"family": "Schiffer",
"given": "Christian"
},
{
"family": "Mangin",
"given": "Jean-Francois"
},
{
"family": "Leprince",
"given": "Yann"
},
{
"family": "Jirsa",
"given": "Viktor"
},
{
"family": "Bjaalie",
"given": "Jan G"
},
{
"family": "Leergaard",
"given": "Trygve B"
},
{
"family": "Bludau",
"given": "Sebastian"
},
{
"family": "Amunts",
"given": "Katrin"
}
],
"container-title-short": "Nat Methods",
"volume": "23",
"issue": "8",
"page": "1647-1658",
"DOI": "10.1038/s41592-026-03159-x",
"PMID": "42477446",
"PMCID": "PMC13441889",
"ISSN": "1548-7091",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s41592-026-03159-x",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
20
]
]
}
}

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.1162/imag.a.1260
Organization, fine structure, and stereotaxic maps of the human Bed nucleus of the Stria terminalis.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: 3 references, 2 authors
[2] doi:10.1038/s41467-026-71428-6 [code]
Binding items to contexts through conjunctive neural representations with the method of loci.
Journal: Nature communications
In common: Nilearn, Plotly, h5py, 7 other tools, 3 references
[3] doi:10.1038/s41467-026-76812-w [code]
Assessing molecular, cellular and transcriptomic bases of laminar perfusion and cytoarchitecture coupling in the human cortex.
Journal: Nature communications
In common: Nilearn, NiBabel, scikit-learn, 4 other tools, structural MRI / diffusion, 5 references
[4] doi:10.1162/imag.a.1198 [code]
MEPrep: A robust pipeline for multi-echo fMRI denoising and preprocessing.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, scikit-image, h5py, 6 other tools, methods / tools, 3 references
[5] doi:10.1038/s41597-026-06869-1 [code]
Individual Brain Charting: fifth release of high-resolution fMRI data for cognitive mapping.
Journal: Scientific data
In common: Nilearn, scikit-image, NiBabel, 6 other tools, methods / tools, 3 references
[6] doi:10.1162/imag.a.1269 [code]
From early to contemporary normative modeling: Mapping individual differences in neurophysiological signals.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: Nilearn, Plotly, NiBabel, 6 other tools, 3 references
[7] doi:10.1162/imag.a.1154 [code]
Evolutionary signatures in deep white matter architecture: A comparative study of humans and chimpanzees.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: scikit-learn, pandas, SciPy, 2 other tools, 2 references, author Yann Leprince
[8] doi:10.1016/j.isci.2026.117180 [code]
Developmental changes in similarity between neural representations of mental arithmetic and artificial neural networks.
Journal: iScience
In common: Nilearn, h5py, NiBabel, 6 other tools, 3 references
[9] doi:10.3389/frai.2026.1771088 [code]
Few-shot deployment of pretrained MRI transformers in brain imaging tasks.
Journal: Frontiers in artificial intelligence
In common: scikit-image, h5py, NiBabel, 6 other tools, methods / tools, structural MRI / diffusion, 2 references
[10] doi:10.1038/s41467-026-75837-5 [code]
A hierarchical framework for cortical and subcortical gray-matter parcellation across rodents, primates, and humans.
Journal: Nature communications
In common: NiBabel, seaborn, pandas, 3 other tools, methods / tools, 4 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.

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: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

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