ENH Visualize HCP fs_LR 32k data and resample to fsaverage (cortex.hcp)#3
ENH Visualize HCP fs_LR 32k data and resample to fsaverage (cortex.hcp)#3mvdoc wants to merge 2 commits into
Conversation
Add a `cortex.hcp` module for working with HCP fs_LR 32k data in pycortex: - download_fs_lr(): fetch the prebuilt 32k_fs_LR pycortex subject (registered with download_subject) for native fs_LR rendering. - cifti_to_surface()/get_cifti_vertex_indices(): expand CIFTI grayordinates to the full 64984-vertex fs_LR surface, filling the medial wall with NaN. - project_fslr_to_fsaverage()/to_fsaverage(): resample fs_LR 32k data to fsaverage5/6/fsaverage via a pure-python spherical-barycentric matrix built from the HCP standard-mesh spheres. This matches `wb_command -metric-resample ... BARYCENTRIC`, so Connectome Workbench is not required at runtime. Also make upsample_to_fsaverage() FreeSurfer-free for the common cases: the per-vertex nearest-neighbor mapping is a fixed property of the fsaverage tessellation, so ship precomputed neighbor tables for fsaverage5/6 (0.68 MB) and fall back to computing from the sphere (and caching) for other orders. Output is numerically identical to the previous sphere-based implementation. Add examples/hcp/plot_hcp_fs_lr.py (renders both the native and the resampled flatmaps) and tests in cortex/tests/test_hcp.py, plus FreeSurfer-free upsample_to_fsaverage tests. The barycentric matrix is validated against wb_command when it is available. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019z37ML4wsmosg5Xvvq6ENN
![gemini-code-assist[bot]](https://avatars.githubusercontent.com/in/956858?s=60&v=4){kind=small}
There was a problem hiding this comment.
Code Review
This pull request introduces a new hcp module to support the visualization and resampling of HCP fs_LR 32k data, along with corresponding tests, examples, and integration updates across the codebase. The code review feedback focuses on enhancing the robustness and performance of this new module. Key recommendations include adding integrity checks for cached sphere files, improving parameter type validation for CIFTI inputs, supporting hemispherically split CIFTI files, ensuring compatibility with older SciPy versions by passing string paths to load_npz and save_npz, optimizing memory usage by preserving float32 precision, and removing a redundant np.abs call on non-negative matrices.
Important
The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.
| if fpath.exists(): | ||
| return fpath |
There was a problem hiding this comment.
If a previous download was interrupted or corrupted, fpath.exists() will still be True, and the function will return a corrupted file, leading to cryptic errors later when trying to load it.
We can add a quick integrity check by attempting to load the cached file with nib.load. If it fails, we log a warning, delete the corrupted file, and proceed to re-download it.
| if fpath.exists(): | |
| return fpath | |
| if fpath.exists(): | |
| try: | |
| nib.load(str(fpath)) | |
| return fpath | |
| except Exception: | |
| logger.warning("Cached sphere file %s is corrupted. Re-downloading...", fpath) | |
| fpath.unlink(missing_ok=True) |
| if not isinstance(cifti, nib.Cifti2Image): | ||
| cifti = nib.load(str(cifti)) |
There was a problem hiding this comment.
Using str(cifti) when cifti is not a Cifti2Image can lead to confusing errors. If cifti is already a loaded image of another type (e.g., a Nifti1Image passed by mistake), str(cifti) will return its string representation (like "<nibabel.nifti1.Nifti1Image object at 0x...>") rather than a file path, causing a FileNotFoundError when nib.load is called.
It is safer to explicitly check if cifti is a path-like object (a string or Path) before loading it, and raise a clear TypeError if it is neither a path nor a Cifti2Image.
| if not isinstance(cifti, nib.Cifti2Image): | |
| cifti = nib.load(str(cifti)) | |
| if isinstance(cifti, (str, Path)): | |
| cifti = nib.load(str(cifti)) | |
| if not isinstance(cifti, nib.Cifti2Image): | |
| raise TypeError("cifti must be a Cifti2Image object or a path to a CIFTI-2 file.") |
| if "L" not in indices or "R" not in indices: | ||
| raise ValueError("CIFTI file is missing a left and/or right cortex structure.") | ||
| return indices["L"], indices["R"] |
There was a problem hiding this comment.
Currently, get_cifti_vertex_indices raises a ValueError if either the left or right cortex structure is missing. However, some CIFTI files are hemispherically split and only contain one hemisphere.
We can make this function much more robust and compatible with hemispherically split files by defaulting any missing hemisphere to an empty array instead of raising an error. cifti_to_surface and project_fslr_to_fsaverage already handle empty index arrays gracefully.
| if "L" not in indices or "R" not in indices: | |
| raise ValueError("CIFTI file is missing a left and/or right cortex structure.") | |
| return indices["L"], indices["R"] | |
| if "L" not in indices and "R" not in indices: | |
| raise ValueError("CIFTI file is missing both left and right cortex structures.") | |
| left_indices = indices.get("L", np.array([], dtype=np.int64)) | |
| right_indices = indices.get("R", np.array([], dtype=np.int64)) | |
| return left_indices, right_indices |
| if cache and cache_path.exists(): | ||
| logger.info("Loading cached matrix from %s", cache_path) | ||
| return load_npz(cache_path) | ||
|
|
There was a problem hiding this comment.
In older versions of SciPy, load_npz does not support pathlib.Path objects and requires string paths. To ensure maximum compatibility across different SciPy versions, we should convert cache_path to a string.
| return load_npz(str(cache_path)) |
|
|
||
| if cache: | ||
| cache_path.parent.mkdir(parents=True, exist_ok=True) | ||
| save_npz(cache_path, matrix) |
There was a problem hiding this comment.
In older versions of SciPy, save_npz does not support pathlib.Path objects and requires string paths. To ensure maximum compatibility across different SciPy versions, we should convert cache_path to a string.
| save_npz(cache_path, matrix) | |
| save_npz(str(cache_path), matrix) |
| projected = (matrix @ hemi_clean.T).T | ||
|
|
||
| # Weight actually contributed by valid (non-NaN) source vertices. | ||
| valid_weight = (np.abs(matrix) @ (~nan_mask).astype(np.float64).T).T |
There was a problem hiding this comment.
The resampling matrix contains only non-negative barycentric weights (clamped to >= 0 and normalized to sum to 1) or fallback weights of 1.0. Therefore, taking the absolute value with np.abs(matrix) is completely redundant and adds unnecessary computation overhead on a large sparse matrix.
We can safely remove np.abs and use matrix directly.
| valid_weight = (np.abs(matrix) @ (~nan_mask).astype(np.float64).T).T | |
| valid_weight = (matrix @ (~nan_mask).astype(np.float64).T).T |
| Data on the target surface. Medial-wall and all-NaN-source vertices are | ||
| NaN; apply ``numpy.nan_to_num`` before handing to pycortex if needed. | ||
| """ | ||
| data = np.asarray(data, dtype=np.float64) |
There was a problem hiding this comment.
Forcing data to np.float64 can double the memory usage for large datasets (e.g., long timeseries), which are extremely common in neuroimaging and typically stored as float32.
We should preserve the input's floating-point precision (or default to float32 for integer inputs) to optimize memory usage and performance.
| data = np.asarray(data, dtype=np.float64) | |
| out_dtype = data.dtype if np.issubdtype(data.dtype, np.floating) else np.float32 | |
| data = np.asarray(data, dtype=out_dtype) |
| ) | ||
|
|
||
| n_tgt_hemi = N_VERTICES_TARGET_HEM[matrix_target] | ||
| result = np.full((data.shape[0], 2 * n_tgt_hemi), np.nan, dtype=np.float64) |
There was a problem hiding this comment.
Initialize result with out_dtype instead of forcing np.float64 to preserve the input's floating-point precision and save memory.
| result = np.full((data.shape[0], 2 * n_tgt_hemi), np.nan, dtype=np.float64) | |
| result = np.full((data.shape[0], 2 * n_tgt_hemi), np.nan, dtype=out_dtype) |
| if not isinstance(cifti, nib.Cifti2Image): | ||
| cifti = nib.load(str(cifti)) |
There was a problem hiding this comment.
Using str(cifti) when cifti is not a Cifti2Image can lead to confusing errors. If cifti is already a loaded image of another type (e.g., a Nifti1Image passed by mistake), str(cifti) will return its string representation (like "<nibabel.nifti1.Nifti1Image object at 0x...>") rather than a file path, causing a FileNotFoundError when nib.load is called.
It is safer to explicitly check if cifti is a path-like object (a string or Path) before loading it, and raise a clear TypeError if it is neither a path nor a Cifti2Image.
| if not isinstance(cifti, nib.Cifti2Image): | |
| cifti = nib.load(str(cifti)) | |
| if isinstance(cifti, (str, Path)): | |
| cifti = nib.load(str(cifti)) | |
| if not isinstance(cifti, nib.Cifti2Image): | |
| raise TypeError("cifti must be a Cifti2Image object or a path to a CIFTI-2 file.") |
- get_cifti_vertex_indices/to_fsaverage: raise TypeError for non-path, non-Cifti2Image inputs instead of a cryptic FileNotFoundError. - get_cifti_vertex_indices: tolerate hemispherically split CIFTI files by returning an empty index array for a missing hemisphere. - project_fslr_to_fsaverage: preserve float32 input precision instead of forcing float64 (halves memory for timeseries). - get_fslr_to_fsaverage_matrix: pass str paths to load_npz/save_npz for older SciPy compatibility. - _project_hemi: drop the redundant np.abs on the non-negative weight matrix. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_019z37ML4wsmosg5Xvvq6ENN
|
Superseded by gallantlab#655 (opened against upstream gallantlab:main). |
1 participant
Adds a
cortex.hcpmodule to visualize HCP fs_LR 32k data with pycortex — natively on the HCP template and resampled to fsaverage.What's new
cortex.hcpdownload_fs_lr()— fetch the prebuilt32k_fs_LRpycortex subject (registered withdownload_subject) for native fs_LR rendering via the usualcortex.Vertex(...)/quickshowpath.cifti_to_surface()/get_cifti_vertex_indices()— expand CIFTI cortical grayordinates (59412) to the full 64984-vertex fs_LR surface, NaN medial wall.project_fslr_to_fsaverage()/to_fsaverage()— resample fs_LR 32k data to fsaverage5/6/fsaverage using a pure-python spherical-barycentric matrix built from the HCP standard-mesh spheres. This matcheswb_command -metric-resample ... BARYCENTRIC, so Connectome Workbench is not required at runtime.FreeSurfer-free
upsample_to_fsaveragecortex/data/upsample_fsaverage_neighbors.npz, 0.68 MB), with a compute-from-sphere (+cache) fallback for other orders. Output is numerically identical to the previous implementation;$SUBJECTS_DIRis no longer needed for the common cases.Example + tests
examples/hcp/plot_hcp_fs_lr.pyrenders both the native and the resampled flatmaps (runs on CI with no FreeSurfer).cortex/tests/test_hcp.py(barycentric matrix math, CIFTI expansion,wb_commandvalidation when available) and new FreeSurfer-freeupsample_to_fsaveragetests.Notes
32k_fs_LRsubject is derived from HCP S1200 group-average Open Access surfaces and redistributed under the WU-Minn HCP Consortium Open Access Data Use Terms (with the required acknowledgment); its download URL is registered indownload_subject..npzships in the built wheel and loads from the installed package.🤖 Generated with Claude Code
https://claude.ai/code/session_019z37ML4wsmosg5Xvvq6ENN