Merge branch 'master' into refactor_stat

affine_transfo.py

@@ -18,10 +18,6 @@
 # About the license: see the file LICENSE
 ###################################################################
 
-# TODO: add feature to re-use transformation if csv file exists v
-# TODO: clarify what i_dir is. v
-# TODO: save image as FLOAT32 v
-# TODO: add flag interp to apply transfo to seg-labeled data --> in that case, interp=0 (nearestneighbor) v
 # TODO (less priority): check padding (seems unecessary)
 
 import glob, os, sys

config_script.yml

@@ -3,12 +3,11 @@
 # process_data.sh
 # ---------------
 # Number of transformations to apply to each subject
-n_transfo: 3
+n_transfo: 2
 # Rescaling factor to apply to each subject. Note: value '1' should always be present.
 rescaling:
   - 1
   - 0.95
-  - 0.8
 # Contrasts on which to perform the analysis. Available values are: 't1', 't2'
 contrast: "t1"
 transfo_file: transfo_values.csv

config_sct_run_batch.yml

@@ -2,7 +2,7 @@
 path_data: spinegeneric_r20200801
 path_output: csa_atrophy_results
 script: process_data.sh
-script_args: --config config_script.yml
+script_args: config_script.yml
 jobs: -1
 batch_log: sct_run_batch_log.txt
 continue_on_error: 1

manual_correction.py

@@ -106,7 +106,7 @@ def get_suffix(task, suffix=''):
     elif task == 'FILES_GMSEG':
         return '_gmseg'+suffix
     elif task == 'FILES_LABEL':
-        return '_labels'+suffix
+        return '_labels-disc'+suffix
     else:
         raise ValueError("This task is not recognized: {}".format(task))
 
@@ -150,20 +150,6 @@ def create_json(fname_nifti, name_rater):
         json.dump(metadata, outfile, indent=4)
 
 
-def get_rescaling(file):
-    """
-    Get subject from BIDS file name
-    :param file:
-    :return: rescaling
-    """
-    if "_r" in file:
-        rescale = "_"+file.split('_')[2]
-        return rescale
-    else:
-        rescale = ""
-        return rescale
-
-
 def check_files_exist(dict_files, path_data):
     """
     Check if all files listed in the input dictionary exist
@@ -175,7 +161,7 @@ def check_files_exist(dict_files, path_data):
     for task, files in dict_files.items():
         if files is not None:
             for file in files:
-                fname = os.path.join(path_data, sg.bids.get_subject(file), sg.bids.get_contrast(file) + get_rescaling(file) , file)
+                fname = os.path.join(path_data, sg.bids.get_subject(file), sg.bids.get_contrast(file), file)
                 if not os.path.exists(fname):
                     missing_files.append(fname)
     if missing_files:
@@ -235,7 +221,7 @@ def main():
                 # build file names
                 subject = sg.bids.get_subject(file)
                 contrast = sg.bids.get_contrast(file)
-                fname = os.path.join(args.path_in, subject, contrast + get_rescaling(file), file)
+                fname = os.path.join(args.path_in, subject, contrast, file)
                 fname_label = os.path.join(
                     path_out_deriv, subject, contrast, sg.utils.add_suffix(file, get_suffix(task, '-manual')))
                 os.makedirs(os.path.join(path_out_deriv, subject, contrast), exist_ok=True)

process_data.sh

@@ -11,7 +11,7 @@
 # Author: Julien Cohen-Adad, Paul Bautin
 ###################################################
 
-# TODO: crop data to save computation time
+# TODO: simplify variable names: no need to call file_or_c_blablabla. Instead: update file
 
 # Uncomment for full verbose
 set -x
@@ -23,6 +23,7 @@ set -e
 trap "echo Caught Keyboard Interrupt within script. Exiting now.; exit" INT
 
 # get starting time:
+start_all=`date +%s`
 start=`date +%s`
 
 # Global variables
@@ -37,18 +38,13 @@ rescaling=$(yaml_parser -o rescaling -i $config_script)
 R_COEFS=$(echo $rescaling | tr '[]' ' ' | tr ',' ' ' | tr "'" ' ')
 contrast=$(yaml_parser -o contrast -i $config_script)
 # TODO: enable to input a list of contrast and loop across contrasts
-if [ $contrast == "t2" ]; then
-  contrast_str="T2w"
-fi
-if [ $contrast == "t1" ]; then
-  contrast_str="T1w"
-fi
 transfo_file=$(yaml_parser -o transfo_file -i $config_script)
 echo "transfo_file: $transfo_file"
 
 
 # FUNCTIONS
 # ==============================================================================
+
 # Check if manual label already exists. If it does, copy it locally. If it does
 # not, perform automatic labeling.
 label_if_does_not_exist(){
@@ -57,38 +53,40 @@ label_if_does_not_exist(){
   local contrast="$3"
   local contrast_str="$4"
   # Update global variable with segmentation file name
-  FILELABEL="${file}_labels"
-  FILELABELMANUAL="${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${FILELABEL}-manual.nii.gz"
-  if [ -e "$FILELABELMANUAL" ]; then
-    echo "manual labeled file was found: $FILELABELMANUAL"
+  FILELABEL="${file}_labels-disc"
+  FILELABELMANUAL="${path_derivatives}/${SUBJECT}_${contrast_str}_labels-disc-manual.nii.gz"
+  if [ -e "${FILELABELMANUAL}" ]; then
+    echo "manual labeled file was found: ${FILELABELMANUAL}"
+    # reorienting and resampling image
+    sct_image -i ${FILELABELMANUAL} -setorient RPI -o "${path_derivatives}/${SUBJECT}_${contrast_str}_RPI_labels-disc-manual.nii.gz"
+    sct_resample -i ${FILELABELMANUAL} -mm $interp -o "${path_derivatives}/${SUBJECT}_${contrast_str}_RPI_r_labels-disc-manual.nii.gz"
     rsync -avzh $FILELABELMANUAL ${FILELABEL}.nii.gz
     # Generate labeled segmentation
     sct_label_vertebrae -i ${file}.nii.gz -s ${file_seg}.nii.gz -c ${contrast} -discfile ${FILELABELMANUAL} -qc ${PATH_QC} -qc-subject ${SUBJECT}
-
-    sct_label_utils -i ${file_seg}_labeled.nii.gz -vert-body 0 -o ${FILELABEL}.nii.gz
   else
     # Generate labeled segmentation
     sct_label_vertebrae -i ${file}.nii.gz -s ${file_seg}.nii.gz -c ${contrast} -qc ${PATH_QC} -qc-subject ${SUBJECT}
-    # Create labels in the Spinal Cord
-    sct_label_utils -i ${file_seg}_labeled.nii.gz -vert-body 0 -o ${FILELABEL}.nii.gz -qc ${PATH_QC} -qc-subject ${SUBJECT}
   fi
+  # Create labels in the Spinal Cord
+  sct_label_utils -i ${file_seg}_labeled.nii.gz -vert-body 0 -o ${FILELABEL}.nii.gz
 }
 
 # Check if manual segmentation already exists. If it does, copy it locally. If
 # it does not, perform seg.
 segment_if_does_not_exist(){
   local file="$1"
   local contrast="$2"
+  local qc=$3
   # Update global variable with segmentation file name
   FILESEG="${file}_seg"
-  FILESEGMANUAL="${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${FILESEG}-manual.nii.gz"
+  FILESEGMANUAL="${path_derivatives}/${FILESEG}-manual.nii.gz"
   if [ -e $FILESEGMANUAL ]; then
     echo "Found! Using manual segmentation."
     rsync -avzh $FILESEGMANUAL ${FILESEG}.nii.gz
-    sct_qc -i ${file}.nii.gz -s ${FILESEG}.nii.gz -p sct_deepseg_sc -qc ${PATH_QC} -qc-subject ${SUBJECT}
+    sct_qc -i ${file}.nii.gz -s ${FILESEG}.nii.gz -p sct_deepseg_sc $qc
   else
     # Segment spinal cord
-    sct_deepseg_sc -i ${file}.nii.gz -c ${contrast} -qc ${PATH_QC} -qc-subject ${SUBJECT}
+    sct_deepseg_sc -i ${file}.nii.gz -c ${contrast} $qc
   fi
 }
 
@@ -97,8 +95,10 @@ segment_if_does_not_exist(){
 # ==============================================================================
 # Display useful info for the log, such as SCT version, RAM and CPU cores available
 sct_check_dependencies -short
-# Copy config files to output results folder
-cp $config_script ${PATH_RESULTS}/
+# copy derivatives directory containing manual corrections to PATH_DATA_PROCESSED
+mkdir -p "${PATH_DATA_PROCESSED}/derivatives/labels/${SUBJECT}/anat/"
+cp -r "${PATH_DATA}/derivatives/labels/${SUBJECT}/anat" "${PATH_DATA_PROCESSED}/derivatives/labels/${SUBJECT}"
+path_derivatives="${PATH_DATA_PROCESSED}/derivatives/labels/${SUBJECT}/anat"
 # Go to results folder, where most of the outputs will be located
 cd $PATH_DATA_PROCESSED
 # Copy source images
@@ -107,31 +107,78 @@ cd $SUBJECT
 # we don't need dwi data, so let's remove it
 rm -r dwi
 
-# Image analysis
+
+# Image analysis in folder anat
 #=======================================================================
-# Segment spinal cord (only if it does not exist) in dir anat
 cd anat
-echo "contrast: $contrast"
-file_c=${SUBJECT}_${contrast_str}
-segment_if_does_not_exist $file_c ${contrast}
+# Reorient to RPI and resample file
+if [ $contrast == "t2" ]; then
+  contrast_str="T2w"
+  interp="0.8x0.8x0.8"
+elif [ $contrast == "t1" ]; then
+  contrast_str="T1w"
+  interp="1x1x1"
+fi
+file=${SUBJECT}_${contrast_str}
+# Reorient image to RPI
+sct_image -i ${file}.nii.gz -setorient RPI -o ${file}_RPI.nii.gz
+file=${file}_RPI
+# Resample image isotropically
+sct_resample -i ${file}.nii.gz -mm $interp -o ${file}_r.nii.gz
+file=${file}_r
+end=`date +%s`
+runtime=$((end-start))
+echo "+++++++++++ TIME: Duration of reorienting and resampling:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
+
+# Segment spinal cord (only if it does not exist) in dir anat
+start=`date +%s`
+segment_if_does_not_exist $file ${contrast} "-qc ${PATH_QC} -qc-subject ${SUBJECT}"
 # name segmented file
-file_c_seg=${FILESEG}
+file_seg=${FILESEG}
+end=`date +%s`
+runtime=$((end-start))
+echo "+++++++++++ TIME: Duration of segmentation on non cropped image:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
+
 # Label spinal cord (only if it does not exist) in dir anat
-label_if_does_not_exist $file_c $file_c_seg $contrast $contrast_str
-file_label=${file_c_seg}_labeled
+start=`date +%s`
+label_if_does_not_exist $file $file_seg $contrast $contrast_str
+file_label=${file_seg}_labeled
+end=`date +%s`
+runtime=$((end-start))
+echo "+++++++++++ TIME: Duration of labelling:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
+
+# dilate segmentation (for cropping)
+sct_maths -i ${file_seg}.nii.gz -dilate 15 -shape cube -o ${file_seg}_dil.nii.gz
+# crop image
+sct_crop_image -i ${file}.nii.gz -m ${file_seg}_dil.nii.gz
+file=${file}_crop
+# crop segmentation
+sct_crop_image -i ${file_seg}.nii.gz -m ${file_seg}_dil.nii.gz
+file_seg=${file_seg}_crop
 cd ../
 
-# iterate across rescaling
+
+# CSA measure iterated across rescaling in folder anat_r{r_coef}
+# ====================================================================
 for r_coef in ${R_COEFS[@]}; do
+  # If directory exists (e.g. 2nd pass after QC and manual correction), we should remove it
+  if [ -d "anat_r${r_coef}" ]; then
+    rm -r "anat_r${r_coef}"
+    echo "anat_r${r_coef} already exists: removing folder"
+  fi
+  if [ -f "$PATH_RESULTS/csa_perlevel_${SUBJECT}_${r_coef}.csv" ]; then
+    rm "$PATH_RESULTS/csa_perlevel_${SUBJECT}_${r_coef}.csv"
+    echo "csa_perlevel_${SUBJECT}_${r_coef}.csv already exists: removing file"
+  fi
   mkdir anat_r$r_coef
   cd anat_r${r_coef}
 
   # Rescale header of native nifti file
-  file_c_r=${file_c}_r${r_coef}
-  affine_rescale -i ../anat/${file_c}.nii.gz -r ${r_coef} -o ${file_c_r}.nii.gz
+  file_r=${file}_r${r_coef}
+  affine_rescale -i ../anat/${file}.nii.gz -r ${r_coef} -o ${file_r}.nii.gz
   # rescale labeled segmentation
-  file_label_c_r=${file_c_r}_seg_labeled
-  affine_rescale -i ../anat/${file_label}.nii.gz -r ${r_coef} -o ${file_label_c_r}.nii.gz
+  file_label_r=${file_r}_seg_labeled
+  affine_rescale -i ../anat/${file_label}.nii.gz -r ${r_coef} -o ${file_label_r}.nii.gz
 
   # create list of array to iterate on (e.g.: seq_transfo = 1 2 3 4 5 if n_transfo=5)
   seq_transfo=$(seq ${n_transfo})
@@ -141,21 +188,40 @@ for r_coef in ${R_COEFS[@]}; do
     # "transfo_values.csv" file if it already exists.
     # We keep a transfo_values.csv file, so that after first pass of the pipeline and QC, if segmentations
     # need to be manually-corrected, we want the transformations to be the same for the 2nd pass of the pipeline.
-    affine_transfo -i ${file_c_r}.nii.gz -transfo ${PATH_RESULTS}/$transfo_file -config ${PATH_RESULTS}/$config_script -o _t${i_transfo}
-    file_c_r_t=${file_c_r}_t${i_transfo}
+    start=`date +%s`
+    affine_transfo -i ${file_r}.nii.gz -transfo ${PATH_RESULTS}/transfo_${file_r} -config $config_script -o _t${i_transfo}
+    file_r_t=${file_r}_t${i_transfo}
+    end=`date +%s`
+    runtime=$((end-start))
+    echo "+++++++++++ TIME: Duration of of image transfo t${i_transfo}:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
     # transform the labeled segmentation with same transfo values
-    affine_transfo -i ${file_label_c_r}.nii.gz -transfo ${PATH_RESULTS}/$transfo_file -config ${PATH_RESULTS}/$config_script -o _t${i_transfo}_seg_labeled -interpolation 0
-    file_label_c_r_t=${file_label_c_r}_t${i_transfo}
+
+    start=`date +%s`
+    affine_transfo -i ${file_label_r}.nii.gz -transfo ${PATH_RESULTS}/transfo_${file_r} -config $config_script -o _t${i_transfo} -interpolation 0
+    file_label_r_t=${file_label_r}_t${i_transfo}
+    end=`date +%s`
+    runtime=$((end-start))
+    echo "+++++++++++ TIME: Duration of labelling transfo t${i_transfo}:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
     # Segment spinal cord (only if it does not exist)
-    segment_if_does_not_exist ${file_c_r_t} ${contrast}
+    start=`date +%s`
+    segment_if_does_not_exist ${file_r_t} ${contrast}
+    end=`date +%s`
+    runtime=$((end-start))
+    echo "+++++++++++ TIME: Duration of segmentation t${i_transfo}:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
     # name segmented file
-    file_c_r_t_seg=${FILESEG}
+    file_r_t_seg=${FILESEG}
     # Compute average CSA between C2 and C5 levels (append across subjects)
-    sct_process_segmentation -i $file_c_r_t_seg.nii.gz -vert 2:5 -perlevel 1 -vertfile $file_label_c_r_t.nii.gz -o $PATH_RESULTS/csa_perlevel_${SUBJECT}_t${i_transfo}_${r_coef}.csv -qc ${PATH_QC}
+    start=`date +%s`
+    sct_process_segmentation -i $file_r_t_seg.nii.gz -vert 3:5 -perlevel 1 -vertfile $file_label_r_t.nii.gz -o $PATH_RESULTS/csa_perlevel_${SUBJECT}_t${i_transfo}_${r_coef}.csv
+    end=`date +%s`
+    runtime=$((end-start))
+    echo "+++++++++++ TIME: Duration of sct_process_segmentation t${i_transfo}:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
     # add files to check
     FILES_TO_CHECK+=(
-    "$PATH_RESULTS/csa_data/csa_perlevel_${SUBJECT}_t${i_transfo}_${r_coef}.csv"
-    "$PATH_DATA_PROCESSED/${SUBJECT}/anat_r${r_coef}/${file_c_seg}.nii.gz"
+    "$PATH_RESULTS/csa_perlevel_${SUBJECT}_t${i_transfo}_${r_coef}.csv"
+    "$PATH_DATA_PROCESSED/${SUBJECT}/anat_r${r_coef}/${file_r_t}.nii.gz"
+    "$PATH_DATA_PROCESSED/${SUBJECT}/anat_r${r_coef}/${file_r_t_seg}.nii.gz"
+    "$PATH_DATA_PROCESSED/${SUBJECT}/anat_r${r_coef}/${file_label_r_t}.nii.gz"
     )
   done
   cd ../
@@ -172,10 +238,10 @@ done
 # Display useful info for the log
 # ===============================
 end=`date +%s`
-runtime=$((end-start))
+runtime=$((end-start_all))
 echo
 echo "~~~"
 echo "SCT version: `sct_version`"
 echo "Ran on:      `uname -nsr`"
 echo "Duration:    $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
-echo "~~~"
+echo "~~~"

requirements.txt

@@ -5,3 +5,9 @@ matplotlib
 scikit-image
 nibabel
 ruamel.yaml
+yaml-1.3
+argparse~=1.4.0
+setuptools~=49.6.0
+PyYAML~=5.3.1
+coloredlogs~=14.0
+