Repository: PennLINC/qsiprep (language: python)

Issue:
[ENH] Warp all reconstruction scalar image outputs to template space
 - [ ] if a T1w -> MNI warp exists, automatically transform scalar outputs
 - [ ] Add documentation once it works


Discussion / hints:

Fix patch:
diff --git a/.circleci/ScalarMapper.sh b/.circleci/ScalarMapper.sh
index 0ab6ecdc..8659d9ec 100644
--- a/.circleci/ScalarMapper.sh
+++ b/.circleci/ScalarMapper.sh
@@ -34,6 +34,7 @@ ${QSIPREP_CMD} \
 	--recon-input ${BIDS_INPUT_DIR} \
 	--sloppy \
     --stop-on-first-crash \
-	--recon-spec bundle_scalar_map \
+	--recon-spec test_scalar_maps \
 	--recon-only \
+	--output-resolution 3.5 \
 	-vv
\ No newline at end of file
diff --git a/.circleci/get_data.sh b/.circleci/get_data.sh
index 985d7ecb..5bb393b3 100644
--- a/.circleci/get_data.sh
+++ b/.circleci/get_data.sh
@@ -42,7 +42,7 @@ run_qsiprep_cmd () {
 
     if [[ -n "${LOCAL_PATCH}" ]]; then
       #echo "Using qsiprep patch: ${LOCAL_PATCH}"
-      QSIPREP_RUN="${QSIPREP_RUN} --patch-qsiprep ${LOCAL_PATCH}"
+      QSIPREP_RUN="${QSIPREP_RUN} --patch-qsiprep ${LOCAL_PATCH} --patch-nipype /Users/mcieslak/projects/Nipype/nipype"
     fi
   fi
   echo "${QSIPREP_RUN} --nthreads ${NTHREADS} --omp-nthreads ${OMP_NTHREADS}"
diff --git a/qsiprep/cli/run.py b/qsiprep/cli/run.py
index f731b14d..a43b13a8 100644
--- a/qsiprep/cli/run.py
+++ b/qsiprep/cli/run.py
@@ -1221,7 +1221,9 @@ def build_recon_workflow(opts, retval):
         b0_threshold=opts.b0_threshold,
         freesurfer_input=opts.freesurfer_input,
         skip_odf_plots=opts.skip_odf_reports,
-        pipeline_source=opts.recon_input_pipeline
+        pipeline_source=opts.recon_input_pipeline,
+        output_resolution=opts.output_resolution,
+        infant_mode=opts.infant
     )
     retval['return_code'] = 0
 
diff --git a/qsiprep/data/pipelines/TORTOISE.json b/qsiprep/data/pipelines/TORTOISE.json
index fe048b6b..0683dc0b 100644
--- a/qsiprep/data/pipelines/TORTOISE.json
+++ b/qsiprep/data/pipelines/TORTOISE.json
@@ -8,7 +8,7 @@
       "name": "tortoise_dtmapmri",
       "software": "TORTOISE",
       "action": "estimate",
-      "output_suffix": "tortoise",
+      "qsirecon_suffix": "TORTOISE",
       "input": "qsiprep",
       "parameters": {
         "estimate_tensor_separately": true,
diff --git a/qsiprep/data/pipelines/amico_noddi.json b/qsiprep/data/pipelines/amico_noddi.json
index 07278800..430c813c 100644
--- a/qsiprep/data/pipelines/amico_noddi.json
+++ b/qsiprep/data/pipelines/amico_noddi.json
@@ -8,7 +8,7 @@
       "action": "fit_noddi",
       "software": "AMICO",
       "input": "qsiprep",
-      "output_suffix": "NODDI",
+      "qsirecon_suffix": "NODDI",
       "parameters": {
         "isExvivo": false,
         "dPar": 1.7E-3,
diff --git a/qsiprep/data/pipelines/csdsi_3dshore.json b/qsiprep/data/pipelines/csdsi_3dshore.json
index aaf246bb..c66c294c 100644
--- a/qsiprep/data/pipelines/csdsi_3dshore.json
+++ b/qsiprep/data/pipelines/csdsi_3dshore.json
@@ -8,7 +8,7 @@
       "software": "Dipy",
       "action": "3dSHORE_reconstruction",
       "input": "qsiprep",
-      "output_suffix": "3dSHORE",
+      "qsirecon_suffix": "3dSHORE",
       "parameters": {
         "extrapolate_scheme": "HCP",
         "radial_order": 8,
@@ -43,13 +43,13 @@
         "name": "controlability",
         "input": "streamline_connectivity",
         "action": "controllability",
-        "output_suffix": "rec-3dshore_control"
+        "qsirecon_suffix": "rec-3dshore_control"
     },
     {
       "name": "msmt_csd",
       "software": "MRTrix3",
       "action": "csd",
-      "output_suffix": "msmtcsd",
+      "qsirecon_suffix": "msmtcsd",
       "input":"csdsi_3dshore",
       "parameters": {
         "mtnormalize": true,
@@ -66,7 +66,7 @@
       "name": "track_ifod2",
       "software": "MRTrix3",
       "action": "tractography",
-      "output_suffix": "ifod2",
+      "qsirecon_suffix": "ifod2",
       "input": "msmt_csd",
       "parameters": {
         "use_5tt": true,
@@ -87,7 +87,7 @@
       "name": "mrtrix_conn",
       "software": "MRTrix3",
       "action": "connectivity",
-      "output_suffix": "msmtconnectome",
+      "qsirecon_suffix": "msmtconnectome",
       "input": "track_ifod2",
       "parameters":{
         "tck2connectome": [
diff --git a/qsiprep/data/pipelines/dipy_3dshore.json b/qsiprep/data/pipelines/dipy_3dshore.json
index 91e9ec90..b8e74850 100644
--- a/qsiprep/data/pipelines/dipy_3dshore.json
+++ b/qsiprep/data/pipelines/dipy_3dshore.json
@@ -8,7 +8,7 @@
       "software": "Dipy",
       "action": "3dSHORE_reconstruction",
       "input": "qsiprep",
-      "output_suffix": "3dSHORE",
+      "qsirecon_suffix": "3dSHORE",
       "parameters": {
         "radial_order": 6,
         "zeta": 700,
@@ -42,7 +42,7 @@
       "software": "DSI Studio",
       "action": "connectivity",
       "input": "tractography",
-      "output_suffix": "3dshorenetwork",
+      "qsirecon_suffix": "3dshorenetwork",
       "parameters": {
         "connectivity_value": "count,ncount,mean_length,gfa",
         "connectivity_type": "pass,end"
diff --git a/qsiprep/data/pipelines/dipy_dki.json b/qsiprep/data/pipelines/dipy_dki.json
index f986c0f2..677dfbfd 100644
--- a/qsiprep/data/pipelines/dipy_dki.json
+++ b/qsiprep/data/pipelines/dipy_dki.json
@@ -8,7 +8,7 @@
       "software": "Dipy",
       "action": "DKI_reconstruction",
       "input": "qsiprep",
-      "output_suffix": "DKI",
+      "qsirecon_suffix": "DKI",
       "parameters": {
         "write_mif": false,
         "write_fibgz": false
diff --git a/qsiprep/data/pipelines/dipy_mapmri.json b/qsiprep/data/pipelines/dipy_mapmri.json
index 79291324..5167946c 100644
--- a/qsiprep/data/pipelines/dipy_mapmri.json
+++ b/qsiprep/data/pipelines/dipy_mapmri.json
@@ -8,7 +8,7 @@
       "software": "Dipy",
       "action": "MAPMRI_reconstruction",
       "input": "qsiprep",
-      "output_suffix": "MAPMRI",
+      "qsirecon_suffix": "DIPYMAPMRI",
       "parameters": {
         "radial_order": 6,
         "laplacian_regularization": true,
diff --git a/qsiprep/data/pipelines/dsi_studio_autotrack.json b/qsiprep/data/pipelines/dsi_studio_autotrack.json
index 4e62e7b2..62d22b08 100644
--- a/qsiprep/data/pipelines/dsi_studio_autotrack.json
+++ b/qsiprep/data/pipelines/dsi_studio_autotrack.json
@@ -8,7 +8,7 @@
       "software": "DSI Studio",
       "action": "reconstruction",
       "input": "qsiprep",
-      "output_suffix": "gqi",
+      "qsirecon_suffix": "DSIStudio",
       "parameters": {"method": "gqi"}
     },
     {
@@ -16,7 +16,7 @@
       "software": "DSI Studio",
       "action": "autotrack",
       "input": "dsistudio_gqi",
-      "output_suffix": "AutoTrackGQI",
+      "qsirecon_suffix": "DSIStudio",
       "parameters": {
         "track_id": "Fasciculus,Cingulum,Aslant,Corticos,Thalamic_R,Reticular,Optic,Fornix,Corpus",
         "tolerance": "22,26,30",
diff --git a/qsiprep/data/pipelines/dsi_studio_gqi.json b/qsiprep/data/pipelines/dsi_studio_gqi.json
index 4dcb4d3d..d96ceaed 100644
--- a/qsiprep/data/pipelines/dsi_studio_gqi.json
+++ b/qsiprep/data/pipelines/dsi_studio_gqi.json
@@ -8,7 +8,7 @@
       "software": "DSI Studio",
       "action": "reconstruction",
       "input": "qsiprep",
-      "output_suffix": "gqi",
+      "qsirecon_suffix": "DSIStudio",
       "parameters": {"method": "gqi"}
     },
     {
@@ -16,7 +16,7 @@
       "software": "DSI Studio",
       "action": "export",
       "input": "dsistudio_gqi",
-      "output_suffix": "gqiscalar"
+      "qsirecon_suffix": "DSIStudio"
     },
     {
       "name": "tractography",
@@ -37,7 +37,7 @@
       "software": "DSI Studio",
       "action": "connectivity",
       "input": "tractography",
-      "output_suffix": "gqinetwork",
+      "qsirecon_suffix": "DSIStudio",
       "parameters": {
         "connectivity_value": "count,ncount,mean_length,gfa",
         "connectivity_type": "pass,end"
diff --git a/qsiprep/data/pipelines/mrtrix_multishell_msmt_ACT-fast.json b/qsiprep/data/pipelines/mrtrix_multishell_msmt_ACT-fast.json
index de1524c8..e9613dc1 100644
--- a/qsiprep/data/pipelines/mrtrix_multishell_msmt_ACT-fast.json
+++ b/qsiprep/data/pipelines/mrtrix_multishell_msmt_ACT-fast.json
@@ -10,7 +10,7 @@
       "name": "msmt_csd",
       "software": "MRTrix3",
       "action": "csd",
-      "output_suffix": "msmtcsd",
+      "qsirecon_suffix": "MRtrix3_act-FAST",
       "input": "qsiprep",
       "parameters": {
         "mtnormalize": true,
@@ -31,7 +31,7 @@
       "name": "track_ifod2",
       "software": "MRTrix3",
       "action": "tractography",
-      "output_suffix": "ifod2",
+      "qsirecon_suffix": "MRtrix3_act-FAST",
       "input": "msmt_csd",
       "parameters": {
         "use_5tt": true,
@@ -54,7 +54,7 @@
       "name": "mrtrix_conn",
       "software": "MRTrix3",
       "action": "connectivity",
-      "output_suffix": "msmtconnectome",
+      "qsirecon_suffix": "MRtrix3_act-FAST",
       "input": "track_ifod2",
       "parameters": {
         "tck2connectome": [
diff --git a/qsiprep/data/pipelines/mrtrix_multishell_msmt_ACT-hsvs.json b/qsiprep/data/pipelines/mrtrix_multishell_msmt_ACT-hsvs.json
index f46fde49..decd310f 100644
--- a/qsiprep/data/pipelines/mrtrix_multishell_msmt_ACT-hsvs.json
+++ b/qsiprep/data/pipelines/mrtrix_multishell_msmt_ACT-hsvs.json
@@ -10,7 +10,7 @@
       "name": "msmt_csd",
       "software": "MRTrix3",
       "action": "csd",
-      "output_suffix": "msmtcsd",
+      "qsirecon_suffix": "MRtrix3_act-HSVS",
       "input": "qsiprep",
       "parameters": {
         "mtnormalize": true,
@@ -31,7 +31,7 @@
       "name": "track_ifod2",
       "software": "MRTrix3",
       "action": "tractography",
-      "output_suffix": "ifod2",
+      "qsirecon_suffix": "MRtrix3_act-HSVS",
       "input": "msmt_csd",
       "parameters": {
         "use_5tt": true,
@@ -54,7 +54,7 @@
       "name": "mrtrix_conn",
       "software": "MRTrix3",
       "action": "connectivity",
-      "output_suffix": "msmtconnectome",
+      "qsirecon_suffix": "MRtrix3_act-HSVS",
       "input": "track_ifod2",
       "parameters": {
         "tck2connectome": [
diff --git a/qsiprep/data/pipelines/mrtrix_multishell_msmt_noACT.json b/qsiprep/data/pipelines/mrtrix_multishell_msmt_noACT.json
index 9f0c17dc..da0abb4c 100644
--- a/qsiprep/data/pipelines/mrtrix_multishell_msmt_noACT.json
+++ b/qsiprep/data/pipelines/mrtrix_multishell_msmt_noACT.json
@@ -8,7 +8,7 @@
       "name": "msmt_csd",
       "software": "MRTrix3",
       "action": "csd",
-      "output_suffix": "msmtcsd",
+      "qsirecon_suffix": "MRtrix3_act-None",
       "input": "qsiprep",
       "parameters": {
         "mtnormalize": true,
@@ -25,7 +25,7 @@
       "name": "track_ifod2",
       "software": "MRTrix3",
       "action": "tractography",
-      "output_suffix": "ifod2",
+      "qsirecon_suffix": "MRtrix3_act-None",
       "input": "msmt_csd",
       "parameters": {
         "use_5tt": false,
@@ -45,7 +45,7 @@
       "name": "mrtrix_conn",
       "software": "MRTrix3",
       "action": "connectivity",
-      "output_suffix": "dhollanderconnectome",
+      "qsirecon_suffix": "MRtrix3_act-None",
       "input": "track_ifod2",
       "parameters": {
         "tck2connectome": [
diff --git a/qsiprep/data/pipelines/mrtrix_multishell_msmt_pyafq_tractometry.json b/qsiprep/data/pipelines/mrtrix_multishell_msmt_pyafq_tractometry.json
index 351baeb1..2582b507 100644
--- a/qsiprep/data/pipelines/mrtrix_multishell_msmt_pyafq_tractometry.json
+++ b/qsiprep/data/pipelines/mrtrix_multishell_msmt_pyafq_tractometry.json
@@ -8,7 +8,7 @@
             "name": "msmt_csd",
             "software": "MRTrix3",
             "action": "csd",
-            "output_suffix": "msmtcsd",
+            "qsirecon_suffix": "MRtrix3",
             "input": "qsiprep",
             "parameters": {
                 "mtnormalize": true,
@@ -25,7 +25,7 @@
             "name": "track_ifod2",
             "software": "MRTrix3",
             "action": "tractography",
-            "output_suffix": "ifod2",
+            "qsirecon_suffix": "MRtrix3",
             "input": "msmt_csd",
             "parameters": {
                 "use_5tt": false,
@@ -45,7 +45,7 @@
             "software": "pyAFQ",
             "action": "pyafq_tractometry",
             "input": "track_ifod2",
-            "output_suffix": "PYAFQ_TRACTOMETRY_ET",
+            "qsirecon_suffix": "PYAFQ",
             "parameters": {
                 "use_external_tracking": true,
                 "export": "all",
diff --git a/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.json b/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.json
index d7786eca..ec49a12b 100644
--- a/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.json
+++ b/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_ACT-fast.json
@@ -8,7 +8,7 @@
       "name": "ss3t_csd",
       "software": "MRTrix3",
       "action": "csd",
-      "output_suffix": "ss3tcsd",
+      "qsirecon_suffix": "MRtrix3_fork-SS3T_act-FAST",
       "input": "qsiprep",
       "parameters": {
         "mtnormalize": true,
@@ -24,7 +24,7 @@
       "name": "track_ifod2",
       "software": "MRTrix3",
       "action": "tractography",
-      "output_suffix": "ifod2",
+      "qsirecon_suffix": "MRtrix3_fork-SS3T_act-FAST",
       "input": "ss3t_csd",
       "parameters": {
         "use_5tt": true,
@@ -47,7 +47,7 @@
       "name": "mrtrix_conn",
       "software": "MRTrix3",
       "action": "connectivity",
-      "output_suffix": "dhollanderconnectome",
+      "qsirecon_suffix": "MRtrix3_fork-SS3T_act-FAST",
       "input": "track_ifod2",
       "parameters":{
         "tck2connectome": [
diff --git a/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_ACT-hsvs.json b/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_ACT-hsvs.json
index c53a77b8..e195728a 100644
--- a/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_ACT-hsvs.json
+++ b/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_ACT-hsvs.json
@@ -8,7 +8,7 @@
       "name": "ss3t_csd",
       "software": "MRTrix3",
       "action": "csd",
-      "output_suffix": "ss3tcsd",
+      "qsirecon_suffix": "MRtrix3_fork-SS3T_act-HSVS",
       "input": "qsiprep",
       "parameters": {
         "mtnormalize": true,
@@ -24,7 +24,7 @@
       "name": "track_ifod2",
       "software": "MRTrix3",
       "action": "tractography",
-      "output_suffix": "ifod2",
+      "qsirecon_suffix": "MRtrix3_fork-SS3T_act-HSVS",
       "input": "ss3t_csd",
       "parameters": {
         "use_5tt": true,
@@ -47,7 +47,7 @@
       "name": "mrtrix_conn",
       "software": "MRTrix3",
       "action": "connectivity",
-      "output_suffix": "dhollanderconnectome",
+      "qsirecon_suffix": "MRtrix3_fork-SS3T_act-HSVS",
       "input": "track_ifod2",
       "parameters":{
         "tck2connectome": [
diff --git a/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_noACT.json b/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_noACT.json
index 3e60ce41..f9c0214b 100644
--- a/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_noACT.json
+++ b/qsiprep/data/pipelines/mrtrix_singleshell_ss3t_noACT.json
@@ -8,7 +8,7 @@
       "name": "ss3t_csd",
       "software": "MRTrix3",
       "action": "csd",
-      "output_suffix": "ss3tcsd",
+      "qsirecon_suffix": "MRtrix3_fork-SS3T_act-None",
       "input": "qsiprep",
       "parameters": {
         "mtnormalize": true,
@@ -24,7 +24,7 @@
       "name": "track_ifod2",
       "software": "MRTrix3",
       "action": "tractography",
-      "output_suffix": "ifod2",
+      "qsirecon_suffix": "MRtrix3_fork-SS3T_act-None",
       "input": "ss3t_csd",
       "parameters": {
         "use_5tt": false,
@@ -44,7 +44,7 @@
       "name": "mrtrix_conn",
       "software": "MRTrix3",
       "action": "connectivity",
-      "output_suffix": "dhollanderconnectome",
+      "qsirecon_suffix": "MRtrix3_fork-SS3T_act-None",
       "input": "track_ifod2",
       "parameters":{
         "tck2connectome": [
diff --git a/qsiprep/data/pipelines/mrtrix_tckglobal.json b/qsiprep/data/pipelines/mrtrix_tckglobal.json
deleted file mode 100644
index 7851f795..00000000
--- a/qsiprep/data/pipelines/mrtrix_tckglobal.json
+++ /dev/null
@@ -1,66 +0,0 @@
-{
-  "name": "mrtrix_msmt",
-  "space": "T1w",
-  "atlases": ["schaefer100", "schaefer200", "schaefer400", "brainnetome246", "aicha384", "gordon333", "aal116"],
-  "anatomical": [],
-  "nodes": [
-    {
-      "name": "msmt_csd",
-      "software": "MRTrix3",
-      "action": "csd",
-      "output_suffix": "msmtcsd",
-      "input": "qsiprep",
-      "parameters": {
-        "mtnormalize": true,
-        "response": {
-          "algorithm": "dhollander"
-        },
-        "fod": {
-          "algorithm": "msmt_csd",
-          "max_sh": [8, 8, 8]
-        }
-      }
-    },
-    {
-      "name": "global_tractography",
-      "software": "MRTrix3",
-      "action": "global_tractography",
-      "output_suffix": "tckglobal",
-      "input": "msmt_csd",
-      "parameters": {
-        "niter": 10000000
-      }
-    },
-    {
-      "name": "mrtrix_conn",
-      "software": "MRTrix3",
-      "action": "connectivity",
-      "output_suffix": "globalconnectome",
-      "input": "global_tractography",
-      "parameters":{
-        "use_sift_weights": false,
-        "tck2connectome": {
-          "zero_diagonal":false,
-          "search_radius": 2.0,
-          "scale_invnodevol":true,
-          "symmetric": true
-        }
-      }
-    },
-    {
-      "name": "mrtrix_invlen_conn",
-      "software": "MRTrix3",
-      "action": "connectivity",
-      "output_suffix": "globalinvlenconnectome",
-      "input": "global_tractography",
-      "parameters":{
-        "use_sift_weights": false,
-        "tck2connectome": {
-          "zero_diagonal":false,
-          "search_radius": 2.0,
-          "symmetric": true
-        }
-      }
-    }
-  ]
-}
diff --git a/qsiprep/data/pipelines/multishell_scalarfest.json b/qsiprep/data/pipelines/multishell_scalarfest.json
index 8f851eb0..eb802476 100644
--- a/qsiprep/data/pipelines/multishell_scalarfest.json
+++ b/qsiprep/data/pipelines/multishell_scalarfest.json
@@ -8,7 +8,7 @@
       "name": "msmt_csd",
       "software": "MRTrix3",
       "action": "csd",
-      "output_suffix": "msmtcsd",
+      "qsirecon_suffix": "msmtcsd",
       "input": "qsiprep",
       "parameters": {
         "mtnormalize": true,
@@ -26,7 +26,7 @@
       "action": "fit_noddi",
       "software": "AMICO",
       "input": "qsiprep",
-      "output_suffix": "wmNODDI",
+      "qsirecon_suffix": "wmNODDI",
       "parameters": {
         "isExvivo": false,
         "dPar": 1.7E-3,
@@ -38,7 +38,7 @@
       "action": "fit_noddi",
       "software": "AMICO",
       "input": "qsiprep",
-      "output_suffix": "gmNODDI",
+      "qsirecon_suffix": "gmNODDI",
       "parameters": {
         "isExvivo": false,
         "dPar": 1.1E-3,
@@ -49,7 +49,7 @@
       "name": "calc_sops",
       "action": "steinhardt_order_parameters",
       "input": "msmt_csd",
-      "output_suffix": "SOP",
+      "qsirecon_suffix": "SOP",
       "parameters": {
         "order": 8
       }
@@ -59,7 +59,7 @@
       "software": "Dipy",
       "action": "MAPMRI_reconstruction",
       "input": "qsiprep",
-      "output_suffix": "MAPMRI",
+      "qsirecon_suffix": "MAPMRI",
       "parameters": {
         "radial_order": 6,
         "laplacian_regularization": true,
@@ -76,7 +76,7 @@
       "software": "Dipy",
       "action": "DKI_reconstruction",
       "input": "qsiprep",
-      "output_suffix": "DKI",
+      "qsirecon_suffix": "DKI",
       "parameters": {
         "write_mif": false,
         "write_fibgz": false
@@ -87,7 +87,7 @@
       "software": "Dipy",
       "action": "3dSHORE_reconstruction",
       "input": "qsiprep",
-      "output_suffix": "3dSHORE",
+      "qsirecon_suffix": "DIPY3dSHORE",
       "parameters": {
         "extrapolate_scheme": "HCP",
         "radial_order": 8,
@@ -102,7 +102,7 @@
       "name": "msmt_csd_fromcs",
       "software": "MRTrix3",
       "action": "csd",
-      "output_suffix": "csmsmtcsd",
+      "qsirecon_suffix": "csmsmtcsd",
       "input":"csdsi_3dshore",
       "parameters": {
         "mtnormalize": true,
@@ -120,7 +120,7 @@
       "software": "DSI Studio",
       "action": "reconstruction",
       "input": "qsiprep",
-      "output_suffix": "gqi",
+      "qsirecon_suffix": "DSIStudio",
       "parameters": {"method": "gqi"}
     },
     {
@@ -128,7 +128,7 @@
       "software": "DSI Studio",
       "action": "export",
       "input": "dsistudio_gqi",
-      "output_suffix": "gqiscalar"
+      "qsirecon_suffix": "DSIStudio"
     }
   ]
 }
diff --git a/qsiprep/data/pipelines/pyafq_tractometry.json b/qsiprep/data/pipelines/pyafq_tractometry.json
index cf50e428..050ec401 100644
--- a/qsiprep/data/pipelines/pyafq_tractometry.json
+++ b/qsiprep/data/pipelines/pyafq_tractometry.json
@@ -9,7 +9,7 @@
             "software": "pyAFQ",
             "action": "pyafq_tractometry",
             "input": "qsiprep",
-            "output_suffix": "PYAFQ_TRACTOMETRY",
+            "qsirecon_suffix": "PYAFQ",
             "parameters": {
                 "use_external_tracking": false,
                 "export": "all",
diff --git a/qsiprep/data/pipelines/reorient_fslstd.json b/qsiprep/data/pipelines/reorient_fslstd.json
index 0d17f78c..f20836d2 100644
--- a/qsiprep/data/pipelines/reorient_fslstd.json
+++ b/qsiprep/data/pipelines/reorient_fslstd.json
@@ -7,7 +7,7 @@
       "name": "to_fslstd",
       "action": "reorient_fslstd",
       "input": "qsiprep",
-      "output_suffix": "fslstd"
+      "qsirecon_suffix": "FSL"
     }
   ]
 }
diff --git a/qsiprep/interfaces/anatomical.py b/qsiprep/interfaces/anatomical.py
index 0d24fb73..790af396 100644
--- a/qsiprep/interfaces/anatomical.py
+++ b/qsiprep/interfaces/anatomical.py
@@ -18,13 +18,14 @@
 import numpy as np
 from dipy.segment.threshold import otsu
 from nipype import logging
-from nipype.interfaces.base import (
-    BaseInterfaceInputSpec,
-    File,
-    SimpleInterface,
-    TraitedSpec,
-    traits,
-)
+from glob import glob
+import nibabel as nb
+from scipy.spatial import distance
+from scipy import ndimage
+from nipype.interfaces.base import (traits, TraitedSpec, BaseInterfaceInputSpec,
+                                    SimpleInterface, File, isdefined)
+from pkg_resources import resource_filename as pkgrf
+
 from nipype.utils.filemanip import fname_presuffix
 from pkg_resources import resource_filename as pkgr
 from scipy import ndimage
@@ -43,6 +44,7 @@ class QsiprepAnatomicalIngressInputSpec(BaseInterfaceInputSpec):
                                        help="directory containing subject results directories")
     subject_id = traits.Str()
     subjects_dir = File(exists=True)
+    infant_mode = traits.Bool(mandatory=True)
 
 
 class QsiprepAnatomicalIngressOutputSpec(TraitedSpec):
@@ -68,6 +70,8 @@ class QsiprepAnatomicalIngressOutputSpec(TraitedSpec):
     t1_2_mni_reverse_transform = File()
     # sub-1_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5
     t1_2_mni_forward_transform = File()
+    # Generic: what template was used?
+    template_image = File()
 
 
 class QsiprepAnatomicalIngress(SimpleInterface):
@@ -128,6 +132,13 @@ def _run_interface(self, runtime):
         self._get_if_exists(
             't1_2_mni_forward_transform',
             "%s/sub-%s*_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5" % (anat_root, sub))
+        if not self.inputs.infant_mode:
+            self._results["template_image"] = pkgrf(
+                'qsiprep', 'data/mni_1mm_t1w_lps_brain.nii.gz')
+        else:
+            self._results["template_image"] = pkgrf(
+                'qsiprep', 'data/mni_1mm_t1w_lps_brain_infant.nii.gz')
+
         return runtime
 
     def _get_if_exists(self, name, pattern, excludes=None):
@@ -195,6 +206,42 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _VoxelSizeChooserInputSpec(BaseInterfaceInputSpec):
+    voxel_size = traits.Float()
+    input_image = File(exists=True)
+    anisotropic_strategy = traits.Enum("min", "max", "mean", usedefault=True)
+
+
+class _VoxelSizeChooserOutputSpec(TraitedSpec):
+    voxel_size = traits.Float()
+
+
+class VoxelSizeChooser(SimpleInterface):
+    input_spec = _VoxelSizeChooserInputSpec
+    output_spec = _VoxelSizeChooserOutputSpec
+
+    def _run_interface(self, runtime):
+        if not isdefined(self.inputs.input_image) and not isdefined(self.inputs.voxel_size):
+            raise Exception("Either voxel_size or input_image need to be defined")
+
+        # A voxel size was specified without an image
+        if isdefined(self.inputs.voxel_size):
+            voxel_size = self.inputs.voxel_size
+        else:
+            # An image was provided
+            img = nb.load(self.inputs.input_image)
+            zooms = img.header.get_zooms()[:3]
+            if self.inputs.anisotropic_strategy == "min":
+                voxel_size = min(zooms)
+            elif self.inputs.anisotropic_strategy == "max":
+                voxel_size = max(zooms)
+            else:
+                voxel_size = np.round(np.mean(zooms), 2)
+
+        self._results["voxel_size"] = voxel_size
+        return runtime
+
+
 class _FakeSegmentationInputSpec(BaseInterfaceInputSpec):
     mask_file = File(exists=True, mandatory=True)
 
diff --git a/qsiprep/interfaces/bids.py b/qsiprep/interfaces/bids.py
index 4a24de98..ccbf05d8 100644
--- a/qsiprep/interfaces/bids.py
+++ b/qsiprep/interfaces/bids.py
@@ -23,9 +23,12 @@
 import os
 import os.path as op
 import re
-from shutil import copyfileobj, copytree, rmtree
 
 import simplejson as json
+import gzip
+from shutil import copytree, rmtree, copyfileobj
+from bids.layout import parse_file_entities
+
 from nipype import logging
 from nipype.interfaces.base import (
     BaseInterfaceInputSpec,
@@ -39,7 +42,7 @@
     isdefined,
     traits,
 )
-from nipype.utils.filemanip import copyfile
+from nipype.utils.filemanip import copyfile, split_filename
 
 LOGGER = logging.getLogger('nipype.interface')
 BIDS_NAME = re.compile(
@@ -253,8 +256,7 @@ class DerivativesDataSink(SimpleInterface):
     >>> dsink.inputs.desc = 'preproc'
     >>> res = dsink.run()
     >>> res.outputs.out_file  # doctest: +ELLIPSIS
-    '.../qsiprep/sub-02/ses-noanat/func/sub-02_ses-noanat_task-rest_run-01_\
-desc-preproc_bold.nii.gz'
+    '.../qsiprep/sub-02/ses-noanat/func/sub-02_ses-noanat_task-rest_run-01_desc-preproc_bold.nii.gz'
 
     """
     input_spec = DerivativesDataSinkInputSpec
@@ -295,7 +297,7 @@ def _run_interface(self, runtime):
         os.makedirs(out_path, exist_ok=True)
         base_fname = op.join(out_path, src_fname)
 
-        formatstr = '{bname}{space}{desc}{bundle}{suffix}{dtype}{ext}'
+        formatstr = '{bname}{space}{desc}{suffix}{dtype}{ext}'
 
         space = '_space-{}'.format(self.inputs.space) if self.inputs.space else ''
         desc = '_desc-{}'.format(self.inputs.desc) if self.inputs.desc else ''
@@ -311,7 +313,6 @@ def _run_interface(self, runtime):
                 bname=base_fname,
                 space=space,
                 desc=desc,
-                bundle=bundle,
                 suffix=suffix,
                 dtype=dtype,
                 ext='')
@@ -321,7 +322,7 @@ def _run_interface(self, runtime):
             return runtime
 
         if len(self.inputs.in_file) > 1 and not isdefined(self.inputs.extra_values):
-            formatstr = '{bname}{space}{desc}{bundle}{suffix}{i:04d}{dtype}{ext}'
+            formatstr = '{bname}{space}{desc}{suffix}{i:04d}{dtype}{ext}'
 
         # Otherwise it's file(s)
         self._results['compression'] = []
@@ -330,7 +331,6 @@ def _run_interface(self, runtime):
                 bname=base_fname,
                 space=space,
                 desc=desc,
-                bundle=bundle,
                 suffix=suffix,
                 i=i,
                 dtype=dtype,
@@ -358,9 +358,149 @@ def _run_interface(self, runtime):
         return super(DerivativesMaybeDataSink, self)._run_interface(runtime)
 
 
+recon_entity_order = [
+    "model",
+    "bundles",
+    "fit",
+    "mdp",
+    "mfp",
+    "bundle",
+    "label"
+]
+
+def get_recon_output_name(base_dir, source_file, derivative_file, qsirecon_suffix, output_bids_entities, use_ext=True):
+
+    source_entities = parse_file_entities(source_file)
+    out_path = op.join(base_dir, f"qsirecon-{qsirecon_suffix}")
+    out_path = op.join(out_path, "sub-" + source_entities["subject"])
+    if "session" in source_entities:
+        out_path += "/ses-{session}".format(**source_entities)
+
+    # Which datatype directory should this go under?
+    # If it's not in the output bids entities use the source_file datatype
+    datatype_dir = output_bids_entities.get(
+        "datatype",
+        op.basename(op.dirname(source_file))
+    )
+    out_path += f"/{datatype_dir}"
+    _, source_fname, _ = split_filename(source_file)
+    source_fname, _ = source_fname.rsplit('_', 1)
+    _, _, extension = split_filename(derivative_file)
+
+    # It may be that the space has changed. Check if it has
+    if "space" in output_bids_entities:
+        source_fname = re.sub(
+            "_space-[a-zA-Z0-9]+",
+            "_space-" + output_bids_entities["space"],
+            source_fname)
+    base_fname = op.join(out_path,source_fname)
+
+    # Add the new bids entities for the output file
+    for entity_name in recon_entity_order:
+        if entity_name in output_bids_entities:
+            base_fname += "_{entity}-{value}".format(
+                entity=entity_name,
+                value=output_bids_entities[entity_name])
+
+    # Add the suffix
+    suffix = output_bids_entities.get("suffix", "dwimap")
+    if use_ext:
+        return f"{base_fname}_{suffix}{extension}"
+
+    return f"{base_fname}_{suffix}"
+
+
+class _ReconDerivativesDataSinkInputSpec(DerivativesDataSinkInputSpec):
+    mdp = traits.Str('', usedefault=True, desc='Label for model derived parameter field')
+    mfp = traits.Str('', usedefault=True, desc='Label for model fit parameter field')
+    model = traits.Str('', usedefault=True, desc='Label for model field')
+    bundle = traits.Str('', usedefault=True, desc='Label for bundle field')
+    bundles = traits.Str('', usedefault=True, desc='Label for bundles field')
+    fit = traits.Str('', usedefault=True, desc='Label for fit field')
+    label = traits.Str('', usedefault=True, desc='Label for label field')
+    qsirecon_suffix = traits.Str('', usedefault=True, desc='name appended to qsirecon- in the derivatives')
+
+
 class ReconDerivativesDataSink(DerivativesDataSink):
+    input_spec = _ReconDerivativesDataSinkInputSpec
     out_path_base = "qsirecon"
 
+    def _run_interface(self, runtime):
+
+        # If there is no qsirecon suffix, then we're not saving this file
+        if not self.inputs.qsirecon_suffix:
+            return runtime
+
+        # Figure out what the extension should be based on the input file and compression
+        src_fname, _ = _splitext(self.inputs.source_file)
+        src_fname, dtype = src_fname.rsplit('_', 1)
+        _, ext = _splitext(self.inputs.in_file[0])
+        if self.inputs.compress is True and not ext.endswith('.gz'):
+            ext += '.gz'
+        elif self.inputs.compress is False and ext.endswith('.gz'):
+            ext = ext[:-3]
+
+        # Prepare the bids entities from the inputs
+        output_bids = {}
+        if self.inputs.bundles:
+            output_bids["bundles"] = self.inputs.bundles
+        if self.inputs.bundle:
+            output_bids["bundle"] = self.inputs.bundle
+        if self.inputs.space:
+            output_bids["space"] = self.inputs.space
+        if self.inputs.model:
+            output_bids["model"] = self.inputs.model
+        if self.inputs.mdp:
+            output_bids["mdp"] = self.inputs.mdp
+        if self.inputs.mfp:
+            output_bids["mfp"] = self.inputs.mfp
+        if self.inputs.fit:
+            output_bids["fit"] = self.inputs.fit
+        if self.inputs.suffix:
+            output_bids["suffix"] = self.inputs.suffix
+        if self.inputs.label:
+            output_bids["label"] = self.inputs.label
+
+        # Get the output name without an extension
+        bname = get_recon_output_name(
+            base_dir=self.inputs.base_directory,
+            source_file=self.inputs.source_file,
+            derivative_file=self.inputs.in_file[0],
+            qsirecon_suffix=self.inputs.qsirecon_suffix,
+            output_bids_entities=output_bids,
+            use_ext=False
+        )
+
+        # Ensure the directory exists
+        os.makedirs(op.dirname(bname), exist_ok=True)
+
+        formatstr = '{bname}{ext}'
+        # If the derivative is a directory, copy it over
+        copy_dir = op.isdir(str(self.inputs.in_file))
+        if copy_dir:
+            out_file = formatstr.format(
+                bname=base_fname,
+                ext='')
+            copytree(str(self.inputs.in_file), out_file, dirs_exist_ok=True)
+            self._results['out_file'] = out_file
+            return runtime
+
+        if len(self.inputs.in_file) > 1 and not isdefined(self.inputs.extra_values):
+            formatstr = '{bname}{i:04d}{ext}'
+
+        # Otherwise it's file(s)
+        self._results['compression'] = []
+        for i, fname in enumerate(self.inputs.in_file):
+            out_file = formatstr.format(
+                bname=bname,
+                i=i,
+                ext=ext)
+            if isdefined(self.inputs.extra_values):
+                out_file = out_file.format(extra_value=self.inputs.extra_values[i])
+            self._results['out_file'].append(out_file)
+            self._results['compression'].append(_copy_any(fname, out_file))
+        return runtime
+
 
 class ReadSidecarJSONInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='the input nifti file')
diff --git a/qsiprep/interfaces/interchange.py b/qsiprep/interfaces/interchange.py
index 52320178..0026d4da 100644
--- a/qsiprep/interfaces/interchange.py
+++ b/qsiprep/interfaces/interchange.py
@@ -29,7 +29,7 @@
 
 # dMRI + registered anatomical fields
 recon_workflow_anatomical_input_fields = anatomical_workflow_outputs + [
-    "dwi_mask", "atlas_configs", "odf_rois"]
+    "dwi_mask", "atlas_configs", "odf_rois", "resampling_template", "mapping_metadata"]
 
 # Check that no conflicts have been introduced
 overlapping_names = set(qsiprep_output_names
diff --git a/qsiprep/interfaces/recon_scalars.py b/qsiprep/interfaces/recon_scalars.py
index a725c560..93478848 100644
--- a/qsiprep/interfaces/recon_scalars.py
+++ b/qsiprep/interfaces/recon_scalars.py
@@ -8,20 +8,25 @@
 
 
 """
+import os
 import os.path as op
-
+import re
+from pkg_resources import resource_filename as pkgr
+import pandas as pd
+from nipype import logging
+from bids.layout import parse_file_entities
+from .bids import _copy_any, get_recon_output_name
+from nipype.utils.filemanip import fname_presuffix, split_filename
 from nipype.interfaces.base import (
-    BaseInterfaceInputSpec,
-    File,
-    SimpleInterface,
-    TraitedSpec,
-    isdefined,
-    traits,
+    InputMultiObject, traits, TraitedSpec, BaseInterfaceInputSpec, SimpleInterface, isdefined, File
 )
 
 
 class ReconScalarsInputSpec(BaseInterfaceInputSpec):
-    workflow_name = traits.Str(mandatory=True)
+    source_file = File(exists=True, mandatory=True)
+    qsirecon_suffix = traits.Str(mandatory=True)
+    model_info = traits.Dict()
+    model_name = traits.Str()
 
 
 class ReconScalarsOutputSpec(TraitedSpec):
@@ -32,7 +37,9 @@ class ReconScalars(SimpleInterface):
     input_spec = ReconScalarsInputSpec
     output_spec = ReconScalarsOutputSpec
     scalar_metadata = {}
-    _ignore_traits = ("workflow_name", "scalar_metadata")
+    _ignore_traits = ("model_name", "qsirecon_suffix", "scalar_metadata",
+                      "model_info", "source_file")
+
     def __init__(self, from_file=None, resource_monitor=None, **inputs):
 
         # Get self._results defined
@@ -49,64 +56,153 @@ def _validate_scalar_metadata(self):
             if not input_key in self.scalar_metadata:
                 raise Exception(f"No entry found for {input_key} in ``scalar_metadata`` in this class.")
 
+            # Check that BIDS attributes are defined
+            if "bids" not in self.scalar_metadata[input_key]:
+                raise Exception(f"Missing BIDS metadata for {input_key}")
+
     def _run_interface(self, runtime):
         results = []
         inputs = self.inputs.get()
 
+        # Get the BIDS info from the source file
+        source_file_bids = parse_file_entities(self.inputs.source_file)
+        del source_file_bids["extension"], source_file_bids["suffix"]
+
         file_traits = [name for name in self.inputs.editable_traits()
                        if name not in self._ignore_traits]
+
         for input_name in file_traits:
             if not isdefined(inputs[input_name]):
                 continue
             result = self.scalar_metadata[input_name].copy()
             result["path"] = op.abspath(inputs[input_name])
-            result["workflow_name"] = self.inputs.workflow_name
+            result["qsirecon_suffix"] = self.inputs.qsirecon_suffix
             result["variable_name"] = input_name
+            result["source_file"] = self.inputs.source_file
+            # Update the BIDS with the source file
+            bids_overlap = set(source_file_bids.keys()).intersection(result['bids'].keys())
+            if bids_overlap:
+                raise Exception(f"BIDS fields for {input_name} conflict with source file BIDS {bids_overlap}")
             results.append(result)
         self._results["scalar_info"] = results
         return runtime
 
 
+class _ReconScalarsDataSinkInputSpec(BaseInterfaceInputSpec):
+    source_file = File()
+    base_directory = File()
+    resampled_files = InputMultiObject(File(exists=True))
+    recon_scalars = InputMultiObject(traits.Any())
+
+
+class ReconScalarsDataSink(SimpleInterface):
+    input_spec = _ReconScalarsDataSinkInputSpec
+    _always_run = True
+
+    def _run_interface(self, runtime):
+
+        for recon_scalar in self.inputs.recon_scalars:
+            output_filename = get_recon_output_name(
+                base_dir=self.inputs.base_directory,
+                source_file=self.inputs.source_file,
+                derivative_file=recon_scalar["path"],
+                qsirecon_suffix=recon_scalar["qsirecon_suffix"],
+                output_bids_entities=recon_scalar["bids"],
+                use_ext=True)
+            output_dir = op.dirname(output_filename)
+            os.makedirs(output_dir, exist_ok=True)
+            _copy_any(recon_scalar["path"], output_filename)
+
+        return runtime
+
+
+class _ReconScalarsTableSplitterDataSinkInputSpec(_ReconScalarsDataSinkInputSpec):
+    summary_tsv = File(exists=True, mandatory=True,
+                       desc="tsv of combined scalar summaries")
+    suffix = traits.Str(mandatory=True)
+
+
+class ReconScalarsTableSplitterDataSink(ReconScalarsDataSink):
+    input_spec = _ReconScalarsTableSplitterDataSinkInputSpec
+    def _run_interface(self, runtime):
+        summary_df = pd.read_csv(self.inputs.summary_tsv, sep="\t")
+        for groupname, group_df in summary_df.groupby("qsirecon_suffix"):
+
+            # reset the index for this df
+            group_df.reset_index(drop=True, inplace=True)
+
+            qsirecon_suffixed_tsv = get_recon_output_name(
+                base_dir=self.inputs.base_directory,
+                source_file=group_df.loc[0, "source_file"],
+                derivative_file=self.inputs.summary_tsv,
+                qsirecon_suffix=group_df.loc[0, "qsirecon_suffix"],
+                output_bids_entities={
+                    "suffix": self.inputs.suffix,
+                    "bundles": group_df.loc[0, "bundle_source"]}
+            )
+            output_dir = op.dirname(qsirecon_suffixed_tsv)
+            os.makedirs(output_dir, exist_ok=True)
+            group_df.to_csv(
+                qsirecon_suffixed_tsv,
+                index=False,
+                sep="\t")
+
+        return runtime
+
+
 # Scalars produced in the TORTOISE recon workflow
 tortoise_scalars = {
     "fa_file": {
-        "desc": "Fractional Anisotropy from a tensor fit"
+        "desc": "Fractional Anisotropy from a tensor fit",
+        "bids":{"mdp": "fa", "model": "tensor"}
     },
     "rd_file": {
-        "desc": "Radial Diffusivity from a tensor fit"
+        "desc": "Radial Diffusivity from a tensor fit",
+        "bids":{"mdp": "rd", "model": "tensor"}
     },
     "ad_file": {
-        "desc": "Apparent Diffusivity from a tensor fit"
+        "desc": "Apparent Diffusivity from a tensor fit",
+        "bids":{"mdp": "ad", "model": "tensor"}
     },
     "li_file": {
-        "desc": "LI from a tensor fit"
+        "desc": "LI from a tensor fit",
+        "bids":{"mdp": "li", "model": "tensor"}
     },
     "am_file": {
-        "desc": "A0 from a tensor fit"
+        "desc": "A0 from a tensor fit",
+        "bids":{"mfp": "AM", "model": "tensor"}
     },
     "pa_file": {
-        "desc": "PA from MAPMRI"
+        "desc": "PA from MAPMRI",
+        "bids":{"mdp": "PA", "model": "mapmri"}
     },
     "path_file": {
-        "desc": "PAth from MAPMRI"
+        "desc": "PAth from MAPMRI",
+        "bids":{"mdp": "PAth", "model": "mapmri"}
     },
     "rtop_file": {
-        "desc": "Return to origin probability from MAPMRI"
+        "desc": "Return to origin probability from MAPMRI",
+        "bids":{"mdp": "RTOP", "model": "mapmri"}
     },
     "rtap_file": {
-        "desc": "Return to axis probability from MAPMRI"
+        "desc": "Return to axis probability from MAPMRI",
+        "bids":{"mdp": "RTAP", "model": "mapmri"}
     },
     "rtpp_file": {
-        "desc": "Return to plane probability from MAPMRI"
+        "desc": "Return to plane probability from MAPMRI",
+        "bids":{"mdp": "RTPP", "model": "mapmri"}
     },
     "ng_file": {
-        "desc": "Non-Gaussianity from MAPMRI"
+        "desc": "Non-Gaussianity from MAPMRI",
+        "bids":{"mdp": "NG", "model": "mapmri"}
     },
     "ngpar_file": {
-        "desc": "Non-Gaussianity parallel from MAPMRI"
+        "desc": "Non-Gaussianity parallel from MAPMRI",
+        "bids":{"mdp": "NGpar", "model": "mapmri"}
     },
     "ngperp_file": {
-        "desc": "Non-Gaussianity perpendicular from MAPMRI"
+        "desc": "Non-Gaussianity perpendicular from MAPMRI",
+        "bids":{"mdp": "NGperp", "model": "mapmri"}
     },
 }
 
@@ -124,13 +220,21 @@ class TORTOISEReconScalars(ReconScalars):
 # Scalars produced in the AMICO recon workflow
 amico_scalars = {
     "icvf_image": {
-        "desc": "Intracellular volume fraction from NODDI"
+        "desc": "Intracellular volume fraction from NODDI",
+        "bids":{"mdp": "icvf", "model": "noddi"}
     },
     "isovf_image": {
-        "desc": "Isotropic volume fraction from NODDI"
+        "desc": "Isotropic volume fraction from NODDI",
+        "bids":{"mdp": "isovf", "model": "noddi"}
     },
     "od_image": {
-        "desc": "OD from NODDI"
+        "desc": "OD from NODDI",
+        "bids":{"mdp": "od", "model": "noddi"}
+    },
+    "directions_image": {
+        "desc": "Peak directions from NODDI",
+        "reorient_on_resample": True,
+        "bids":{"mfp": "direction", "model": "noddi"}
     }
 }
 
@@ -148,64 +252,84 @@ class AMICOReconScalars(ReconScalars):
 # Scalars produced by DSI Studio
 dsistudio_scalars = {
     "qa_file": {
-        "desc": "Fractional Anisotropy from a tensor fit"
+        "desc": "Fractional Anisotropy from a tensor fit",
+        "bids":{"mdp": "qa", "model": "GQI"}
     },
     "dti_fa_file": {
-        "desc": "Radial Diffusivity from a tensor fit"
+        "desc": "Radial Diffusivity from a tensor fit",
+        "bids":{"mdp": "qa", "model": "tensor"}
     },
     "txx_file": {
-        "desc": "Tensor fit txx"
+        "desc": "Tensor fit txx",
+        "bids":{"mfp": "txx", "model": "tensor"}
     },
     "txy_file": {
-        "desc": "Tensor fit txy"
+        "desc": "Tensor fit txy",
+        "bids":{"mfp": "txx", "model": "tensor"}
     },
     "txz_file": {
-        "desc": "Tensor fit txz"
+        "desc": "Tensor fit txz",
+        "bids":{"mfp": "txx", "model": "tensor"}
     },
     "tyy_file": {
-        "desc": "Tensor fit tyy"
+        "desc": "Tensor fit tyy",
+        "bids":{"mfp": "txx", "model": "tensor"}
     },
     "tyz_file": {
-        "desc": "Tensor fit tyz"
+        "desc": "Tensor fit tyz",
+        "bids":{"mfp": "txx", "model": "tensor"}
     },
     "tzz_file": {
-        "desc": "Tensor fit tzz"
+        "desc": "Tensor fit tzz",
+        "bids":{"mfp": "txx", "model": "tensor"}
     },
     "rd1_file": {
-        "desc": "RD1"
+        "desc": "RD1",
+        "bids":{"mdp": "rd1", "model": "RDI"}
     },
     "rd2_file": {
-        "desc": "RD2"
+        "desc": "RD2",
+        "bids":{"mdp": "rd2", "model": "RDI"}
     },
     "ha_file": {
-        "desc": "HA"
+        "desc": "HA",
+        "bids":{"mdp": "ha", "model": "RDI"}
     },
     "md_file": {
-        "desc": "Mean Diffusivity"
+        "desc": "Mean Diffusivity",
+        "bids":{"mdp": "rd2", "model": "RDI"}
     },
     "ad_file": {
-        "desc": "AD"
+        "desc": "AD",
+        "bids":{"mdp": "ad", "model": "tensor"}
     },
     "rd_file": {
-        "desc": "Radial Diffusivity"
+        "desc": "Radial Diffusivity",
+        "bids":{"mdp": "rd", "model": "tensor"}
     },
     "gfa_file": {
-        "desc": "Generalized Fractional Anisotropy"
+        "desc": "Generalized Fractional Anisotropy",
+        "bids": {"mdp": "gfa", "fit": "GQI"}
     },
     "iso_file": {
-        "desc": "Isotropic Diffusion"
+        "desc": "Isotropic Diffusion",
+        "bids": {"mdp": "iso", "fit": "GQI"}
     },
     "rdi_file": {
-        "desc": "RDI"
+        "desc": "RDI",
+        "bids":{"mdp": "rdi", "model": "RDI"}
     },
     "nrdi02L_file": {
-        "desc": "NRDI at 02L"
+        "desc": "NRDI at 02L",
+        "bids":{"mdp": "nrdi02L", "model": "RDI"}
     },
     "nrdi04L_file": {
-        "desc": "NRDI at 04L"
+        "desc": "NRDI at 04L",
+        "bids":{"mdp": "nrdi04L", "model": "RDI"}
     },
     "nrdi06L_file": {
-        "desc": "NRDI at 06L"
+        "desc": "NRDI at 06L",
+        "bids":{"mdp": "nrdi06L", "model": "RDI"}
     },
 }
 
@@ -222,33 +346,43 @@ class DSIStudioReconScalars(ReconScalars):
 
 dipy_dki_scalars = {
     'dki_fa': {
-        "desc": "DKI FA"
+        "desc": "DKI FA",
+        "bids": {"mdp": "FA", "model": "tensor"}
     },
     'dki_md': {
-        "desc": "DKI MD"
+        "desc": "DKI MD",
+        "bids": {"mdp": "MD", "model": "dki"}
     },
     'dki_rd': {
-        "desc": "DKI RD"
+        "desc": "DKI RD",
+        "bids": {"mdp": "RD", "model": "dki"}
     },
     'dki_ad': {
-        "desc": "DKI AD"
+        "desc": "DKI AD",
+        "bids": {"mdp": "AD", "model": "dki"}
     },
     'dki_kfa': {
-        "desc": "DKI KFA"
+        "desc": "DKI KFA",
+        "bids": {"mdp": "KFA", "model": "dki"}
     },
     'dki_mk': {
-        "desc": "DKI MK"
+        "desc": "DKI MK",
+        "bids": {"mdp": "MK", "model": "dki"}
     },
     'dki_ak': {
-        "desc": "DKI AK"
+        "desc": "DKI AK",
+        "bids": {"mdp": "AK", "model": "dki"}
     },
     'dki_rk': {
-        "desc": "DKI RK"
+        "desc": "DKI RK",
+        "bids": {"mdp": "RK", "model": "dki"}
     },
     'dki_mkt': {
-        "desc": "DKI MKT"
+        "desc": "DKI MKT",
+        "bids": {"mdp": "MKT", "model": "dki"}
     }
 }
+
 class _DIPYDKIReconScalarInputSpec(ReconScalarsInputSpec):
     pass
 
@@ -263,31 +397,40 @@ class DIPYDKIReconScalars(ReconScalars):
 # DIPY implementation of MAPMRI
 dipy_mapmri_scalars = {
     "qiv_file": {
-        "desc": "q-space inverse variance from MAPMRI"
+        "desc": "q-space inverse variance from MAPMRI",
+        "bids":{"mdp": "QIV", "model": "mapmri"}
     },
     "msd_file": {
-        "desc": "mean square displacement from MAPMRI"
+        "desc": "mean square displacement from MAPMRI",
+        "bids":{"mdp": "MSD", "model": "mapmri"}
     },
     "lapnorm_file": {
-        "desc": "Laplacian norm from regularized MAPMRI (MAPL)"
+        "desc": "Laplacian norm from regularized MAPMRI (MAPL)",
+        "bids":{"mfp": "lapnorm", "model": "mapmri"}
     },
     "rtop_file": {
-        "desc": "Return to origin probability from MAPMRI"
+        "desc": "Return to origin probability from MAPMRI",
+        "bids":{"mdp": "RTOP", "model": "mapmri"}
     },
     "rtap_file": {
-        "desc": "Return to axis probability from MAPMRI"
+        "desc": "Return to axis probability from MAPMRI",
+        "bids":{"mdp": "RTAP", "model": "mapmri"}
     },
     "rtpp_file": {
-        "desc": "Return to plane probability from MAPMRI"
+        "desc": "Return to plane probability from MAPMRI",
+        "bids":{"mdp": "RTPP", "model": "mapmri"}
     },
     "ng_file": {
-        "desc": "Non-Gaussianity from MAPMRI"
+        "desc": "Non-Gaussianity from MAPMRI",
+        "bids":{"mdp": "NG", "model": "mapmri"}
     },
     "ngpar_file": {
-        "desc": "Non-Gaussianity parallel from MAPMRI"
+        "desc": "Non-Gaussianity parallel from MAPMRI",
+        "bids":{"mdp": "NGpar", "model": "mapmri"}
     },
     "ngperp_file": {
-        "desc": "Non-Gaussianity perpendicular from MAPMRI"
+        "desc": "Non-Gaussianity perpendicular from MAPMRI",
+        "bids":{"mdp": "NGperp", "model": "mapmri"}
     },
 }
 
diff --git a/qsiprep/interfaces/scalar_mapping.py b/qsiprep/interfaces/scalar_mapping.py
index 40518de9..4b70ccbc 100644
--- a/qsiprep/interfaces/scalar_mapping.py
+++ b/qsiprep/interfaces/scalar_mapping.py
@@ -10,14 +10,14 @@
 import pandas as pd
 from nilearn.maskers import NiftiMasker
 from nipype.interfaces.base import (
-    BaseInterfaceInputSpec,
-    File,
-    InputMultiObject,
-    SimpleInterface,
-    TraitedSpec,
-    traits,
+    traits, TraitedSpec, BaseInterfaceInputSpec, File, SimpleInterface, isdefined,
+    InputMultiObject, OutputMultiObject
 )
+from nipype.interfaces import ants
 from nipype.utils.filemanip import fname_presuffix
+import nilearn.image as nim
+from nilearn.maskers import NiftiMasker
+import pandas as pd
 
 from .bids import get_bids_params
 
@@ -26,6 +26,8 @@ class ScalarMapperInputSpec(BaseInterfaceInputSpec):
     scalars_from = InputMultiObject(traits.Str())
     recon_scalars = InputMultiObject(traits.Any())
     dwiref_image = File(exists=True)
+    mapping_metadata = traits.Dict(
+        desc="Info about the upstream workflow that created the anatomical mapping units")
 
 
 class ScalarMapperOutputSpec(TraitedSpec):
@@ -47,8 +49,12 @@ def _update_with_bids_info(self, summary_row_list):
         for summary_row in summary_row_list:
             summary_row.update(bids_info)
 
+    def _unload_scalars(self):
+        for recon_scalar in self._results["mapped_scalars"]:
+            if "image" in recon_scalar:
+                del recon_scalar["image"]
+
     def _run_interface(self, runtime):
-        self._load_scalars()
         self._do_mapping(runtime)
         return runtime
 
@@ -63,7 +69,7 @@ class _BundleMapperInputSpec(ScalarMapperInputSpec):
 
 class _BundleMapperOutputSpec(ScalarMapperOutputSpec):
     bundle_summary = File(exists=True)
-    bundle_profiles = File(exists=True)
+    tdi_stats = File(exists=True)
 
 
 def _get_tdi_img(dwiref_image, tck_file, output_tdi_file):
@@ -73,12 +79,16 @@ def _get_tdi_img(dwiref_image, tck_file, output_tdi_file):
         check=True)
     return nim.load_img(output_tdi_file)
 
+
 class BundleMapper(ScalarMapper):
     input_spec = _BundleMapperInputSpec
     output_spec = _BundleMapperOutputSpec
 
     def _do_mapping(self, runtime):
+        self._load_scalars()
         bundle_dfs = []
+        tdi_dfs = []
+        source_suffix = self.inputs.mapping_metadata.get("qsirecon_suffix", "QSIPrep")
         for tck_name, tck_file in zip(self.inputs.bundle_names, self.inputs.tck_files):
             output_tdi_file = fname_presuffix(
                 tck_file,
@@ -101,28 +111,46 @@ def _do_mapping(self, runtime):
             tdi_weights = tdi_weights / tdi_weights.sum()
 
             # Start gathering stats with the TDI first
-            bundle_dfs.append(
+            tdi_dfs.append(
                 calculate_mask_stats(
                     bundle_masker,
                     tck_name,
                     "bundle",
                     {"image": tdi_img,
                      "variable_name": "tdi",
-                     "workflow_name": "scalar_to_bundle",
+                     # Check that this is ok:
+                     "source_file": self.inputs.recon_scalars[0]["source_file"],
+                     "qsirecon_suffix": source_suffix,
                      "desc": "Streamline counts per voxel"})
             )
 
             # Then get the same stats for the scalars
-            for recon_scalar in self.recon_scalars:
+            for recon_scalar in self.inputs.recon_scalars:
                 bundle_dfs.append(
                     calculate_mask_stats(bundle_masker, tck_name, "bundle",
                                          recon_scalar, tdi_weights))
+
+        # Write the scalar summary df
         self._update_with_bids_info(bundle_dfs)
         summary_file = op.join(runtime.cwd, "bundle_stats.tsv")
         summary_df = pd.DataFrame(bundle_dfs)
-
+        # Add information about which bundle workflow created these summaries
+        if isdefined(self.inputs.mapping_metadata):
+            summary_df["bundle_source"] = self.inputs.mapping_metadata.get("qsirecon_suffix")
+            summary_df["bundle_params_id"] = self.inputs.mapping_metadata.get("name")
         summary_df.to_csv(summary_file, index=False, sep="\t")
-        self._results['bundle_summary'] = summary_file
+        self._results["bundle_summary"] = summary_file
+
+        # Write the TDI df
+        self._update_with_bids_info(tdi_dfs)
+        tdi_file = op.join(runtime.cwd, "tdi_stats.tsv")
+        tdi_summary_df = pd.DataFrame(tdi_dfs)
+        # Add information about which bundle workflow created these summaries
+        if isdefined(self.inputs.mapping_metadata):
+            tdi_summary_df["bundle_source"] = self.inputs.mapping_metadata.get("qsirecon_suffix")
+            tdi_summary_df["bundle_params_id"] = self.inputs.mapping_metadata.get("name")
+        tdi_summary_df.to_csv(tdi_file, index=False, sep="\t")
+        self._results['tdi_stats'] = tdi_file
 
 
 # For mapping to atlases
@@ -138,6 +166,7 @@ class AtlasMapper(ScalarMapper):
     input_spec = _AtlasMapperInputSpec
     output_spec = _AtlasMapperOutputSpec
 
+
 def calculate_mask_stats(masker, mask_name, mask_variable_name, recon_scalar, weighting_vector=None):
 
     # Get the scalar data in the masked region
@@ -153,7 +182,8 @@ def calculate_mask_stats(masker, mask_name, mask_variable_name, recon_scalar, we
     results = {
         mask_variable_name: mask_name,
         "variable_name": variable_name,
-        "workflow": recon_scalar["workflow_name"],
+        "qsirecon_suffix": recon_scalar["qsirecon_suffix"],
+        "source_file": recon_scalar["source_file"],
         "zero_proportion": np.sum(np.isnan(nz_voxel_data)) / voxel_data.shape[0],
         "mean": np.mean(voxel_data),
         "stdev": np.std(voxel_data),
@@ -171,4 +201,54 @@ def calculate_mask_stats(masker, mask_name, mask_variable_name, recon_scalar, we
         results["masked_weighted_mean"] = np.nansum(
             nz_voxel_data * nz_weighting_vector )
 
-    return results
\ No newline at end of file
+    return results
+
+
+class _TemplateMapperInputSpec(ScalarMapperInputSpec):
+    template_reference_image = File(exists=True, mandatory=True)
+    to_template_transform = File(exists=True, mandatory=True)
+    interpolation = traits.Str("NearestNeighbor", usedefault=True)
+
+
+class _TemplateMapperOutputSpec(ScalarMapperOutputSpec):
+    template_space_scalars = OutputMultiObject(traits.Any())
+    template_space_scalar_info = OutputMultiObject(traits.Any())
+
+
+class TemplateMapper(ScalarMapper):
+    input_spec = _TemplateMapperInputSpec
+    output_spec = _TemplateMapperOutputSpec
+
+    def _do_mapping(self, runtime):
+        resampled_images = []
+        resampled_image_metadata = []
+        # Then get the same stats for the scalars
+        for recon_scalar in self.inputs.recon_scalars:
+            if recon_scalar.get("reorient_on_resample", False):
+                # LOGGER.info(f"Skipping {recon_scalar}")
+                continue
+            new_metadata = recon_scalar.copy()
+            output_fname = op.split(recon_scalar["path"])[1]
+            output_fname = output_fname.replace("_space-T1w_", "_transformed_")
+            output_fname = op.join(runtime.cwd, output_fname)
+            transform = ants.ApplyTransforms(
+                input_image=recon_scalar["path"],
+                dimension=3,
+                transforms=[self.inputs.to_template_transform],
+                reference_image=self.inputs.template_reference_image,
+                output_image=output_fname,
+                interpolation=self.inputs.interpolation)
+            transform.terminal_output = 'allatonce'
+            transform.resource_monitor = False
+            transform.run()
+            resampled_images.append(output_fname)
+
+            # Create new metadata for the resampled image
+            new_metadata["path"] = output_fname
+            if "bids" not in new_metadata:
+                raise Exception(f"incomplete metadata spec {new_metadata}")
+            new_metadata["bids"]["space"] = "MNI152NLin2009cAsym"
+            resampled_image_metadata.append(new_metadata)
+
+        self._results['template_space_scalars'] = resampled_images
+        self._results['template_space_scalar_info'] = resampled_image_metadata
diff --git a/qsiprep/workflows/anatomical/volume.py b/qsiprep/workflows/anatomical/volume.py
index 2e8fa25c..2604005a 100644
--- a/qsiprep/workflows/anatomical/volume.py
+++ b/qsiprep/workflows/anatomical/volume.py
@@ -29,7 +29,11 @@
 from ...niworkflows.interfaces.masks import ROIsPlot
 from ...niworkflows.interfaces.registration import RobustMNINormalizationRPT
 from ...utils.misc import fix_multi_source_name
-from qsiprep.interfaces import Conform
+from ...interfaces.freesurfer import (
+        PrepareSynthStripGrid, FixHeaderSynthStrip, SynthSeg)
+from ...interfaces.anatomical import DesaturateSkull, GetTemplate, VoxelSizeChooser
+from ...interfaces.itk import DisassembleTransform, AffineToRigid
+from ...interfaces import Conform
 
 LOGGER = logging.getLogger('nipype.workflow')
 
@@ -963,26 +967,35 @@ def init_synthseg_wf(omp_nthreads, sloppy, name="synthseg_wf"):
 def init_output_grid_wf(voxel_size, padding, name='output_grid_wf'):
     """Generate a non-oblique, uniform voxel-size grid around a brain."""
     workflow = Workflow(name=name)
-    inputnode = pe.Node(niu.IdentityInterface(fields=['template_image']), name='inputnode')
+    inputnode = pe.Node(niu.IdentityInterface(fields=['template_image', 'input_image']), name='inputnode')
     outputnode = pe.Node(niu.IdentityInterface(fields=['grid_image']), name='outputnode')
     autobox_template = pe.Node(afni.Autobox(outputtype="NIFTI_GZ", padding=padding),
                                name='autobox_template')
     deoblique_autobox = pe.Node(afni.Warp(outputtype="NIFTI_GZ", deoblique=True),
                                 name="deoblique_autobox")
+    voxel_size_chooser = pe.Node(
+        VoxelSizeChooser(voxel_size=voxel_size),
+        name="voxel_size_chooser")
     resample_to_voxel_size = pe.Node(afni.Resample(outputtype="NIFTI_GZ"),
                                      name="resample_to_voxel_size")
-    resample_to_voxel_size.inputs.voxel_size = (voxel_size, voxel_size, voxel_size)
 
     workflow.connect([
         (inputnode, autobox_template, [('template_image', 'in_file')]),
         (autobox_template, deoblique_autobox, [('out_file', 'in_file')]),
         (deoblique_autobox, resample_to_voxel_size, [('out_file', 'in_file')]),
-        (resample_to_voxel_size, outputnode, [('out_file', 'grid_image')])
+        (resample_to_voxel_size, outputnode, [('out_file', 'grid_image')]),
+        (inputnode, voxel_size_chooser, [('input_image', 'input_image')]),
+        (voxel_size_chooser, resample_to_voxel_size, [(('voxel_size', _tupleize), 'voxel_size')])
     ])
 
     return workflow
 
 
+def _tupleize(value):
+    # Nipype did not like having a Tuple output trait
+    return (value, value, value)
+
+
 def init_anat_reports_wf(reportlets_dir, nonlinear_register_to_template,
                          name='anat_reports_wf'):
     """
diff --git a/qsiprep/workflows/recon/amico.py b/qsiprep/workflows/recon/amico.py
index fb992ef4..f252a7c8 100644
--- a/qsiprep/workflows/recon/amico.py
+++ b/qsiprep/workflows/recon/amico.py
@@ -14,17 +14,19 @@
 from ...engine import Workflow
 from ...interfaces.amico import NODDI
 from ...interfaces.converters import NODDItoFIBGZ
+from ...interfaces.recon_scalars import AMICOReconScalars, ReconScalarsDataSink
 from ...interfaces.interchange import recon_workflow_input_fields
 from ...interfaces.recon_scalars import AMICOReconScalars
 from ...interfaces.reports import CLIReconPeaksReport
 from qsiprep.interfaces.bids import ReconDerivativesDataSink
 
+
 LOGGER = logging.getLogger('nipype.interface')
 
 
 def init_amico_noddi_fit_wf(omp_nthreads, available_anatomical_data,
                             name="amico_noddi_recon",
-                            output_suffix="", params={}):
+                            qsirecon_suffix="", params={}):
     """Reconstruct EAPs, ODFs, using 3dSHORE (brainsuite-style basis set).
 
     Inputs
@@ -56,7 +58,7 @@ def init_amico_noddi_fit_wf(omp_nthreads, available_anatomical_data,
         name="outputnode")
 
     workflow = Workflow(name=name)
-    recon_scalars = pe.Node(AMICOReconScalars(workflow_name=name),
+    recon_scalars = pe.Node(AMICOReconScalars(qsirecon_suffix=qsirecon_suffix),
                             name="recon_scalars",
                             run_without_submitting=True)
     plot_reports = params.pop("plot_reports", True)
@@ -92,6 +94,7 @@ def init_amico_noddi_fit_wf(omp_nthreads, available_anatomical_data,
             ('icvf_image', 'icvf_image'),
             ('od_image', 'od_image'),
             ('isovf_image', 'isovf_image'),
+            ('directions_image', 'directions_image')
         ]),
         (recon_scalars, outputnode, [("scalar_info", "recon_scalars")]),
         (noddi_fit, convert_to_fibgz, [
@@ -119,56 +122,29 @@ def init_amico_noddi_fit_wf(omp_nthreads, available_anatomical_data,
             (noddi_fit, plot_peaks, [('icvf_image', 'background_image')]),
             (plot_peaks, ds_report_peaks, [('peak_report', 'in_file')])])
 
-    if output_suffix:
+    if qsirecon_suffix:
         ds_fibgz = pe.Node(
         ReconDerivativesDataSink(extension='.fib.gz',
-                                    suffix=output_suffix,
-                                    compress=True),
-        name='ds_{}_fibgz'.format(output_suffix),
+                                 qsirecon_suffix=qsirecon_suffix,
+                                 compress=True),
+        name='ds_{}_fibgz'.format(qsirecon_suffix),
         run_without_submitting=True)
         workflow.connect(outputnode, 'fibgz', ds_fibgz, 'in_file')
 
-        # Niftis from AMICO
-        ds_directions = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                     desc="directions",
-                                     suffix=output_suffix,
-                                     compress=True),
-            name='ds_noddi_directions',
-            run_without_submitting=True)
-        workflow.connect(outputnode, 'directions_image', ds_directions, 'in_file')
-
-        ds_icvf = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                     desc="ICVF",
-                                     suffix=output_suffix,
-                                     compress=True),
-            name='ds_noddi_icvf',
-            run_without_submitting=True)
-        workflow.connect(outputnode, 'icvf_image', ds_icvf, 'in_file')
-
-        ds_isovf = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                     desc="ISOVF",
-                                     suffix=output_suffix,
-                                     compress=True),
-            name='ds_noddi_isovf',
-            run_without_submitting=True)
-        workflow.connect(outputnode, 'isovf_image', ds_isovf, 'in_file')
-
-        ds_od = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                     desc="OD",
-                                     suffix=output_suffix,
-                                     compress=True),
-            name='ds_noddi_od',
+        ds_recon_scalars = pe.Node(
+            ReconScalarsDataSink(),
+            name="ds_recon_scalars",
             run_without_submitting=True)
-        workflow.connect(outputnode, 'od_image', ds_od, 'in_file')
+        workflow.connect(
+            recon_scalars,
+            "scalar_info",
+            ds_recon_scalars,
+            "recon_scalars")
 
         ds_config = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                     desc="AMICOconfig",
-                                     suffix=output_suffix,
+            ReconDerivativesDataSink(mfp="AMICOconfig",
+                                     model="NODDI",
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_noddi_config',
             run_without_submitting=True)
diff --git a/qsiprep/workflows/recon/anatomical.py b/qsiprep/workflows/recon/anatomical.py
index 192956b8..f422b920 100644
--- a/qsiprep/workflows/recon/anatomical.py
+++ b/qsiprep/workflows/recon/anatomical.py
@@ -13,7 +13,12 @@
 
 import logging
 from pathlib import Path
-
+from pkg_resources import resource_filename as pkgrf
+from nipype.pipeline import engine as pe
+from nipype.interfaces import utility as niu
+from nipype.interfaces import mrtrix3, ants, afni
+from nipype.interfaces.base import traits
+from ...niworkflows.interfaces.registration import RobustMNINormalizationRPT
 import nipype.interfaces.io as nio
 from nipype.interfaces import afni, ants, mrtrix3
 from nipype.interfaces import utility as niu
@@ -36,6 +41,9 @@
 from ...niworkflows.interfaces.registration import RobustMNINormalizationRPT
 from qsiprep.interfaces.utils import GetConnectivityAtlases
 
+from ..anatomical.volume import init_output_grid_wf
+
+
 LOGGER = logging.getLogger('nipype.workflow')
 
 # Required freesurfer files for mrtrix's HSV 5tt generation
@@ -68,7 +76,7 @@
 
 def init_highres_recon_anatomical_wf(
     subject_id, recon_input_dir, extras_to_make, freesurfer_dir,
-    needs_t1w_transform, pipeline_source, name='recon_anatomical_wf'):
+    needs_t1w_transform, pipeline_source, infant_mode, name='recon_anatomical_wf'):
     """Gather any high-res anatomical data (images, transforms, segmentations) to use in recon workflows.
 
     This workflow searches through input data to see what anatomical data is available. The anatomical data
@@ -94,7 +102,7 @@ def init_highres_recon_anatomical_wf(
             subject_id, recon_input_dir, name='gather_ukb_anatomical_wf')
     else:
         raise Exception(f"Unknown pipeline source '{pipeline_source}'")
-
+    anat_ingress_node.inputs.infant_mode = infant_mode
     if needs_t1w_transform and not status["has_qsiprep_t1w_transforms"]:
         raise Exception("Cannot compute to-template")
 
@@ -199,7 +207,7 @@ def init_highres_recon_anatomical_wf(
 
 
 def gather_ukb_anatomical_data(
-        subject_id, recon_input_dir, name='gather_ukb_anatomical_wf'):
+        subject_id, recon_input_dir, infant_mode, name='gather_ukb_anatomical_wf'):
     """
     Check a UKB directory for the necessary files for recon workflows.
 
@@ -420,7 +428,7 @@ def init_dwi_recon_anatomical_workflow(
     atlas_names, omp_nthreads, has_qsiprep_5tt_hsvs, needs_t1w_transform,
     has_freesurfer_5tt_hsvs, has_qsiprep_t1w, has_qsiprep_t1w_transforms,
     infant_mode, has_freesurfer, extras_to_make, freesurfer_dir, b0_threshold,
-    sloppy=False, prefer_dwi_mask=False, name="qsirecon_anat_wf"):
+    output_resolution, sloppy=False, prefer_dwi_mask=False, name="qsirecon_anat_wf"):
     """Ensure that anatomical data is available for the reconstruction workflows.
 
     This workflow calculates images/transforms that require a DWI spatial reference.
@@ -460,7 +468,8 @@ def _exchange_fields(fields):
         connect_from_inputnode.difference_update(fields)
         connect_from_buffernode.update(fields)
 
-    _exchange_fields(['dwi_mask', 'atlas_configs', 'odf_rois'])
+    # These are always created here
+    _exchange_fields(['dwi_mask', 'atlas_configs', 'odf_rois', 'resampling_template'])
 
     outputnode = pe.Node(
         niu.IdentityInterface(
@@ -475,6 +484,18 @@ def _get_status():
                 "has_qsiprep_t1w": has_qsiprep_t1w,
                 "has_qsiprep_t1w_transforms": has_qsiprep_t1w_transforms}
 
+    # Create the output reference grid_image
+    if output_resolution is None:
+        output_resolution = traits.Undefined
+
+    reference_grid_wf = init_output_grid_wf(voxel_size=output_resolution,
+                                            padding=4 if infant_mode else 8)
+    workflow.connect([
+        (inputnode, reference_grid_wf, [
+            ('template_image', 'inputnode.template_image'),
+            ('dwi_ref', 'inputnode.input_image')]),
+        (reference_grid_wf, buffernode, [('outputnode.grid_image', 'resampling_template')])])
+
     # Missing Freesurfer AND QSIPrep T1ws, or the user wants a DWI-based mask
     if not (has_qsiprep_t1w or has_freesurfer) or prefer_dwi_mask:
         desc += "No T1w weighted images were available for masking, so a mask " \
diff --git a/qsiprep/workflows/recon/base.py b/qsiprep/workflows/recon/base.py
index 6be667cc..2d42a69d 100644
--- a/qsiprep/workflows/recon/base.py
+++ b/qsiprep/workflows/recon/base.py
@@ -49,7 +49,8 @@
 
 def init_qsirecon_wf(subject_list, run_uuid, work_dir, output_dir, recon_input,
                      recon_spec, low_mem, omp_nthreads, sloppy, freesurfer_input,
-                     b0_threshold, skip_odf_plots, pipeline_source, name="qsirecon_wf"):
+                     b0_threshold, skip_odf_plots, pipeline_source, infant_mode,
+                     output_resolution, name="qsirecon_wf"):
     """
     This workflow organizes the execution of qsiprep, with a sub-workflow for
     each subject.
@@ -70,7 +71,9 @@ def init_qsirecon_wf(subject_list, run_uuid, work_dir, output_dir, recon_input,
                               sloppy=False,
                               omp_nthreads=1,
                               skip_odf_plots=False,
-                              pipeline_source="qsiprep"
+                              pipeline_source="qsiprep",
+                              output_resolution=0.0,
+                              infant_mode=False
                               )
 
 
@@ -98,6 +101,11 @@ def init_qsirecon_wf(subject_list, run_uuid, work_dir, output_dir, recon_input,
             Path to the directory containing subject freesurfer outputs ($SUBJECTS_DIR)
         sloppy : bool
             If True, replace reconstruction options with fast but bad options.
+        infant_mode : bool
+            Use MNI Infant templates
+        output_resolution : float
+            resolution to resample template space outputs. if 0, the native resolution
+            will be used
 
     """
     qsiprep_wf = Workflow(name=name)
@@ -118,7 +126,9 @@ def init_qsirecon_wf(subject_list, run_uuid, work_dir, output_dir, recon_input,
             sloppy=sloppy,
             b0_threshold=b0_threshold,
             freesurfer_input=freesurfer_input,
-            skip_odf_plots=skip_odf_plots
+            skip_odf_plots=skip_odf_plots,
+            infant_mode=infant_mode,
+            output_resolution=output_resolution
             )
 
         single_subject_wf.config['execution']['crashdump_dir'] = (os.path.join(
@@ -132,8 +142,9 @@ def init_qsirecon_wf(subject_list, run_uuid, work_dir, output_dir, recon_input,
 
 
 def init_single_subject_wf(
-        subject_id, name, reportlets_dir, output_dir, freesurfer_input, skip_odf_plots,
-        low_mem, omp_nthreads, recon_input, recon_spec, sloppy, b0_threshold, pipeline_source):
+        subject_id, name, reportlets_dir, output_dir, freesurfer_input,
+        skip_odf_plots, infant_mode, output_resolution, low_mem, omp_nthreads,
+        recon_input, recon_spec, sloppy, b0_threshold, pipeline_source):
     """
     This workflow organizes the reconstruction pipeline for a single subject.
     Reconstruction is performed using a separate workflow for each dwi series.
@@ -162,6 +173,11 @@ def init_single_subject_wf(
             Use bad parameters for reconstruction to make the workflow faster.
         pipeline_source : str
             Which pipeline was used to process the input data
+        infant_mode : bool
+            Use MNI Infant templates
+        output_resolution : float
+            resolution to resample template space outputs. if 0, the native resolution
+            will be used
     """
     if name in ('single_subject_wf', 'single_subject_test_recon_wf'):
         # a fake spec
@@ -221,6 +237,7 @@ def init_single_subject_wf(
             freesurfer_dir=freesurfer_input,
             needs_t1w_transform=needs_t1w_transform,
             pipeline_source="qsiprep",
+            infant_mode=infant_mode,
             name='anat_ingress_wf')
 
         # Connect the anatomical-only inputs. NOTE this is not to the inputnode!
@@ -255,6 +272,7 @@ def init_single_subject_wf(
                 freesurfer_dir=freesurfer_input,
                 pipeline_source="ukb",
                 needs_t1w_transform=needs_t1w_transform,
+                infant_mode=infant_mode,
                 name=wf_name + "_ingressed_ukb_anat_data")
 
         # Create scan-specific anatomical data (mask, atlas configs, odf ROIs for reports)
@@ -263,7 +281,7 @@ def init_single_subject_wf(
             init_dwi_recon_anatomical_workflow(
                 atlas_names=atlas_names,
                 omp_nthreads=omp_nthreads,
-                infant_mode=False,
+                infant_mode=infant_mode,
                 prefer_dwi_mask=False,
                 sloppy=sloppy,
                 needs_t1w_transform=needs_t1w_transform,
@@ -271,6 +289,7 @@ def init_single_subject_wf(
                 freesurfer_dir=freesurfer_input,
                 extras_to_make=spec.get('anatomical', []),
                 name=wf_name + "_dwi_specific_anat_wf",
+                output_resolution=output_resolution,
                 **available_anatomical_data)
 
         # This node holds all the inputs that will go to the recon workflow.
diff --git a/qsiprep/workflows/recon/build_workflow.py b/qsiprep/workflows/recon/build_workflow.py
index e4dc2952..046cd324 100644
--- a/qsiprep/workflows/recon/build_workflow.py
+++ b/qsiprep/workflows/recon/build_workflow.py
@@ -20,11 +20,19 @@
     init_dsi_studio_tractography_wf,
 )
 from .dynamics import init_controllability_wf
+from .utils import init_conform_dwi_wf, init_discard_repeated_samples_wf
+from .steinhardt import init_steinhardt_order_param_wf
+from .scalar_mapping import init_scalar_to_bundle_wf, init_scalar_to_template_wf
+from ...engine import Workflow
+from ...interfaces.interchange import (
+    default_input_set, recon_workflow_input_fields
+)
 from .mrtrix import (
     init_global_tractography_wf,
     init_mrtrix_connectivity_wf,
     init_mrtrix_csd_recon_wf,
     init_mrtrix_tractography_wf,
+
 )
 from .scalar_mapping import init_scalar_to_bundle_wf
 from .steinhardt import init_steinhardt_order_param_wf
@@ -54,6 +62,7 @@ def init_dwi_recon_workflow(workflow_spec, output_dir,
         name='inputnode')
     # Read nodes from workflow spec, make sure we can implement them
     nodes_to_add = []
+    workflow_metadata_nodes = {}
     for node_spec in workflow_spec['nodes']:
         if not node_spec['name']:
             raise Exception("Node has no name [{}]".format(node_spec))
@@ -65,6 +74,14 @@ def init_dwi_recon_workflow(workflow_spec, output_dir,
         if new_node is None:
             raise Exception("Unable to create a node for %s" % node_spec)
         nodes_to_add.append(new_node)
+
+        # Make an identity interface that just has the info of this node
+        workflow_metadata_nodes[node_spec["name"]] = pe.Node(
+            niu.IdentityInterface(fields=["input_metadata"]),
+            name=node_spec["name"] + "_spec")
+        workflow_metadata_nodes[node_spec["name"]].inputs.input_metadata = node_spec
+        nodes_to_add.append(workflow_metadata_nodes[node_spec["name"]])
+
     workflow.add_nodes(nodes_to_add)
     _check_repeats(workflow.list_node_names())
 
@@ -110,7 +127,9 @@ def init_dwi_recon_workflow(workflow_spec, output_dir,
             #             inputnode, node)
             workflow.connect([
                 (inputnode, node,
-                 _as_connections(connect_from_qsiprep, dest_prefix='inputnode.'))])
+                 _as_connections(
+                    connect_from_qsiprep - set(("mapping_metadata",)),
+                    dest_prefix='inputnode.'))])
             # for qp_connection in connect_from_qsiprep:
             #    workflow.connect(inputnode, qp_connection, node, 'inputnode.' + qp_connection)
             _check_repeats(workflow.list_node_names())
@@ -120,21 +139,32 @@ def init_dwi_recon_workflow(workflow_spec, output_dir,
             workflow.connect([
                 (upstream_node, node,
                  _as_connections(
-                    connect_from_upstream, src_prefix='outputnode.', dest_prefix='inputnode.'))])
+                    connect_from_upstream - set(("mapping_metadata",)),
+                    src_prefix='outputnode.',
+                    dest_prefix='inputnode.'))])
             # for upstream_connection in connect_from_upstream:
             #     workflow.connect(upstream_node, "outputnode." + upstream_connection,
             #                      node, 'inputnode.' + upstream_connection)
             _check_repeats(workflow.list_node_names())
 
+            # Send metadata about the upstream node into the downstream node
+            workflow.connect(
+                workflow_metadata_nodes[node_spec['input']],
+                "input_metadata",
+                node,
+                "inputnode.mapping_metadata")
+
     # Fill-in datasinks and reportlet datasinks seen so far
     for node in workflow.list_node_names():
         node_suffix = node.split('.')[-1]
-        if node_suffix.startswith('ds'):
+        if node_suffix.startswith('ds_') or node_suffix.startswith("recon_scalars"):
             base_dir = reportlets_dir if "report" in node_suffix else output_dir
+            workflow.connect(inputnode, 'dwi_file', workflow.get_node(node), 'source_file')
             # LOGGER.info("setting %s base dir to %s", node_suffix, base_dir )
-            workflow.get_node(node).inputs.base_directory = base_dir
-            if node_suffix.startswith('ds_'):
-                workflow.connect(inputnode, 'dwi_file', workflow.get_node(node), 'source_file')
+            if node_suffix.startswith("ds"):
+                workflow.get_node(node).inputs.base_directory = base_dir
+            # if node_suffix.startswith('ds_') or node_suffix.startswith("recon_scalars"):
+            #     workflow.connect(inputnode, 'dwi_file', workflow.get_node(node), 'source_file')
 
     return workflow
 
@@ -143,7 +173,7 @@ def workflow_from_spec(omp_nthreads, available_anatomical_data, node_spec,
                        skip_odf_plots):
     """Build a nipype workflow based on a json file."""
     software = node_spec.get("software", "qsiprep")
-    output_suffix = node_spec.get("output_suffix", "")
+    qsirecon_suffix = node_spec.get("qsirecon_suffix", "")
     node_name = node_spec.get("name", None)
     parameters = node_spec.get("parameters", {})
 
@@ -165,11 +195,9 @@ def workflow_from_spec(omp_nthreads, available_anatomical_data, node_spec,
         "omp_nthreads": omp_nthreads,
         "available_anatomical_data": available_anatomical_data,
         "name": node_name,
-        "output_suffix": output_suffix,
+        "qsirecon_suffix": qsirecon_suffix,
         "params": parameters}
 
-
-
     # DSI Studio operations
     if software == "DSI Studio":
         if node_spec["action"] == "reconstruction":
@@ -213,7 +241,7 @@ def workflow_from_spec(omp_nthreads, available_anatomical_data, node_spec,
         if node_spec["action"] == "pyafq_tractometry":
             return init_pyafq_wf(**kwargs)
 
-    elif software ==  "TORTOISE":
+    elif software == "TORTOISE":
         if node_spec["action"] == "estimate":
             return init_tortoise_estimator_wf(**kwargs)
 
@@ -233,6 +261,8 @@ def workflow_from_spec(omp_nthreads, available_anatomical_data, node_spec,
             return init_steinhardt_order_param_wf(**kwargs)
         if node_spec['action'] == 'bundle_map':
             return init_scalar_to_bundle_wf(**kwargs)
+        if node_spec['action'] == 'template_map':
+            return init_scalar_to_template_wf(**kwargs)
 
     raise Exception("Unknown node %s" % node_spec)
 
diff --git a/qsiprep/workflows/recon/converters.py b/qsiprep/workflows/recon/converters.py
index 98055977..26e634e5 100644
--- a/qsiprep/workflows/recon/converters.py
+++ b/qsiprep/workflows/recon/converters.py
@@ -21,7 +21,7 @@
 
 
 def init_mif_to_fibgz_wf(omp_nthreads, available_anatomical_data,
-                         name="mif_to_fibgz", output_suffix="", params={}):
+                         name="mif_to_fibgz", qsirecon_suffix="", params={}):
     """Converts a MRTrix mif file to DSI Studio fib file.
 
     This workflow uses ``sh2amp`` to sample the FODs on the standard DSI Studio
@@ -52,7 +52,7 @@ def init_mif_to_fibgz_wf(omp_nthreads, available_anatomical_data,
     return workflow
 
 
-def init_fibgz_to_mif_wf(name="fibgz_to_mif", output_suffix="", params={}):
+def init_fibgz_to_mif_wf(name="fibgz_to_mif", qsirecon_suffix="", params={}):
     """Needs Documentation"""
     inputnode = pe.Node(niu.IdentityInterface(fields=recon_workflow_input_fields + ["mif_file"]),
                         name="inputnode")
@@ -68,7 +68,7 @@ def init_fibgz_to_mif_wf(name="fibgz_to_mif", output_suffix="", params={}):
 
 
 def init_qsiprep_to_fsl_wf(omp_nthreads, available_anatomical_data,
-                           name="qsiprep_to_fsl", output_suffix="", params={}):
+                           name="qsiprep_to_fsl", qsirecon_suffix="", params={}):
     """Converts QSIPrep outputs (images, bval, bvec) to fsl standard orientation"""
     inputnode = pe.Node(niu.IdentityInterface(fields=recon_workflow_input_fields),
                         name="inputnode")
@@ -97,28 +97,34 @@ def init_qsiprep_to_fsl_wf(omp_nthreads, available_anatomical_data,
         (convert_mask_to_fsl, outputnode, [('dwi_file', 'mask_file')])
     ])
 
-    if output_suffix:
+    if qsirecon_suffix:
         # Save the output in the outputs directory
         ds_dwi_file = pe.Node(
             ReconDerivativesDataSink(
-                suffix=output_suffix + "_dwi"),
-                name='ds_dwi_' + name,
-                run_without_submitting=True)
+                qsirecon_suffix=qsirecon_suffix,
+                suffix="dwi",
+                extension="bvec"),
+            name='ds_dwi_' + name,
+            run_without_submitting=True)
         ds_bval_file = pe.Node(
             ReconDerivativesDataSink(
-                suffix=output_suffix + "_dwi"),
-                name='ds_bval_' + name,
-                run_without_submitting=True)
+                qsirecon_suffix=qsirecon_suffix,
+                suffix="dwi",
+                extension="bval"),
+            name='ds_bval_' + name,
+            run_without_submitting=True)
         ds_bvec_file = pe.Node(
             ReconDerivativesDataSink(
-                suffix=output_suffix + "_dwi"),
-                name='ds_bvec_' + name,
-                run_without_submitting=True)
+                qsirecon_suffix=qsirecon_suffix,
+                suffix="dwi"),
+            name='ds_bvec_' + name,
+            run_without_submitting=True)
         ds_mask_file = pe.Node(
             ReconDerivativesDataSink(
-                suffix=output_suffix + "_mask"),
-                name='ds_mask_' + name,
-                run_without_submitting=True)
+                qsirecon_suffix=qsirecon_suffix),
+                suffix="mask",
+            name='ds_mask_' + name,
+            run_without_submitting=True)
         workflow.connect([
             (convert_dwi_to_fsl, ds_bval_file, [('bval_file', 'in_file')]),
             (convert_dwi_to_fsl, ds_bvec_file, [('bvec_file', 'in_file')]),
diff --git a/qsiprep/workflows/recon/dipy.py b/qsiprep/workflows/recon/dipy.py
index c9ecf52c..743189af 100644
--- a/qsiprep/workflows/recon/dipy.py
+++ b/qsiprep/workflows/recon/dipy.py
@@ -18,39 +18,41 @@
 )
 from ...interfaces.interchange import recon_workflow_input_fields
 from ...interfaces.recon_scalars import (
+    ReconScalarsDataSink, 
+    DIPYDKIReconScalars, 
     BrainSuite3dSHOREReconScalars,
-    DIPYDKIReconScalars,
     DIPYMAPMRIReconScalars,
 )
 from ...interfaces.reports import CLIReconPeaksReport
 from qsiprep.interfaces.bids import ReconDerivativesDataSink
 
+
 LOGGER = logging.getLogger('nipype.interface')
 
 
-def external_format_datasinks(output_suffix, params, wf):
+def external_format_datasinks(qsirecon_suffix, params, wf):
     """Add datasinks for Dipy Reconstructions in other formats."""
     outputnode = wf.get_node("outputnode")
     if params["write_fibgz"]:
         ds_fibgz = pe.Node(
             ReconDerivativesDataSink(extension='.fib.gz',
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
-            name='ds_{}_fibgz'.format(output_suffix),
+            name='ds_{}_fibgz'.format(qsirecon_suffix),
             run_without_submitting=True)
         wf.connect(outputnode, 'fibgz', ds_fibgz, 'in_file')
     if params["write_mif"]:
         ds_mif = pe.Node(
             ReconDerivativesDataSink(extension='.mif',
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=False),
-            name='ds_{}_mif'.format(output_suffix),
+            name='ds_{}_mif'.format(qsirecon_suffix),
             run_without_submitting=True)
         wf.connect(outputnode, 'fod_sh_mif', ds_mif, 'in_file')
 
 
 def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data, name="dipy_3dshore_recon",
-                                        output_suffix="", params={}):
+                                        qsirecon_suffix="", params={}):
     """Reconstruct EAPs, ODFs, using 3dSHORE (brainsuite-style basis set).
 
     Inputs
@@ -115,14 +117,12 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
         name="outputnode")
 
     workflow = Workflow(name=name)
-    desc = """Dipy Reconstruction
-
-: """
+    desc = "Dipy Reconstruction\n\n: "
     plot_reports = params.pop("plot_reports", True)
     recon_shore = pe.Node(
         BrainSuiteShoreReconstruction(**params), name="recon_shore")
     recon_scalars = pe.Node(
-        BrainSuite3dSHOREReconScalars(workflow_name="name"),
+        BrainSuite3dSHOREReconScalars(qsirecon_suffix="name"),
         name="recon_scalars",
         run_without_submitting=True)
     doing_extrapolation = params.get("extrapolate_scheme") in ("HCP", "ABCD")
@@ -183,13 +183,13 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
             run_without_submitting=True)
         workflow.connect(plot_peaks, 'odf_report', ds_report_odfs, 'in_file')
 
-    if output_suffix:
-        external_format_datasinks(output_suffix, params, workflow)
+    if qsirecon_suffix:
+        external_format_datasinks(qsirecon_suffix, params, workflow)
 
         ds_rtop = pe.Node(
             ReconDerivativesDataSink(extension='.nii.gz',
                                      desc="rtop",
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_bsshore_rtop',
             run_without_submitting=True)
@@ -198,7 +198,7 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
         ds_coeff = pe.Node(
             ReconDerivativesDataSink(extension='.nii.gz',
                                      desc="SHOREcoeff",
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_bsshore_coeff',
             run_without_submitting=True)
@@ -207,7 +207,7 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
         ds_alpha = pe.Node(
             ReconDerivativesDataSink(extension='.nii.gz',
                                      desc="L1alpha",
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_bsshore_alpha',
             run_without_submitting=True)
@@ -216,7 +216,7 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
         ds_r2 = pe.Node(
             ReconDerivativesDataSink(extension='.nii.gz',
                                      desc="r2",
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_bsshore_r2',
             run_without_submitting=True)
@@ -225,7 +225,7 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
         ds_cnr = pe.Node(
             ReconDerivativesDataSink(extension='.nii.gz',
                                      desc="CNR",
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_bsshore_cnr',
             run_without_submitting=True)
@@ -234,7 +234,7 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
         ds_regl = pe.Node(
             ReconDerivativesDataSink(extension='.nii.gz',
                                      desc="regularization",
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_bsshore_regl',
             run_without_submitting=True)
@@ -243,7 +243,7 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
             ds_extrap_dwi = pe.Node(
                 ReconDerivativesDataSink(extension='.nii.gz',
                                          desc="extrapolated",
-                                         suffix=output_suffix,
+                                         qsirecon_suffix=qsirecon_suffix,
                                          compress=True),
                 name='ds_extrap_dwi',
                 run_without_submitting=True)
@@ -251,21 +251,21 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
             ds_extrap_bval = pe.Node(
                 ReconDerivativesDataSink(extension='.bval',
                                          desc="extrapolated",
-                                         suffix=output_suffix),
+                                         qsirecon_suffix=qsirecon_suffix),
                 name='ds_extrap_bval',
                 run_without_submitting=True)
             workflow.connect(outputnode, 'bval_file', ds_extrap_bval, 'in_file')
             ds_extrap_bvec = pe.Node(
                 ReconDerivativesDataSink(extension='.bvec',
                                          desc="extrapolated",
-                                         suffix=output_suffix),
+                                         qsirecon_suffix=qsirecon_suffix),
                 name='ds_extrap_bvec',
                 run_without_submitting=True)
             workflow.connect(outputnode, 'bvec_file', ds_extrap_bvec, 'in_file')
             ds_extrap_b = pe.Node(
                 ReconDerivativesDataSink(extension='.b',
                                          desc="extrapolated",
-                                         suffix=output_suffix),
+                                         qsirecon_suffix=qsirecon_suffix),
                 name='ds_extrap_b',
                 run_without_submitting=True)
             workflow.connect(outputnode, 'b_file', ds_extrap_b, 'in_file')
@@ -274,7 +274,7 @@ def init_dipy_brainsuite_shore_recon_wf(omp_nthreads, available_anatomical_data,
 
 
 def init_dipy_mapmri_recon_wf(omp_nthreads, available_anatomical_data, name="dipy_mapmri_recon",
-                              output_suffix="", params={}):
+                              qsirecon_suffix="", params={}):
     """Reconstruct EAPs, ODFs, using 3dSHORE (brainsuite-style basis set).
 
     Inputs
@@ -358,13 +358,11 @@ def init_dipy_mapmri_recon_wf(omp_nthreads, available_anatomical_data, name="dip
         name="outputnode")
 
     workflow = Workflow(name=name)
-    desc = """Dipy Reconstruction
-
-: """
+    desc = "Dipy Reconstruction\n\n: "
     plot_reports = params.pop("plot_reports", True)
     recon_map = pe.Node(MAPMRIReconstruction(**params), name="recon_map")
     recon_scalars = pe.Node(
-        DIPYMAPMRIReconScalars(workflow_name=name),
+        DIPYMAPMRIReconScalars(qsirecon_suffix=name),
         name="recon_scalars",
         run_without_submitting=True)
     workflow.connect([
@@ -424,23 +422,24 @@ def init_dipy_mapmri_recon_wf(omp_nthreads, available_anatomical_data, name="dip
             run_without_submitting=True)
         workflow.connect(plot_peaks, 'odf_report', ds_report_odfs, 'in_file')
 
-    if output_suffix:
-        external_format_datasinks(output_suffix, params, workflow)
-        connections = []
-        for scalar_name in ['rtop', 'rtap', 'rtpp', 'qiv', 'msd', 'lapnorm']:
-            connections += [(outputnode,
-                             pe.Node(
-                                 ReconDerivativesDataSink(desc=scalar_name,
-                                                          suffix=output_suffix),
-                                 name='ds_%s_%s' % (name, scalar_name)),
-                             [(scalar_name, 'in_file')])]
-        workflow.connect(connections)
+    if qsirecon_suffix:
+        external_format_datasinks(qsirecon_suffix, params, workflow)
+        ds_recon_scalars = pe.Node(
+            ReconScalarsDataSink(),
+            name="ds_recon_scalars",
+            run_without_submitting=True)
+        workflow.connect(
+            recon_scalars,
+            "scalar_info",
+            ds_recon_scalars,
+            "recon_scalars")
+
     workflow.__desc__ = desc
     return workflow
 
 
 def init_dipy_dki_recon_wf(omp_nthreads, available_anatomical_data, name="dipy_dki_recon",
-                              output_suffix="", params={}):
+                              qsirecon_suffix="", params={}):
     """Fit DKI
 
     Inputs
@@ -480,13 +479,11 @@ def init_dipy_dki_recon_wf(omp_nthreads, available_anatomical_data, name="dipy_d
                     'colorFA', 'kfa', 'mk', 'ak', 'rk',
                     'mkt', 'recon_scalars']),
         name="outputnode")
-    recon_scalars = pe.Node(DIPYDKIReconScalars(workflow_name=name),
+    recon_scalars = pe.Node(DIPYDKIReconScalars(qsirecon_suffix=qsirecon_suffix),
                             run_without_submitting=True,
                             name="recon_scalars")
     workflow = Workflow(name=name)
-    desc = """Dipy Reconstruction
-
-: """
+    desc = "Dipy Reconstruction\n\n: "
     plot_reports = params.pop("plot_reports", True)
     recon_dki = pe.Node(KurtosisReconstruction(**params), name="recon_dki")
 
@@ -541,19 +538,17 @@ def init_dipy_dki_recon_wf(omp_nthreads, available_anatomical_data, name="dipy_d
                                      ('odf_amplitudes', 'odf_file')]),
             (plot_peaks, ds_report_peaks, [('peak_report', 'in_file')])])
 
-    if output_suffix:
-        external_format_datasinks(output_suffix, params, workflow)
-        connections = []
-        for scalar_name in ['tensor', 'fa', 'md', 'rd', 'ad',
-                    'colorFA', 'kfa', 'mk', 'ak', 'rk',
-                    'mkt']:
-            connections += [(outputnode,
-                             pe.Node(
-                                 ReconDerivativesDataSink(desc=scalar_name,
-                                                          suffix=output_suffix),
-                                 name='ds_%s_%s' % (name, scalar_name)),
-                             [(scalar_name, 'in_file')])]
-        workflow.connect(connections)
+    if qsirecon_suffix:
+        external_format_datasinks(qsirecon_suffix, params, workflow)
+        ds_recon_scalars = pe.Node(
+            ReconScalarsDataSink(),
+            name="ds_recon_scalars",
+            run_without_submitting=True)
+        workflow.connect(
+            recon_scalars,
+            "scalar_info",
+            ds_recon_scalars,
+            "recon_scalars")
 
     workflow.__desc__ = desc
     return workflow
\ No newline at end of file
diff --git a/qsiprep/workflows/recon/dsi_studio.py b/qsiprep/workflows/recon/dsi_studio.py
index a5c71d8e..15f16421 100644
--- a/qsiprep/workflows/recon/dsi_studio.py
+++ b/qsiprep/workflows/recon/dsi_studio.py
@@ -14,6 +14,8 @@
 
 from ...engine import Workflow
 from ...interfaces.bids import ReconDerivativesDataSink
+from ...interfaces.recon_scalars import DSIStudioReconScalars, ReconScalarsDataSink
+
 from ...interfaces.converters import DSIStudioTrkToTck
 from ...interfaces.interchange import recon_workflow_input_fields
 from ...interfaces.recon_scalars import DSIStudioReconScalars
@@ -35,7 +37,7 @@
 
 
 def init_dsi_studio_recon_wf(omp_nthreads, available_anatomical_data, name="dsi_studio_recon",
-                             output_suffix="", params={}):
+                             qsirecon_suffix="", params={}):
     """Reconstructs diffusion data using DSI Studio.
 
     This workflow creates a ``.src.gz`` file from the input dwi, bvals and bvecs,
@@ -119,12 +121,12 @@ def init_dsi_studio_recon_wf(omp_nthreads, available_anatomical_data, name="dsi_
             workflow.connect(plot_peaks, 'odf_report', ds_report_odfs, 'in_file')
 
 
-    if output_suffix:
+    if qsirecon_suffix:
         # Save the output in the outputs directory
         ds_gqi_fibgz = pe.Node(
             ReconDerivativesDataSink(
                 extension='.fib.gz',
-                suffix=output_suffix,
+                qsirecon_suffix=qsirecon_suffix,
                 compress=True),
             name='ds_gqi_fibgz',
             run_without_submitting=True)
@@ -134,7 +136,7 @@ def init_dsi_studio_recon_wf(omp_nthreads, available_anatomical_data, name="dsi_
 
 
 def init_dsi_studio_tractography_wf(omp_nthreads, available_anatomical_data, name="dsi_studio_tractography",
-                                    params={}, output_suffix=""):
+                                    params={}, qsirecon_suffix=""):
     """Calculate streamline-based connectivity matrices using DSI Studio.
 
     DSI Studio has a deterministic tractography algorithm that can be used to
@@ -201,10 +203,9 @@ def init_dsi_studio_tractography_wf(omp_nthreads, available_anatomical_data, nam
     outputnode.inputs.recon_scalars = []
     plot_reports = params.pop("plot_reports", True)
     workflow = Workflow(name=name)
-    desc = """DSI Studio Tractography
-
-: Tractography was run in DSI Studio (version %s) using a deterministic algorithm
-[@yeh2013deterministic]. """ %  DSI_STUDIO_VERSION
+    desc = "DSI Studio Tractography\n\n: Tractography was run in DSI Studio " \
+        "(version %s) using a deterministic algorithm " \
+        "[@yeh2013deterministic]. " %  DSI_STUDIO_VERSION
     tracking = pe.Node(
         DSIStudioTracking(
             num_threads=omp_nthreads,
@@ -216,17 +217,21 @@ def init_dsi_studio_tractography_wf(omp_nthreads, available_anatomical_data, nam
         (tracking, outputnode, [('output_trk', 'trk_file')]),
         (inputnode, outputnode, [('fibgz', 'fibgz')])
     ])
-    if output_suffix:
+    if qsirecon_suffix:
         # Save the output in the outputs directory
-        ds_tracking = pe.Node(ReconDerivativesDataSink(suffix=output_suffix),
-                              name='ds_' + name,
-                              run_without_submitting=True)
+        ds_tracking = pe.Node(
+            ReconDerivativesDataSink(
+                model="QAthresh",
+                suffix="streamlines",
+                qsirecon_suffix=qsirecon_suffix),
+            name='ds_' + name,
+            run_without_submitting=True)
         workflow.connect(tracking, 'output_trk', ds_tracking, 'in_file')
     return workflow
 
 
 def init_dsi_studio_autotrack_wf(omp_nthreads, available_anatomical_data,
-                                 params={}, output_suffix="", name="dsi_studio_autotrack_wf"):
+                                 params={}, qsirecon_suffix="", name="dsi_studio_autotrack_wf"):
     """Run DSI Studio's AutoTrack method to produce bundles and bundle stats.
 
     Inputs
@@ -278,10 +283,8 @@ def init_dsi_studio_autotrack_wf(omp_nthreads, available_anatomical_data,
     outputnode = pe.Node(niu.IdentityInterface(fields=['tck_files', "bundle_names", "recon_scalars"]),
                          name="outputnode")
     outputnode.inputs.recon_scalars = []
-    desc = """DSI Studio Automatic Tractography
-
-: Automatic Tractography was run in DSI Studio (version %s) and
-bundle shape statistics were calculated [@autotrack]. """ %  DSI_STUDIO_VERSION
+    desc = "DSI Studio Automatic Tractography\n\n: Automatic Tractography was run in " \
+        "DSI Studio (version %s) and bundle shape statistics were calculated [@autotrack]. " %  DSI_STUDIO_VERSION
     plot_reports = params.pop("plot_reports", True)
     bundle_names = _get_dsi_studio_bundles(params.get("track_id", ""))
     bundle_desc = "AutoTrack attempted to reconstruct the following bundles:\n  * " \
@@ -311,21 +314,24 @@ def init_dsi_studio_autotrack_wf(omp_nthreads, available_anatomical_data,
     # Save tck files of the bundles into the outputs
     ds_tckfiles = pe.MapNode(
         ReconDerivativesDataSink(
-            suffix=output_suffix),
+            suffix="streamlines",
+            qsirecon_suffix=qsirecon_suffix),
         iterfield=["in_file", "bundle"],
         name="ds_tckfiles")
 
     # Save the bundle csv
     ds_bundle_csv = pe.Node(
         ReconDerivativesDataSink(
-            suffix=output_suffix),
+            suffix="bundlestats",
+            qsirecon_suffix=qsirecon_suffix),
         name="ds_bundle_csv",
         run_without_submitting=True)
 
     # Save the mapping file
     ds_mapping = pe.Node(
         ReconDerivativesDataSink(
-            suffix=output_suffix),
+            suffix="mapping",
+            qsirecon_suffix=qsirecon_suffix),
         name="ds_mapping",
         run_without_submitting=True)
 
@@ -352,7 +358,7 @@ def init_dsi_studio_autotrack_wf(omp_nthreads, available_anatomical_data,
 
 
 def init_dsi_studio_connectivity_wf(omp_nthreads, available_anatomical_data, name="dsi_studio_connectivity",
-                                    params={}, output_suffix=""):
+                                    params={}, qsirecon_suffix=""):
     """Calculate streamline-based connectivity matrices using DSI Studio.
 
     DSI Studio has a deterministic tractography algorithm that can be used to
@@ -441,17 +447,20 @@ def init_dsi_studio_connectivity_wf(omp_nthreads, available_anatomical_data, nam
                 ('connectivity_matfile', 'connectivity_matfile')]),
             (plot_connectivity, ds_report_connectivity, [('out_report', 'in_file')]),
         ])
-    if output_suffix:
+    if qsirecon_suffix:
         # Save the output in the outputs directory
-        ds_connectivity = pe.Node(ReconDerivativesDataSink(suffix=output_suffix),
-                                  name='ds_' + name,
-                                  run_without_submitting=True)
+        ds_connectivity = pe.Node(
+            ReconDerivativesDataSink(
+                qsirecon_suffix=qsirecon_suffix,
+                suffix="connectivity"),
+            name='ds_' + name,
+            run_without_submitting=True)
         workflow.connect(calc_connectivity, 'connectivity_matfile', ds_connectivity, 'in_file')
     return workflow
 
 
 def init_dsi_studio_export_wf(omp_nthreads, available_anatomical_data, name="dsi_studio_export",
-                              params={}, output_suffix=""):
+                              params={}, qsirecon_suffix=""):
     """Export scalar maps from a DSI Studio fib file into NIfTI files with correct headers.
 
     This workflow exports gfa, fa0, fa1, fa2 and iso.
@@ -489,7 +498,7 @@ def init_dsi_studio_export_wf(omp_nthreads, available_anatomical_data, name="dsi
     workflow = pe.Workflow(name=name)
     export = pe.Node(DSIStudioExport(to_export=",".join(scalar_names)), name='export')
     recon_scalars = pe.Node(
-        DSIStudioReconScalars(workflow_name=name),
+        DSIStudioReconScalars(qsirecon_suffix=qsirecon_suffix),
         name="recon_scalars",
         n_procs=1)
     fixhdr_nodes = {}
@@ -501,15 +510,19 @@ def init_dsi_studio_export_wf(omp_nthreads, available_anatomical_data, name="dsi
                                                                 'correct_header_nifti')]),
                        (fixhdr_nodes[scalar_name], outputnode, [('out_file', output_name)]),
                        (fixhdr_nodes[scalar_name], recon_scalars, [("out_file", output_name)])]
-        if output_suffix:
-            connections += [(fixhdr_nodes[scalar_name],
-                             pe.Node(
-                                 ReconDerivativesDataSink(desc=scalar_name,
-                                                          suffix=output_suffix),
-                                 name='ds_%s_%s' % (name, scalar_name)),
-                             [('out_file', 'in_file')])]
         workflow.connect(connections)
 
+    if qsirecon_suffix:
+        ds_recon_scalars = pe.Node(
+            ReconScalarsDataSink(),
+            name="ds_recon_scalars",
+            run_without_submitting=True)
+        workflow.connect(
+            recon_scalars,
+            "scalar_info",
+            ds_recon_scalars,
+            "recon_scalars")
+
     workflow.connect([
         (inputnode, export, [('fibgz', 'input_file')]),
         (recon_scalars, outputnode, [("scalar_info", "recon_scalars")])])
diff --git a/qsiprep/workflows/recon/dynamics.py b/qsiprep/workflows/recon/dynamics.py
index 14bb7706..5718e140 100644
--- a/qsiprep/workflows/recon/dynamics.py
+++ b/qsiprep/workflows/recon/dynamics.py
@@ -17,7 +17,7 @@
 LOGGER = logging.getLogger('nipype.workflow')
 
 
-def init_controllability_wf(name="controllability", output_suffix="", params={}, **kwargs):
+def init_controllability_wf(name="controllability", qsirecon_suffix="", params={}, **kwargs):
     """Calculates network controllability from connectivity matrices.
 
     Calculates modal and average controllability using the method of Gu et al. 2015.
@@ -49,10 +49,13 @@ def init_controllability_wf(name="controllability", output_suffix="", params={},
         (inputnode, calc_control, [('matfile', 'matfile')]),
         (calc_control, outputnode, [('controllability', 'matfile')])
     ])
-    if output_suffix:
+    if qsirecon_suffix:
         # Save the output in the outputs directory
-        ds_control = pe.Node(ReconDerivativesDataSink(suffix=output_suffix),
-                             name='ds_' + name,
-                             run_without_submitting=True)
+        ds_control = pe.Node(
+            ReconDerivativesDataSink(
+                qsirecon_suffix=qsirecon_suffix,
+                suffix="control"),
+            name='ds_' + name,
+            run_without_submitting=True)
         workflow.connect(calc_control, 'controllability', ds_control, 'in_file')
     return workflow
diff --git a/qsiprep/workflows/recon/mrtrix.py b/qsiprep/workflows/recon/mrtrix.py
index 18d8daa1..de85b606 100644
--- a/qsiprep/workflows/recon/mrtrix.py
+++ b/qsiprep/workflows/recon/mrtrix.py
@@ -46,7 +46,7 @@
 
 
 def init_mrtrix_csd_recon_wf(omp_nthreads, available_anatomical_data, name="mrtrix_recon",
-                             output_suffix="", params={}):
+                             qsirecon_suffix="", params={}):
     """Create FOD images for WM, GM and CSF.
 
     This workflow uses mrtrix tools to run csd on multishell data. At the end,
@@ -103,9 +103,7 @@ def init_mrtrix_csd_recon_wf(omp_nthreads, available_anatomical_data, name="mrtr
     workflow = Workflow(name=name)
     outputnode.inputs.recon_scalars = []
     plot_reports = params.pop("plot_reports", True)
-    desc = """MRtrix3 Reconstruction
-
-: """
+    desc = """MRtrix3 Reconstruction\n\n: """
 
     # Response estimation
     response = params.get('response', {})
@@ -122,11 +120,12 @@ def init_mrtrix_csd_recon_wf(omp_nthreads, available_anatomical_data, name="mrtr
         desc += 'Multi-tissue '
     LOGGER.info("Response configuration: %s", response)
 
-    desc += """\
-fiber response functions were estimated using the {} algorithm.
-FODs were estimated via constrained spherical deconvolution
-(CSD, @originalcsd, @tournier2008csd) \
-""".format(response_algorithm)
+    desc += "\n".join([
+        "fiber response functions were estimated using the {} algorithm. ",
+        "FODs were estimated via constrained spherical deconvolution ",
+        "(CSD, @originalcsd, @tournier2008csd) "
+    ]).format(response_algorithm)
+
     if response_algorithm == 'msmt_5tt':
         desc += 'using a T1w-based segmentation {}.'.format(CITATIONS[response_algorithm])
     else:
@@ -170,9 +169,9 @@ def init_mrtrix_csd_recon_wf(omp_nthreads, available_anatomical_data, name="mrtr
             SS3TEstimateFOD(**fod),
             name='estimate_fod',
             n_procs=omp_nthreads)
-        desc += """ \
-A single-shell-optimized multi-tissue CSD was performed using MRtrix3Tissue
-(https://3Tissue.github.io), a fork of MRtrix3 (@mrtrix3)"""
+        desc += "\n".join([
+            "A single-shell-optimized multi-tissue CSD was performed using MRtrix3Tissue",
+            "(https://3Tissue.github.io), a fork of MRtrix3 (@mrtrix3)"])
 
     workflow.connect([
         (estimate_response, estimate_fod, [('wm_file', 'wm_txt'),
@@ -249,20 +248,23 @@ def init_mrtrix_csd_recon_wf(omp_nthreads, available_anatomical_data, name="mrtr
             else (intensity_norm, "wm_normed_odf")
         workflow.connect(fod_source, fod_key, plot_peaks, "mif_file")
 
-    if output_suffix:
-        normed = '' if not run_mtnormalize else 'mtnormed'
+    if qsirecon_suffix:
         ds_wm_odf = pe.Node(
             ReconDerivativesDataSink(extension='.mif.gz',
-                                     desc="wmFOD" + normed,
-                                     suffix=output_suffix,
+                                     model=fod_algorithm,
+                                     mfp="FOD",
+                                     label="WM",
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_wm_odf',
             run_without_submitting=True)
         workflow.connect(outputnode, 'wm_odf', ds_wm_odf, 'in_file')
         ds_wm_txt = pe.Node(
             ReconDerivativesDataSink(extension='.txt',
-                                     desc="wmFOD",
-                                     suffix=output_suffix),
+                                     model=fod_algorithm,
+                                     mfp="FOD",
+                                     label="WM",
+                                     qsirecon_suffix=qsirecon_suffix),
             name='ds_wm_txt',
             run_without_submitting=True)
         workflow.connect(outputnode, 'wm_txt', ds_wm_txt, 'in_file')
@@ -271,32 +273,40 @@ def init_mrtrix_csd_recon_wf(omp_nthreads, available_anatomical_data, name="mrtr
         if using_multitissue:
             ds_gm_odf = pe.Node(
                 ReconDerivativesDataSink(extension='.mif.gz',
-                                         desc="gmFOD" + normed,
-                                         suffix=output_suffix,
+                                         model=fod_algorithm,
+                                         mfp="FOD",
+                                         label="GM",
+                                         qsirecon_suffix=qsirecon_suffix,
                                          compress=True),
                 name='ds_gm_odf',
                 run_without_submitting=True)
             workflow.connect(outputnode, 'gm_odf', ds_gm_odf, 'in_file')
             ds_gm_txt = pe.Node(
                 ReconDerivativesDataSink(extension='.txt',
-                                         desc="gmFOD",
-                                         suffix=output_suffix),
+                                         model=fod_algorithm,
+                                         mfp="FOD",
+                                         label="GM",
+                                         qsirecon_suffix=qsirecon_suffix),
                 name='ds_gm_txt',
                 run_without_submitting=True)
             workflow.connect(outputnode, 'gm_txt', ds_gm_txt, 'in_file')
 
             ds_csf_odf = pe.Node(
                 ReconDerivativesDataSink(extension='.mif.gz',
-                                         desc="csfFOD" + normed,
-                                         suffix=output_suffix,
+                                         model=fod_algorithm,
+                                         mfp="FOD",
+                                         label="CSF",
+                                         qsirecon_suffix=qsirecon_suffix,
                                          compress=True),
                 name='ds_csf_odf',
                 run_without_submitting=True)
             workflow.connect(outputnode, 'csf_odf', ds_csf_odf, 'in_file')
             ds_csf_txt = pe.Node(
                 ReconDerivativesDataSink(extension='.txt',
-                                         desc="csfFOD",
-                                         suffix=output_suffix),
+                                         model=fod_algorithm,
+                                         mfp="FOD",
+                                         label="CSF",
+                                         qsirecon_suffix=qsirecon_suffix),
                 name='ds_csf_txt',
                 run_without_submitting=True)
             workflow.connect(outputnode, 'csf_txt', ds_csf_txt, 'in_file')
@@ -304,16 +314,18 @@ def init_mrtrix_csd_recon_wf(omp_nthreads, available_anatomical_data, name="mrtr
             if run_mtnormalize:
                 ds_mt_norm = pe.Node(
                     ReconDerivativesDataSink(extension='.mif.gz',
-                                             desc="mtnorm",
-                                             suffix=output_suffix,
+                                             model="mtnorm",
+                                             mfp="norm",
+                                             qsirecon_suffix=qsirecon_suffix,
                                              compress=True),
                     name='ds_mt_norm',
                     run_without_submitting=True)
                 workflow.connect(intensity_norm, 'norm_image', ds_mt_norm, 'in_file')
                 ds_inlier_mask = pe.Node(
                     ReconDerivativesDataSink(extension='.mif.gz',
-                                             desc="mtinliermask",
-                                             suffix=output_suffix,
+                                             model="mtnorm",
+                                             mfp="inliermask",
+                                             qsirecon_suffix=qsirecon_suffix,
                                              compress=True),
                     name='ds_inlier_mask',
                     run_without_submitting=True)
@@ -324,7 +336,7 @@ def init_mrtrix_csd_recon_wf(omp_nthreads, available_anatomical_data, name="mrtr
 
 
 def init_global_tractography_wf(omp_nthreads, available_anatomical_data, name="mrtrix_recon",
-                                output_suffix="", params={}):
+                                qsirecon_suffix="", params={}):
     """Run multi-shell, multi-tissue global tractography
 
     This workflow uses mrtrix tools to run csd on multishell data.
@@ -385,11 +397,11 @@ def init_global_tractography_wf(omp_nthreads, available_anatomical_data, name="m
                                   ("wm_odf", "fod_sh_mif")])
         ])
 
-    if output_suffix:
+    if qsirecon_suffix:
         ds_globalwm_odf = pe.Node(
             ReconDerivativesDataSink(extension='.mif.gz',
                                      desc="globalwmFOD",
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_globalwm_odf',
             run_without_submitting=True)
@@ -398,7 +410,7 @@ def init_global_tractography_wf(omp_nthreads, available_anatomical_data, name="m
         ds_isotropic_fraction = pe.Node(
             ReconDerivativesDataSink(extension='.mif.gz',
                                      desc="ISOfraction",
-                                     suffix=output_suffix),
+                                     qsirecon_suffix=qsirecon_suffix),
             name='ds_isotropic_fraction',
             run_without_submitting=True)
         workflow.connect(outputnode, 'isotropic_fraction', ds_isotropic_fraction, 'in_file')
@@ -406,7 +418,7 @@ def init_global_tractography_wf(omp_nthreads, available_anatomical_data, name="m
         ds_tck_file = pe.Node(
             ReconDerivativesDataSink(extension='.tck.gz',
                                      desc="global",
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_tck_file',
             run_without_submitting=True)
@@ -415,7 +427,7 @@ def init_global_tractography_wf(omp_nthreads, available_anatomical_data, name="m
         ds_residual_energy = pe.Node(
             ReconDerivativesDataSink(extension='.tck.gz',
                                      desc="residualEnergy",
-                                     suffix=output_suffix,
+                                     qsirecon_suffix=qsirecon_suffix,
                                      compress=True),
             name='ds_residual_energy',
             run_without_submitting=True)
@@ -425,7 +437,7 @@ def init_global_tractography_wf(omp_nthreads, available_anatomical_data, name="m
 
 
 def init_mrtrix_tractography_wf(omp_nthreads, available_anatomical_data, name="mrtrix_tracking",
-                                output_suffix="", params={}):
+                                qsirecon_suffix="", params={}):
     """Run tractography
 
     This workflow uses mrtrix tools to run csd on multishell data.
@@ -498,22 +510,24 @@ def init_mrtrix_tractography_wf(omp_nthreads, available_anatomical_data, name="m
                 ('out_mu', 'mu'),
                 ('out_weights', 'sift_weights')])
         ])
-        if output_suffix:
+        if qsirecon_suffix:
             ds_sift_weights = pe.Node(
                 ReconDerivativesDataSink(extension='.csv',
-                                         desc="siftweights",
-                                         suffix=output_suffix),
+                                         model="sift2",
+                                         suffix="streamlineweights",
+                                         qsirecon_suffix=qsirecon_suffix),
                 name='ds_sift_weights',
                 run_without_submitting=True)
             workflow.connect(outputnode, 'sift_weights', ds_sift_weights, 'in_file')
         if use_5tt:
             workflow.connect(inputnode, connect_5tt, tck_sift2, "act_file")
 
-    if output_suffix:
+    if qsirecon_suffix:
         ds_tck_file = pe.Node(
             ReconDerivativesDataSink(extension='.tck',
-                                     desc="tracks",
-                                     suffix=output_suffix),
+                                     model=tracking_params["algorithm"],
+                                     suffix="streamlines",
+                                     qsirecon_suffix=qsirecon_suffix),
             name='ds_tck_file',
             run_without_submitting=True)
         workflow.connect(outputnode, 'tck_file', ds_tck_file, 'in_file')
@@ -522,7 +536,7 @@ def init_mrtrix_tractography_wf(omp_nthreads, available_anatomical_data, name="m
 
 
 def init_mrtrix_connectivity_wf(omp_nthreads, available_anatomical_data, name="mrtrix_connectiity",
-                                params={}, output_suffix=""):
+                                params={}, qsirecon_suffix=""):
     """Runs ``tck2connectome`` on a ``tck`` file.
 
     Inputs
@@ -576,15 +590,18 @@ def init_mrtrix_connectivity_wf(omp_nthreads, available_anatomical_data, name="m
                 ('connectivity_matfile', 'connectivity_matfile')]),
             (plot_connectivity, ds_report_connectivity, [('out_report', 'in_file')])])
 
-    if output_suffix:
+    if qsirecon_suffix:
         # Save the output in the outputs directory
-        ds_connectivity = pe.Node(ReconDerivativesDataSink(suffix=output_suffix),
-                                  name='ds_' + name,
-                                  run_without_submitting=True)
+        ds_connectivity = pe.Node(
+            ReconDerivativesDataSink(
+                suffix="connectivity",
+                qsirecon_suffix=qsirecon_suffix),
+            name='ds_' + name,
+            run_without_submitting=True)
         ds_exemplars = pe.Node(
             ReconDerivativesDataSink(
-                suffix=output_suffix,
-                desc="exemplarbundles"),
+                qsirecon_suffix=qsirecon_suffix,
+                suffix="exemplarbundles"),
             name='ds_exemplars',
             run_without_submitting=True)
         workflow.connect([
diff --git a/qsiprep/workflows/recon/pyafq.py b/qsiprep/workflows/recon/pyafq.py
index 4275a72a..88e4bf17 100644
--- a/qsiprep/workflows/recon/pyafq.py
+++ b/qsiprep/workflows/recon/pyafq.py
@@ -49,7 +49,7 @@ def _parse_qsiprep_params_dict(params_dict):
 
 
 def init_pyafq_wf(omp_nthreads, available_anatomical_data,
-                  name="afq", output_suffix="", params={}):
+                  name="afq", qsirecon_suffix="", params={}):
     """Run PyAFQ on some qsiprep outputs
 
     Inputs
@@ -86,7 +86,7 @@ def init_pyafq_wf(omp_nthreads, available_anatomical_data,
             ('t1_2_mni_reverse_transform', 'itk_file')]),
         (run_afq, outputnode, [('afq_dir', 'afq_dir')])
     ])
-    if output_suffix:
+    if qsirecon_suffix:
         # Save the output in the outputs directory
         ds_afq = pe.Node(
             ReconDerivativesDataSink(),
diff --git a/qsiprep/workflows/recon/scalar_mapping.py b/qsiprep/workflows/recon/scalar_mapping.py
index 60a84179..a2d27b49 100644
--- a/qsiprep/workflows/recon/scalar_mapping.py
+++ b/qsiprep/workflows/recon/scalar_mapping.py
@@ -1,9 +1,8 @@
 """
-Converting between file formats
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Summarize and Transform recon outputs
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
 
-.. autofunction:: init_mif_to_fibgz_wf
-.. autofunction:: init_fibgz_to_mif_wf
 
 """
 import logging
@@ -14,17 +13,17 @@
 from ...engine import Workflow
 from ...interfaces.bids import ReconDerivativesDataSink
 from ...interfaces.interchange import recon_workflow_input_fields
-from ...interfaces.scalar_mapping import BundleMapper
+from ...engine import Workflow
+from ...interfaces.recon_scalars import ReconScalarsDataSink, ReconScalarsTableSplitterDataSink
+from ...interfaces.scalar_mapping import BundleMapper, TemplateMapper
 
 LOGGER = logging.getLogger('nipype.workflow')
 
 
 def init_scalar_to_bundle_wf(omp_nthreads, available_anatomical_data,
-                             name="scalar_to_bundle", output_suffix="", params={}):
+                             name="scalar_to_bundle", qsirecon_suffix="", params={}):
     """Map scalar images to bundles
 
-
-
     Inputs
         tck_files
             MRtrix3 format tck files for each bundle
@@ -36,7 +35,7 @@ def init_scalar_to_bundle_wf(omp_nthreads, available_anatomical_data,
     Outputs
 
         bundle_summaries
-            summary statistics in tsv formar
+            summary statistics in tsv format
 
     """
     inputnode = pe.Node(
@@ -50,31 +49,36 @@ def init_scalar_to_bundle_wf(omp_nthreads, available_anatomical_data,
     bundle_mapper = pe.Node(
         BundleMapper(**params),
         name="bundle_mapper")
+    ds_bundle_mapper = pe.Node(
+        ReconScalarsTableSplitterDataSink(
+            suffix="scalarstats"),
+        name="ds_bundle_mapper",
+        run_without_submitting=True)
+    ds_tdi_summary = pe.Node(
+        ReconScalarsTableSplitterDataSink(
+            suffix="tdistats"),
+        name='ds_tdi_summary',
+        run_without_submitting=True)
     workflow.connect([
         (inputnode, bundle_mapper, [
             ("collected_scalars", "recon_scalars"),
             ("tck_files", "tck_files"),
             ("dwi_ref", "dwiref_image"),
+            ("mapping_metadata", "mapping_metadata"),
             ("bundle_names", "bundle_names")]),
+        (bundle_mapper, ds_bundle_mapper, [
+            ("bundle_summary", "summary_tsv")]),
         (bundle_mapper, outputnode, [
-            ("bundle_summary", "bundle_summary")])
-    ])
-    if output_suffix:
-
-        ds_bundle_summaries = pe.Node(
-            ReconDerivativesDataSink(desc="bundlemap",
-                                     suffix=output_suffix),
-            name='ds_bundle_summaries',
-            run_without_submitting=True)
-        workflow.connect([
-            (bundle_mapper, ds_bundle_summaries, [("bundle_summary", "in_file")])
+            ("bundle_summary", "bundle_summary")]),
+        (bundle_mapper, ds_tdi_summary, [
+            ("tdi_stats", "summary_tsv")])
         ])
 
     return workflow
 
 
 def init_scalar_to_atlas_wf(omp_nthreads, available_anatomical_data,
-                            name="scalar_to_template", output_suffix="", params={}):
+                            name="scalar_to_template", qsirecon_suffix="", params={}):
     """Map scalar images to atlas regions
 
     Inputs
@@ -93,10 +97,10 @@ def init_scalar_to_atlas_wf(omp_nthreads, available_anatomical_data,
     inputnode = pe.Node(
         niu.IdentityInterface(
             fields=recon_workflow_input_fields +
-            ["tck_files", "bundle_names", "recon_scalars"]),
+            ["recon_scalars", "collected_scalars"]),
         name="inputnode")
     outputnode = pe.Node(
-        niu.IdentityInterface(fields=['bundle_summaries']), name="outputnode")
+        niu.IdentityInterface(fields=['atlas_summaries']), name="outputnode")
     workflow = Workflow(name=name)
     bundle_mapper = pe.Node(
         BundleMapper(**params),
@@ -107,11 +111,11 @@ def init_scalar_to_atlas_wf(omp_nthreads, available_anatomical_data,
             ("tck_files", "tck_files"),
             ("dwi_ref", "dwiref_image")])
     ])
-    if output_suffix:
+    if qsirecon_suffix:
 
         ds_bundle_summaries = pe.Node(
             ReconDerivativesDataSink(desc="bundlemap",
-                                     suffix=output_suffix),
+                                     qsirecon_suffix=qsirecon_suffix),
             name='ds_bundle_summaries',
             run_without_submitting=True)
         workflow.connect([
@@ -121,14 +125,56 @@ def init_scalar_to_atlas_wf(omp_nthreads, available_anatomical_data,
 
 
 def init_scalar_to_template_wf(omp_nthreads, available_anatomical_data,
-                               name="scalar_to_template", output_suffix="", params={}):
-    """Maps scalar data to a volumetric template"""
+                               name="scalar_to_template", qsirecon_suffix="", params={}):
+    """Maps scalar data to a volumetric template
+
+
+    Inputs
+        recon_scalars
+            List of dictionaries containing scalar info
+
+    Outputs
+
+        template_scalars
+            List of transformed files
+
+    """
+    inputnode = pe.Node(
+        niu.IdentityInterface(
+            fields=recon_workflow_input_fields +
+            ["recon_scalars", "collected_scalars"]),
+        name="inputnode")
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=['template_scalars', 'template_scalar_sidecars']),
+        name="outputnode")
+    workflow = Workflow(name=name)
+    template_mapper = pe.Node(
+        TemplateMapper(**params),
+        name="template_mapper")
+    # Datasink will always be used, and the qsirecon-suffix is determined when
+    # the scalars are originally calculated
+    ds_template_scalars = pe.Node(
+        ReconScalarsDataSink(),
+        name='ds_template_scalars',
+        run_without_submitting=True)
+    workflow.connect([
+        (inputnode, template_mapper, [
+            ("collected_scalars", "recon_scalars"),
+            ("t1_2_mni_forward_transform", "to_template_transform"),
+            ("resampling_template", "template_reference_image")]),
+        (template_mapper, outputnode, [
+            ("template_space_scalars", "template_scalars")]),
+        (template_mapper, ds_template_scalars, [
+            # Send the resampled files (without metadata) so they don't get cleaned up
+            ("template_space_scalars", "resampled_files"),
+            ("template_space_scalar_info", "recon_scalars")])
+        ])
 
     return workflow
 
 
 def init_scalar_to_surface_wf(omp_nthreads, available_anatomical_data,
-                              name="scalar_to_surface", output_suffix="", params={}):
+                              name="scalar_to_surface", qsirecon_suffix="", params={}):
     """Maps scalar data to a surface."""
     raise NotImplementedError()
     return workflow
diff --git a/qsiprep/workflows/recon/steinhardt.py b/qsiprep/workflows/recon/steinhardt.py
index d2116b35..6247b59c 100644
--- a/qsiprep/workflows/recon/steinhardt.py
+++ b/qsiprep/workflows/recon/steinhardt.py
@@ -22,7 +22,7 @@
 
 def init_steinhardt_order_param_wf(omp_nthreads, available_anatomical_data,
                                    name="sop_recon",
-                                   output_suffix="", params={}):
+                                   qsirecon_suffix="", params={}):
     """Compute Steinhardt order parameters based on ODFs or FODs
 
     Inputs
@@ -78,14 +78,16 @@ def init_steinhardt_order_param_wf(omp_nthreads, available_anatomical_data,
         ])
 
     sop_sinks = {}
-    if output_suffix:
+    if qsirecon_suffix:
         for sop_order in range(2, sop_order + 1, 2):
             key = 'q%d_file' % sop_order
             sop_sinks[key] = pe.Node(
-                ReconDerivativesDataSink(suffix="q-%d_SOP" % sop_order),
-                    name='ds_sop_q%d' % sop_order,
-                    compress=True,
-                    run_without_submitting=True)
+                ReconDerivativesDataSink(
+                    model="steinhardt",
+                    mfp=f"q{sop_order}",
+                    compress=True),
+                name='ds_sop_q%d' % sop_order,
+                run_without_submitting=True)
             workflow.connect(outputnode, key, sop_sinks[key], 'in_file')
 
     workflow.__desc__ = desc
diff --git a/qsiprep/workflows/recon/tortoise.py b/qsiprep/workflows/recon/tortoise.py
index c24f3879..df25fa8b 100644
--- a/qsiprep/workflows/recon/tortoise.py
+++ b/qsiprep/workflows/recon/tortoise.py
@@ -15,6 +15,9 @@
 from ...engine import Workflow
 from ...interfaces.bids import ReconDerivativesDataSink
 from ...interfaces.interchange import recon_workflow_input_fields
+from ...interfaces.recon_scalars import TORTOISEReconScalars, ReconScalarsDataSink
+from ...engine import Workflow
+
 from ...interfaces.recon_scalars import TORTOISEReconScalars
 from qsiprep.interfaces.tortoise import (
     ComputeADMap,
@@ -26,9 +29,10 @@
     ComputeRDMap,
     EstimateMAPMRI,
     EstimateTensor,
-    TORTOISEConvert,
+    TORTOISEConvert
 )
 
+
 LOGGER = logging.getLogger('nipype.interface')
 
 CITATIONS = {
@@ -43,7 +47,7 @@ def init_tortoise_estimator_wf(
     omp_nthreads,
     available_anatomical_data,
     name="tortoise_recon",
-    output_suffix="", params={}):
+    qsirecon_suffix="", params={}):
     """Run estimators from TORTOISE.
 
     This workflow may run ``EstimateTensor`` and/or ``EstimateMAPMRI``
@@ -85,13 +89,12 @@ def init_tortoise_estimator_wf(
                 'csf_txt', 'scalar_image_info', 'recon_scalars']),
         name="outputnode")
     workflow = Workflow(name=name)
-    recon_scalars = pe.Node(TORTOISEReconScalars(workflow_name=name), name="recon_scalars")
+    recon_scalars = pe.Node(
+        TORTOISEReconScalars(qsirecon_suffix=qsirecon_suffix),
+        name="recon_scalars")
     plot_reports = params.pop("plot_reports", True)
-    desc = """TORTOISE Reconstruction
-
-:
-
-Methods implemented in TORTOISE (@tortoisev3) were used for reconstruction. """
+    desc = "TORTOISE Reconstruction\n\n: " + \
+        "Methods implemented in TORTOISE (@tortoisev3) were used for reconstruction. "
 
     tensor_opts = params.get("estimate_tensor", {})
     estimate_tensor_separately = params.get("estimate_tensor_separately", False)
@@ -151,58 +154,6 @@ def init_tortoise_estimator_wf(
             (compute_dt_ad, recon_scalars, [("ad_file", "ad_file")]),
             (compute_dt_li, recon_scalars, [("li_file", "li_file")])
         ])
-        if output_suffix:
-            ds_dt = pe.Node(
-                ReconDerivativesDataSink(extension='.nii.gz',
-                                        desc="DT",
-                                        suffix=output_suffix,
-                                        compress=True),
-                name='ds_dt',
-                run_without_submitting=True)
-            ds_dt_am = pe.Node(
-                ReconDerivativesDataSink(extension='.nii.gz',
-                                        desc="DTAM",
-                                        suffix=output_suffix,
-                                        compress=True),
-                name='ds_dt_am',
-                run_without_submitting=True)
-            ds_dt_fa = pe.Node(
-                ReconDerivativesDataSink(extension='.nii.gz',
-                                         desc="DTFA",
-                                         suffix=output_suffix,
-                                         compress=True),
-                name='ds_dt_fa',
-                run_without_submitting=True)
-            ds_dt_rd = pe.Node(
-                ReconDerivativesDataSink(extension='.nii.gz',
-                                         desc="DTRD",
-                                         suffix=output_suffix,
-                                         compress=True),
-                name='ds_dt_rd',
-                run_without_submitting=True)
-            ds_dt_ad = pe.Node(
-                ReconDerivativesDataSink(extension='.nii.gz',
-                                         desc="DTAD",
-                                         suffix=output_suffix,
-                                         compress=True),
-                name='ds_dt_ad',
-                run_without_submitting=True)
-            ds_dt_li = pe.Node(
-                ReconDerivativesDataSink(extension='.nii.gz',
-                                         desc="DTLI",
-                                         suffix=output_suffix,
-                                         compress=True),
-                name='ds_dt_li',
-                run_without_submitting=True)
-            workflow.connect([
-                (estimate_tensor, ds_dt, [("dt_file", "in_file")]),
-                (estimate_tensor, ds_dt_am, [("am_file", "in_file")]),
-                (compute_dt_fa, ds_dt_fa, [("fa_file", "in_file")]),
-                (compute_dt_rd, ds_dt_rd, [("rd_file", "in_file")]),
-                (compute_dt_ad, ds_dt_ad, [("ad_file", "in_file")]),
-                (compute_dt_li, ds_dt_li, [("li_file", "in_file")])
-            ])
-
 
     # EstimateMAPMRI
     mapmri_opts = params.get("estimate_mapmri", {})
@@ -212,7 +163,7 @@ def init_tortoise_estimator_wf(
     # Set deltas if we have them. Prevent only one from being defined
     if approximate_deltas:
         LOGGER.warning('Both "big_delta" and "small_delta" are required for precise MAPMRI')
-        big_delta = little_delta = traits.undefined
+        big_delta = little_delta = traits.Undefined
     else:
         mapmri_opts["big_delta"], mapmri_opts["small_delta"] = deltas
     mapmri_opts["num_threads"] = omp_nthreads
@@ -243,7 +194,6 @@ def init_tortoise_estimator_wf(
                 ("dt_file", "dt_file"),
                 ("am_file", "a0_file")])])
 
-
     workflow.connect([
         (tortoise_convert, estimate_mapmri,[
             ("bmtxt_file", "bmtxt_file"),
@@ -270,90 +220,15 @@ def init_tortoise_estimator_wf(
             ("ngpar_file", "ngpar_file"),
             ("ngperp_file", "ngperp_file")]),
     ])
-    if output_suffix:
-        ds_map_coeffs = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmri",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_coeffs',
-            run_without_submitting=True)
-        ds_map_uvec = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmriuvec",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_uvec',
-            run_without_submitting=True)
-        ds_map_pa = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmriPA",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_pa',
-            run_without_submitting=True)
-        ds_map_path = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmriPAth",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_path',
+    if qsirecon_suffix:
+        ds_recon_scalars = pe.Node(
+            ReconScalarsDataSink(),
+            name="ds_recon_scalars",
             run_without_submitting=True)
-        ds_map_rtop = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmriRTOP",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_rtop',
-            run_without_submitting=True)
-        ds_map_rtap = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmriRTAP",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_rtap',
-            run_without_submitting=True)
-        ds_map_rtpp = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmriRTPP",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_rtpp',
-            run_without_submitting=True)
-        ds_map_ng = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmriNG",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_ng',
-            run_without_submitting=True)
-        ds_map_ngpar = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmriNGpar",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_ngpar',
-            run_without_submitting=True)
-        ds_map_ngperp = pe.Node(
-            ReconDerivativesDataSink(extension='.nii.gz',
-                                    desc="mapmriNGperp",
-                                    suffix=output_suffix,
-                                    compress=True),
-            name='ds_map_ngperp',
-            run_without_submitting=True)
-        workflow.connect([
-            (estimate_mapmri, ds_map_coeffs, [("coeffs_file", "in_file")]),
-            (estimate_mapmri, ds_map_uvec, [("uvec_file", "in_file")]),
-            (compute_mapmri_pa, ds_map_pa, [("pa_file", "in_file")]),
-            (compute_mapmri_pa, ds_map_path, [("path_file", "in_file")]),
-            (compute_mapmri_rtop, ds_map_rtop, [("rtop_file", "in_file")]),
-            (compute_mapmri_rtop, ds_map_rtap, [("rtap_file", "in_file")]),
-            (compute_mapmri_rtop, ds_map_rtpp, [("rtpp_file", "in_file")]),
-            (compute_mapmri_ng, ds_map_ng, [("ng_file", "in_file")]),
-            (compute_mapmri_ng, ds_map_ngpar, [("ngpar_file", "in_file")]),
-            (compute_mapmri_ng, ds_map_ngperp, [("ngperp_file", "in_file")]),
-        ])
-
-
+        workflow.connect(
+            recon_scalars,
+            "scalar_info",
+            ds_recon_scalars,
+            "recon_scalars")
     workflow.__desc__ = desc
     return workflow
\ No newline at end of file
diff --git a/qsiprep/workflows/recon/utils.py b/qsiprep/workflows/recon/utils.py
index c92ea61a..99ac4d7e 100644
--- a/qsiprep/workflows/recon/utils.py
+++ b/qsiprep/workflows/recon/utils.py
@@ -21,7 +21,7 @@
 
 
 def init_conform_dwi_wf(omp_nthreads, available_anatomical_data,
-                        name="conform_dwi", output_suffix="", params={}):
+                        name="conform_dwi", qsirecon_suffix="", params={}):
     """If data were preprocessed elsewhere, ensure the gradients and images
     conform to LPS+ before running other parts of the pipeline."""
     inputnode = pe.Node(niu.IdentityInterface(fields=recon_workflow_input_fields),
@@ -45,7 +45,7 @@ def init_conform_dwi_wf(omp_nthreads, available_anatomical_data,
 
 
 def init_discard_repeated_samples_wf(omp_nthreads, available_anatomical_data,
-                                     name="discard_repeats", output_suffix="",
+                                     name="discard_repeats", qsirecon_suffix="",
                                      space="T1w", params={}):
     """Remove a sample if a similar direction/gradient has already been sampled."""
     inputnode = pe.Node(niu.IdentityInterface(fields=recon_workflow_input_fields),


Test patch:
diff --git a/qsiprep/data/pipelines/bundle_scalar_map.json b/qsiprep/data/pipelines/test_scalar_maps.json
similarity index 69%
rename from qsiprep/data/pipelines/bundle_scalar_map.json
rename to qsiprep/data/pipelines/test_scalar_maps.json
index 2c342804..dbb9d735 100644
--- a/qsiprep/data/pipelines/bundle_scalar_map.json
+++ b/qsiprep/data/pipelines/test_scalar_maps.json
@@ -5,22 +5,32 @@
   "anatomical": [ ],
   "nodes": [
     {
-      "name": "DIPYdki",
+      "name": "dipy_dki",
       "software": "Dipy",
       "action": "DKI_reconstruction",
       "input": "qsiprep",
-      "output_suffix": "DKI",
+      "qsirecon_suffix": "DIPYDKI",
       "parameters": {
         "write_mif": false,
         "write_fibgz": false
       }
     },
+    {
+      "name": "template_map",
+      "software": "qsiprep",
+      "action": "template_map",
+      "input": "qsiprep",
+      "scalars_from": ["dipy_dki", "gqi_scalars"],
+      "parameters": {
+        "interpolation": "NearestNeighbor"
+      }
+    },
     {
       "name": "dsistudio_gqi",
       "software": "DSI Studio",
       "action": "reconstruction",
       "input": "qsiprep",
-      "output_suffix": "gqi",
+      "qsirecon_suffix": "DSIStudio",
       "parameters": {"method": "gqi"}
     },
     {
@@ -28,7 +38,7 @@
       "software": "DSI Studio",
       "action": "autotrack",
       "input": "dsistudio_gqi",
-      "output_suffix": "AutoTrackGQI",
+      "qsirecon_suffix": "DSIStudio",
       "parameters": {
         "track_id": "Fasciculus,Cingulum,Aslant,Corticos,Thalamic_R,Reticular,Optic,Fornix,Corpus",
         "tolerance": "22,26,30",
@@ -41,15 +51,14 @@
       "software": "DSI Studio",
       "action": "export",
       "input": "dsistudio_gqi",
-      "output_suffix": "gqiscalar"
+      "qsirecon_suffix": "DSIStudio"
     },
     {
       "name": "bundle_means",
       "software": "qsiprep",
       "action": "bundle_map",
       "input": "autotrackgqi",
-      "scalars_from": ["gqi_scalars", "DIPYdki"],
-      "output_suffix": "bundlemap"
+      "scalars_from": ["gqi_scalars", "dipy_dki"]
     }
   ]
 }