diff --git a/.config/b0_to_template_affine_config.py b/.config/b0_to_template_affine_config.py
index 7efa503..bc7f160 100644
--- a/.config/b0_to_template_affine_config.py
+++ b/.config/b0_to_template_affine_config.py
@@ -14,7 +14,7 @@
 
 c.AntsRegistration.log_level = 30
 
-c.AntsRegistration.init_with_ants_ai = True
+c.AntsRegistration.init_with_ants_ai = False
 
 c.AntsRegistration.base_config_file = ""
 
@@ -24,6 +24,10 @@
 # AntsConfiguration(mrHARDIConfigurable) configuration
 # -----------------------------------------------------------------------------
 
+c.AntsConfiguration.coarse_angular_split = 3
+
+c.AntsConfiguration.fine_angular_split = 4
+
 c.AntsConfiguration.accross_modalities = True
 
 c.AntsConfiguration.dimension = 3
@@ -39,20 +43,31 @@
 c.AntsConfiguration.match_histogram = False
 
 c.AntsConfiguration.passes = [{
-    "conv_eps": 1e-5,
+    "conv_eps": 1e-7,
     "conv_max_iter": [400, 200, 100, 50],
     "conv_win": 20,
-    "grad_step": 0.1,
+    "grad_step": 0.05,
     "klass": "mrHARDI.traits.ants.AntsRigid",
     "metrics": [
+        {
+            "target_index": 0,
+            "moving_index": 0,
+            "args": [
+                0.75,
+                128,
+                "Regular",
+                1.
+            ],
+            "klass": "mrHARDI.traits.ants.MetricMI"
+        },
         {
             "target_index": 0,
             "moving_index": 1,
             "args": [
-                1.,
-                64,
+                0.25,
+                128,
                 "Regular",
-                0.7,
+                1.,
                 True
             ],
             "klass": "mrHARDI.traits.ants.MetricMI"
@@ -71,20 +86,31 @@
         0
     ]
 }, {
-    "conv_eps": 1e-5,
+    "conv_eps": 1e-7,
     "conv_max_iter": [500, 300, 150, 75],
     "conv_win": 10,
-    "grad_step": 0.1,
+    "grad_step": 0.05,
     "klass": "mrHARDI.traits.ants.AntsAffine",
     "metrics": [
+        {
+            "target_index": 0,
+            "moving_index": 0,
+            "args": [
+                0.75,
+                128,
+                "Regular",
+                1.
+            ],
+            "klass": "mrHARDI.traits.ants.MetricMI"
+        },
         {
             "target_index": 0,
             "moving_index": 1,
             "args": [
-                1.,
-                64,
+                0.25,
+                128,
                 "Regular",
-                0.8,
+                1.,
                 True
             ],
             "klass": "mrHARDI.traits.ants.MetricMI"
diff --git a/.config/b0_to_template_syn_config.py b/.config/b0_to_template_syn_config.py
index f6662e6..cc239f6 100644
--- a/.config/b0_to_template_syn_config.py
+++ b/.config/b0_to_template_syn_config.py
@@ -35,17 +35,17 @@
 c.AntsConfiguration.match_histogram = False
 
 c.AntsConfiguration.passes = [{
-    "conv_eps": 1e-5,
-    "conv_max_iter": [100, 50, 20, 10, 10],
+    "conv_eps": 1e-7,
+    "conv_max_iter": [200, 200, 140, 100, 40],
     "conv_win": 10,
-    "grad_step": 0.2,
+    "grad_step": 0.05,
     "var_penality": 3,
-    "var_total": 0,
+    "var_total": 3,
     "klass": "mrHARDI.traits.ants.AntsSyN",
     "metrics": [
         {
             "target_index": 0,
-            "moving_index": 1,
+            "moving_index": 0,
             "args": [
                 0.65,
                 64,
@@ -54,49 +54,15 @@
             ],
             "klass": "mrHARDI.traits.ants.MetricMI"
         },
-        {
-            "target_index": 0,
-            "moving_index": 2,
-            "args": [
-                0.35,
-                2,
-                "Regular",
-                1.
-            ],
-            "klass": "mrHARDI.traits.ants.MetricCC"
-        }
-    ],
-    "shrinks": [
-        12,
-        8,
-        4,
-        2,
-        1
-    ],
-    "smoothing": [
-        3.,
-        1.5,
-        1.,
-        0.5,
-        0
-    ]
-}, {
-    "conv_eps": 1e-5,
-    "conv_max_iter": [50, 20, 20],
-    "conv_win": 20,
-    "grad_step": 0.1,
-    "var_penality": 3,
-    "var_total": 0,
-    "klass": "mrHARDI.traits.ants.AntsSyN",
-    "metrics": [
         {
             "target_index": 0,
             "moving_index": 1,
             "args": [
-                0.65,
+                0.15,
                 64,
                 "Regular",
-                1.
+                1.,
+                True
             ],
             "klass": "mrHARDI.traits.ants.MetricMI"
         },
@@ -104,8 +70,8 @@
             "target_index": 0,
             "moving_index": 2,
             "args": [
-                0.35,
-                2,
+                0.20,
+                4,
                 "Regular",
                 1.
             ],
@@ -113,14 +79,18 @@
         }
     ],
     "shrinks": [
-        8,
+        10,
+        6,
         4,
+        2,
         1
     ],
     "smoothing": [
-        1.5,
+        5.,
+        3.,
+        2.,
         1.,
-        0
+        0.
     ]
 }]
 
diff --git a/.config/dwi_n4_normalization_config.py b/.config/dwi_n4_normalization_config.py
index 05c51a6..a734bde 100644
--- a/.config/dwi_n4_normalization_config.py
+++ b/.config/dwi_n4_normalization_config.py
@@ -23,21 +23,21 @@
 
 # -----------------------------------------------------------------------------
 # N4BiasCorrectionConfiguration(mrHARDIConfigurable) configuration
-c.N4BiasCorrectionConfiguration.bins = 150
+c.N4BiasCorrectionConfiguration.bins = 200
 
-c.N4BiasCorrectionConfiguration.filter_width = 0.075
+c.N4BiasCorrectionConfiguration.filter_width = 0.60
 
 c.N4BiasCorrectionConfiguration.spline_order = 3
 
-c.N4BiasCorrectionConfiguration.knot_distance = 64.0
+c.N4BiasCorrectionConfiguration.nvox_between_knots = 7.0
 
-c.N4BiasCorrectionConfiguration.iterations = [100, 100, 50, 30]
+c.N4BiasCorrectionConfiguration.iterations = [300, 300, 150, 75, 40]
 
 c.N4BiasCorrectionConfiguration.noise = 0.01
 
 c.N4BiasCorrectionConfiguration.rescale = True
 
-c.N4BiasCorrectionConfiguration.shrink = 2
+c.N4BiasCorrectionConfiguration.shrink = 4
 
 c.N4BiasCorrectionConfiguration.threshold = 1E-8
 
diff --git a/.config/dwi_n4_normalization_quick_config.py b/.config/dwi_n4_normalization_quick_config.py
index 1f293e0..4f7e55f 100644
--- a/.config/dwi_n4_normalization_quick_config.py
+++ b/.config/dwi_n4_normalization_quick_config.py
@@ -25,11 +25,13 @@
 # N4BiasCorrectionConfiguration(mrHARDIConfigurable) configuration
 c.N4BiasCorrectionConfiguration.bins = 150
 
-c.N4BiasCorrectionConfiguration.filter_width = 0.075
+c.N4BiasCorrectionConfiguration.filter_width = 0.6
 
-c.N4BiasCorrectionConfiguration.knot_distance = 4.0
+c.N4BiasCorrectionConfiguration.spline_order = 3
 
-c.N4BiasCorrectionConfiguration.iterations = [300, 150, 75, 40]
+c.N4BiasCorrectionConfiguration.nvox_between_knots = 7.0
+
+c.N4BiasCorrectionConfiguration.iterations = [300, 300, 150, 75, 40]
 
 c.N4BiasCorrectionConfiguration.noise = 0.01
 
diff --git a/.config/prepare_topup_base_config.py b/.config/prepare_topup_base_config.py
index 484ba06..615f94f 100644
--- a/.config/prepare_topup_base_config.py
+++ b/.config/prepare_topup_base_config.py
@@ -28,74 +28,74 @@
 c.TopupConfiguration.klass = "mrHARDI.config.topup.TopupConfiguration"
 
 c.TopupConfiguration.passes = [{
-    "warp_resolution": 8.8,
-    "subsampling": 2,
-    "blur_fwhm": 3.5,
-    "n_iter": 100,
+    "warp_resolution": 25.1,
+    "subsampling": 4,
+    "blur_fwhm": 10,
+    "n_iter": 10,
     "estimate_motion": 1,
     "minimizer": 0,
     "w_reg": 5E-3
 }, {
-    "warp_resolution": 7.1,
-    "subsampling": 2,
-    "blur_fwhm": 2.6,
-    "n_iter": 100,
+    "warp_resolution": 20.3,
+    "subsampling": 4,
+    "blur_fwhm": 7.4,
+    "n_iter": 10,
     "estimate_motion": 1,
     "minimizer": 0,
     "w_reg": 1e-3
 }, {
-    "warp_resolution": 6.2,
-    "subsampling": 2,
-    "blur_fwhm": 1.8,
-    "n_iter": 100,
+    "warp_resolution": 17.7,
+    "subsampling": 4,
+    "blur_fwhm": 5.1,
+    "n_iter": 10,
     "estimate_motion": 1,
     "minimizer": 0,
     "w_reg": 1e-4
 }, {
-    "warp_resolution": 5.3,
-    "subsampling": 2,
-    "blur_fwhm": 1.3,
-    "n_iter": 100,
+    "warp_resolution": 15.1,
+    "subsampling": 4,
+    "blur_fwhm": 3.7,
+    "n_iter": 10,
     "estimate_motion": 1,
     "minimizer": 0,
     "w_reg": 1.5e-5
 }, {
-    "warp_resolution": 4.4,
+    "warp_resolution": 12.6,
     "subsampling": 2,
-    "blur_fwhm": 1.3,
-    "n_iter": 100,
+    "blur_fwhm": 3.7,
+    "n_iter": 15,
     "estimate_motion": 1,
     "minimizer": 0,
     "w_reg": 5e-6
 }, {
-    "warp_resolution": 2.7,
-    "subsampling": 1,
-    "blur_fwhm": 0.9,
-    "n_iter": 100,
+    "warp_resolution": 7.7,
+    "subsampling": 2,
+    "blur_fwhm": 2.6,
+    "n_iter": 15,
     "estimate_motion": 0,
     "minimizer": 1,
     "w_reg": 5e-7
 }, {
-    "warp_resolution": 1.7,
-    "subsampling": 1,
-    "blur_fwhm": 0.4,
-    "n_iter": 150,
+    "warp_resolution": 6.1,
+    "subsampling": 2,
+    "blur_fwhm": 1.1,
+    "n_iter": 20,
     "estimate_motion": 0,
     "minimizer": 1,
     "w_reg": 5e-8
 }, {
-    "warp_resolution": 1.7,
+    "warp_resolution": 6.1,
     "subsampling": 1,
     "blur_fwhm": 0.,
-    "n_iter": 200,
+    "n_iter": 20,
     "estimate_motion": 0,
     "minimizer": 1,
     "w_reg": 5e-10
 }, {
-    "warp_resolution": 1.5,
+    "warp_resolution": 6.1,
     "subsampling": 1,
     "blur_fwhm": 0.,
-    "n_iter": 400,
+    "n_iter": 40,
     "estimate_motion": 0,
     "minimizer": 1,
     "w_reg": 1e-11
@@ -105,7 +105,7 @@
 
 c.TopupConfiguration.reg_model = "bending_energy"
 
-c.TopupConfiguration.scale_intensities = False
+c.TopupConfiguration.scale_intensities = True
 
 c.TopupConfiguration.spl_order = "cubic"
 
diff --git a/.config/t1_n4_normalization_config.py b/.config/t1_n4_normalization_config.py
index a77e744..edc2195 100644
--- a/.config/t1_n4_normalization_config.py
+++ b/.config/t1_n4_normalization_config.py
@@ -25,21 +25,21 @@
 # N4BiasCorrectionConfiguration(mrHARDIConfigurable) configuration
 c.N4BiasCorrectionConfiguration.bins = 200
 
-c.N4BiasCorrectionConfiguration.filter_width = 0.075
+c.N4BiasCorrectionConfiguration.filter_width = 0.60
 
 c.N4BiasCorrectionConfiguration.spline_order = 3
 
-c.N4BiasCorrectionConfiguration.knot_distance = 16.0
+c.N4BiasCorrectionConfiguration.nvox_between_knots = 14.0
 
-c.N4BiasCorrectionConfiguration.iterations = [80, 80, 40, 30]
+c.N4BiasCorrectionConfiguration.iterations = [300, 300, 150, 70]
 
 c.N4BiasCorrectionConfiguration.noise = 0.01
 
 c.N4BiasCorrectionConfiguration.rescale = True
 
-c.N4BiasCorrectionConfiguration.shrink = 2
+c.N4BiasCorrectionConfiguration.shrink = 4
 
-c.N4BiasCorrectionConfiguration.threshold = 1E-6
+c.N4BiasCorrectionConfiguration.threshold = 1E-8
 
 # Base traits configuration
 
diff --git a/.config/t1_to_b0_registration_affine_config.py b/.config/t1_to_b0_registration_affine_config.py
index bf41d20..9ecae3b 100644
--- a/.config/t1_to_b0_registration_affine_config.py
+++ b/.config/t1_to_b0_registration_affine_config.py
@@ -14,7 +14,7 @@
 
 c.AntsRegistration.log_level = 30
 
-c.AntsRegistration.init_with_ants_ai = True
+c.AntsRegistration.init_with_ants_ai = False
 
 c.AntsRegistration.base_config_file = ""
 
@@ -42,20 +42,32 @@
 c.AntsConfiguration.match_histogram = False
 
 c.AntsConfiguration.passes = [{
-    "conv_eps": 1e-07,
+    "conv_eps": 1e-6,
     "conv_max_iter": [500, 300, 150, 50],
-    "conv_win": 10,
+    "conv_win": 20,
     "grad_step": 0.1,
     "klass": "mrHARDI.traits.ants.AntsRigid",
     "metrics": [
+        {
+            "target_index": 0,
+            "moving_index": 0,
+            "args": [
+                0.75,
+                128,
+                "Regular",
+                1.0
+            ],
+            "klass": "mrHARDI.traits.ants.MetricMI"
+        },
         {
             "target_index": 1,
             "moving_index": 0,
             "args": [
-                1.0,
-                64,
+                0.25,
+                128,
                 "Regular",
-                0.5
+                1.0,
+                True
             ],
             "klass": "mrHARDI.traits.ants.MetricMI"
         }
@@ -67,26 +79,38 @@
         1
     ],
     "smoothing": [
-        1.5,
+        3,
+        2,
         1,
-        0.5,
         0
     ]
 }, {
-    "conv_eps": 1e-07,
+    "conv_eps": 1e-6,
     "conv_max_iter": [300, 150, 100, 50],
-    "conv_win": 10,
+    "conv_win": 20,
     "grad_step": 0.1,
     "klass": "mrHARDI.traits.ants.AntsAffine",
     "metrics": [
+        {
+            "target_index": 0,
+            "moving_index": 0,
+            "args": [
+                0.75,
+                128,
+                "Regular",
+                1.0
+            ],
+            "klass": "mrHARDI.traits.ants.MetricMI"
+        },
         {
             "target_index": 1,
             "moving_index": 0,
             "args": [
-                1.0,
-                64,
+                0.25,
+                128,
                 "Regular",
-                0.5
+                1.0,
+                True
             ],
             "klass": "mrHARDI.traits.ants.MetricMI"
         }
@@ -98,9 +122,9 @@
         1
     ],
     "smoothing": [
-        1.5,
+        3,
+        2,
         1,
-        0.5,
         0
     ]
 }]
diff --git a/.config/t1_to_b0_syn_config.py b/.config/t1_to_b0_syn_config.py
index 7d27dc4..e80c9ba 100644
--- a/.config/t1_to_b0_syn_config.py
+++ b/.config/t1_to_b0_syn_config.py
@@ -36,9 +36,9 @@
 
 c.AntsConfiguration.passes = [{
     "conv_eps": 1e-6,
-    "conv_max_iter": [30, 30],
+    "conv_max_iter": [200, 100, 40, 20, 10],
     "conv_win": 10,
-    "grad_step": 0.1,
+    "grad_step": 0.05,
     "var_penality": 3,
     "var_total": 0,
     "klass": "mrHARDI.traits.ants.AntsSyN",
@@ -47,19 +47,19 @@
             "target_index": 0,
             "moving_index": 0,
             "args": [
-                0.65,
-                2,
+                0.75,
+                128,
                 "Regular",
                 1.
             ],
-            "klass": "mrHARDI.traits.ants.MetricCC"
+            "klass": "mrHARDI.traits.ants.MetricMI"
         },
         {
             "target_index": 1,
             "moving_index": 0,
             "args": [
-                0.35,
-                2,
+                0.25,
+                8,
                 "Regular",
                 1.
             ],
@@ -67,12 +67,18 @@
         }
     ],
     "shrinks": [
+        10,
+        6,
+        4,
         2,
         1
     ],
     "smoothing": [
-        0.5,
-        0
+        5.,
+        3.,
+        2.,
+        1.,
+        0.
     ]
 }]
 
diff --git a/.config/t1_to_template_affine_config.py b/.config/t1_to_template_affine_config.py
index 271581a..6e211b8 100644
--- a/.config/t1_to_template_affine_config.py
+++ b/.config/t1_to_template_affine_config.py
@@ -14,11 +14,9 @@
 
 c.AntsRegistration.log_level = 30
 
-c.AntsRegistration.init_with_ants_ai = True
-
 c.AntsRegistration.base_config_file = ""
 
-c.AntsRegistration.init_with_ants_ai = True
+c.AntsRegistration.init_with_ants_ai = False
 
 c.AntsRegistration.verbose = True
 
@@ -28,7 +26,7 @@
 
 c.AntsConfiguration.coarse_angular_split = 3
 
-c.AntsConfiguration.coarse_linear_split = 3
+c.AntsConfiguration.fine_angular_split = 4
 
 c.AntsConfiguration.accross_modalities = False
 
@@ -45,64 +43,45 @@
 c.AntsConfiguration.match_histogram = True
 
 c.AntsConfiguration.passes = [{
-    "conv_eps": 1e-6,
-    "conv_max_iter": [400, 100, 50, 20],
-    "conv_win": 10,
-    "grad_step": 0.1,
-    "klass": "mrHARDI.traits.ants.AntsRigid",
+    "conv_eps": 1e-5,
+    "conv_max_iter": [1000, 400, 200, 100, 100],
+    "conv_win": 20,
+    "grad_step": 0.05,
+    "klass": "mrHARDI.traits.ants.AntsAffine",
     "metrics": [
         {
             "target_index": 0,
             "moving_index": 0,
             "args": [
-                1.0,
-                64,
+                0.75,
+                2,
                 "Regular",
-                0.8,
-                True
+                1.0
             ],
-            "klass": "mrHARDI.traits.ants.MetricMI"
-        }
-    ],
-    "shrinks": [
-        8,
-        4,
-        2,
-        1
-    ],
-    "smoothing": [
-        3.,
-        2.,
-        1.,
-        0
-    ]
-}, {
-    "conv_eps": 1e-6,
-    "conv_max_iter": [200, 100, 50, 20],
-    "conv_win": 20,
-    "grad_step": 0.1,
-    "klass": "mrHARDI.traits.ants.AntsAffine",
-    "metrics": [
+            "klass": "mrHARDI.traits.ants.MetricCC"
+        },
         {
             "target_index": 0,
             "moving_index": 0,
             "args": [
-                1.0,
-                64,
+                0.25,
+                32,
                 "Regular",
-                0.8,
+                1.0,
                 True
             ],
             "klass": "mrHARDI.traits.ants.MetricMI"
         }
     ],
     "shrinks": [
-        8,
+        10,
+        6,
         4,
         2,
         1
     ],
     "smoothing": [
+        5.,
         3.,
         2.,
         1.,
diff --git a/.config/t1_to_template_syn_config.py b/.config/t1_to_template_syn_config.py
index 3f35c0b..f91e774 100644
--- a/.config/t1_to_template_syn_config.py
+++ b/.config/t1_to_template_syn_config.py
@@ -36,9 +36,9 @@
 
 c.AntsConfiguration.passes = [{
     "conv_eps": 1e-7,
-    "conv_max_iter": [50, 50, 20, 20],
+    "conv_max_iter": [200, 200, 140, 100, 40],
     "conv_win": 10,
-    "grad_step": 0.2,
+    "grad_step": 0.1,
     "var_penality": 3,
     "var_total": 0,
     "klass": "mrHARDI.traits.ants.AntsSyN",
@@ -47,53 +47,38 @@
             "target_index": 0,
             "moving_index": 0,
             "args": [
-                1.0,
-                64,
+                0.6,
+                4,
                 "Regular",
                 1.0
             ],
-            "klass": "mrHARDI.traits.ants.MetricMI"
-        }
-    ],
-    "shrinks": [
-        8,
-        4,
-        2,
-        1
-    ],
-    "smoothing": [
-        4.,
-        2.,
-        1.,
-        0.
-    ]
-}, {
-    "conv_eps": 1e-7,
-    "conv_max_iter": [20, 10],
-    "conv_win": 10,
-    "grad_step": 0.1,
-    "var_penality": 3,
-    "var_total": 0,
-    "klass": "mrHARDI.traits.ants.AntsSyN",
-    "metrics": [
+            "klass": "mrHARDI.traits.ants.MetricCC"
+        },
         {
             "target_index": 0,
             "moving_index": 0,
             "args": [
-                1.0,
-                4,
+                0.4,
+                64,
                 "Regular",
-                1.0
+                1.0,
+                True
             ],
-            "klass": "mrHARDI.traits.ants.MetricCC"
+            "klass": "mrHARDI.traits.ants.MetricMI"
         }
     ],
     "shrinks": [
+        10,
+        6,
         4,
+        2,
         1
     ],
     "smoothing": [
+        5.,
+        3.,
         2.,
+        1.,
         0.
     ]
 }]
diff --git a/.data/maccaca_mulatta/tissue_segmentation/tissue_segmentation_mask_no_bv_dilated.nii.gz b/.data/maccaca_mulatta/tissue_segmentation/tissue_segmentation_mask_no_bv_dilated.nii.gz
index 4f9ef3f..0ec15d3 100644
Binary files a/.data/maccaca_mulatta/tissue_segmentation/tissue_segmentation_mask_no_bv_dilated.nii.gz and b/.data/maccaca_mulatta/tissue_segmentation/tissue_segmentation_mask_no_bv_dilated.nii.gz differ
diff --git a/README.md b/README.md
index 7ecaf97..f783469 100644
--- a/README.md
+++ b/README.md
@@ -27,24 +27,20 @@ Frontiers in Neuroinformatics, 10.3389/fninf.2023.1191200, doi.org/10.3389/fninf
 To run the pipeline, the following tools must be installed :
 
 - [Nextflow](https://www.nextflow.io/docs/latest/getstarted.html) (we recommend installing the last available release of version 21)
-- [Singularity](https://docs.sylabs.io/guides/3.0/user-guide/installation.html) 3.7.1 or higher, or [Docker](https://docs.docker.com/engine/install/)
+- [Apptainer](https://apptainer.org/docs/admin/main/installation.html) 1.3.2 or higher, or [Docker](https://docs.docker.com/engine/install/) (we recommend [Docker Desktop](https://docs.docker.com/desktop/))
 
 ## Dependencies
 
-All dependencies required to run the pipeline have been packaged in singularity 
+All dependencies required to run the pipeline have been packaged in apptainer 
 and docker images. The `latest` versions come pre-packaged with the CUDA runtime 
 and require a Nvidia GPU to execute. For usage on a machine without a Nvidia GPU, 
 use the images tagged `nogpu`.
 
 - Docker : [Docker Hub](https://hub.docker.com/r/avcaron/versa)
   - `docker pull avcaron/versa:latest`
-- Singularity :
-  - Singularity images are no longer produced in house. To build your singularity, 
-    use a docker tag and the following command : `singularity build <image.sif> docker://avcaron/versa:<tag>`.
-  - Old versions of the singularity images can still be access through the ORAS 
-    container registry :
-    - Nvidia GPU  : `singularity pull oras://mrhardi.azurecr.io/mrHARDI/mrhardi:latest`
-    - Without GPU : `singularity pull oras://mrhardi.azurecr.io/mrHARDI/mrhardi:nogpu`
+- Apptainer : apptainer images must be built locally, using the following command :
+  - `apptainer build <image.sif> docker://avcaron/versa:latest .`
+
 
 # Data input format
 
@@ -153,7 +149,7 @@ nextflow run \
 Additional parameters to this command can be supplied, such as :
 
 - Running with docker image : `-with-docker <docker image>`
-- Running with singularity image : `-with-singularity <singularity image>`
+- Running with apptainer image : `-with-apptainer <apptainer image>`
 - Changing output directory : `--output_root <output folder>`
 - Display help and usage : `--help`
 
diff --git a/main.nf b/main.nf
index f8eb7bc..9481057 100644
--- a/main.nf
+++ b/main.nf
@@ -3,6 +3,10 @@
 nextflow.enable.dsl=2
 
 params.help = false
+params.preprocess_images = true
+params.reconstruct_diffusion = true
+params.measure_diffusion = true
+params.pft_tracking = true
 
 //include { print_channel } from "./modules/debug.nf" // For debugging purpose only
 include { load_dataset } from "./workflows/io.nf"
@@ -18,34 +22,40 @@ workflow {
         display_run_info()
         dataloader = load_dataset()
 
-        preprocess_wkf(
-            dataloader.dwi,
-            dataloader.rev,
-            dataloader.anat,
-            dataloader.pvf,
-            dataloader.metadata,
-            dataloader.rev_metadata,
-            dataloader.dwi_mask,
-            dataloader.anat_mask
-        )
+        if ( params.preprocess_images ) {
+            preprocess_wkf(
+                dataloader.dwi,
+                dataloader.rev,
+                dataloader.anat,
+                dataloader.pvf,
+                dataloader.metadata,
+                dataloader.rev_metadata,
+                dataloader.dwi_mask,
+                dataloader.anat_mask
+            )
+        }
 
-        reconstruct_wkf(
-            preprocess_wkf.out.dwi,
-            preprocess_wkf.out.mask,
-            preprocess_wkf.out.tissue_masks,
-            preprocess_wkf.out.safe_wm_mask,
-            preprocess_wkf.out.metadata
-        )
+        if ( params.reconstruct_diffusion || params.measure_diffusion || params.pft_tracking ) {
+            reconstruct_wkf(
+                preprocess_wkf.out.dwi,
+                preprocess_wkf.out.mask,
+                preprocess_wkf.out.tissue_masks,
+                preprocess_wkf.out.safe_wm_mask,
+                preprocess_wkf.out.metadata
+            )
+        }
 
-        measure_wkf(
-            preprocess_wkf.out.dwi,
-            preprocess_wkf.out.mask,
-            preprocess_wkf.out.tissue_masks,
-            reconstruct_wkf.out.csd,
-            reconstruct_wkf.out.diamond,
-            reconstruct_wkf.out.diamond_summary,
-            preprocess_wkf.out.metadata
-        )
+        if ( params.measure_diffusion || params.pft_tracking ) {
+            measure_wkf(
+                preprocess_wkf.out.dwi,
+                preprocess_wkf.out.mask,
+                preprocess_wkf.out.tissue_masks,
+                reconstruct_wkf.out.csd,
+                reconstruct_wkf.out.diamond,
+                reconstruct_wkf.out.diamond_summary,
+                preprocess_wkf.out.metadata
+            )
+        }
 
         if ( params.pft_tracking ) {
             tracking_wkf(reconstruct_wkf.out.csd, preprocess_wkf.out.pvf)
diff --git a/modules/functions.nf b/modules/functions.nf
index 180f7e4..c500742 100644
--- a/modules/functions.nf
+++ b/modules/functions.nf
@@ -205,7 +205,7 @@ def is_path_list ( pth ) {
 
 def collect_paths ( data_channel ) {
     return data_channel.map{
-        [it[0], (it.size() > 2 ? it[1..-1] : [it[1]] ).findAll{ it }]
+        [it[0], (it.size() > 2 ? it[1..-1] : [it[1]].flatten() ).findAll{ it }]
     }
 }
 
diff --git a/modules/processes/denoise.nf b/modules/processes/denoise.nf
index c6fbd31..5394a3b 100644
--- a/modules/processes/denoise.nf
+++ b/modules/processes/denoise.nf
@@ -60,7 +60,7 @@ process nlmeans_denoise {
         tuple val(sid), path("${image.simpleName}__nlmeans_denoised_metadata.*"), optional: true, emit: metadata
     script:
         def args = ""
-        if ( !mask.empty() ) args += "--mask $mask"
+        // if ( !mask.empty() ) args += "--mask $mask"
         def after_script = ""
         if ( !metadata.empty() ) after_script += "cp $metadata ${image.simpleName}__nlmeans_denoised_metadata.py"
         """
@@ -70,6 +70,7 @@ process nlmeans_denoise {
         mrhardi nlmeans \
             --in $image \
             --out ${image.simpleName}__nlmeans_denoised \
+            --coils 64 \
             --processes $task.cpus $args
         $after_script
         """
@@ -106,68 +107,74 @@ process ants_gaussian_denoise {
 }
 
 process n4_denoise {
-    label "N4_CORRECTION"
-    label params.conservative_resources ? "res_conservative_cpu" : "res_max_cpu"
 
-    publishDir "${params.output_root}/all/${sid}/$caller_name/${task.process.replaceAll(":", "/")}", mode: "$params.publish_all_mode", enabled: params.publish_all, overwrite: true
-    publishDir "${params.output_root}/${sid}", saveAs: {
-        f -> ("$publish" == "true") ? f.contains("metadata") || f.contains("bias_field") ? null 
-                                                                                         : remove_alg_suffixes(f)
-                                    : null
-    }, mode: params.publish_mode, overwrite: true
+label "N4_CORRECTION"
+label params.conservative_resources ? "res_conservative_cpu" : "res_max_cpu"
+
+publishDir "${params.output_root}/all/${sid}/$caller_name/${task.process.replaceAll(":", "/")}", mode: "$params.publish_all_mode", enabled: params.publish_all, overwrite: true
+publishDir "${params.output_root}/${sid}", saveAs: {
+    f -> ("$publish" == "true") ? f.contains("metadata") || f.contains("bias_field") ? null 
+                                                                                        : remove_alg_suffixes(f)
+                                : null
+}, mode: params.publish_mode, overwrite: true
+
+input:
+    tuple val(sid), path(image), file(anat), file(mask), file(metadata)
+    val(caller_name)
+    val(publish)
+    path(config)
+output:
+    tuple val(sid), path("${image.simpleName}__n4denoised.nii.gz"), emit: image
+    tuple val(sid), path("${image.simpleName}__n4denoised_metadata.*"), optional: true, emit: metadata
+    tuple val(sid), path("${image.simpleName}_n4_bias_field.nii.gz"), emit: bias_field
+script:
+def before_denoise =""
+def after_denoise = ""
+def args = ""
+if ( anat.empty() ) {
+    before_denoise += "fslmaths $image -thr 0 image4n4.nii.gz\n"
+    after_denoise += "mv n4denoise_bias_field.nii.gz ${image.simpleName}_n4_bias_field.nii.gz\n"
+}
+else {
+    args += "--apply $image"
+    before_denoise += "fslmaths $anat -thr 0 image4n4.nii.gz\n"
+    after_denoise += "mv tmp_n4denoised_bias_field.nii.gz ${image.simpleName}_n4_bias_field.nii.gz\n"
+}
 
-    input:
-        tuple val(sid), path(image), file(anat), file(mask), file(metadata)
-        val(caller_name)
-        val(publish)
-        path(config)
-    output:
-        tuple val(sid), path("${image.simpleName}__n4denoised.nii.gz"), emit: image
-        tuple val(sid), path("${image.simpleName}__n4denoised_metadata.*"), optional: true, emit: metadata
-        tuple val(sid), path("${image.simpleName}_n4_bias_field.nii.gz"), emit: bias_field
-    script:
-        def before_denoise =""
-        def after_denoise = ""
-        def args = ""
-        if ( anat.empty() ) {
-            args += "--in $image"
-            after_denoise += "mv n4denoise_bias_field.nii.gz ${image.simpleName}_n4_bias_field.nii.gz\n"
-        }
-        else {
-            args += "--in $anat --apply $image"
-            after_denoise += "mv tmp_n4denoised_bias_field.nii.gz ${image.simpleName}_n4_bias_field.nii.gz\n"
-        }
+if ( !metadata.empty() ) {
+    after_denoise += "mv n4denoise_metadata.py ${image.simpleName}__n4denoised_metadata.py\n"
+    args += " --metadata $metadata"
+}
 
-        if ( !metadata.empty() ) {
-            after_denoise += "mv n4denoise_metadata.py ${image.simpleName}__n4denoised_metadata.py\n"
-            args += " --metadata $metadata"
-        }
-        after_denoise += "fslmaths n4denoise.nii.gz -thr 0 ${image.simpleName}__n4denoised.nii.gz\n"
+if ( !mask.empty() ) {
+    // args += " --mask mask_for_n4.nii.gz"
+    if ( !anat.empty() ) {
+        before_denoise += "antsApplyTransforms -n NearestNeighbor -e 0 -r $anat -t identity -i $mask -o mask_for_n4.nii.gz\n"
+        before_denoise += "scil_image_math.py convert mask_for_n4.nii.gz mask_for_n4.nii.gz -f --data_type uint8\n"
+    }
+    else {
+        before_denoise += "antsApplyTransforms -n NearestNeighbor -e 0 -r $image -t identity -i $mask -o mask_for_n4.nii.gz\n"
+        before_denoise += "scil_image_math.py convert mask_for_n4.nii.gz mask_for_n4.nii.gz -f --data_type uint8\n"
+    }
+}
 
-        if ( !mask.empty() ) {
-            args += " --mask mask_for_n4.nii.gz"
-            if ( !anat.empty() ) {
-                before_denoise += "antsApplyTransforms -n NearestNeighbor -e 0 -r $anat -t identity -i $mask -o mask_for_n4.nii.gz\n"
-                before_denoise += "scil_image_math.py convert mask_for_n4.nii.gz mask_for_n4.nii.gz -f --data_type uint8\n"
-            }
-            else {
-                before_denoise += "antsApplyTransforms -n NearestNeighbor -e 0 -r $image -t identity -i $mask -o mask_for_n4.nii.gz\n"
-                before_denoise += "scil_image_math.py convert mask_for_n4.nii.gz mask_for_n4.nii.gz -f --data_type uint8\n"
-            }
-        }
+// after_denoise += "scil_apply_bias_field_on_dwi.py $image ${image.simpleName}_n4_bias_field.nii.gz n4denoised.nii.gz -f\n"
+after_denoise += "fslmaths n4denoise.nii.gz -thr 0 ${image.simpleName}__n4denoised.nii.gz\n"
 
-        """
-        export OMP_NUM_THREADS=$task.cpus
-        export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=$task.cpus
-        export OPENBLAS_NUM_THREADS=1
-        export ANTS_RANDOM_SEED=$params.random_seed
+"""
+export OMP_NUM_THREADS=$task.cpus
+export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=$task.cpus
+export OPENBLAS_NUM_THREADS=1
+export ANTS_RANDOM_SEED=$params.random_seed
+
+$before_denoise
+
+mrhardi n4 --in image4n4.nii.gz $args \
+    --out n4denoise \
+    --config $config
+$after_denoise
+"""
 
-        $before_denoise
-        mrhardi n4 $args \
-            --out n4denoise \
-            --config $config
-        $after_denoise
-        """
 }
 
 process apply_n4_bias_field {
diff --git a/modules/processes/reconstruct.nf b/modules/processes/reconstruct.nf
index 4f55eed..e5dbbfd 100644
--- a/modules/processes/reconstruct.nf
+++ b/modules/processes/reconstruct.nf
@@ -256,7 +256,7 @@ process scilpy_csd {
         mrhardi sh_order --bvals $bval --bvecs $bvec --out sh_order.txt $args
         scil_compute_ssst_fodf.py $dwi $bval $bvec $response ${sid}_fodf.nii.gz \
             --mask mask4scil.nii.gz \
-            --force_b0_threshold \
+            --skip_b0_check \
             --sh_order \$(cat sh_order.txt) \
             --processes $task.cpus
         """
@@ -294,7 +294,7 @@ process scilpy_msmt_csd {
             --vf ${sid}_vf.nii.gz \
             --vf_rgb ${sid}_vf_rgb.nii.gz \
             --mask mask4scil.nii.gz \
-            --force_b0_threshold \
+            --skip_b0_check \
             --sh_order \$(cat sh_order.txt) \
             --processes $task.cpus
         """
diff --git a/modules/processes/register.nf b/modules/processes/register.nf
index 3f93349..688c034 100644
--- a/modules/processes/register.nf
+++ b/modules/processes/register.nf
@@ -18,6 +18,7 @@ process ants_register {
 
     input:
         tuple val(sid), path(moving), path(target), val(reference), file(mask), file(metadata)
+        val(register_without_masks)
         val(caller_name)
         val(additional_publish_path)
         val(publish)
@@ -31,19 +32,47 @@ process ants_register {
         tuple val(sid), path("${moving[0].simpleName}__registration_warped_metadata.*"), optional: true, emit: metadata
     script:
         def mask_arg = ""
-        if ( !mask.iterator().inject(false) { c, i -> c || i.empty() } ) {
+        def before_cmd = ""
+        def after_cmd = ""
+        def tg = []
+        def mv = []
+        def has_masks = mask.iterator().inject(false) { c, i -> c || i.empty() }
+        if ( !(register_without_masks && has_masks) ) {
             mask_arg = "--mask ${mask.iterator().collect{ it.name }.join(',')}"
         }
-
+        if ( mask.size() == 2 ) {
+            def i = 0
+            before_cmd += "cp ${file(reference).name} reference.nii.gz\n"
+            before_cmd += "scil_volume_crop.py ${mask[0]} m.nii.gz --output_bbox bbox.nii.gz -f\n"
+            for (f in target) {
+                before_cmd += "scil_volume_crop.py $f target${i}.nii.gz --input_bbox bbox.nii.gz -f\n"
+                tg += ["target${i}.nii.gz"]
+                i += 1
+            }
+            before_cmd += "scil_volume_crop.py ${mask[1]} m.nii.gz --output_bbox bbox.nii.gz -f\n"
+            i = 0
+            for (f in moving) {
+                before_cmd += "scil_volume_crop.py $f moving${i}.nii.gz --input_bbox bbox.nii.gz -f\n"
+                mv += ["moving${i}.nii.gz"]
+                i += 1
+            }
+            after_cmd += "cp reference.nii.gz ${moving[0].simpleName}__registration_ref.nii.gz\n"
+        }
+        else {
+            tg = target
+            mv = moving
+        }
         """
         export OMP_NUM_THREADS=$task.cpus
         export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=$task.cpus
         export OPENBLAS_NUM_THREADS=1
         export ANTS_RANDOM_SEED=$params.random_seed
 
+        $before_cmd
+
         mrhardi ants_registration \
-            --moving ${moving.join(",")} \
-            --target ${target.join(",")} \
+            --moving ${mv.join(",")} \
+            --target ${tg.join(",")} \
             --out ${moving[0].simpleName}__registration $mask_arg \
             ${params.verbose_outputs ? "--verbose" : ""} \
             --config $config
@@ -88,7 +117,9 @@ process ants_register {
                 break
             fi
         done
-            
+
+        $after_cmd
+
         """
 }
 
@@ -155,6 +186,37 @@ process ants_transform {
         """
 }
 
+process identity {
+    label "FAST"
+    label "res_single_cpu"
+
+    publishDir "${params.output_root}/all/${sid}/$caller_name/${task.process.replaceAll(":", "/")}", mode: "$params.publish_all_mode", enabled: params.publish_all, overwrite: true
+    publishDir "${["${params.output_root}/${sid}", additional_publish_path].findAll({ it }).join("/")}", saveAs: { f -> ("$publish" == "true") ? f.contains("metadata") ? null : publish_suffix ? "${sid}_${publish_suffix}.nii.gz" : remove_alg_suffixes(f) : null }, mode: params.publish_mode, overwrite: true
+
+    input:
+        tuple val(sid), path(img), path(ref), file(metadata)
+        val(caller_name)
+        val(additional_publish_path)
+        val(publish)
+        val(publish_suffix)
+        val(intype)
+        val(outype)
+        val(interp)
+        path(config)
+    output:
+        tuple val(sid), path("${img.simpleName}__transformed.nii.gz"), emit: image
+        tuple val(sid), path("${img.simpleName}__transformed.bvec"), optional: true, emit: bvec
+        tuple val(sid), path("${img.simpleName}__transformed_metadata.*"), optional: true, emit: metadata
+    script:
+        """
+        export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
+        export OMP_NUM_THREADS=1
+        export OPENBLAS_NUM_THREADS=1
+        export ANTS_RANDOM_SEED=$params.random_seed
+        antsApplyTransforms -e $intype -u $outype -n $interp -t identity -i $img -r $ref -o ${img.simpleName}__transformed.nii.gz
+        """
+}
+
 process align_to_closest {
     label "ALIGN"
     input:
diff --git a/modules/processes/upsample.nf b/modules/processes/upsample.nf
index 8743635..16b9c9b 100644
--- a/modules/processes/upsample.nf
+++ b/modules/processes/upsample.nf
@@ -17,9 +17,11 @@ process scilpy_resample {
     publishDir "${["${params.output_root}/${sid}", additional_publish_path].findAll({ it }).join("/")}", saveAs: { f -> f.contains("metadata") ? null : f.contains("${mask.simpleName}") ? ("$publish_mask" == "true") ? mask_prefix ? "${sid}_${mask_prefix}.nii.gz" : remove_alg_suffixes(f) : null : remove_alg_suffixes(f) }, mode: params.publish_mode, overwrite: true
 
     input:
-        tuple val(sid), path(image), file(mask), file(metadata)
+        tuple val(sid), path(image), file(mask), file(metadata), file(reference)
         val(caller_name)
         val(interpolation)
+        val(antstype)
+        val(antsinterp)
         val(publish_mask)
         val(mask_prefix)
         val(additional_publish_path)
@@ -35,6 +37,9 @@ process scilpy_resample {
         }
         if ( !metadata.empty() )
             after_script += "mrhardi metadata --in ${image.getSimpleName()}__resampled.nii.gz --update_affine --metadata $metadata\n"
+        if ( !reference.empty() ) {
+            after_script += "antsApplyTransforms -v -d 3 -e $antstype -n $antsinterp -r $reference -i ${image.simpleName}__resampled.nii.gz -o ${image.simpleName}__resampled.nii.gz -t identity\n"
+        }
         """
         export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
         export OMP_NUM_THREADS=1
@@ -71,7 +76,7 @@ process scilpy_resample_to_reference {
     script:
         def after_script = ""
         if ( !mask.empty() ) {
-            after_script += "scil_resample_volume.py $mask mask_resampled.nii.gz --ref $reference --interp nn\n"
+            after_script += "scil_resample_volume.py $mask mask_resampled.nii.gz --ref $reference --interp nn --enforce_dimensions\n"
             after_script += "scil_image_math.py floor mask_resampled.nii.gz ${mask.getSimpleName()}__resampled.nii.gz --data_type uint8 -f\n"
         }
         if ( !metadata.empty() )
@@ -80,7 +85,7 @@ process scilpy_resample_to_reference {
         export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
         export OMP_NUM_THREADS=1
         export OPENBLAS_NUM_THREADS=1
-        scil_resample_volume.py $image resampled.nii.gz --ref $reference --interp $interpolation
+        scil_resample_volume.py $image resampled.nii.gz --ref $reference --interp $interpolation --enforce_dimensions
         fslmaths resampled.nii.gz -thr 0 ${image.simpleName}__resampled.nii.gz
         if [ "\$(mrinfo -datatype $image)" != "\$(mrinfo -datatype ${image.simpleName}__resampled.nii.gz)" ]
         then
diff --git a/modules/processes/utils.nf b/modules/processes/utils.nf
index 5821c6b..93114e6 100644
--- a/modules/processes/utils.nf
+++ b/modules/processes/utils.nf
@@ -277,8 +277,8 @@ process check_dwi_conformity {
         """
 }
 
-def get_tissue_file (tissue_class, image_list) {
-    return image_list[params.segmentation_classes.indexOf(tissue_class)]
+def get_tissue_file (tissue_class, image_list, classes) {
+    return image_list[classes.indexOf(tissue_class)]
 }
 
 def pvf_space_to_masks ( sid, brain_mask ) {
@@ -327,6 +327,7 @@ process pvf_to_mask {
 
     input:
         tuple val(sid), path(pvf_images), path(brain_mask)
+        val(classes)
         val(caller_name)
         val(additional_publish_path)
     output:
@@ -335,8 +336,8 @@ process pvf_to_mask {
         tuple val(sid), path("${sid}_3t_{${params.tissue_masks_mapping.collect{ it.key }.join(',')}}_pvf.nii.gz"), emit: pvf_3t
     script:
         def map_initializer = params.tissue_masks_mapping.collect{ 
-            it.value.size() == 1 ? "cp ${get_tissue_file(it.value[0], pvf_images)} ${sid}_3t_${it.key}_pvf.nii.gz"
-                                 : "scil_image_math.py addition ${it.value.collect{ i -> get_tissue_file(i, pvf_images) }.join(' ')} ${sid}_3t_${it.key}_pvf.nii.gz --data_type float32 -f"
+            it.value.size() == 1 ? "cp ${get_tissue_file(it.value[0], pvf_images, classes)} ${sid}_3t_${it.key}_pvf.nii.gz"
+                                 : "scil_image_math.py addition ${it.value.collect{ i -> get_tissue_file(i, pvf_images, classes) }.join(' ')} ${sid}_3t_${it.key}_pvf.nii.gz --data_type float32 -f"
         }
         def safe_wm_initializer = params.tissue_masks_mapping.collect{ 
             compute_safe_maps(sid, it.key, dilation_factor = it.key == "csf" ? params.safe_csf_mask_dilation : it.key == "gm" ? params.safe_gm_mask_dilation : false)
@@ -744,7 +745,7 @@ process validate_gradients {
         tuple val(sid), path("${bvec.simpleName}__validated.bvec"), emit: bvecs
     script:
         def args = ""
-        if ( !mask.empty() ) args += " --mask $mask"
+        // if ( !mask.empty() ) args += " --mask $mask"
         if ( !peaks_vals.empty() ) args += " --peaks_vals $peaks_vals"
         """
         export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=1
diff --git a/modules/workflows/preprocess.nf b/modules/workflows/preprocess.nf
index e09185e..93ea839 100644
--- a/modules/workflows/preprocess.nf
+++ b/modules/workflows/preprocess.nf
@@ -83,6 +83,7 @@ workflow registration_wkf {
         publish
         publish_suffix
         trans_publish_suffix
+        register_without_masks
         registration_parameters
         transformation_parameters
     main:
@@ -96,6 +97,7 @@ workflow registration_wkf {
         into_register = join_optional(into_register, mask_channel)
         ants_register(
             join_optional(into_register, reg_metadata),
+            register_without_masks,
             "preprocess",
             additional_publish_path,
             publish, publish_suffix,
@@ -185,17 +187,6 @@ workflow epi_correction_wkf {
             .join(rev_b0_metadata)
             .map{ [it[0], it[1..-1]] }
 
-        ec_align_b0_wkf(
-            b0_channel,
-            reverse_b0_channel,
-            b0_metadata,
-            rev_b0_metadata
-        )
-
-        b0_channel = ec_align_b0_wkf.out.b0
-        reverse_b0_channel = ec_align_b0_wkf.out.rev_b0
-        b0_meta_with_reverse_channel = ec_align_b0_wkf.out.metadata
-
         acq_channel = dwi_with_reverse_channel
             .map{ [it[0], [it[2]]] }
             .join(existing_rev.map{ [it[0], [it[2]]] })
@@ -216,7 +207,7 @@ workflow epi_correction_wkf {
 
         metadata_channel = ec_concatenate_b0.out.metadata
             .join(meta_with_reverse_channel)
-            .join(b0_meta_with_reverse_channel)
+            .join(b0_metadata.map{ [it[0], it[1..-1]] })
             .map{ [it[0], [it[1]] + it[2] + it[3]] }
 
         generate_b0_bval(
@@ -296,9 +287,6 @@ workflow epi_correction_wkf {
         excluded_dwi_metadata = excluded_indexes
             .join(metadata_channel)
             .map{ it[0..1] }
-        forward_transform = ec_align_b0_wkf.out.b0_transform
-        reverse_transform = ec_align_b0_wkf.out.rev_transform
-        transform_reference = ec_align_b0_wkf.out.reference
 }
 
 workflow ec_align_b0_wkf {
diff --git a/modules/workflows/segment.nf b/modules/workflows/segment.nf
index 8dcf112..05bb05c 100644
--- a/modules/workflows/segment.nf
+++ b/modules/workflows/segment.nf
@@ -23,15 +23,16 @@ include {
 include {
     resampling_reference;
     resampling_reference as resampling_reference_fa;
-    scilpy_resample_to_reference as resample_template;
-    scilpy_resample_to_reference as resample_segmentation;
-    scilpy_resample_to_reference as resample_segmentation_pre_transform;
+    scilpy_resample as resample_template;
+    scilpy_resample as resample_template_mask;
+    scilpy_resample as resample_segmentation;
+    scilpy_resample as resample_segmentation_pre_transform;
     scilpy_resample_to_reference as resample_template_fa;
     scilpy_resample_to_reference as resample_wm_atlas;
-    scilpy_resample_to_reference as resample_d99;
-    scilpy_resample_to_reference as resample_charm;
-    scilpy_resample_to_reference as resample_sarm;
-    scilpy_resample_to_reference as resample_inia19
+    scilpy_resample as resample_d99;
+    scilpy_resample as resample_charm;
+    scilpy_resample as resample_sarm;
+    scilpy_resample as resample_inia19
 } from '../processes/upsample.nf'
 include {
     ants_transform as transform_d99;
@@ -70,11 +71,10 @@ workflow segment_nmt_wkf {
         if ( template_to_t1_transform && template_resampling_reference ) {
             resample_segmentation_pre_transform(
                 segmentation_channel
-                    .join(template_resampling_reference)
-                    .map{ it + ["", ""] },
+                    .map{ it + ["", ""] }
+                    .join(template_resampling_reference),
                 "segmentation",
-                "nn",
-                false,
+                "nn", "0", "NearestNeighbor",
                 false,
                 "", ""
             )
@@ -96,7 +96,7 @@ workflow segment_nmt_wkf {
             template_mask_channel = prepend_sid_template_mask(t1_channel.map{ [it[0], file("${params.tissue_segmentation_root}/tissue_segmentation_mask_whole_no_bv.nii.gz")] })
 
             resampling_reference(
-                t1_channel.join(template_channel).map{ [it[0], it[1..-1]] },
+                t1_channel.map{ [it[0], it[1..-1]] },
                 "segmentation",
                 params.resampling_subdivision,
                 params.resampling_min_resolution,
@@ -104,22 +104,27 @@ workflow segment_nmt_wkf {
             )
             resample_template(
                 template_channel
-                    .join(resampling_reference.out.reference)
-                    .join(template_mask_channel)
-                    .map{ it + [""] },
+                    .map{ it + ["", "", ""] },
                 "segmentation",
-                "lin",
+                "lin", "", "",
                 false,
+                "", ""
+            )
+            resample_template_mask(
+                template_mask_channel
+                    .map{ it + ["", ""] }
+                    .join(resample_template.out.image),
+                "segmentation",
+                "lin", "0", "NearestNeighbor",
                 false,
                 "", ""
             )
             resample_segmentation(
                 segmentation_channel
-                    .join(resampling_reference.out.reference)
-                    .map{ it + ["", ""] },
+                    .map{ it + ["", ""] }
+                    .join(resample_template.out.image),
                 "segmentation",
-                "nn",
-                false,
+                "nn", "0", "NearestNeighbor",
                 false,
                 "", ""
             )
@@ -128,12 +133,13 @@ workflow segment_nmt_wkf {
                 t1_channel.map{ [it[0], [it[1]]] },
                 resample_template.out.image.map{ [it[0], [it[1]]] },
                 resample_segmentation.out.image.map{ [it[0], [it[1]]] },
-                mask_channel.join(resample_template.out.mask).map{ [it[0], [it[1]]] },
+                mask_channel.join(resample_template_mask.out.mask).map{ [it[0], [it[1]]] },
                 null,
                 null,
                 "segmentation",
                 false,
                 "", "",
+                false,
                 params.segmentation_registration_config,
                 params.ants_transform_segmentation_config
             )
@@ -148,6 +154,7 @@ workflow segment_nmt_wkf {
             atropos.out.vol_fractions
                 .map{ [it[0], it[1].sort{ a, b -> params.segmentation_classes.findIndexOf{ i -> a.simpleName.contains("_$i") } <=> params.segmentation_classes.findIndexOf{ i -> b.simpleName.contains("_$i") } }] }
                 .join(mask_channel),
+            params.segmentation_classes,
             "segmentation",
             "segmentation"
         )
@@ -186,11 +193,9 @@ workflow transform_atlases_wkf {
             
             resample_d99(
                 d99_channel
-                    .join(template_resampling_reference)
-                    .map{ it + ["", ""] },
+                    .map{ it + ["", "", ""] },
                 "segmentation",
-                "nn",
-                false,
+                "nn", "", "",
                 false,
                 "", ""
             )
@@ -213,11 +218,9 @@ workflow transform_atlases_wkf {
             
             resample_charm(
                 charm_channel
-                    .join(template_resampling_reference)
-                    .map{ it + ["", ""] },
+                    .map{ it + ["", "", ""] },
                 "segmentation",
-                "nn",
-                false,
+                "nn", "", "",
                 false,
                 "", ""
             )
@@ -240,11 +243,9 @@ workflow transform_atlases_wkf {
             
             resample_sarm(
                 sarm_channel
-                    .join(template_resampling_reference)
-                    .map{ it + ["", ""] },
+                    .map{ it + ["", "", ""] },
                 "segmentation",
-                "nn",
-                false,
+                "nn", "", "",
                 false,
                 "", ""
             )
@@ -267,11 +268,9 @@ workflow transform_atlases_wkf {
             
             resample_inia19(
                 inia19_channel
-                    .join(template_resampling_reference)
-                    .map{ it + ["", ""] },
+                    .map{ it + ["", "", ""] },
                 "segmentation",
-                "nn",
-                false,
+                "nn", "", "",
                 false,
                 "", ""
             )
@@ -343,6 +342,7 @@ workflow segment_wm_wkf {
             "segmentation",
             true,
             "", "",
+            false,
             params.segmentation_registration_config,
             params.ants_transform_segmentation_config
         )
diff --git a/modules/workflows/t1_registration.nf b/modules/workflows/t1_registration.nf
index 4d8e25f..0d7bae1 100644
--- a/modules/workflows/t1_registration.nf
+++ b/modules/workflows/t1_registration.nf
@@ -72,9 +72,9 @@ include {
 } from '../processes/utils.nf'
 include {
     resampling_reference;
-    scilpy_resample_to_reference as resample_template;
-    scilpy_resample_to_reference as resample_dilated_mask;
-    scilpy_resample_to_reference as resample_whole_mask
+    scilpy_resample as resample_template;
+    scilpy_resample as resample_dilated_mask;
+    scilpy_resample as resample_whole_mask
 } from '../processes/upsample.nf'
 include {
     get_data_path;
@@ -86,7 +86,8 @@ include {
     rename_sequentially as reorder_template_to_b0;
     rename_sequentially as reorder_b0_to_template;
     rename_sequentially as reorder_t1_to_b0;
-    rename_sequentially as reorder_b0_to_t1
+    rename_sequentially as reorder_b0_to_t1;
+    change_name as rename_reference_mask;
 } from "../processes/io.nf"
 
 params.use_quick = false
@@ -145,10 +146,7 @@ workflow t12b0_registration {
         }
         else {
             resampling_reference(
-                dwi_channel
-                    .map{ it[0..1] }
-                    .join(t1_channel)
-                    .join(template_channel)
+                t1_channel
                     .map{ [it[0], it[1..-1]] },
                 "preprocess",
                 params.resampling_subdivision,
@@ -160,32 +158,27 @@ workflow t12b0_registration {
 
         resample_template(
             template_channel
-                .join(registration_reference)
-                .join(template_mask_channel)
-                .map{ it + [""] },
+                .map{ it + ["", "", ""] },
             "preprocess",
-            "lin",
-            false,
+            "lin", "", "",
             false,
             "", ""
         )
         resample_dilated_mask(
             template_dilated_mask_channel
-                .join(registration_reference)
-                .map{ it + ["", ""] },
+                .map{ it + ["", ""] }
+                .join(resample_template.out.image),
             "preprocess",
-            "nn",
-            false,
+            "nn", "0", "NearestNeighbor",
             false,
             "", ""
         )
         resample_whole_mask(
             template_whole_mask_channel
-                .join(registration_reference)
-                .map{ it + ["", ""] },
+                .map{ it + ["", ""] }
+                .join(resample_template.out.image),
             "preprocess",
-            "nn",
-            false,
+            "nn", "0", "NearestNeighbor",
             false,
             "", ""
         )
@@ -216,7 +209,7 @@ workflow t12b0_registration {
             t1_to_b0_affine.out.t1,
             t1_to_b0_affine.out.dwi_mask,
             t1_to_b0_affine.out.t1_mask,
-            mask_template.out.image,
+            resample_template.out.image,
             resample_template.out.mask,
             resample_dilated_mask.out.image,
             dwi_metadata_channel,
@@ -390,6 +383,7 @@ workflow t1_to_b0_affine {
             false,
             "",
             "",
+            true,
             t1_affine_config,
             params.ants_transform_mask_config
         )
@@ -405,6 +399,7 @@ workflow t1_to_b0_affine {
             false,
             "",
             "",
+            true,
             b0_affine_config,
             params.ants_transform_mask_config
         )
@@ -539,16 +534,19 @@ workflow t1_to_b0_syn {
             false,
             "",
             "",
+            false,
             t1_syn_config,
             params.ants_transform_mask_config
         )
 
+        dwi_mask_channel = rename_reference_mask(dilated_reference_mask_channel, "dwi_mask")
+
         b0_to_reference_syn(
             reference_fixed_channel,
             b0_moving_channel,
             dwi_mask_channel,
             dilated_reference_mask_channel
-                .join(difference_masks.out.mask)
+                .join(dwi_mask_channel)
                 .map{ [it[0], it[1..-1]] },
             null,
             null,
@@ -556,6 +554,7 @@ workflow t1_to_b0_syn {
             false,
             "",
             "",
+            true,
             b0_syn_config,
             params.ants_transform_mask_config
         )
@@ -586,6 +585,7 @@ workflow t1_to_b0_syn {
             false,
             "",
             "",
+            true,
             t1_to_b0_syn_config,
             params.ants_transform_mask_config
         )
@@ -686,30 +686,34 @@ workflow t1_mask_to_b0 {
         bet_mask(extract_target_b0.out.b0, "preprocess", "false", "")
         dwi_mask_channel = bet_mask.out.mask
 
-        mask_target_b0(extract_target_b0.out.b0.join(dwi_mask_channel).map{ it + [""] }, "preprocess", "false")
+        //mask_target_b0(extract_target_b0.out.b0.join(dwi_mask_channel).map{ it + [""] }, "preprocess", "false")
 
         compute_target_pdavg(dwi_channel.map{ it[0..2] }.map{ it + ["", ""] }, "preprocess", "false")
-        mask_target_pdavg(compute_target_pdavg.out.image.join(dwi_mask_channel).map{ it + [""] }, "preprocess", "false")
+        //mask_target_pdavg(compute_target_pdavg.out.image.join(dwi_mask_channel).map{ it + [""] }, "preprocess", "false")
 
-        mask_moving_t1(t1_channel.join(t1_mask_channel).map{ it + [""] }, "preprocess", "false")
+        //mask_moving_t1(t1_channel.join(t1_mask_channel).map{ it + [""] }, "preprocess", "false")
 
-        target_channel = mask_target_b0.out.image
-            .join(mask_target_pdavg.out.image)
+        target_channel = extract_target_b0.out.b0
+            .join(compute_target_pdavg.out.image)
             .map{ [it[0], it[1..-1]] }
-        moving_channel = mask_moving_t1.out.image
+        moving_channel = t1_channel
             .map{ [it[0], [it[1]]] }
+        mask_channel = bet_mask.out.mask
+            .join(t1_mask_channel)
+            .map{ [it[0], it[1..-1]] }
 
         t1_to_b0_registration_wkf(
             target_channel,
             moving_channel,
             null,
-            null,
+            mask_channel,
             null,
             null,
             "",
             false,
             "",
             "",
+            true,
             params.t1_to_b0_registration_config,
             ""
         )
diff --git a/nextflow.config b/nextflow.config
index 6c8ed5e..8ad5a7e 100644
--- a/nextflow.config
+++ b/nextflow.config
@@ -90,7 +90,7 @@ params.use_cuda = false
 params.conservative_resources = true
 params.free_processes = 1
 params.memory_buffer_gb = 10
-params.max_cpu_per_process = null
+params.max_cpu_per_process = 10
 params.max_attempts = 3
 params.check_memory_requirements = false
 
diff --git a/workflows/preprocess.nf b/workflows/preprocess.nf
index 09af1aa..8906403 100644
--- a/workflows/preprocess.nf
+++ b/workflows/preprocess.nf
@@ -16,7 +16,8 @@ include {
     extract_b0 as dwi_b0;
     extract_b0 as extract_epi_corrected_b0;
     extract_b0 as extract_b0_preprocessed;
-    extract_b0 as extract_b0_reference
+    extract_b0 as extract_b0_reference;
+    extract_b0 as extract_b0_identity
 } from '../modules/processes/preprocess.nf'
 include {
     scil_compute_dti_fa
@@ -41,7 +42,17 @@ include {
     ants_transform as ants_transform_raw_t1_mask;
     ants_transform as apply_transform_epi_dwi;
     ants_transform as apply_transform_epi_rev;
-    ants_transform as apply_transform_epi_field
+    ants_transform as apply_transform_epi_field;
+    identity;
+    identity as identity_t1;
+    identity as identity_t1_mask;
+    identity as identity_safe_wm_mask;
+    identity as identity_wm;
+    identity as identity_wm_mask;
+    identity as identity_gm;
+    identity as identity_gm_mask;
+    identity as identity_csf;
+    identity as identity_csf_mask
 } from '../modules/processes/register.nf'
 include {
     convert_float_to_integer as convert_wm_segmentation;
@@ -81,13 +92,15 @@ include {
     normalize_inter_b0
 } from '../modules/processes/denoise.nf'
 include {
-    scilpy_resample_to_reference as resample_wm;
-    scilpy_resample_to_reference as resample_gm;
-    scilpy_resample_to_reference as resample_csf;
-    scilpy_resample_to_reference as resample_t1;
-    scilpy_resample_to_reference as resample_dwi;
-    scilpy_resample_to_reference as resample_raw_dwi;
-    scilpy_resample_to_reference as resample_raw_t1;
+    scilpy_resample as resample_wm;
+    scilpy_resample as resample_gm;
+    scilpy_resample as resample_csf;
+    scilpy_resample as resample_t1;
+    scilpy_resample as resample_dwi;
+    scilpy_resample as resample_dwi_mask;
+    scilpy_resample as resample_raw_dwi;
+    scilpy_resample as resample_raw_t1;
+    scilpy_resample as resample_t1_mask;
     resampling_reference
 } from '../modules/processes/upsample.nf'
 include {
@@ -139,6 +152,7 @@ include {
 params.gaussian_noise_correction = true
 params.gibbs_ringing_correction = true
 params.dwi_mask_from_t1_mask = true
+params.dwi_use_native_t1_mask = false
 params.epi_correction = true
 params.epi_algorithm = "topup"
 params.eddy_correction = true
@@ -152,7 +166,11 @@ params.resampling_subdivision = 2
 params.resampling_min_resolution = false
 params.force_resampling_resolution = false
 
-// T1 preprocess workflow parameters
+// Preprocess workflow parameters
+params.preprocess_dwi = true
+params.preprocess_t1 = true
+params.resample_dwi = true
+params.resample_t1 = true
 params.denoise_t1 = true
 params.nlmeans_t1 = true
 params.t1_intensity_normalization = true
@@ -194,576 +212,592 @@ workflow preprocess_wkf {
         ).map{ [it[0]] }
 
         // T1 preprocessing
-        t1_preprocess_wkf(t1_channel, t1_mask_channel)
-        t1_channel = t1_preprocess_wkf.out.t1
-        t1_mask_channel = t1_preprocess_wkf.out.mask
-
-        // Fix odd number of slices in phase direction for EPI correction
-        check_odd_dimensions(
-            dwi_channel
-                .join(rev_channel)
-                .join(dwi_mask_channel)
-                .join(collect_paths(meta_channel.join(rev_meta_channel))),
-            "preprocess"
-        )
-
-        rev_bval_bvec_channel = fill_missing_datapoints(
-            check_odd_dimensions.out.rev_bval_bvec,
-            ref_id_channel,
-            1, ["", ""]
-        )
+        if ( params.preprocess_t1 ) {
+            t1_preprocess_wkf(t1_channel, t1_mask_channel)
+            t1_channel = t1_preprocess_wkf.out.t1
+            t1_mask_channel = t1_preprocess_wkf.out.mask
+        }
 
-        dwi_channel = check_odd_dimensions.out.dwi
-        rev_channel = fill_missing_datapoints(
-            check_odd_dimensions.out.rev,
-            ref_id_channel,
-            1, [""]
-        ).join(rev_bval_bvec_channel)
-        dwi_mask_channel = fill_missing_datapoints(
-            check_odd_dimensions.out.mask,
-            ref_id_channel,
-            1, [""]
-        )
+        if (params.preprocess_dwi) {
+            // Fix odd number of slices in phase direction for EPI correction
+            check_odd_dimensions(
+                dwi_channel
+                    .join(rev_channel)
+                    .join(dwi_mask_channel)
+                    .join(collect_paths(meta_channel.join(rev_meta_channel))),
+                "preprocess"
+            )
 
-        checked_meta_channel = exclude_missing_datapoints(
-            check_odd_dimensions.out.metadata
-                .map{ it.flatten() }
-                .map{ [it[0], it[1..-1]] }
-                .transpose(),
-             1, ""
-        ).groupTuple()
-
-        meta_channel = checked_meta_channel
-            .map{ [it[0], it[1].findAll{ i -> !i.simpleName.contains("_rev") }].flatten() }
-        rev_meta_channel = fill_missing_datapoints(
-            checked_meta_channel
-                .map{ [it[0], it[1].findAll{ i -> i.simpleName.contains("_rev") }].flatten() },
-            ref_id_channel,
-            1, [""]
-        )
+            rev_bval_bvec_channel = fill_missing_datapoints(
+                check_odd_dimensions.out.rev_bval_bvec,
+                ref_id_channel,
+                1, ["", ""]
+            )
 
-        // Copy input channels for later
-        dwi_channel.tap{ raw_dwi_channel }
-        rev_channel.tap{ raw_rev_channel }
-        t1_channel.tap{ raw_t1_channel }
-        t1_mask_channel.tap{ raw_t1_mask_channel }
-        meta_channel.tap{ raw_meta_channel }
-        rev_meta_channel.tap{ raw_rev_meta_channel }
-
-        // Perform DWI and b0 denoising
-        if ( params.gaussian_noise_correction ) {
-            dwi_denoise_wkf(dwi_channel, dwi_mask_channel, meta_channel, "true")
-            dwi_channel = replace_dwi_file(dwi_channel, dwi_denoise_wkf.out.image)
-            meta_channel = dwi_denoise_wkf.out.metadata
-
-            rev_denoise_wkf(rev_channel, dwi_mask_channel, rev_meta_channel, "false")
-            rev_channel = replace_dwi_file(rev_channel, rev_denoise_wkf.out.image)
-            rev_meta_channel = fill_missing_datapoints(
-                rev_denoise_wkf.out.metadata,
+            dwi_channel = check_odd_dimensions.out.dwi
+            rev_channel = fill_missing_datapoints(
+                check_odd_dimensions.out.rev,
+                ref_id_channel,
+                1, [""]
+            ).join(rev_bval_bvec_channel)
+            dwi_mask_channel = fill_missing_datapoints(
+                check_odd_dimensions.out.mask,
                 ref_id_channel,
                 1, [""]
             )
-        }
 
-        // Perform DWI and b0 gibbs correction
-        if ( params.gibbs_ringing_correction ) {
-            dwi_gibbs_removal(dwi_channel.map{ it[0..1] }.join(meta_channel), "preprocess", "true")
-            dwi_channel = replace_dwi_file(dwi_channel, dwi_gibbs_removal.out.image)
-            meta_channel = dwi_gibbs_removal.out.metadata
-
-            rev_gibbs_removal(
-                exclude_missing_datapoints(rev_channel.map{ it[0..1] }.join(rev_meta_channel), 1, ""),
-                "preprocess", "false"
-            )
-            rev_channel = replace_dwi_file(
-                rev_channel,
-                fill_missing_datapoints(rev_gibbs_removal.out.image, ref_id_channel, 1, [""])
-            )
-            rev_meta_channel = fill_missing_datapoints(rev_gibbs_removal.out.metadata, ref_id_channel, 1, [""])
-        }
+            checked_meta_channel = exclude_missing_datapoints(
+                check_odd_dimensions.out.metadata
+                    .map{ it.flatten() }
+                    .map{ [it[0], it[1..-1]] }
+                    .transpose(),
+                1, ""
+            ).groupTuple()
 
-        // Perform DWI signal normalization between b0 volumes
-        if ( params.normalize_inter_b0 ) {
-            normalize_inter_b0(
-                dwi_channel
-                    .map{ it[0..2] }
-                    .join(rev_channel.map{ it[0..2] })
-                    .join(meta_channel)
-                    .join(rev_meta_channel),
-                "preprocess",
-                params.b0_to_b0_normalization_config
-            )
-            dwi_channel = replace_dwi_file(dwi_channel, normalize_inter_b0.out.dwi)
-            meta_channel = normalize_inter_b0.out.dwi_metadata
-            rev_channel = replace_dwi_file(
-                rev_channel,
-                fill_missing_datapoints(normalize_inter_b0.out.rev, ref_id_channel, 1, [""])
-            )
+            meta_channel = checked_meta_channel
+                .map{ [it[0], it[1].findAll{ i -> !i.simpleName.contains("_rev") }].flatten() }
             rev_meta_channel = fill_missing_datapoints(
-                normalize_inter_b0.out.rev_metadata,
+                checked_meta_channel
+                    .map{ [it[0], it[1].findAll{ i -> i.simpleName.contains("_rev") }].flatten() },
                 ref_id_channel,
                 1, [""]
             )
-        }
-
-        // Average consecutive b0 volumes just like for EPI correction
-        squash_wkf(
-            dwi_channel,
-            rev_channel,
-            meta_channel.join(rev_meta_channel),
-            ""
-        )
 
-        dwi_channel = squash_wkf.out.dwi
-        rev_channel = squash_wkf.out.rev
+            // Copy input channels for later
+            dwi_channel.tap{ raw_dwi_channel }
+            rev_channel.tap{ raw_rev_channel }
+            t1_channel.tap{ raw_t1_channel }
+            t1_mask_channel.tap{ raw_t1_mask_channel }
+            meta_channel.tap{ raw_meta_channel }
+            rev_meta_channel.tap{ raw_rev_meta_channel }
+
+            // Perform DWI and b0 denoising
+            if ( params.gaussian_noise_correction ) {
+                dwi_denoise_wkf(dwi_channel, dwi_mask_channel, meta_channel, "true")
+                dwi_channel = replace_dwi_file(dwi_channel, dwi_denoise_wkf.out.image)
+                meta_channel = dwi_denoise_wkf.out.metadata
+
+                rev_denoise_wkf(rev_channel, dwi_mask_channel, rev_meta_channel, "false")
+                rev_channel = replace_dwi_file(rev_channel, rev_denoise_wkf.out.image)
+                rev_meta_channel = fill_missing_datapoints(
+                    rev_denoise_wkf.out.metadata,
+                    ref_id_channel,
+                    1, [""]
+                )
+            }
 
-        squashed_meta_channel = exclude_missing_datapoints(
-            squash_wkf.out.metadata
-                .map{ it.flatten() }
-                .map{ [it[0], it[1..-1]] }
-                .transpose(),
-             1, ""
-        ).groupTuple()
-
-        meta_channel = squashed_meta_channel
-            .map{ [it[0], it[1].findAll{ i -> !i.simpleName.contains("_rev") }].flatten() }
-
-        rev_meta_channel = fill_missing_datapoints(
-            squashed_meta_channel
-                .map{ [it[0], it[1].findAll{ i -> i.simpleName.contains("_rev") }].flatten() },
-            ref_id_channel,
-            1, [""]
-        )
+            // Perform DWI and b0 gibbs correction
+            if ( params.gibbs_ringing_correction ) {
+                dwi_gibbs_removal(dwi_channel.map{ it[0..1] }.join(meta_channel), "preprocess", "true")
+                dwi_channel = replace_dwi_file(dwi_channel, dwi_gibbs_removal.out.image)
+                meta_channel = dwi_gibbs_removal.out.metadata
 
-        // Extract mean b0
-        dwi_b0(
-            dwi_channel
-                .map{ it[0..2] }
-                .join(meta_channel.map{ [it[0], it[1..-1]] }),
-            "preprocess",
-            "false",
-            params.extract_mean_b0_base_config
-        )
-
-        b0_channel = dwi_b0.out.b0
-        b0_metadata_channel = dwi_b0.out.metadata
+                rev_gibbs_removal(
+                    exclude_missing_datapoints(rev_channel.map{ it[0..1] }.join(rev_meta_channel), 1, ""),
+                    "preprocess", "false"
+                )
+                rev_channel = replace_dwi_file(
+                    rev_channel,
+                    fill_missing_datapoints(rev_gibbs_removal.out.image, ref_id_channel, 1, [""])
+                )
+                rev_meta_channel = fill_missing_datapoints(rev_gibbs_removal.out.metadata, ref_id_channel, 1, [""])
+            }
 
-        // EPI correction
-        ec2eddy_channel = ref_id_channel.map{ it + ["", "", []] }
-        epi_corrected_dwi_channel = dwi_channel
-        epi_corrected_meta_channel = meta_channel
-        if ( params.epi_correction ) {
-            ref_rev_id_channel = exclude_missing_datapoints(
-                raw_rev_channel, 1, ""
-            ).map{ [it[0]] }
-            excluded_id_channel = filter_datapoints(
-                raw_rev_channel, { it[1] == "" }
-            ).map{ [it[0]] }
+            // Perform DWI signal normalization between b0 volumes
+            if ( params.normalize_inter_b0 ) {
+                normalize_inter_b0(
+                    dwi_channel
+                        .map{ it[0..2] }
+                        .join(rev_channel.map{ it[0..2] })
+                        .join(meta_channel)
+                        .join(rev_meta_channel),
+                    "preprocess",
+                    params.b0_to_b0_normalization_config
+                )
+                dwi_channel = replace_dwi_file(dwi_channel, normalize_inter_b0.out.dwi)
+                meta_channel = normalize_inter_b0.out.dwi_metadata
+                rev_channel = replace_dwi_file(
+                    rev_channel,
+                    fill_missing_datapoints(normalize_inter_b0.out.rev, ref_id_channel, 1, [""])
+                )
+                rev_meta_channel = fill_missing_datapoints(
+                    normalize_inter_b0.out.rev_metadata,
+                    ref_id_channel,
+                    1, [""]
+                )
+            }
 
-            // Average consecutive b0 volumes to speed up EPI correction
-            squash_for_epi_correction_wkf(
-                ref_rev_id_channel.join(raw_dwi_channel),
-                ref_rev_id_channel.join(raw_rev_channel),
-                ref_rev_id_channel.join(raw_meta_channel).join(raw_rev_meta_channel),
+            // Average consecutive b0 volumes just like for EPI correction
+            squash_wkf(
+                dwi_channel,
+                rev_channel,
+                meta_channel.join(rev_meta_channel),
                 ""
             )
 
-            // Run EPI correction sub-workflow
-            epi_correction_wkf(
-                squash_for_epi_correction_wkf.out.dwi,
-                squash_for_epi_correction_wkf.out.rev,
-                squash_for_epi_correction_wkf.out.metadata.map{ it.flatten() }
-            )        
+            dwi_channel = squash_wkf.out.dwi
+            rev_channel = squash_wkf.out.rev
 
-            ec_input_dwi_channel = rename_ec_input_dwi(
-                collect_paths(epi_correction_wkf.out.corrected_indexes.join(dwi_channel)),
-                "ec_input_dwi"
-            ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
-
-
-            ec_input_rev_channel = rename_ec_input_rev(
-                collect_paths(epi_correction_wkf.out.corrected_indexes.join(rev_channel)),
-                "ec_input_rev"
-            ).map{ it.flatten() }.map{ it.size() == 4 ? [it[0], it[3], it[1], it[2]] : it + ["", ""] }
-
-            ec_input_dwi_meta_channel = rename_ec_input_dwi_meta(
-                epi_correction_wkf.out.in_metadata_w_epi_correction
-                    .map{ [it[0], it[1].find{m -> m.simpleName.contains("_dwi")}, "dwi__ec_input_dwi_metadata"] }
-            ).map{ it.flatten() }
+            squashed_meta_channel = exclude_missing_datapoints(
+                squash_wkf.out.metadata
+                    .map{ it.flatten() }
+                    .map{ [it[0], it[1..-1]] }
+                    .transpose(),
+                1, ""
+            ).groupTuple()
 
-            ec_input_rev_meta_channel = rename_ec_input_rev_meta(
-                epi_correction_wkf.out.in_metadata_w_epi_correction
-                    .map{ [it[0], it[1].find{m -> m.simpleName.contains("_rev")}, "rev__ec_input_rev_metadata"] }
-            ).map{ it.flatten() }
+            meta_channel = squashed_meta_channel
+                .map{ [it[0], it[1].findAll{ i -> !i.simpleName.contains("_rev") }].flatten() }
 
-            extract_b0_reference(
-                ec_input_dwi_channel.map{ it[0..2] + [""] },
-                "preprocess", "false", params.extract_mean_b0_base_config
+            rev_meta_channel = fill_missing_datapoints(
+                squashed_meta_channel
+                    .map{ [it[0], it[1].findAll{ i -> i.simpleName.contains("_rev") }].flatten() },
+                ref_id_channel,
+                1, [""]
             )
 
-            b0_reference_for_registration = extract_b0_reference.out.b0
-            apply_transform_epi_rev(
-                ec_input_rev_channel.map{ it[0..1] }
-                    .join(b0_reference_for_registration)
-                    .join(epi_correction_wkf.out.forward_transform)
-                    .join(epi_correction_wkf.out.reverse_transform)
-                    .map{ it[0..-3] + [it[-2] + it[-1]] }
-                    .map{ it + [["true", "false"], "", ""] },
+            // Extract mean b0
+            dwi_b0(
+                dwi_channel
+                    .map{ it[0..2] }
+                    .join(meta_channel.map{ [it[0], it[1..-1]] }),
                 "preprocess",
-                "",
                 "false",
-                "",
-                params.ants_transform_base_config
+                params.extract_mean_b0_base_config
             )
 
-            rev_channel = replace_dwi_file(rev_channel, apply_transform_epi_rev.out.image)
-
-            // Applied estimated susceptibility correction to DWI
-            ec2eddy_channel = Channel.empty()
-            epi_fieldmap_channel = Channel.empty()
-            epi_displacement_field_channel = Channel.empty()
-            if ( params.epi_algorithm == "topup" ) {
-                ec2eddy_channel = epi_correction_wkf.out.param
-                    .join(epi_correction_wkf.out.prefix)
-                    .join(epi_correction_wkf.out.topup.map{ [it[0], it[1..-1]] })
-                
-                // Applied estimated susceptibility correction to DWI
-                apply_topup_wkf(
-                    ec_input_dwi_channel,
-                    apply_transform_epi_rev.out.image,
-                    ec2eddy_channel,
-                    ec_input_dwi_meta_channel
-                        .join(ec_input_rev_meta_channel)
-                        .map{ [it[0], it[1..-1]] },
+            b0_channel = dwi_b0.out.b0
+            b0_metadata_channel = dwi_b0.out.metadata
+
+            // EPI correction
+            ec2eddy_channel = ref_id_channel.map{ it + ["", "", []] }
+            epi_corrected_dwi_channel = dwi_channel
+            epi_corrected_meta_channel = meta_channel
+            if ( params.epi_correction ) {
+                ref_rev_id_channel = exclude_missing_datapoints(
+                    raw_rev_channel, 1, ""
+                ).map{ [it[0]] }
+                excluded_id_channel = filter_datapoints(
+                    raw_rev_channel, { it[1] == "" }
+                ).map{ [it[0]] }
+
+                // Average consecutive b0 volumes to speed up EPI correction
+                squash_for_epi_correction_wkf(
+                    ref_rev_id_channel.join(raw_dwi_channel),
+                    ref_rev_id_channel.join(raw_rev_channel),
+                    ref_rev_id_channel.join(raw_meta_channel).join(raw_rev_meta_channel),
                     ""
                 )
 
-                epi_corrected_dwi_channel = rename_epi_corrected_dwi(
-                    collect_paths(apply_topup_wkf.out.dwi),
-                    "topup_corrected"
+                // Run EPI correction sub-workflow
+                epi_correction_wkf(
+                    squash_for_epi_correction_wkf.out.dwi,
+                    squash_for_epi_correction_wkf.out.rev,
+                    squash_for_epi_correction_wkf.out.metadata.map{ it.flatten() }
+                )        
+
+                ec_input_dwi_channel = rename_ec_input_dwi(
+                    collect_paths(epi_correction_wkf.out.corrected_indexes.join(dwi_channel)),
+                    "ec_input_dwi"
                 ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
-                epi_corrected_meta_channel = rename_epi_corrected_meta(
-                    collect_paths(apply_topup_wkf.out.metadata),
-                    "topup_corrected_metadata"
-                ).map{ it.flatten() }
-            }
-            else {
-                epi_displacement_field_channel = epi_correction_wkf.out.field
-                epi_fieldmap_channel = epi_correction_wkf.out.fieldmap
-
-                apply_transform_epi_field(
-                    epi_displacement_field_channel
-                        .join(b0_reference_for_registration)
-                        .join(epi_correction_wkf.out.forward_transform)
-                        .map{ it[0..-2] + [it[-1]] }
-                        .map{ it + [["true"], "", ""] },
-                    "preprocess",
-                    "",
-                    "false",
-                    "",
-                    params.ants_transform_base_config
-                )
 
-                apply_epi_field_wkf(
-                    ec_input_dwi_channel.map{ it[0..1] },
-                    apply_transform_epi_rev.out.image,
-                    apply_transform_epi_field.out.image,
-                    ec_input_dwi_meta_channel
-                        .map{ [it[0], it[1]] },
-                    ""
-                )
 
-                epi_corrected_dwi_channel = rename_epi_corrected_dwi(
-                    collect_paths(replace_dwi_file(ec_input_dwi_channel, apply_epi_field_wkf.out.dwi)),
-                    "epi_corrected"
-                ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
-                epi_corrected_meta_channel = rename_epi_corrected_meta(
-                    collect_paths(ec_input_dwi_meta_channel),
-                    "epi_corrected_metadata"
+                ec_input_rev_channel = rename_ec_input_rev(
+                    collect_paths(epi_correction_wkf.out.corrected_indexes.join(rev_channel)),
+                    "ec_input_rev"
+                ).map{ it.flatten() }.map{ it.size() == 4 ? [it[0], it[3], it[1], it[2]] : it + ["", ""] }
+
+                ec_input_dwi_meta_channel = rename_ec_input_dwi_meta(
+                    epi_correction_wkf.out.in_metadata_w_epi_correction
+                        .map{ [it[0], it[1].find{m -> m.simpleName.contains("_dwi")}, "dwi__ec_input_dwi_metadata"] }
                 ).map{ it.flatten() }
-            }
 
-            epi_displacement_field_channel = fill_missing_datapoints(
-                epi_displacement_field_channel,
-                ref_id_channel,
-                1, [""]
-            )
-            epi_fieldmap_channel = fill_missing_datapoints(
-                epi_fieldmap_channel,
-                ref_id_channel,
-                1, [""]
-            )
+                ec_input_rev_meta_channel = rename_ec_input_rev_meta(
+                    epi_correction_wkf.out.in_metadata_w_epi_correction
+                        .map{ [it[0], it[1].find{m -> m.simpleName.contains("_rev")}, "rev__ec_input_rev_metadata"] }
+                ).map{ it.flatten() }
 
-            ec2eddy_channel = fill_missing_datapoints(
-                ec2eddy_channel,
-                ref_id_channel,
-                1, ["", "", []]
-            )
+                extract_b0_reference(
+                    ec_input_dwi_channel.map{ it[0..2] + [""] },
+                    "preprocess", "false", params.extract_mean_b0_base_config
+                )
 
-            epi_corrected_dwi_channel = excluded_id_channel
-                .join(dwi_channel)
-                .mix(epi_corrected_dwi_channel)
-            epi_corrected_meta_channel = excluded_id_channel
-                .join(meta_channel)
-                .mix(epi_corrected_meta_channel)
+                // b0_reference_for_registration = extract_b0_reference.out.b0
+                // apply_transform_epi_rev(
+                //     ec_input_rev_channel.map{ it[0..1] }
+                //         .join(b0_reference_for_registration)
+                //         .join(epi_correction_wkf.out.forward_transform)
+                //         .join(epi_correction_wkf.out.reverse_transform)
+                //         .map{ it[0..-3] + [it[-2] + it[-1]] }
+                //         .map{ it + [["true", "false"], "", ""] },
+                //     "preprocess",
+                //     "",
+                //     "false",
+                //     "",
+                //     params.ants_transform_base_config
+                // )
+
+                // rev_channel = replace_dwi_file(rev_channel, apply_transform_epi_rev.out.image)
 
-            // Get average susceptibility corrected b0
-            extract_epi_corrected_b0(
-                epi_corrected_dwi_channel
-                    .map{ it[0..2] }
-                    .join(collect_paths(epi_corrected_meta_channel)),
-                "preprocess",
-                "true",
-                params.extract_mean_b0_base_config
-            )
+                // Applied estimated susceptibility correction to DWI
+                ec2eddy_channel = Channel.empty()
+                epi_fieldmap_channel = Channel.empty()
+                epi_displacement_field_channel = Channel.empty()
+                if ( params.epi_algorithm == "topup" ) {
+                    ec2eddy_channel = epi_correction_wkf.out.param
+                        .join(epi_correction_wkf.out.prefix)
+                        .join(epi_correction_wkf.out.topup.map{ [it[0], it[1..-1]] })
+                    
+                    // Applied estimated susceptibility correction to DWI
+                    apply_topup_wkf(
+                        ec_input_dwi_channel,
+                        rev_channel.map{ it[0..1] },
+                        ec2eddy_channel,
+                        ec_input_dwi_meta_channel
+                            .join(ec_input_rev_meta_channel)
+                            .map{ [it[0], it[1..-1]] },
+                        ""
+                    )
 
-            b0_channel = excluded_id_channel
-                .join(b0_channel)
-                .mix(extract_epi_corrected_b0.out.b0)
-            b0_metadata_channel = excluded_id_channel
-                .join(b0_metadata_channel)
-                .mix(extract_epi_corrected_b0.out.metadata)
+                    epi_corrected_dwi_channel = rename_epi_corrected_dwi(
+                        collect_paths(apply_topup_wkf.out.dwi),
+                        "topup_corrected"
+                    ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
+                    epi_corrected_meta_channel = rename_epi_corrected_meta(
+                        collect_paths(apply_topup_wkf.out.metadata),
+                        "topup_corrected_metadata"
+                    ).map{ it.flatten() }
+                }
+                else {
+                    epi_displacement_field_channel = epi_correction_wkf.out.field
+                    epi_fieldmap_channel = epi_correction_wkf.out.fieldmap
+
+                    apply_epi_field_wkf(
+                        ec_input_dwi_channel.map{ it[0..1] },
+                        rev_channel.map{ it[0..1] },
+                        epi_correction_wkf.out.field,
+                        ec_input_dwi_meta_channel
+                            .map{ [it[0], it[1]] },
+                        ""
+                    )
 
-            if ( !params.eddy_correction ) {
-                dwi_channel = epi_corrected_dwi_channel
-                meta_channel = epi_corrected_meta_channel
-            }
+                    epi_corrected_dwi_channel = rename_epi_corrected_dwi(
+                        collect_paths(replace_dwi_file(ec_input_dwi_channel, apply_epi_field_wkf.out.dwi)),
+                        "epi_corrected"
+                    ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
+                    epi_corrected_meta_channel = rename_epi_corrected_meta(
+                        collect_paths(ec_input_dwi_meta_channel),
+                        "epi_corrected_metadata"
+                    ).map{ it.flatten() }
+                }
 
-            // Apply susceptibility corrections to raw images (for comparison)
-            if ( params.raw_to_processed_space ) {
-                raw_dwi_channel = rename_transformed_raw_dwi(
-                    collect_paths(raw_dwi_channel),
-                    "raw"
-                ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
-                raw_rev_channel = rename_transformed_raw_rev(
-                    collect_paths(raw_rev_channel).filter{ it[1] },
-                    "raw"
+                epi_displacement_field_channel = fill_missing_datapoints(
+                    epi_displacement_field_channel,
+                    ref_id_channel,
+                    1, [""]
+                )
+                epi_fieldmap_channel = fill_missing_datapoints(
+                    epi_fieldmap_channel,
+                    ref_id_channel,
+                    1, [""]
                 )
-                raw_rev_channel = raw_rev_channel
-                    .map{ it.flatten() }
-                    .map{ it.size() == 4 ? [it[0], it[3], it[1], it[2]] : it + ["", ""] }
 
-                raw_meta_channel = excluded_id_channel.join(raw_meta_channel)
-                    .mix(
-                        rename_transformed_raw_metadata(
-                            collect_paths(ec_input_dwi_meta_channel).filter{ it[1] },
-                            "raw_metadata"
-                        ).map{ it.flatten() }
-                    )
+                ec2eddy_channel = fill_missing_datapoints(
+                    ec2eddy_channel,
+                    ref_id_channel,
+                    1, ["", "", []]
+                )
 
-                raw_rev_meta_channel = excluded_id_channel.join(raw_meta_channel)
-                    .mix(
-                        rename_transformed_raw_rev_metadata(
-                            collect_paths(ec_input_rev_meta_channel),
-                            "raw_metadata"
-                        ).map{ it.flatten() }
-                    )
+                epi_corrected_dwi_channel = excluded_id_channel
+                    .join(dwi_channel)
+                    .mix(epi_corrected_dwi_channel)
+                epi_corrected_meta_channel = excluded_id_channel
+                    .join(meta_channel)
+                    .mix(epi_corrected_meta_channel)
 
-                if ( params.epi_algorithm == "topup" ) {
-                    raw_apply_topup_wkf(
-                        epi_correction_wkf.out.corrected_indexes.join(raw_dwi_channel),
-                        epi_correction_wkf.out.corrected_indexes
-                            .join(raw_rev_channel)
-                            .map{ it[0..1] },
-                        ec2eddy_channel,
-                        collect_paths(raw_meta_channel.join(raw_rev_meta_channel)),
-                        "raw"
-                    )
+                // Get average susceptibility corrected b0
+                extract_epi_corrected_b0(
+                    epi_corrected_dwi_channel
+                        .map{ it[0..2] }
+                        .join(collect_paths(epi_corrected_meta_channel)),
+                    "preprocess",
+                    "true",
+                    params.extract_mean_b0_base_config
+                )
 
-                    raw_dwi_channel = excluded_id_channel
-                        .join(raw_dwi_channel)
-                        .mix(raw_apply_topup_wkf.out.dwi)
+                b0_channel = excluded_id_channel
+                    .join(b0_channel)
+                    .mix(extract_epi_corrected_b0.out.b0)
+                b0_metadata_channel = excluded_id_channel
+                    .join(b0_metadata_channel)
+                    .mix(extract_epi_corrected_b0.out.metadata)
 
-                    raw_meta_channel = excluded_id_channel
-                        .join(raw_meta_channel)
-                        .mix(raw_apply_topup_wkf.out.metadata)
+                if ( !params.eddy_correction ) {
+                    dwi_channel = epi_corrected_dwi_channel
+                    meta_channel = epi_corrected_meta_channel
                 }
-                else {
-                    raw_apply_epi_field_wkf(
-                        epi_correction_wkf.out.corrected_indexes.join(raw_dwi_channel),
-                        epi_correction_wkf.out.corrected_indexes
-                            .join(raw_rev_channel)
-                            .map{ it[0..1] },
-                        epi_displacement_field_channel,
-                        collect_paths(raw_meta_channel.join(raw_rev_meta_channel)),
+
+                // Apply susceptibility corrections to raw images (for comparison)
+                if ( params.raw_to_processed_space ) {
+                    raw_dwi_channel = rename_transformed_raw_dwi(
+                        collect_paths(raw_dwi_channel),
+                        "raw"
+                    ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
+                    raw_rev_channel = rename_transformed_raw_rev(
+                        collect_paths(raw_rev_channel).filter{ it[1] },
                         "raw"
                     )
-
-                    raw_dwi_channel = excluded_id_channel
-                        .join(raw_dwi_channel)
-                        .mix(raw_apply_epi_field_wkf.out.dwi)
-
-                    raw_meta_channel = excluded_id_channel
-                        .join(raw_meta_channel)
-                        .mix(raw_apply_epi_field_wkf.out.metadata)
+                    raw_rev_channel = raw_rev_channel
+                        .map{ it.flatten() }
+                        .map{ it.size() == 4 ? [it[0], it[3], it[1], it[2]] : it + ["", ""] }
+
+                    raw_meta_channel = excluded_id_channel.join(raw_meta_channel)
+                        .mix(
+                            rename_transformed_raw_metadata(
+                                collect_paths(ec_input_dwi_meta_channel).filter{ it[1] },
+                                "raw_metadata"
+                            ).map{ it.flatten() }
+                        )
+
+                    raw_rev_meta_channel = excluded_id_channel.join(raw_meta_channel)
+                        .mix(
+                            rename_transformed_raw_rev_metadata(
+                                collect_paths(ec_input_rev_meta_channel),
+                                "raw_metadata"
+                            ).map{ it.flatten() }
+                        )
+
+                    if ( params.epi_algorithm == "topup" ) {
+                        raw_apply_topup_wkf(
+                            epi_correction_wkf.out.corrected_indexes.join(raw_dwi_channel),
+                            epi_correction_wkf.out.corrected_indexes
+                                .join(raw_rev_channel)
+                                .map{ it[0..1] },
+                            ec2eddy_channel,
+                            collect_paths(raw_meta_channel.join(raw_rev_meta_channel)),
+                            "raw"
+                        )
+
+                        raw_dwi_channel = excluded_id_channel
+                            .join(raw_dwi_channel)
+                            .mix(raw_apply_topup_wkf.out.dwi)
+
+                        raw_meta_channel = excluded_id_channel
+                            .join(raw_meta_channel)
+                            .mix(raw_apply_topup_wkf.out.metadata)
+                    }
+                    else {
+                        raw_apply_epi_field_wkf(
+                            epi_correction_wkf.out.corrected_indexes.join(raw_dwi_channel),
+                            epi_correction_wkf.out.corrected_indexes
+                                .join(raw_rev_channel)
+                                .map{ it[0..1] },
+                            epi_displacement_field_channel,
+                            collect_paths(raw_meta_channel.join(raw_rev_meta_channel)),
+                            "raw"
+                        )
+
+                        raw_dwi_channel = excluded_id_channel
+                            .join(raw_dwi_channel)
+                            .mix(raw_apply_epi_field_wkf.out.dwi)
+
+                        raw_meta_channel = excluded_id_channel
+                            .join(raw_meta_channel)
+                            .mix(raw_apply_epi_field_wkf.out.metadata)
+                    }
                 }
             }
-        }
-
-        empty_dwi_mask_id_channel = filter_datapoints(
-            dwi_mask_channel,
-            { it[1] == "" }
-        ).map{ [it[0]] }
-        dwi_mask_channel = exclude_missing_datapoints(dwi_mask_channel, 1, "")
-
-        // Compute brain mask for the DWI (when missing)
-        bet_mask(
-            empty_dwi_mask_id_channel.join(b0_channel),
-            "preprocess",
-            "${!params.dwi_mask_from_t1_mask}",
-            "dwi_mask"
-        )
-        dwi_mask_channel = dwi_mask_channel.mix(bet_mask.out.mask)
 
-        // Get better mask for the DWI from the T1 (when missing and if present)
-        if ( params.dwi_mask_from_t1_mask ) {
-            existing_t1_mask_id_channel = exclude_missing_datapoints(
-                t1_mask_channel,
-                1, ""
-            ).map{ [it[0]] }
-            absent_t1_mask_id_channel = filter_datapoints(
-                t1_mask_channel,
+            empty_dwi_mask_id_channel = filter_datapoints(
+                dwi_mask_channel,
                 { it[1] == "" }
             ).map{ [it[0]] }
+            dwi_mask_channel = exclude_missing_datapoints(dwi_mask_channel, 1, "")
 
-            n4_denoise_t1_to_b0_wkf(
-                existing_t1_mask_id_channel.join(epi_corrected_dwi_channel.map{ it[0..1] }),
-                existing_t1_mask_id_channel.join(b0_channel),
-                existing_t1_mask_id_channel.join(dwi_mask_channel),
-                existing_t1_mask_id_channel.join(epi_corrected_meta_channel),
-                params.dwi_n4_normalization_quick_config,
-                false
+            // Compute brain mask for the DWI (when missing)
+            bet_mask(
+                empty_dwi_mask_id_channel.join(b0_channel),
+                "preprocess",
+                "${!params.dwi_mask_from_t1_mask}",
+                "dwi_mask"
             )
+            dwi_mask_channel = dwi_mask_channel.mix(bet_mask.out.mask)
+
+            // Get better mask for the DWI from the T1 (when missing and if present)
+            if ( params.dwi_mask_from_t1_mask ) {
+                if ( params.dwi_use_native_t1_mask )
+                {
+                    identity(
+                        t1_mask_channel
+                            .join(b0_channel)
+                            .map{ it + [""] },
+                        "preprocess",
+                        "", false, "dwi_mask",
+                        0, "uchar", "NearestNeighbor",
+                        params.ants_transform_base_config
+                    )
+                    dwi_mask_channel = identity.out.image
+                }
+                else if ( !params.register_t1_to_dwi ) {
+                    existing_t1_mask_id_channel = exclude_missing_datapoints(
+                        t1_mask_channel,
+                        1, ""
+                    ).map{ [it[0]] }
+                    absent_t1_mask_id_channel = filter_datapoints(
+                        t1_mask_channel,
+                        { it[1] == "" }
+                    ).map{ [it[0]] }
+
+                    n4_denoise_t1_to_b0_wkf(
+                        existing_t1_mask_id_channel.join(epi_corrected_dwi_channel.map{ it[0..1] }),
+                        existing_t1_mask_id_channel.join(b0_channel),
+                        Channel.empty(),
+                        existing_t1_mask_id_channel.join(epi_corrected_meta_channel),
+                        params.dwi_n4_normalization_quick_config,
+                        false
+                    )
 
-            t1_mask_to_b0(
-                replace_dwi_file(epi_corrected_dwi_channel, n4_denoise_t1_to_b0_wkf.out.image),
-                existing_t1_mask_id_channel.join(t1_channel),
-                existing_t1_mask_id_channel.join(t1_mask_channel),
-                "false"
-            )
+                    t1_mask_to_b0(
+                        replace_dwi_file(epi_corrected_dwi_channel, n4_denoise_t1_to_b0_wkf.out.image),
+                        existing_t1_mask_id_channel.join(t1_channel),
+                        existing_t1_mask_id_channel.join(t1_mask_channel),
+                        "false"
+                    )
 
-            t1_mask_convert_datatype(
-                t1_mask_to_b0.out.mask,
-                "uint8", "preprocess",
-                !params.register_t1_to_dwi,
-                "dwi_mask", ""
-            )
+                    t1_mask_convert_datatype(
+                        t1_mask_to_b0.out.mask,
+                        "uint8", "preprocess",
+                        !params.register_t1_to_dwi,
+                        "dwi_mask", ""
+                    )
 
-            dwi_mask_channel = t1_mask_convert_datatype.out.image
-                .mix(absent_t1_mask_id_channel.join(dwi_mask_channel))
-        }
+                    dwi_mask_channel = t1_mask_convert_datatype.out.image
+                        .mix(absent_t1_mask_id_channel.join(dwi_mask_channel))
+                }
+            }
 
-        // Perform Eddy currents and motion correction on the DWI
-        if ( params.eddy_correction ) {
-            dwi_channel = rename_dwi_for_eddy(
-                collect_paths(dwi_channel),
-                "to_eddy"
-            ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
+            // Perform Eddy currents and motion correction on the DWI
+            if ( params.eddy_correction ) {
+                dwi_channel = rename_dwi_for_eddy(
+                    collect_paths(dwi_channel),
+                    "to_eddy"
+                ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
 
-            rev_channel = rename_rev_for_eddy(
-                collect_paths(rev_channel).filter{ it[1] },
-                "to_eddy"
-            )
-            rev_channel = rev_channel
-                .map{ it.flatten() }
-                .map{ it.size() == 4 ? [it[0], it[3], it[1], it[2]] : it + ["", ""] }
+                rev_channel = rename_rev_for_eddy(
+                    collect_paths(rev_channel).filter{ it[1] },
+                    "to_eddy"
+                )
+                rev_channel = rev_channel
+                    .map{ it.flatten() }
+                    .map{ it.size() == 4 ? [it[0], it[3], it[1], it[2]] : it + ["", ""] }
 
-            rev_channel = fill_missing_datapoints(
-                rev_channel,
-                ref_id_channel,
-                1, ["", "", ""]
-            )
+                rev_channel = fill_missing_datapoints(
+                    rev_channel,
+                    ref_id_channel,
+                    1, ["", "", ""]
+                )
 
-            meta_channel = rename_dwi_metadata_for_eddy(
-                collect_paths(meta_channel).filter{ it[1] },
-                "to_eddy_metadata"
-            ).map{ it.flatten() }
+                meta_channel = rename_dwi_metadata_for_eddy(
+                    collect_paths(meta_channel).filter{ it[1] },
+                    "to_eddy_metadata"
+                ).map{ it.flatten() }
 
-            rev_meta_channel = rename_rev_metadata_for_eddy(
-                collect_paths(rev_meta_channel).filter{ it[1] },
-                "to_eddy_metadata"
-            ).map{ it.flatten() }
+                rev_meta_channel = rename_rev_metadata_for_eddy(
+                    collect_paths(rev_meta_channel).filter{ it[1] },
+                    "to_eddy_metadata"
+                ).map{ it.flatten() }
 
-            rev_meta_channel = fill_missing_datapoints(
-                rev_meta_channel,
-                ref_id_channel,
-                1, [""]
-            )
+                rev_meta_channel = fill_missing_datapoints(
+                    rev_meta_channel,
+                    ref_id_channel,
+                    1, [""]
+                )
 
-            // Run Eddy sub-workflow
-            eddy_wkf(
-                dwi_channel,
-                dwi_mask_channel,
-                ec2eddy_channel,
-                epi_fieldmap_channel,
-                epi_displacement_field_channel,
-                rev_channel,
-                meta_channel.join(rev_meta_channel)
-            )
+                // Run Eddy sub-workflow
+                eddy_wkf(
+                    dwi_channel,
+                    dwi_mask_channel,
+                    ec2eddy_channel,
+                    epi_fieldmap_channel,
+                    epi_displacement_field_channel,
+                    rev_channel,
+                    meta_channel.join(rev_meta_channel)
+                )
 
-            dwi_channel = eddy_wkf.out.dwi
-                .join(eddy_wkf.out.bval)
-                .join(eddy_wkf.out.bvec)
-            meta_channel = eddy_wkf.out.metadata
-        }
+                dwi_channel = eddy_wkf.out.dwi
+                    .join(eddy_wkf.out.bval)
+                    .join(eddy_wkf.out.bvec)
+                meta_channel = eddy_wkf.out.metadata
+            }
 
-        // Perform intensity normalization on the DWI
-        if ( params.dwi_intensity_normalization ) {
-            // Run N4 sub-workflow
-            n4_denoise_wkf(
-                dwi_channel.map{ it[0..1] },
-                b0_channel,
-                dwi_mask_channel,
-                meta_channel,
-                params.dwi_n4_normalization_config,
-                true
-            )
+            // Perform intensity normalization on the DWI
+            if ( params.dwi_intensity_normalization ) {
+                // Run N4 sub-workflow
+                n4_denoise_wkf(
+                    dwi_channel.map{ it[0..1] },
+                    b0_channel,
+                    Channel.empty(),
+                    meta_channel,
+                    params.dwi_n4_normalization_config,
+                    true
+                )
 
-            dwi_channel = replace_dwi_file(dwi_channel, n4_denoise_wkf.out.image)
-            meta_channel = n4_denoise_wkf.out.metadata
-        }
+                dwi_channel = replace_dwi_file(dwi_channel, n4_denoise_wkf.out.image)
+                meta_channel = n4_denoise_wkf.out.metadata
+            }
 
-        absent_t1_mask_id_channel = filter_datapoints(
-            t1_mask_channel,
-            { it[1] == "" }
-        ).map{ [it[0]] }
-        existing_t1_mask_id_channel = exclude_missing_datapoints(
-            t1_mask_channel,
-            1, ""
-        ).map{ [it[0]] }
+            absent_t1_mask_id_channel = filter_datapoints(
+                t1_mask_channel,
+                { it[1] == "" }
+            ).map{ [it[0]] }
+            existing_t1_mask_id_channel = exclude_missing_datapoints(
+                t1_mask_channel,
+                1, ""
+            ).map{ [it[0]] }
 
-        // Get DWI mask in T1 space (for missing T1 masks)
-        dwi_mask_registration_wkf(
-            absent_t1_mask_id_channel.join(t1_channel).map{ [it[0], [it[1]]] },
-            absent_t1_mask_id_channel.join(b0_channel).map{ [it[0], [it[1]]] },
-            absent_t1_mask_id_channel.join(dwi_mask_channel),
-            null,
-            null,
-            absent_t1_mask_id_channel
-                .join(b0_metadata_channel)
-                .map{ it[0..1] + [""] },
-            "",
-            false, "", "",
-            params.b02t1_mask_registration_config,
-            null
-        )
+            // Get DWI mask in T1 space (for missing T1 masks)
+            dwi_mask_registration_wkf(
+                absent_t1_mask_id_channel.join(t1_channel).map{ [it[0], [it[1]]] },
+                absent_t1_mask_id_channel.join(b0_channel).map{ [it[0], [it[1]]] },
+                absent_t1_mask_id_channel.join(dwi_mask_channel),
+                null,
+                null,
+                absent_t1_mask_id_channel
+                    .join(b0_metadata_channel)
+                    .map{ it[0..1] + [""] },
+                "",
+                false, "", "",
+                true,
+                params.b02t1_mask_registration_config,
+                null
+            )
 
-        dwi_mask_convert_datatype(
-            dwi_mask_registration_wkf.out.image,
-            "uint8", "preprocess",
-            true,
-            "t1_mask", ""
-        )
+            dwi_mask_convert_datatype(
+                dwi_mask_registration_wkf.out.image,
+                "uint8", "preprocess",
+                true,
+                "t1_mask", ""
+            )
 
-        t1_mask_channel = existing_t1_mask_id_channel
-            .join(t1_mask_channel)
-            .mix(dwi_mask_convert_datatype.out.image)
-        raw_t1_mask_channel = t1_mask_channel
-        raw_dwi_mask_channel = dwi_mask_channel
+            t1_mask_channel = existing_t1_mask_id_channel
+                .join(t1_mask_channel)
+                .mix(dwi_mask_convert_datatype.out.image)
+            raw_t1_mask_channel = t1_mask_channel
+            raw_dwi_mask_channel = dwi_mask_channel
+        }
+        else {
+            // Copy input channels for later
+            dwi_channel.tap{ raw_dwi_channel }
+            rev_channel.tap{ raw_rev_channel }
+            t1_channel.tap{ raw_t1_channel }
+            t1_mask_channel.tap{ raw_t1_mask_channel }
+            meta_channel.tap{ raw_meta_channel }
+            rev_meta_channel.tap{ raw_rev_meta_channel }
+        }
 
         // Compute best resampling reference
         resampling_reference(
-            collect_paths(dwi_channel.map{ it[0..1] }.join(t1_channel)),
+            collect_paths(t1_channel),
             "preprocess",
             params.resample_data ? params.resampling_subdivision : "1",
             params.resample_data ? params.resampling_min_resolution : "",
@@ -774,81 +808,102 @@ workflow preprocess_wkf {
         pvf_to_resample_channel = pvf_channel.filter{ !it[1].isEmpty() }
 
         // Resample all volumes
-        resample_dwi(
-            dwi_channel
-                .map{ it[0..1] }
-                .join(reference_channel)
-                .join(dwi_mask_channel)
-                .join(meta_channel),
-            "preprocess", "lin",
-            true, true,
-            "dwi_mask", ""
-        )
+        if (params.resample_dwi) {
+            resample_dwi(
+                dwi_channel
+                    .map{ it[0..1] + [""] }
+                    .join(meta_channel)
+                    .map{ it + [""] },
+                "preprocess", "lin", 3, "",
+                false, "", ""
+            )
 
-        resample_t1(
-            t1_channel
-                .join(reference_channel)
-                .join(t1_mask_channel)
-                .map{ it + [""] },
-            "preprocess", "lin",
-            true, true,
-            "t1_mask", ""
-        )
+            dwi_channel = replace_dwi_file(dwi_channel, resample_dwi.out.image)
+            meta_channel = resample_dwi.out.metadata
 
-        resample_wm(
-            pvf_to_resample_channel
-                .map{ [it[0], it[1][0]] }
-                .join(reference_channel)
-                .join(t1_mask_channel)
-                .map{ it + [""] },
-            "preprocess", "nn",
-            true, false,
-            "", "segmentation"
-        )
-        resample_gm(
-            pvf_to_resample_channel
-                .map{ [it[0], it[1][1]] }
-                .join(reference_channel)
-                .join(t1_mask_channel)
-                .map{ it + [""] },
-            "preprocess", "nn",
-            true, false,
-            "", "segmentation"
-        )
-        resample_csf(
-            pvf_to_resample_channel
-                .map{ [it[0], it[1][2]] }
-                .join(reference_channel)
-                .join(t1_mask_channel)
-                .map{ it + [""] },
-            "preprocess", "nn",
-            true, false,
-            "", "segmentation"
-        )
+            if ( params.dwi_mask_from_t1_mask ){
+                resample_dwi_mask(
+                    dwi_mask_channel
+                        .map{ it + ["", "", ""] },
+                    "preprocess", "lin", "", "NearestNeighbor",
+                    false, "", ""
+                )
+                dwi_mask_channel = resample_dwi_mask.out.image
+            }
+        }
 
-        dwi_channel = replace_dwi_file(dwi_channel, resample_dwi.out.image)
-        dwi_mask_channel = resample_dwi.out.mask
-        meta_channel = resample_dwi.out.metadata
-        t1_channel = resample_t1.out.image
-        t1_mask_channel = resample_t1.out.mask
+        if ( params.resample_t1 ) {
+            resample_t1(
+                t1_channel
+                    .map{ it + ["", "", ""] },
+                "preprocess", "lin", "", "",
+                false,
+                "", ""
+            )
 
-        pvf_channel = resample_wm.out.image
-            .join(resample_gm.out.image)
-            .join(resample_csf.out.image)
-            .map{ [it[0], it[1..-1]] }
-            .mix(pvf_channel.filter{ it[1].isEmpty() })
+            resample_t1_mask(
+                t1_mask_channel
+                    .map{ it + ["", ""] }
+                    .join(resample_t1.out.image),
+                "preprocess", "nn", "0", "NearestNeighbor",
+                false,
+                "", ""
+            )
 
-        if ( params.raw_to_processed_space ) {
-            resample_raw_dwi(
-                raw_dwi_channel
-                    .map{ it[0..1] }
-                    .join(reference_channel)
-                    .join(raw_dwi_mask_channel)
-                    .join(raw_meta_channel),
-                "preprocess", "lin",
-                true, true,
-                "dwi_mask", "raw"
+            resample_wm(
+                pvf_to_resample_channel
+                    .map{ [it[0], it[1][0]] }
+                    .map{ it + ["", ""] }
+                    .join(resample_t1.out.image),
+                "preprocess", "nn", "0", "NearestNeighbor",
+                true,
+                "", "segmentation"
+            )
+            resample_gm(
+                pvf_to_resample_channel
+                    .map{ [it[0], it[1][1]] }
+                    .map{ it + ["", ""] }
+                    .join(resample_t1.out.image),
+                "preprocess", "nn", "0", "NearestNeighbor",
+                true,
+                "", "segmentation"
             )
+            resample_csf(
+                pvf_to_resample_channel
+                    .map{ [it[0], it[1][2]] }
+                    .map{ it + ["", ""] }
+                    .join(resample_t1.out.image),
+                "preprocess", "nn", "0", "NearestNeighbor",
+                true,
+                "", "segmentation"
+            )
+
+            t1_channel = resample_t1.out.image
+            t1_mask_channel = resample_t1_mask.out.image
+
+            pvf_channel = resample_wm.out.image
+                .join(resample_gm.out.image)
+                .join(resample_csf.out.image)
+                .map{ [it[0], it[1..-1]] }
+                .mix(pvf_channel.filter{ it[1].isEmpty() })
+        }
+
+        if ( params.raw_to_processed_space ) {
+            if (params.resample_dwi) {
+                resample_raw_dwi(
+                    raw_dwi_channel
+                        .map{ it[0..1] }
+                        .join(reference_channel)
+                        .join(raw_dwi_mask_channel)
+                        .join(raw_meta_channel),
+                    "preprocess", "lin",
+                    true, true,
+                    "dwi_mask", "raw"
+                )
+                raw_dwi_channel = replace_dwi_file(raw_dwi_channel, resample_raw_dwi.out.image)
+                raw_meta_channel = resample_raw_dwi.out.metadata
+                raw_dwi_mask_channel = resample_raw_dwi.out.mask
+            }
 
             resample_raw_t1(
                 raw_t1_channel
@@ -860,18 +915,16 @@ workflow preprocess_wkf {
                 "t1_mask", "raw"
             )
 
-            raw_dwi_channel = replace_dwi_file(raw_dwi_channel, resample_raw_dwi.out.image)
-            raw_meta_channel = resample_raw_dwi.out.metadata
             raw_t1_channel = resample_raw_t1.out.image
             raw_t1_mask_channel = resample_raw_t1.out.mask
-            raw_dwi_mask_channel = resample_raw_dwi.out.mask
         }
 
         // Get tissue masks from PVF
         pvf_to_mask(
              pvf_channel
                 .filter{ !it[1].isEmpty() }
-                .join(dwi_mask_channel),
+                .join(t1_mask_channel),
+            params.segmentation_classes,
             "preprocess",
             "segmentation"
         )
@@ -1026,6 +1079,130 @@ workflow preprocess_wkf {
                 raw_t1_mask_channel = ants_transform_raw_t1_mask.out.image
             }
         }
+        else
+        {
+            extract_b0_identity(
+                dwi_channel
+                    .map{ it[0..2] }
+                    .join(meta_channel),
+                "preprocess", "false",
+                params.extract_mean_b0_base_config
+            )
+            identity_t1(
+                t1_channel
+                    .join(extract_b0_identity.out.b0)
+                    .map{ it + [""] },
+                "preprocess",
+                "", false, "t1",
+                0, "float", "Linear",
+                params.ants_transform_base_config
+            )
+            identity_t1_mask(
+                t1_mask_channel
+                    .join(extract_b0_identity.out.b0)
+                    .map{ it + [""] },
+                "preprocess",
+                "", false, "t1_mask",
+                0, "uchar", "NearestNeighbor",
+                params.ants_transform_base_config
+            )
+            identity_wm(
+                pvf_channel
+                    .filter{ it[1] }
+                    .map{ [it[0], it[1][0]] }
+                    .join(extract_b0_identity.out.b0)
+                    .map{ it + [""] },
+                "preprocess",
+                "", false, "wm",
+                0, "float", "Linear",
+                params.ants_transform_base_config
+            )
+            identity_safe_wm_mask(
+                safe_wm_mask_channel
+                    .filter{ it[1] }
+                    .join(extract_b0_identity.out.b0)
+                    .map{ it + [""] },
+                "preprocess",
+                "", false, "wm_mask",
+                0, "uchar", "NearestNeighbor",
+                params.ants_transform_base_config
+            )
+            identity_wm_mask(
+                tissue_mask_channel
+                    .filter{ it[1] }
+                    .map{ [it[0], it[1][0]] }
+                    .join(extract_b0_identity.out.b0)
+                    .map{ it + [""] },
+                "preprocess",
+                "", false, "wm_mask",
+                0, "uchar", "NearestNeighbor",
+                params.ants_transform_base_config
+            )
+            identity_gm(
+                pvf_channel
+                    .filter{ it[1] }
+                    .map{ [it[0], it[1][1]] }
+                    .join(extract_b0_identity.out.b0)
+                    .map{ it + [""] },
+                "preprocess",
+                "", false, "gm",
+                0, "float", "Linear",
+                params.ants_transform_base_config
+            )
+            identity_gm_mask(
+                tissue_mask_channel
+                    .filter{ it[1] }
+                    .map{ [it[0], it[1][1]] }
+                    .join(extract_b0_identity.out.b0)
+                    .map{ it + [""] },
+                "preprocess",
+                "", false, "gm_mask",
+                0, "uchar", "NearestNeighbor",
+                params.ants_transform_base_config
+            )
+            identity_csf(
+                pvf_channel
+                    .filter{ it[1] }
+                    .map{ [it[0], it[1][2]] }
+                    .join(extract_b0_identity.out.b0)
+                    .map{ it + [""] },
+                "preprocess",
+                "", false, "csf",
+                0, "float", "Linear",
+                params.ants_transform_base_config
+            )
+            identity_csf_mask(
+                tissue_mask_channel
+                    .filter{ it[1] }
+                    .map{ [it[0], it[1][2]] }
+                    .join(extract_b0_identity.out.b0)
+                    .map{ it + [""] },
+                "preprocess",
+                "", false, "csf_mask",
+                0, "uchar", "NearestNeighbor",
+                params.ants_transform_base_config
+            )
+
+            t1_channel = identity_t1.out.image
+            t1_mask_channel = identity_t1_mask.out.image
+
+            pvf_channel = collect_paths(
+                identity_wm.out.image
+                    .join(identity_gm.out.image)
+                    .join(identity_csf.out.image)
+            ).mix(pvf_channel.filter{ it[1].isEmpty() })
+
+            tissue_mask_channel = collect_paths(
+                identity_wm_mask.out.image
+                    .join(identity_gm_mask.out.image)
+                    .join(identity_csf_mask.out.image)
+            ).mix(pvf_channel.filter{ it[1].isEmpty() })
+
+            safe_wm_mask_channel = pvf_channel
+                .filter{ it[1].isEmpty() }
+                .map{ [it[0], ""] }
+                .mix(identity_safe_wm_mask.out.image)
+        }
 
         // Generate tissue segmentation from T1
         d99_atlas_channel = Channel.empty()
@@ -1093,7 +1270,7 @@ workflow preprocess_wkf {
 
         crop_wm(
             pvf_to_crop_channel
-                .map{ [it[0], it[1][2]] }
+                .map{ [it[0], it[1][0]] }
                 .join(t1_mask_channel)
                 .join(t1_bbox_channel)
                 .map{ it + [""] },
@@ -1109,7 +1286,7 @@ workflow preprocess_wkf {
         )
         crop_csf(
             pvf_to_crop_channel
-                .map{ [it[0], it[1][0]] }
+                .map{ [it[0], it[1][2]] }
                 .join(t1_mask_channel)
                 .join(t1_bbox_channel)
                 .map{ it + [""] },
@@ -1241,23 +1418,26 @@ workflow preprocess_wkf {
             params.extract_mean_b0_base_config
         )
 
-        validate_gradients_wkf(dwi_channel, dwi_mask_channel)
+        if ( params.preprocess_dwi ) {
+            validate_gradients_wkf(dwi_channel, dwi_mask_channel)
+            dwi_channel = rename_processed_dwi(
+                collect_paths(validate_gradients_wkf.out.dwi),
+                "dwi_preprocessed"
+            ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
+
+            meta_channel = rename_processed_dwi_metadata(
+                collect_paths(meta_channel),
+                "dwi_preprocessed_metadata"
+            ).map{ it.flatten() }
+            dwi_mask_channel = rename_processed_dwi_mask(
+                collect_paths(dwi_mask_channel),
+                "mask_preprocessed"
+            ).map{ it.flatten() }
+        }
+
         mask_t1(t1_channel.join(t1_mask_channel).map{ it + [""] }, "preprocess", true)
         t1_channel = mask_t1.out.image
 
-        dwi_channel = rename_processed_dwi(
-            collect_paths(validate_gradients_wkf.out.dwi),
-            "dwi_preprocessed"
-        ).map{ [it[0], it[1][2], it[1][0], it[1][1]] }
-        meta_channel = rename_processed_dwi_metadata(
-            collect_paths(meta_channel),
-            "dwi_preprocessed_metadata"
-        ).map{ it.flatten() }
-        dwi_mask_channel = rename_processed_dwi_mask(
-            collect_paths(dwi_mask_channel),
-            "mask_preprocessed"
-        ).map{ it.flatten() }
-
         wm_segmentation_channel = Channel.empty()
         if ( params.generate_wm_segmentation ) {
             segment_wm_wkf(dwi_channel, dwi_mask_channel)