Skip to content

Commit 92f635e

Browse files
valosekjvalosekj
authored
Register PSIR/STIR to PAM50 and generate lesion frequency maps (#67)
* add `01_register_to_pam50.sh` * add `02_generate_lesion_frequency_maps.py` * add `03_lfm_to_png.py`
1 parent e58f977 commit 92f635e

File tree

6 files changed

+879
-6
lines changed

6 files changed

+879
-6
lines changed

‎lesion-mapping/01_register_to_pam50.sh

+202
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,202 @@
1+
#!/bin/bash
2+
#
3+
# Register PSIR/STIR contrasts to PAM50 space using the spinal cord segmentation and disc labels
4+
# Then, use the transformation to bring the GT lesion to PAM50 space
5+
#
6+
# NOTE: the commands below assume that GT SC and lesion segmentation and manual disc labels are located under
7+
# derivatives/labels
8+
#
9+
# Usage:
10+
# sct_run_batch -config config.json
11+
#
12+
#
13+
# Example of config.json:
14+
# {
15+
# "path_data" : "<PATH_TO_DATASET>/canproco",
16+
# "path_output" : "<PATH_TO_DATASET>/canproco_register_to_PAM50_2023-10-21",
17+
# "script" : "<PATH_TO_REPO>/canproco/segment_sc_contrast-agnostic/01_register_to_pam50.sh",
18+
# "jobs" : 16,
19+
# "exclude_list" : "sub-mon118 sub-mon006 ..."
20+
# }
21+
#
22+
# The following global variables are retrieved from the caller sct_run_batch
23+
# but could be overwritten by uncommenting the lines below:
24+
# PATH_DATA_PROCESSED="~/data_processed"
25+
# PATH_RESULTS="~/results"
26+
# PATH_LOG="~/log"
27+
# PATH_QC="~/qc"
28+
#
29+
# Author: Jan Valosek
30+
#
31+
32+
# Uncomment for full verbose
33+
set -x
34+
35+
# Immediately exit if error
36+
set -e -o pipefail
37+
38+
# Exit if user presses CTRL+C (Linux) or CMD+C (OSX)
39+
trap "echo Caught Keyboard Interrupt within script. Exiting now.; exit" INT
40+
41+
# Print retrieved variables from the sct_run_batch script to the log (to allow easier debug)
42+
echo "Retrieved variables from from the caller sct_run_batch:"
43+
echo "PATH_DATA: ${PATH_DATA}"
44+
echo "PATH_DATA_PROCESSED: ${PATH_DATA_PROCESSED}"
45+
echo "PATH_RESULTS: ${PATH_RESULTS}"
46+
echo "PATH_LOG: ${PATH_LOG}"
47+
echo "PATH_QC: ${PATH_QC}"
48+
49+
SUBJECT=$1
50+
PATH_SCRIPT=$2
51+
PATH_MODEL=$3
52+
53+
echo "SUBJECT: ${SUBJECT}"
54+
echo "PATH_SCRIPT: ${PATH_SCRIPT}"
55+
echo "PATH_MODEL: ${PATH_MODEL}"
56+
57+
# ------------------------------------------------------------------------------
58+
# CONVENIENCE FUNCTIONS
59+
# ------------------------------------------------------------------------------
60+
61+
# Copy GT lesion segmentation
62+
copy_gt(){
63+
local file="$1"
64+
local type="$2" # seg or lesion
65+
# Construct file name to GT lesion segmentation located under derivatives/labels
66+
FILESEGMANUAL="${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${file}_${type}-manual.nii.gz"
67+
echo ""
68+
echo "Looking for manual segmentation: $FILESEGMANUAL"
69+
if [[ -e $FILESEGMANUAL ]]; then
70+
echo "Found! Copying ..."
71+
rsync -avzh $FILESEGMANUAL ${file}_${type}-manual.nii.gz
72+
else
73+
echo "File ${FILESEGMANUAL}.nii.gz does not exist" >> ${PATH_LOG}/missing_files.log
74+
echo "ERROR: Manual GT segmentation ${FILESEGMANUAL}.nii.gz does not exist. Exiting."
75+
exit 1
76+
fi
77+
}
78+
79+
# Copy manually created disc labels from derivatives/labels
80+
copy_disc_labels(){
81+
local file="$1"
82+
# Update global variable with disc labels file name
83+
FILELABEL="${file}_labels-disc"
84+
FILELABELMANUAL="${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${FILELABEL}.nii.gz"
85+
echo "Looking for manual disc labels: $FILELABELMANUAL"
86+
if [[ -e $FILELABELMANUAL ]]; then
87+
echo "Found! Copying manual disc labels."
88+
rsync -avzh $FILELABELMANUAL ${FILELABEL}.nii.gz
89+
else
90+
echo "File ${FILESEGMANUAL}.nii.gz does not exist" >> ${PATH_LOG}/missing_files.log
91+
echo "ERROR: Manual disc labels ${FILESEGMANUAL}.nii.gz does not exist. Exiting."
92+
exit 1
93+
fi
94+
}
95+
96+
# ------------------------------------------------------------------------------
97+
# SCRIPT STARTS HERE
98+
# ------------------------------------------------------------------------------
99+
# get starting time:
100+
start=`date +%s`
101+
102+
# Display useful info for the log, such as SCT version, RAM and CPU cores available
103+
sct_check_dependencies -short
104+
105+
# Go to folder where data will be copied and processed
106+
cd $PATH_DATA_PROCESSED
107+
108+
# Copy source PSIR/STIR images
109+
# Note: we use '/./' in order to include the sub-folder 'ses-M0'
110+
# We do a substitution '/' --> '_' in case there is a subfolder 'ses-M0/'
111+
rsync -Ravzh ${PATH_DATA}/./${SUBJECT}/anat/${SUBJECT//[\/]/_}_*IR.* .
112+
113+
# Go to subject folder for source images
114+
cd ${SUBJECT}/anat
115+
116+
# ------------------------------------------------------------------------------
117+
# PSIR/STIR
118+
# ------------------------------------------------------------------------------
119+
120+
# We do a substitution '/' --> '_' in case there is a subfolder 'ses-M0/'
121+
# Calgary has STIR
122+
if [[ ${SUBJECT} =~ "cal" ]]; then
123+
file="${SUBJECT//[\/]/_}"_STIR
124+
# Other sites have PSIR
125+
else
126+
file="${SUBJECT//[\/]/_}"_PSIR
127+
fi
128+
129+
# Check if file exists
130+
if [[ ! -e ${file}.nii.gz ]]; then
131+
echo "File ${file}.nii.gz does not exist" >> ${PATH_LOG}/missing_files.log
132+
echo "ERROR: File ${file}.nii.gz does not exist. Exiting."
133+
exit 1
134+
else
135+
136+
# Create a variable with fname corresponding to the files under derivatives/labels (because we modify the variable
137+
# for PSIR images)
138+
file_gt=${file}
139+
140+
if [[ ${file} =~ "PSIR" ]]; then
141+
# For PSIR, swap contrast from light cord and dark CSF to dark cord and light CSF to use -c t2 during template registration
142+
# Context: https://github.com/ivadomed/canproco/issues/46#issuecomment-1752142304
143+
# TODO: this line will be deleted once the SCT script will include the flag for contrast swapping
144+
sct_maths -i ${file}.nii.gz -mul -1 -o ${file}_mul.nii.gz
145+
file=${file}_mul
146+
fi
147+
148+
# NOTE: the commands below assume that GT SC, GT lesion segmentation and manual disc labels are located under
149+
# derivatives/labels
150+
# Copy GT SC seg
151+
copy_gt "${file_gt}" "seg"
152+
file_seg="${file_gt}_seg-manual"
153+
154+
# Copy GT lesion seg
155+
copy_gt "${file_gt}" "lesion"
156+
file_lesion="${file_gt}_lesion-manual"
157+
158+
# Copy manually created disc labels from derivatives/labels
159+
copy_disc_labels "${file_gt}"
160+
file_disc="${file_gt}_labels-disc"
161+
162+
# Binarize GT lesion segmentation
163+
sct_maths -i ${file_lesion}.nii.gz -bin 0 -o ${file_lesion}_bin.nii.gz
164+
165+
# Register to template
166+
# The flag '-param' is inspired by: Eden, D. et al. Brain, 2019. https://doi.org/10.1093/brain/awy352
167+
# https://github.com/neuropoly/lesion-mapping/blob/497584281526bbcee2f56bfa9ec076982c7f1503/spinalcord/1_register_data.py#L22
168+
# Note: we use '-ldisc' to use all disc labels
169+
sct_register_to_template -i ${file}.nii.gz -s ${file_seg}.nii.gz -ldisc ${file_disc}.nii.gz -c t2 -param step=1,type=seg,algo=centermass,metric=MeanSquares,slicewise=1:step=2,type=seg,algo=bsplinesyn,metric=MeanSquares,slicewise=1,iter=3 -qc ${PATH_QC} -qc-subject ${SUBJECT}
170+
171+
# Generate additional QC to check the registration
172+
# Native image in PAM50 space
173+
# https://github.com/spinalcordtoolbox/spinalcordtoolbox/issues/4166#issuecomment-1773808040
174+
sct_qc -i ${SCT_DIR}/data/PAM50/template/PAM50_t2.nii.gz -s ${SCT_DIR}/data/PAM50/template/PAM50_cord.nii.gz -d anat2template.nii.gz -p sct_image_stitch -qc ${PATH_QC} -qc-subject ${SUBJECT}
175+
# Native image in PAM50 space overlaid with PAM50_levels
176+
# https://github.com/spinalcordtoolbox/spinalcordtoolbox/issues/4166#issuecomment-1773810021
177+
sct_qc -i anat2template.nii.gz -s ${SCT_DIR}/data/PAM50/template/PAM50_levels.nii.gz -p sct_label_vertebrae -qc ${PATH_QC} -qc-subject ${SUBJECT}
178+
179+
# Bring GT lesion segmentation to template space
180+
sct_apply_transfo -i ${file_lesion}_bin.nii.gz -d ${SCT_DIR}/data/PAM50/template/PAM50_t2.nii.gz -w warp_anat2template.nii.gz -x linear -o ${file_lesion}_bin_reg.nii.gz
181+
# Generate QC (native image in PAM50 space overlaid with GT lesion segmentation in PAM50 space)
182+
sct_qc -i anat2template.nii.gz -d ${file_lesion}_bin_reg.nii.gz -s ${SCT_DIR}/data/PAM50/template/PAM50_cord.nii.gz -p sct_deepseg_lesion -plane sagittal -qc ${PATH_QC} -qc-subject ${SUBJECT}
183+
184+
# Bring GT SC segmentation to template space
185+
sct_apply_transfo -i ${file_seg}.nii.gz -d ${SCT_DIR}/data/PAM50/template/PAM50_t2.nii.gz -w warp_anat2template.nii.gz -x linear -o ${file_seg}_reg.nii.gz
186+
# Generate QC (native image in PAM50 space overlaid with GT SC segmentation in PAM50 space)
187+
sct_qc -i anat2template.nii.gz -d ${file_seg}_reg.nii.gz -s ${SCT_DIR}/data/PAM50/template/PAM50_cord.nii.gz -p sct_deepseg_lesion -plane sagittal
188+
fi
189+
190+
# ------------------------------------------------------------------------------
191+
# End
192+
# ------------------------------------------------------------------------------
193+
194+
# Display results (to easily compare integrity across SCT versions)
195+
end=`date +%s`
196+
runtime=$((end-start))
197+
echo
198+
echo "~~~"
199+
echo "SCT version: `sct_version`"
200+
echo "Ran on: `uname -nsr`"
201+
echo "Duration: $(($runtime / 3600))hrs $((($runtime / 60) % 60))min $(($runtime % 60))sec"
202+
echo "~~~"

0 commit comments

Comments
 (0)