Changes made to sct_preprocess.py

Created Diff never expires
8 removals
920 lines
25 additions
932 lines
#!/usr/bin/env python3
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# -*- coding: utf-8 -*-
"""
"""
Created on Wednesday 24 may 2023
Created on Wednesday 24 may 2023


@author: evefra
@author: evefra
"""
"""


# import libraries
# import libraries
import os
import os
import shutil
import shutil
import subprocess
import subprocess
import fnmatch
import fnmatch
import nibabel as nib
import nibabel as nib
import numpy as np
import numpy as np


# work in conda env, where SCT is installed. For me: conda activate venv_sct
# work in conda env, where SCT is installed. For me: conda activate venv_sct


######################################################################################
######################################################################################
# INDICATE WHAT TO DO#
# INDICATE WHAT TO DO#
######################################################################################
######################################################################################


# FUNCTIONAL
# FUNCTIONAL
remove_volumes = 0
remove_volumes = 0
motion_correction = 0
motion_correction = 0
# MATLAB: slice time correction
# MATLAB: slice time correction
# MATLAB: mean slice time corrected image
# MATLAB: mean slice time corrected image
segment_mean = 0
segment_mean = 0


# ANATOMICAL T1
# ANATOMICAL T1
correct_qform = 0
correct_qform = 0
first_segment_t1 = 0
first_segment_t1 = 0
smooth = 0
smooth = 0
second_segment_t1 = 0
second_segment_t1 = 0
label_t1 = 0
label_t1 = 0
register_t1_temp = 0
register_t1_temp = 0


# ANATOMICAL T2*
# ANATOMICAL T2*
# MATLAB: transfer from dicom to nii
# MATLAB: transfer from dicom to nii
first_segment_t2 = 0
first_segment_t2 = 0


# COREGISTERING
# COREGISTERING
register_t2_t1 = 0
register_t2_t1 = 0
register_fmri_t2 = 0
register_fmri_t2 = 0
t2_to_t1_mult_t1_to_temp = 0
t2_to_t1_mult_t1_to_temp = 0
correct_qform_sc = 0
correct_qform_sc = 0
gm_segment = 0
gm_segment = 0
gm_segment_add = 0
gm_segment_add = 1
epi_to_t2_mult_t2_to_temp = 0
epi_to_t2_mult_t2_to_temp = 0
final_warp_apply = 0
final_warp_apply = 0


# LOOP OVER SUB NUMBERS
# LOOP OVER SUB NUMBERS
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------


sub_list = [10]
sub_list = [10]
# [1, 2, 4, 5, 6, 8, 9, 10, 11, 13, 14, 16, 17, 18, 24]
# [1, 2, 4, 5, 6, 8, 9, 10, 11, 13, 14, 16, 17, 18, 24]
for sub in sub_list:
for sub in sub_list:


# data path pre subject unspecific
# data path pre subject unspecific
data_pre = "/project/3023009.09/evefra"
data_pre = "./data/share"
subdirectory = f"sub-{sub:03d}"
subdirectory = f"sub-{sub:03d}"


# random testing
# random testing


# warp = '/project/3023009.09/evefra/sub-010/spinalcord/concat_t2_t1_temp/gm_warp_t22temp.nii.gz'
# warp = '/project/3023009.09/evefra/sub-010/spinalcord/concat_t2_t1_temp/gm_warp_t22temp.nii.gz'
# des_pam = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz'
# des_pam = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz'
# des_t1 = '/project/3023009.09/evefra/sub-013/anat/s302300909_sub013-0017-00001-000192-01.qform.nii'
# des_t1 = '/project/3023009.09/evefra/sub-013/anat/s302300909_sub013-0017-00001-000192-01.qform.nii'
# des_t2 = '/project/3023009.09/evefra/sub-013/anat_sc/s302300909_sub013-0009-00001-000001-01.nii'
# des_t2 = '/project/3023009.09/evefra/sub-013/anat_sc/s302300909_sub013-0009-00001-000001-01.nii'
# input_im = '/project/3023009.09/evefra/sub-013/anat_sc/s302300909_sub013-0009-00001-000001-01.nii'
# input_im = '/project/3023009.09/evefra/sub-013/anat_sc/s302300909_sub013-0009-00001-000001-01.nii'
# des_mean = '/project/3023009.09/evefra/sub-013/spinalcord/slice_time_corrected/spinacord_slc_mean.nii'
# des_mean = '/project/3023009.09/evefra/sub-013/spinalcord/slice_time_corrected/spinacord_slc_mean.nii'


# command = f"sct_apply_transfo -i {des_t2} -d {des_pam} -w {warp}"
# command = f"sct_apply_transfo -i {des_t2} -d {des_pam} -w {warp}"
# process = subprocess.Popen(command, shell=True)
# process = subprocess.Popen(command, shell=True)
# process.wait()
# process.wait()


#################################################################################
#################################################################################
# FUNCTIONAL prep #
# FUNCTIONAL prep #
#################################################################################
#################################################################################


# try:
# try:
# REMOVE FIRST TWO VOLUMES
# REMOVE FIRST TWO VOLUMES
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
if remove_volumes == 1:
if remove_volumes == 1:
directory = os.path.join(
directory = os.path.join(
data_pre, subdirectory, "spinalcord", "niftie")
data_pre, subdirectory, "spinalcord", "niftie")
removed_directory = os.path.join(
removed_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "removed_volumes")
data_pre, subdirectory, "spinalcord", "removed_volumes")
file_extension = ".nii"
file_extension = ".nii"


# Check if the removed volumes directory exists, create it if necessary
# Check if the removed volumes directory exists, create it if necessary
if not os.path.exists(removed_directory):
if not os.path.exists(removed_directory):
os.makedirs(removed_directory)
os.makedirs(removed_directory)


# Get a list of all files in the directory with the specified extension
# Get a list of all files in the directory with the specified extension
file_list = [f for f in os.listdir(
file_list = [f for f in os.listdir(
directory) if f.endswith(file_extension)]
directory) if f.endswith(file_extension)]


# Sort the file list alphabetically
# Sort the file list alphabetically
file_list.sort()
file_list.sort()


# Move the first two files to the removed volumes directory
# Move the first two files to the removed volumes directory
files_to_move = file_list[:2]
files_to_move = file_list[:2]


# Move the files one by one
# Move the files one by one
for file_name in files_to_move:
for file_name in files_to_move:
source_path = os.path.join(directory, file_name)
source_path = os.path.join(directory, file_name)
destination_path = os.path.join(removed_directory, file_name)
destination_path = os.path.join(removed_directory, file_name)


# Move the file to the destination
# Move the file to the destination
shutil.move(source_path, destination_path)
shutil.move(source_path, destination_path)


# Check if only two files in the removed-directory
# Check if only two files in the removed-directory
file_list = os.listdir(removed_directory)
file_list = os.listdir(removed_directory)


if len(file_list) != 2:
if len(file_list) != 2:
raise ValueError(
raise ValueError(
"The directory does not contain exactly two files.")
"The directory does not contain exactly two files.")
print('remove volumes completed for ', subdirectory)
print('remove volumes completed for ', subdirectory)


# MOTION CORRECTION
# MOTION CORRECTION
# ----------------------------------------------------------------------------------------
# ----------------------------------------------------------------------------------------
if motion_correction == 1:
if motion_correction == 1:


input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "niftie")
data_pre, subdirectory, "spinalcord", "niftie")
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "motion_corrected")
data_pre, subdirectory, "spinalcord", "motion_corrected")


# Check if the removed volumes directory exists, create it if necessary
# Check if the removed volumes directory exists, create it if necessary
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)


# Get a list of all NIfTI files in the input directory
# Get a list of all NIfTI files in the input directory
file_list_before = os.listdir(input_directory)
file_list_before = os.listdir(input_directory)
file_list = []
file_list = []
for names in file_list_before:
for names in file_list_before:
if names.endswith(".nii"):
if names.endswith(".nii"):
file_list.append(names)
file_list.append(names)


# Perform motion correction for each input file
# Perform motion correction for each input file
for input_file in file_list:
for input_file in file_list:
# Construct input and output file paths
# Construct input and output file paths
input_path = os.path.join(input_directory, input_file)
input_path = os.path.join(input_directory, input_file)
output_file = input_file.replace(".nii.gz", "_moco.nii.gz")
output_file = input_file.replace(".nii.gz", "_moco.nii.gz")
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)


# Run motion correction using sct_fmri_moco
# Run motion correction using sct_fmri_moco
command = f"sct_fmri_moco -i {input_path} -o {output_path}"
command = f"sct_fmri_moco -i {input_path} -o {output_path}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()


# print if failed
# print if failed
if process.returncode != 0:
if process.returncode != 0:
print(f"Motion correction failed for {input_file}.")
print(f"Motion correction failed for {input_file}.")
print('moco completed for ', subdirectory)
print('moco completed for ', subdirectory)


# CHEKC MOTION CORRECTION AND CORRECT
# CHEKC MOTION CORRECTION AND CORRECT
# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------


# check FSLeyes
# check FSLeyes


# SEGMENT MEAN slc IMAGE
# SEGMENT MEAN slc IMAGE
# -----------------------------------------------------------------
# -----------------------------------------------------------------


if segment_mean == 1:
if segment_mean == 1:
# Define the input
# Define the input
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "spinalcord", 'slice_time_corrected')
data_pre, subdirectory, "spinalcord", 'slice_time_corrected')
wildcard_pattern = 'spinacord_slc_mean.nii'
wildcard_pattern = 'spinacord_slc_mean.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path_fmri = os.path.join(input_directory, input_file[0])
input_path_fmri = os.path.join(input_directory, input_file[0])


# define output
# define output
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "segment_slc_mean")
data_pre, subdirectory, "spinalcord", "segment_slc_mean")
# Create the output directory if it doesn't exist
# Create the output directory if it doesn't exist
os.makedirs(output_directory, exist_ok=True)
os.makedirs(output_directory, exist_ok=True)


# Define the quality control directory
# Define the quality control directory
qc_directory = os.path.join(output_directory, 'qc')
qc_directory = os.path.join(output_directory, 'qc')


# Create the output directory if it doesn't exist
# Create the output directory if it doesn't exist
os.makedirs(output_directory, exist_ok=True)
os.makedirs(output_directory, exist_ok=True)


# segmentation
# segmentation
command = f'sct_deepseg_sc -i {input_path_fmri} -c t2 -ofolder {output_directory} -qc {qc_directory} -centerline viewer'
command = f'sct_deepseg_sc -i {input_path_fmri} -c t2 -ofolder {output_directory} -qc {qc_directory} -centerline viewer'
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()
print('segment mean func completed for ', subdirectory)
print('segment mean func completed for ', subdirectory)


##############################################################################
##############################################################################
# T2 prep #
# T2 prep #
#################################################################################
#################################################################################


# SEGMENTAION OF T2 (including manual C2 and T1 labelling)
# SEGMENTAION OF T2 (including manual C2 and T1 labelling)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------


if first_segment_t2 == 1:
if first_segment_t2 == 1:


# define input and output
# define input and output
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "anat_sc")
data_pre, subdirectory, "anat_sc")
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "segment_anat_t2")
data_pre, subdirectory, "spinalcord", "segment_anat_t2")


# Check if the removed volumes directory exists, create it if necessary
# Check if the removed volumes directory exists, create it if necessary
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)
qc_directory = os.path.join(output_directory, 'qc')
qc_directory = os.path.join(output_directory, 'qc')


# Find the input file in the "anat" folder
# Find the input file in the "anat" folder
import fnmatch
import fnmatch


# Define the wildcard pattern
# Define the wildcard pattern
wildcard_pattern = 's*.nii'
wildcard_pattern = 's*.nii'


# Get the list of files matching the wildcard pattern in the source directory
# Get the list of files matching the wildcard pattern in the source directory
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)


# Construct input and output file paths
# Construct input and output file paths
input_path = os.path.join(input_directory, input_file[0])
input_path = os.path.join(input_directory, input_file[0])
output_file = input_file[0].replace("nii", "segmentation.nii")
output_file = input_file[0].replace("nii", "segmentation.nii")
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)


command = f"sct_deepseg_sc -i {input_path} -c t2 -o {output_path} -qc {qc_directory} -centerline viewer"
command = f"sct_deepseg_sc -i {input_path} -c t2 -o {output_path} -qc {qc_directory} -centerline viewer"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()
print('first segment t2 completed for ', subdirectory)
print('first segment t2 completed for ', subdirectory)


# CHECK SEGMENTATION AND CORRECT
# CHECK SEGMENTATION AND CORRECT
# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------


# in FSLeyes
# in FSLeyes


#################################################################################
#################################################################################
# T1 prep #
# T1 prep #
#################################################################################
#################################################################################


# CORRECT QFORM ANATOMICAL
# CORRECT QFORM ANATOMICAL
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------


if correct_qform == 1:
if correct_qform == 1:
# Define the input
# Define the input
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "anat")
data_pre, subdirectory, "anat")
wildcard_pattern = 's*.nii'
wildcard_pattern = 's*.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path = os.path.join(input_directory, input_file[0])
input_path = os.path.join(input_directory, input_file[0])


# define output
# define output
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "anat")
data_pre, subdirectory, "anat")
output_file = input_file[0].replace("nii", "qform.nii")
output_file = input_file[0].replace("nii", "qform.nii")
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)


# run command
# run command
command = f"sct_image -i {input_path} -set-sform-to-qform -o {output_path}"
command = f"sct_image -i {input_path} -set-sform-to-qform -o {output_path}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()


print('qform fix anat completed for ', subdirectory)
print('qform fix anat completed for ', subdirectory)


# FIRST SEGMENTAION OF T1 (including manual C2 and T1 labelling)
# FIRST SEGMENTAION OF T1 (including manual C2 and T1 labelling)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------


if first_segment_t1 == 1:
if first_segment_t1 == 1:


# define input
# define input
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "anat")
data_pre, subdirectory, "anat")
wildcard_pattern = 's*qform.nii'
wildcard_pattern = 's*qform.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path = os.path.join(input_directory, input_file[0])
input_path = os.path.join(input_directory, input_file[0])


# define output
# define output
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "segment_anat")
data_pre, subdirectory, "spinalcord", "segment_anat")
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)
output_file = input_file[0].replace("nii", "segmentation.nii")
output_file = input_file[0].replace("nii", "segmentation.nii")
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)


# define QC directory
# define QC directory
qc_directory = os.path.join(output_directory, 'qc')
qc_directory = os.path.join(output_directory, 'qc')


# run first segmentation
# run first segmentation
# label centerline C2 and C7
# label centerline C2 and C7
command = f"sct_deepseg_sc -i {input_path} -c t1 -o {output_path} -qc {qc_directory} -centerline viewer"
command = f"sct_deepseg_sc -i {input_path} -c t1 -o {output_path} -qc {qc_directory} -centerline viewer"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()
print('first segment t1 completed for ', subdirectory)
print('first segment t1 completed for ', subdirectory)


# SMOOTH ALONG SC TO IMPROVE SEGMENTATION
# SMOOTH ALONG SC TO IMPROVE SEGMENTATION
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
if smooth == 1:
if smooth == 1:
# define input
# define input
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "anat")
data_pre, subdirectory, "anat")
wildcard_pattern = 's*qform.nii'
wildcard_pattern = 's*qform.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path = os.path.join(input_directory, input_file[0])
input_path = os.path.join(input_directory, input_file[0])


# select segmentation image
# select segmentation image
input_directory_seg = os.path.join(
input_directory_seg = os.path.join(
data_pre, subdirectory, 'spinalcord', "segment_anat")
data_pre, subdirectory, 'spinalcord', "segment_anat")
wildcard_pattern = '*.segmentation.nii'
wildcard_pattern = '*.segmentation.nii'
input_file_seg = fnmatch.filter(os.listdir(
input_file_seg = fnmatch.filter(os.listdir(
input_directory_seg), wildcard_pattern)
input_directory_seg), wildcard_pattern)
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])


# define output
# define output
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "segment_anat")
data_pre, subdirectory, "spinalcord", "segment_anat")
output_file_seg = input_file_seg[0].replace(
output_file_seg = input_file_seg[0].replace(
"segmentation.nii", "smooth.nii")
"segmentation.nii", "smooth.nii")
output_path_seg = os.path.join(output_directory, output_file_seg)
output_path_seg = os.path.join(output_directory, output_file_seg)


# RUN SMOOOTHING
# RUN SMOOOTHING
command = f"sct_smooth_spinalcord -i {input_path} -o {output_path_seg} -s {input_path_seg} -smooth 0,0,5"
command = f"sct_smooth_spinalcord -i {input_path} -o {output_path_seg} -s {input_path_seg} -smooth 0,0,5"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()
print('smooth t1 completed for ', subdirectory)
print('smooth t1 completed for ', subdirectory)


# DO SECOND SEGMENTATION
# DO SECOND SEGMENTATION
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------


if second_segment_t1 == 1:
if second_segment_t1 == 1:


# define input
# define input
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "segment_anat")
data_pre, subdirectory, "spinalcord", "segment_anat")
wildcard_pattern = '*smooth.nii'
wildcard_pattern = '*smooth.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path = os.path.join(input_directory, input_file[0])
input_path = os.path.join(input_directory, input_file[0])


# define output
# define output
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "segment_anat_second")
data_pre, subdirectory, "spinalcord", "segment_anat_second")
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)
output_file = input_file[0].replace(
output_file = input_file[0].replace(
"smooth.nii", "second_segmentation.nii")
"smooth.nii", "second_segmentation.nii")
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)
# define qcdirectory
# define qcdirectory
qc_directory = os.path.join(output_directory, 'qc')
qc_directory = os.path.join(output_directory, 'qc')


# select segmentation image
# select segmentation image
input_directory_seg = os.path.join(
input_directory_seg = os.path.join(
data_pre, subdirectory, 'spinalcord', "segment_anat")
data_pre, subdirectory, 'spinalcord', "segment_anat")
wildcard_pattern = '*.segmentation.nii'
wildcard_pattern = '*.segmentation.nii'
input_file_seg = fnmatch.filter(os.listdir(
input_file_seg = fnmatch.filter(os.listdir(
input_directory_seg), wildcard_pattern)
input_directory_seg), wildcard_pattern)
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])


# run second segmentation
# run second segmentation
command = f"sct_deepseg_sc -i {input_path} -c t1 -o {output_path} -qc {qc_directory} -centerline file -file_centerline {input_path_seg}"
command = f"sct_deepseg_sc -i {input_path} -c t1 -o {output_path} -qc {qc_directory} -centerline file -file_centerline {input_path_seg}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()
print('second segment t1 completed for ', subdirectory)
print('second segment t1 completed for ', subdirectory)


# CHECK SEGMENTATION AND CORRECT
# CHECK SEGMENTATION AND CORRECT
# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------


# in FSLeyes
# in FSLeyes


# LABEL T1
# LABEL T1
# -------------------------------------------------------------------------------
# -------------------------------------------------------------------------------
if label_t1 == 1:
if label_t1 == 1:


# define input anatomical image
# define input anatomical image
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "anat")
data_pre, subdirectory, "anat")
wildcard_pattern = 's3*qform.nii'
wildcard_pattern = 's3*qform.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path = os.path.join(input_directory, input_file[0])
input_path = os.path.join(input_directory, input_file[0])


# define input segmentation
# define input segmentation
input_directory_seg = os.path.join(
input_directory_seg = os.path.join(
data_pre, subdirectory, 'spinalcord', "segment_anat_second")
data_pre, subdirectory, 'spinalcord', "segment_anat_second")
wildcard_pattern = '*segmentation.nii'
wildcard_pattern = '*segmentation.nii'
input_file_seg = fnmatch.filter(os.listdir(
input_file_seg = fnmatch.filter(os.listdir(
input_directory_seg), wildcard_pattern)
input_directory_seg), wildcard_pattern)
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])


# define output paths and name
# define output paths and name
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "label_anat")
data_pre, subdirectory, "spinalcord", "label_anat")
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)
output_file_manual = input_file[0].replace(".nii", "manlabel.nii")
output_file_manual = input_file[0].replace(".nii", "manlabel.nii")
output_path_manual = os.path.join(output_directory, output_file_manual)
output_path_manual = os.path.join(output_directory, output_file_manual)
output_file = input_file[0].replace(".nii", "label.nii")
output_file = input_file[0].replace(".nii", "label.nii")
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)


# define qc directory
# define qc directory
qc_directory_man = os.path.join(output_directory, 'qc_manual')
qc_directory_man = os.path.join(output_directory, 'qc_manual')
qc_directory = os.path.join(output_directory, 'qc')
qc_directory = os.path.join(output_directory, 'qc')


# manual label C2 and C7
# manual label C2 and C7
command = f"sct_label_utils -i {input_path} -qc {qc_directory_man} -create-viewer 3,7 -o {output_path_manual} -msg 'Please manually label the inter-vertebral disc '"
command = f"sct_label_utils -i {input_path} -qc {qc_directory_man} -create-viewer 3,7 -o {output_path_manual} -msg 'Please manually label the inter-vertebral disc '"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()


# define input man label
# define input man label
# wildcard_pattern = '*manlabel.nii'
# wildcard_pattern = '*manlabel.nii'
# input_file_man = fnmatch.filter(os.listdir(
# input_file_man = fnmatch.filter(os.listdir(
# output_directory), wildcard_pattern)
# output_directory), wildcard_pattern)
# input_path_man = os.path.join(output_directory, input_file_man[0])
# input_path_man = os.path.join(output_directory, input_file_man[0])


# correcting s or qform if needed
# correcting s or qform if needed
# sct_image -set-sform-to-qform -i {
# sct_image -set-sform-to-qform -i {


# run labelling
# run labelling
# command = f"sct_label_vertebrae -i {input_path} -c t1 -s {input_path_seg} -qc {qc_directory} -o {output_path} -initlabel {input_path_man}"
# command = f"sct_label_vertebrae -i {input_path} -c t1 -s {input_path_seg} -qc {qc_directory} -o {output_path} -initlabel {input_path_man}"
# process = subprocess.Popen(command, shell=True)
# process = subprocess.Popen(command, shell=True)
# process.wait()
# process.wait()
print('label t1 completed for ', subdirectory)
print('label t1 completed for ', subdirectory)


# CHECK LABELLING
# CHECK LABELLING
# --------------------------------------------------------------------------------
# --------------------------------------------------------------------------------


# manual label more vertebra
# manual label more vertebra
# repeat steps
# repeat steps


# command = f"sct_label_utils -i {input_path} -qc {qc_directory_man} -create-viewer 3,4,5,6,7 -o {output_path_manual} -msg 'Please manually label the inter-vertebral disc 3 and 7.'"
# command = f"sct_label_utils -i {input_path} -qc {qc_directory_man} -create-viewer 3,4,5,6,7 -o {output_path_manual} -msg 'Please manually label the inter-vertebral disc 3 and 7.'"
# process = subprocess.Popen(command, shell=True)
# process = subprocess.Popen(command, shell=True)
# process.wait()
# process.wait()


# REGISTER FROM T1 to TEMPLATE
# REGISTER FROM T1 to TEMPLATE
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------


if register_t1_temp == 1:
if register_t1_temp == 1:
# define input image (T1)
# define input image (T1)
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "anat")
data_pre, subdirectory, "anat")
wildcard_pattern = 's*qform.nii'
wildcard_pattern = 's*qform.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path = os.path.join(input_directory, input_file[0])
input_path = os.path.join(input_directory, input_file[0])


# define input segmentation
# define input segmentation
input_directory_seg = os.path.join(
input_directory_seg = os.path.join(
data_pre, subdirectory, 'spinalcord', "segment_anat_second")
data_pre, subdirectory, 'spinalcord', "segment_anat_second")
wildcard_pattern = '*segmentation.nii'
wildcard_pattern = '*segmentation.nii'
input_file_seg = fnmatch.filter(os.listdir(
input_file_seg = fnmatch.filter(os.listdir(
input_directory_seg), wildcard_pattern)
input_directory_seg), wildcard_pattern)
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])


# define labels input image
# define labels input image
input_directory_lab = os.path.join(
input_directory_lab = os.path.join(
data_pre, subdirectory, "spinalcord", "label_anat")
data_pre, subdirectory, "spinalcord", "label_anat")
wildcard_pattern = 's*qformmanlabel.nii' # name if not manual labelled?
wildcard_pattern = 's*qformmanlabel.nii' # name if not manual labelled?
input_file_lab = fnmatch.filter(os.listdir(
input_file_lab = fnmatch.filter(os.listdir(
input_directory_lab), wildcard_pattern)
input_directory_lab), wildcard_pattern)
input_path_lab = os.path.join(input_directory_lab, input_file_lab[0])
input_path_lab = os.path.join(input_directory_lab, input_file_lab[0])


# define output directory
# define output directory
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "register_t1_temp")
data_pre, subdirectory, "spinalcord", "register_t1_temp")
# define qc
# define qc
qc_directory = os.path.join(output_directory, 'qc_t1_temp')
qc_directory = os.path.join(output_directory, 'qc_t1_temp')


# run registering T1 --> template
# run registering T1 --> template
command = f"sct_register_to_template -i {input_path} -s {input_path_seg} -l {input_path_lab} -c t1 -qc {qc_directory} -ofolder {output_directory}"
command = f"sct_register_to_template -i {input_path} -s {input_path_seg} -l {input_path_lab} -c t1 -qc {qc_directory} -ofolder {output_directory}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()
print('register calc t1 - temp completed for ', subdirectory)
print('register calc t1 - temp completed for ', subdirectory)


#################################################################################
#################################################################################
# COREGISTER #
# COREGISTER #
#################################################################################
#################################################################################


# REGISTER FROM T2 to T1
# REGISTER FROM T2 to T1
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------


if register_t2_t1 == 1:
if register_t2_t1 == 1:
# define input image (T2)
# define input image (T2)
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "anat_sc")
data_pre, subdirectory, "anat_sc")
wildcard_pattern = 's3*01.nii'
wildcard_pattern = 's3*01.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path = os.path.join(input_directory, input_file[0])
input_path = os.path.join(input_directory, input_file[0])


# define input segmentation of T2
# define input segmentation of T2
input_directory_seg = os.path.join(
input_directory_seg = os.path.join(
data_pre, subdirectory, 'spinalcord', "segment_anat_t2")
data_pre, subdirectory, 'spinalcord', "segment_anat_t2")
wildcard_pattern = '*segmentation.nii'
wildcard_pattern = '*segmentation.nii'
input_file_seg = fnmatch.filter(os.listdir(
input_file_seg = fnmatch.filter(os.listdir(
input_directory_seg), wildcard_pattern)
input_directory_seg), wildcard_pattern)
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])
input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])


# define input segmentation of T1
# define input segmentation of T1
input_directory_seg = os.path.join(
input_directory_seg = os.path.join(
data_pre, subdirectory, 'spinalcord', "segment_anat_second")
data_pre, subdirectory, 'spinalcord', "segment_anat_second")
wildcard_pattern = '*segmentation.nii'
wildcard_pattern = '*segmentation.nii'
input_file_seg = fnmatch.filter(os.listdir(
input_file_seg = fnmatch.filter(os.listdir(
input_directory_seg), wildcard_pattern)
input_directory_seg), wildcard_pattern)
input_path_seg_t1 = os.path.join(
input_path_seg_t1 = os.path.join(
input_directory_seg, input_file_seg[0])
input_directory_seg, input_file_seg[0])


# define labels input image (from t1?)
# define labels input image (from t1?)
input_directory_lab = os.path.join(
input_directory_lab = os.path.join(
data_pre, subdirectory, "spinalcord", "label_anat")
data_pre, subdirectory, "spinalcord", "label_anat")
wildcard_pattern = 's*manlabel.nii' # name if not manual labelled?
wildcard_pattern = 's*manlabel.nii' # name if not manual labelled?
input_file_lab = fnmatch.filter(os.listdir(
input_file_lab = fnmatch.filter(os.listdir(
input_directory_lab), wildcard_pattern)
input_directory_lab), wildcard_pattern)
input_path_lab = os.path.join(input_directory_lab, input_file_lab[0])
input_path_lab = os.path.join(input_directory_lab, input_file_lab[0])


# define to-be-warped image (T1)
# define to-be-warped image (T1)
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "anat")
data_pre, subdirectory, "anat")
wildcard_pattern = 's*qform.nii'
wildcard_pattern = 's*qform.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path_warped = os.path.join(input_directory, input_file[0])
input_path_warped = os.path.join(input_directory, input_file[0])


# define output directory
# define output directory
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "register_t2_t1")
data_pre, subdirectory, "spinalcord", "register_t2_t1")
# define qc
# define qc
qc_directory = os.path.join(output_directory, 'qc_t1_temp')
qc_directory = os.path.join(output_directory, 'qc_t1_temp')


# run coregistering T2 --> T1
# run coregistering T2 --> T1
command = f"sct_register_multimodal -i {input_path} -d {input_path_warped} -iseg {input_path_seg} -dseg {input_path_seg_t1} -qc {qc_directory} -ofolder {output_directory}"
command = f"sct_register_multimodal -i {input_path} -d {input_path_warped} -iseg {input_path_seg} -dseg {input_path_seg_t1} -qc {qc_directory} -ofolder {output_directory}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()
print('Register calc t2 - t1 completed for ', subdirectory)
print('Register calc t2 - t1 completed for ', subdirectory)


# CALCULATE warping field from fmri to T2*
# CALCULATE warping field from fmri to T2*
# ---------------------------------------------------------------------------------
# ---------------------------------------------------------------------------------
if register_fmri_t2 == 1:
if register_fmri_t2 == 1:
# SOURCE define input image (slc fmri mean)
# SOURCE define input image (slc fmri mean)
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "spinalcord", 'slice_time_corrected')
data_pre, subdirectory, "spinalcord", 'slice_time_corrected')
wildcard_pattern = 'spinacord_slc_mean.nii'
wildcard_pattern = 'spinacord_slc_mean.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path_fmri = os.path.join(input_directory, input_file[0])
input_path_fmri = os.path.join(input_directory, input_file[0])


# Source SEGMENT define input segmentation of fmri mean
# Source SEGMENT define input segmentation of fmri mean
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "segment_slc_mean")
data_pre, subdirectory, "spinalcord", "segment_slc_mean")
wildcard_pattern = 'spinacord_slc_mean_seg.nii'
wildcard_pattern = 'spinacord_slc_mean_seg.nii'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)
input_path_fmri_seg = os.path.join(input_directory, input_file[0])
input_path_fmri_seg = os.path.join(input_directory, input_file[0])


# DESTINATION define input of T2
# DESTINATION define input of T2
input_directory_d = os.path.join(
input_directory_d = os.path.join(
data_pre, subdirectory, "anat_sc")
data_pre, subdirectory, "anat_sc")
wildcard_pattern_d = 's3*01.nii'
wildcard_pattern_d = 's3*01.nii'
input_file_d = fnmatch.filter(os.listdir(
input_file_d = fnmatch.filter(os.listdir(
input_directory_d), wildcard_pattern_d)
input_directory_d), wildcard_pattern_d)
input_path_des = os.path.join(input_directory_d, input_file_d[0])
input_path_des = os.path.join(input_directory_d, input_file_d[0])


# DESTINATION SEGMENT define input segmentation of T2
# DESTINATION SEGMENT define input segmentation of T2
input_directory_seg = os.path.join(
input_directory_seg = os.path.join(
data_pre, subdirectory, 'spinalcord', "segment_anat_t2")
data_pre, subdirectory, 'spinalcord', "segment_anat_t2")
wildcard_pattern = '*segmentation.nii'
wildcard_pattern = '*segmentation.nii'
input_file_seg = fnmatch.filter(os.listdir(
input_file_seg = fnmatch.filter(os.listdir(
input_directory_seg), wildcard_pattern)
input_directory_seg), wildcard_pattern)
input_path_des_seg = os.path.join(
input_path_des_seg = os.path.join(
input_directory_seg, input_file_seg[0])
input_directory_seg, input_file_seg[0])


# define output directory
# define output directory
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "registration_fmri_t2")
data_pre, subdirectory, "spinalcord", "registration_fmri_t2")
# define qc
# define qc
qc_directory = os.path.join(output_directory, 'qc')
qc_directory = os.path.join(output_directory, 'qc')


# run coregistering T2 --> T1
# run coregistering T2 --> T1
command = f"sct_register_multimodal -i {input_path_fmri} -d {input_path_des} -dseg {input_path_des_seg} -qc {qc_directory} -ofolder {output_directory}"
command = f"sct_register_multimodal -i {input_path_fmri} -d {input_path_des} -dseg {input_path_des_seg} -qc {qc_directory} -ofolder {output_directory}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()
print('Register calc fmri - t2 completed for ', subdirectory)
print('Register calc fmri - t2 completed for ', subdirectory)


# mult t2-t1 to t1-temp
# mult t2-t1 to t1-temp
# ------------------------------------------------------------------------------------
# ------------------------------------------------------------------------------------
if t2_to_t1_mult_t1_to_temp == 1:
if t2_to_t1_mult_t1_to_temp == 1:


# Load the first warping field
# Load the first warping field
input_directory_w2 = os.path.join(
input_directory_w2 = os.path.join(
data_pre, subdirectory, "spinalcord", "register_t2_t1")
data_pre, subdirectory, "spinalcord", "register_t2_t1")
wildcard_pattern_w2 = '*qform2s*'
wildcard_pattern_w2 = '*qform2s*'
input_file_w2 = fnmatch.filter(os.listdir(
input_file_w2 = fnmatch.filter(os.listdir(
input_directory_w2), wildcard_pattern_w2)
input_directory_w2), wildcard_pattern_w2)
input_path_w2 = os.path.join(input_directory_w2, input_file_w2[0])
input_path_w2 = os.path.join(input_directory_w2, input_file_w2[0])


# define warping field t1 - temp
# define warping field t1 - temp
input_directory_w3 = os.path.join(
input_directory_w3 = os.path.join(
data_pre, subdirectory, "spinalcord", "register_t1_temp")
data_pre, subdirectory, "spinalcord", "register_t1_temp")
wildcard_pattern_w3 = 'warp_anat2template.nii.gz'
wildcard_pattern_w3 = 'warp_anat2template.nii.gz'
input_file_w3 = fnmatch.filter(os.listdir(
input_file_w3 = fnmatch.filter(os.listdir(
input_directory_w3), wildcard_pattern_w3)
input_directory_w3), wildcard_pattern_w3)
input_path_w3 = os.path.join(input_directory_w3, input_file_w3[0])
input_path_w3 = os.path.join(input_directory_w3, input_file_w3[0])


# destination image hardcoded
# destination image hardcoded
des_image = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz'
des_image = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz'


# output file
# output file
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)
output_file = 'warp_t22temp.nii.gz'
output_file = 'warp_t22temp.nii.gz'
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)


# run command
# run command
command = f"sct_concat_transfo -d {des_image} -o {output_path} -w {input_path_w2} {input_path_w3}"
command = f"sct_concat_transfo -d {des_image} -o {output_path} -w {input_path_w2} {input_path_w3}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()


################################################################################################
################################################################################################
# COREGISTER: ADD GM
# COREGISTER: ADD GM
#################################################################################################
#################################################################################################


# CORRECT QFORM ANATOMICAL
# CORRECT QFORM ANATOMICAL
# -----------------------------------------------------------------------------
# -----------------------------------------------------------------------------


if correct_qform_sc == 1:
if correct_qform_sc == 1:
# Define the input
# Define the input
input_directory_d = os.path.join(
input_directory_d = os.path.join(
data_pre, subdirectory, "anat_sc")
data_pre, subdirectory, "anat_sc")
wildcard_pattern_d = 's3*01.nii'
wildcard_pattern_d = 's3*01.nii'
input_file_d = fnmatch.filter(os.listdir(
input_file_d = fnmatch.filter(os.listdir(
input_directory_d), wildcard_pattern_d)
input_directory_d), wildcard_pattern_d)
input_path_des = os.path.join(input_directory_d, input_file_d[0])
input_path_des = os.path.join(input_directory_d, input_file_d[0])


# define output
# define output
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "anat_sc")
data_pre, subdirectory, "anat_sc")
output_file = input_file_d[0].replace("nii", "qform.nii")
output_file = input_file_d[0].replace("nii", "qform.nii")
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)


# run command
# run command
command = f"sct_image -i {input_path_des} -set-sform-to-qform -o {output_path}"
command = f"sct_image -i {input_path_des} -set-sform-to-qform -o {output_path}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()


print('qform fix anat sc completed for ', subdirectory)
print('qform fix anat sc completed for ', subdirectory)


# CALC SEGMENTATION GM
# CALC SEGMENTATION GM
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------


if gm_segment == 1:
if gm_segment == 1:


# define input
# define input
input_directory_d = os.path.join(
input_directory_d = os.path.join(
data_pre, subdirectory, "anat_sc")
data_pre, subdirectory, "anat_sc")
wildcard_pattern_d = 's3*01.qform.nii'
wildcard_pattern_d = 's3*01.qform.nii'
input_file_d = fnmatch.filter(os.listdir(
input_file_d = fnmatch.filter(os.listdir(
input_directory_d), wildcard_pattern_d)
input_directory_d), wildcard_pattern_d)
input_path_des = os.path.join(input_directory_d, input_file_d[0])
input_path_des = os.path.join(input_directory_d, input_file_d[0])


# define output
# define output
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "gm_segment_t2_anat")
data_pre, subdirectory, "spinalcord", "gm_segment_t2_anat")
output_file_seg = 'gm_t2.nii'
output_file_seg = 'gm_t2.nii'
output_path = os.path.join(output_directory, output_file_seg)
output_path = os.path.join(output_directory, output_file_seg)


# Check if the removed volumes directory exists, create it if necessary
# Check if the removed volumes directory exists, create it if necessary
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)


# define qc
# define qc
qc_directory = os.path.join(output_directory, 'qc')
qc_directory = os.path.join(output_directory, 'qc')


# Run gm segment
# Run gm segment
command = f"sct_deepseg_gm -i {input_path_des} -qc {qc_directory} -o {output_path}"
print('segmenting gm', subdirectory)
process = subprocess.Popen(command, shell=True)
from spinalcordtoolbox.scripts import sct_deepseg_gm
process.wait()
sct_deepseg_gm.main(['-i', input_path_des, '-qc', qc_directory, '-o', output_path])
# process = subprocess.Popen(command, shell=True)
# process.wait()
print('segment gm completed for ', subdirectory)
print('segment gm completed for ', subdirectory)


# ADD SEGMENTATION GM
# ADD SEGMENTATION GM
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
if gm_segment_add == 1:
if gm_segment_add == 1:


# define input template (always the same)
# define input template (always the same)
input_directory = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz'
input_directory = './data/PAM50/template/PAM50_t2.nii.gz'


# define input gm segmentation (always the same)
# define input gm segmentation (always the same)
input_directory_iseg = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_gm.nii.gz'
input_directory_iseg = './data/PAM50/template/PAM50_wm.nii.gz'


# define destination (T2)
# define destination (T2)
input_directory_d = os.path.join(
input_directory_d = os.path.join(
data_pre, subdirectory, "anat_sc")
data_pre, subdirectory, "anat_sc")
wildcard_pattern_d = 's3*01.nii'
wildcard_pattern_d = 's3*01.qform.nii'
input_file_d = fnmatch.filter(os.listdir(
input_file_d = fnmatch.filter(os.listdir(
input_directory_d), wildcard_pattern_d)
input_directory_d), wildcard_pattern_d)
input_path_des = os.path.join(input_directory_d, input_file_d[0])
input_path_des = os.path.join(input_directory_d, input_file_d[0])


# define destination segmentation (T2 GM)
# define destination segmentation (T2 GM)
input_directory_des_seg = os.path.join(
input_directory_des_seg = os.path.join(
data_pre, subdirectory, "spinalcord", "gm_segment_t2_anat")
data_pre, subdirectory, "spinalcord", "gm_segment_t2_anat")
wildcard_pattern_des_seg = 'gm*.nii'
wildcard_pattern_des_seg = 'wm*.nii'
input_file_des_seg = fnmatch.filter(os.listdir(
input_file_des_seg = fnmatch.filter(os.listdir(
input_directory_des_seg), wildcard_pattern_des_seg)
input_directory_des_seg), wildcard_pattern_des_seg)
input_path_des_seg = os.path.join(
input_path_des_seg = os.path.join(
input_directory_des_seg, input_file_des_seg[0])
input_directory_des_seg, input_file_des_seg[0])


# define init warp
# define init warp
input_directory_warp = os.path.join(
input_directory_warp = os.path.join(
data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
wildcard_pattern_warp = 'warp_t22temp.nii.gz'
wildcard_pattern_warp = 'warp_t22temp.nii.gz'
input_file_warp = fnmatch.filter(os.listdir(
input_file_warp = fnmatch.filter(os.listdir(
input_directory_warp), wildcard_pattern_warp)
input_directory_warp), wildcard_pattern_warp)
input_path_des_warp = os.path.join(
input_path_des_warp = os.path.join(
input_directory_warp, input_file_warp[0])
input_directory_warp, input_file_warp[0])


# define output warp
# define output warp
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)
output_file = 'gm_warp_t22temp.nii.gz'
output_file = 'gm_warp_t22temp.nii.gz'
output_path_warp = os.path.join(output_directory, output_file)
output_path_warp = os.path.join(output_directory, output_file)


print(input_directory)
print(input_directory)
print(input_directory_iseg)
print(input_directory_iseg)
print(input_path_des)
print(input_path_des)
print(input_path_des_seg)
print(input_path_des_seg)
print(input_path_des_warp)
print(input_path_des_warp)
print(output_path_warp)
print(output_path_warp)


# run warping calculation
# run warping calculation
command = f"sct_register_multimodal -i {input_directory} -iseg {input_directory_iseg} -d {input_path_des} -dseg {input_path_des_seg} -initwarp {input_path_des_warp} -owarp {output_path_warp}"
print('add gm to t2-temp warp for ', subdirectory)
process = subprocess.Popen(command, shell=True)
command = [
process.wait()
'-i', input_path_des,
'-iseg', input_path_des_seg,
'-d', input_directory,
'-dseg', input_directory_iseg,
'-initwarp', input_path_des_warp,
'-param', 'step=1,type=seg,algo=rigid:step=2,type=seg,algo=bsplinesyn,slicewise=1,iter=3'
]
from spinalcordtoolbox.scripts import sct_register_multimodal
sct_register_multimodal.main(command)
# process = subprocess.Popen(command, shell=True)
# process.wait()
print('add gm to t2-temp warp completed for ', subdirectory)
print('add gm to t2-temp warp completed for ', subdirectory)


# mult epi-t2 to t2-temp
# mult epi-t2 to t2-temp
# ------------------------------------------------------------------------------------
# ------------------------------------------------------------------------------------
if epi_to_t2_mult_t2_to_temp == 1:
if epi_to_t2_mult_t2_to_temp == 1:


# Load the first warping field
# Load the first warping field
input_directory_w2 = os.path.join(
input_directory_w2 = os.path.join(
data_pre, subdirectory, "spinalcord", "registration_fmri_t2")
data_pre, subdirectory, "spinalcord", "registration_fmri_t2")
wildcard_pattern_w2 = '*mean2s30*'
wildcard_pattern_w2 = '*mean2s30*'
input_file_w2 = fnmatch.filter(os.listdir(
input_file_w2 = fnmatch.filter(os.listdir(
input_directory_w2), wildcard_pattern_w2)
input_directory_w2), wildcard_pattern_w2)
input_path_w2 = os.path.join(input_directory_w2, input_file_w2[0])
input_path_w2 = os.path.join(input_directory_w2, input_file_w2[0])
print(input_path_w2)
print(input_path_w2)


# define warping field t2 - temp
# define warping field t2 - temp
input_directory_w3 = os.path.join(
input_directory_w3 = os.path.join(
data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
data_pre, subdirectory, "spinalcord", "concat_t2_t1_temp")
wildcard_pattern_w3 = 'gm_warp_t22temp.nii.gz'
wildcard_pattern_w3 = 'gm_warp_t22temp.nii.gz'
input_file_w3 = fnmatch.filter(os.listdir(
input_file_w3 = fnmatch.filter(os.listdir(
input_directory_w3), wildcard_pattern_w3)
input_directory_w3), wildcard_pattern_w3)
input_path_w3 = os.path.join(input_directory_w3, input_file_w3[0])
input_path_w3 = os.path.join(input_directory_w3, input_file_w3[0])
print(input_path_w3)
print(input_path_w3)


# destination image hardcoded
# destination image hardcoded
des_image = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz'
des_image = '/home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz'


# output file
# output file
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "concat_epi_temp")
data_pre, subdirectory, "spinalcord", "concat_epi_temp")
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)
output_file = 'warp_epi2temp.nii.gz'
output_file = 'warp_epi2temp.nii.gz'
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)
print(output_path)
print(output_path)


# run command
# run command
command = f"sct_concat_transfo -d {des_image} -o {output_path} -w {input_path_w2} {input_path_w3}"
command = f"sct_concat_transfo -d {des_image} -o {output_path} -w {input_path_w2} {input_path_w3}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()


# apply WARPING EPI to TEMP
# apply WARPING EPI to TEMP
# ------------------------------------------------------------------------------------------
# ------------------------------------------------------------------------------------------
if final_warp_apply == 1:
if final_warp_apply == 1:


# define input EPI images
# define input EPI images
input_directory = os.path.join(
input_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "slice_time_corrected/")
data_pre, subdirectory, "spinalcord", "slice_time_corrected/")
wildcard_pattern = '*moco*'
wildcard_pattern = '*moco*'
input_file = fnmatch.filter(os.listdir(
input_file = fnmatch.filter(os.listdir(
input_directory), wildcard_pattern)
input_directory), wildcard_pattern)


for input_filename in input_file:
for input_filename in input_file:
input_path = os.path.join(input_directory, input_filename)
input_path = os.path.join(input_directory, input_filename)


# define output
# define output
output_directory = os.path.join(
output_directory = os.path.join(
data_pre, subdirectory, "spinalcord", "epi_warped_gm")
data_pre, subdirectory, "spinalcord", "epi_warped_gm")
# Check if excist, otherwise make
# Check if excist, otherwise make
if not os.path.exists(output_directory):
if not os.path.exists(output_directory):
os.makedirs(output_directory)
os.makedirs(output_directory)
output_file = input_filename.replace(".nii", "_gm_warped.nii")
output_file = input_filename.replace(".nii", "_gm_warped.nii")
output_path = os.path.join(output_directory, output_file)
output_path = os.path.join(output_directory, output_file)


# define warping field epi - temp with gm
# define warping field epi - temp with gm
input_directory_warp = os.path.join(
input_directory_warp = os.path.join(
data_pre, subdirectory, "spinalcord", "concat_epi_temp")
data_pre, subdirectory, "spinalcord", "concat_epi_temp")
wildcard_pattern_warp = 'warp_epi2temp.nii.gz'
wildcard_pattern_warp = 'warp_epi2temp.nii.gz'
input_file_warp = fnmatch.filter(os.listdir(
input_file_warp = fnmatch.filter(os.listdir(
input_directory_warp), wildcard_pattern_warp)
input_directory_warp), wildcard_pattern_warp)
input_path_warp = os.path.join(
input_path_warp = os.path.join(
input_directory_warp, input_file_warp[0])
input_directory_warp, input_file_warp[0])


# define warping field epi - t2
# define warping field epi - t2
# input_directory_w1 = os.path.join(
# input_directory_w1 = os.path.join(
# data_pre, subdirectory, "spinalcord", "registration_fmri_t2")
# data_pre, subdirectory, "spinalcord", "registration_fmri_t2")
# wildcard_pattern_w1 = '*mean2s*'
# wildcard_pattern_w1 = '*mean2s*'
# input_file_w1 = fnmatch.filter(os.listdir(
# input_file_w1 = fnmatch.filter(os.listdir(
# input_directory_w1), wildcard_pattern_w1)
# input_directory_w1), wildcard_pattern_w1)
# input_path_w1 = os.path.join(input_directory_w1, input_file_w1[0])
# input_path_w1 = os.path.join(input_directory_w1, input_file_w1[0])


# define warping field t2- t1
# define warping field t2- t1
# input_directory_w2 = os.path.join(
# input_directory_w2 = os.path.join(
# data_pre, subdirectory, "spinalcord", "register_t2_t1")
# data_pre, subdirectory, "spinalcord", "register_t2_t1")
# wildcard_pattern_w2 = '*qform2s*'
# wildcard_pattern_w2 = '*qform2s*'
# input_file_w2 = fnmatch.filter(os.listdir(
# input_file_w2 = fnmatch.filter(os.listdir(
# input_directory_w2), wildcard_pattern_w2)
# input_directory_w2), wildcard_pattern_w2)
# input_path_w2 = os.path.join(input_directory_w2, input_file_w2[0])
# input_path_w2 = os.path.join(input_directory_w2, input_file_w2[0])


# define warping field t1 - temp
# define warping field t1 - temp
# input_directory_w3 = os.path.join(
# input_directory_w3 = os.path.join(
# data_pre, subdirectory, "spinalcord", "register_t1_temp")
# data_pre, subdirectory, "spinalcord", "register_t1_temp")
# wildcard_pattern_w3 = 'warp_anat2template.nii.gz'
# wildcard_pattern_w3 = 'warp_anat2template.nii.gz'
# input_file_w3 = fnmatch.filter(os.listdir(
# input_file_w3 = fnmatch.filter(os.listdir(
# input_directory_w3), wildcard_pattern_w3)
# input_directory_w3), wildcard_pattern_w3)
# input_path_w3 = os.path.join(input_directory_w3, input_file_w3[0])
# input_path_w3 = os.path.join(input_directory_w3, input_file_w3[0])


# perform warping
# perform warping
# command = f"sct_apply_transfo -i {input_path} -d {output_path} -w {input_path_w3}"
# command = f"sct_apply_transfo -i {input_path} -d {output_path} -w {input_path_w3}"
command = f"sct_apply_transfo -i {input_path} -d /home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz -o {output_path} -w {input_path_warp}"
command = f"sct_apply_transfo -i {input_path} -d /home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz -o {output_path} -w {input_path_warp}"
process = subprocess.Popen(command, shell=True)
process = subprocess.Popen(command, shell=True)
process.wait()
process.wait()
print('Apply warping completed for ', subdirectory)
print('Apply warping completed for ', subdirectory)


# except:
# except:
# print('Error for ', subdirectory)
# print('Error for ', subdirectory)


##################################################################################
##################################################################################
# REST #
# REST #
##################################################################################
##################################################################################
# -----------------------------------------------------
# -----------------------------------------------------
# Run centerlining
# Run centerlining
# command = f"sct_deepseg_sc -i {input_path} -c t2 -qc {qc_directory} -o {output_path} -centerline svm"
# command = f"sct_deepseg_sc -i {input_path} -c t2 -qc {qc_directory} -o {output_path} -centerline svm"
# process = subprocess.Popen(command, shell=True)
# process = subprocess.Popen(command, shell=True)
# process.wait()
# process.wait()


# run segmentation
# run segmentation
# input_directory_centerline = os.path.join(
# input_directory_centerline = os.path.join(
# data_pre, subdirectory, 'spinalcord', 'segment_anat')
# data_pre, subdirectory, 'spinalcord', 'segment_anat')


# Define the wildcard pattern
# Define the wildcard pattern
# wildcard_pattern_centerline = '*centerline.nii'
# wildcard_pattern_centerline = '*centerline.nii'


# Get the list of files matching the wildcard pattern in the source directory
# Get the list of files matching the wildcard pattern in the source directory
# input_file_centerline = fnmatch.filter(os.listdir(
# input_file_centerline = fnmatch.filter(os.listdir(
# input_directory_centerline), wildcard_pattern_centerline)
# input_directory_centerline), wildcard_pattern_centerline)
# input_path_centerline = os.path.join(
# input_path_centerline = os.path.join(
# input_directory_centerline, input_file_centerline[0])
# input_directory_centerline, input_file_centerline[0])
# output_file = input_file[0].replace(".nii", "segmentation.nii")
# output_file = input_file[0].replace(".nii", "segmentation.nii")
# output_path = os.path.join(output_directory, output_file)
# output_path = os.path.join(output_directory, output_file)
# csf_directory = os.path.join(output_directory, 'csf')
# csf_directory = os.path.join(output_directory, 'csf')


# define smoothed input_path
# define smoothed input_path
# input_directory_seg = os.path.join(
# input_directory_seg = os.path.join(
# data_pre, subdirectory, 'spinalcord', "segment_anat")
# data_pre, subdirectory, 'spinalcord', "segment_anat")


# Define the wildcard pattern
# Define the wildcard pattern
# wildcard_pattern = '*smooth_segmentation.nii'
# wildcard_pattern = '*smooth_segmentation.nii'
# Get the list of files matching the wildcard pattern in the source directory
# Get the list of files matching the wildcard pattern in the source directory
# input_file_smooth = fnmatch.filter(os.listdir(
# input_file_smooth = fnmatch.filter(os.listdir(
# input_directory_seg), wildcard_pattern)
# input_directory_seg), wildcard_pattern)


# Construct input and output file paths
# Construct input and output file paths
# input_path_smooth = os.path.join(input_directory_seg, input_file_smooth[0])
# input_path_smooth = os.path.join(input_directory_seg, input_file_smooth[0])


# outputpath
# outputpath
# output_file_smooth = input_file_smooth[0].replace(
# output_file_smooth = input_file_smooth[0].replace(
# "smooth_segmentation.nii", "smooth_segmentation_twice.nii")
# "smooth_segmentation.nii", "smooth_segmentation_twice.nii")
# output_path_smooth = os.path.join(output_directory, output_file_smooth)
# output_path_smooth = os.path.join(output_directory, output_file_smooth)




# RUN SEGMENT AGAIN
# RUN SEGMENT AGAIN
# command = f"sct_propseg -i {input_path_smooth} -c t1 -init-centerline #{input_path_seg} -qc {qc_directory} -o {output_path_smooth}"
# command = f"sct_propseg -i {input_path_smooth} -c t1 -init-centerline #{input_path_seg} -qc {qc_directory} -o {output_path_smooth}"
# process = subprocess.Popen(command, shell=True)
# process = subprocess.Popen(command, shell=True)
# process.wait()
# process.wait()


# input_directory = os.path.join(
# input_directory = os.path.join(
# data_pre, subdirectory, "spinalcord", "slice_time_corrected")
# data_pre, subdirectory, "spinalcord", "slice_time_corrected")
# output_directory = os.path.join(
# output_directory = os.path.join(
# data_pre, subdirectory, "spinalcord", "gm_segment_t2")
# data_pre, subdirectory, "spinalcord", "gm_segment_t2")


# Check if the removed volumes directory exists, create it if necessary
# Check if the removed volumes directory exists, create it if necessary
# if not os.path.exists(output_directory):
# if not os.path.exists(output_directory):
# os.makedirs(output_directory)
# os.makedirs(output_directory)
# define qc
# define qc
# qc_directory = os.path.join(output_directory, 'qc')
# qc_directory = os.path.join(output_directory, 'qc')


# Get a list of all NIfTI files in the input directory
# Get a list of all NIfTI files in the input directory
# wildcard_pattern = 'slc*mean.nii'
# wildcard_pattern = 'slc*mean.nii'
# file_list = fnmatch.filter(os.listdir(
# file_list = fnmatch.filter(os.listdir(
# input_directory), wildcard_pattern) # cell array
# input_directory), wildcard_pattern) # cell array


# Perform motion correction for each input file
# Perform motion correction for each input file
# for input_file in file_list:
# for input_file in file_list:
# Construct input and output file paths
# Construct input and output file paths
# input_path = os.path.join(input_directory, input_file)
# input_path = os.path.join(input_directory, input_file)
# output_file = input_file.replace("mean.nii", "gm_segment.nii")
# output_file = input_file.replace("mean.nii", "gm_segment.nii")
# output_path = os.path.join(output_directory, output_file)
# output_path = os.path.join(output_directory, output_file)
# print(output_path)
# print(output_path)


# Run gm segment
# Run gm segment
# command = f"sct_deepseg_gm -i {input_path} -qc {qc_directory} -o {output_path} -thr 0.9999"
# command = f"sct_deepseg_gm -i {input_path} -qc {qc_directory} -o {output_path} -thr 0.9999"
# process = subprocess.Popen(command, shell=True)
# process = subprocess.Popen(command, shell=True)
# process.wait()
# process.wait()
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wednesday 24 may 2023

@author: evefra
"""

# import libraries
import os
import shutil
import subprocess
import fnmatch
import nibabel as nib
import numpy as np

# work in conda env, where SCT is installed. For me: conda activate venv_sct

######################################################################################
# INDICATE WHAT TO DO#
######################################################################################

# FUNCTIONAL
remove_volumes = 0
motion_correction = 0
# MATLAB: slice time correction
# MATLAB: mean slice time corrected image
segment_mean = 0

# ANATOMICAL T1
correct_qform = 0
first_segment_t1 = 0
smooth = 0
second_segment_t1 = 0
label_t1 = 0
register_t1_temp = 0
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wednesday 24 may 2023

@author: evefra
"""

# import libraries
import os
import shutil
import subprocess
import fnmatch
import nibabel as nib
import numpy as np

# work in conda env, where SCT is installed. For me: conda activate venv_sct

######################################################################################
# INDICATE WHAT TO DO#
######################################################################################

# FUNCTIONAL
remove_volumes = 0
motion_correction = 0
# MATLAB: slice time correction
# MATLAB: mean slice time corrected image
segment_mean = 0

# ANATOMICAL T1
correct_qform = 0
first_segment_t1 = 0
smooth = 0
second_segment_t1 = 0
label_t1 = 0
register_t1_temp = 0