From 89e5f82031d432751566a8fcc3e3d594d7dd240b Mon Sep 17 00:00:00 2001 From: Gab-D-G Date: Wed, 8 Nov 2023 10:35:45 -0500 Subject: [PATCH 1/4] The resampling interpolator was changed from BSpline to Linear for all resampling operations involving 4D images. Only the resample_image_spacing function is still on BSpline. --- rabies/preprocess_pkg/hmc.py | 4 ++-- rabies/utils.py | 6 ++---- 2 files changed, 4 insertions(+), 6 deletions(-) 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/utils.py b/rabies/utils.py index 4c1d8c6..c4c59df 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) @@ -242,7 +240,7 @@ def exec_applyTransforms(transforms, inverses, input_image, ref_image, output_im if mask: interpolation = 'GenericLabel' else: - interpolation = 'BSpline[5]' + interpolation = 'Linear' command = f'antsApplyTransforms -i {input_image} {transform_string}-n {interpolation} -r {ref_image} -o {output_image}' rc,c_out = run_command(command) From 1593ca597b17745073271ef01cbaf1d69e14a114 Mon Sep 17 00:00:00 2001 From: Gab-D-G Date: Fri, 15 Dec 2023 14:40:35 -0500 Subject: [PATCH 2/4] Exposed the --keep-mask-after-extract registration parameter, which was always turned on with brain extraction. It is now enforced that masking should be set to true if using brain extraction, since brain extraction isn't doing anything otherwise. --- rabies/parser.py | 38 +++++++++++++++++++----- rabies/preprocess_pkg/commonspace_reg.py | 15 +++++++--- rabies/preprocess_pkg/inho_correction.py | 10 +++++-- rabies/preprocess_pkg/main_wf.py | 2 +- rabies/preprocess_pkg/registration.py | 15 +++++++--- scripts/error_check_rabies.py | 14 ++++----- 6 files changed, 67 insertions(+), 27 deletions(-) diff --git a/rabies/parser.py b/rabies/parser.py index 0c9c35e..85200fc 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" @@ -1052,24 +1064,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..1c676c8 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} @@ -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 diff --git a/rabies/preprocess_pkg/inho_correction.py b/rabies/preprocess_pkg/inho_correction.py index ed4e10b..907ea78 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 @@ -153,7 +157,7 @@ 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'], + 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'], keep_mask_after_extract=opts.bold_robust_inho_cor['keep_mask_after_extract'], template_reg=opts.bold_robust_inho_cor['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) 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/scripts/error_check_rabies.py b/scripts/error_check_rabies.py index 86cb3a2..d4a6c37 100755 --- a/scripts/error_check_rabies.py +++ b/scripts/error_check_rabies.py @@ -58,8 +58,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 +77,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 +85,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 +115,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 +125,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 +197,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, From 0e7381a7dc2cafcf1e3ba53b609347340aaedacd Mon Sep 17 00:00:00 2001 From: Gab-D-G Date: Fri, 15 Dec 2023 14:47:24 -0500 Subject: [PATCH 3/4] Bug: the robust_anat_inho_cor parameters were not fed properly. The ones from robust_bold_inho_cor were always inputed. --- rabies/preprocess_pkg/inho_correction.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/rabies/preprocess_pkg/inho_correction.py b/rabies/preprocess_pkg/inho_correction.py index 907ea78..5627dbd 100644 --- a/rabies/preprocess_pkg/inho_correction.py +++ b/rabies/preprocess_pkg/inho_correction.py @@ -131,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'] @@ -144,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: @@ -157,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'], keep_mask_after_extract=opts.bold_robust_inho_cor['keep_mask_after_extract'], - 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([ From 96b7be5860ed56559d5ac789c9e26030555b1702 Mon Sep 17 00:00:00 2001 From: Gab-D-G Date: Fri, 15 Dec 2023 15:09:06 -0500 Subject: [PATCH 4/4] Exposed the interpolator for antsApplyTransform for resampling of timeseries. --- rabies/analysis_pkg/utils.py | 2 +- rabies/parser.py | 8 ++++++++ rabies/preprocess_pkg/commonspace_reg.py | 4 ++-- rabies/preprocess_pkg/resampling.py | 3 ++- rabies/utils.py | 11 ++++------- scripts/error_check_rabies.py | 1 + 6 files changed, 18 insertions(+), 11 deletions(-) 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 85200fc..5c5d875 100644 --- a/rabies/parser.py +++ b/rabies/parser.py @@ -479,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', diff --git a/rabies/preprocess_pkg/commonspace_reg.py b/rabies/preprocess_pkg/commonspace_reg.py index 1c676c8..21c6d34 100644 --- a/rabies/preprocess_pkg/commonspace_reg.py +++ b/rabies/preprocess_pkg/commonspace_reg.py @@ -267,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'], @@ -551,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/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 c4c59df..8e96dcb 100644 --- a/rabies/utils.py +++ b/rabies/utils.py @@ -155,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.") @@ -214,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) @@ -226,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): @@ -237,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 = 'Linear' - 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 d4a6c37..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