diff --git a/macapype/nodes/prepare.py b/macapype/nodes/prepare.py index 6f2a4cfa8..19feeb597 100644 --- a/macapype/nodes/prepare.py +++ b/macapype/nodes/prepare.py @@ -294,8 +294,17 @@ def average_align(list_img, reorient=False): img_0 = nib.load(list_img[0]) path, fname, ext = split_f(list_img[0]) - list_data = [img_0.get_data()] + data_0 = img_0.get_fdata() + + data_0_shape = data_0.shape + list_data = [data_0] + for i, img in enumerate(list_img[1:]): + data_orig_shape = nib.load(img).get_fdata().shape + if data_orig_shape != data_0_shape: + print(f"**** Warning , original image shape \ + {data_orig_shape} \ + have a different shape than ref {data_0_shape}") print("running flirt on {}".format(img)) flirt = fsl.FLIRT(dof=6) @@ -306,11 +315,21 @@ def average_align(list_img, reorient=False): out_file = flirt.run().outputs.out_file print(out_file) - data = nib.load(out_file).get_data() + data = nib.load(out_file).get_fdata() + assert data.shape == data_0_shape, \ + f"**** Error, flirt results {out_file} \ + have a different shape {data.shape} \ + than ref {data_0_shape}" + list_data.append(data) avg_data = np.mean(np.array(list_data), axis=0) - print(avg_data.shape) + + assert data_0_shape == avg_data.shape, \ + f"Error with average shape {avg_data.shape} != {data_0_shape}" + + # after casting type + avg_data = avg_data.astype(np.int16) av_img_file = os.path.abspath("avg_" + fname + ext) nib.save(nib.Nifti1Image(avg_data, header=img_0.header, diff --git a/macapype/nodes/register.py b/macapype/nodes/register.py index eba3cdf08..3b1fa9d94 100644 --- a/macapype/nodes/register.py +++ b/macapype/nodes/register.py @@ -190,15 +190,14 @@ class IterREGBETInputSpec(CommandLineInputSpec): mandatory=False, argstr="-refw %s") k = traits.Bool( - False, usedefault=True, - position=3, argstr="-k", + False, usedefault=True, argstr="-k", desc="Will keep temporary files", mandatory=False) p = traits.String( desc="Prefix for running FSL functions\ (can be a path or just a prefix)", - position=3, argstr="-p %s") + argstr="-p %s") class IterREGBETOutputSpec(TraitedSpec): @@ -749,11 +748,26 @@ def pad_zero_mri(img_file, pad_val=10, const=0): img_arr = np.array(img.dataobj) img_arr_copy = np.copy(img_arr) - img_arr_copy_padded = np.pad( - img_arr_copy, - pad_width=pad_val, - mode='constant', - constant_values=const) + if isinstance(pad_val, int): + print("running with pad_val") + + img_arr_copy_padded = np.pad( + img_arr_copy, + pad_width=pad_val, + mode='constant', + constant_values=const) + + elif pad_val == 'double': + + pad_val = np.max(img_arr.shape) + + print(f"running with double pad_val = {pad_val}") + + img_arr_copy_padded = np.pad( + img_arr_copy, + pad_width=pad_val, + mode='constant', + constant_values=const) img_padded = nib.Nifti1Image(img_arr_copy_padded, header=img.header, diff --git a/macapype/pipelines/full_pipelines.py b/macapype/pipelines/full_pipelines.py index 482a6dc35..95140c091 100644 --- a/macapype/pipelines/full_pipelines.py +++ b/macapype/pipelines/full_pipelines.py @@ -514,111 +514,118 @@ def create_brain_segment_from_mask_pipe( "prob_csf", "gen_5tt"]), name='outputnode') - if "register_NMT_pipe" in params: + if "use_priors" in params["segment_atropos_pipe"].keys(): - # register NMT template, template mask and priors to subject T1 - register_NMT_pipe = create_register_NMT_pipe( - params_template=params_template, - params=parse_key(params, "register_NMT_pipe")) - - brain_segment_pipe.connect( - inputnode, 'masked_debiased_T1', - register_NMT_pipe, "inputnode.T1") - - brain_segment_pipe.connect( - inputnode, 'indiv_params', - register_NMT_pipe, "inputnode.indiv_params") + if "register_NMT_pipe" in params: - elif "reg" in params: - # Iterative registration to the INIA19 template - reg = NodeParams( - IterREGBET(), - params=parse_key(params, "reg"), - name='reg') + # register NMT template, template mask and priors to subject T1 + register_NMT_pipe = create_register_NMT_pipe( + params_template=params_template, + params=parse_key(params, "register_NMT_pipe")) - reg.inputs.refb_file = params_template["template_brain"] + brain_segment_pipe.connect( + inputnode, 'masked_debiased_T1', + register_NMT_pipe, "inputnode.T1") - brain_segment_pipe.connect( - inputnode, 'debiased_T1', - reg, 'inw_file') + brain_segment_pipe.connect( + inputnode, 'indiv_params', + register_NMT_pipe, "inputnode.indiv_params") - brain_segment_pipe.connect( - inputnode, 'masked_debiased_T1', - reg, 'inb_file') + elif "reg" in params: + # Iterative registration to the template + reg = NodeParams( + IterREGBET(), + params=parse_key(params, "reg"), + name='reg') - brain_segment_pipe.connect( - inputnode, ('indiv_params', parse_key, "reg"), - reg, 'indiv_params') + reg.inputs.refb_file = params_template["template_brain"] - if "template_seg" in params_template.keys(): + reg.inputs.refw_file = params_template["template_head"] + reg.inputs.k = True - # seg - register_seg_to_nat = pe.Node( - fsl.ApplyXFM(), name="register_seg_to_nat") - register_seg_to_nat.inputs.interp = "nearestneighbour" + brain_segment_pipe.connect( + inputnode, 'debiased_T1', + reg, 'inw_file') - register_seg_to_nat.inputs.in_file = params_template[ - "template_seg"] brain_segment_pipe.connect( inputnode, 'masked_debiased_T1', - register_seg_to_nat, 'reference') + reg, 'inb_file') brain_segment_pipe.connect( - reg, 'inv_transfo_file', - register_seg_to_nat, "in_matrix_file") + inputnode, ('indiv_params', parse_key, "reg"), + reg, 'indiv_params') - else: - # gm - register_gm_to_nat = pe.Node( - fsl.ApplyXFM(), name="register_gm_to_nat") - register_gm_to_nat.inputs.output_type = "NIFTI" # for SPM segment - register_gm_to_nat.inputs.interp = "nearestneighbour" + if "template_seg" in params_template.keys(): - register_gm_to_nat.inputs.in_file = params_template["template_gm"] + # seg + register_seg_to_nat = pe.Node( + fsl.ApplyXFM(), name="register_seg_to_nat") + register_seg_to_nat.inputs.interp = "nearestneighbour" - brain_segment_pipe.connect( - inputnode, 'masked_debiased_T1', - register_gm_to_nat, 'reference') + register_seg_to_nat.inputs.in_file = params_template[ + "template_seg"] + brain_segment_pipe.connect( + inputnode, 'masked_debiased_T1', + register_seg_to_nat, 'reference') - brain_segment_pipe.connect( - reg, 'inv_transfo_file', - register_gm_to_nat, "in_matrix_file") + brain_segment_pipe.connect( + reg, 'inv_transfo_file', + register_seg_to_nat, "in_matrix_file") - # wm - register_wm_to_nat = pe.Node( - fsl.ApplyXFM(), name="register_wm_to_nat") - register_wm_to_nat.inputs.output_type = "NIFTI" # for SPM segment - register_wm_to_nat.inputs.interp = "nearestneighbour" + else: + # gm + register_gm_to_nat = pe.Node( + fsl.ApplyXFM(), name="register_gm_to_nat") + register_gm_to_nat.inputs.output_type = "NIFTI" # for SPM + register_gm_to_nat.inputs.interp = "nearestneighbour" - register_wm_to_nat.inputs.in_file = params_template["template_wm"] + register_gm_to_nat.inputs.in_file = \ + params_template["template_gm"] - brain_segment_pipe.connect( - inputnode, 'masked_debiased_T1', - register_wm_to_nat, 'reference') + brain_segment_pipe.connect( + inputnode, 'masked_debiased_T1', + register_gm_to_nat, 'reference') - brain_segment_pipe.connect( - reg, 'inv_transfo_file', - register_wm_to_nat, "in_matrix_file") + brain_segment_pipe.connect( + reg, 'inv_transfo_file', + register_gm_to_nat, "in_matrix_file") - # csf - register_csf_to_nat = pe.Node( - fsl.ApplyXFM(), name="register_csf_to_nat") - register_csf_to_nat.inputs.output_type = "NIFTI" # for SPM segment - register_csf_to_nat.inputs.interp = "nearestneighbour" + # wm + register_wm_to_nat = pe.Node( + fsl.ApplyXFM(), name="register_wm_to_nat") + register_wm_to_nat.inputs.output_type = "NIFTI" # for SPM + register_wm_to_nat.inputs.interp = "nearestneighbour" - register_csf_to_nat.inputs.in_file = params_template[ - "template_csf"] + register_wm_to_nat.inputs.in_file = \ + params_template["template_wm"] - brain_segment_pipe.connect( - inputnode, 'masked_debiased_T1', - register_csf_to_nat, 'reference') + brain_segment_pipe.connect( + inputnode, 'masked_debiased_T1', + register_wm_to_nat, 'reference') - brain_segment_pipe.connect( - reg, 'inv_transfo_file', - register_csf_to_nat, "in_matrix_file") - else: - print("##### Error, no coregistration method is defined") - return brain_segment_pipe + brain_segment_pipe.connect( + reg, 'inv_transfo_file', + register_wm_to_nat, "in_matrix_file") + + # csf + register_csf_to_nat = pe.Node( + fsl.ApplyXFM(), name="register_csf_to_nat") + register_csf_to_nat.inputs.output_type = "NIFTI" # for SPM + register_csf_to_nat.inputs.interp = "nearestneighbour" + + register_csf_to_nat.inputs.in_file = \ + params_template["template_csf"] + + brain_segment_pipe.connect( + inputnode, 'masked_debiased_T1', + register_csf_to_nat, 'reference') + + brain_segment_pipe.connect( + reg, 'inv_transfo_file', + register_csf_to_nat, "in_matrix_file") + else: + print("##### Error, no coregistration method is defined") + return brain_segment_pipe # ants Atropos if "template_seg" in params_template.keys(): diff --git a/macapype/pipelines/prepare.py b/macapype/pipelines/prepare.py index 67cf369b5..ee97debf8 100644 --- a/macapype/pipelines/prepare.py +++ b/macapype/pipelines/prepare.py @@ -4,9 +4,9 @@ import nipype.interfaces.fsl as fsl -from nipype.interfaces.niftyreg import reg +from nipype.interfaces.niftyreg import reg, regutils -from nipype.interfaces.niftyreg import regutils +from nipype.interfaces.ants.utils import ImageMath from nipype.interfaces.ants.segmentation import DenoiseImage @@ -395,6 +395,48 @@ def create_short_preparation_pipe(params, params_template={}, data_preparation_pipe.connect(av_T2, 'avg_img', align_T2_on_T1, 'flo_file') + # pad image T2 + pad_image_T2 = pe.Node( + ImageMath(), + name="pad_image_T2") + + pad_image_T2.inputs.copy_header = True + pad_image_T2.inputs.operation = "PadImage" + pad_image_T2.inputs.op2 = '200' + + if "avg_reorient_pipe" in params.keys(): + data_preparation_pipe.connect( + av_T2, 'outputnode.std_img', + pad_image_T2, "op1") + + else: + data_preparation_pipe.connect( + av_T2, 'avg_img', + pad_image_T2, "op1") + + # resampling using transfo on much bigger image + reg_resample_T2 = pe.Node( + regutils.RegResample(pad_val=0.0), + name="reg_resample_T2") + + data_preparation_pipe.connect( + align_T2_on_T1, 'aff_file', + reg_resample_T2, 'trans_file') + + data_preparation_pipe.connect( + pad_image_T2, 'output_image', + reg_resample_T2, "flo_file") + + if "avg_reorient_pipe" in params.keys(): + data_preparation_pipe.connect( + av_T1, 'outputnode.std_img', + reg_resample_T2, "ref_file") + + else: + data_preparation_pipe.connect( + av_T1, 'avg_img', + reg_resample_T2, "ref_file") + # align avg T2 on avg T1 else: @@ -429,7 +471,7 @@ def create_short_preparation_pipe(params, params_template={}, if 'aladin_T2_on_T1' in params.keys(): data_preparation_pipe.connect( - align_T2_on_T1, "res_file", + reg_resample_T2, 'out_file', outputnode, 'native_T1') else: data_preparation_pipe.connect( @@ -449,7 +491,7 @@ def create_short_preparation_pipe(params, params_template={}, if 'aladin_T2_on_T1' in params.keys(): data_preparation_pipe.connect( - align_T2_on_T1, "res_file", + reg_resample_T2, 'out_file', outputnode, 'native_T2') else: data_preparation_pipe.connect( @@ -490,8 +532,9 @@ def create_short_preparation_pipe(params, params_template={}, crop_T1, 'in_file') if 'aladin_T2_on_T1' in params.keys(): - data_preparation_pipe.connect(align_T2_on_T1, "res_file", - crop_T2, 'in_file') + data_preparation_pipe.connect( + reg_resample_T2, 'out_file', + crop_T2, 'in_file') else: data_preparation_pipe.connect(align_T2_on_T1, "out_file", @@ -527,8 +570,9 @@ def create_short_preparation_pipe(params, params_template={}, name='apply_crop_aladin_T2') if 'aladin_T2_on_T1' in params.keys(): - data_preparation_pipe.connect(align_T2_on_T1, "res_file", - apply_crop_aladin_T2, 'flo_file') + data_preparation_pipe.connect( + reg_resample_T2, 'out_file', + apply_crop_aladin_T2, 'flo_file') else: data_preparation_pipe.connect(align_T2_on_T1, "out_file", @@ -678,7 +722,7 @@ def create_short_preparation_pipe(params, params_template={}, # T2 to pad if 'aladin_T2_on_T1' in params.keys(): data_preparation_pipe.connect( - align_T2_on_T1, "res_file", + reg_resample_T2, 'out_file', resample_T2_pad, "flo_file") else: @@ -687,6 +731,8 @@ def create_short_preparation_pipe(params, params_template={}, resample_T2_pad, "flo_file") if "padded_template_head" in params_template.keys(): + print("Found padded_template_head in template_stereo") + resample_T1_pad.inputs.ref_file = \ params_template["padded_template_head"] @@ -694,6 +740,10 @@ def create_short_preparation_pipe(params, params_template={}, params_template["padded_template_head"] elif "template_head" in params_template.keys(): + if "pad_val" in params["resample_T1_pad"].keys(): + pad_val = params["resample_T1_pad"]["pad_val"] + print("Found pad_val {}".format(pad_val)) + # padding versio of the template pad_template = NodeParams( niu.Function( diff --git a/macapype/pipelines/register.py b/macapype/pipelines/register.py index 8e98019cc..707f40a46 100644 --- a/macapype/pipelines/register.py +++ b/macapype/pipelines/register.py @@ -590,29 +590,50 @@ def create_crop_aladin_pipe(name="crop_aladin_pipe", params={}): reg_T1_on_template, 'res_file', crop_z_T1, 'in_file') - # remove nans - remove_nans = pe.Node( - fsl.maths.MathsCommand(nan2zeros=True), - name="remove_nans") + pad_image_T1 = pe.Node( + ants.utils.ImageMath(), + name="pad_image_T1") - if "crop_z_T1" in params.keys(): + pad_image_T1.inputs.copy_header = True + pad_image_T1.inputs.operation = "PadImage" + pad_image_T1.inputs.op2 = '200' + + if "pre_crop_z_T1" in params.keys(): reg_pipe.connect( - crop_z_T1, 'out_roi', - remove_nans, "in_file") + pre_crop_z_T1, 'out_roi', + pad_image_T1, "op1") else: - if "reg_T1_on_template2" in params.keys(): - reg_pipe.connect( - reg_T1_on_template2, 'res_file', - remove_nans, "in_file") + reg_pipe.connect( + inputnode, 'native_T1', + pad_image_T1, "op1") - else: - reg_pipe.connect( - reg_T1_on_template, 'res_file', - remove_nans, "in_file") + # resampling using transfo on much bigger image + reg_resample_T1 = pe.Node( + regutils.RegResample(pad_val=0.0), + name="reg_resample_T1") + + # transfo + if "reg_T1_on_template2" in params.keys(): + reg_pipe.connect( + compose_transfo, 'out_file', + reg_resample_T1, 'trans_file') + + else: + reg_pipe.connect( + reg_T1_on_template, 'aff_file', + reg_resample_T1, 'trans_file') + + reg_pipe.connect( + pad_image_T1, 'output_image', + reg_resample_T1, "flo_file") + + reg_pipe.connect( + inputnode, 'stereo_template_T1', + reg_resample_T1, "ref_file") # outputnode - reg_pipe.connect(remove_nans, 'out_file', + reg_pipe.connect(reg_resample_T1, 'out_file', outputnode, "stereo_T1") if "reg_T1_on_template2" in params.keys(): diff --git a/macapype/utils/subs.json b/macapype/utils/subs.json index 737cc0037..88679ab8b 100644 --- a/macapype/utils/subs.json +++ b/macapype/utils/subs.json @@ -24,6 +24,8 @@ "T1w_":"", "T2w_":"", "PDw_":"", + "_aff_":"_", + "_trans":"", "_roi":"", "_ROI":"", "_restore_debiased_brain":"", diff --git a/workflows/segment_pnh.py b/workflows/segment_pnh.py index b584535da..1a4e4d8d4 100644 --- a/workflows/segment_pnh.py +++ b/workflows/segment_pnh.py @@ -282,11 +282,21 @@ def create_main_workflow(cmd, data_dir, process_dir, soft, species, datatypes, assert os.path.exists(template_head), "Could not find template_head {}".format(template_head) params_template["template_head"] = template_head + template_brain = os.path.join(template_path,template_files[1]) assert os.path.exists(template_brain), "Could not find template_brain {}".format(template_brain) params_template["template_brain"] = template_brain - if len(template_files) == 3: + if len(template_files) == 2: + + print("Only two files (template_head and template_brain) have been specified, segmentation will be without priors") + + if "brain_segment_pipe" in params.keys(): + if "segment_atropos_pipe" in params["brain_segment_pipe"].keys(): + if "use_priors" in params["brain_segment_pipe"]["segment_atropos_pipe"].keys(): + del params["brain_segment_pipe"]["segment_atropos_pipe"]["use_priors"] + + elif len(template_files) == 3: template_seg = os.path.join(template_path,template_files[2]) assert os.path.exists(template_seg), "Could not find template_seg {}".format(template_seg)