From ccf5af816e60603edd7a3998524181dfade5339b Mon Sep 17 00:00:00 2001 From: Enning Yang Date: Mon, 4 May 2026 13:30:03 -0400 Subject: [PATCH] Implement issues #274, #343 (Destrieux), #341 #274: Invert BigBrain histology profiles by default read_histology_profile now flips the 8-bit intensity values so that higher numbers correspond to darker, more cell-dense regions, matching the convention used in Paquola et al. 2021 and most downstream BigBrain analyses (the raw BigBrainWarp data uses the opposite convention). Pass `invert=False` to keep the previous behaviour. MPC and gradient computations are correlation-based and therefore invariant to the flip; only direct inspection of profile intensities is affected. #343: Add the Destrieux atlas to fetch_parcellation `fetch_parcellation(template, "destrieux", n_regions, ...)` now returns the FreeSurfer aparc.a2009s parcellation (75 cortical labels per hemisphere). Sourced from nilearn.datasets.fetch_atlas_surf_destrieux on fsaverage5; for other fsaverage* templates we resample with nearest- neighbour interpolation. fslr32k is rejected (FreeSurfer-specific). n_regions is accepted for API consistency but ignored. This is the closest off-the-shelf FreeSurfer atlas BrainStat can ship without hosting new data files; a hosted aparc/Desikan-Killiany option will need a separate download URL and is left as a follow-up. #341: Resilient Neurosynth download meta_analytic_decoder>fetch_neurosynth_data is now retry-safe: - The downloaded zip is validated via java.util.zip.ZipFile before unzip, so truncated downloads surface a clear error instead of MATLAB's terse "Invalid ZIP file" later on. - websave + verify + unzip is wrapped in a small retry loop (3 attempts) and any partial zip is cleaned up between attempts. - After extraction we verify the expected file count and surface a remediation hint if it doesn't match. - find_files now handles a missing data_dir gracefully. --- brainstat/context/histology.py | 29 ++++++-- brainstat/datasets/base.py | 51 ++++++++++++- .../context/meta_analytic_decoder.m | 73 +++++++++++++++++-- 3 files changed, 140 insertions(+), 13 deletions(-) diff --git a/brainstat/context/histology.py b/brainstat/context/histology.py index 095ae6f5..e70ca4d3 100644 --- a/brainstat/context/histology.py +++ b/brainstat/context/histology.py @@ -114,6 +114,7 @@ def read_histology_profile( data_dir: Optional[Union[str, Path]] = None, template: str = "fsaverage", overwrite: bool = False, + invert: bool = True, ) -> np.ndarray: """Reads BigBrain histology profiles. @@ -127,6 +128,15 @@ def read_histology_profile( default 'fsaverage'. overwrite : bool, optional If true, existing data will be overwrriten, by default False. + invert : bool, optional + If True (default), invert the 8-bit intensity values so that higher + numbers correspond to darker, more cell-dense regions, matching the + convention used in Paquola et al. 2021 (eLife) and most downstream + BigBrain analyses. The raw BigBrainWarp data uses the opposite + convention (0 = black). Set to False to return the raw values. MPC + and gradient computations (correlation-based) are invariant to this + flip; only direct inspection of profile intensities is affected. + See https://github.com/MICA-MNI/BrainStat/issues/274 for context. Returns ------- @@ -169,12 +179,19 @@ def read_histology_profile( profiles = h5_file.get("fs_LR_64k")[...] else: profiles = h5_file.get(template)[...] - if civet_template: - fsaverage_surface = fetch_template_surface("fsaverage") - civet_surface = fetch_template_surface(civet_template) - return _surf2surf(fsaverage_surface, civet_surface, profiles.T).T - else: - return profiles + + if invert: + # Map from BigBrainWarp convention (low = dark/dense) to the + # literature convention (high = dark/dense). The data is unsigned + # 8-bit, so the inversion is `255 - x`. + max_value = np.iinfo(profiles.dtype).max if np.issubdtype(profiles.dtype, np.integer) else 255 + profiles = max_value - profiles + + if civet_template: + fsaverage_surface = fetch_template_surface("fsaverage") + civet_surface = fetch_template_surface(civet_template) + return _surf2surf(fsaverage_surface, civet_surface, profiles.T).T + return profiles def download_histology_profiles( diff --git a/brainstat/datasets/base.py b/brainstat/datasets/base.py index f63118de..af8232e1 100644 --- a/brainstat/datasets/base.py +++ b/brainstat/datasets/base.py @@ -111,12 +111,15 @@ def fetch_parcellation( The surface template. Valid values are "fsaverage", "fsaverage5", "fsaverage6", "fslr32k", "civet41k", "civet164k", by default "fsaverage5". atlas : str - Name of the atlas. Valid names are "cammoun", "glasser", "schaefer", "yeo". + Name of the atlas. Valid names are "cammoun", "destrieux", "glasser", + "schaefer", "yeo". n_regions : int Number of regions of the requested atlas. Valid values for the cammoun atlas are 33, 60, 125, 250, 500. Valid values for the glasser atlas are 360. Valid values for the "schaefer" atlas are 100, 200, 300, 400, 500, - 600, 800, 1000. Valid values for "yeo" are 7 and 17. + 600, 800, 1000. Valid values for "yeo" are 7 and 17. The + "destrieux" atlas (a2009s, FreeSurfer's 75 cortical labels per + hemisphere) ignores ``n_regions``. join : bool, optional If true, returns parcellation as a single array, if false, returns an array per hemisphere, by default True. @@ -154,6 +157,8 @@ def fetch_parcellation( parcellations = _fetch_glasser_parcellation(template, data_dir) elif atlas == "yeo": parcellations = _fetch_yeo_parcellation(template, n_regions, data_dir) + elif atlas == "destrieux": + parcellations = _fetch_destrieux_parcellation(template, data_dir) else: raise ValueError(f"Invalid atlas: {atlas}") @@ -537,6 +542,48 @@ def _fetch_glasser_parcellation(template: str, data_dir: Path) -> List[np.ndarra return parcellations +def _fetch_destrieux_parcellation( + template: str, data_dir: Path +) -> List[np.ndarray]: + """Fetches the Destrieux 2009 (FreeSurfer ``aparc.a2009s``) parcellation. + + Uses :func:`nilearn.datasets.fetch_atlas_surf_destrieux`, which provides + the parcellation on the ``fsaverage5`` surface. For other ``fsaverage*`` + templates we resample with nearest-neighbour interpolation; ``fslr32k`` + is not supported because the destrieux atlas is FreeSurfer-specific. + """ + from nilearn.datasets import fetch_atlas_surf_destrieux + + if template == "fslr32k": + raise ValueError( + "The destrieux atlas is FreeSurfer-specific; fslr32k is not " + "supported. Use a fsaverage* template instead." + ) + + bunch = fetch_atlas_surf_destrieux(data_dir=str(data_dir)) + # nilearn returns int32 arrays for fsaverage5 (10242 vertices/hemi). + parcellations = [np.asarray(bunch["map_left"]), np.asarray(bunch["map_right"])] + + if template != "fsaverage5": + # Nearest-neighbour resample fsaverage5 -> requested fsaverage*. + from brainstat.mesh.interpolate import _surf2surf + + src_left, src_right = fetch_template_surface( + "fsaverage5", layer="white", join=False + ) + dst_left, dst_right = fetch_template_surface( + template, layer="white", join=False + ) + parcellations[0] = _surf2surf( + src_left, dst_left, parcellations[0], interpolation="nearest" + ).astype(parcellations[0].dtype) + parcellations[1] = _surf2surf( + src_right, dst_right, parcellations[1], interpolation="nearest" + ).astype(parcellations[1].dtype) + + return parcellations + + def _fetch_yeo_parcellation( template: str, n_regions: int, data_dir: Path ) -> List[np.ndarray]: diff --git a/brainstat_matlab/context/meta_analytic_decoder.m b/brainstat_matlab/context/meta_analytic_decoder.m index a177b391..b911cfc5 100644 --- a/brainstat_matlab/context/meta_analytic_decoder.m +++ b/brainstat_matlab/context/meta_analytic_decoder.m @@ -79,29 +79,92 @@ function neurosynth_files = fetch_neurosynth_data(data_dir) % Fetches neurosynth data files. +% +% Downloads the precomputed Neurosynth archive and extracts it into +% data_dir. Robust to interrupted downloads (which previously left a +% truncated .zip on disk and broke subsequent calls; see #341): the zip is +% verified before unzip, and the download is retried up to a small number +% of times. The zip and any partial extracted state are cleaned up if the +% download or extraction fails. json = brainstat_utils.read_data_fetcher_json(); neurosynth_files = find_files(data_dir); if numel(neurosynth_files) ~= json.neurosynth_precomputed.n_files + if ~exist(data_dir, 'dir'); mkdir(data_dir); end disp('Downloading Neurosynth files. This may take several minutes.') - zip_file = tempname(data_dir) + ".zip"; - cleanup = onCleanup(@() delete(zip_file)); - websave(zip_file, json.neurosynth_precomputed.url); - unzip(zip_file, data_dir) + + max_attempts = 3; + last_err = []; + for attempt = 1:max_attempts + zip_file = tempname(data_dir) + ".zip"; + cleanup = onCleanup(@() delete_if_exists(zip_file)); %#ok + try + websave(zip_file, json.neurosynth_precomputed.url); + verify_zip_integrity(zip_file); + unzip(zip_file, data_dir); + last_err = []; + break + catch ME + last_err = ME; + warning('BrainStat:NeurosynthDownload', ... + 'Attempt %d/%d failed: %s', attempt, max_attempts, ME.message); + end + end + + if ~isempty(last_err) + rethrow(last_err); + end + neurosynth_files = find_files(data_dir); + if numel(neurosynth_files) ~= json.neurosynth_precomputed.n_files + error('BrainStat:NeurosynthDownload', ... + ['Neurosynth archive extracted but %d files were found (expected %d). ' ... + 'Try removing %s and re-running.'], ... + numel(neurosynth_files), json.neurosynth_precomputed.n_files, data_dir); + end end function neurosynth_files = find_files(data_dir) + if ~exist(data_dir, 'dir') + neurosynth_files = strings(0, 1); + return + end data_dir_contents = dir(data_dir); data_dir_filenames = {data_dir_contents.name}; neurosynth_files = regexp(data_dir_filenames, "Neurosynth_TFIDF.*_z_.*consistency.nii.gz", ... 'match', 'once'); - neurosynth_files(cellfun(@isempty, neurosynth_files)) = []; + neurosynth_files(cellfun(@isempty, neurosynth_files)) = []; neurosynth_files = data_dir + filesep + neurosynth_files; end end +function verify_zip_integrity(zip_file) +% Sanity-checks a downloaded zip before handing it to MATLAB's unzip, +% which is unhelpfully terse when given a truncated archive. We use +% Java's ZipFile to validate the central directory without extracting +% the (multi-GB) contents. +info = dir(zip_file); +if isempty(info) || info.bytes < 1024 + error('BrainStat:NeurosynthDownload', ... + 'Downloaded zip is missing or implausibly small (%d bytes).', ... + max(0, info.bytes)); +end +try + zf = java.util.zip.ZipFile(zip_file); + zf.close(); +catch ME + error('BrainStat:NeurosynthDownload', ... + 'Downloaded zip failed integrity check: %s', ME.message); +end +end + +function delete_if_exists(filepath) +if exist(filepath, 'file') + delete(filepath); +end +end + function interpolated_volume = nearest_neighbor_interpolation(template, labels) % Performs nearest neighbor interpolation using precomputed volumes. precomputed_data_dir = brainstat_utils.get_brainstat_directories('brainstat_precomputed_data');