From 8876d945ec8daaf75c1b94ae048178df9e6bf26c Mon Sep 17 00:00:00 2001 From: Karl Koschutnig Date: Thu, 6 Aug 2026 12:26:27 +0200 Subject: [PATCH] Warn about coarse input voxel size and clarify resolution requirements Segmentations at ~1 mm rather than the ~0.33 mm this method is designed for fail checkSurface with "Surface contains holes", which sends users toward manual editing even though the segmentation is usually fine. The mask filters are parameterised in voxels, not mm, so their physical strength scales with the input voxel size. - Warn at startup when the input voxel edge length exceeds a threshold, and for --lut freesurfer additionally point out the .FSvoxelSpace. variant, which is the common way to hit this. - Add the voxel size to the "Surface contains holes" message when the input is coarse, so the error itself suggests checking the resolution before editing the segmentation. - Log when --upsample is a no-op, which it silently is for isotropic input, since reaching for it is the natural response to the error. - Document the two FreeSurfer output files and the resolution dependency in the tutorial and documentation. Refs #25 --- TUTORIAL.md | 28 +++++++++++++++++++++++++++- hipsta/cfg/config.py | 4 ++++ hipsta/doc/DOCUMENTATION.md | 6 +++++- hipsta/hipsta.py | 21 +++++++++++++++++++++ hipsta/processImage.py | 12 ++++++++++++ hipsta/utils/check_surface.py | 14 ++++++++++++++ 6 files changed, 83 insertions(+), 2 deletions(-) diff --git a/TUTORIAL.md b/TUTORIAL.md index 86e9ebd..15278e2 100644 --- a/TUTORIAL.md +++ b/TUTORIAL.md @@ -228,6 +228,17 @@ hemisphere in a FreeSurfer segmentation): run_hipsta --filename /path/to/my/segmentation/image --hemi lh --lut freesurfer --outputdir /path/to/my/output/directory ``` +Note that the FreeSurfer subfield module writes two copies of every +segmentation: one at its native resolution of about 0.33 mm +(e.g. `lh.hippoAmygLabels.mgz`), and one that has been resampled to the +conformed space of about 1 mm (e.g. `lh.hippoAmygLabels.FSvoxelSpace.mgz`). +Use the native-resolution file, i.e. the one **without** `FSvoxelSpace` in its +name. The resampled file is too coarse for this method and will frequently fail +the surface check with [holes in the surface](#holes-in-the-surface), even when +the underlying segmentation is perfectly fine. Hipsta will issue a warning at +startup if the voxel size of the input image is substantially larger than +0.33 mm. + ### Additional arguments for ASHS segmentations Segmentation images that are produced by the ASHS require a set of additional @@ -403,7 +414,22 @@ We outline two strategies to mitigate these issues below. #### Holes in the surface Holes in the surfaces can be the result of less than optimal image -preprocessing. We recommend to try adjusting the width and threshold for the +preprocessing. + +Before adjusting any processing parameters, check the voxel size of the input +image. The mask filtering operations are parametrised in voxels rather than in +mm, so their physical strength scales with the voxel size of the input. At +around 1 mm the hippocampal ribbon is only about two voxels thick, and the +default filtering can erode through it, which produces holes even for an +otherwise correct segmentation. If you are working with FreeSurfer +segmentations, make sure you are using the native-resolution file rather than +the `FSvoxelSpace` one (see [above](#mandatory-arguments)). Note that +resampling a coarse segmentation to a smaller voxel size after the fact does +not recover the lost detail, because nearest-neighbour interpolation preserves +the original staircase pattern. + +If the input resolution is appropriate, we recommend to try adjusting the width +and threshold for the gaussian filter using the `--gauss-filter-size` argument, and optionally to employ additional smoothing along the longitudinal axis using the `--long-filter` flag (which can be fine-tuned using the `--long-filter-size` argument). It may diff --git a/hipsta/cfg/config.py b/hipsta/cfg/config.py index 6c30ce4..cd872d4 100644 --- a/hipsta/cfg/config.py +++ b/hipsta/cfg/config.py @@ -44,6 +44,10 @@ def get_defaults(x): # deprecated options no_orient=False, # internal options (not set during parsing or class definition, but during args evaluation) + # voxel edge length (in mm) above which a warning about the input + # resolution will be issued; the method is designed for segmentations + # at around 0.33 mm, so this is set slightly above that value + voxel_size_threshold=0.4, map_values_integrate="none", map_values_select=None, map_values_interp="nearest", diff --git a/hipsta/doc/DOCUMENTATION.md b/hipsta/doc/DOCUMENTATION.md index 02c95fc..fdaa579 100644 --- a/hipsta/doc/DOCUMENTATION.md +++ b/hipsta/doc/DOCUMENTATION.md @@ -231,7 +231,11 @@ Also the intermediate volume and surface files can be useful for quality control ## Supported segmentations: - A hippocampal subfields segmentation based on FreeSurfer 7.1.1 or later should - work as-is. Use `--lut freesurfer`. + work as-is. Use `--lut freesurfer`. Use the native-resolution segmentation + (about 0.33 mm), not the one that was resampled to the conformed space, which + has `FSvoxelSpace` in its filename: the latter is too coarse for this method + and will frequently fail the surface check with holes even when the + segmentation itself is fine. - Two ASHS atlases are currently supported, the Penn ABC-3T ASHS atlas for T2-weighted MRI and the UMC Utrecht 7T atlas. Use `--lut ashs-penn_abc_3t_t2` or `--lut ashs-umcutrecht_7t`. If additional labels for the hippocampal head diff --git a/hipsta/hipsta.py b/hipsta/hipsta.py index 3d3d46e..143f565 100644 --- a/hipsta/hipsta.py +++ b/hipsta/hipsta.py @@ -11,6 +11,7 @@ import sys import time +import nibabel as nb import numpy as np from .cfg.atlases import get_atlases @@ -667,6 +668,26 @@ def _check_params(params): else: raise RuntimeError("Could not find " + params.FILENAME) + # check voxel size; this is stored on params so that later stages (e.g. + # checkSurface) can refer to it without re-reading the input image + + params.internal.VOXEL_SIZE = tuple(float(x) for x in nb.load(params.FILENAME).header.get_zooms()[0:3]) + + if max(params.internal.VOXEL_SIZE) > get_defaults("voxel_size_threshold"): + LOGGER.warning( + "The input image has a voxel size of %s mm. This method is designed for segmentations with a " + "voxel edge length of about 0.33 mm; with substantially coarser images, the surface check will " + "frequently report holes even though the segmentation itself is fine.", + " x ".join(format(x, ".3f") for x in params.internal.VOXEL_SIZE), + ) + + if params.LUT == "freesurfer": + LOGGER.warning( + "For FreeSurfer input, please make sure to use the native-resolution segmentation " + "(.hippoAmygLabels*.mgz) rather than the version that was resampled to the conformed " + "space (.hippoAmygLabels*.FSvoxelSpace.mgz)." + ) + # check hemisphere if params.HEMI != "lh" and params.HEMI != "rh": diff --git a/hipsta/processImage.py b/hipsta/processImage.py index fbaac43..08590d8 100644 --- a/hipsta/processImage.py +++ b/hipsta/processImage.py @@ -9,6 +9,7 @@ import subprocess import nibabel as nb +import numpy as np from nilearn import image as nli # ============================================================================== @@ -131,6 +132,17 @@ def upsampleImage(params): target_affn[0:4, 1] *= min(vxsz) / vxsz[1] target_affn[0:4, 2] *= min(vxsz) / vxsz[2] + # upsampling to the smallest voxel edge length has no effect if the + # image already has isotropic voxels; report this, as it is not + # otherwise apparent from the output + if np.allclose(target_affn, affn): + LOGGER.warning( + "The image already has isotropic voxels of %s mm, so upsampling to the smallest voxel " + "edge length leaves it unchanged. An explicit target size can be requested with the " + "--upsample-size argument.", + format(float(min(vxsz)), ".3f"), + ) + # resample img_int = nli.resample_img(img, target_affine=target_affn, interpolation="nearest") diff --git a/hipsta/utils/check_surface.py b/hipsta/utils/check_surface.py index a29973a..443025a 100644 --- a/hipsta/utils/check_surface.py +++ b/hipsta/utils/check_surface.py @@ -8,6 +8,8 @@ from lapy import TriaMesh +from ..cfg.config import get_defaults + # ============================================================================== # LOGGING @@ -36,6 +38,18 @@ def checkSurface(params, stage=None): if euler != 2: LOGGER.info("Surface contains holes. Please edit the corresponding hippocampal segmentation and re-run.") + + voxel_size = getattr(params.internal, "VOXEL_SIZE", None) + + if voxel_size is not None and max(voxel_size) > get_defaults("voxel_size_threshold"): + LOGGER.info( + "Note that the input image has a voxel size of %s mm, which is coarse for this method. " + "Holes are common at this resolution even for an otherwise correct segmentation, so " + "please check whether a higher-resolution version of the segmentation is available " + "before editing it.", + " x ".join(format(x, ".3f") for x in voxel_size), + ) + continue_program = False else: