ENH: improve pet mask generation

Author: mnoergaard

fmriprep/workflows/pet/confounds.py

@@ -580,12 +580,12 @@ def _binary_union(mask1, mask2):
 
 
 def _smooth_binarize(in_file, fwhm=10.0, thresh=0.2):
-    """Smooth ``in_file`` with a Gaussian kernel and binarize at ``thresh``."""
+    """Smooth ``in_file`` with a Gaussian kernel, binarize and keep largest cluster."""
     from pathlib import Path
 
     import nibabel as nb
     import numpy as np
-    from scipy.ndimage import gaussian_filter
+    from scipy.ndimage import gaussian_filter, label
 
     img = nb.load(in_file)
     data = img.get_fdata(dtype=np.float32)
@@ -594,6 +594,12 @@ def _smooth_binarize(in_file, fwhm=10.0, thresh=0.2):
     smoothed = gaussian_filter(data, sigma=sigma)
     mask = smoothed > (thresh * smoothed.max())
 
+    labeled, n_labels = label(mask)
+    if n_labels > 1:
+        sizes = np.bincount(labeled.ravel())
+        sizes[0] = 0  # ignore background
+        mask = labeled == sizes.argmax()
+
     out_img = img.__class__(mask.astype('uint8'), img.affine, img.header)
     out_img.set_data_dtype('uint8')
     out_name = Path('smoothed_bin_mask.nii.gz').absolute()

fmriprep/workflows/pet/tests/test_smooth_binarize.py

@@ -0,0 +1,19 @@
+import nibabel as nb
+import numpy as np
+from scipy.ndimage import label
+
+from ..confounds import _smooth_binarize
+
+
+def test_smooth_binarize_largest(tmp_path):
+    data = np.zeros((5, 5, 5))
+    data[1:3, 1:3, 1:3] = 1
+    data[4, 4, 4] = 1
+    img = nb.Nifti1Image(data, np.eye(4))
+    src = tmp_path / 'input.nii.gz'
+    img.to_filename(src)
+
+    out = _smooth_binarize(str(src), fwhm=0.0, thresh=0.5)
+    result = nb.load(out).get_fdata()
+    _, num = label(result > 0)
+    assert num == 1