diff --git a/xcp_d/interfaces/censoring.py b/xcp_d/interfaces/censoring.py
index d6d66c236..929b77a3b 100644
--- a/xcp_d/interfaces/censoring.py
+++ b/xcp_d/interfaces/censoring.py
@@ -425,6 +425,7 @@ class _GenerateConfoundsOutputSpec(TraitedSpec):
             "It will also always have the linear trend and a constant column."
         ),
     )
+    confounds_metadata = traits.Dict(desc="Metadata associated with the confounds_file output.")
     motion_file = File(
         exists=True,
         desc="The filtered motion parameters.",
@@ -497,7 +498,7 @@ def _run_interface(self, runtime):
         )
 
         # Load nuisance regressors, but use filtered motion parameters.
-        confounds_df = load_confound_matrix(
+        confounds_df, confounds_metadata = load_confound_matrix(
             params=self.inputs.params,
             img_file=self.inputs.in_file,
             confounds_file=self._results["filtered_confounds_file"],
@@ -505,6 +506,57 @@ def _run_interface(self, runtime):
             custom_confounds=self.inputs.custom_confounds_file,
         )
 
+        if confounds_df is not None:
+            # Orthogonalize full nuisance regressors w.r.t. any signal regressors
+            signal_columns = [c for c in confounds_df.columns if c.startswith("signal__")]
+            if signal_columns:
+                LOGGER.warning(
+                    "Signal columns detected. "
+                    "Orthogonalizing nuisance columns w.r.t. the following signal columns: "
+                    f"{', '.join(signal_columns)}"
+                )
+                noise_columns = [c for c in confounds_df.columns if not c.startswith("signal__")]
+
+                # Don't orthogonalize the intercept or linear trend regressors
+                untouched_cols = ["linear_trend", "intercept"]
+                cols_to_orth = [c for c in noise_columns if c not in untouched_cols]
+                orth_confounds_df = confounds_df[noise_columns].copy()
+                orth_cols = [f"{c}_orth" for c in cols_to_orth]
+                orth_confounds_df = pd.DataFrame(
+                    index=confounds_df.index,
+                    columns=orth_cols + untouched_cols,
+                )
+                orth_confounds_df.loc[:, untouched_cols] = confounds_df[untouched_cols]
+
+                # Do the orthogonalization
+                signal_regressors = confounds_df[signal_columns].to_numpy()
+                noise_regressors = confounds_df[cols_to_orth].to_numpy()
+                signal_betas = np.linalg.lstsq(signal_regressors, noise_regressors, rcond=None)[0]
+                pred_noise_regressors = np.dot(signal_regressors, signal_betas)
+                orth_noise_regressors = noise_regressors - pred_noise_regressors
+
+                # Replace the old data
+                orth_confounds_df.loc[:, orth_cols] = orth_noise_regressors
+                confounds_df = orth_confounds_df
+
+                for col in cols_to_orth:
+                    desc_str = (
+                        "This regressor is orthogonalized with respect to the 'signal' regressors "
+                        f"({', '.join(signal_columns)}) after dummy scan removal, "
+                        "but prior to any censoring."
+                    )
+
+                    col_metadata = {}
+                    if col in confounds_metadata.keys():
+                        col_metadata = confounds_metadata.pop(col)
+                        if "Description" in col_metadata.keys():
+                            desc_str = f"{col_metadata['Description']} {desc_str}"
+
+                    col_metadata["Description"] = desc_str
+                    confounds_metadata[f"{col}_orth"] = col_metadata
+
+        self._results["confounds_metadata"] = confounds_metadata
+
         # get the output
         self._results["motion_file"] = fname_presuffix(
             self.inputs.fmriprep_confounds_file,
diff --git a/xcp_d/tests/data/test_ds001419_cifti_outputs.txt b/xcp_d/tests/data/test_ds001419_cifti_outputs.txt
index e4bd18a37..f65d6c054 100644
--- a/xcp_d/tests/data/test_ds001419_cifti_outputs.txt
+++ b/xcp_d/tests/data/test_ds001419_cifti_outputs.txt
@@ -57,6 +57,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-dcan_qc.hdf5
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_design.tsv
+xcp_d/sub-01/func/sub-01_task-imagery_run-01_design.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_outliers.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_space-fsLR_atlas-4S1052Parcels_coverage.tsv
@@ -169,6 +170,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-dcan_qc.hdf5
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_design.tsv
+xcp_d/sub-01/func/sub-01_task-imagery_run-02_design.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_outliers.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_space-fsLR_atlas-4S1052Parcels_coverage.tsv
diff --git a/xcp_d/tests/data/test_ds001419_nifti_outputs.txt b/xcp_d/tests/data/test_ds001419_nifti_outputs.txt
index d320b2e11..2f78e96bc 100644
--- a/xcp_d/tests/data/test_ds001419_nifti_outputs.txt
+++ b/xcp_d/tests/data/test_ds001419_nifti_outputs.txt
@@ -55,6 +55,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-imagery_run-01_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_outliers.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-01_space-MNI152NLin2009cAsym_atlas-4S1052Parcels_coverage.tsv
@@ -151,6 +152,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_run-01_space-MNI152NLin2009cAsym_res-2_reh
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-imagery_run-02_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_outliers.json
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-imagery_run-02_space-MNI152NLin2009cAsym_atlas-4S1052Parcels_coverage.tsv
@@ -265,6 +267,7 @@ xcp_d/sub-01/func/sub-01_task-imagery_space-MNI152NLin2009cAsym_res-2_desc-denoi
 xcp_d/sub-01/func/sub-01_task-rest_desc-filtered_motion.json
 xcp_d/sub-01/func/sub-01_task-rest_desc-filtered_motion.tsv
 xcp_d/sub-01/func/sub-01_task-rest_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-rest_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-rest_outliers.json
 xcp_d/sub-01/func/sub-01_task-rest_outliers.tsv
 xcp_d/sub-01/func/sub-01_task-rest_space-MNI152NLin2009cAsym_atlas-4S1052Parcels_coverage.tsv
diff --git a/xcp_d/tests/data/test_fmriprep_without_freesurfer_outputs.txt b/xcp_d/tests/data/test_fmriprep_without_freesurfer_outputs.txt
index 096b1ca0d..a608fa821 100644
--- a/xcp_d/tests/data/test_fmriprep_without_freesurfer_outputs.txt
+++ b/xcp_d/tests/data/test_fmriprep_without_freesurfer_outputs.txt
@@ -55,6 +55,7 @@ xcp_d/sub-01/anat/sub-01_space-MNI152NLin2009cAsym_dseg.nii.gz
 xcp_d/sub-01/func
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_desc-dcan_qc.hdf5
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_motion.json
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_motion.tsv
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_outliers.json
@@ -126,6 +127,7 @@ xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_space-MNI152NLin2009cAsym_r
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-1_space-MNI152NLin2009cAsym_reho.nii.gz
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_desc-dcan_qc.hdf5
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_desc-preproc_design.tsv
+xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_desc-preproc_design.json
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_motion.json
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_motion.tsv
 xcp_d/sub-01/func/sub-01_task-mixedgamblestask_run-2_outliers.json
diff --git a/xcp_d/tests/data/test_nibabies_outputs.txt b/xcp_d/tests/data/test_nibabies_outputs.txt
index 9d7091e9e..5e1ecb352 100644
--- a/xcp_d/tests/data/test_nibabies_outputs.txt
+++ b/xcp_d/tests/data/test_nibabies_outputs.txt
@@ -53,6 +53,7 @@ xcp_d/sub-01/ses-1mo/anat/sub-01_ses-1mo_run-001_space-MNIInfant_cohort-1_dseg.n
 xcp_d/sub-01/ses-1mo/func
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-dcan_qc.hdf5
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-preproc_design.tsv
+xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_desc-preproc_design.json
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_motion.json
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_motion.tsv
 xcp_d/sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_outliers.json
diff --git a/xcp_d/tests/data/test_pnc_cifti_outputs.txt b/xcp_d/tests/data/test_pnc_cifti_outputs.txt
index 769cb840f..162f64f0e 100644
--- a/xcp_d/tests/data/test_pnc_cifti_outputs.txt
+++ b/xcp_d/tests/data/test_pnc_cifti_outputs.txt
@@ -104,6 +104,7 @@ xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleb
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-filtered_motion.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-filtered_motion.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_design.tsv
+xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_design.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_outliers.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_outliers.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_space-fsLR_atlas-4S1052Parcels_coverage.tsv
diff --git a/xcp_d/tests/data/test_pnc_nifti_outputs.txt b/xcp_d/tests/data/test_pnc_nifti_outputs.txt
index f11536ea1..1deb436ec 100644
--- a/xcp_d/tests/data/test_pnc_nifti_outputs.txt
+++ b/xcp_d/tests/data/test_pnc_nifti_outputs.txt
@@ -54,6 +54,7 @@ xcp_d/sub-1648798153/ses-PNC1/func
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_desc-filtered_motion.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_desc-filtered_motion.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_desc-preproc_design.tsv
+xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_desc-preproc_design.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_outliers.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_outliers.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_space-MNI152NLin6Asym_atlas-4S1052Parcels_coverage.tsv
@@ -150,6 +151,7 @@ xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-frac2back_space-
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-filtered_motion.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-filtered_motion.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-preproc_design.tsv
+xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_desc-preproc_design.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_outliers.json
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_outliers.tsv
 xcp_d/sub-1648798153/ses-PNC1/func/sub-1648798153_ses-PNC1_task-rest_acq-singleband_space-MNI152NLin6Asym_atlas-4S1052Parcels_coverage.tsv
diff --git a/xcp_d/tests/test_confounds.py b/xcp_d/tests/test_confounds.py
index aa5a1f631..22a71765f 100644
--- a/xcp_d/tests/test_confounds.py
+++ b/xcp_d/tests/test_confounds.py
@@ -44,7 +44,7 @@ def test_custom_confounds(ds001419_data, tmp_path_factory):
     )
     custom_confounds.to_csv(custom_confounds_file, sep="\t", index=False)
 
-    combined_confounds = load_confound_matrix(
+    combined_confounds, _ = load_confound_matrix(
         params="24P",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -56,7 +56,7 @@ def test_custom_confounds(ds001419_data, tmp_path_factory):
     assert "condition01" in combined_confounds.columns
     assert "condition02" in combined_confounds.columns
 
-    custom_confounds = load_confound_matrix(
+    custom_confounds, _ = load_confound_matrix(
         params="custom",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -126,7 +126,7 @@ def test_load_confounds(ds001419_data):
 
     N_VOLUMES = 60
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="24P",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -134,7 +134,7 @@ def test_load_confounds(ds001419_data):
     )
     assert confounds_df.shape == (N_VOLUMES, 26)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="27P",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -142,7 +142,7 @@ def test_load_confounds(ds001419_data):
     )
     assert confounds_df.shape == (N_VOLUMES, 29)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="36P",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -150,7 +150,7 @@ def test_load_confounds(ds001419_data):
     )
     assert confounds_df.shape == (N_VOLUMES, 38)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="acompcor",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -158,7 +158,7 @@ def test_load_confounds(ds001419_data):
     )
     assert confounds_df.shape == (N_VOLUMES, 30)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="acompcor_gsr",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -166,7 +166,7 @@ def test_load_confounds(ds001419_data):
     )
     assert confounds_df.shape == (N_VOLUMES, 31)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="aroma",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -174,7 +174,7 @@ def test_load_confounds(ds001419_data):
     )
     assert confounds_df.shape == (N_VOLUMES, 50)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="aroma_gsr",
         img_file=bold_file,
         confounds_file=confounds_file,
@@ -182,7 +182,7 @@ def test_load_confounds(ds001419_data):
     )
     assert confounds_df.shape == (N_VOLUMES, 51)
 
-    confounds_df = load_confound_matrix(
+    confounds_df, _ = load_confound_matrix(
         params="none",
         img_file=bold_file,
         confounds_file=confounds_file,
diff --git a/xcp_d/tests/test_interfaces_censoring.py b/xcp_d/tests/test_interfaces_censoring.py
index be8f96dc9..98cab51b7 100644
--- a/xcp_d/tests/test_interfaces_censoring.py
+++ b/xcp_d/tests/test_interfaces_censoring.py
@@ -1,4 +1,5 @@
 """Tests for framewise displacement calculation."""
+import json
 import os
 
 import nibabel as nb
@@ -16,16 +17,28 @@ def test_generate_confounds(ds001419_data, tmp_path_factory):
     confounds_json = ds001419_data["confounds_json"]
 
     df = pd.read_table(confounds_file)
+    with open(confounds_json, "r") as fo:
+        metadata = json.load(fo)
 
     # Replace confounds tsv values with values that should be omitted
     df.loc[1:3, "trans_x"] = [6, 8, 9]
     df.loc[4:6, "trans_y"] = [7, 8, 9]
     df.loc[7:9, "trans_z"] = [12, 8, 9]
 
+    # Modify JSON file
+    metadata["trans_x"] = {"test": "hello"}
+    confounds_json = os.path.join(tmpdir, "edited_confounds.json")
+    with open(confounds_json, "w") as fo:
+        json.dump(metadata, fo)
+
     # Rename with same convention as initial confounds tsv
     confounds_tsv = os.path.join(tmpdir, "edited_confounds.tsv")
     df.to_csv(confounds_tsv, sep="\t", index=False, header=True)
 
+    custom_confounds_file = os.path.join(tmpdir, "custom_confounds.tsv")
+    df2 = pd.DataFrame(columns=["signal__test"], data=np.random.random((df.shape[0], 1)))
+    df2.to_csv(custom_confounds_file, sep="\t", index=False, header=True)
+
     # Run workflow
     interface = censoring.GenerateConfounds(
         in_file=in_file,
@@ -35,7 +48,7 @@ def test_generate_confounds(ds001419_data, tmp_path_factory):
         head_radius=50,
         fmriprep_confounds_file=confounds_tsv,
         fmriprep_confounds_json=confounds_json,
-        custom_confounds_file=None,
+        custom_confounds_file=custom_confounds_file,
         motion_filter_type=None,
         motion_filter_order=4,
         band_stop_min=0,
@@ -47,6 +60,10 @@ def test_generate_confounds(ds001419_data, tmp_path_factory):
     assert os.path.isfile(results.outputs.confounds_file)
     assert os.path.isfile(results.outputs.motion_file)
     assert os.path.isfile(results.outputs.temporal_mask)
+    out_confounds_file = results.outputs.confounds_file
+    out_df = pd.read_table(out_confounds_file)
+    assert out_df.shape[1] == 26  # 24(P) + linear trend + constant
+    assert sum(out_df.columns.str.endswith("_orth")) == 24  # all 24(P), not linear trend/constant
 
 
 def test_random_censor(tmp_path_factory):
diff --git a/xcp_d/tests/test_utils_utils.py b/xcp_d/tests/test_utils_utils.py
index 3b3d50a01..9754f9f01 100644
--- a/xcp_d/tests/test_utils_utils.py
+++ b/xcp_d/tests/test_utils_utils.py
@@ -112,30 +112,6 @@ def test_denoise_with_nilearn(ds001419_data, tmp_path_factory):
     assert uncensored_denoised_bold.shape == (n_volumes, n_voxels)
     assert interpolated_filtered_bold.shape == (n_volumes, n_voxels)
 
-    # Next, do the orthogonalization
-    reduced_confounds_df["signal__test"] = confounds_df["global_signal"]
-
-    # Move intercept to end of dataframe
-    reduced_confounds_df = reduced_confounds_df[
-        [c for c in reduced_confounds_df.columns if c not in ["intercept"]] + ["intercept"]
-    ]
-    orth_confounds_file = os.path.join(tmpdir, "orth_confounds.tsv")
-    reduced_confounds_df.to_csv(orth_confounds_file, sep="\t", index=False)
-    (
-        uncensored_denoised_bold,
-        interpolated_filtered_bold,
-    ) = utils.denoise_with_nilearn(
-        preprocessed_bold=preprocessed_bold_arr,
-        confounds_file=orth_confounds_file,
-        temporal_mask=temporal_mask,
-        low_pass=low_pass,
-        high_pass=high_pass,
-        filter_order=filter_order,
-        TR=TR,
-    )
-    assert uncensored_denoised_bold.shape == (n_volumes, n_voxels)
-    assert interpolated_filtered_bold.shape == (n_volumes, n_voxels)
-
     # Finally, run without denoising
     (
         uncensored_denoised_bold,
diff --git a/xcp_d/utils/confounds.py b/xcp_d/utils/confounds.py
index fb081e7a6..8506dcc24 100644
--- a/xcp_d/utils/confounds.py
+++ b/xcp_d/utils/confounds.py
@@ -1,6 +1,7 @@
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """Confound matrix selection based on Ciric et al. 2007."""
+import json
 import os
 import warnings
 
@@ -371,6 +372,8 @@ def load_confound_matrix(
         If "AROMA" is requested, then this DataFrame will include signal components as well.
         These will be named something like "signal_[XX]".
         If ``params`` is "none", ``confounds_df`` will be None.
+    confounds_metadata : :obj:`dict`
+        Metadata for the columns in the confounds file.
     """
     PARAM_KWARGS = {
         # Get rot and trans values, as well as derivatives and square
@@ -427,7 +430,7 @@ def load_confound_matrix(
     }
 
     if params == "none":
-        return None
+        return None, {}
 
     if params in PARAM_KWARGS:
         kwargs = PARAM_KWARGS[params]
@@ -462,10 +465,24 @@ def load_confound_matrix(
         custom_confounds_df = pd.read_table(custom_confounds, sep="\t")
         confounds_df = pd.concat([custom_confounds_df, confounds_df], axis=1)
 
+    with open(confounds_json_file, "r") as fo:
+        full_confounds_metadata = json.load(fo)
+
+    confounds_metadata = {
+        k: v for k, v in full_confounds_metadata.items() if k in confounds_df.columns
+    }
+    confounds_metadata["linear_trend"] = {
+        "Description": (
+            "Linear trend regressor, generated after dummy scan removal, "
+            "but before any censoring."
+        ),
+    }
+    confounds_metadata["intercept"] = {"Description": "Intercept regressor."}
+
     confounds_df["linear_trend"] = np.arange(confounds_df.shape[0])
     confounds_df["intercept"] = np.ones(confounds_df.shape[0])
 
-    return confounds_df
+    return confounds_df, confounds_metadata
 
 
 def _get_mixing_matrix(img_file):
diff --git a/xcp_d/utils/utils.py b/xcp_d/utils/utils.py
index 14c690618..9360faeca 100644
--- a/xcp_d/utils/utils.py
+++ b/xcp_d/utils/utils.py
@@ -1,8 +1,6 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 """Miscellaneous utility functions for xcp_d."""
-import warnings
-
 import nibabel as nb
 import numpy as np
 from nipype import logging
@@ -357,14 +355,13 @@ def denoise_with_nilearn(
 
     This step does the following:
 
-        1. Orthogonalize nuisance regressors w.r.t. any signal regressors.
-        2. Censor the data and associated confounds.
-        3. Mean-center the censored and uncensored confounds, based on the censored confounds.
-        4. Estimate betas using only the censored data.
-        5. Apply the betas to denoise the *full* (uncensored) BOLD data.
-        6. Apply the betas to denoise the censored BOLD data.
-        7. Interpolate the censored, denoised data.
-        8. Bandpass filter the interpolated, denoised data.
+        1. Censor the data and associated confounds.
+        2. Mean-center the censored and uncensored confounds, based on the censored confounds.
+        3. Estimate betas using only the censored data.
+        4. Apply the betas to denoise the *full* (uncensored) BOLD data.
+        5. Apply the betas to denoise the censored BOLD data.
+        6. Interpolate the censored, denoised data.
+        7. Bandpass filter the interpolated, denoised data.
 
     Parameters
     ----------
@@ -398,7 +395,6 @@ def denoise_with_nilearn(
     # Only remove high-motion outliers in this step (not the random volumes for trimming).
     sample_mask = ~censoring_df["framewise_displacement"].to_numpy().astype(bool)
 
-    signal_columns = None
     denoise = bool(confounds_file)
     if denoise:
         confounds_df = pd.read_table(confounds_file)
@@ -407,28 +403,6 @@ def denoise_with_nilearn(
         assert "linear_trend" in confounds_df.columns
         assert confounds_df.columns[-1] == "intercept"
 
-        signal_columns = [c for c in confounds_df.columns if c.startswith("signal__")]
-
-    # Orthogonalize full nuisance regressors w.r.t. any signal regressors
-    if signal_columns:
-        warnings.warn(
-            "Signal columns detected. "
-            "Orthogonalizing nuisance columns w.r.t. the following signal columns: "
-            f"{', '.join(signal_columns)}"
-        )
-        noise_columns = [c for c in confounds_df.columns if not c.startswith("signal__")]
-        # Don't orthogonalize the intercept or linear trend regressors
-        columns_to_denoise = [c for c in noise_columns if c not in ["linear_trend", "intercept"]]
-        temp_confounds_df = confounds_df[noise_columns].copy()
-
-        signal_regressors = confounds_df[signal_columns].to_numpy()
-        noise_regressors = confounds_df[columns_to_denoise].to_numpy()
-        signal_betas = np.linalg.lstsq(signal_regressors, noise_regressors, rcond=None)[0]
-        pred_noise_regressors = np.dot(signal_regressors, signal_betas)
-        orth_noise_regressors = noise_regressors - pred_noise_regressors
-        temp_confounds_df.loc[:, columns_to_denoise] = orth_noise_regressors
-        confounds_df = temp_confounds_df
-
     # Censor the data and confounds
     preprocessed_bold_censored = preprocessed_bold[sample_mask, :]
     if denoise:
diff --git a/xcp_d/workflows/bold.py b/xcp_d/workflows/bold.py
index 2d824ccd6..0ce607bfc 100644
--- a/xcp_d/workflows/bold.py
+++ b/xcp_d/workflows/bold.py
@@ -534,6 +534,7 @@ def init_postprocess_nifti_wf(
         (inputnode, postproc_derivatives_wf, [("atlas_names", "inputnode.atlas_names")]),
         (prepare_confounds_wf, postproc_derivatives_wf, [
             ("outputnode.confounds_file", "inputnode.confounds_file"),
+            ("outputnode.confounds_metadata", "inputnode.confounds_metadata"),
             ("outputnode.filtered_motion", "inputnode.filtered_motion"),
             ("outputnode.motion_metadata", "inputnode.motion_metadata"),
             ("outputnode.temporal_mask", "inputnode.temporal_mask"),
diff --git a/xcp_d/workflows/cifti.py b/xcp_d/workflows/cifti.py
index 942908cbd..41d08ffbf 100644
--- a/xcp_d/workflows/cifti.py
+++ b/xcp_d/workflows/cifti.py
@@ -512,6 +512,7 @@ def init_postprocess_cifti_wf(
         (qc_report_wf, postproc_derivatives_wf, [("outputnode.qc_file", "inputnode.qc_file")]),
         (prepare_confounds_wf, postproc_derivatives_wf, [
             ("outputnode.confounds_file", "inputnode.confounds_file"),
+            ("outputnode.confounds_metadata", "inputnode.confounds_metadata"),
             ("outputnode.filtered_motion", "inputnode.filtered_motion"),
             ("outputnode.motion_metadata", "inputnode.motion_metadata"),
             ("outputnode.temporal_mask", "inputnode.temporal_mask"),
diff --git a/xcp_d/workflows/outputs.py b/xcp_d/workflows/outputs.py
index d307f3b65..553897f63 100644
--- a/xcp_d/workflows/outputs.py
+++ b/xcp_d/workflows/outputs.py
@@ -208,6 +208,7 @@ def init_postproc_derivatives_wf(
     reho
     parcellated_reho
     confounds_file
+    confounds_metadata
     %(filtered_motion)s
     motion_metadata
     %(temporal_mask)s
@@ -221,6 +222,7 @@ def init_postproc_derivatives_wf(
             fields=[
                 "atlas_names",
                 "confounds_file",
+                "confounds_metadata",
                 "coverage",
                 "timeseries",
                 "correlations",
@@ -329,7 +331,14 @@ def init_postproc_derivatives_wf(
             name="ds_confounds",
             run_without_submitting=False,
         )
-        workflow.connect([(inputnode, ds_confounds, [("confounds_file", "in_file")])])
+        # fmt:off
+        workflow.connect([
+            (inputnode, ds_confounds, [
+                ("confounds_file", "in_file"),
+                ("confounds_metadata", "meta_dict"),
+            ]),
+        ])
+        # fmt:on
 
     ds_coverage_files = pe.MapNode(
         DerivativesDataSink(
diff --git a/xcp_d/workflows/postprocessing.py b/xcp_d/workflows/postprocessing.py
index 240484881..6798b471b 100644
--- a/xcp_d/workflows/postprocessing.py
+++ b/xcp_d/workflows/postprocessing.py
@@ -112,6 +112,7 @@ def init_prepare_confounds_wf(
     %(fmriprep_confounds_file)s
     confounds_file : :obj:`str`
         The selected confounds, potentially including custom confounds, after dummy scan removal.
+    confounds_metadata : :obj:`dict`
     %(dummy_scans)s
         If originally set to "auto", this output will have the actual number of dummy volumes.
     %(filtered_motion)s
@@ -181,6 +182,7 @@ def init_prepare_confounds_wf(
                 "preprocessed_bold",
                 "fmriprep_confounds_file",  # used to calculate motion in concatenation workflow
                 "confounds_file",
+                "confounds_metadata",
                 "dummy_scans",
                 "filtered_motion",
                 "motion_metadata",
@@ -216,7 +218,10 @@ def init_prepare_confounds_wf(
             ("fmriprep_confounds_json", "fmriprep_confounds_json"),
             ("custom_confounds_file", "custom_confounds_file"),
         ]),
-        (generate_confounds, outputnode, [("motion_metadata", "motion_metadata")]),
+        (generate_confounds, outputnode, [
+            ("motion_metadata", "motion_metadata"),
+            ("confounds_metadata", "confounds_metadata"),
+        ]),
     ])
     # fmt:on