Merge pull request #18 from ivadomed/jv/refactor_utilities_folder

Move scripts from `utilities` to `inter-rater_variability`

inter-rater_variability/01_combine_segmentations_from_different_raters.py

@@ -0,0 +1,156 @@
+#
+# Combine several multi-class (i.e., not binary) segmentations (obtained manually by expert raters) into a reference
+# segmentation using the STAPLE algorithm.
+#
+# Note: since the segmentations are multi-class (i.e., not binary), we need to binaries them. We do this for each
+# spinal level separately.
+#
+# Example:
+#   python 01_combine_segmentations_from_different_raters.py
+#       -i  sub-001_T2w_label-rootlet_rater1.nii.gz
+#           sub-001_T2w_label-rootlet_rater2.nii.gz
+#           sub-001_T2w_label-rootlet_rater3.nii.gz
+#           sub-001_T2w_label-rootlet_rater4.nii.gz
+#       -o  sub-001_T2w_label-rootlet_staple.nii.gz
+#
+# Authors: Jan Valosek
+#
+
+import os
+import argparse
+
+import numpy as np
+import SimpleITK as sitk
+
+
+def get_parser():
+    """
+    parser function
+    """
+
+    parser = argparse.ArgumentParser(
+        description='Combine several multi-class (i.e., not binary) segmentations (obtained manually by expert raters) '
+                    'into a reference segmentation using the STAPLE algorithm.',
+        prog=os.path.basename(__file__).strip('.py')
+    )
+    parser.add_argument(
+        '-i',
+        required=True,
+        nargs='+',
+        help='Paths to the segmentation files separated by spaces. Example: sub-001_T2w_label-rootlet_rater1.nii.gz '
+             'sub-001_T2w_label-rootlet_rater2.nii.gz sub-001_T2w_label-rootlet_rater3.nii.gz'
+    )
+    parser.add_argument(
+        '-o',
+        required=True,
+        type=str,
+        help='Path to the output. Example: sub-001_T2w_label-rootlet_staple.nii.gz'
+    )
+
+    return parser
+
+
+def combine_staple(segmentations, fname_out):
+    """
+    Combine several segmentations into a reference segmentation using the STAPLE algorithm.
+
+    Note: since the segmentations are multi-class (i.e., not binary), we need to binaries them. We do this for each
+    spinal level separately.
+
+    Inspiration: https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1STAPLEImageFilter.html#details
+
+    :param segmentations: list of segmentations
+    :param fname_out: output file name
+    Returns:
+    """
+
+    # Initialize the list for reference segmentations for each level
+    reference_segmentations_all_levels = list()
+
+    # Loop across spinal levels
+    for level in np.unique(sitk.GetArrayFromImage(segmentations[0])):
+        # Skip zero (background)
+        if level == 0:
+            continue
+
+        print(f'Processing level {level} ...')
+        segmentations_current_level = list()
+        # Binarize the segmentations for a specific level
+        for i in range(len(segmentations)):
+            segmentations_current_level.append(sitk.BinaryThreshold(
+                    segmentations[i],
+                    lowerThreshold=int(level),
+                    upperThreshold=int(level),
+                    insideValue=1,
+                    outsideValue=0
+                )
+            )
+
+        # Combine the segmentations for a specific level using STAPLE
+        foregroundValue = 1
+        threshold = 0.95
+        reference_segmentation_STAPLE_probabilities = sitk.STAPLE(
+            segmentations_current_level, foregroundValue
+        )
+
+        # Threshold the probabilities
+        reference_segmentation_current_level = reference_segmentation_STAPLE_probabilities > threshold
+
+        # Change the reference_segmentation_current_level value to the current level to be able to recreate the
+        # multi-class segmentation
+        reference_segmentation_current_level = sitk.Cast(
+            reference_segmentation_current_level,
+            sitk.sitkUInt8
+        ) * level
+
+        # Store the combined segmentation for each level into a list
+        reference_segmentations_all_levels.append(reference_segmentation_current_level)
+
+    # Now, combine the segmentations for each level back into one multi-class file
+    # Initialize the final segmentation
+    final_segmentation = sitk.Image(
+        segmentations[0].GetSize(),
+        sitk.sitkUInt8
+    )
+    final_segmentation.CopyInformation(segmentations[0])
+
+    for i in range(len(reference_segmentations_all_levels)):
+        final_segmentation = sitk.Add(
+            final_segmentation,
+            reference_segmentations_all_levels[i]
+        )
+
+    # Save the reference segmentation
+    sitk.WriteImage(
+        final_segmentation,
+        fname_out
+    )
+    print(f'Reference segmentation saved as {fname_out}.')
+
+
+def main():
+    # Parse the command line arguments
+    parser = get_parser()
+    args = parser.parse_args()
+
+    # Parse paths
+    full_paths = [os.path.join(os.getcwd(), path) for path in args.i]
+
+    # Check if the input path exists
+    for fname in full_paths:
+        if not os.path.exists(fname):
+            raise ValueError('Input segmentation does not exist: {}'.format(fname))
+
+    print('Input segmentations:')
+    print('\n'.join(full_paths))
+
+    segmentations = [
+        sitk.ReadImage(file_name, sitk.sitkUInt8)
+        for file_name in full_paths
+    ]
+
+    combine_staple(segmentations, args.o)
+
+
+if __name__ == '__main__':
+    main()

inter-rater_variability/02_run_batch_inter_rater_variability.sh

@@ -0,0 +1,184 @@
+#!/bin/bash
+#
+# The script performs inter-rater variability analysis
+# Namely:
+#   - segment spinal cord
+#   - detect PMJ label
+#   - run 02a_rootlets_to_spinal_levels.py to project the nerve rootlets on the spinal cord segmentation to obtain spinal
+#   levels, and compute the distance between the pontomedullary junction (PMJ) and the start and end of the spinal
+#   level. The 02a_rootlets_to_spinal_levels.py script outputs .nii.gz file with spinal levels and saves the results in CSV files.
+#   - run 02b_compute_f1_and_dice.py to compute the F1 and Dice scores between the reference and GT segmentations. The
+#   02b_compute_f1_and_dice.py saves the results in CSV files.
+#
+# Usage:
+#       sct_run_batch -script 02_run_batch_inter_rater_variability.sh -path-data <DATA> -path-output <DATA>_202X-XX-XX -jobs 16 -script-args <PATH_TO_PYTHON_SCRIPTS>
+#
+# The following global variables are retrieved from the caller sct_run_batch
+# but could be overwritten by uncommenting the lines below:
+# PATH_DATA_PROCESSED="~/data_processed"
+# PATH_RESULTS="~/results"
+# PATH_LOG="~/log"
+# PATH_QC="~/qc"
+#
+# Authors: Jan Valosek
+#
+
+# Uncomment for full verbose
+set -x
+
+# Immediately exit if error
+set -e -o pipefail
+
+# Exit if user presses CTRL+C (Linux) or CMD+C (OSX)
+trap "echo Caught Keyboard Interrupt within script. Exiting now.; exit" INT
+
+# Print retrieved variables from the sct_run_batch script to the log (to allow easier debug)
+echo "Retrieved variables from from the caller sct_run_batch:"
+echo "PATH_DATA: ${PATH_DATA}"
+echo "PATH_DATA_PROCESSED: ${PATH_DATA_PROCESSED}"
+echo "PATH_RESULTS: ${PATH_RESULTS}"
+echo "PATH_LOG: ${PATH_LOG}"
+echo "PATH_QC: ${PATH_QC}"
+
+# ------------------------------------------------------------------------------
+# CONVENIENCE FUNCTIONS
+# ------------------------------------------------------------------------------
+# Check if manual spinal cord segmentation file already exists. If it does, copy it locally.
+# If it doesn't, perform automatic spinal cord segmentation
+segment_if_does_not_exist() {
+  local file="$1"
+  local contrast="$2"
+  # Update global variable with segmentation file name
+  FILESEG="${file}_seg"
+  FILESEGMANUAL="${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${FILESEG}-manual.nii.gz"
+  echo
+  echo "Looking for manual segmentation: $FILESEGMANUAL"
+  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}
+  else
+    echo "Not found. Proceeding with automatic segmentation."
+    # Segment spinal cord
+    sct_deepseg_sc -i ${file}.nii.gz -o ${FILESEG}.nii.gz -c ${contrast} -qc ${PATH_QC} -qc-subject ${SUBJECT}
+  fi
+}
+
+# Check if manual PMJ label file already exists. If it does, copy it locally.
+# If it doesn't, perform automatic PMJ labeling
+detect_pmj() {
+  local file="$1"
+  local contrast="$2"
+  # Update global variable with segmentation file name
+  FILESEG="${file}_pmj"
+  FILESEGMANUAL="${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${FILESEG}-manual.nii.gz"
+  echo
+  echo "Looking for manual PMJ label: $FILESEGMANUAL"
+  if [[ -e $FILESEGMANUAL ]]; then
+    echo "Found! Using manual PMJ label."
+    rsync -avzh $FILESEGMANUAL ${FILESEG}.nii.gz
+    sct_qc -i ${file}.nii.gz -s ${FILESEG}.nii.gz -p sct_label_pmj -qc ${PATH_QC} -qc-subject ${SUBJECT}
+  else
+    echo "Not found. Proceeding with automatic segmentation."
+    # Segment spinal cord
+    sct_detect_pmj -i ${file}.nii.gz -o ${FILESEG}.nii.gz -c ${contrast} -qc ${PATH_QC} -qc-subject ${SUBJECT}
+  fi
+}
+
+# Copy GT segmentation
+copy_gt(){
+  local file="$1"
+  local suffix="$2"
+  # Construct file name to GT segmentation located under derivatives/labels
+  FILESEGMANUAL="${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${file}_${suffix}.nii.gz"
+  echo ""
+  echo "Looking for manual segmentation: $FILESEGMANUAL"
+  if [[ -e $FILESEGMANUAL ]]; then
+      echo "Found! Copying ..."
+      rsync -avzh $FILESEGMANUAL ${file}_${suffix}.nii.gz
+  else
+      echo "File ${FILESEGMANUAL} does not exist" >> ${PATH_LOG}/missing_files.log
+      echo "ERROR: Manual GT segmentation ${FILESEGMANUAL} does not exist. Exiting."
+      exit 1
+  fi
+}
+
+
+# Retrieve input params and other params
+SUBJECT=$1
+# Path to the directory with 02a_rootlets_to_spinal_levels.py and 02b_compute_f1_and_dice.py scripts
+PATH_PYTHON_SCRIPTS=$2
+
+# get starting time:
+start=`date +%s`
+
+# ------------------------------------------------------------------------------
+# SCRIPT STARTS HERE
+# ------------------------------------------------------------------------------
+# Display useful info for the log, such as SCT version, RAM and CPU cores available
+sct_check_dependencies -short
+
+# Go to folder where data will be copied and processed
+cd $PATH_DATA_PROCESSED
+
+# Copy source images
+# Note: we use '/./' in order to include the sub-folder 'ses-0X'
+rsync -Ravzh $PATH_DATA/./$SUBJECT .
+
+# Go to subject folder for source images
+cd ${SUBJECT}/anat
+
+# Define variables
+# We do a substitution '/' --> '_' in case there is a subfolder 'ses-0X/'
+file="${SUBJECT//[\/]/_}"
+
+# -------------------------------------------------------------------------
+# T2w
+# -------------------------------------------------------------------------
+# Add suffix corresponding to contrast
+file_t2w=${file}_T2w
+# Check if T2w image exists
+if [[ -f ${file_t2w}.nii.gz ]];then
+
+    # Segment spinal cord
+    segment_if_does_not_exist $file_t2w t2
+
+    # Label PMJ
+    detect_pmj $file_t2w t2
+
+    # Copy manual GT rootlets segmentation for each rater
+    copy_gt $file_t2w label-rootlet_rater1
+    copy_gt $file_t2w label-rootlet_rater2
+    copy_gt $file_t2w label-rootlet_rater3
+    copy_gt $file_t2w label-rootlet_rater4
+
+    # Project the nerve rootlets on the spinal cord segmentation to obtain spinal levels and compute the distance
+    # between the pontomedullary junction (PMJ) and the start and end of the spinal level
+    python ${PATH_PYTHON_SCRIPTS}/02a_rootlets_to_spinal_levels.py -i ${file_t2w}_label-rootlet_rater1.nii.gz -s ${file_t2w}_seg.nii.gz -pmj ${file_t2w}_pmj.nii.gz
+    python ${PATH_PYTHON_SCRIPTS}/02a_rootlets_to_spinal_levels.py -i ${file_t2w}_label-rootlet_rater2.nii.gz -s ${file_t2w}_seg.nii.gz -pmj ${file_t2w}_pmj.nii.gz
+    python ${PATH_PYTHON_SCRIPTS}/02a_rootlets_to_spinal_levels.py -i ${file_t2w}_label-rootlet_rater3.nii.gz -s ${file_t2w}_seg.nii.gz -pmj ${file_t2w}_pmj.nii.gz
+    python ${PATH_PYTHON_SCRIPTS}/02a_rootlets_to_spinal_levels.py -i ${file_t2w}_label-rootlet_rater4.nii.gz -s ${file_t2w}_seg.nii.gz -pmj ${file_t2w}_pmj.nii.gz
+
+    # Copy the reference segmentation generated using the STAPLE algorithm
+    copy_gt $file_t2w label-rootlet_staple
+
+    # Compute f1 and dice scores for each level between the reference segmentation and GT segmentations
+    python ${PATH_PYTHON_SCRIPTS}/02b_compute_f1_and_dice.py -gt ${file_t2w}_label-rootlet_staple.nii.gz -pr ${file_t2w}_label-rootlet_rater1.nii.gz -im ${file_t2w}.nii.gz -o ${file_t2w}_label-rootlet_rater1
+    python ${PATH_PYTHON_SCRIPTS}/02b_compute_f1_and_dice.py -gt ${file_t2w}_label-rootlet_staple.nii.gz -pr ${file_t2w}_label-rootlet_rater2.nii.gz -im ${file_t2w}.nii.gz -o ${file_t2w}_label-rootlet_rater2
+    python ${PATH_PYTHON_SCRIPTS}/02b_compute_f1_and_dice.py -gt ${file_t2w}_label-rootlet_staple.nii.gz -pr ${file_t2w}_label-rootlet_rater3.nii.gz -im ${file_t2w}.nii.gz -o ${file_t2w}_label-rootlet_rater3
+    python ${PATH_PYTHON_SCRIPTS}/02b_compute_f1_and_dice.py -gt ${file_t2w}_label-rootlet_staple.nii.gz -pr ${file_t2w}_label-rootlet_rater4.nii.gz -im ${file_t2w}.nii.gz -o ${file_t2w}_label-rootlet_rater4
+
+fi
+# ------------------------------------------------------------------------------
+# End
+# ------------------------------------------------------------------------------
+
+# Display results (to easily compare integrity across SCT versions)
+end=`date +%s`
+runtime=$((end-start))
+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 "~~~"