Warn about coarse input voxel size and clarify resolution requirements

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 

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",

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

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 "
+                "(<lh|rh>.hippoAmygLabels*.mgz) rather than the version that was resampled to the conformed "
+                "space (<lh|rh>.hippoAmygLabels*.FSvoxelSpace.mgz)."
+            )
+
     # check hemisphere
 
     if params.HEMI != "lh" and params.HEMI != "rh":

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")
 

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: