From 1f0fd3d262c4fc48fe954d4ea19fd4e059eecd0f Mon Sep 17 00:00:00 2001 From: valosekj Date: Fri, 18 Jul 2025 14:48:52 +0200 Subject: [PATCH 01/17] Add scripts to compare AP and RL diameters computed using 'skimage.regionprops' and 'HOG' SCT PR: https://github.com/spinalcordtoolbox/spinalcordtoolbox/pull/4958 --- figures/AP_RL_diameters_PAM50.py | 192 ++++++++++++++++++++++++++++ figures/AP_RL_diameters_plotting.py | 116 +++++++++++++++++ 2 files changed, 308 insertions(+) create mode 100644 figures/AP_RL_diameters_PAM50.py create mode 100644 figures/AP_RL_diameters_plotting.py diff --git a/figures/AP_RL_diameters_PAM50.py b/figures/AP_RL_diameters_PAM50.py new file mode 100644 index 0000000..009443b --- /dev/null +++ b/figures/AP_RL_diameters_PAM50.py @@ -0,0 +1,192 @@ +# +# Plot a single subject morphometric metrics in the PAM50 space per slice and vertebral levels +# +# You can use SCT's conda environment to run this script: +# # Go to the SCT directory +# cd $SCT_DIR +# # Activate SCT conda environment +# source ./python/etc/profile.d/conda.sh +# conda activate venv_sct +# +# Example usage on a single subject: +# python AP_RL_diameters_PAM50.py -i /path/to/morphometrics_PAM50.csv -o /path/to/output/morphometrics_PAM50.png +# +# Authors: Jan Valosek +# + +import sys +import argparse +import numpy as np +import pandas as pd +import seaborn as sns +import matplotlib as mpl +import matplotlib.pyplot as plt + +METRICS = [ + 'MEAN(diameter_AP)', + 'MEAN(diameter_RL)', +] + +LABELS_FONT_SIZE = 14 +TICKS_FONT_SIZE = 12 + + +def get_parser(): + parser = argparse.ArgumentParser( + description="Plot AP and RL diameters measured using different methods in the PAM50 space.") + parser.add_argument('-i', required=True, type=str, + help="Path to the CSV file with diameter measurements.") + parser.add_argument('-o', required=True, type=str, + help="Path to save the output figure.") + + return parser + + +# Apply smoothing to the metric data +def smooth(y: np.ndarray, box_pts: int) -> np.ndarray: + """ + Smooths a 1D array using a simple moving average (box filter). + Inspired by: https://github.com/sct-pipeline/rootlets-informed-reg2template/blob/main/csa_analysis.py#L199 + Args: + y (np.ndarray): Input 1D array to be smoothed. + box_pts (int): Number of points in the moving average window. Needs to be >= 1. + Returns: + np.ndarray: Smoothed array of the same length as the input. + """ + box = np.ones(box_pts) / box_pts + y_smooth = np.convolve(y, box, mode='same') + return y_smooth + + +def load_single_subject_data(path_single_subject, df_spine_generic_min, df_spine_generic_max): + """ + Load single subject data + :param path_single_subject: path to single subject CSV file (from session1 or session2) + :param df_spine_generic_min: minimum slice number from spine-generic dataset + :param df_spine_generic_max: maximum slice number from spine-generic dataset + :return: + """ + df_single_subject = pd.read_csv(path_single_subject) + # Compute compression ratio (CR) as MEAN(diameter_AP) / MEAN(diameter_RL) + df_single_subject['MEAN(compression_ratio)'] = df_single_subject['MEAN(diameter_AP)'] / \ + df_single_subject['MEAN(diameter_RL)'] + + # Keep only slices from C1 to Th1 to match the slices of the spine-generic normative values + df_single_subject = df_single_subject[(df_single_subject['Slice (I->S)'] >= df_spine_generic_min) & + (df_single_subject['Slice (I->S)'] <= df_spine_generic_max)] + + return df_single_subject + +def get_vert_indices(df): + """ + Get indices of slices corresponding to mid-vertebrae + Args: + df (pd.dataFrame): dataframe with CSA values + Returns: + vert (pd.Series): vertebrae levels across slices + ind_vert (np.array): indices of slices corresponding to the beginning of each level (=intervertebral disc) + ind_vert_mid (np.array): indices of slices corresponding to mid-levels + """ + # Get unique participant IDs + subjects = df['Filename'].unique() + # Get vert levels for one certain subject + vert = df[df['Filename'] == subjects[0]]['VertLevel'] + # Get indexes of where array changes value + ind_vert = vert.diff()[vert.diff() != 0].index.values + # Get the beginning of C1 + ind_vert = np.append(ind_vert, vert.index.values[-1]) + ind_vert_mid = [] + # Get indexes of mid-vertebrae + for i in range(len(ind_vert)-1): + ind_vert_mid.append(int(ind_vert[i:i+2].mean())) + + return vert, ind_vert, ind_vert_mid + + +def create_lineplot(df, figure_path): + """ + Create lineplot for individual metrics per vertebral levels. + Args: + df (pd.DataFrame): dataframe with single subject values + figure_path (str): path to save the figure + """ + mpl.rcParams['font.family'] = 'Arial' + + # Plot figures + fig, axes = plt.subplots(1, 2, figsize=(10, 5)) + axs = axes.ravel() + + # Loop across metrics + for index, metric in enumerate(METRICS): + + # Smooth the data to improve visualization + df[metric] = smooth(df[metric].values, 5) + df[f'{metric.replace(")", "_hog)")}'] = smooth(df[f'{metric.replace(")", "_hog)")}'].values, 5) + sns.lineplot(ax=axs[index], x="Slice (I->S)", y=metric, + data=df, linewidth=2, + label=metric) + sns.lineplot(ax=axs[index], x="Slice (I->S)", y=f'{metric.replace(")", "_hog)")}', + data=df, linewidth=2, + label=f'{metric.replace(")", "_hog)")}') + + # Tweak y-axis limits + # ymin, ymax = METRICS_YLIMITS[metric] + # axes.set_ylim(ymin, ymax) + # Remove first and last 4 slices from the x-axis to remove smoothing artifacts + axs[index].set_xlim(df['Slice (I->S)'].iloc[4], df['Slice (I->S)'].iloc[-4]) + + # Remove xticks to hide PAM50 Axial Slice numbers + axs[index].set_xticks([]) + axs[index].tick_params(axis='both', which='major', labelsize=TICKS_FONT_SIZE) + axs[index].spines['right'].set_visible(False) + axs[index].spines['left'].set_visible(False) + axs[index].spines['top'].set_visible(False) + axs[index].spines['bottom'].set_visible(True) + + # Get ymin and ymax for the y-axis + ymin, ymax = axs[index].get_ylim() + + vert, ind_vert, ind_vert_mid = get_vert_indices(df) + for idx, x in enumerate(ind_vert[1:-1]): + axs[index].axvline(df.loc[x, 'Slice (I->S)'], color='black', linestyle='--', alpha=0.3, zorder=0) + for idx, x in enumerate(ind_vert_mid, 0): + level = f'T{vert[x] - 7}' if vert[x] > 7 else f'C{vert[x]}' + axs[index].text(df.loc[ind_vert_mid[idx], 'Slice (I->S)'], + ymin - (ymax-ymin)*0.05, level, horizontalalignment='center', + verticalalignment='bottom', color='black', fontsize=TICKS_FONT_SIZE) + + axs[index].invert_xaxis() + axs[index].yaxis.grid(True) + axs[index].set_axisbelow(True) + axs[index].set_ylabel(f'{metric}', fontsize=LABELS_FONT_SIZE) + axs[index].set_xlabel('') + + # Decrease the font size of the legend + axs[index].legend(loc='lower left', fontsize=TICKS_FONT_SIZE) + + plt.tight_layout() + plt.savefig(figure_path, dpi=300, bbox_inches='tight') + print(f'Figure saved: {figure_path}') + plt.show() + # plt.close() + + +def main(): + parser = get_parser() + args = parser.parse_args() + + # Get the file path from args + csv_path = args.i + out_path = args.o + + # Create a list of dataframes for each session file + df = load_single_subject_data(csv_path, 700, 980) + if df.empty: + print('WARNING: No slices found in the range C1-Th1 in the single subject data. Exiting...') + sys.exit(1) + + create_lineplot(df, out_path) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/figures/AP_RL_diameters_plotting.py b/figures/AP_RL_diameters_plotting.py new file mode 100644 index 0000000..d6c0a56 --- /dev/null +++ b/figures/AP_RL_diameters_plotting.py @@ -0,0 +1,116 @@ +# +# Script for comparing AP and RL diameters measured using different methods (skimage.regionprops vs. HOG). +# Context: https://github.com/spinalcordtoolbox/spinalcordtoolbox/pull/4958/ +# This script generates scatter plots with regression lines to visualize the correlation +# between different measurement techniques. +# +# You can use SCT's conda environment to run this script: +# # Go to the SCT directory +# cd $SCT_DIR +# # Activate SCT conda environment +# source ./python/etc/profile.d/conda.sh +# conda activate venv_sct +# +# Example usage on a single subject: +# python AP_RL_diameters_plotting.py -i /path/to/morphometrics_plotting.csv -o /path/to/output/morphometrics.png +# +# Authors: Jan Valosek +# + +import argparse +import pandas as pd +import matplotlib as mpl +import matplotlib.pyplot as plt +from scipy.stats import linregress + + +def get_parser(): + """ + Parse command line arguments. + """ + parser = argparse.ArgumentParser( + description="Create scatter plots with regression lines for AP and RL diameters measured using different methods.") + parser.add_argument('-i', required=True, type=str, + help="Path to the CSV file with diameter measurements.") + parser.add_argument('-o', required=True, type=str, + help="Path to save the output figure.") + return parser + + +def plot_with_regression(x: pd.Series, y: pd.Series, xlabel: str, ylabel: str, ax: plt.Axes) -> None: + """Plot scatter with regression line""" + # Drop NaNs + df = pd.DataFrame({xlabel: x, ylabel: y}).dropna() + x_clean, y_clean = df[xlabel], df[ylabel] + + # Regression: + # slope -- Slope of the regression line + # intercept -- Intercept of the regression line + # r -- The Pearson correlation coefficient + slope, intercept, r, _, _ = linregress(x_clean, y_clean) + reg_line = slope * x_clean + intercept + + mpl.rcParams['font.family'] = 'Arial' + + # Scatter + ax.scatter(x_clean, y_clean, s=20, edgecolor='k', facecolor='gray', alpha=0.7) + ax.plot(x_clean, reg_line, 'r-', linewidth=2) + + # Identity line + min_val, max_val = min(x_clean.min(), y_clean.min()), max(x_clean.max(), y_clean.max()) + ax.plot([min_val, max_val], [min_val, max_val], 'k--', linewidth=1) + + # Labels and text + ax.set_xlabel(xlabel, fontsize=10) + ax.set_ylabel(ylabel, fontsize=10) + ax.set_aspect('equal', adjustable='box') + ax.set_xlim(min_val*0.95, max_val*1.05) + ax.set_ylim(min_val*0.95, max_val*1.05) + ax.text(0.05, 0.95, + f'$r$ = {r:.2f}\n$y$ = {slope:.2f}$x$ + {intercept:.2f}', + transform=ax.transAxes, + verticalalignment='top', + fontsize=9) + + # Style + ax.spines[['top', 'right']].set_visible(False) + ax.tick_params(labelsize=8) + + +def main() -> None: + """Read CSV and generate comparison plots for AP and RL diameters.""" + parser = get_parser() + args = parser.parse_args() + + # Use command line arguments + csv_path = args.i + out_path = args.o + + df = pd.read_csv(csv_path) + + fig, axes = plt.subplots(1, 2, figsize=(10, 5), dpi=300) + plot_with_regression( + df['MEAN(diameter_AP)'], + df['MEAN(diameter_AP_hog)'], + 'AP diameter (skimage.regionprops)', + 'AP diameter (HOG)', + axes[0] + ) + plot_with_regression( + df['MEAN(diameter_RL)'], + df['MEAN(diameter_RL_hog)'], + 'RL diameter (skimage.regionprops)', + 'RL diameter (HOG)', + axes[1] + ) + plt.tight_layout() + + # Save figure + plt.savefig(out_path, bbox_inches='tight', dpi=300) + print(f"Figure saved: {out_path}") + plt.show() + # plt.close() + + +if __name__ == "__main__": + main() From 4eef5adc75218c6cc8548e457a9ef7c2e3bb58f0 Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 10:17:34 +0200 Subject: [PATCH 02/17] improve description and comments --- figures/AP_RL_diameters_PAM50.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/figures/AP_RL_diameters_PAM50.py b/figures/AP_RL_diameters_PAM50.py index 009443b..48a0e1f 100644 --- a/figures/AP_RL_diameters_PAM50.py +++ b/figures/AP_RL_diameters_PAM50.py @@ -1,5 +1,6 @@ # # Plot a single subject morphometric metrics in the PAM50 space per slice and vertebral levels +# Original vs HOG-based AP and RL diameters # # You can use SCT's conda environment to run this script: # # Go to the SCT directory @@ -122,9 +123,11 @@ def create_lineplot(df, figure_path): # Smooth the data to improve visualization df[metric] = smooth(df[metric].values, 5) df[f'{metric.replace(")", "_hog)")}'] = smooth(df[f'{metric.replace(")", "_hog)")}'].values, 5) + # Original sns.lineplot(ax=axs[index], x="Slice (I->S)", y=metric, data=df, linewidth=2, label=metric) + # HOG-based sns.lineplot(ax=axs[index], x="Slice (I->S)", y=f'{metric.replace(")", "_hog)")}', data=df, linewidth=2, label=f'{metric.replace(")", "_hog)")}') From 8e6b3074c292a4138c8047edb3007f597af68114 Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 10:17:56 +0200 Subject: [PATCH 03/17] comment out plt.show() and uncomment plt.close() --- figures/AP_RL_diameters_PAM50.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/figures/AP_RL_diameters_PAM50.py b/figures/AP_RL_diameters_PAM50.py index 48a0e1f..5d255f3 100644 --- a/figures/AP_RL_diameters_PAM50.py +++ b/figures/AP_RL_diameters_PAM50.py @@ -170,8 +170,8 @@ def create_lineplot(df, figure_path): plt.tight_layout() plt.savefig(figure_path, dpi=300, bbox_inches='tight') print(f'Figure saved: {figure_path}') - plt.show() - # plt.close() + #plt.show() + plt.close() def main(): From 74062e365fa564effa9225fc17370dddec780535 Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 10:18:13 +0200 Subject: [PATCH 04/17] add compute_diameters.sh to compute AP and RL diameters using "traditional" and HOG-based methods --- figures/compute_diameters.sh | 159 +++++++++++++++++++++++++++++++++++ 1 file changed, 159 insertions(+) create mode 100644 figures/compute_diameters.sh diff --git a/figures/compute_diameters.sh b/figures/compute_diameters.sh new file mode 100644 index 0000000..f9e62b4 --- /dev/null +++ b/figures/compute_diameters.sh @@ -0,0 +1,159 @@ +#!/bin/bash +# +# Compute AP and RL diameters using "traditional" and HOG-based methods +# +# Requirements: +# SCT PR #4958: https://github.com/spinalcordtoolbox/spinalcordtoolbox/pull/4958 +# +# Usage: +# sct_run_batch -config config.json +# +# Example of config.json: +# { +# "script" : "~/code/dcm-metric-normalization/figures/compute_diameters.sh", +# "jobs" : 8 +# } +# +# Example of dataset.txt: +# path/to/dataset1 path/to/output1 +# path/to/dataset2 path/to/output1 +# +# Manual segmentations and disc labels should be located under: +# PATH_DATA/derivatives/labels/SUBJECT// +# +# Author: 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}" + +SUBJECT=$1 + +echo "SUBJECT: ${SUBJECT}" + +# 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 --no-g $PATH_DATA/./$SUBJECT . +# • -R → Preserves relative path structure. +# • -a → Archive mode (preserves most attributes except group ownership). +# • -v → Verbose output. +# • -z → Compresses data during transfer. +# • -h → Human-readable file sizes. +# • --no-g → Skips changing group ownership. + +# Remove dwi folder if it exists (because we do not process DWI data in this script) +if [[ -d ${SUBJECT}/dwi ]]; then + echo "Removing DWI folder: ${SUBJECT}/dwi" + rm -rf ${SUBJECT}/dwi +fi + +# 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//[\/]/_}" +echo "file: ${file}" + +# ------------------------------------------------------------------------- +# T2w +# ------------------------------------------------------------------------- +# Steps: +# - find T2w image +# - remove all other files except the selected T2w and its derivatives to save space +# - copy SC segmentation from derivatives/labels folder (we assume it exists) +# - copy disc labels from derivatives/labels folder (we assume it exists) +# - generate labeled segmentation using init disc labels +# - generate metrics using sct_process_segmentation +# - generate figures using custom script + +# Find T2w image +if [[ -e "${file}_T2w.nii.gz" ]]; then + file_t2="${file}_T2w" +else + echo "No suitable T2w file found for subject ${SUBJECT}" >&2 + exit 1 +fi + +# Remove all other files except the selected T2w and its derivatives to save space +find . -maxdepth 1 -type f \ + ! -name "${file_t2}.*" \ + ! -name "${file_t2}_*" \ + -exec rm -f {} + +echo "file_t2: ${file_t2}" + +echo "👉 Processing: ${file_t2}" + +# Copy SC segmentation from derivatives/labels folder (we assume it exists) +FILESEG="${file_t2}_label-SC_seg" +rsync -avzh --no-g "${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${FILESEG}.nii.gz" "${FILESEG}.nii.gz" +echo "✅ [$(date '+%Y-%m-%d %H:%M:%S')] ${FILESEG}.nii.gz found under derivatives/labels --> copying it" >> "${PATH_LOG}/T2w_sc_segmentation.log" + +# Copy disc labels from derivatives/labels folder (we assume it exists) as we need them for vertebral labeling +FILEDISCS="${file_t2}_labels-disc-manual" +# whole-spine +if [[ -e "${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${file_t2}_labels-disc-manual.nii.gz" ]]; then + rsync -avzh --no-g "${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${file_t2}_labels-disc-manual.nii.gz" "${FILEDISCS}.nii.gz" + echo "✅ [$(date '+%Y-%m-%d %H:%M:%S')] ${file_t2}_labels-disc-manual.nii.gz found under derivatives/labels --> copying it" >> "${PATH_LOG}/T2w_disc_labels.log" +# spine-generic +elif [[ -e "${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${file_t2}_label-discs_dlabel.nii.gz" ]]; then + rsync -avzh --no-g "${PATH_DATA}/derivatives/labels/${SUBJECT}/anat/${file_t2}_label-discs_dlabel.nii.gz" "${FILEDISCS}.nii.gz" + echo "✅ [$(date '+%Y-%m-%d %H:%M:%S')] ${file_t2}_label-discs_dlabel.nii.gz found under derivatives/labels --> copying it" >> "${PATH_LOG}/T2w_disc_labels.log" +fi + +# Generate labeled segmentation using init disc labels +sct_label_vertebrae -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -discfile ${FILEDISCS}.nii.gz -c t2 -qc ${PATH_QC} -qc-subject ${file} + +sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -o ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr0.csv -qc qc -angle-corr 0 +sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -o ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr1.csv -qc qc -angle-corr 1 + +# Create figures folder under ${PATH_RESULTS} if it does not exist +mkdir -p ${PATH_RESULTS}/figures + +# Generate figure using custom script (assuming this repo is cloned in `~/code/dcm-metric-normalization`) +${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/figures/AP_RL_diameters_PAM50.py \ + -i ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr0.csv \ + -o ${PATH_RESULTS}/figures/${file}_AP_RL_diameters_PAM50_angle_corr0.png \ + +${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/figures/AP_RL_diameters_PAM50.py \ + -i ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr1.csv \ + -o ${PATH_RESULTS}/figures/${file}_AP_RL_diameters_PAM50_angle_corr1.png + +# ------------------------------------------------------------------------------ +# End +# ------------------------------------------------------------------------------ +# 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 "~~~" From 671f19fb00740a342795f970d2cfdb919233aa4d Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 10:22:09 +0200 Subject: [PATCH 05/17] rename `figures` folder to `hog_angle` --- {figures => hog_angle}/AP_RL_diameters_PAM50.py | 0 {figures => hog_angle}/AP_RL_diameters_plotting.py | 0 {figures => hog_angle}/compute_diameters.sh | 8 ++++---- 3 files changed, 4 insertions(+), 4 deletions(-) rename {figures => hog_angle}/AP_RL_diameters_PAM50.py (100%) rename {figures => hog_angle}/AP_RL_diameters_plotting.py (100%) rename {figures => hog_angle}/compute_diameters.sh (96%) diff --git a/figures/AP_RL_diameters_PAM50.py b/hog_angle/AP_RL_diameters_PAM50.py similarity index 100% rename from figures/AP_RL_diameters_PAM50.py rename to hog_angle/AP_RL_diameters_PAM50.py diff --git a/figures/AP_RL_diameters_plotting.py b/hog_angle/AP_RL_diameters_plotting.py similarity index 100% rename from figures/AP_RL_diameters_plotting.py rename to hog_angle/AP_RL_diameters_plotting.py diff --git a/figures/compute_diameters.sh b/hog_angle/compute_diameters.sh similarity index 96% rename from figures/compute_diameters.sh rename to hog_angle/compute_diameters.sh index f9e62b4..76bbab0 100644 --- a/figures/compute_diameters.sh +++ b/hog_angle/compute_diameters.sh @@ -10,7 +10,7 @@ # # Example of config.json: # { -# "script" : "~/code/dcm-metric-normalization/figures/compute_diameters.sh", +# "script" : "~/code/dcm-metric-normalization/hog_angle/compute_diameters.sh", # "jobs" : 8 # } # @@ -91,7 +91,7 @@ echo "file: ${file}" # - copy disc labels from derivatives/labels folder (we assume it exists) # - generate labeled segmentation using init disc labels # - generate metrics using sct_process_segmentation -# - generate figures using custom script +# - generate hog_angle using custom script # Find T2w image if [[ -e "${file}_T2w.nii.gz" ]]; then @@ -137,11 +137,11 @@ sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${F mkdir -p ${PATH_RESULTS}/figures # Generate figure using custom script (assuming this repo is cloned in `~/code/dcm-metric-normalization`) -${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/figures/AP_RL_diameters_PAM50.py \ +${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/hog_angle/AP_RL_diameters_PAM50.py \ -i ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr0.csv \ -o ${PATH_RESULTS}/figures/${file}_AP_RL_diameters_PAM50_angle_corr0.png \ -${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/figures/AP_RL_diameters_PAM50.py \ +${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/hog_angle/AP_RL_diameters_PAM50.py \ -i ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr1.csv \ -o ${PATH_RESULTS}/figures/${file}_AP_RL_diameters_PAM50_angle_corr1.png From 77c8cfe68f925635b37b08ea4b003a9c400c93a1 Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 15:19:02 +0200 Subject: [PATCH 06/17] rename `AP_RL_diameters_plotting.py` folder to `AP_RL_diameters.py` --- hog_angle/{AP_RL_diameters_plotting.py => AP_RL_diameters.py} | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename hog_angle/{AP_RL_diameters_plotting.py => AP_RL_diameters.py} (96%) diff --git a/hog_angle/AP_RL_diameters_plotting.py b/hog_angle/AP_RL_diameters.py similarity index 96% rename from hog_angle/AP_RL_diameters_plotting.py rename to hog_angle/AP_RL_diameters.py index d6c0a56..e0a5627 100644 --- a/hog_angle/AP_RL_diameters_plotting.py +++ b/hog_angle/AP_RL_diameters.py @@ -12,7 +12,7 @@ # conda activate venv_sct # # Example usage on a single subject: -# python AP_RL_diameters_plotting.py -i /path/to/morphometrics_plotting.csv -o /path/to/output/morphometrics.png +# python AP_RL_diameters.py -i /path/to/morphometrics_plotting.csv -o /path/to/output/morphometrics.png # # Authors: Jan Valosek # From d01c4ede78731657ab45dc86c920c12f7867394a Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 15:19:56 +0200 Subject: [PATCH 07/17] Allow to specify smoothing as an input arg --- hog_angle/AP_RL_diameters_PAM50.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/hog_angle/AP_RL_diameters_PAM50.py b/hog_angle/AP_RL_diameters_PAM50.py index 5d255f3..1627266 100644 --- a/hog_angle/AP_RL_diameters_PAM50.py +++ b/hog_angle/AP_RL_diameters_PAM50.py @@ -39,12 +39,15 @@ def get_parser(): help="Path to the CSV file with diameter measurements.") parser.add_argument('-o', required=True, type=str, help="Path to save the output figure.") + parser.add_argument('-smooth', required=False, type=int, default=0, metavar='INT', + help="Smooth measurements before plotting. Number of points in the moving average window." + "Examples: 0: no smoothing, 5: window of 5") return parser # Apply smoothing to the metric data -def smooth(y: np.ndarray, box_pts: int) -> np.ndarray: +def smooth_data(y: np.ndarray, box_pts: int) -> np.ndarray: """ Smooths a 1D array using a simple moving average (box filter). Inspired by: https://github.com/sct-pipeline/rootlets-informed-reg2template/blob/main/csa_analysis.py#L199 @@ -104,12 +107,13 @@ def get_vert_indices(df): return vert, ind_vert, ind_vert_mid -def create_lineplot(df, figure_path): +def create_lineplot(df, figure_path, smooth): """ Create lineplot for individual metrics per vertebral levels. Args: df (pd.DataFrame): dataframe with single subject values figure_path (str): path to save the figure + smooth (int): Smooth measurements before plotting. 0: no smoothing; 1: smoothing """ mpl.rcParams['font.family'] = 'Arial' @@ -120,9 +124,10 @@ def create_lineplot(df, figure_path): # Loop across metrics for index, metric in enumerate(METRICS): - # Smooth the data to improve visualization - df[metric] = smooth(df[metric].values, 5) - df[f'{metric.replace(")", "_hog)")}'] = smooth(df[f'{metric.replace(")", "_hog)")}'].values, 5) + if smooth > 0: + # Smooth the data to improve visualization + df[metric] = smooth_data(df[metric].values, smooth) + df[f'{metric.replace(")", "_hog)")}'] = smooth_data(df[f'{metric.replace(")", "_hog)")}'].values, smooth) # Original sns.lineplot(ax=axs[index], x="Slice (I->S)", y=metric, data=df, linewidth=2, @@ -188,7 +193,7 @@ def main(): print('WARNING: No slices found in the range C1-Th1 in the single subject data. Exiting...') sys.exit(1) - create_lineplot(df, out_path) + create_lineplot(df, out_path, args.smooth) if __name__ == '__main__': From 8f2dbd95565b9a32462b739c2aee9cf0dbcc66a5 Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 15:20:29 +0200 Subject: [PATCH 08/17] Use the same y-axis for all subjects --- hog_angle/AP_RL_diameters_PAM50.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/hog_angle/AP_RL_diameters_PAM50.py b/hog_angle/AP_RL_diameters_PAM50.py index 1627266..c442193 100644 --- a/hog_angle/AP_RL_diameters_PAM50.py +++ b/hog_angle/AP_RL_diameters_PAM50.py @@ -28,6 +28,12 @@ 'MEAN(diameter_RL)', ] +# Set ylim to do not overlap horizontal grid with vertebrae labels +METRICS_TO_YLIM = { + 'MEAN(diameter_AP)': (5.5, 9.5), + 'MEAN(diameter_RL)': (8.5, 14.5), +} + LABELS_FONT_SIZE = 14 TICKS_FONT_SIZE = 12 @@ -152,6 +158,7 @@ def create_lineplot(df, figure_path, smooth): axs[index].spines['bottom'].set_visible(True) # Get ymin and ymax for the y-axis + axs[index].set_ylim(METRICS_TO_YLIM[metric][0], METRICS_TO_YLIM[metric][1]) ymin, ymax = axs[index].get_ylim() vert, ind_vert, ind_vert_mid = get_vert_indices(df) From 9d19adefff2a802e6e446b881ada5f6c377d2471 Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 15:22:29 +0200 Subject: [PATCH 09/17] remove `-R` from rsync to prevent copying the full path Context: https://github.com/spinalcordtoolbox/spinalcordtoolbox/issues/4930#issuecomment-3132524965 --- hog_angle/compute_diameters.sh | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/hog_angle/compute_diameters.sh b/hog_angle/compute_diameters.sh index 76bbab0..153202a 100644 --- a/hog_angle/compute_diameters.sh +++ b/hog_angle/compute_diameters.sh @@ -59,8 +59,7 @@ cd $PATH_DATA_PROCESSED # Copy source images # Note: we use '/./' in order to include the sub-folder 'ses-0X' -rsync -Ravzh --no-g $PATH_DATA/./$SUBJECT . -# • -R → Preserves relative path structure. +rsync -avzh $PATH_DATA/./$SUBJECT . # • -a → Archive mode (preserves most attributes except group ownership). # • -v → Verbose output. # • -z → Compresses data during transfer. From ab689f6409ded9a6252c0d3973b60071e3c5d533 Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 15:22:46 +0200 Subject: [PATCH 10/17] make sure images are gzip files --- hog_angle/compute_diameters.sh | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/hog_angle/compute_diameters.sh b/hog_angle/compute_diameters.sh index 153202a..0b43f06 100644 --- a/hog_angle/compute_diameters.sh +++ b/hog_angle/compute_diameters.sh @@ -54,6 +54,11 @@ start=`date +%s` # Display useful info for the log, such as SCT version, RAM and CPU cores available sct_check_dependencies -short +cd $PATH_DATA/$SUBJECT/anat +git annex get *T2w* +cd $PATH_DATA/derivatives/labels/$SUBJECT/anat +git annex get *T2w* + # Go to folder where data will be copied and processed cd $PATH_DATA_PROCESSED From 602890f7911f507589d4ab439b329b9a1162352f Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 15:23:01 +0200 Subject: [PATCH 11/17] add example usage --- hog_angle/compute_diameters.sh | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/hog_angle/compute_diameters.sh b/hog_angle/compute_diameters.sh index 0b43f06..ebc65fe 100644 --- a/hog_angle/compute_diameters.sh +++ b/hog_angle/compute_diameters.sh @@ -6,7 +6,7 @@ # SCT PR #4958: https://github.com/spinalcordtoolbox/spinalcordtoolbox/pull/4958 # # Usage: -# sct_run_batch -config config.json +# xargs -n2 sh -c 'sct_run_batch -path-data $1 -path-output $2 -config config_compute_diameters.json' sh < datasets.txt # # Example of config.json: # { @@ -15,8 +15,10 @@ # } # # Example of dataset.txt: -# path/to/dataset1 path/to/output1 -# path/to/dataset2 path/to/output1 +# ~/data/data.neuro.polymtl.ca/data-multi-subject/ ~/results/hog_angle/compute_diameters_2025-07-29 +# ~/data/data.neuro.polymtl.ca/whole-spine ~/results/hog_angle/compute_diameters_2025-07-29 +# +# NOTE: we use the same results folder for all datasets to store results from all datasets in a single folder. # # Manual segmentations and disc labels should be located under: # PATH_DATA/derivatives/labels/SUBJECT// From 83aac0ead26392837e8a3c91a06dd8a26ed3e4b7 Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 15:23:18 +0200 Subject: [PATCH 12/17] add PAM50 normalization --- hog_angle/compute_diameters.sh | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/hog_angle/compute_diameters.sh b/hog_angle/compute_diameters.sh index ebc65fe..2c438cb 100644 --- a/hog_angle/compute_diameters.sh +++ b/hog_angle/compute_diameters.sh @@ -136,8 +136,13 @@ fi # Generate labeled segmentation using init disc labels sct_label_vertebrae -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -discfile ${FILEDISCS}.nii.gz -c t2 -qc ${PATH_QC} -qc-subject ${file} -sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -o ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr0.csv -qc qc -angle-corr 0 -sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -o ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr1.csv -qc qc -angle-corr 1 +# Normalize to PAM50, with and without the angle correction +sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -normalize-PAM50 1 -o ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr0.csv -qc ${PATH_QC} -angle-corr 0 +sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -normalize-PAM50 1 -o ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr1.csv -qc ${PATH_QC} -angle-corr 1 + +# No normalization to PAM50, with and without the angle correction +sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -o ${PATH_RESULTS}/${file}_metrics_angle_corr0.csv -qc ${PATH_QC} -angle-corr 0 +sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -o ${PATH_RESULTS}/${file}_metrics_angle_corr1.csv -qc ${PATH_QC} -angle-corr 1 # Create figures folder under ${PATH_RESULTS} if it does not exist mkdir -p ${PATH_RESULTS}/figures From f6b913fdd85c45b3b8e2b096edcc7a7e2e366538 Mon Sep 17 00:00:00 2001 From: valosekj Date: Tue, 29 Jul 2025 15:23:43 +0200 Subject: [PATCH 13/17] add plotting in the native space --- hog_angle/compute_diameters.sh | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/hog_angle/compute_diameters.sh b/hog_angle/compute_diameters.sh index 2c438cb..96b65ba 100644 --- a/hog_angle/compute_diameters.sh +++ b/hog_angle/compute_diameters.sh @@ -151,10 +151,20 @@ mkdir -p ${PATH_RESULTS}/figures ${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/hog_angle/AP_RL_diameters_PAM50.py \ -i ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr0.csv \ -o ${PATH_RESULTS}/figures/${file}_AP_RL_diameters_PAM50_angle_corr0.png \ + -smooth 0 ${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/hog_angle/AP_RL_diameters_PAM50.py \ -i ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr1.csv \ - -o ${PATH_RESULTS}/figures/${file}_AP_RL_diameters_PAM50_angle_corr1.png + -o ${PATH_RESULTS}/figures/${file}_AP_RL_diameters_PAM50_angle_corr1.png \ + -smooth 0 + +${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/hog_angle/AP_RL_diameters.py \ + -i ${PATH_RESULTS}/${file}_metrics_angle_corr0.csv \ + -o ${PATH_RESULTS}/figures/${file}_AP_RL_diameters_angle_corr0.png \ + +${SCT_DIR}/python/envs/venv_sct/bin/python ~/code/dcm-metric-normalization/hog_angle/AP_RL_diameters.py \ + -i ${PATH_RESULTS}/${file}_metrics_angle_corr1.csv \ + -o ${PATH_RESULTS}/figures/${file}_AP_RL_diameters_angle_corr1.png \ # ------------------------------------------------------------------------------ # End From 16dfd193d63f7303a276bd8872b0c865fc3415c5 Mon Sep 17 00:00:00 2001 From: valosekj Date: Mon, 11 Aug 2025 19:21:34 +0200 Subject: [PATCH 14/17] tweak METRICS_TO_YLIM --- hog_angle/AP_RL_diameters_PAM50.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hog_angle/AP_RL_diameters_PAM50.py b/hog_angle/AP_RL_diameters_PAM50.py index c442193..fbf1c38 100644 --- a/hog_angle/AP_RL_diameters_PAM50.py +++ b/hog_angle/AP_RL_diameters_PAM50.py @@ -30,8 +30,8 @@ # Set ylim to do not overlap horizontal grid with vertebrae labels METRICS_TO_YLIM = { - 'MEAN(diameter_AP)': (5.5, 9.5), - 'MEAN(diameter_RL)': (8.5, 14.5), + 'MEAN(diameter_AP)': (5, 10), + 'MEAN(diameter_RL)': (8, 16), } LABELS_FONT_SIZE = 14 From 07b4d4d1c87a860704cd9d50a54753ff3187f1c2 Mon Sep 17 00:00:00 2001 From: valosekj Date: Mon, 11 Aug 2025 19:21:46 +0200 Subject: [PATCH 15/17] add alpha to HOG-based curves --- hog_angle/AP_RL_diameters_PAM50.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hog_angle/AP_RL_diameters_PAM50.py b/hog_angle/AP_RL_diameters_PAM50.py index fbf1c38..3f207af 100644 --- a/hog_angle/AP_RL_diameters_PAM50.py +++ b/hog_angle/AP_RL_diameters_PAM50.py @@ -140,7 +140,7 @@ def create_lineplot(df, figure_path, smooth): label=metric) # HOG-based sns.lineplot(ax=axs[index], x="Slice (I->S)", y=f'{metric.replace(")", "_hog)")}', - data=df, linewidth=2, + data=df, linewidth=2, alpha=0.5, label=f'{metric.replace(")", "_hog)")}') # Tweak y-axis limits From d41ceee1fb20f3362df638ab87b3007c32fe251f Mon Sep 17 00:00:00 2001 From: valosekj Date: Mon, 11 Aug 2025 19:22:05 +0200 Subject: [PATCH 16/17] add `-v 2` to get debug figures for the HOG angle --- hog_angle/compute_diameters.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hog_angle/compute_diameters.sh b/hog_angle/compute_diameters.sh index 96b65ba..a6d627d 100644 --- a/hog_angle/compute_diameters.sh +++ b/hog_angle/compute_diameters.sh @@ -137,7 +137,7 @@ fi sct_label_vertebrae -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -discfile ${FILEDISCS}.nii.gz -c t2 -qc ${PATH_QC} -qc-subject ${file} # Normalize to PAM50, with and without the angle correction -sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -normalize-PAM50 1 -o ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr0.csv -qc ${PATH_QC} -angle-corr 0 +sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -normalize-PAM50 1 -o ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr0.csv -qc ${PATH_QC} -angle-corr 0 -v 2 sct_process_segmentation -i ${file_t2}.nii.gz -s ${FILESEG}.nii.gz -vertfile ${FILESEG}_labeled.nii.gz -perslice 1 -normalize-PAM50 1 -o ${PATH_RESULTS}/${file}_metrics_PAM50_angle_corr1.csv -qc ${PATH_QC} -angle-corr 1 # No normalization to PAM50, with and without the angle correction From 2dc657f42a4231e632e238e4477a414febe4686a Mon Sep 17 00:00:00 2001 From: valosekj Date: Wed, 13 Aug 2025 09:05:48 +0200 Subject: [PATCH 17/17] do not show the figure --- hog_angle/AP_RL_diameters.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hog_angle/AP_RL_diameters.py b/hog_angle/AP_RL_diameters.py index e0a5627..f5d8f42 100644 --- a/hog_angle/AP_RL_diameters.py +++ b/hog_angle/AP_RL_diameters.py @@ -108,8 +108,8 @@ def main() -> None: # Save figure plt.savefig(out_path, bbox_inches='tight', dpi=300) print(f"Figure saved: {out_path}") - plt.show() - # plt.close() + #plt.show() + plt.close() if __name__ == "__main__":