[ENH] Efforts towards keeping memory low

docs/workflows.rst

@@ -180,7 +180,7 @@ packages, including FreeSurfer and the `Connectome Workbench`_.
 
 BOLD preprocessing
 ------------------
-:mod:`fmriprep.workflows.bold.init_func_preproc_wf`
+:mod:`fmriprep.workflows.bold.base.init_func_preproc_wf`
 
 .. workflow::
     :graph2use: orig
@@ -215,7 +215,7 @@ Preprocessing of BOLD files is split into multiple sub-workflows decribed below.
 
 BOLD reference image estimation
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-:mod:`fmriprep.workflows.bold.init_bold_reference_wf`
+:mod:`fmriprep.workflows.bold.util.init_bold_reference_wf`
 
 .. workflow::
     :graph2use: orig
@@ -243,7 +243,7 @@ Skullstripping of the reference image is performed using Nilearn.
 
 Slice time correction
 ~~~~~~~~~~~~~~~~~~~~~
-:mod:`fmriprep.workflows.bold.init_bold_stc_wf`
+:mod:`fmriprep.workflows.bold.stc.init_bold_stc_wf`
 
 .. workflow::
     :graph2use: colored
@@ -270,14 +270,16 @@ correction will be disabled.
 
 Head-motion estimation
 ~~~~~~~~~~~~~~~~~~~~~~
-:mod:`fmriprep.workflows.bold.init_bold_hmc_wf`
+:mod:`fmriprep.workflows.bold.hmc.init_bold_hmc_wf`
 
 .. workflow::
     :graph2use: colored
     :simple_form: yes
 
     from fmriprep.workflows.bold import init_bold_hmc_wf
-    wf = init_bold_hmc_wf(bold_file_size_gb=3, omp_nthreads=1)
+    wf = init_bold_hmc_wf(
+        mem_gb=1,
+        omp_nthreads=1)
 
 FSL MCFLIRT is used to estimate motion transformations using an automatically
 estimated reference scan (see :ref:`bold_ref`).
@@ -291,22 +293,35 @@ Susceptibility Distortion Correction (SDC)
     :show-inheritance:
 
 
+Pre-processed BOLD in native space
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+:mod:`fmriprep.workflows.bold.resampling.init_bold_preproc_trans_wf`
+
+.. workflow::
+    :graph2use: colored
+    :simple_form: yes
+
+    from fmriprep.workflows.bold import init_bold_preproc_trans_wf
+    wf = init_bold_preproc_trans_wf(mem_gb=3, omp_nthreads=1)
+
+
 .. _bold_reg:
 
 EPI to T1w registration
 ~~~~~~~~~~~~~~~~~~~~~~~
-:mod:`fmriprep.workflows.bold.init_bold_reg_wf`
+:mod:`fmriprep.workflows.bold.registration.init_bold_reg_wf`
 
 .. workflow::
     :graph2use: orig
     :simple_form: yes
 
     from fmriprep.workflows.bold import init_bold_reg_wf
-    wf = init_bold_reg_wf(freesurfer=True,
-                          bold_file_size_gb=3,
-                          omp_nthreads=1,
-                          use_bbr=True,
-                          bold2t1w_dof=9)
+    wf = init_bold_reg_wf(
+        freesurfer=True,
+        mem_gb=1,
+        omp_nthreads=1,
+        use_bbr=True,
+        bold2t1w_dof=9)
 
 The reference EPI image of each run is aligned by the ``bbregister`` routine to the
 reconstructed subject using the gray/white matter boundary (FreeSurfer's
@@ -323,17 +338,18 @@ boundary.
 
 EPI to MNI transformation
 ~~~~~~~~~~~~~~~~~~~~~~~~~
-:mod:`fmriprep.workflows.bold.init_bold_mni_trans_wf`
+:mod:`fmriprep.workflows.bold.resampling.init_bold_mni_trans_wf`
 
 .. workflow::
     :graph2use: colored
     :simple_form: yes
 
     from fmriprep.workflows.bold import init_bold_mni_trans_wf
-    wf = init_bold_mni_trans_wf(template='MNI152NLin2009cAsym',
-                                bold_file_size_gb=3,
-                                omp_nthreads=1,
-                                output_grid_ref=None)
+    wf = init_bold_mni_trans_wf(
+        template='MNI152NLin2009cAsym',
+        mem_gb=1,
+        omp_nthreads=1,
+        output_grid_ref=None)
 
 This sub-workflow concatenates the transforms calculated upstream (see
 `Head-motion estimation`_, `Susceptibility Distortion Correction (SDC)`_ (if
@@ -347,17 +363,18 @@ step, so as little information is lost as possible.
 
 EPI sampled to FreeSurfer surfaces
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-:mod:`fmriprep.workflows.bold.init_bold_surf_wf`
+:mod:`fmriprep.workflows.bold.resampling.init_bold_surf_wf`
 
 .. workflow::
     :graph2use: colored
     :simple_form: yes
 
     from fmriprep.workflows.bold import init_bold_surf_wf
-    wf = init_bold_surf_wf(bold_file_size_gb=0.1,
-                           output_spaces=['T1w', 'fsnative',
-                                         'template', 'fsaverage5'],
-                           medial_surface_nan=False)
+    wf = init_bold_surf_wf(
+        mem_gb=1,
+        output_spaces=['T1w', 'fsnative',
+                       'template', 'fsaverage5'],
+        medial_surface_nan=False)
 
 If FreeSurfer processing is enabled, the motion-corrected functional series
 (after single shot resampling to T1w space) is sampled to the
@@ -380,7 +397,9 @@ Confounds estimation
     from fmriprep.workflows.bold.confounds import init_bold_confs_wf
     wf = init_bold_confs_wf(
         name="discover_wf",
-        use_aroma=False, ignore_aroma_err=False, bold_file_size_gb=3,
+        use_aroma=False,
+        ignore_aroma_err=False,
+        mem_gb=1,
         metadata={"RepetitionTime": 2.0,
                   "SliceTiming": [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]})
 

fmriprep/interfaces/__init__.py

@@ -13,7 +13,7 @@
 )
 from .surf import NormalizeSurf, GiftiNameSource, GiftiSetAnatomicalStructure
 from .reports import SubjectSummary, FunctionalSummary, AboutSummary
-from .utils import TPM2ROI, ConcatROIs, CombineROIs, AddTSVHeader
+from .utils import TPM2ROI, AddTPMs, AddTSVHeader
 from .fmap import FieldEnhance
 from .confounds import GatherConfounds, ICAConfounds
 from .itk import MCFLIRT2ITK, MultiApplyTransforms

fmriprep/interfaces/utils.py

@@ -14,17 +14,19 @@
 import scipy.ndimage as nd
 from nilearn.image import resample_to_img
 
+from niworkflows.nipype import logging
 from niworkflows.nipype.utils.filemanip import fname_presuffix
 from niworkflows.nipype.interfaces.base import (
     traits, isdefined, TraitedSpec, BaseInterfaceInputSpec, File, InputMultiPath,
     SimpleInterface
 )
 
+IFLOGGER = logging.getLogger('interfaces')
+
 
 class TPM2ROIInputSpec(BaseInterfaceInputSpec):
-    t1_tpm = File(exists=True, mandatory=True, desc='Tissue probability map file in T1 space')
-    t1_mask = File(exists=True, mandatory=True, desc='Binary mask of skull-stripped T1w image')
-    bold_mask = File(exists=True, mandatory=True, desc='Binary mask of skull-stripped BOLD image')
+    in_tpm = File(exists=True, mandatory=True, desc='Tissue probability map file in T1 space')
+    in_mask = File(exists=True, mandatory=True, desc='Binary mask of skull-stripped T1w image')
     mask_erode_mm = traits.Float(xor=['mask_erode_prop'],
                                  desc='erode input mask (kernel width in mm)')
     erode_mm = traits.Float(xor=['erode_prop'],
@@ -47,13 +49,10 @@ class TPM2ROI(SimpleInterface):
 
     This interface follows the following logic:
 
-    #. Erode ``t1_mask`` by ``mask_erode_mm`` and apply to ``t1_tpm``
+    #. Erode ``in_mask`` by ``mask_erode_mm`` and apply to ``in_tpm``
     #. Threshold masked TPM at ``prob_thresh``
     #. Erode resulting mask by ``erode_mm``
 
-    Both the eroded brain mask and eroded ROI mask are then resampled to BOLD
-    resolution, and masked by ``bold_mask``.
-
     """
 
     input_spec = TPM2ROIInputSpec
@@ -73,9 +72,8 @@ def _run_interface(self, runtime):
         if not isdefined(erode_prop):
             erode_prop = None
         roi_file, eroded_mask = _tpm2roi(
-            self.inputs.t1_tpm,
-            self.inputs.t1_mask,
-            self.inputs.bold_mask,
+            self.inputs.in_tpm,
+            self.inputs.in_mask,
             mask_erode_mm,
             erode_mm,
             mask_erode_prop,
@@ -87,49 +85,49 @@ def _run_interface(self, runtime):
         return runtime
 
 
-class CombineROIsInputSpec(BaseInterfaceInputSpec):
+class AddTPMsInputSpec(BaseInterfaceInputSpec):
     in_files = InputMultiPath(File(exists=True), mandatory=True, desc='input list of ROIs')
-    ref_header = File(exists=True, mandatory=True,
-                      desc='reference NIfTI file with desired output header/affine')
+    indexes = traits.List(traits.Int, desc='select specific maps')
 
 
-class CombineROIsOutputSpec(TraitedSpec):
+class AddTPMsOutputSpec(TraitedSpec):
     out_file = File(exists=True, desc='union of binarized input files')
 
 
-class CombineROIs(SimpleInterface):
+class AddTPMs(SimpleInterface):
     """Generate the union (logical or) of a series of ROI masks"""
-    input_spec = CombineROIsInputSpec
-    output_spec = CombineROIsOutputSpec
+    input_spec = AddTPMsInputSpec
+    output_spec = AddTPMsOutputSpec
 
     def _run_interface(self, runtime):
-        self._results['out_file'] = _combine_rois(self.inputs.in_files, self.inputs.ref_header)
-        return runtime
+        in_files = self.inputs.in_files
 
+        indexes = list(range(len(in_files)))
+        if isdefined(self.inputs.indexes):
+            indexes = self.inputs.indexes
 
-class ConcatROIsInputSpec(BaseInterfaceInputSpec):
-    in_file = File(exists=True, mandatory=True, desc='input file')
-    in_mask = File(exists=True, mandatory=True, desc='input mask')
-    ref_header = File(exists=True, mandatory=True,
-                      desc='reference NIfTI file with desired output header/affine')
+        if len(self.inputs.in_files) < 2:
+            self._results['out_file'] = in_files[0]
+            return runtime
 
+        first_fname = in_files[indexes[0]]
+        if len(indexes) == 1:
+            self._results['out_file'] = first_fname
+            return runtime
 
-class ConcatROIsOutputSpec(TraitedSpec):
-    out_file = File(exists=True, desc='output average file')
+        im = nb.load(first_fname)
+        data = im.get_data()
 
+        for idx in indexes[1:]:
+            data += nb.load(in_files[idx]).get_data()
 
-class ConcatROIs(SimpleInterface):
-    """Concatenate two ROI files along time axis
+        data[data > 1] = 1.0
+        data[data < 0] = 0.0
 
-    ``in_file`` is resampled to match ``in_mask``, and the two concatenated
-    datasets are saved with the affine and header of ``ref_header``.
-    """
-    input_spec = ConcatROIsInputSpec
-    output_spec = ConcatROIsOutputSpec
+        out_file = fname_presuffix(first_fname, suffix='_tpmsum')
+        im.__class__(data, im.affine, im.header).to_filename(out_file)
+        self._results['out_file'] = out_file
 
-    def _run_interface(self, runtime):
-        self._results['out_file'] = _concat_rois(
-            self.inputs.in_file, self.inputs.in_mask, self.inputs.ref_header)
         return runtime
 
 
@@ -196,40 +194,36 @@ def _run_interface(self, runtime):
         return runtime
 
 
-def _tpm2roi(t1_tpm, t1_mask, bold_mask, mask_erosion_mm=None, erosion_mm=None,
+def _tpm2roi(in_tpm, in_mask, mask_erosion_mm=None, erosion_mm=None,
              mask_erosion_prop=None, erosion_prop=None, pthres=0.95):
     """
     Generate a mask from a tissue probability map
     """
-    tpm_img = nb.load(t1_tpm)
-    bold_mask_img = nb.load(bold_mask)
-
-    probability_map_data = (tpm_img.get_data() >= pthres).astype(np.uint8)
+    tpm_img = nb.load(in_tpm)
+    roi_mask = (tpm_img.get_data() >= pthres).astype(np.uint8)
 
-    eroded_mask_file = bold_mask
+    eroded_mask_file = None
     erode_in = (mask_erosion_mm is not None and mask_erosion_mm > 0 or
                 mask_erosion_prop is not None and mask_erosion_prop < 1)
     if erode_in:
-        eroded_mask_file = os.path.abspath("eroded_mask.nii.gz")
-
-        mask_img = nb.load(t1_mask)
+        eroded_mask_file = fname_presuffix(in_mask, suffix='_eroded')
+        mask_img = nb.load(in_mask)
+        mask_data = mask_img.get_data().astype(np.uint8)
         if mask_erosion_mm:
             iter_n = max(int(mask_erosion_mm / max(mask_img.header.get_zooms())), 1)
-            mask_data = nd.binary_erosion(mask_img.get_data().astype(np.uint8), iterations=iter_n)
+            mask_data = nd.binary_erosion(mask_data, iterations=iter_n)
         else:
-            mask_data = mask_img.get_data().astype(np.uint8)
             orig_vol = np.sum(mask_data > 0)
             while np.sum(mask_data > 0) / orig_vol > mask_erosion_prop:
                 mask_data = nd.binary_erosion(mask_data, iterations=1)
 
-        # Resample to BOLD and save
-        eroded_t1 = nb.Nifti1Image(mask_data, mask_img.affine, mask_img.header)
-        eroded_t1.set_data_dtype(np.uint8)
-        eroded_bold = _resample_and_mask(eroded_t1, bold_mask_img, interpolation='nearest')
-        eroded_bold.to_filename(eroded_mask_file)
+        # Store mask
+        eroded = nb.Nifti1Image(mask_data, mask_img.affine, mask_img.header)
+        eroded.set_data_dtype(np.uint8)
+        eroded.to_filename(eroded_mask_file)
 
         # Mask TPM data (no effect if not eroded)
-        probability_map_data[~mask_data] = 0
+        roi_mask[~mask_data] = 0
 
     # shrinking
     erode_out = (erosion_mm is not None and erosion_mm > 0 or
@@ -238,19 +232,18 @@ def _tpm2roi(t1_tpm, t1_mask, bold_mask, mask_erosion_mm=None, erosion_mm=None,
         if erosion_mm:
             iter_n = max(int(erosion_mm / max(tpm_img.header.get_zooms())), 1)
             iter_n = int(erosion_mm / max(tpm_img.header.get_zooms()))
-            probability_map_data = nd.binary_erosion(probability_map_data, iterations=iter_n)
+            roi_mask = nd.binary_erosion(roi_mask, iterations=iter_n)
         else:
-            orig_vol = np.sum(probability_map_data > 0)
-            while np.sum(probability_map_data > 0) / orig_vol > erosion_prop:
-                probability_map_data = nd.binary_erosion(probability_map_data, iterations=1)
+            orig_vol = np.sum(roi_mask > 0)
+            while np.sum(roi_mask > 0) / orig_vol > erosion_prop:
+                roi_mask = nd.binary_erosion(roi_mask, iterations=1)
 
     # Create image to resample
-    img_t1 = nb.Nifti1Image(probability_map_data, tpm_img.affine, tpm_img.header)
-    img_t1.set_data_dtype(np.uint8)
-    roi_img = _resample_and_mask(img_t1, bold_mask_img, interpolation='nearest')
-
-    roi_img.to_filename("roi.nii.gz")
-    return os.path.abspath("roi.nii.gz"), eroded_mask_file
+    roi_fname = fname_presuffix(in_tpm, suffix='_roi')
+    roi_img = nb.Nifti1Image(roi_mask, tpm_img.affine, tpm_img.header)
+    roi_img.set_data_dtype(np.uint8)
+    roi_img.to_filename(roi_fname)
+    return roi_fname, eroded_mask_file or in_mask
 
 
 def _resample_and_mask(source_img, mask_img, interpolation='continuous'):

fmriprep/workflows/base.py

@@ -445,6 +445,6 @@ def init_single_subject_wf(subject_id, task_id, name,
                'inputnode.t1_2_fsnative_forward_transform'),
               ('outputnode.t1_2_fsnative_reverse_transform',
                'inputnode.t1_2_fsnative_reverse_transform')]),
-            ])
+        ])
 
     return workflow

fmriprep/workflows/bold/__init__.py

@@ -26,6 +26,7 @@
 from .resampling import (
     init_bold_mni_trans_wf,
     init_bold_surf_wf,
+    init_bold_preproc_trans_wf,
 )
 
 from .confounds import (