Calculate basic CBF in whole FOV (#600)

Author: tsalo

aslprep/interfaces/cbf.py

@@ -95,7 +95,6 @@ class _ExtractCBFInputSpec(BaseInterfaceInputSpec):
         mandatory=True,
         desc="metadata for M0 scan. Only defined if M0Type is 'Separate'.",
     )
-    in_mask = File(exists=True, mandatory=True, desc='mask')
     fwhm = traits.Float(default_value=5, usedefault=True, mandatory=False, desc='fwhm')
 
 
@@ -129,12 +128,12 @@ def _run_interface(self, runtime):
         aslcontext = pd.read_table(self.inputs.aslcontext)
         metadata = self.inputs.metadata.copy()
 
-        mask_data = nb.load(self.inputs.in_mask).get_fdata()
-
         # read the preprocessed ASL data
         asl_img = nb.load(self.inputs.asl_file)
         asl_data = asl_img.get_fdata()
 
+        dims = asl_data.shape[:3]
+
         if aslcontext.shape[0] != asl_img.shape[3]:
             raise ValueError(
                 f'Number of rows in aslcontext ({aslcontext.shape[0]}) != '
@@ -152,11 +151,9 @@ def _run_interface(self, runtime):
         # extract m0 file and register it to ASL if separate
         if metadata['M0Type'] == 'Separate':
             m0file = self.inputs.m0scan
-            m0data_smooth = smooth_image(nb.load(m0file), fwhm=self.inputs.fwhm).get_fdata()
-            if len(m0data_smooth.shape) > 3:
-                m0data = mask_data * np.mean(m0data_smooth, axis=3)
-            else:
-                m0data = mask_data * m0data_smooth
+            m0data = smooth_image(nb.load(m0file), fwhm=self.inputs.fwhm).get_fdata()
+            if len(m0data.shape) > 3:
+                m0data = np.mean(m0data, axis=3)
 
             m0tr = self.inputs.m0scan_metadata['RepetitionTimePreparation']
             if np.array(m0tr).size > 1 and np.std(m0tr) > 0:
@@ -165,8 +162,8 @@ def _run_interface(self, runtime):
         elif metadata['M0Type'] == 'Included':
             m0data = asl_data[:, :, :, m0_volume_idx]
             m0img = nb.Nifti1Image(m0data, asl_img.affine, asl_img.header)
-            m0data_smooth = smooth_image(m0img, fwhm=self.inputs.fwhm).get_fdata()
-            m0data = mask_data * np.mean(m0data_smooth, axis=3)
+            m0data = smooth_image(m0img, fwhm=self.inputs.fwhm).get_fdata()
+            m0data = np.mean(m0data, axis=3)
 
             if np.array(metadata['RepetitionTimePreparation']).size > 1:
                 m0tr = np.array(metadata['RepetitionTimePreparation'])[m0_volume_idx]
@@ -177,7 +174,7 @@ def _run_interface(self, runtime):
                 raise ValueError('M0 scans have variable TR. ASLPrep does not support this.')
 
         elif metadata['M0Type'] == 'Estimate':
-            m0data = metadata['M0Estimate'] * mask_data
+            m0data = np.full(dims, metadata['M0Estimate'])
 
             m0tr = None
 
@@ -194,7 +191,7 @@ def _run_interface(self, runtime):
                 control_data = asl_data[:, :, :, control_volume_idx]
                 control_img = nb.Nifti1Image(control_data, asl_img.affine, asl_img.header)
                 control_img = smooth_image(control_img, fwhm=self.inputs.fwhm).get_fdata()
-                m0data = mask_data * np.mean(control_img, axis=3)
+                m0data = np.mean(control_img, axis=3)
 
                 # Use the control volumes' TR as the M0 TR.
                 if np.array(metadata['RepetitionTimePreparation']).size > 1:
@@ -204,7 +201,7 @@ def _run_interface(self, runtime):
 
             elif cbf_volume_idx:
                 # If we have precalculated CBF data, we don't need M0, so we'll just use the mask.
-                m0data = mask_data
+                m0data = np.ones(dims)
 
                 m0tr = None
 
@@ -314,7 +311,6 @@ class _ComputeCBFInputSpec(BaseInterfaceInputSpec):
     )
     m0_file = File(exists=True, mandatory=True, desc='M0 nifti file')
     m0tr = traits.Float(mandatory=False, desc='M0 TR, in seconds.')
-    mask = File(exists=True, mandatory=True, desc='Mask nifti file')
     cbf_only = traits.Bool(
         mandatory=True,
         desc='Whether data are deltam (False) or CBF (True).',
@@ -392,9 +388,12 @@ def _run_interface(self, runtime):
         metadata = self.inputs.metadata
         m0_file = self.inputs.m0_file
         m0_scale = self.inputs.m0_scale
-        mask_file = self.inputs.mask
         deltam_file = self.inputs.deltam  # control - label signal intensities
 
+        deltam_img = nb.load(deltam_file)
+        dims = deltam_img.shape[:3]
+        mask_img = nb.Nifti1Image(np.ones(dims), deltam_img.affine, deltam_img.header)
+
         if self.inputs.cbf_only:
             config.loggers.interface.debug('CBF data detected. Skipping CBF estimation.')
             self._results['cbf_ts'] = fname_presuffix(
@@ -434,7 +433,7 @@ def _run_interface(self, runtime):
 
         # NOTE: Nilearn will still add a singleton time dimension for 3D imgs with
         # NiftiMasker.transform, until 0.12.0, so the arrays will currently be 2D no matter what.
-        masker = maskers.NiftiMasker(mask_img=mask_file)
+        masker = maskers.NiftiMasker(mask_img=mask_img)
         deltam_arr = masker.fit_transform(deltam_file).T  # Transpose to SxT
         if deltam_arr.ndim != 2:
             raise ValueError(f'deltam is {deltam_arr.ndim}')

aslprep/tests/test_interfaces_cbf.py

@@ -27,7 +27,6 @@ def test_computecbf_casl(datasets, tmp_path_factory):
     asl_file = _save_img(asl_data, tmpdir, 'asl.nii.gz')
     asl_mask = np.zeros((30, 30, 30), dtype=np.uint8)
     asl_mask[10:20, 10:20, 10:20] = 1
-    mask_file = _save_img(asl_mask, tmpdir, 'mask.nii.gz')
     m0_file = _save_img(asl_mask, tmpdir, 'm0.nii.gz')
 
     single_pld = 1.5
@@ -67,7 +66,6 @@ def test_computecbf_casl(datasets, tmp_path_factory):
                 metadata=metadata,
                 m0_scale=1,
                 m0_file=m0_file,
-                mask=mask_file,
             )
             results = interface.run(cwd=tmpdir)
             assert os.path.isfile(results.outputs.cbf_ts)
@@ -93,7 +91,6 @@ def test_computecbf_casl(datasets, tmp_path_factory):
                 metadata=metadata,
                 m0_scale=1,
                 m0_file=m0_file,
-                mask=mask_file,
             )
             with pytest.raises(ValueError, match='Number of PostLabelingDelays'):
                 results = interface.run(cwd=tmpdir)
@@ -113,7 +110,6 @@ def test_computecbf_casl(datasets, tmp_path_factory):
                 metadata=metadata,
                 m0_scale=1,
                 m0_file=m0_file,
-                mask=mask_file,
             )
             results = interface.run(cwd=tmpdir)
             assert results.outputs.cbf_ts is None
@@ -142,7 +138,6 @@ def test_computecbf_pasl(datasets, tmp_path_factory):
     asl_file = _save_img(asl_data, tmpdir, 'asl.nii.gz')
     asl_mask = np.zeros((30, 30, 30), dtype=np.uint8)
     asl_mask[10:20, 10:20, 10:20] = 1
-    mask_file = _save_img(asl_mask, tmpdir, 'mask.nii.gz')
     m0_file = _save_img(asl_mask, tmpdir, 'm0.nii.gz')
 
     single_pld = 1.5
@@ -180,7 +175,6 @@ def test_computecbf_pasl(datasets, tmp_path_factory):
             metadata=metadata,
             m0_scale=1,
             m0_file=m0_file,
-            mask=mask_file,
         )
         with pytest.raises(ValueError, match='not supported in ASLPrep'):
             results = interface.run(cwd=tmpdir)
@@ -201,7 +195,6 @@ def test_computecbf_pasl(datasets, tmp_path_factory):
             metadata=metadata,
             m0_scale=1,
             m0_file=m0_file,
-            mask=mask_file,
         )
         results = interface.run(cwd=tmpdir)
         assert os.path.isfile(results.outputs.cbf_ts)
@@ -227,7 +220,6 @@ def test_computecbf_pasl(datasets, tmp_path_factory):
             metadata=metadata,
             m0_scale=1,
             m0_file=m0_file,
-            mask=mask_file,
         )
         with pytest.raises(ValueError, match=r'Unsupported BolusCutOffTechnique \(QUIPSS\)'):
             results = interface.run(cwd=tmpdir)
@@ -247,7 +239,6 @@ def test_computecbf_pasl(datasets, tmp_path_factory):
             metadata=metadata,
             m0_scale=1,
             m0_file=m0_file,
-            mask=mask_file,
         )
         with pytest.raises(ValueError, match=r'Unsupported BolusCutOffTechnique \(QUIPSS\)'):
             results = interface.run(cwd=tmpdir)
@@ -268,7 +259,6 @@ def test_computecbf_pasl(datasets, tmp_path_factory):
             metadata=metadata,
             m0_scale=1,
             m0_file=m0_file,
-            mask=mask_file,
         )
         results = interface.run(cwd=tmpdir)
         assert os.path.isfile(results.outputs.cbf_ts)
@@ -294,7 +284,6 @@ def test_computecbf_pasl(datasets, tmp_path_factory):
             metadata=metadata,
             m0_scale=1,
             m0_file=m0_file,
-            mask=mask_file,
         )
         results = interface.run(cwd=tmpdir)
         assert os.path.isfile(results.outputs.mean_cbf)
@@ -317,7 +306,6 @@ def test_computecbf_pasl(datasets, tmp_path_factory):
             metadata=metadata,
             m0_scale=1,
             m0_file=m0_file,
-            mask=mask_file,
         )
         results = interface.run(cwd=tmpdir)
         assert os.path.isfile(results.outputs.cbf_ts)
@@ -343,7 +331,6 @@ def test_computecbf_pasl(datasets, tmp_path_factory):
             metadata=metadata,
             m0_scale=1,
             m0_file=m0_file,
-            mask=mask_file,
         )
         results = interface.run(cwd=tmpdir)
         assert os.path.isfile(results.outputs.mean_cbf)
@@ -370,7 +357,6 @@ def test_compare_slicetiming(datasets, tmp_path_factory):
     asl_file = _save_img(asl_data, tmpdir, 'asl.nii.gz')
     asl_mask = np.zeros((30, 30, 30), dtype=np.uint8)
     asl_mask[10:20, 10:20, 10:20] = 1
-    mask_file = _save_img(asl_mask, tmpdir, 'mask.nii.gz')
     m0_file = _save_img(asl_mask, tmpdir, 'm0.nii.gz')
 
     ACQ_DICTS = [
@@ -396,7 +382,6 @@ def test_compare_slicetiming(datasets, tmp_path_factory):
             metadata=metadata,
             m0_scale=1,
             m0_file=m0_file,
-            mask=mask_file,
         )
         results = interface.run(cwd=tmpdir)
         cbf_data.append(nb.load(results.outputs.cbf_ts).get_fdata())

aslprep/workflows/asl/base.py

@@ -400,7 +400,6 @@ def init_asl_wf(
     )
     workflow.connect([
         (inputnode, cbf_wf, [
-            ('t1w_mask', 'inputnode.t1w_mask'),
             ('t1w_tpms', 'inputnode.t1w_tpms'),
             ('m0scan_metadata', 'inputnode.m0scan_metadata'),
         ]),

aslprep/workflows/asl/cbf.py

@@ -18,7 +18,6 @@
     BASILCBF,
     ComputeCBF,
     ExtractCBF,
-    RefineMask,
     ScoreAndScrubCBF,
 )
 from aslprep.interfaces.parcellation import NiftiParcellate
@@ -31,7 +30,6 @@
     pcasl_or_pasl,
     prepare_basil_kwargs,
 )
-from aslprep.workflows.asl.reference import init_synthstrip_aslref_wf
 
 
 def init_cbf_wf(
@@ -104,8 +102,6 @@ def init_cbf_wf(
         asl mask NIFTI file
     t1w_tpms
         t1w probability maps
-    t1w_mask
-        t1w mask Nifti
     aslref2anat_xfm
         asl to t1w transformation file
 
@@ -226,7 +222,6 @@ def init_cbf_wf(
                 'm0scan_metadata',
                 'asl_mask',
                 't1w_tpms',
-                't1w_mask',
                 'aslref2anat_xfm',
             ],
         ),
@@ -257,37 +252,6 @@ def init_cbf_wf(
         name='outputnode',
     )
 
-    warp_t1w_mask_to_asl = pe.Node(
-        ApplyTransforms(
-            dimension=3,
-            float=True,
-            interpolation='NearestNeighbor',
-            invert_transform_flags=[True],
-            input_image_type=3,
-            args='-v',
-        ),
-        name='warp_t1w_mask_to_asl',
-    )
-    workflow.connect([
-        (inputnode, warp_t1w_mask_to_asl, [
-            ('asl_mask', 'reference_image'),
-            ('t1w_mask', 'input_image'),
-            ('aslref2anat_xfm', 'transforms'),
-        ]),
-    ])  # fmt:skip
-
-    reduce_mask = pe.Node(
-        RefineMask(),
-        mem_gb=0.2,
-        run_without_submitting=True,
-        name='reduce_mask',
-    )
-
-    workflow.connect([
-        (inputnode, reduce_mask, [('asl_mask', 'asl_mask')]),
-        (warp_t1w_mask_to_asl, reduce_mask, [('output_image', 't1w_mask')]),
-    ])  # fmt:skip
-
     # Warp tissue probability maps to ASL space
     def _pick_gm(files):
         if not isinstance(files, list):
@@ -383,7 +347,6 @@ def _getfiledir(file):
             ('aslcontext', 'aslcontext'),
             ('m0scan_metadata', 'm0scan_metadata'),
         ]),
-        (reduce_mask, extract_deltam, [('out_mask', 'in_mask')]),
     ])  # fmt:skip
 
     if metadata['M0Type'] == 'Separate':
@@ -393,15 +356,6 @@ def _getfiledir(file):
             (mean_m0, extract_deltam, [('out_file', 'm0scan')]),
         ])  # fmt:skip
 
-        enhance_and_skullstrip_m0scan_wf = init_synthstrip_aslref_wf(
-            disable_n4=config.workflow.disable_n4,
-            name='enhance_and_skullstrip_m0scan_wf',
-        )
-        workflow.connect([
-            (mean_m0, enhance_and_skullstrip_m0scan_wf, [('out_file', 'inputnode.in_file')]),
-            (enhance_and_skullstrip_m0scan_wf, reduce_mask, [('outputnode.mask_file', 'm0_mask')]),
-        ])  # fmt:skip
-
     compute_cbf = pe.Node(
         ComputeCBF(
             cbf_only=processing_target == 'cbf',
@@ -413,7 +367,6 @@ def _getfiledir(file):
     )
 
     workflow.connect([
-        (reduce_mask, compute_cbf, [('out_mask', 'mask')]),
         (extract_deltam, compute_cbf, [
             ('out_file', 'deltam'),
             ('m0_file', 'm0_file'),
@@ -447,7 +400,7 @@ def _getfiledir(file):
 [@dolui2017structural;@dolui2016scrub].
 """
         workflow.connect([
-            (reduce_mask, score_and_scrub_cbf, [('out_mask', 'mask')]),
+            (inputnode, score_and_scrub_cbf, [('asl_mask', 'mask')]),
             (compute_cbf, score_and_scrub_cbf, [('cbf_ts', 'cbf_ts')]),
             (gm_tfm, score_and_scrub_cbf, [('output_image', 'gm_tpm')]),
             (wm_tfm, score_and_scrub_cbf, [('output_image', 'wm_tpm')]),
@@ -520,7 +473,7 @@ def _getfiledir(file):
         )
 
         workflow.connect([
-            (reduce_mask, basilcbf, [('out_mask', 'mask')]),
+            (inputnode, basilcbf, [('asl_mask', 'mask')]),
             (extract_deltam, basilcbf, [
                 ('out_file', 'deltam'),
                 ('m0_file', 'mzero'),