Comparer le texte

Trouver la différence entre deux fichiers texte

Éditeur live

Cacher identiques

Sans retour à la ligne

Vue

Niveau de précision

Coloration syntaxique

Diffchecker Desktop La façon la plus sécurisée d'utiliser Diffchecker. Obtenez l'application Diffchecker Desktop : vos diffs ne quittent jamais votre ordinateur !Obtenir Desktop

sct_preproces.py: Before/after @sandrinebedard's changes

Créé il y a 3 ansLe diff n'expire jamais

37 suppressions

Lignes
Total
Supprimé

Caractères
Total
Supprimé

Pour continuer à utiliser cette fonctionnalité, passez à Diffchecker Pro Voir les prix

603 lignes

43 ajouts

Lignes
Total
Ajouté

Caractères
Total
Ajouté

Pour continuer à utiliser cette fonctionnalité, passez à Diffchecker Pro Voir les prix

595 lignes

#!/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

# ANATOMICAL T2*

# MATLAB: transfer from dicom to nii

first_segment_t2 = 0

# COREGISTERING

register_t2_t1 = 0

register_fmri_t2 = 0

t2_to_t1_mult_t1_to_temp = 0

correct_qform_sc = 0

gm_wm_segment = 0

gm_segment_add = 1

epi_to_t2_mult_t2_to_temp = 1

epi_to_t2_mult_t2_to_temp = 0

final_warp_apply = 1

final_warp_apply = 0

final_warp_apply_old = 1

final_warp_apply_old = 0

# LOOP OVER SUB NUMBERS

# ------------------------------------------------------------------------------

sub_list = [18, 20, 21, 24]

sub_list = [10]

for sub in sub_list:

# data path pre subject unspecific

data_pre = "/project/3023009.09/evefra"

data_pre = "./data/"

subdirectory = f"sub-{sub:03d}"

# random testing

#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_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'

# 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'

#command = f"sct_apply_transfo -i {des_t2} -d {des_pam} -w {warp}"

#process = subprocess.Popen(command, shell=True)

#process.wait()

#################################################################################

# FUNCTIONAL prep #

#################################################################################

# try:

# REMOVE FIRST TWO VOLUMES

# ------------------------------------------------------------------------------

if remove_volumes == 1:

directory = os.path.join(

data_pre, subdirectory, "spinalcord", "niftie")

removed_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "removed_volumes")

file_extension = ".nii"

# Check if the removed volumes directory exists, create it if necessary

if not os.path.exists(removed_directory):

os.makedirs(removed_directory)

# Get a list of all files in the directory with the specified extension

file_list = [f for f in os.listdir(

directory) if f.endswith(file_extension)]

# Sort the file list alphabetically

file_list.sort()

# Move the first two files to the removed volumes directory

files_to_move = file_list[:2]

# Move the files one by one

for file_name in files_to_move:

source_path = os.path.join(directory, file_name)

destination_path = os.path.join(removed_directory, file_name)

# Move the file to the destination

shutil.move(source_path, destination_path)

# Check if only two files in the removed-directory

file_list = os.listdir(removed_directory)

if len(file_list) != 2:

raise ValueError(

"The directory does not contain exactly two files.")

print('remove volumes completed for ', subdirectory)

# MOTION CORRECTION

# ----------------------------------------------------------------------------------------

if motion_correction == 1:

input_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "niftie")

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "motion_corrected")

# Check if the removed volumes directory exists, create it if necessary

if not os.path.exists(output_directory):

os.makedirs(output_directory)

# Get a list of all NIfTI files in the input directory

file_list_before = os.listdir(input_directory)

file_list = []

for names in file_list_before:

if names.endswith(".nii"):

file_list.append(names)

# Perform motion correction for each input file

for input_file in file_list:

# Construct input and output file paths

input_path = os.path.join(input_directory, input_file)

output_file = input_file.replace(".nii.gz", "_moco.nii.gz")

output_path = os.path.join(output_directory, output_file)

# Run motion correction using sct_fmri_moco

command = f"sct_fmri_moco -i {input_path} -o {output_path}"

process = subprocess.Popen(command, shell=True)

process.wait()

# print if failed

if process.returncode != 0:

print(f"Motion correction failed for {input_file}.")

print('moco completed for ', subdirectory)

# CHEKC MOTION CORRECTION AND CORRECT

# --------------------------------------------------------------------------------

# check FSLeyes

# SEGMENT MEAN slc IMAGE

# -----------------------------------------------------------------

if segment_mean == 1:

# Define the input

input_directory = os.path.join(

data_pre, subdirectory, "spinalcord", 'slice_time_corrected')

wildcard_pattern = 'spinacord_slc_mean.nii'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path_fmri = os.path.join(input_directory, input_file[0])

# define output

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "segment_slc_mean")

# Create the output directory if it doesn't exist

os.makedirs(output_directory, exist_ok=True)

# Define the quality control directory

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

# Create the output directory if it doesn't exist

os.makedirs(output_directory, exist_ok=True)

# 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 t2s -ofolder {output_directory} -qc {qc_directory} -centerline viewer' #test with t2star

command = f'sct_deepseg_sc -i {input_path_fmri} -c t2s -ofolder {output_directory} -qc {qc_directory}' #test with t2star

process = subprocess.Popen(command, shell=True)

process.wait()

print('segment mean func completed for ', subdirectory)

##############################################################################

# T2 prep #

#################################################################################

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

# ------------------------------------------------------------------------

if first_segment_t2 == 1:

# define input and output

input_directory = os.path.join(

data_pre, subdirectory, "anat_sc")

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "segment_anat_t2")

# Check if the removed volumes directory exists, create it if necessary

if not os.path.exists(output_directory):

os.makedirs(output_directory)

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

# Find the input file in the "anat" folder

import fnmatch

# Define the wildcard pattern

wildcard_pattern = 's*.nii'

wildcard_pattern = '*t2s*.nii'

# Get the list of files matching the wildcard pattern in the source directory

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

# Construct input and output file paths

input_path = os.path.join(input_directory, input_file[0])

output_file = input_file[0].replace("nii", "segmentation.nii.gz")

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 t2s -o {output_path} -qc {qc_directory} -centerline viewer" # TODO: t2 to t2star

process = subprocess.Popen(command, shell=True)

process.wait()

print('first segment t2 completed for ', subdirectory)

# CHECK SEGMENTATION AND CORRECT

# --------------------------------------------------------------------------------

# in FSLeyes

#################################################################################

# T1 prep #

#################################################################################

# CORRECT QFORM ANATOMICAL

# -----------------------------------------------------------------------------

if correct_qform == 1:

# Define the input

input_directory = os.path.join(

data_pre, subdirectory, "anat")

data_pre, subdirectory, "anat_sc")

wildcard_pattern = 's*.nii'

wildcard_pattern = '*t1*'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path = os.path.join(input_directory, input_file[0])

# define output

output_directory = os.path.join(

data_pre, subdirectory, "anat")

data_pre, subdirectory, "anat_sc")

output_file = input_file[0].replace("nii", "qform.nii")

output_path = os.path.join(output_directory, output_file)

# run command

command = f"sct_image -i {input_path} -set-sform-to-qform -o {output_path}"

process = subprocess.Popen(command, shell=True)

process.wait()

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

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

# ------------------------------------------------------------------------

if first_segment_t1 == 1:

# define input

input_directory = os.path.join(

data_pre, subdirectory, "anat")

data_pre, subdirectory, "anat_sc")

wildcard_pattern = 's*qform.nii'

wildcard_pattern = '*t1*qform.nii'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path = os.path.join(input_directory, input_file[0])

# define output

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "segment_anat")

if not os.path.exists(output_directory):

os.makedirs(output_directory)

output_file = input_file[0].replace("nii", "segmentation.nii")

output_path = os.path.join(output_directory, output_file)

# define QC directory

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

# run first segmentation

# label centerline C2 and C7

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.wait()

print('first segment t1 completed for ', subdirectory)

# SMOOTH ALONG SC TO IMPROVE SEGMENTATION

# -------------------------------------------------------------------------------

if smooth == 1:

# define input

input_directory = os.path.join(

data_pre, subdirectory, "anat")

data_pre, subdirectory, "anat_sc")

wildcard_pattern = 's*qform.nii'

wildcard_pattern = '*t1*qform.nii'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path = os.path.join(input_directory, input_file[0])

# select segmentation image

input_directory_seg = os.path.join(

data_pre, subdirectory, 'spinalcord', "segment_anat")

wildcard_pattern = '*.segmentation.nii'

input_file_seg = fnmatch.filter(os.listdir(

input_directory_seg), wildcard_pattern)

input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

# define output

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "segment_anat")

output_file_seg = input_file_seg[0].replace(

"segmentation.nii", "smooth.nii")

output_path_seg = os.path.join(output_directory, output_file_seg)

# RUN SMOOOTHING

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.wait()

print('smooth t1 completed for ', subdirectory)

# DO SECOND SEGMENTATION

# -----------------------------------------------------------------------------

if second_segment_t1 == 1:

# define input

input_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "segment_anat")

wildcard_pattern = '*smooth.nii'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path = os.path.join(input_directory, input_file[0])

# define output

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "segment_anat_second")

if not os.path.exists(output_directory):

os.makedirs(output_directory)

output_file = input_file[0].replace(

"smooth.nii", "second_segmentation.nii")

output_path = os.path.join(output_directory, output_file)

# define qcdirectory

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

# select segmentation image

input_directory_seg = os.path.join(

data_pre, subdirectory, 'spinalcord', "segment_anat")

wildcard_pattern = '*.segmentation.nii'

input_file_seg = fnmatch.filter(os.listdir(

input_directory_seg), wildcard_pattern)

input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

# 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}"

process = subprocess.Popen(command, shell=True)

process.wait()

print('second segment t1 completed for ', subdirectory)

# CHECK SEGMENTATION AND CORRECT

# --------------------------------------------------------------------------------

# in FSLeyes

# LABEL T1

# -------------------------------------------------------------------------------

if label_t1 == 1:

# define input anatomical image

input_directory = os.path.join(

data_pre, subdirectory, "anat")

data_pre, subdirectory, "anat_sc")

wildcard_pattern = 's3*qform.nii'

wildcard_pattern = '*t1*qform.nii'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path = os.path.join(input_directory, input_file[0])

# define input segmentation

input_directory_seg = os.path.join(

data_pre, subdirectory, 'spinalcord', "segment_anat_second")

wildcard_pattern = '*segmentation.nii'

input_file_seg = fnmatch.filter(os.listdir(

input_directory_seg), wildcard_pattern)

input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

# define output paths and name

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "label_anat")

if not os.path.exists(output_directory):

os.makedirs(output_directory)

output_file_manual = input_file[0].replace(".nii", "manlabel.nii")

output_path_manual = os.path.join(output_directory, output_file_manual)

output_file = input_file[0].replace(".nii", "label.nii")

output_path = os.path.join(output_directory, output_file)

# define qc directory

qc_directory_man = os.path.join(output_directory, 'qc_manual')

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

# 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 '"

process = subprocess.Popen(command, shell=True)

process.wait()

# define input man label

# wildcard_pattern = '*manlabel.nii'

# input_file_man = fnmatch.filter(os.listdir(

# output_directory), wildcard_pattern)

# input_path_man = os.path.join(output_directory, input_file_man[0])

# correcting s or qform if needed

# sct_image -set-sform-to-qform -i {

# 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}"

# process = subprocess.Popen(command, shell=True)

# process.wait()

print('label t1 completed for ', subdirectory)

# CHECK LABELLING

# --------------------------------------------------------------------------------

# manual label more vertebra

# 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.'"

# process = subprocess.Popen(command, shell=True)

# process.wait()

# REGISTER FROM T1 to TEMPLATE

# -----------------------------------------------------------------------------

if register_t1_temp == 1:

# define input image (T1)

input_directory = os.path.join(

data_pre, subdirectory, "anat")

data_pre, subdirectory, "anat_sc")

wildcard_pattern = 's*qform.nii'

wildcard_pattern = '*t1*qform.nii'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path = os.path.join(input_directory, input_file[0])

# define input segmentation

input_directory_seg = os.path.join(

data_pre, subdirectory, 'spinalcord', "segment_anat_second")

wildcard_pattern = '*segmentation.nii'

input_file_seg = fnmatch.filter(os.listdir(

input_directory_seg), wildcard_pattern)

input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

# define labels input image

input_directory_lab = os.path.join(

data_pre, subdirectory, "spinalcord", "label_anat")

wildcard_pattern = 's*qformmanlabel.nii' # name if not manual labelled?

input_file_lab = fnmatch.filter(os.listdir(

input_directory_lab), wildcard_pattern)

input_path_lab = os.path.join(input_directory_lab, input_file_lab[0])

# define output directory

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "register_t1_temp")

# define qc

qc_directory = os.path.join(output_directory, 'qc_t1_temp')

# 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} -ldisc {input_path_lab} -c t1 -qc {qc_directory} -ofolder {output_directory}"

process = subprocess.Popen(command, shell=True)

process.wait()

print('register calc t1 - temp completed for ', subdirectory)

#################################################################################

# COREGISTER #

#################################################################################

# REGISTER FROM T2 to T1

# -----------------------------------------------------------------------------

if register_t2_t1 == 1:

# define input image (T2)

input_directory = os.path.join(

data_pre, subdirectory, "anat_sc")

wildcard_pattern = 's3*01.nii'

wildcard_pattern = '*t2s*'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path = os.path.join(input_directory, input_file[0])

# define input segmentation of T2

input_directory_seg = os.path.join(

data_pre, subdirectory, 'spinalcord', "segment_anat_t2")

wildcard_pattern = '*segmentation.nii'

wildcard_pattern = '*segmentation*'

input_file_seg = fnmatch.filter(os.listdir(

input_directory_seg), wildcard_pattern)

input_path_seg = os.path.join(input_directory_seg, input_file_seg[0])

# define input segmentation of T1

input_directory_seg = os.path.join(

data_pre, subdirectory, 'spinalcord', "segment_anat_second")

wildcard_pattern = '*segmentation.nii'

input_file_seg = fnmatch.filter(os.listdir(

input_directory_seg), wildcard_pattern)

input_path_seg_t1 = os.path.join(

input_directory_seg, input_file_seg[0])

# define labels input image (from t1?)

input_directory_lab = os.path.join(

data_pre, subdirectory, "spinalcord", "label_anat")

wildcard_pattern = 's*manlabel.nii' # name if not manual labelled?

input_file_lab = fnmatch.filter(os.listdir(

input_directory_lab), wildcard_pattern)

input_path_lab = os.path.join(input_directory_lab, input_file_lab[0])

# define to-be-warped image (T1)

input_directory = os.path.join(

data_pre, subdirectory, "anat")

data_pre, subdirectory, "anat_sc")

wildcard_pattern = 's*qform.nii'

wildcard_pattern = '*t1*qform.nii'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path_warped = os.path.join(input_directory, input_file[0])

# define output directory

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "register_t2_t1")

# define qc

qc_directory = os.path.join(output_directory, 'qc_t1_temp')

# 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}"

command = f"sct_register_multimodal -i {input_path} -d {input_path_warped} -iseg {input_path_seg} -dseg {input_path_seg_t1} -param step=1,type=seg,algo=centermass -qc {qc_directory} -ofolder {output_directory}"

process = subprocess.Popen(command, shell=True)

process.wait()

print('Register calc t2 - t1 completed for ', subdirectory)

# CALCULATE warping field from fmri to T2*

# ---------------------------------------------------------------------------------

if register_fmri_t2 == 1:

# SOURCE define input image (slc fmri mean)

input_directory = os.path.join(

data_pre, subdirectory, "spinalcord", 'slice_time_corrected')

wildcard_pattern = 'spinacord_slc_mean.nii'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path_fmri = os.path.join(input_directory, input_file[0])

# Source SEGMENT define input segmentation of fmri mean

input_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "segment_slc_mean")

wildcard_pattern = 'spinacord_slc_mean_seg.nii'

input_file = fnmatch.filter(os.listdir(

input_directory), wildcard_pattern)

input_path_fmri_seg = os.path.join(input_directory, input_file[0])

# DESTINATION define input of T2

input_directory_d = os.path.join(

data_pre, subdirectory, "anat_sc")

wildcard_pattern_d = 's3*01.nii'

wildcard_pattern_d = '*t2s*qform.nii'

input_file_d = fnmatch.filter(os.listdir(

input_directory_d), wildcard_pattern_d)

input_path_des = os.path.join(input_directory_d, input_file_d[0])

# DESTINATION SEGMENT define input segmentation of T2

input_directory_seg = os.path.join(

data_pre, subdirectory, 'spinalcord', "segment_anat_t2")

wildcard_pattern = '*segmentation.nii'

wildcard_pattern = '*segmentation*'

input_file_seg = fnmatch.filter(os.listdir(

input_directory_seg), wildcard_pattern)

input_path_des_seg = os.path.join(

input_directory_seg, input_file_seg[0])

# define output directory

output_directory = os.path.join(

data_pre, subdirectory, "spinalcord", "registration_fmri_t2")

# define qc

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

# 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}"

# run coregistering FMRI --> T2

command = f"sct_register_multimodal -i {input_path_fmri} -d {input_path_des} -iseg {input_path_fmri_seg} -dseg {input_path_des_seg} -param step=1,type=seg,algo=centermass -qc {qc_directory} -ofolder {output_directory}"

process = subprocess.Popen(command, shell=True)

process.wait()

print('Register calc fmri - t2 completed for ', subdirectory)

# mult t2-t1 to t1-temp

# ------------------------------------------------------------------------------------

if t2_to_t1_mult_t1_to_temp == 1:

# Load the first warping field

# Input T2

input_directory_w2 = os.path.join(

input_directory_t2 = os.path.joi

data_pre, subdirectory, "spinalcord", "register_t2_t1")

wildcard_pattern_w2 = '*qform2s*'

input_file_w2 = fnmatch.filter(os.listdir(

input_directory_w2), wildcard_pattern_w2)

input_path_w2 = os.path.join(input_directory_w2, input_file_w2[0])

# define warping field t1 - temp

input_directory_w3 = os.path.join(

data_pre, subdirectory, "spinalcord", "register_t1_temp")

wildcard_pattern_w3 = 'warp_anat2template.nii.gz'

input_file_w3 = fnmatch.filter(os.listdir(

input_directory_w3), wildcard_pattern_w3)

input_path_w3 = os.path.join(input_directory_

Différences enregistrées

Texte d'origine

Ouvrir un fichier

#!/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

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

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

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

sub_list = [18, 20, 21, 24]
for sub in sub_list:

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

# random testing

#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_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'
    # 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'

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

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

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

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

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

# Sort the file list alphabetically
        file_list.sort()

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

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

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

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

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

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

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

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

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

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

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

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

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

# check FSLeyes

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

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

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

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

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

# segmentation
        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.wait()
        print('segment mean func completed for ', subdirectory)

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

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

if first_segment_t2 == 1:

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

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

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

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

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

# Construct input and output file paths
        input_path = os.path.join(input_directory, input_file[0])
        output_file = input_file[0].replace("nii", "segmentation.nii.gz")
        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"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('first segment t2 completed for ', subdirectory)

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

# in FSLeyes

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

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

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

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

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

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

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

if first_segment_t1 == 1:

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

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

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

# run first segmentation
        # label centerline C2 and C7
        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.wait()
        print('first segment t1 completed for ', subdirectory)

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

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

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

# RUN SMOOOTHING
        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.wait()
        print('smooth t1 completed for ', subdirectory)

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

if second_segment_t1 == 1:

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

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

# 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}"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('second segment t1 completed for ', subdirectory)

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

# in FSLeyes

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

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

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

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

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

# 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 '"
        process = subprocess.Popen(command, shell=True)
        process.wait()

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

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

# 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}"
        # process = subprocess.Popen(command, shell=True)
        # process.wait()
        print('label t1 completed for ', subdirectory)

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

# manual label more vertebra
    # 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.'"
    # process = subprocess.Popen(command, shell=True)
    # process.wait()

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

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

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

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

# 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}"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('register calc t1 - temp completed for ', subdirectory)

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

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

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

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

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

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

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

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

# 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}"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('Register calc t2 - t1 completed for ', subdirectory)

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

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

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

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

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

# 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}"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('Register calc fmri - t2 completed for ', subdirectory)

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

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

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

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

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

# run command
        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.wait()

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

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

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

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

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

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

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

if gm_wm_segment == 1:

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

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

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

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

# Run gm segment
        from spinalcordtoolbox.scripts import sct_deepseg_gm
        sct_deepseg_gm.main(['-i', input_path_des, '-qc', qc_directory, '-o', output_path])
        #process = subprocess.Popen(command, shell=True)
        #process.wait()

# segment wm
        # input segment
        input_directory_seg = os.path.join(
            data_pre, subdirectory, 'spinalcord', "segment_anat_t2")
        wildcard_pattern = '*segmentation.nii'
        input_file_seg = fnmatch.filter(os.listdir(
            input_directory_seg), wildcard_pattern)
        input_path_des_seg = os.path.join(
            input_directory_seg, input_file_seg[0])
        
        # segmentation gm
        input_directory_gm = os.path.join(
            data_pre, subdirectory, "spinalcord", "gm_segment_t2_anat")
        wildcard_pattern_gm = 'gm*.nii.gz'
        input_file_gm = fnmatch.filter(os.listdir(
            input_directory_gm), wildcard_pattern_gm)
        input_path_gm = os.path.join(
            input_directory_gm, input_file_gm[0])

# output name
        output_directory = os.path.join(
            data_pre, subdirectory, "spinalcord", "wm_segment_t2_anat")
        output_file_seg = 'wm_t2.nii.gz'
        output_path = os.path.join(output_directory, output_file_seg)

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

# run
        command = f"sct_maths -i {input_path_des_seg} -sub {input_path_gm} -o {output_path}"
        process = subprocess.Popen(command, shell=True)
        process.wait()
        print('segment gm and wm completed for ', subdirectory)
        
        
        
        
    # ADD SEGMENTATION GM
    # ----------------------------------------------------------------------------
    if gm_segment_add == 1:

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

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

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

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

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

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

# 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}"
        #process = subprocess.Popen(command, shell=True)
        #process.wait()
        
        print('add gm to t2-temp warp for ', subdirectory)
        command = [
            '-i', input_path_des,
            '-iseg', input_path_des_seg,
            '-d', input_directory,
            '-dseg', input_directory_iseg,
            '-initwarp', input_path_des_warp,
            '-owarp', output_path_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)

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

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

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

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

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

# run command
        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.wait()

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

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

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

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

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

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

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

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

# 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 /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.wait()
            print('Apply warping completed for ', subdirectory)
            
# Simpel approach warping
#----------------------------------------------------------------------------
if final_warp_apply_old == 1:

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

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

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

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

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

# 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 /home/affneu/evefra/sct_5.8/data/PAM50/template/PAM50_t1.nii.gz -o {output_path} -w {input_path_w1} {input_path_w2} {input_path_w3}"
            process = subprocess.Popen(command, shell=True)
            process.wait()
            print('Apply warping completed for ', subdirectory)

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

##################################################################################
# REST #
##################################################################################
# -----------------------------------------------------
# Run centerlining
# 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.wait()

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

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

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

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

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

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

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

# 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}"
# process = subprocess.Popen(command, shell=True)
# process.wait()

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

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

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

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

# Run gm segment
# command = f"sct_deepseg_gm -i {input_path} -qc {qc_directory} -o {output_path} -thr 0.9999"
# process = subprocess.Popen(command, shell=True)
# process.wait()

Texte modifié

Ouvrir un fichier

#!/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

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

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

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

sub_list = [10]
for sub in sub_list:

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