First prototype of processing script

batch_processing.sh

@@ -0,0 +1,201 @@
+#!/bin/bash
+# NB: This file is best displayed with 120 col.
+# 
+# Processing of dataset of MS patients. Data organized with BIDS. This script is designed to be run 
+# across multiple subjects in parallel using 'sct_run_batch', but it can also be used to run processing on a single 
+# subject. The input data is assumed to be in BIDS format.
+# 
+# IMPORTANT: This script MUST be run from the root folder of the repository, because it relies on Python scripts located 
+#  in the root folder.
+#
+# Usage:
+#   ./batch_processing.sh <SUBJECT>
+#
+# Example:
+#   ./batch_processing.sh sub-003
+#
+# Author: Julien Cohen-Adad and Pierre-Louis Benveniste
+
+# Parameters
+# =====================================================================================================================
+# 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"
+
+# Uncomment for full verbose
+# set -v
+
+# Immediately exit if error
+set -e
+
+# Exit if user presses CTRL+C (Linux) or CMD+C (OSX)
+trap "echo Caught Keyboard Interrupt within script. Exiting now.; exit" INT
+
+# Save script path
+PATH_SCRIPT=$PWD
+
+
+# CONVENIENCE FUNCTIONS
+# =====================================================================================================================
+
+# label_vertebrae_if_does_not_exist() {
+#   # This function checks if a disc manual label file already exists, then:
+#   #   - If it does, copy it locally.
+#   #   - If it doesn't, perform automatic labeling.
+#   # This allows you to add manual labels on a subject-by-subject basis without disrupting the pipeline.
+
+#   local file="${1}"
+#   local file_seg="${2}"
+#   # Update global variable with segmentation file name
+#   # TODO: modify the name below to _label_disc when Nathan's TotalSpineSeg is ready (the current method outputs a file called FILE_seg_labeled_discs.nii.gz)
+#   FILELABEL="${file}"_seg_labeled_discs
+#   FILELABELMANUAL="${PATH_DATA}"/derivatives/labels/"${SUBJECT}"/anat/"${FILELABEL}".nii.gz
+#   echo "Looking for manual label: ${FILELABELMANUAL}"
+#   if [[ -e "${FILELABELMANUAL}" ]]; then
+#     echo "Found! Copying manual labels."
+#     rsync -avzh "${FILELABELMANUAL}" "${FILELABEL}".nii.gz
+#   else
+#     echo "Not found. Proceeding with automatic labeling."
+#     # Generate labeled segmentation
+#     # TODO: replace with Nathan's TotalSpineSeg
+#     sct_label_vertebrae -i "${file}".nii.gz -s "${file_seg}".nii.gz -c t2 -qc "${PATH_QC}" -qc-subject "${SUBJECT}"
+#   fi
+#   # Generate labeled segmentation based on disc labels
+#   sct_label_vertebrae -i "${file}".nii.gz -s "${file_seg}".nii.gz -discfile "${FILELABEL}".nii.gz -c t2 -qc "${PATH_QC}" -qc-subject "${SUBJECT}"
+#   FILELABELVERTEBRAE="${file}"_seg_labeled
+# }
+
+segment_sc_if_does_not_exist() {
+  # This function checks if a manual spinal cord segmentation file already exists, then:
+  #   - If it does, copy it locally.
+  #   - If it doesn't, perform automatic spinal cord segmentation.
+  # This allows you to add manual segmentations on a subject-by-subject basis without disrupting the pipeline.
+
+  local file="${1}"
+  # Update global variable with segmentation file name
+  FILESEG="${file}"_seg
+  FILESEGMANUAL="${PATH_DATA}"/derivatives/labels/"${SUBJECT}"/anat/"${FILESEG}".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 spinalcord -i "${file}".nii.gz -qc "${PATH_QC}" -qc-subject "${SUBJECT}"
+  fi
+}
+
+segment_lesion_if_does_not_exist() {
+  # This function checks if a manual MS lesion segmentation file already exists, then:
+  #   - If it does, copy it locally.
+  #   - If it doesn't, perform automatic MS lesion segmentation.
+  # This allows you to add manual segmentations on a subject-by-subject basis without disrupting the pipeline.
+
+  local file="${1}"
+  local file_sc_seg="${2}"
+  # Update global variable with segmentation file name
+  FILELESIONSEG="${file}"_label-lesion_seg
+  FILELESIONSEGMANUAL="${PATH_DATA}"/derivatives/labels/"${SUBJECT}"/anat/"${FILELESIONSEG}".nii.gz
+  echo
+  echo "Looking for manual segmentations: ${FILELESIONSEGMANUAL}"
+  if [[ -e "${FILELESIONSEGMANUAL}" ]]; then
+    echo "Found! Using manual segmentation."
+    rsync -avzh "${FILELESIONSEGMANUAL}" "${FILELESIONSEG}".nii.gz
+    sct_qc -i "${file}".nii.gz -p sct_deepseg_lesion -d "${FILELESIONSEG}".nii.gz -s "${file_sc_seg}".nii.gz -qc "${PATH_QC}" -plane axial -qc-subject "${SUBJECT}"
+  else
+    echo "Not found. Proceeding with automatic segmentation."
+    # Segment MS lesion
+    sct_deepseg lesion_ms -i "${file}".nii.gz  -o "${FILELESIONSEG}".nii.gz -qc "${PATH_QC}" -qc-subject "${SUBJECT}" -qc-seg "${file_sc_seg}".nii.gz
+  fi
+}
+
+compare_mayo_segmentation() {
+  # This function compares the Mayo Clinic segmentations with our segmentations.
+  # It requires the Mayo Clinic segmentations to be present in the data folder.
+
+  local file="${1}"
+  local file_sc_seg="${2}"
+  # Update global variable with segmentation file name
+  FILEMAYOSEG="${file}"_label-criticalLesion_dseg
+  FILEMAYOSEGMANUAL="${PATH_DATA}"/derivatives/labels/"${SUBJECT}"/anat/"${FILEMAYOSEG}".nii.gz
+  rsync -avzh "${FILEMAYOSEGMANUAL}" "${FILEMAYOSEG}".nii.gz
+  # The file contains values 1 for lesion and values 2 for segmentation
+  ## We first change the file to remove the values 2, so that we can compare it with our segmentation
+  sct_maths -i "${FILEMAYOSEG}".nii.gz -uthr 1.1 -o "${FILEMAYOSEG}".nii.gz
+  # Generate QC
+  echo "Generating QC for comparison of Mayo Clinic segmentations with our segmentations."
+  sct_qc -i "${file}".nii.gz -p sct_deepseg_lesion -d "${FILEMAYOSEG}".nii.gz -s "${file_sc_seg}".nii.gz -qc "${PATH_QC}" -plane axial -qc-subject "${SUBJECT}"
+}
+
+# SCRIPT STARTS HERE
+# =====================================================================================================================
+
+# Retrieve input params
+SUBJECT="${1}"
+
+# get starting time:
+start="$(date +%s)"
+
+# 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
+rsync -avzh "${PATH_DATA}"/"${SUBJECT}" .
+
+
+# For each subject, we process the T2w image
+# =====================================================================================================================
+cd "${SUBJECT}"/anat/
+file_t2="${SUBJECT}"_acq-ax_T2w
+echo "👉 Processing: ${file_t2}"
+# Segment spinal cord (only if it does not exist)
+# Note: we output the soft segmentation for better CSA precision
+segment_sc_if_does_not_exist "${file_t2}"
+file_t2_sc_seg="${FILESEG}"
+echo SC seg file is here "${file_t2_sc_seg}".nii.gz
+# Segment MS lesions (only if it does not exist)
+segment_lesion_if_does_not_exist "${file_t2}" "${file_t2_sc_seg}"
+file_t2_lesion_seg="${FILELESIONSEG}"
+# Compare mayo segmentations with our segmentations by adding them to the QC just after our predicted segmentation
+compare_mayo_segmentation "${file_t2}" "${file_t2_sc_seg}"
+
+# # Create labels in the cord at mid-vertebral levels
+# label_vertebrae_if_does_not_exist "${file_t2}" "${file_t2_seg}"
+# file_label_vert="${FILELABELVERTEBRAE}"
+# # Compute average CSA as defined by variable 'vertebral_levels'
+# sct_process_segmentation -i "${file_t2_seg}".nii.gz -vertfile "${file_label_vert}".nii.gz \
+#                          -perslice 1 -o "${PATH_RESULTS}"/"${SUBJECT}"_CSA.csv -append 1 -normalize-PAM50 1
+
+
+# # Go back to parent folder
+# cd ..
+
+
+# Verify presence of output files and write log file if error
+# ======================================================================================================================
+FILES_TO_CHECK=(
+  "anat/${file_t2_sc_seg}".nii.gz
+)
+for file in "${FILES_TO_CHECK[@]}"; do
+  if [ ! -e "${file}" ]; then
+    echo "${SUBJECT}/${file} does not exist" >> "${PATH_LOG}/error.log"
+  fi
+done
+
+# Display useful info for the log
+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 "~~~"

detect_atrophy.py

@@ -0,0 +1,62 @@
+"""
+This script detects atrophy in the spinal cord using the output of sct_process_segmentation with PAM-50 normalized CSA. 
+It considers that there is atrophy in the spinal cord if the CSA, if locally the CSA is less than 70% of the average CSA
+of the above and below vertebral level.
+
+Input: 
+    -csa-file: csv file with the CSA values of the spinal cord of the subject.
+
+Output:
+    -vertebral-levels: a list containing the vertebral levels where atrophy was detected. It is empty if no atrophy was detected.
+
+Example of usage:
+    python detect_atrophy.py -csa-file csa_values.csv
+
+Author: Pierre-Louis Benveniste
+"""
+import argparse
+import pandas as pd
+import matplotlib.pyplot as plt
+
+
+def parse_arguments():
+    parser = argparse.ArgumentParser(description='Detect atrophy in the spinal cord from CSA values')
+    parser.add_argument('-csv-file', type=str, required=True, help='CSV file with CSA values of the spinal cord')
+    return parser.parse_args()
+
+
+def main(): 
+    # Parse arguments
+    args = parse_arguments()
+    csv_file = args.csv_file
+
+    # Read csv file
+    data = pd.read_csv(csv_file)
+    # Keep only some columns
+    data = data[['Slice (I->S)', 'VertLevel', 'MEAN(area)', 'MEAN(diameter_RL)']]
+    # Convert VertLevel to string
+    data['VertLevel'] = data['VertLevel'].astype(str)
+
+    # We first plot the mean CSA and mean diameter across the spinal cord
+    fig, ax1 = plt.subplots()
+    ax2 = ax1.twinx()
+    ax1.plot(data['Slice (I->S)'], data['MEAN(area)'], 'g-')
+    ax2.plot(data['Slice (I->S)'], data['MEAN(diameter_RL)'], 'b-')
+    ax1.set_xlabel('Slice (I->S)')
+    ax1.set_ylabel('CSA (mm^2)', color='g')
+    ax2.set_ylabel('Diameter (mm)', color='b')
+    plt.show()
+
+    # We look at the same plot but with the vertebral levels
+    fig, ax1 = plt.subplots()
+    ax2 = ax1.twinx()
+    ax1.plot(data['VertLevel'], data['MEAN(area)'], 'g-')
+    ax2.plot(data['VertLevel'], data['MEAN(diameter_RL)'], 'b-')
+    ax1.set_xlabel('VertLevel')
+    ax1.set_ylabel('CSA (mm^2)', color='g')
+    ax2.set_ylabel('Diameter (mm)', color='b')
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()