[ENH/WIP] Robust EPI Mask for heavily nonuniform image

.circleci/config.yml

@@ -260,7 +260,13 @@ jobs:
               -e FMRIPREP_REGRESSION_REPORTS=/tmp/data/reports \
               --entrypoint="py.test" poldracklab/fmriprep:latest \
               /root/src/fmriprep/ \
-              --doctest-modules --ignore=/root/src/fmriprep/docs --ignore=setup.py
+              -svx --doctest-modules --ignore=/root/src/fmriprep/docs --ignore=setup.py
+      - run:
+          name: Package new masks
+          no_output_timeout: 10m
+          working_directory: /tmp/data/reports
+          command: |
+            tar cfz fmriprep_bold_mask.tar.gz fmriprep_bold_mask/*/*.nii.gz
       - run:
           name: Test fmriprep-wrapper (Python 2)
           command: |

fmriprep/interfaces/nilearn.py

@@ -82,7 +82,7 @@ def _run_interface(self, runtime):
 
         if self.inputs.fill_holes:
             filled = binary_fill_holes(masknii.get_data().astype(
-                np.uint8), sim.ball(6)).astype(np.uint8)
+                np.uint8), sim.ball(1)).astype(np.uint8)
             masknii = masknii.__class__(filled, masknii.affine,
                                         masknii.header)
 

fmriprep/workflows/bold/tests/test_util.py

@@ -1,10 +1,11 @@
 ''' Testing module for fmriprep.workflows.bold.util '''
 import pytest
 import os
+from pathlib import Path
 
 import numpy as np
 from nipype.pipeline import engine as pe
-from nipype.utils.filemanip import fname_presuffix
+from nipype.utils.filemanip import fname_presuffix, copyfile
 from nilearn.image import load_img
 
 from niworkflows.interfaces.masks import ROIsPlot
@@ -54,25 +55,33 @@ def test_masking(input_fname, expected_fname):
 
     # Reconstruct base_fname from above
     dirname, basename = os.path.split(input_fname)
-    newpath = os.path.join(os.getenv('FMRIPREP_REGRESSION_REPORTS', '.'),
-                           os.path.basename(dirname))
+    dsname = os.path.basename(dirname)
+    reports_dir = Path(os.getenv('FMRIPREP_REGRESSION_REPORTS', ''))
+    newpath = reports_dir / dirname
     out_fname = fname_presuffix(basename, suffix='_masks.svg', use_ext=False,
-                                newpath=newpath)
-    os.makedirs(newpath, exist_ok=True)
+                                newpath=str(newpath))
+    newpath.mkdir(parents=True, exist_ok=True)
 
     mask_diff_plot = pe.Node(ROIsPlot(), name='mask_diff_plot')
     mask_diff_plot.inputs.in_mask = expected_fname
     mask_diff_plot.inputs.out_report = out_fname
 
+    outputnode = bold_reference_wf.get_node('outputnode')
     bold_reference_wf.connect([
-        (bold_reference_wf.get_node('outputnode'), mask_diff_plot, [
-            ('ref_image', 'in_file'),
-            ('bold_mask', 'in_rois'),
-            ])])
+        (outputnode, mask_diff_plot, [('ref_image', 'in_file'),
+                                      ('bold_mask', 'in_rois')])
+        ])
     res = bold_reference_wf.run(plugin='MultiProc')
 
-    combine_masks = [node for node in res.nodes if node.name.endswith('combine_masks')][0]
+    combine_masks = [node for node in res.nodes if node.name.endswith('combine_masks2')][0]
     overlap = symmetric_overlap(expected_fname,
                                 combine_masks.result.outputs.out_file)
 
     assert overlap > 0.95, input_fname
+
+    mask_dir = reports_dir / 'fmriprep_bold_mask' / dsname
+    mask_dir.mkdir(parents=True, exist_ok=True)
+    copyfile(combine_masks.result.outputs.out_file,
+             fname_presuffix(basename, suffix='_mask',
+                             use_ext=True, newpath=str(mask_dir)),
+             copy=True)

fmriprep/workflows/bold/util.py

@@ -148,19 +148,19 @@ def init_enhance_and_skullstrip_bold_wf(name='enhance_and_skullstrip_bold_wf',
 
     Steps of this workflow are:
 
-
-      1. Calculate a conservative mask using Nilearn's ``create_epi_mask``.
+      1. Calculate a loose mask using FSL's ``bet``.
       2. Run ANTs' ``N4BiasFieldCorrection`` on the input
          :abbr:`BOLD (blood-oxygen level-dependant)` average, using the
          mask generated in 1) instead of the internal Otsu thresholding.
-      3. Calculate a loose mask using FSL's ``bet``, with one mathematical morphology
-         dilation of one iteration and a sphere of 6mm as structuring element.
+      3. Calculate a mask using Nilearn's ``create_epi_mask`` on corrected image,
+         with one mathematical morphology dilation of one iteration and
+         a sphere of 3mm as structuring element.
       4. Mask the :abbr:`INU (intensity non-uniformity)`-corrected image
          with the latest mask calculated in 3), then use AFNI's ``3dUnifize``
          to *standardize* the T2* contrast distribution.
       5. Calculate a mask using AFNI's ``3dAutomask`` after the contrast
          enhancement of 4).
-      6. Calculate a final mask as the intersection of 3) and 5).
+      6. Calculate a final mask as the intersection of 1), 3) and 5).
       7. Apply final mask on the enhanced reference.
 
 
@@ -208,20 +208,21 @@ def init_enhance_and_skullstrip_bold_wf(name='enhance_and_skullstrip_bold_wf',
         'mask_file', 'skull_stripped_file', 'bias_corrected_file']), name='outputnode')
 
     # Create a loose mask to avoid N4 internal's Otsu mask
-    n4_mask = pe.Node(MaskEPI(upper_cutoff=0.75, enhance_t2=enhance_t2, opening=1,
-                      no_sanitize=True), name='n4_mask')
+    n4_mask = pe.Node(fsl.BET(frac=0.3, mask=True, robust=True),
+                      name='n4_mask')
+    fixhdr_n4_mask = pe.Node(CopyXForm(), name='fixhdr_n4_mask', mem_gb=0.1)
 
     # Run N4 normally, force num_threads=1 for stability (images are small, no need for >1)
     n4_correct = pe.Node(ants.N4BiasFieldCorrection(dimension=3, copy_header=True),
                          name='n4_correct', n_procs=1)
 
-    # Create a generous BET mask out of the bias-corrected EPI
-    skullstrip_first_pass = pe.Node(fsl.BET(frac=0.2, mask=True),
-                                    name='skullstrip_first_pass')
-    bet_dilate = pe.Node(fsl.DilateImage(
-        operation='max', kernel_shape='sphere', kernel_size=6.0,
+    # Create a generous mask out of the bias-corrected EPI
+    skullstrip_first_pass = pe.Node(MaskEPI(upper_cutoff=0.75, enhance_t2=enhance_t2, opening=1,
+                                    no_sanitize=True), name='skullstrip_first_pass')
+    dilate_mask_first_pass = pe.Node(fsl.DilateImage(
+        operation='max', kernel_shape='sphere', kernel_size=3.0,
         internal_datatype='char'), name='skullstrip_first_dilate')
-    bet_mask = pe.Node(fsl.ApplyMask(), name='skullstrip_first_mask')
+    apply_mask_first_pass = pe.Node(fsl.ApplyMask(), name='skullstrip_first_mask')
 
     # Use AFNI's unifize for T2 constrast & fix header
     unifize = pe.Node(afni.Unifize(
@@ -238,34 +239,40 @@ def init_enhance_and_skullstrip_bold_wf(name='enhance_and_skullstrip_bold_wf',
                                      name='skullstrip_second_pass')
     fixhdr_skullstrip2 = pe.Node(CopyXForm(), name='fixhdr_skullstrip2', mem_gb=0.1)
 
-    # Take intersection of both masks
-    combine_masks = pe.Node(fsl.BinaryMaths(operation='mul'),
-                            name='combine_masks')
+    # Take intersection of all 3 masks
+    combine_masks1 = pe.Node(fsl.BinaryMaths(operation='mul'),
+                             name='combine_masks1')
+    combine_masks2 = pe.Node(fsl.BinaryMaths(operation='mul'),
+                             name='combine_masks2')
 
     # Compute masked brain
     apply_mask = pe.Node(fsl.ApplyMask(), name='apply_mask')
 
     workflow.connect([
-        (inputnode, n4_mask, [('in_file', 'in_files')]),
+        (inputnode, n4_mask, [('in_file', 'in_file')]),
         (inputnode, n4_correct, [('in_file', 'input_image')]),
+        (inputnode, fixhdr_n4_mask, [('in_file', 'hdr_file')]),
         (inputnode, fixhdr_unifize, [('in_file', 'hdr_file')]),
         (inputnode, fixhdr_skullstrip2, [('in_file', 'hdr_file')]),
-        (n4_mask, n4_correct, [('out_mask', 'mask_image')]),
-        (n4_correct, skullstrip_first_pass, [('output_image', 'in_file')]),
-        (skullstrip_first_pass, bet_dilate, [('mask_file', 'in_file')]),
-        (bet_dilate, bet_mask, [('out_file', 'mask_file')]),
-        (skullstrip_first_pass, bet_mask, [('out_file', 'in_file')]),
-        (bet_mask, unifize, [('out_file', 'in_file')]),
+        (n4_mask, fixhdr_n4_mask, [('mask_file', 'in_file')]),
+        (fixhdr_n4_mask, n4_correct, [('out_file', 'mask_image')]),
+        (fixhdr_n4_mask, combine_masks2, [('out_file', 'operand_file')]),
+        (n4_correct, skullstrip_first_pass, [('output_image', 'in_files')]),
+        (n4_correct, apply_mask_first_pass, [('output_image', 'in_file')]),
+        (skullstrip_first_pass, dilate_mask_first_pass, [('out_mask', 'in_file')]),
+        (dilate_mask_first_pass, apply_mask_first_pass, [('out_file', 'mask_file')]),
+        (dilate_mask_first_pass, combine_masks1, [('out_file', 'in_file')]),
+        (apply_mask_first_pass, unifize, [('out_file', 'in_file')]),
         (unifize, fixhdr_unifize, [('out_file', 'in_file')]),
         (fixhdr_unifize, skullstrip_second_pass, [('out_file', 'in_file')]),
-        (skullstrip_first_pass, combine_masks, [('mask_file', 'in_file')]),
-        (skullstrip_second_pass, fixhdr_skullstrip2, [('out_file', 'in_file')]),
-        (fixhdr_skullstrip2, combine_masks, [('out_file', 'operand_file')]),
         (fixhdr_unifize, apply_mask, [('out_file', 'in_file')]),
-        (combine_masks, apply_mask, [('out_file', 'mask_file')]),
-        (combine_masks, outputnode, [('out_file', 'mask_file')]),
+        (skullstrip_second_pass, fixhdr_skullstrip2, [('out_file', 'in_file')]),
+        (fixhdr_skullstrip2, combine_masks1, [('out_file', 'operand_file')]),
+        (combine_masks1, combine_masks2, [('out_file', 'in_file')]),
+        (combine_masks2, apply_mask, [('out_file', 'mask_file')]),
+        (combine_masks2, outputnode, [('out_file', 'mask_file')]),
         (apply_mask, outputnode, [('out_file', 'skull_stripped_file')]),
-        (n4_correct, outputnode, [('output_image', 'bias_corrected_file')]),
+        (n4_correct, outputnode, [('output_image', 'bias_corrected_file')])
     ])
 
     return workflow