diff --git a/src/fmripost_aroma/cli/parser.py b/src/fmripost_aroma/cli/parser.py
index 805f965..59afc84 100644
--- a/src/fmripost_aroma/cli/parser.py
+++ b/src/fmripost_aroma/cli/parser.py
@@ -356,8 +356,9 @@ def _bids_filter(value, parser):
     g_aroma.add_argument(
         "--denoising-method",
         action="store",
-        choices=["aggr", "nonaggr", "both", "none"],
-        default="nonaggr",
+        nargs="+",
+        choices=["aggr", "nonaggr", "orthaggr"],
+        default=None,
         help="Denoising method to apply, if any.",
     )
     g_aroma.add_argument(
diff --git a/src/fmripost_aroma/config.py b/src/fmripost_aroma/config.py
index f207334..9cfc476 100644
--- a/src/fmripost_aroma/config.py
+++ b/src/fmripost_aroma/config.py
@@ -534,7 +534,7 @@ class workflow(_Config):
     melodic_dim = None
     """Number of ICA components to be estimated by MELODIC
     (positive = exact, negative = maximum)."""
-    denoising_method = "nonaggr"
+    denoising_method = None
     """Denoising strategy to be used."""
     orthogonalize = False
     """Orthogonalize the AROMA-flagged components w.r.t. the non-flagged components."""
diff --git a/src/fmripost_aroma/interfaces/confounds.py b/src/fmripost_aroma/interfaces/confounds.py
index 961373d..99645d4 100644
--- a/src/fmripost_aroma/interfaces/confounds.py
+++ b/src/fmripost_aroma/interfaces/confounds.py
@@ -131,6 +131,7 @@ def _get_ica_confounds(ica_out_dir, skip_vols, newpath=None):
 class _ICADenoiseInputSpec(BaseInterfaceInputSpec):
     method = traits.Enum("aggr", "nonaggr", "orthaggr", mandatory=True, desc="denoising method")
     bold_file = File(exists=True, mandatory=True, desc="input file to denoise")
+    mask_file = File(exists=True, mandatory=True, desc="mask file")
     confounds = File(exists=True, mandatory=True, desc="confounds file")
 
 
@@ -144,17 +145,16 @@ class ICADenoise(SimpleInterface):
     output_spec = _ICADenoiseOutputSpec
 
     def _run_interface(self, runtime):
-        import nibabel as nb
         import numpy as np
         import pandas as pd
+        from nilearn.maskers import NiftiMasker
+        from nilearn.masking import apply_mask, unmask
 
         method = self.inputs.method
         bold_file = self.inputs.bold_file
         confounds_file = self.inputs.confounds
         metrics_file = self.inputs.metrics_file
 
-        bold_img = nb.load(bold_file)
-        bold_data = bold_img.get_fdata()
         confounds_df = pd.read_table(confounds_file)
 
         # Split up component time series into accepted and rejected components
@@ -167,7 +167,7 @@ def _run_interface(self, runtime):
         if method == "aggr":
             # Denoise the data with the motion components
             masker = NiftiMasker(
-                mask_img=mask_file,
+                mask_img=self.inputs.mask_file,
                 standardize_confounds=True,
                 standardize=False,
                 smoothing_fwhm=None,
@@ -179,15 +179,10 @@ def _run_interface(self, runtime):
             )
 
             # Denoise the data by fitting and transforming the data file using the masker
-            denoised_img_2d = masker.fit_transform(data_file, confounds=rejected_components)
+            denoised_img_2d = masker.fit_transform(bold_file, confounds=rejected_components)
 
             # Transform denoised data back into 4D space
-            denoised_img_4d = masker.inverse_transform(denoised_img_2d)
-
-            # Save to file
-            denoised_img_4d.to_filename(
-                "sub-01_task-rest_space-MNI152NLin2009cAsym_desc-aggrDenoised_bold.nii.gz"
-            )
+            denoised_img = masker.inverse_transform(denoised_img_2d)
         elif method == "orthaggr":
             # Regress the good components out of the bad time series to get "pure evil" regressors
             betas = np.linalg.lstsq(accepted_components, rejected_components, rcond=None)[0]
@@ -196,7 +191,7 @@ def _run_interface(self, runtime):
 
             # Once you have these "pure evil" components, you can denoise the data
             masker = NiftiMasker(
-                mask_img=mask_file,
+                mask_img=self.inputs.mask_file,
                 standardize_confounds=True,
                 standardize=False,
                 smoothing_fwhm=None,
@@ -208,26 +203,21 @@ def _run_interface(self, runtime):
             )
 
             # Denoise the data by fitting and transforming the data file using the masker
-            denoised_img_2d = masker.fit_transform(data_file, confounds=orth_bad_timeseries)
+            denoised_img_2d = masker.fit_transform(bold_file, confounds=orth_bad_timeseries)
 
             # Transform denoised data back into 4D space
-            denoised_img_4d = masker.inverse_transform(denoised_img_2d)
-
-            # Save to file
-            denoised_img_4d.to_filename(
-                "sub-01_task-rest_space-MNI152NLin2009cAsym_desc-orthAggrDenoised_bold.nii.gz"
-            )
+            denoised_img = masker.inverse_transform(denoised_img_2d)
         else:
             # Apply the mask to the data image to get a 2d array
-            data = apply_mask(data_file, mask_file)
+            data = apply_mask(bold_file, self.inputs.mask_file)
 
-            # Fit GLM to accepted components, rejected components and nuisance regressors
+            # Fit GLM to accepted components and rejected components
             # (after adding a constant term)
             regressors = np.hstack(
                 (
                     rejected_components,
                     accepted_components,
-                    np.ones((mixing_df.shape[0], 1)),
+                    np.ones((confounds_df.shape[0], 1)),
                 ),
             )
             betas = np.linalg.lstsq(regressors, data, rcond=None)[0][:-1]
@@ -238,9 +228,9 @@ def _run_interface(self, runtime):
             data_denoised = data - pred_data
 
             # Save to file
-            denoised_img = unmask(data_denoised, mask_file)
-            denoised_img.to_filename(
-                "sub-01_task-rest_space-MNI152NLin2009cAsym_desc-nonaggrDenoised_bold.nii.gz"
-            )
+            denoised_img = unmask(data_denoised, self.inputs.mask_file)
+
+        self._results["denoised_file"] = os.path.abspath("denoised.nii.gz")
+        denoised_img.to_filename(self._results["denoised_file"])
 
         return runtime
diff --git a/src/fmripost_aroma/workflows/aroma.py b/src/fmripost_aroma/workflows/aroma.py
index ae46d21..d9e8f8e 100644
--- a/src/fmripost_aroma/workflows/aroma.py
+++ b/src/fmripost_aroma/workflows/aroma.py
@@ -313,13 +313,11 @@ def init_denoise_wf(bold_file):
     """Build a workflow that denoises a BOLD series using AROMA confounds.
 
     This workflow performs the denoising in the requested output space(s).
+    It doesn't currently work on CIFTIs.
     """
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 
-    if config.workflow.denoise_method == "all":
-        denoise_methods = ["nonaggr", "aggr", "orthaggr"]
-    else:
-        denoise_methods = [config.workflow.denoise_method]
+    from fmripost_aroma.interfaces.confounds import ICADenoise
 
     workflow = Workflow(name=_get_wf_name(bold_file, "denoise"))
 
@@ -327,6 +325,7 @@ def init_denoise_wf(bold_file):
         niu.IdentityInterface(
             fields=[
                 "bold_file",
+                "bold_mask_std",
                 "confounds",
                 "name_source",
                 "skip_vols",
@@ -347,7 +346,7 @@ def init_denoise_wf(bold_file):
         ]),
     ])  # fmt:skip
 
-    for denoise_method in denoise_methods:
+    for denoise_method in config.workflow.denoise_method:
         denoise = pe.Node(
             ICADenoise(method=denoise_method),
             name=f"denoise_{denoise_method}",
@@ -356,6 +355,7 @@ def init_denoise_wf(bold_file):
             (inputnode, denoise, [
                 ("confounds", "confounds_file"),
                 ("skip_vols", "skip_vols"),
+                ("bold_mask_std", "mask_file"),
             ]),
             (rm_non_steady_state, denoise, [("bold_cut", "bold_file")]),
         ])  # fmt:skip