Merge pull request #32 from sct-pipeline/jca/modif-reg-algo

Changed intercontrast registration algo

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

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()