diff --git a/process_data.sh b/process_data.sh index 09f195e6..33d728fd 100644 --- a/process_data.sh +++ b/process_data.sh @@ -307,14 +307,19 @@ for file_path in "${inc_contrasts[@]}";do # Find contrast to do segmentation if [[ $file_path == *"T1w"* ]];then contrast_seg="t1" + contrast="t1w" elif [[ $file_path == *"T2star"* ]];then contrast_seg="t2s" + contrast=contrast_seg elif [[ $file_path == *"T1w_MTS"* ]];then contrast_seg="t1" + contrast="T1w_MTS" elif [[ $file_path == *"MTon_MTS"* ]];then contrast_seg="t2s" + contrast="MTon_MTS" elif [[ $file_path == *"dwi"* ]];then contrast_seg="dwi" + contrast=contrast_seg fi type=$(find_contrast $file_path) @@ -326,7 +331,7 @@ for file_path in "${inc_contrasts[@]}";do python ${PATH_SCRIPT}/pad_seg.py -i ${fileseg}.nii.gz -o ${fileseg}_pad.nii.gz # Registration # ------------------------------------------------------------------------------ - sct_register_multimodal -i ${file_path}.nii.gz -d ./anat/${file_t2}.nii.gz -iseg ${fileseg}_pad.nii.gz -dseg ./anat/${file_t2_seg}.nii.gz -param step=1,type=seg,algo=slicereg,metric=MeanSquares,iter=10,poly=2 -qc ${PATH_QC} -qc-subject ${SUBJECT} -o ${file_path}_reg.nii.gz + sct_register_multimodal -i ${file_path}.nii.gz -d ./anat/${file_t2}.nii.gz -iseg ${fileseg}_pad.nii.gz -dseg ./anat/${file_t2_seg}.nii.gz -param step=1,type=seg,algo=centermass -qc ${PATH_QC} -qc-subject ${SUBJECT} -o ${file_path}_reg.nii.gz warping_field=${type}warp_${file}2${file_t2} # Generate SC segmentation coverage and register to T2w @@ -336,6 +341,13 @@ for file_path in "${inc_contrasts[@]}";do sct_apply_transfo -i ${file_path}_ones.nii.gz -d ./anat/${file_t2}.nii.gz -w ${warping_field}.nii.gz -x linear -o ${file_path}_ones_reg.nii.gz # Bring SC segmentation to T2w space sct_apply_transfo -i ${fileseg}.nii.gz -d ./anat/${file_t2_seg}.nii.gz -w ${warping_field}.nii.gz -x linear -o ${fileseg}_reg.nii.gz + # Remove 8 to 10 slices before adding all segmentation because of partial slices (except for T1w) + if [[ $contrast != "t1w" ]];then + mv ${fileseg}_reg.nii.gz ${fileseg}_reg_no_crop.nii.gz + mv ${file_path}_ones_reg.nii.gz ${file_path}_ones_reg_no_crop.nii.gz + python ${PATH_SCRIPT}/remove_slices_seg.py -i ${fileseg}_reg_no_crop.nii.gz -c contrast -coverage-map ${file_path}_ones_reg_no_crop.nii.gz -o ${fileseg}_reg.nii.gz -o-coverage-map ${file_path}_ones_reg.nii.gz + fi + done # Create coverage mask for T2w sct_create_mask -i ./anat/${file_t2}.nii.gz -o ./anat/${file_t2}_ones.nii.gz -size 500 -p centerline,./anat/${file_t2_seg}.nii.gz diff --git a/remove_slices_seg.py b/remove_slices_seg.py new file mode 100644 index 00000000..598da3ef --- /dev/null +++ b/remove_slices_seg.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python +# -*- coding: utf-8 +# Removes top and bottom slices of segmentation and coverage map +# +# For usage, type: python remove_slices_seg.py -h + +# Authors: Sandrine Bédard + +import argparse +import numpy as np +import nibabel as nib + +NEAR_ZERO_THRESHOLD = 1e-6 + + +def get_parser(): + parser = argparse.ArgumentParser( + description="Removes top and bottom slices of segmentation and coverage map to avoid propagation of registration errors when summing segmentations and coverage maps.") + parser.add_argument('-i', required=True, type=str, + help="Input segmentation in T2w space.") + parser.add_argument('-coverage-map', required=True, type=str, + help="Coverage map in T2w space.") + parser.add_argument('-o', required=True, type=str, + help="Output name of crop segmentation") + parser.add_argument('-o-coverage-map', required=True, type=str, + help="Output name of crop coverage map") + parser.add_argument('-c', required=True, type=str, + help="Contrast.") + + return parser + + +def save_Nifti1(data, original_image, filename): + empty_header = nib.Nifti1Header() + image = nib.Nifti1Image(data, original_image.affine, empty_header) + nib.save(image, filename) + + +def remove_slices(image, max_z_index, min_z_index, nb_slices): + image_seg_crop = image.copy() + # Remove top slices + image_seg_crop[:, :, (max_z_index - nb_slices)::] = 0 + # Remove bottom slice + image_seg_crop[:, :, 0:(min_z_index + nb_slices + 1)] = 0 + return image_seg_crop + + +def main(): + parser = get_parser() + args = parser.parse_args() + + contrast_seg = args.c + image = nib.load(args.i) + image_np = image.get_fdata() + coverage_map = nib.load(args.coverage_map) + coverage_map_np = coverage_map.get_fdata() + + # Get max and min index of the segmentation + _, _, Z = (image_np > NEAR_ZERO_THRESHOLD).nonzero() + # Find index of top and bottom slice of the segmentation + min_z_index, max_z_index = min(Z), max(Z) + + # Number of slices to remove: + # MTS: 10 to and bottom + # DWI & T2star : 8 slices + if contrast_seg == "t2s" or contrast_seg == "dwi": + nb_slices = 7 + else: + nb_slices = 9 + image_np_crop = remove_slices(image_np, max_z_index, min_z_index, nb_slices) + coverage_map_np_crop = remove_slices(coverage_map_np, max_z_index, min_z_index, nb_slices) + + # save new nifti for both segmentation and coverage maps + filename_seg_crop = args.o + filename_coverage_crop = args.o_coverage_map + save_Nifti1(image_np_crop, image, filename_seg_crop) + save_Nifti1(coverage_map_np_crop, coverage_map, filename_coverage_crop) + + +if __name__ == '__main__': + main() \ No newline at end of file