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: