Merge pull request #380 from MICA-MNI/fix/issues-274-341-343

Implement #274 (histology bit-flip), #343 (Destrieux atlas), #341 (resilient Neurosynth fetch)

Author: zihuaihuai

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(

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]:

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<NASGU>
+        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');