diff --git a/rabies/analysis_pkg/utils.py b/rabies/analysis_pkg/utils.py
index 2494370..b9235fb 100644
--- a/rabies/analysis_pkg/utils.py
+++ b/rabies/analysis_pkg/utils.py
@@ -17,7 +17,7 @@ def resample_prior_maps(in_file, ref_file, transforms = [], inverses = []):
         warped_vol_fname = os.path.abspath(
             "deformed_volume" + str(x) + ".nii.gz")
         warped_volumes.append(warped_vol_fname)
-        exec_applyTransforms(transforms=transforms, inverses=inverses, input_image=volumes_list[x], ref_image=ref_file, output_image=warped_vol_fname, mask=False)
+        exec_applyTransforms(transforms=transforms, inverses=inverses, input_image=volumes_list[x], ref_image=ref_file, output_image=warped_vol_fname, interpolation='Linear')
 
     sample_volume = sitk.ReadImage(
         warped_volumes[0])
diff --git a/rabies/parser.py b/rabies/parser.py
index 0c9c35e..5c5d875 100644
--- a/rabies/parser.py
+++ b/rabies/parser.py
@@ -324,7 +324,7 @@ def get_parser():
         )
     g_registration.add_argument(
         '--anat_robust_inho_cor', type=str,
-        default='apply=false,masking=false,brain_extraction=false,template_registration=SyN',
+        default='apply=false,masking=false,brain_extraction=false,keep_mask_after_extract=false,template_registration=SyN',
         help=
             "When selecting this option, inhomogeneity correction is executed twice to optimize \n"
             "outcomes. After completing an initial inhomogeneity correction step, the corrected outputs \n"
@@ -342,6 +342,10 @@ def get_parser():
             " combined masks from inhomogeneity correction. This will enhance brain edge-matching, but \n"
             " requires good quality masks. This must be selected along the 'masking' option.\n"
             " *** Specify 'true' or 'false'. \n"
+            "* keep_mask_after_extract: If using brain_extraction, use the mask to compute the registration metric \n"
+            " within the mask only. Choose to prevent stretching of the images beyond the limit of the brain mask \n"
+            " (e.g. if the moving and target images don't have the same brain coverage).\n"
+            "*** Specify 'true' or 'false'. \n"
             "* template_registration: Specify a registration script for the alignment of the \n"
             " dataset-generated unbiased template to a reference template for masking.\n"
             "*** Rigid: conducts only rigid registration.\n"
@@ -361,7 +365,7 @@ def get_parser():
         )
     g_registration.add_argument(
         '--bold_robust_inho_cor', type=str,
-        default='apply=false,masking=false,brain_extraction=false,template_registration=SyN',
+        default='apply=false,masking=false,brain_extraction=false,keep_mask_after_extract=false,template_registration=SyN',
         help=
             "Same as --anat_robust_inho_cor, but for the EPI images.\n"
             "(default: %(default)s)\n"
@@ -369,7 +373,7 @@ def get_parser():
         )
     g_registration.add_argument(
         '--commonspace_reg', type=str,
-        default='masking=false,brain_extraction=false,template_registration=SyN,fast_commonspace=false',
+        default='masking=false,brain_extraction=false,keep_mask_after_extract=false,template_registration=SyN,fast_commonspace=false',
         help=
             "Specify registration options for the commonspace registration.\n"
             "* masking: Combine masks derived from the inhomogeneity correction step to support \n"
@@ -380,6 +384,10 @@ def get_parser():
             " combined masks from inhomogeneity correction. This will enhance brain edge-matching, but \n"
             " requires good quality masks. This must be selected along the 'masking' option.\n"
             "*** Specify 'true' or 'false'. \n"
+            "* keep_mask_after_extract: If using brain_extraction, use the mask to compute the registration metric \n"
+            " within the mask only. Choose to prevent stretching of the images beyond the limit of the brain mask \n"
+            " (e.g. if the moving and target images don't have the same brain coverage).\n"
+            "*** Specify 'true' or 'false'. \n"
             "* template_registration: Specify a registration script for the alignment of the \n"
             " dataset-generated unbiased template to the commonspace atlas.\n"
             "*** Rigid: conducts only rigid registration.\n"
@@ -411,7 +419,7 @@ def get_parser():
         )
     g_registration.add_argument(
         "--bold2anat_coreg", type=str, 
-        default='masking=false,brain_extraction=false,registration=SyN',
+        default='masking=false,brain_extraction=false,keep_mask_after_extract=false,registration=SyN',
         help=
             "Specify the registration script for cross-modal alignment between the EPI and structural\n"
             "images. This operation is responsible for correcting EPI susceptibility distortions.\n"
@@ -422,6 +430,10 @@ def get_parser():
             " inhomogeneity correction. This will enhance brain edge-matching, but requires good quality \n"
             " masks. This must be selected along the 'masking' option.\n"
             "*** Specify 'true' or 'false'. \n"
+            "* keep_mask_after_extract: If using brain_extraction, use the mask to compute the registration metric \n"
+            " within the mask only. Choose to prevent stretching of the images beyond the limit of the brain mask \n"
+            " (e.g. if the moving and target images don't have the same brain coverage).\n"
+            "*** Specify 'true' or 'false'. \n"
             "* registration: Specify a registration script.\n"
             "*** Rigid: conducts only rigid registration.\n"
             "*** Affine: conducts Rigid then Affine registration.\n"
@@ -467,6 +479,14 @@ def get_parser():
             "(default: %(default)s)\n"
             "\n"
         )
+    g_resampling.add_argument(
+        "--interpolation", type=str, default='Linear',
+        help=
+            "Select the interpolator which will be used by antsApplyTransforms (e.g. 'Linear' or 'BSpline[5]') \n"
+            "to resample preprocessed timeseries. \n"
+            "(default: %(default)s)\n"
+            "\n"
+        )    
 
     g_stc = preprocess.add_argument_group(
         title='STC Options', 
@@ -1052,24 +1072,34 @@ def read_parser(parser, args):
             name='bold_inho_cor')
 
         opts.commonspace_reg = parse_argument(opt=opts.commonspace_reg, 
-            key_value_pairs = {'masking':['true', 'false'], 'brain_extraction':['true', 'false'], 
+            key_value_pairs = {'masking':['true', 'false'], 'brain_extraction':['true', 'false'], 'keep_mask_after_extract':['true', 'false'], 
                 'template_registration':['Rigid', 'Affine', 'SyN', 'no_reg'], 'fast_commonspace':['true', 'false']},
             name='commonspace_reg')
 
         opts.bold2anat_coreg = parse_argument(opt=opts.bold2anat_coreg, 
-            key_value_pairs = {'masking':['true', 'false'], 'brain_extraction':['true', 'false'], 
+            key_value_pairs = {'masking':['true', 'false'], 'brain_extraction':['true', 'false'], 'keep_mask_after_extract':['true', 'false'], 
                 'registration':['Rigid', 'Affine', 'SyN', 'no_reg']},
             name='bold2anat_coreg')
 
         opts.anat_robust_inho_cor = parse_argument(opt=opts.anat_robust_inho_cor, 
-            key_value_pairs = {'apply':['true', 'false'], 'masking':['true', 'false'], 'brain_extraction':['true', 'false'], 
+            key_value_pairs = {'apply':['true', 'false'], 'masking':['true', 'false'], 'brain_extraction':['true', 'false'], 'keep_mask_after_extract':['true', 'false'], 
                 'template_registration':['Rigid', 'Affine', 'SyN', 'no_reg']},
             name='anat_robust_inho_cor')
 
         opts.bold_robust_inho_cor = parse_argument(opt=opts.bold_robust_inho_cor, 
-            key_value_pairs = {'apply':['true', 'false'], 'masking':['true', 'false'], 'brain_extraction':['true', 'false'], 
+            key_value_pairs = {'apply':['true', 'false'], 'masking':['true', 'false'], 'brain_extraction':['true', 'false'], 'keep_mask_after_extract':['true', 'false'], 
                 'template_registration':['Rigid', 'Affine', 'SyN', 'no_reg']},
             name='bold_robust_inho_cor')
+        
+        # check that masking/extraction options are well set
+        for name,opt in zip(['--commonspace_reg','--bold2anat_coreg','--anat_robust_inho_cor','--bold_robust_inho_cor'],
+                            [opts.commonspace_reg, opts.bold2anat_coreg, opts.anat_robust_inho_cor, opts.bold_robust_inho_cor]):
+            if opt['brain_extraction']:
+                if not opt['masking']:
+                    raise ValueError(f"For {name}, 'masking' must be set to 'true' if 'brain_extraction' is set to 'true'.")
+            if opt['keep_mask_after_extract']:
+                if not opt['brain_extraction']:
+                    raise ValueError(f"For {name}, 'brain_extraction' must be set to 'true' if 'keep_mask_after_extract' is set to 'true'.")
 
     elif opts.rabies_stage == 'confound_correction':
         opts.frame_censoring = parse_argument(opt=opts.frame_censoring, 
diff --git a/rabies/preprocess_pkg/commonspace_reg.py b/rabies/preprocess_pkg/commonspace_reg.py
index 6857fa4..21c6d34 100644
--- a/rabies/preprocess_pkg/commonspace_reg.py
+++ b/rabies/preprocess_pkg/commonspace_reg.py
@@ -11,7 +11,7 @@
 from .preprocess_visual_QC import PlotOverlap,template_masking
 
 
-def init_commonspace_reg_wf(opts, commonspace_masking, brain_extraction, template_reg, fast_commonspace, inherit_unbiased, output_folder, transforms_datasink, num_procs, output_datasinks, joinsource_list, name='commonspace_reg_wf'):
+def init_commonspace_reg_wf(opts, commonspace_masking, brain_extraction, keep_mask_after_extract, template_reg, fast_commonspace, inherit_unbiased, output_folder, transforms_datasink, num_procs, output_datasinks, joinsource_list, name='commonspace_reg_wf'):
     # commonspace_wf_head_start
     """
     This workflow handles the alignment of all MRI sessions to a common space. This is conducted first by generating
@@ -42,6 +42,10 @@ def init_commonspace_reg_wf(opts, commonspace_masking, brain_extraction, templat
                                 combined masks from inhomogeneity correction. This will enhance brain edge-matching, but 
                                 requires good quality masks. This should be selected along the 'masking' option.
                                 *** Specify 'true' or 'false'. 
+                                * keep_mask_after_extract: If using brain_extraction, use the mask to compute the registration metric
+                                within the mask only. Choose to prevent stretching of the images beyond the limit of the brain mask
+                                (e.g. if the moving and target images don't have the same brain coverage).
+                                *** Specify 'true' or 'false'.
                                 * template_registration: Specify a registration script for the alignment of the 
                                 dataset-generated unbiased template to the commonspace atlas.
                                 *** Rigid: conducts only rigid registration.
@@ -53,13 +57,14 @@ def init_commonspace_reg_wf(opts, commonspace_masking, brain_extraction, templat
                                 template_registration. This option can be faster, but may decrease the quality of 
                                 alignment between subjects. 
                                 *** Specify 'true' or 'false'. 
-                                (default: masking=false,brain_extraction=false,template_registration=SyN,fast_commonspace=false)
+                                (default: masking=false,brain_extraction=false,keep_mask_after_extract=false,template_registration=SyN,fast_commonspace=false)
                         
     Workflow:
         parameters
             opts: command line interface parameters
             commonspace_masking: whether masking is applied during template generation and registration
             brain_extraction: whether brain extraction is applied for template registration
+            keep_mask_after_extract: whether to keep using mask to delineate metric computation after brain extraction.
             template_reg: registration method
             fast_commonspace: whether the template generation step is skipped and instead each scan is registered directly in commonspace
             output_folder: specify a folder to execute the workflow and store important outputs
@@ -124,7 +129,7 @@ def init_commonspace_reg_wf(opts, commonspace_masking, brain_extraction, templat
                                                     'inverse_warp', 'warped_image']),
                                             name="atlas_reg_inherited")
     else:
-        atlas_reg = pe.Node(Function(input_names=['reg_method', 'brain_extraction', 'moving_image', 'moving_mask', 'fixed_image', 'fixed_mask', 'rabies_data_type'],
+        atlas_reg = pe.Node(Function(input_names=['reg_method', 'brain_extraction', 'keep_mask_after_extract', 'moving_image', 'moving_mask', 'fixed_image', 'fixed_mask', 'rabies_data_type'],
                                         output_names=['affine', 'warp',
                                                     'inverse_warp', 'warped_image'],
                                         function=run_antsRegistration),
@@ -136,6 +141,7 @@ def init_commonspace_reg_wf(opts, commonspace_masking, brain_extraction, templat
                 brain_extraction=False
 
         atlas_reg.inputs.brain_extraction = brain_extraction
+        atlas_reg.inputs.keep_mask_after_extract = keep_mask_after_extract
         atlas_reg.inputs.rabies_data_type = opts.data_type
         atlas_reg.plugin_args = {
             'qsub_args': f'-pe smp {str(3*opts.min_proc)}', 'overwrite': True}
@@ -261,7 +267,7 @@ def resample_unbiased_mask(unbiased_template, template_mask,to_atlas_affine,to_a
             transform_list=[to_atlas_affine,to_atlas_inverse_warp]
             inverse_list=[1,0]
             exec_applyTransforms(transforms = transform_list, inverses = inverse_list, 
-                input_image = template_mask, ref_image = unbiased_template, output_image = unbiased_mask, mask=True)
+                input_image = template_mask, ref_image = unbiased_template, output_image = unbiased_mask, interpolation='GenericLabel')
             return unbiased_mask
 
         resample_unbiased_mask_node = pe.Node(Function(input_names=['unbiased_template', 'template_mask','to_atlas_affine','to_atlas_inverse_warp'],
@@ -361,13 +367,14 @@ def resample_unbiased_mask(unbiased_template, template_mask,to_atlas_affine,to_a
                                                                             'unbiased_to_atlas_inverse_warp', 'warped_unbiased']),
                                                 name="inherit_unbiased_inputnode")
 
-            inherit_unbiased_reg_node = pe.Node(Function(input_names=['reg_method', 'brain_extraction', 'moving_image', 'moving_mask', 'fixed_image', 'fixed_mask', 'rabies_data_type'],
+            inherit_unbiased_reg_node = pe.Node(Function(input_names=['reg_method', 'brain_extraction', 'keep_mask_after_extract', 'moving_image', 'moving_mask', 'fixed_image', 'fixed_mask', 'rabies_data_type'],
                                             output_names=['affine', 'warp',
                                                         'inverse_warp', 'warped_image'],
                                             function=run_antsRegistration),
                                 name='inherit_unbiased_reg', mem_gb=2*opts.scale_min_memory)
             inherit_unbiased_reg_node.inputs.reg_method = 'Rigid' # this is the registration modelbuild conducts
             inherit_unbiased_reg_node.inputs.brain_extraction = False # brain extraction is not applied during template building
+            inherit_unbiased_reg_node.inputs.keep_mask_after_extract = False # brain extraction is not applied during template building
             inherit_unbiased_reg_node.inputs.rabies_data_type = opts.data_type
 
             if commonspace_masking: # this parameter should be inherited from previous run
@@ -544,7 +551,7 @@ def prep_commonspace_transform(native_ref, atlas_mask, native_to_unbiased_affine
     from rabies.utils import exec_applyTransforms
     # resample the atlas brain mask to native space
     exec_applyTransforms(transforms = commonspace_to_native_transform_list, inverses = commonspace_to_native_inverse_list, 
-        input_image = atlas_mask, ref_image = native_ref, output_image = native_mask, mask=True)
+        input_image = atlas_mask, ref_image = native_ref, output_image = native_mask, interpolation='GenericLabel')
 
     return native_mask, native_to_commonspace_transform_list,native_to_commonspace_inverse_list,commonspace_to_native_transform_list,commonspace_to_native_inverse_list
 
diff --git a/rabies/preprocess_pkg/hmc.py b/rabies/preprocess_pkg/hmc.py
index 18960a2..ce81ca3 100644
--- a/rabies/preprocess_pkg/hmc.py
+++ b/rabies/preprocess_pkg/hmc.py
@@ -239,7 +239,7 @@ def register_slice(fixed_image, moving_image):
     registration_method.SetMetricSamplingStrategy(registration_method.NONE)
     # registration_method.SetMetricSamplingPercentage(0.01)
 
-    registration_method.SetInterpolator(sitk.sitkBSpline)
+    registration_method.SetInterpolator(sitk.sitkLinear)
 
     # Optimizer settings.
     registration_method.SetOptimizerAsGradientDescent(
@@ -278,7 +278,7 @@ def slice_specific_registration(i, ref_file, timeseries_file):
 
         final_transform = register_slice(fixed_image, moving_image)
         moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform,
-                                         sitk.sitkBSplineResamplerOrder4, 0.0, moving_image.GetPixelID())
+                                         sitk.sitkLinear, 0.0, moving_image.GetPixelID())
 
         resampled_slice = sitk.GetArrayFromImage(moving_resampled)
         volume_array[:, j, :] = resampled_slice
diff --git a/rabies/preprocess_pkg/inho_correction.py b/rabies/preprocess_pkg/inho_correction.py
index ed4e10b..5627dbd 100644
--- a/rabies/preprocess_pkg/inho_correction.py
+++ b/rabies/preprocess_pkg/inho_correction.py
@@ -59,13 +59,17 @@ def init_inho_correction_wf(opts, image_type, output_folder, num_procs, name='in
                                 combined masks from inhomogeneity correction. This will enhance brain edge-matching, but 
                                 requires good quality masks. This should be selected along the 'masking' option.
                                 *** Specify 'true' or 'false'. 
+                                * keep_mask_after_extract: If using brain_extraction, use the mask to compute the registration metric
+                                within the mask only. Choose to prevent stretching of the images beyond the limit of the brain mask
+                                (e.g. if the moving and target images don't have the same brain coverage).
+                                *** Specify 'true' or 'false'.
                                 * template_registration: Specify a registration script for the alignment of the 
                                 dataset-generated unbiased template to a reference template for masking.
                                 *** Rigid: conducts only rigid registration.
                                 *** Affine: conducts Rigid then Affine registration.
                                 *** SyN: conducts Rigid, Affine then non-linear registration.
                                 *** no_reg: skip registration.
-                                (default: apply=false,masking=false,brain_extraction=false,template_registration=SyN)
+                                (default: apply=false,masking=false,brain_extraction=false,keep_mask_after_extract=false,template_registration=SyN)
                                 
         --bold_inho_cor BOLD_INHO_COR
                                 Same as --anat_inho_cor, but for the EPI images.
@@ -73,7 +77,7 @@ def init_inho_correction_wf(opts, image_type, output_folder, num_procs, name='in
                                 
         --bold_robust_inho_cor BOLD_ROBUST_INHO_COR
                                 Same as --anat_robust_inho_cor, but for the EPI images.
-                                (default: apply=false,masking=false,brain_extraction=false,template_registration=SyN)
+                                (default: apply=false,masking=false,brain_extraction=false,keep_mask_after_extract=false,template_registration=SyN)
                                 
     Workflow:
         parameters
@@ -127,6 +131,7 @@ def init_inho_correction_wf(opts, image_type, output_folder, num_procs, name='in
             multistage_otsu='true'
         if opts.bold_robust_inho_cor['apply']:
             robust_inho_cor=True
+            robust_inho_cor_opts=opts.bold_robust_inho_cor
             commonspace_wf_name='bold_robust_inho_cor_template'
             if not opts.bold_only:
                 joinsource_list = ['run_split', 'main_split']
@@ -140,6 +145,7 @@ def init_inho_correction_wf(opts, image_type, output_folder, num_procs, name='in
             multistage_otsu='true'
         if opts.anat_robust_inho_cor['apply']:
             robust_inho_cor=True
+            robust_inho_cor_opts=opts.anat_robust_inho_cor
             commonspace_wf_name='anat_robust_inho_cor_template'
             joinsource_list = ['main_split']
     else:
@@ -153,8 +159,8 @@ def init_inho_correction_wf(opts, image_type, output_folder, num_procs, name='in
                             name='init_InhoCorrection', mem_gb=0.6*opts.scale_min_memory)
 
         from .commonspace_reg import init_commonspace_reg_wf
-        commonspace_reg_wf = init_commonspace_reg_wf(opts=opts, commonspace_masking=opts.bold_robust_inho_cor['masking'], brain_extraction=opts.bold_robust_inho_cor['brain_extraction'], 
-                                                     template_reg=opts.bold_robust_inho_cor['template_registration'], fast_commonspace=False, inherit_unbiased=False, output_folder=output_folder, 
+        commonspace_reg_wf = init_commonspace_reg_wf(opts=opts, commonspace_masking=robust_inho_cor_opts['masking'], brain_extraction=robust_inho_cor_opts['brain_extraction'], keep_mask_after_extract=robust_inho_cor_opts['keep_mask_after_extract'], 
+                                                     template_reg=robust_inho_cor_opts['template_registration'], fast_commonspace=False, inherit_unbiased=False, output_folder=output_folder, 
                                                      transforms_datasink=None, num_procs=num_procs, output_datasinks=False, joinsource_list=joinsource_list, name=commonspace_wf_name)
 
         workflow.connect([
diff --git a/rabies/preprocess_pkg/main_wf.py b/rabies/preprocess_pkg/main_wf.py
index c073818..cafcbdb 100644
--- a/rabies/preprocess_pkg/main_wf.py
+++ b/rabies/preprocess_pkg/main_wf.py
@@ -245,7 +245,7 @@ def init_main_wf(data_dir_path, output_folder, opts, name='main_wf'):
     EPI_target_buffer = pe.Node(niu.IdentityInterface(fields=['EPI_template', 'EPI_mask']),
                                         name="EPI_target_buffer")
 
-    commonspace_reg_wf = init_commonspace_reg_wf(opts=opts, commonspace_masking=opts.commonspace_reg['masking'], brain_extraction=opts.commonspace_reg['brain_extraction'], 
+    commonspace_reg_wf = init_commonspace_reg_wf(opts=opts, commonspace_masking=opts.commonspace_reg['masking'], brain_extraction=opts.commonspace_reg['brain_extraction'], keep_mask_after_extract=opts.commonspace_reg['keep_mask_after_extract'], 
                                                  template_reg=opts.commonspace_reg['template_registration'], fast_commonspace=opts.commonspace_reg['fast_commonspace'], inherit_unbiased=inherit_unbiased,
                                                  output_folder=output_folder, transforms_datasink=transforms_datasink, num_procs=num_procs, output_datasinks=True, 
                                                  joinsource_list=['main_split'], name='commonspace_reg_wf')
diff --git a/rabies/preprocess_pkg/registration.py b/rabies/preprocess_pkg/registration.py
index f70f29b..9d1b37c 100644
--- a/rabies/preprocess_pkg/registration.py
+++ b/rabies/preprocess_pkg/registration.py
@@ -25,12 +25,16 @@ def init_cross_modal_reg_wf(opts, name='cross_modal_reg_wf'):
                             inhomogeneity correction. This will enhance brain edge-matching, but requires good quality 
                             masks. This should be selected along the 'masking' option.
                             *** Specify 'true' or 'false'. 
+                            * keep_mask_after_extract: If using brain_extraction, use the mask to compute the registration metric
+                            within the mask only. Choose to prevent stretching of the images beyond the limit of the brain mask
+                            (e.g. if the moving and target images don't have the same brain coverage).
+                            *** Specify 'true' or 'false'.
                             * registration: Specify a registration script.
                             *** Rigid: conducts only rigid registration.
                             *** Affine: conducts Rigid then Affine registration.
                             *** SyN: conducts Rigid, Affine then non-linear registration.
                             *** no_reg: skip registration.
-                            (default: masking=false,brain_extraction=false,registration=SyN)
+                            (default: masking=false,brain_extraction=false,keep_mask_after_extract=false,registration=SyN)
 
     Workflow:
         parameters
@@ -63,7 +67,7 @@ def init_cross_modal_reg_wf(opts, name='cross_modal_reg_wf'):
         name='outputnode'
     )
 
-    run_reg = pe.Node(Function(input_names=["reg_method", "brain_extraction", "moving_image", "moving_mask", "fixed_image",
+    run_reg = pe.Node(Function(input_names=["reg_method", "brain_extraction", "keep_mask_after_extract", "moving_image", "moving_mask", "fixed_image",
                                             "fixed_mask", "rabies_data_type"],
                                output_names=['bold_to_anat_affine', 'bold_to_anat_warp',
                                              'bold_to_anat_inverse_warp', 'output_warped_bold'],
@@ -72,6 +76,7 @@ def init_cross_modal_reg_wf(opts, name='cross_modal_reg_wf'):
     # don't use brain extraction without a moving mask
     run_reg.inputs.reg_method = opts.bold2anat_coreg['registration']
     run_reg.inputs.brain_extraction = opts.bold2anat_coreg['brain_extraction']
+    run_reg.inputs.keep_mask_after_extract = opts.bold2anat_coreg['keep_mask_after_extract']
     run_reg.inputs.rabies_data_type = opts.data_type
     run_reg.plugin_args = {
         'qsub_args': f'-pe smp {str(3*opts.min_proc)}', 'overwrite': True}
@@ -99,7 +104,7 @@ def init_cross_modal_reg_wf(opts, name='cross_modal_reg_wf'):
     return workflow
 
 
-def run_antsRegistration(reg_method, brain_extraction=False, moving_image='NULL', moving_mask='NULL', fixed_image='NULL', fixed_mask='NULL', rabies_data_type=8):
+def run_antsRegistration(reg_method, brain_extraction=False, keep_mask_after_extract=False, moving_image='NULL', moving_mask='NULL', fixed_image='NULL', fixed_mask='NULL', rabies_data_type=8):
     import os
     import pathlib  # Better path manipulation
     filename_split = pathlib.Path(moving_image).name.rsplit(".nii")
@@ -109,7 +114,9 @@ def run_antsRegistration(reg_method, brain_extraction=False, moving_image='NULL'
 
     if reg_method == 'Rigid' or reg_method == 'Affine' or reg_method == 'SyN':
         if brain_extraction:
-            reg_call+=" --mask-extract --keep-mask-after-extract"
+            reg_call+=" --mask-extract"
+            if keep_mask_after_extract:
+                reg_call+=" --keep-mask-after-extract"
         command = f"{reg_call} --moving-mask {moving_mask} --fixed-mask {fixed_mask} --resampled-output {filename_split[0]}_output_warped_image.nii.gz {moving_image} {fixed_image} {filename_split[0]}_output_"
     else:
         command = f'{reg_call} {moving_image} {moving_mask} {fixed_image} {fixed_mask} {filename_split[0]}'
diff --git a/rabies/preprocess_pkg/resampling.py b/rabies/preprocess_pkg/resampling.py
index 6064589..4d3c2f4 100644
--- a/rabies/preprocess_pkg/resampling.py
+++ b/rabies/preprocess_pkg/resampling.py
@@ -84,6 +84,7 @@ def init_bold_preproc_trans_wf(opts, resampling_dim, name='bold_native_trans_wf'
         rabies_data_type=opts.data_type), name='bold_transform', mem_gb=1*opts.scale_min_memory)
     bold_transform.inputs.apply_motcorr = (not opts.apply_slice_mc)
     bold_transform.inputs.resampling_dim = resampling_dim
+    bold_transform.inputs.interpolation = opts.interpolation
 
     merge = pe.Node(Merge(rabies_data_type=opts.data_type, clip_negative=True), name='merge', mem_gb=4*opts.scale_min_memory)
     merge.plugin_args = {
@@ -224,7 +225,7 @@ def _run_interface(self, runtime):
         else:
             new_mask_path = os.path.abspath(f'{filename_split[0]}_{self.inputs.name_spec}.nii.gz')
 
-        exec_applyTransforms(self.inputs.transforms, self.inputs.inverses, self.inputs.mask, self.inputs.ref_EPI, new_mask_path, mask=True)
+        exec_applyTransforms(self.inputs.transforms, self.inputs.inverses, self.inputs.mask, self.inputs.ref_EPI, new_mask_path, interpolation='GenericLabel')
         sitk.WriteImage(sitk.ReadImage(
             new_mask_path, sitk.sitkInt16), new_mask_path)
 
diff --git a/rabies/utils.py b/rabies/utils.py
index 4c1d8c6..8e96dcb 100644
--- a/rabies/utils.py
+++ b/rabies/utils.py
@@ -85,9 +85,7 @@ def resample_image_spacing_4d(image_4d, output_spacing, clip_negative=True):
     sitk.ProcessObject_SetGlobalDefaultThreader('Platform')
     resampled_list=[]
     for i in range(input_size[3]):
-        # the resampling is done with BSpline, the interpolator order is 3
-        # other BSpline options are blurry, see https://discourse.itk.org/t/resample-produces-blurry-results-when-just-cropping/4473
-        resampled_image = sitk.Resample(image_4d[:,:,:,i], output_size, identity, sitk.sitkBSpline,
+        resampled_image = sitk.Resample(image_4d[:,:,:,i], output_size, identity, sitk.sitkLinear,
                                         origin, output_spacing, direction_3d)
         resampled_list.append(resampled_image)
     combined = sitk.JoinSeries(resampled_list) 
@@ -157,6 +155,8 @@ class slice_applyTransformsInputSpec(BaseInterfaceInputSpec):
         exists=True, desc="xforms from head motion estimation .csv file")
     resampling_dim = traits.Str(
         desc="Specification for the dimension of resampling.")
+    interpolation = traits.Str(
+        desc="Select the interpolator for antsApplyTransform.")
     rabies_data_type = traits.Int(mandatory=True,
                                   desc="Integer specifying SimpleITK data type.")
 
@@ -216,7 +216,7 @@ def _run_interface(self, runtime):
                 transforms = orig_transforms
                 inverses = orig_inverses
 
-            exec_applyTransforms(transforms, inverses, bold_volumes[x], ref_img, warped_vol_fname, mask=False)
+            exec_applyTransforms(transforms, inverses, bold_volumes[x], ref_img, warped_vol_fname, interpolation=self.inputs.interpolation)
             # change image to specified data type
             sitk.WriteImage(sitk.ReadImage(warped_vol_fname,
                                            self.inputs.rabies_data_type), warped_vol_fname)
@@ -228,7 +228,7 @@ def _list_outputs(self):
         return {'out_files': getattr(self, 'out_files')}
 
 
-def exec_applyTransforms(transforms, inverses, input_image, ref_image, output_image, mask=False):
+def exec_applyTransforms(transforms, inverses, input_image, ref_image, output_image, interpolation):
     # tranforms is a list of transform files, set in order of call within antsApplyTransforms
     transform_string = ""
     for transform, inverse in zip(transforms, inverses):
@@ -239,11 +239,6 @@ def exec_applyTransforms(transforms, inverses, input_image, ref_image, output_im
         else:
             transform_string += f"-t {transform} "
 
-    if mask:
-        interpolation = 'GenericLabel'
-    else:
-        interpolation = 'BSpline[5]'
-
     command = f'antsApplyTransforms -i {input_image} {transform_string}-n {interpolation} -r {ref_image} -o {output_image}'
     rc,c_out = run_command(command)
     if not os.path.isfile(output_image):
diff --git a/scripts/error_check_rabies.py b/scripts/error_check_rabies.py
index 86cb3a2..a46f44d 100755
--- a/scripts/error_check_rabies.py
+++ b/scripts/error_check_rabies.py
@@ -14,6 +14,7 @@
     --anat_robust_inho_cor/bold_robust_inho_cor: requires MINC/ANTs, which doesn't run on token data
     --commonspace_reg: requires MINC/ANTs, which doesn't run on token data
     --bold2anat_coreg: requires MINC/ANTs, which doesn't run on token data
+    --interpolation: is not tested.
 
 confound_correction:
     --highpass/lowpass/edge_cutoff; since filtering doesn't work with 3 timepoints
@@ -58,8 +59,8 @@ def get_parser():
             "   --anat_template {tmppath}/inputs/sub-token1_T1w.nii.gz --brain_mask {tmppath}/inputs/token_mask.nii.gz  \ \n"
             "   --WM_mask {tmppath}/inputs/token_mask.nii.gz --CSF_mask {tmppath}/inputs/token_mask.nii.gz  \ \n"
             "   --vascular_mask {tmppath}/inputs/token_mask.nii.gz --labels {tmppath}/inputs/token_mask.nii.gz \ \n"
-            "   --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false \ \n"
-            "   --commonspace_reg masking=false,brain_extraction=false,fast_commonspace=true,template_registration=no_reg --data_type int16 \n"
+            "   --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false,keep_mask_after_extract=false \ \n"
+            "   --commonspace_reg masking=false,brain_extraction=false,keep_mask_after_extract=false,fast_commonspace=true,template_registration=no_reg --data_type int16 \n"
         )
     
     return parser
@@ -77,7 +78,7 @@ def get_parser():
 if not opts.custom is None:
     minimal_preproc = f"rabies --inclusion_ids {tmppath}/inputs/sub-token1_bold.nii.gz --verbose 1 preprocess {tmppath}/inputs {tmppath}/outputs --anat_inho_cor method=disable,otsu_thresh=2,multiotsu=false --bold_inho_cor method=disable,otsu_thresh=2,multiotsu=false \
         --anat_template {tmppath}/inputs/sub-token1_T1w.nii.gz --brain_mask {tmppath}/inputs/token_mask.nii.gz --WM_mask {tmppath}/inputs/token_mask.nii.gz --CSF_mask {tmppath}/inputs/token_mask.nii.gz --vascular_mask {tmppath}/inputs/token_mask.nii.gz --labels {tmppath}/inputs/token_mask.nii.gz \
-        --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false --commonspace_reg masking=false,brain_extraction=false,fast_commonspace=true,template_registration=no_reg --data_type int16"
+        --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false,keep_mask_after_extract=false --commonspace_reg masking=false,brain_extraction=false,keep_mask_after_extract=false,fast_commonspace=true,template_registration=no_reg --data_type int16"
     minimal_cc = f"rabies --verbose 1 confound_correction {tmppath}/outputs {tmppath}/outputs"
 
     import sys
@@ -85,7 +86,7 @@ def get_parser():
     if 'preprocess' in command:
         command += f" --anat_inho_cor method=disable,otsu_thresh=2,multiotsu=false --bold_inho_cor method=disable,otsu_thresh=2,multiotsu=false \
     --anat_template {tmppath}/inputs/sub-token1_T1w.nii.gz --brain_mask {tmppath}/inputs/token_mask.nii.gz --WM_mask {tmppath}/inputs/token_mask.nii.gz --CSF_mask {tmppath}/inputs/token_mask.nii.gz --vascular_mask {tmppath}/inputs/token_mask.nii.gz --labels {tmppath}/inputs/token_mask.nii.gz \
-    --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false --commonspace_reg masking=false,brain_extraction=false,fast_commonspace=true,template_registration=no_reg --data_type int16"
+    --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false,keep_mask_after_extract=false --commonspace_reg masking=false,brain_extraction=false,keep_mask_after_extract=false,fast_commonspace=true,template_registration=no_reg --data_type int16"
         command += f" {tmppath}/inputs {tmppath}/outputs"
         
     if 'confound_correction' in command or 'analysis' in command:
@@ -115,7 +116,7 @@ def get_parser():
 
 command = f"rabies --exclusion_ids {tmppath}/inputs/sub-token2_bold.nii.gz {tmppath}/inputs/sub-token3_bold.nii.gz --force --verbose 1 preprocess {tmppath}/inputs {tmppath}/outputs --anat_inho_cor method=disable,otsu_thresh=2,multiotsu=false --bold_inho_cor method=disable,otsu_thresh=2,multiotsu=false \
     --anat_template {tmppath}/inputs/sub-token1_T1w.nii.gz --brain_mask {tmppath}/inputs/token_mask.nii.gz --WM_mask {tmppath}/inputs/token_mask.nii.gz --CSF_mask {tmppath}/inputs/token_mask.nii.gz --vascular_mask {tmppath}/inputs/token_mask.nii.gz --labels {tmppath}/inputs/token_mask.nii.gz \
-    --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false --commonspace_reg masking=false,brain_extraction=false,fast_commonspace=true,template_registration=no_reg --data_type int16 --bold_only --detect_dummy \
+    --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false,keep_mask_after_extract=false --commonspace_reg masking=false,brain_extraction=false,keep_mask_after_extract=false,fast_commonspace=true,template_registration=no_reg --data_type int16 --bold_only --detect_dummy \
     --tpattern seq-z --apply_STC --voxelwise_motion --isotropic_HMC --interp_method linear --nativespace_resampling 1x1x1 --commonspace_resampling 1x1x1 --anatomical_resampling 1x1x1 --oblique2card 3dWarp"
 process = subprocess.run(
     command,
@@ -125,7 +126,7 @@ def get_parser():
 
 command = f"rabies --inclusion_ids {tmppath}/inputs/sub-token1_bold.nii.gz --verbose 1 --force preprocess {tmppath}/inputs {tmppath}/outputs --anat_inho_cor method=disable,otsu_thresh=2,multiotsu=false --bold_inho_cor method=disable,otsu_thresh=2,multiotsu=false \
     --anat_template {tmppath}/inputs/sub-token1_T1w.nii.gz --brain_mask {tmppath}/inputs/token_mask.nii.gz --WM_mask {tmppath}/inputs/token_mask.nii.gz --CSF_mask {tmppath}/inputs/token_mask.nii.gz --vascular_mask {tmppath}/inputs/token_mask.nii.gz --labels {tmppath}/inputs/token_mask.nii.gz \
-    --bold2anat_coreg registration=no_reg,masking=true,brain_extraction=true --commonspace_reg masking=true,brain_extraction=true,fast_commonspace=true,template_registration=no_reg --data_type int16  \
+    --bold2anat_coreg registration=no_reg,masking=true,brain_extraction=true,keep_mask_after_extract=false --commonspace_reg masking=true,brain_extraction=true,keep_mask_after_extract=false,fast_commonspace=true,template_registration=no_reg --data_type int16  \
     --HMC_option 0 --apply_despiking --anat_autobox --bold_autobox --oblique2card affine"
 process = subprocess.run(
     command,
@@ -197,7 +198,7 @@ def get_parser():
     ####GROUP LEVEL, RUNNING ALL 3 SCANS####
     command = f"rabies --force --verbose 1 preprocess {tmppath}/inputs {tmppath}/outputs --anat_inho_cor method=disable,otsu_thresh=2,multiotsu=false --bold_inho_cor method=disable,otsu_thresh=2,multiotsu=false \
         --anat_template {tmppath}/inputs/sub-token1_T1w.nii.gz --brain_mask {tmppath}/inputs/token_mask.nii.gz --WM_mask {tmppath}/inputs/token_mask_half.nii.gz --CSF_mask {tmppath}/inputs/token_mask_half.nii.gz --vascular_mask {tmppath}/inputs/token_mask_half.nii.gz --labels {tmppath}/inputs/token_mask.nii.gz \
-        --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false --commonspace_reg masking=false,brain_extraction=false,fast_commonspace=true,template_registration=no_reg --data_type int16  \
+        --bold2anat_coreg registration=no_reg,masking=false,brain_extraction=false,keep_mask_after_extract=false --commonspace_reg masking=false,brain_extraction=false,keep_mask_after_extract=false,fast_commonspace=true,template_registration=no_reg --data_type int16  \
         --HMC_option 0"
     process = subprocess.run(
         command,