<reponame>noxtoby/anatomically_constrained_tractography
#!/usr/bin/env python3
# 
# Network spreading models using personalised (age/gender-matched) connectomes 
# from PPMI healthy controls for anatomically constrained tractography
#
# ACT prep (MRTrix3 pipeline)
#
# Prerequisite: walk through act_ppmi_0_prep.py
#
# <NAME>, February 2020

import os,subprocess,glob
import pandas as pd

#* Name of CSV file containing image file paths and IDs: generated by act_ppmi_0_prep.py
csv_input_file = 'act_ppmi_filepaths.csv'

#* File locations
top_folder = '~/ox_tools/anatomically_constrained_tractography'
data_folder = os.path.join(top_folder,'PPMI')
act_folder  = os.path.normpath(os.path.join(data_folder,'act'))
dMRI_folder = os.path.normpath(os.path.join(act_folder,'dMRI'))
sMRI_folder = os.path.normpath(os.path.join(act_folder,'sMRI'))
json_folder = os.path.normpath(os.path.join(act_folder,'json'))
csv_input_path = os.path.join(top_folder,csv_input_file)
act_csv_path = os.path.join(act_folder,'act_ppmi.csv')

#* Preprocessing functions copied from my connectome pipeline
wd = os.path.join(top_folder,'scripts')
os.chdir(wd)
from act_utilities import *

#* Preprocessing output folders
denoise_tail = '_denoised'
preproc_tail = '_preproc'
biascorrect_tail = '_biascorrect'
dMRI_tails = [denoise_tail,preproc_tail,biascorrect_tail]
dMRI_preproc_folders = [os.path.join(act_folder,'dMRI%s' % d) for d in dMRI_tails]
reg_tail = '_reg'
sMRI_reg_folder      =  os.path.join(act_folder,'sMRI%s' % reg_tail)
make_sure_path_exists(sMRI_reg_folder)
for f in dMRI_preproc_folders:
    make_sure_path_exists(f)

os.chdir(act_folder)


#####################################################################################################
############ Prep: check for dependencies, convert dMRI to MIF, prepare output DataFrame ############
#* Check for dependencies: ANTS, etc.
depCheckVal = checkForDependencies()
if depCheckVal!=0:
    print('ERROR: missing dependency {0}'.format(depCheckVal))
    print('see checkForDependencies() function for details')
    sys.exit(depCheckVal)
#* DataFrame for ACT files
if os.path.exists(act_csv_path):
    df = pd.read_csv(act_csv_path)
else:
    df_act = pd.read_csv(csv_input_path)
    df_act_sMRI = df_act.loc[df_act["Modality"]=='sMRI'].sort_values(by='Subject ID').reset_index(drop=True)
    df_act_dMRI = df_act.loc[df_act["Modality"]=='dMRI'].sort_values(by='Subject ID').reset_index(drop=True)
    col_preproc = "Path preprocessed"
    #* Match sMRI to dMRI using Subject ID
    cols = ["Subject ID","Path"]
    df_act_sMRI = df_act_sMRI[cols].rename(columns={"Path": "Path_sMRI"})
    df_act_dMRI = df_act_dMRI[cols].rename(columns={"Path": "Path_dMRI"})
    df = pd.merge(df_act_sMRI,df_act_dMRI)
    df["Preprocessed"] = False
#* dMRI: convert NII to MIF, with embedded b-values and b-vectors (and other metadata from json, including Phase Encoding Direction)
for i in df.index.values:
    sMRI = df.at[i,'Path_sMRI']
    dMRI = df.at[i,'Path_dMRI']
    dMRI_mif = dMRI.replace('.nii.gz','.mif')
    json_file = os.path.join(json_folder,os.path.basename(dMRI).replace('.nii.gz','.json'))
    #* Convert to MRtrix image format (MIF), with embedded b-values and b-vectors
    if not os.path.isfile(dMRI_mif):
        bvec = dMRI.replace('.nii.gz','.bvec')
        bval = dMRI.replace('.nii.gz','.bval')
        mrconvert_check = subprocess.run(['mrconvert','-json_import',json_file,'-fslgrad',bvec,bval,dMRI,dMRI_mif])
        #rm_check = subprocess.run(['rm',bvec,bval,dMRI])
    df.at[i,"Path_dMRI"] = dMRI_mif
#* Write ACT preprocessing status and results (file locations) to CSV
df.to_csv(act_csv_path,index=False)
#>>>>>>>>>>>> I fixed a missing file. Keeping the code here for future reference
# #* One DWI file was missing bvals/bvecs files, so I manually downloaded an alternative (same participant, earlier visit)
# patno = "3855"
# i = df.loc[df["Subject ID"]==3855].index.values[0]
# #* FIXME: Remove existing files (did it manually)
# # o = subprocess.run(['rm',os.path.join(sMRI_folder,patno+"*")])
# # o = subprocess.run(['rm',os.path.join(dMRI_folder,patno+"*")])
# zipfile = "/Users/noxtoby/Documents/EuroPOND-WP8-NetworkModels/tractography/PPMI/zips/PPMI-for-CJD-missing_grad.zip"
# o = subprocess.run(['open',zipfile])
# new_dicom_dMRI = "/Users/noxtoby/Documents/EuroPOND-WP8-NetworkModels/tractography/PPMI/zips/PPMI/3855/DTI_gated/2012-07-24_08_47_28.0/S204392"
# new_dicom_sMRI = "/Users/noxtoby/Documents/EuroPOND-WP8-NetworkModels/tractography/PPMI/zips/PPMI/3855/MPRAGE_GRAPPA/2012-07-24_08_38_32.0/S204393"
# nifti_folder = os.path.join(data_folder,'nifti')
# o = subprocess.run(['dcm2niix','-o',nifti_folder,'-f','%i_%t_%s_%p','-z','y',new_dicom_sMRI])
# o = subprocess.run(['dcm2niix','-o',nifti_folder,'-f','%i_%t_%s_%p','-z','y',new_dicom_dMRI])
# #* Convert dMRI to MIF, move to dMRI folder, delete old file
# old_dMRI = df.at[i,'Path_dMRI']
# new_dMRI_nifti = os.path.join(nifti_folder,"3855_20120724083832_2_DTI_gated.nii.gz")
# bval = new_dMRI_nifti.replace('.nii.gz','.bval')
# bvec = new_dMRI_nifti.replace('.nii.gz','.bvec')
# json_file = new_dMRI_nifti.replace('.nii.gz','.json')
# new_dMRI_mif = os.path.join(dMRI_folder,os.path.basename(new_dMRI_nifti)).replace('.nii.gz','.mif')
# mrconvert_check = subprocess.run(['mrconvert','-json_import',json_file,'-fslgrad',bvec,bval,new_dMRI_nifti,new_dMRI_mif])
# #* Move sMRI
# old_sMRI = df.at[i,'Path_sMRI']
# new_sMRI = os.path.join(nifti_folder,"3855_20120724083832_9_MPRAGE_GRAPPA.nii.gz")
# new_sMRI_loc = os.path.join(sMRI_folder,os.path.basename(new_sMRI))
# o = subprocess.run(['mv',new_sMRI,new_sMRI_loc])
# #* Fix paths for the new files
# df.at[i,"Path_dMRI"] = new_dMRI_mif
# df.at[i,"Path_sMRI"] = new_sMRI_loc
# df["Path_sMRI"] = df["Path_sMRI"].map(lambda x: x.split('.')[0] + '.nii.gz')
#<<<<<<<<<<<<
#####################################################################################################




#####################################################################################################
############ Preprocessing ############
#* Steps: 1. preprocess_dMRI => bias corrected dMRI
#*        2. register sMRI to dMRI_biascorrected
#* Cluster:
#*        3. run GIF3 on sMRI_registered
#*        4. ACT stuff: 5TT, nodes / WM_FODs (response functions & CSD)
#*        5. mtnormalise
#*        6. Tractography
#*        7. Connectomes
preproc_bool = False

# It might be prudent to use the cluster for this (DWI bias correction takes ages!)
#   => a) copy dMRI and sMRI folders onto cluster
#      b) run act_utilities.preprocess_individual() OR a version of act_utilities.preprocess_dMRI() and .preprocess_sMRI()

new_column_dMRI = 'Path_dMRI_biascorrected'
new_column_sMRI = 'Path_sMRI_registered'

#****** Local batch preprocessing using multicore CPU
n_cpus = 4
#* mrtrix3: dwidenoise => dwipreproc (-rpe_header)
batch_cmd = [
    'foreach', '-%i' % n_cpus,
    'dMRI/*.mif',':',
    'dwidenoise','IN','dMRI_{0}/PRE_{0}.mif'.format(denoise_tail),
    '"&&"',
    'dwipreproc','-rpe_header','dMRI_{0}/PRE_{0}.mif'.format(denoise_tail),'dMRI_{0}/PRE_{0}.mif'.format(preproc_tail)
]
if preproc_bool:
    foreach_check = os.system(' '.join(batch_cmd)) # Subprocess doesn't work! subprocess.run(batch_cmd) FIXME???
#* mrtrix3: dwibiascorrect
batch_cmd = [
    'foreach','-%i' % n_cpus,
    'dMRI%s/*.mif' % preproc_tail,':',
    'dwibiascorrect','-ants','IN','dMRI{0}/PRE{0}.mif'.format(biascorrect_tail)
]
if preproc_bool:
    foreach_check = os.system(' '.join(batch_cmd))

#* niftyreg: register sMRI to biascorrected dMRI
indices = df.index.values
cnt = 0
for i in indices:
    cnt += 1
    if not df.at[i,"Preprocessed"]:
        sMRI = df.at[i,'Path_sMRI']
        dMRI = df.at[i,'Path_dMRI']
        dMRI_denoised    = dMRI.replace('dMRI','dMRI%s' % denoise_tail).replace('.mif','%s.mif' % denoise_tail)
        dMRI_preproc     = dMRI.replace('dMRI','dMRI%s' % preproc_tail).replace('.mif','%s.mif' % preproc_tail)
        dMRI_biascorrect = dMRI.replace('dMRI','dMRI%s' % biascorrect_tail).replace('.mif','%s%s.mif' % (preproc_tail,biascorrect_tail))
        sMRI_reg = sMRI.replace('sMRI','sMRI%s' % reg_tail)
        print('************\n%i of %i - Preprocessing %s and %s ' % (cnt,len(indices),os.path.basename(dMRI),os.path.basename(sMRI)))
        #* Serial preprocessing: (dMRI_denoised,dMRI_preproc,dMRI_biascorrected,sMRI_reg,reg_aladin_output_transform) = preprocess_individual(dMRI_mif,sMRI)
        
        #* register sMRI to biascorrected dMRI
        if os.path.exists(dMRI_biascorrect):
            (reg_aladin_check,sMRI_reg,reg_aladin_output_transform) = reg_sMRI_to_dMRI(sMRI,dMRI_biascorrect)
            #*** Save to DataFrame
            df.at[i,new_column_dMRI] = dMRI_biascorrect
            df.at[i,new_column_sMRI] = sMRI_reg
            df.at[i,"Preprocessed"] = True
        else:
            print("Could not find bias-corrected file: %s " % dMRI_biascorrect)
            print(" >>> registration not performed")
            df.at[i,"Preprocessed"] = False
        #* Update CSV
        df.to_csv(act_csv_path,index=False)
    else:
        print("Files already preprocessed: %s, %s " % (os.path.basename(dMRI),os.path.basename(sMRI)))

#*** 5. run GIF3 / FreeSurfer on sMRI_registered
df_ = pd.read_csv(act_csv_path)
print("\nChecking whether you've run GIF on the sMRI (you'll have to change this for FreeSurfer, or remove: see below)\n")
algorithm = 'gif3'
sMRI_parc_txt = 'sMRI_%s' % algorithm
gif_parc_tag = "_NeuroMorph_Parcellation.nii.gz"
gif_folder = os.path.join(act_folder,sMRI_parc_txt)
gif3_parcellations_path = glob.glob(os.path.join(gif_folder,"*/*%s" % gif_parc_tag))
gif3_parcellations = [os.path.basename(g) for g in gif3_parcellations_path]
df_gif3 = pd.DataFrame(columns=[sMRI_parc_txt],data=gif3_parcellations)
df_gif3['Subject ID'] = df_gif3[sMRI_parc_txt].map(lambda x: x.split('_')[0]).astype(int)
df = df_.merge(df_gif3,on='Subject ID')
print(" > > > Found %i GIF3 parcellations out of %i" % (len(gif3_parcellations),df.shape[0]))

print("\nChecking whether you've run FreeSurfer on the sMRI\n")
algorithm = 'freesurfer'
sMRI_parc_txt = 'sMRI_%s' % algorithm
freesurfer_parc = "aparc+aseg.mgz"
freesurfer_folder = os.path.join(act_folder,sMRI_parc_txt)
freesurfer_parcellations_path = glob.glob(os.path.join(freesurfer_folder,"*/mri/%s" % freesurfer_parc))
#freesurfer_parcellations = freesurfer_parcellations_path # Could rename the files (to approximately match GIF filenames?)
df_freesurfer = pd.DataFrame(columns=[sMRI_parc_txt],data=freesurfer_parcellations_path)
df_freesurfer['Subject ID'] = df_freesurfer[sMRI_parc_txt].map(lambda x: x.split('/')[-3].split('_')[0]).astype(int)
df = df.merge(df_freesurfer,on='Subject ID')
print(" > > > Found %i FreeSurfer parcellations out of %i" % (len(freesurfer_parcellations_path),df.shape[0]))


#*** 5.2 DWI masks
print(">>> DWI masks for dwi2fod")
masks_folder = os.path.join(act_folder,'dMRI_mask')
oot = os.system('mkdir -p %s' % masks_folder)
#* mrtrix3: dwi2mask
mask_tail = '_mask'
batch_cmd = [
    'foreach','-%i' % n_cpus,
    'dMRI%s/*.mif' % biascorrect_tail,':',
    'dwi2mask','IN','dMRI{0}/PRE{0}.mif'.format(mask_tail)
]
if preproc_bool:
    foreach_check = os.system(' '.join(batch_cmd))

#*** 6. ACT stuff: 5TT, nodes / WM_FODs (response functions & CSD)
algorithm = 'freesurfer'
sMRI_parc_txt = 'sMRI_%s' % algorithm
bizzi_lobes = False
if bizzi_lobes:
    algorithm_ = algorithm + "_bizzi2lobes"
else:
    algorithm_ = algorithm
sMRI_5tt_txt = 'Path_5tt_%s' % algorithm_
sMRI_nodes_txt = 'Path_nodes_%s' % algorithm_
#* labelconvert: use bizzi2lobes_cingsep_convert.txt (c.f. MRTRIX3:path/to/labelconvertLUTs/f2lobes....txt)
print(">>> ACT prep: 5TT, nodes, CSD (WM_FODs and response functions)")
print("    *** Chosen algorithm: %s (%s)" % (algorithm,algorithm_))
fivett_folder = os.path.join(act_folder,'sMRI_5tt')
oot = os.system('mkdir -p %s' % fivett_folder)
fttrc = -1
# df.drop(columns=[sMRI_nodes_txt],inplace=True)
for i in df.index.values:
    participant_id = df.at[i,'Subject ID']
    sMRI_reg = os.path.join(act_folder,'sMRI_reg',os.path.basename(df.at[i,'Path_sMRI_registered']))
    if algorithm == 'gif3':
        sMRI_gif = os.path.join(act_folder,sMRI_parc_txt,os.path.basename(df.at[i,sMRI_parc_txt]).replace(gif_parc_tag,'_share_apps_cmic_GIF_bin'),os.path.basename(df.at[i,sMRI_parc_txt]))
        sMRI_parc = sMRI_gif
        sMRI_5TT = os.path.join(fivett_folder,os.path.basename(sMRI_gif).replace(gif_parc_tag,'_5tt_%s.mif' % algorithm))
    elif algorithm == 'freesurfer':
        sMRI_freesurfer = df.at[i,sMRI_parc_txt]
        sMRI_parc = sMRI_freesurfer
        sMRI_5TT = os.path.join(fivett_folder,os.path.basename(sMRI_reg).replace('.nii.gz','_5tt_%s.mif' % algorithm))
    if not os.path.isfile(sMRI_5TT):
        fiveTT_check = subprocess.run(['5ttgen',algorithm,sMRI_parc,sMRI_5TT])
        fttrc = fiveTT_check.returncode
    else:
        print('5TT file found')
        fttrc = 0
    if fttrc == 0:
        df.at[i,sMRI_5tt_txt] = sMRI_5TT
        print('\t> > > \n\t> > > 5ttgen done. Proceeding with labelconvert\n\t> > >')
        (labelconvert_check,nodes) = convert_sMRI_labels(sMRI_parc,algorithm=algorithm,cluster=False,bizzi=bizzi_lobes)
        if labelconvert_check==0:
            df.at[i,sMRI_nodes_txt] = nodes
        else:
            df.at[i,sMRI_nodes_txt] = np.nan
    else:
        df.at[i,sMRI_5tt_txt] = np.nan
#* Move nodes files
nodes_folder = 'sMRI_nodes'
oot = os.system('mkdir -p %s' % os.path.join(act_folder,nodes_folder))
for i in df.index.values:
    nodes = df.at[i,sMRI_nodes_txt]
    nodes_new = nodes.replace("/"+ sMRI_parc_txt +"/","/"+ nodes_folder +"/")
    oot = subprocess.run(['mv',nodes,nodes_new])
    if oot.returncode==0:
        df.at[i,"Path_nodes"] = nodes_new

#* WMFODs: instead could do on the cluster via ppmi_tckgen_cluster
rf_algo = 'tournier' #'msmt_5tt','dhollander'
fods_folder = os.path.join(act_folder,'dMRI_fods'); os.system('mkdir -p %s' % fods_folder)
#* Re-run these after fixing the masks
indices_to_rerun = df.index[df['Subject ID'].isin([3114, 3765, 3767, 3768, 3857])]

for i in df.index.tolist(): #indices_to_rerun: #
    sid = df.at[i,'Subject ID']
    dwi = df.at[i,'Path_dMRI_biascorrected']
    dwi_mask = dwi.replace('/dMRI_biascorrect/','/dMRI_mask/').replace('.mif','_mask.mif')
    fivett = df.at[i,sMRI_5tt_txt]
    wm_response_txt  = os.path.join(fods_folder,'%s_wm_response_%s.txt' % (sid,rf_algo))
    gm_response_txt  = os.path.join(fods_folder,'%s_gm_response_%s.txt' % (sid,rf_algo))
    csf_response_txt = os.path.join(fods_folder,'%s_csf_response_%s.txt' % (sid,rf_algo))
    voxels_mif       = os.path.join(fods_folder,'%s_voxels_%s.mif' % (sid,rf_algo))
    #* dwi2response
    # > PPMI has only two b-values, which precludes multi-shell multi-tissue CSD, 
    # > so use 'dwi2response tournier' rather than msmt_5tt or dhollander 
    # > (which will output matrices rather than vectors in the RF txt file)
    #cmd = ['dwi2response','dhollander','-shells','0,1000',dwi,wm_response_txt,gm_response_txt,csf_response_txt,'-voxels',voxels_mif]
    #cmd = ['dwi2response','msmt_5tt',dwi,fivett,wm_response_txt,gm_response_txt,csf_response_txt,'-voxels',voxels_mif]
    cmd = ['dwi2response',rf_algo,dwi,wm_response_txt]
    if not os.path.exists(wm_response_txt):
        oot = subprocess.run(cmd)
        orc = oot.returncode
    else:
        print(' Subject %s: %s exists' % (sid,wm_response_txt))
        orc = 0
    if orc==0:
        print(' Subject %s: dwi2response completed' % sid)
    else:
        print(' Subject %s: dwi2response failed (return code %i)' % (sid,orc))
    #* dwi2fod
    wm_fods_mif = os.path.join(fods_folder,'%s_wm_fods_%s.mif' % (sid,rf_algo))
    gm_fods_mif = os.path.join(fods_folder,'%s_gm_fods_%s.mif' % (sid,rf_algo))
    csf_fods_mif = os.path.join(fods_folder,'%s_csf_fods_%s.mif' % (sid,rf_algo))
    cmd = ['dwi2fod','csd','-mask',dwi_mask,dwi,wm_response_txt,wm_fods_mif] 
    #cmd = ['dwi2fod','csd',dwi,gm_response_txt,gm_fods_mif,csf_response_txt,csf_fods_mif]
    if not os.path.exists(wm_fods_mif):
        oot = subprocess.run(cmd)
        #os.system(' '.join(cmd))
        orc = oot.returncode
    else:
        print(' Subject %s: %s exists' % (sid,wm_fods_mif))
    if orc==0:
        print(' Subject %s: dwi2fod completed' % sid)
        df.at[i,"Path_wm_fod"] = wm_fods_mif
        #df.at[i,"Path_gm_fod"] = gm_fods_mif
        #df.at[i,"Path_csf_fod"] = csf_fods_mif
    else:
        print(' Subject %s: dwi2fod failed (return code %i)' % (sid,orc))

#* Check the response functions:
# shview wm_response.txt

# #* MSMT-only: Check the FODs: convert to RGB image
# for i in df.index.tolist():
#     sid = df.at[i,'Subject ID']
#     wm_fods_mif = os.path.join(fods_folder,'%s_wm_fods_%s.mif' % (sid,rf_algo))
#     tissue_rgb_mif = wm_fods_mif.replace('wm_fods','tissueRGB')
#     cmd = 'mrconvert %s - -coord 3 0 | mrcat CSF.mif GM.mif - %s -axis 3' % (wm_fods_mif,tissue_rgb_mif)
#     if os.path.exists(wm_fods_mif):
#         df.at[i,"Path_wm_fod"] = wm_fods_mif

# mrconvert WM_FODs.mif - -coord 3 0 | mrcat CSF.mif GM.mif - tissueRGB.mif -axis 3
# mrview tissueRGB.mif -odf.load_sh WM_FODs.mif



#*** 7. mtnormalise
print('Running mtnormalise on WM FODs')
for i in df.index.tolist(): #indices_to_rerun:
    sid = df.at[i,'Subject ID']
    dwi = df.at[i,'Path_dMRI_biascorrected']
    dwi_mask = dwi.replace('/dMRI_biascorrect/','/dMRI_mask/').replace('.mif','_mask.mif')
    wm_fods_mif = os.path.join(fods_folder,'%s_wm_fods_%s.mif' % (sid,rf_algo))
    wm_fods_normed = wm_fods_mif.replace('.mif','_normed.mif')
    
    cmd = ['mtnormalise',wm_fods_mif,wm_fods_normed,'-mask',dwi_mask]
    if not os.path.exists(wm_fods_normed):
        oot = subprocess.run(cmd)
        orc = oot.returncode
    else:
        print(' Normalise WN FOD exists for %s' % (sid))
        orc = 0
    if orc==0:
        print(' Subject %s: mtnormalise completed' % sid)
    else:
        print(' Subject %s: mtnormalise failed (return code %i)' % (sid,orc))


#*** 8. Tractography
print("      ...then tractography")
cmd_tckgen_cluster_freesurfer = 'cd /SAN/medic/Net_Mod_MS/cjd/ppmi_tractography; python ~/scripts/ppmi_tckgen_cluster.py freesurfer/5tt wm_fods freesurfer'

#*** 9. Connectomes
print("         ...then form the connectomes")
nodes_uint32_bool = True
batch_cmd = [
    'foreach', '-%i' % n_cpus,
    'sMRI_nodes/*nodes*.mif',':',
    'mrconvert','-datatype','uint32','IN','sMRI_nodes/PRE_uint32.mif'
]
if nodes_uint32_bool:
    foreach_check = os.system(' '.join(batch_cmd))
tck2connectome_bool = True
#* foreach fail...
# batch_cmd = [
#     'foreach', '-%i' % n_cpus,
#     'tck_freesurfer/*tcksift4M.tck',':', #'sMRI_freesurfer/*nodes*uint32.mif',':',
#     #'echo','sMRI_freesurfer/UNI*uint32.mif'
#     'tck2connectome','-assignment_all_voxels','IN','sMRI_freesurfer/UNI_uint32.mif','tck_freesurfer/PRE_connectome.tsv'
# ]
# if tck2connectome_bool:
#     foreach_check = os.system(' '.join(batch_cmd))
tail = '_tckgen32M_tckedit16M_tcksift4M.tck'
fold = os.path.join(act_folder,'tck_freesurfer')
label_txt = sMRI_5tt_txt.replace('Path_5tt_','')
for i in df.index.tolist(): #indices_to_rerun:
    sid = df.at[i,'Subject ID']
    tcksift = os.path.join(fold,os.path.basename(df.at[i,sMRI_5tt_txt]).split('.')[0] + tail)
    if os.path.exists(tcksift):
        nodes_uint32 = df.at[i,sMRI_nodes_txt].replace('.mif','_uint32.mif').replace('sMRI_freesurfer','sMRI_nodes')
        connectome_tsv = tcksift.replace('.tck',label_txt + '_connectome.tsv')
        cmd = ['tck2connectome','-assignment_all_voxels',tcksift,nodes_uint32,connectome_tsv]
        if not os.path.exists(connectome_tsv):
            oot = subprocess.run(cmd)
            orc = oot.returncode
        else:
            print(' Connectome exists for %s' % (sid))
            orc = 0
        if orc==0:
            print(' Subject %s: tck2connectome completed' % sid)
            df.at[i,'Path_connectome_' + label_txt] = connectome_tsv
        else:
            print(' Subject %s: tck2connectome failed (return code %i)' % (sid,orc))
        


print('Apparent errors in the following connectomes: 3114, 3767, 3768, 3857 \n All because dwi2mask produced a weird mask focussed on the spinal cord. Fixed the masks using dwi2mak on the raw images, then re-ran dwi2fod etc..')
print('Still awaiting the tractography for this one: 3765')




#!/usr/bin/env python3
# 
# Network spreading models using personalised (age/gender-matched) connectomes 
# from PPMI healthy controls for anatomically constrained tractography
#
# ACT prep (MRTrix3 pipeline)
#
# Prerequisite: walk through act_ppmi_0_prep.py
#
# <NAME>, February 2020

import os,subprocess,glob
import pandas as pd

#* Name of CSV file containing image file paths and IDs: generated by act_ppmi_0_prep.py
csv_input_file = 'act_ppmi_filepaths.csv'

#* File locations
top_folder = '~/ox_tools/anatomically_constrained_tractography'
data_folder = os.path.join(top_folder,'PPMI')
act_folder  = os.path.normpath(os.path.join(data_folder,'act'))
dMRI_folder = os.path.normpath(os.path.join(act_folder,'dMRI'))
sMRI_folder = os.path.normpath(os.path.join(act_folder,'sMRI'))
json_folder = os.path.normpath(os.path.join(act_folder,'json'))
csv_input_path = os.path.join(top_folder,csv_input_file)
act_csv_path = os.path.join(act_folder,'act_ppmi.csv')

#* Preprocessing functions copied from my connectome pipeline
wd = os.path.join(top_folder,'scripts')
os.chdir(wd)
from act_utilities import *

#* Preprocessing output folders
denoise_tail = '_denoised'
preproc_tail = '_preproc'
biascorrect_tail = '_biascorrect'
dMRI_tails = [denoise_tail,preproc_tail,biascorrect_tail]
dMRI_preproc_folders = [os.path.join(act_folder,'dMRI%s' % d) for d in dMRI_tails]
reg_tail = '_reg'
sMRI_reg_folder      =  os.path.join(act_folder,'sMRI%s' % reg_tail)
make_sure_path_exists(sMRI_reg_folder)
for f in dMRI_preproc_folders:
    make_sure_path_exists(f)

os.chdir(act_folder)


#####################################################################################################
############ Prep: check for dependencies, convert dMRI to MIF, prepare output DataFrame ############
#* Check for dependencies: ANTS, etc.
depCheckVal = checkForDependencies()
if depCheckVal!=0:
    print('ERROR: missing dependency {0}'.format(depCheckVal))
    print('see checkForDependencies() function for details')
    sys.exit(depCheckVal)
#* DataFrame for ACT files
if os.path.exists(act_csv_path):
    df = pd.read_csv(act_csv_path)
else:
    df_act = pd.read_csv(csv_input_path)
    df_act_sMRI = df_act.loc[df_act["Modality"]=='sMRI'].sort_values(by='Subject ID').reset_index(drop=True)
    df_act_dMRI = df_act.loc[df_act["Modality"]=='dMRI'].sort_values(by='Subject ID').reset_index(drop=True)
    col_preproc = "Path preprocessed"
    #* Match sMRI to dMRI using Subject ID
    cols = ["Subject ID","Path"]
    df_act_sMRI = df_act_sMRI[cols].rename(columns={"Path": "Path_sMRI"})
    df_act_dMRI = df_act_dMRI[cols].rename(columns={"Path": "Path_dMRI"})
    df = pd.merge(df_act_sMRI,df_act_dMRI)
    df["Preprocessed"] = False
#* dMRI: convert NII to MIF, with embedded b-values and b-vectors (and other metadata from json, including Phase Encoding Direction)
for i in df.index.values:
    sMRI = df.at[i,'Path_sMRI']
    dMRI = df.at[i,'Path_dMRI']
    dMRI_mif = dMRI.replace('.nii.gz','.mif')
    json_file = os.path.join(json_folder,os.path.basename(dMRI).replace('.nii.gz','.json'))
    #* Convert to MRtrix image format (MIF), with embedded b-values and b-vectors
    if not os.path.isfile(dMRI_mif):
        bvec = dMRI.replace('.nii.gz','.bvec')
        bval = dMRI.replace('.nii.gz','.bval')
        mrconvert_check = subprocess.run(['mrconvert','-json_import',json_file,'-fslgrad',bvec,bval,dMRI,dMRI_mif])
        #rm_check = subprocess.run(['rm',bvec,bval,dMRI])
    df.at[i,"Path_dMRI"] = dMRI_mif
#* Write ACT preprocessing status and results (file locations) to CSV
df.to_csv(act_csv_path,index=False)
#>>>>>>>>>>>> I fixed a missing file. Keeping the code here for future reference
# #* One DWI file was missing bvals/bvecs files, so I manually downloaded an alternative (same participant, earlier visit)
# patno = "3855"
# i = df.loc[df["Subject ID"]==3855].index.values[0]
# #* FIXME: Remove existing files (did it manually)
# # o = subprocess.run(['rm',os.path.join(sMRI_folder,patno+"*")])
# # o = subprocess.run(['rm',os.path.join(dMRI_folder,patno+"*")])
# zipfile = "/Users/noxtoby/Documents/EuroPOND-WP8-NetworkModels/tractography/PPMI/zips/PPMI-for-CJD-missing_grad.zip"
# o = subprocess.run(['open',zipfile])
# new_dicom_dMRI = "/Users/noxtoby/Documents/EuroPOND-WP8-NetworkModels/tractography/PPMI/zips/PPMI/3855/DTI_gated/2012-07-24_08_47_28.0/S204392"
# new_dicom_sMRI = "/Users/noxtoby/Documents/EuroPOND-WP8-NetworkModels/tractography/PPMI/zips/PPMI/3855/MPRAGE_GRAPPA/2012-07-24_08_38_32.0/S204393"
# nifti_folder = os.path.join(data_folder,'nifti')
# o = subprocess.run(['dcm2niix','-o',nifti_folder,'-f','%i_%t_%s_%p','-z','y',new_dicom_sMRI])
# o = subprocess.run(['dcm2niix','-o',nifti_folder,'-f','%i_%t_%s_%p','-z','y',new_dicom_dMRI])
# #* Convert dMRI to MIF, move to dMRI folder, delete old file
# old_dMRI = df.at[i,'Path_dMRI']
# new_dMRI_nifti = os.path.join(nifti_folder,"3855_20120724083832_2_DTI_gated.nii.gz")
# bval = new_dMRI_nifti.replace('.nii.gz','.bval')
# bvec = new_dMRI_nifti.replace('.nii.gz','.bvec')
# json_file = new_dMRI_nifti.replace('.nii.gz','.json')
# new_dMRI_mif = os.path.join(dMRI_folder,os.path.basename(new_dMRI_nifti)).replace('.nii.gz','.mif')
# mrconvert_check = subprocess.run(['mrconvert','-json_import',json_file,'-fslgrad',bvec,bval,new_dMRI_nifti,new_dMRI_mif])
# #* Move sMRI
# old_sMRI = df.at[i,'Path_sMRI']
# new_sMRI = os.path.join(nifti_folder,"3855_20120724083832_9_MPRAGE_GRAPPA.nii.gz")
# new_sMRI_loc = os.path.join(sMRI_folder,os.path.basename(new_sMRI))
# o = subprocess.run(['mv',new_sMRI,new_sMRI_loc])
# #* Fix paths for the new files
# df.at[i,"Path_dMRI"] = new_dMRI_mif
# df.at[i,"Path_sMRI"] = new_sMRI_loc
# df["Path_sMRI"] = df["Path_sMRI"].map(lambda x: x.split('.')[0] + '.nii.gz')
#<<<<<<<<<<<<
#####################################################################################################




#####################################################################################################
############ Preprocessing ############
#* Steps: 1. preprocess_dMRI => bias corrected dMRI
#*        2. register sMRI to dMRI_biascorrected
#* Cluster:
#*        3. run GIF3 on sMRI_registered
#*        4. ACT stuff: 5TT, nodes / WM_FODs (response functions & CSD)
#*        5. mtnormalise
#*        6. Tractography
#*        7. Connectomes
preproc_bool = False

# It might be prudent to use the cluster for this (DWI bias correction takes ages!)
#   => a) copy dMRI and sMRI folders onto cluster
#      b) run act_utilities.preprocess_individual() OR a version of act_utilities.preprocess_dMRI() and .preprocess_sMRI()

new_column_dMRI = 'Path_dMRI_biascorrected'
new_column_sMRI = 'Path_sMRI_registered'

#****** Local batch preprocessing using multicore CPU
n_cpus = 4
#* mrtrix3: dwidenoise => dwipreproc (-rpe_header)
batch_cmd = [
    'foreach', '-%i' % n_cpus,
    'dMRI/*.mif',':',
    'dwidenoise','IN','dMRI_{0}/PRE_{0}.mif'.format(denoise_tail),
    '"&&"',
    'dwipreproc','-rpe_header','dMRI_{0}/PRE_{0}.mif'.format(denoise_tail),'dMRI_{0}/PRE_{0}.mif'.format(preproc_tail)
]
if preproc_bool:
    foreach_check = os.system(' '.join(batch_cmd)) # Subprocess doesn't work! subprocess.run(batch_cmd) FIXME???
#* mrtrix3: dwibiascorrect
batch_cmd = [
    'foreach','-%i' % n_cpus,
    'dMRI%s/*.mif' % preproc_tail,':',
    'dwibiascorrect','-ants','IN','dMRI{0}/PRE{0}.mif'.format(biascorrect_tail)
]
if preproc_bool:
    foreach_check = os.system(' '.join(batch_cmd))

#* niftyreg: register sMRI to biascorrected dMRI
indices = df.index.values
cnt = 0
for i in indices:
    cnt += 1
    if not df.at[i,"Preprocessed"]:
        sMRI = df.at[i,'Path_sMRI']
        dMRI = df.at[i,'Path_dMRI']
        dMRI_denoised    = dMRI.replace('dMRI','dMRI%s' % denoise_tail).replace('.mif','%s.mif' % denoise_tail)
        dMRI_preproc     = dMRI.replace('dMRI','dMRI%s' % preproc_tail).replace('.mif','%s.mif' % preproc_tail)
        dMRI_biascorrect = dMRI.replace('dMRI','dMRI%s' % biascorrect_tail).replace('.mif','%s%s.mif' % (preproc_tail,biascorrect_tail))
        sMRI_reg = sMRI.replace('sMRI','sMRI%s' % reg_tail)
        print('************\n%i of %i - Preprocessing %s and %s ' % (cnt,len(indices),os.path.basename(dMRI),os.path.basename(sMRI)))
        #* Serial preprocessing: (dMRI_denoised,dMRI_preproc,dMRI_biascorrected,sMRI_reg,reg_aladin_output_transform) = preprocess_individual(dMRI_mif,sMRI)
        
        #* register sMRI to biascorrected dMRI
        if os.path.exists(dMRI_biascorrect):
            (reg_aladin_check,sMRI_reg,reg_aladin_output_transform) = reg_sMRI_to_dMRI(sMRI,dMRI_biascorrect)
            #*** Save to DataFrame
            df.at[i,new_column_dMRI] = dMRI_biascorrect
            df.at[i,new_column_sMRI] = sMRI_reg
            df.at[i,"Preprocessed"] = True
        else:
            print("Could not find bias-corrected file: %s " % dMRI_biascorrect)
            print(" >>> registration not performed")
            df.at[i,"Preprocessed"] = False
        #* Update CSV
        df.to_csv(act_csv_path,index=False)
    else:
        print("Files already preprocessed: %s, %s " % (os.path.basename(dMRI),os.path.basename(sMRI)))

#*** 5. run GIF3 / FreeSurfer on sMRI_registered
df_ = pd.read_csv(act_csv_path)
print("\nChecking whether you've run GIF on the sMRI (you'll have to change this for FreeSurfer, or remove: see below)\n")
algorithm = 'gif3'
sMRI_parc_txt = 'sMRI_%s' % algorithm
gif_parc_tag = "_NeuroMorph_Parcellation.nii.gz"
gif_folder = os.path.join(act_folder,sMRI_parc_txt)
gif3_parcellations_path = glob.glob(os.path.join(gif_folder,"*/*%s" % gif_parc_tag))
gif3_parcellations = [os.path.basename(g) for g in gif3_parcellations_path]
df_gif3 = pd.DataFrame(columns=[sMRI_parc_txt],data=gif3_parcellations)
df_gif3['Subject ID'] = df_gif3[sMRI_parc_txt].map(lambda x: x.split('_')[0]).astype(int)
df = df_.merge(df_gif3,on='Subject ID')
print(" > > > Found %i GIF3 parcellations out of %i" % (len(gif3_parcellations),df.shape[0]))

print("\nChecking whether you've run FreeSurfer on the sMRI\n")
algorithm = 'freesurfer'
sMRI_parc_txt = 'sMRI_%s' % algorithm
freesurfer_parc = "aparc+aseg.mgz"
freesurfer_folder = os.path.join(act_folder,sMRI_parc_txt)
freesurfer_parcellations_path = glob.glob(os.path.join(freesurfer_folder,"*/mri/%s" % freesurfer_parc))
#freesurfer_parcellations = freesurfer_parcellations_path # Could rename the files (to approximately match GIF filenames?)
df_freesurfer = pd.DataFrame(columns=[sMRI_parc_txt],data=freesurfer_parcellations_path)
df_freesurfer['Subject ID'] = df_freesurfer[sMRI_parc_txt].map(lambda x: x.split('/')[-3].split('_')[0]).astype(int)
df = df.merge(df_freesurfer,on='Subject ID')
print(" > > > Found %i FreeSurfer parcellations out of %i" % (len(freesurfer_parcellations_path),df.shape[0]))


#*** 5.2 DWI masks
print(">>> DWI masks for dwi2fod")
masks_folder = os.path.join(act_folder,'dMRI_mask')
oot = os.system('mkdir -p %s' % masks_folder)
#* mrtrix3: dwi2mask
mask_tail = '_mask'
batch_cmd = [
    'foreach','-%i' % n_cpus,
    'dMRI%s/*.mif' % biascorrect_tail,':',
    'dwi2mask','IN','dMRI{0}/PRE{0}.mif'.format(mask_tail)
]
if preproc_bool:
    foreach_check = os.system(' '.join(batch_cmd))

#*** 6. ACT stuff: 5TT, nodes / WM_FODs (response functions & CSD)
algorithm = 'freesurfer'
sMRI_parc_txt = 'sMRI_%s' % algorithm
bizzi_lobes = False
if bizzi_lobes:
    algorithm_ = algorithm + "_bizzi2lobes"
else:
    algorithm_ = algorithm
sMRI_5tt_txt = 'Path_5tt_%s' % algorithm_
sMRI_nodes_txt = 'Path_nodes_%s' % algorithm_
#* labelconvert: use bizzi2lobes_cingsep_convert.txt (c.f. MRTRIX3:path/to/labelconvertLUTs/f2lobes....txt)
print(">>> ACT prep: 5TT, nodes, CSD (WM_FODs and response functions)")
print("    *** Chosen algorithm: %s (%s)" % (algorithm,algorithm_))
fivett_folder = os.path.join(act_folder,'sMRI_5tt')
oot = os.system('mkdir -p %s' % fivett_folder)
fttrc = -1
# df.drop(columns=[sMRI_nodes_txt],inplace=True)
for i in df.index.values:
    participant_id = df.at[i,'Subject ID']
    sMRI_reg = os.path.join(act_folder,'sMRI_reg',os.path.basename(df.at[i,'Path_sMRI_registered']))
    if algorithm == 'gif3':
        sMRI_gif = os.path.join(act_folder,sMRI_parc_txt,os.path.basename(df.at[i,sMRI_parc_txt]).replace(gif_parc_tag,'_share_apps_cmic_GIF_bin'),os.path.basename(df.at[i,sMRI_parc_txt]))
        sMRI_parc = sMRI_gif
        sMRI_5TT = os.path.join(fivett_folder,os.path.basename(sMRI_gif).replace(gif_parc_tag,'_5tt_%s.mif' % algorithm))
    elif algorithm == 'freesurfer':
        sMRI_freesurfer = df.at[i,sMRI_parc_txt]
        sMRI_parc = sMRI_freesurfer
        sMRI_5TT = os.path.join(fivett_folder,os.path.basename(sMRI_reg).replace('.nii.gz','_5tt_%s.mif' % algorithm))
    if not os.path.isfile(sMRI_5TT):
        fiveTT_check = subprocess.run(['5ttgen',algorithm,sMRI_parc,sMRI_5TT])
        fttrc = fiveTT_check.returncode
    else:
        print('5TT file found')
        fttrc = 0
    if fttrc == 0:
        df.at[i,sMRI_5tt_txt] = sMRI_5TT
        print('\t> > > \n\t> > > 5ttgen done. Proceeding with labelconvert\n\t> > >')
        (labelconvert_check,nodes) = convert_sMRI_labels(sMRI_parc,algorithm=algorithm,cluster=False,bizzi=bizzi_lobes)
        if labelconvert_check==0:
            df.at[i,sMRI_nodes_txt] = nodes
        else:
            df.at[i,sMRI_nodes_txt] = np.nan
    else:
        df.at[i,sMRI_5tt_txt] = np.nan
#* Move nodes files
nodes_folder = 'sMRI_nodes'
oot = os.system('mkdir -p %s' % os.path.join(act_folder,nodes_folder))
for i in df.index.values:
    nodes = df.at[i,sMRI_nodes_txt]
    nodes_new = nodes.replace("/"+ sMRI_parc_txt +"/","/"+ nodes_folder +"/")
    oot = subprocess.run(['mv',nodes,nodes_new])
    if oot.returncode==0:
        df.at[i,"Path_nodes"] = nodes_new

#* WMFODs: instead could do on the cluster via ppmi_tckgen_cluster
rf_algo = 'tournier' #'msmt_5tt','dhollander'
fods_folder = os.path.join(act_folder,'dMRI_fods'); os.system('mkdir -p %s' % fods_folder)
#* Re-run these after fixing the masks
indices_to_rerun = df.index[df['Subject ID'].isin([3114, 3765, 3767, 3768, 3857])]

for i in df.index.tolist(): #indices_to_rerun: #
    sid = df.at[i,'Subject ID']
    dwi = df.at[i,'Path_dMRI_biascorrected']
    dwi_mask = dwi.replace('/dMRI_biascorrect/','/dMRI_mask/').replace('.mif','_mask.mif')
    fivett = df.at[i,sMRI_5tt_txt]
    wm_response_txt  = os.path.join(fods_folder,'%s_wm_response_%s.txt' % (sid,rf_algo))
    gm_response_txt  = os.path.join(fods_folder,'%s_gm_response_%s.txt' % (sid,rf_algo))
    csf_response_txt = os.path.join(fods_folder,'%s_csf_response_%s.txt' % (sid,rf_algo))
    voxels_mif       = os.path.join(fods_folder,'%s_voxels_%s.mif' % (sid,rf_algo))
    #* dwi2response
    # > PPMI has only two b-values, which precludes multi-shell multi-tissue CSD, 
    # > so use 'dwi2response tournier' rather than msmt_5tt or dhollander 
    # > (which will output matrices rather than vectors in the RF txt file)
    #cmd = ['dwi2response','dhollander','-shells','0,1000',dwi,wm_response_txt,gm_response_txt,csf_response_txt,'-voxels',voxels_mif]
    #cmd = ['dwi2response','msmt_5tt',dwi,fivett,wm_response_txt,gm_response_txt,csf_response_txt,'-voxels',voxels_mif]
    cmd = ['dwi2response',rf_algo,dwi,wm_response_txt]
    if not os.path.exists(wm_response_txt):
        oot = subprocess.run(cmd)
        orc = oot.returncode
    else:
        print(' Subject %s: %s exists' % (sid,wm_response_txt))
        orc = 0
    if orc==0:
        print(' Subject %s: dwi2response completed' % sid)
    else:
        print(' Subject %s: dwi2response failed (return code %i)' % (sid,orc))
    #* dwi2fod
    wm_fods_mif = os.path.join(fods_folder,'%s_wm_fods_%s.mif' % (sid,rf_algo))
    gm_fods_mif = os.path.join(fods_folder,'%s_gm_fods_%s.mif' % (sid,rf_algo))
    csf_fods_mif = os.path.join(fods_folder,'%s_csf_fods_%s.mif' % (sid,rf_algo))
    cmd = ['dwi2fod','csd','-mask',dwi_mask,dwi,wm_response_txt,wm_fods_mif] 
    #cmd = ['dwi2fod','csd',dwi,gm_response_txt,gm_fods_mif,csf_response_txt,csf_fods_mif]
    if not os.path.exists(wm_fods_mif):
        oot = subprocess.run(cmd)
        #os.system(' '.join(cmd))
        orc = oot.returncode
    else:
        print(' Subject %s: %s exists' % (sid,wm_fods_mif))
    if orc==0:
        print(' Subject %s: dwi2fod completed' % sid)
        df.at[i,"Path_wm_fod"] = wm_fods_mif
        #df.at[i,"Path_gm_fod"] = gm_fods_mif
        #df.at[i,"Path_csf_fod"] = csf_fods_mif
    else:
        print(' Subject %s: dwi2fod failed (return code %i)' % (sid,orc))

#* Check the response functions:
# shview wm_response.txt

# #* MSMT-only: Check the FODs: convert to RGB image
# for i in df.index.tolist():
#     sid = df.at[i,'Subject ID']
#     wm_fods_mif = os.path.join(fods_folder,'%s_wm_fods_%s.mif' % (sid,rf_algo))
#     tissue_rgb_mif = wm_fods_mif.replace('wm_fods','tissueRGB')
#     cmd = 'mrconvert %s - -coord 3 0 | mrcat CSF.mif GM.mif - %s -axis 3' % (wm_fods_mif,tissue_rgb_mif)
#     if os.path.exists(wm_fods_mif):
#         df.at[i,"Path_wm_fod"] = wm_fods_mif

# mrconvert WM_FODs.mif - -coord 3 0 | mrcat CSF.mif GM.mif - tissueRGB.mif -axis 3
# mrview tissueRGB.mif -odf.load_sh WM_FODs.mif



#*** 7. mtnormalise
print('Running mtnormalise on WM FODs')
for i in df.index.tolist(): #indices_to_rerun:
    sid = df.at[i,'Subject ID']
    dwi = df.at[i,'Path_dMRI_biascorrected']
    dwi_mask = dwi.replace('/dMRI_biascorrect/','/dMRI_mask/').replace('.mif','_mask.mif')
    wm_fods_mif = os.path.join(fods_folder,'%s_wm_fods_%s.mif' % (sid,rf_algo))
    wm_fods_normed = wm_fods_mif.replace('.mif','_normed.mif')
    
    cmd = ['mtnormalise',wm_fods_mif,wm_fods_normed,'-mask',dwi_mask]
    if not os.path.exists(wm_fods_normed):
        oot = subprocess.run(cmd)
        orc = oot.returncode
    else:
        print(' Normalise WN FOD exists for %s' % (sid))
        orc = 0
    if orc==0:
        print(' Subject %s: mtnormalise completed' % sid)
    else:
        print(' Subject %s: mtnormalise failed (return code %i)' % (sid,orc))


#*** 8. Tractography
print("      ...then tractography")
cmd_tckgen_cluster_freesurfer = 'cd /SAN/medic/Net_Mod_MS/cjd/ppmi_tractography; python ~/scripts/ppmi_tckgen_cluster.py freesurfer/5tt wm_fods freesurfer'

#*** 9. Connectomes
print("         ...then form the connectomes")
nodes_uint32_bool = True
batch_cmd = [
    'foreach', '-%i' % n_cpus,
    'sMRI_nodes/*nodes*.mif',':',
    'mrconvert','-datatype','uint32','IN','sMRI_nodes/PRE_uint32.mif'
]
if nodes_uint32_bool:
    foreach_check = os.system(' '.join(batch_cmd))
tck2connectome_bool = True
#* foreach fail...
# batch_cmd = [
#     'foreach', '-%i' % n_cpus,
#     'tck_freesurfer/*tcksift4M.tck',':', #'sMRI_freesurfer/*nodes*uint32.mif',':',
#     #'echo','sMRI_freesurfer/UNI*uint32.mif'
#     'tck2connectome','-assignment_all_voxels','IN','sMRI_freesurfer/UNI_uint32.mif','tck_freesurfer/PRE_connectome.tsv'
# ]
# if tck2connectome_bool:
#     foreach_check = os.system(' '.join(batch_cmd))
tail = '_tckgen32M_tckedit16M_tcksift4M.tck'
fold = os.path.join(act_folder,'tck_freesurfer')
label_txt = sMRI_5tt_txt.replace('Path_5tt_','')
for i in df.index.tolist(): #indices_to_rerun:
    sid = df.at[i,'Subject ID']
    tcksift = os.path.join(fold,os.path.basename(df.at[i,sMRI_5tt_txt]).split('.')[0] + tail)
    if os.path.exists(tcksift):
        nodes_uint32 = df.at[i,sMRI_nodes_txt].replace('.mif','_uint32.mif').replace('sMRI_freesurfer','sMRI_nodes')
        connectome_tsv = tcksift.replace('.tck',label_txt + '_connectome.tsv')
        cmd = ['tck2connectome','-assignment_all_voxels',tcksift,nodes_uint32,connectome_tsv]
        if not os.path.exists(connectome_tsv):
            oot = subprocess.run(cmd)
            orc = oot.returncode
        else:
            print(' Connectome exists for %s' % (sid))
            orc = 0
        if orc==0:
            print(' Subject %s: tck2connectome completed' % sid)
            df.at[i,'Path_connectome_' + label_txt] = connectome_tsv
        else:
            print(' Subject %s: tck2connectome failed (return code %i)' % (sid,orc))
        


print('Apparent errors in the following connectomes: 3114, 3767, 3768, 3857 \n All because dwi2mask produced a weird mask focussed on the spinal cord. Fixed the masks using dwi2mak on the raw images, then re-ran dwi2fod etc..')
print('Still awaiting the tractography for this one: 3765')
#!/usr/bin/env python3 # # Network spreading models using personalised (age/gender-matched) connectomes # from PPMI healthy controls for anatomically constrained tractography # # ACT prep (MRTrix3 pipeline) # # Prerequisite: walk through act_ppmi_0_prep.py # # <NAME>, February 2020 #* Name of CSV file containing image file paths and IDs: generated by act_ppmi_0_prep.py #* File locations #* Preprocessing functions copied from my connectome pipeline #* Preprocessing output folders ##################################################################################################### ############ Prep: check for dependencies, convert dMRI to MIF, prepare output DataFrame ############ #* Check for dependencies: ANTS, etc. #* DataFrame for ACT files #* Match sMRI to dMRI using Subject ID #* dMRI: convert NII to MIF, with embedded b-values and b-vectors (and other metadata from json, including Phase Encoding Direction) #* Convert to MRtrix image format (MIF), with embedded b-values and b-vectors #rm_check = subprocess.run(['rm',bvec,bval,dMRI]) #* Write ACT preprocessing status and results (file locations) to CSV #>>>>>>>>>>>> I fixed a missing file. Keeping the code here for future reference # #* One DWI file was missing bvals/bvecs files, so I manually downloaded an alternative (same participant, earlier visit) # patno = "3855" # i = df.loc[df["Subject ID"]==3855].index.values[0] # #* FIXME: Remove existing files (did it manually) # # o = subprocess.run(['rm',os.path.join(sMRI_folder,patno+"*")]) # # o = subprocess.run(['rm',os.path.join(dMRI_folder,patno+"*")]) # zipfile = "/Users/noxtoby/Documents/EuroPOND-WP8-NetworkModels/tractography/PPMI/zips/PPMI-for-CJD-missing_grad.zip" # o = subprocess.run(['open',zipfile]) # new_dicom_dMRI = "/Users/noxtoby/Documents/EuroPOND-WP8-NetworkModels/tractography/PPMI/zips/PPMI/3855/DTI_gated/2012-07-24_08_47_28.0/S204392" # new_dicom_sMRI = "/Users/noxtoby/Documents/EuroPOND-WP8-NetworkModels/tractography/PPMI/zips/PPMI/3855/MPRAGE_GRAPPA/2012-07-24_08_38_32.0/S204393" # nifti_folder = os.path.join(data_folder,'nifti') # o = subprocess.run(['dcm2niix','-o',nifti_folder,'-f','%i_%t_%s_%p','-z','y',new_dicom_sMRI]) # o = subprocess.run(['dcm2niix','-o',nifti_folder,'-f','%i_%t_%s_%p','-z','y',new_dicom_dMRI]) # #* Convert dMRI to MIF, move to dMRI folder, delete old file # old_dMRI = df.at[i,'Path_dMRI'] # new_dMRI_nifti = os.path.join(nifti_folder,"3855_20120724083832_2_DTI_gated.nii.gz") # bval = new_dMRI_nifti.replace('.nii.gz','.bval') # bvec = new_dMRI_nifti.replace('.nii.gz','.bvec') # json_file = new_dMRI_nifti.replace('.nii.gz','.json') # new_dMRI_mif = os.path.join(dMRI_folder,os.path.basename(new_dMRI_nifti)).replace('.nii.gz','.mif') # mrconvert_check = subprocess.run(['mrconvert','-json_import',json_file,'-fslgrad',bvec,bval,new_dMRI_nifti,new_dMRI_mif]) # #* Move sMRI # old_sMRI = df.at[i,'Path_sMRI'] # new_sMRI = os.path.join(nifti_folder,"3855_20120724083832_9_MPRAGE_GRAPPA.nii.gz") # new_sMRI_loc = os.path.join(sMRI_folder,os.path.basename(new_sMRI)) # o = subprocess.run(['mv',new_sMRI,new_sMRI_loc]) # #* Fix paths for the new files # df.at[i,"Path_dMRI"] = new_dMRI_mif # df.at[i,"Path_sMRI"] = new_sMRI_loc # df["Path_sMRI"] = df["Path_sMRI"].map(lambda x: x.split('.')[0] + '.nii.gz') #<<<<<<<<<<<< ##################################################################################################### ##################################################################################################### ############ Preprocessing ############ #* Steps: 1. preprocess_dMRI => bias corrected dMRI #*        2. register sMRI to dMRI_biascorrected #* Cluster: #*        3. run GIF3 on sMRI_registered #*        4. ACT stuff: 5TT, nodes / WM_FODs (response functions & CSD) #*        5. mtnormalise #*        6. Tractography #*        7. Connectomes # It might be prudent to use the cluster for this (DWI bias correction takes ages!) #   => a) copy dMRI and sMRI folders onto cluster #      b) run act_utilities.preprocess_individual() OR a version of act_utilities.preprocess_dMRI() and .preprocess_sMRI() #****** Local batch preprocessing using multicore CPU #* mrtrix3: dwidenoise => dwipreproc (-rpe_header) # Subprocess doesn't work! subprocess.run(batch_cmd) FIXME??? #* mrtrix3: dwibiascorrect #* niftyreg: register sMRI to biascorrected dMRI #* Serial preprocessing: (dMRI_denoised,dMRI_preproc,dMRI_biascorrected,sMRI_reg,reg_aladin_output_transform) = preprocess_individual(dMRI_mif,sMRI) #* register sMRI to biascorrected dMRI #*** Save to DataFrame #* Update CSV #*** 5. run GIF3 / FreeSurfer on sMRI_registered #freesurfer_parcellations = freesurfer_parcellations_path # Could rename the files (to approximately match GIF filenames?) #*** 5.2 DWI masks #* mrtrix3: dwi2mask #*** 6. ACT stuff: 5TT, nodes / WM_FODs (response functions & CSD) #* labelconvert: use bizzi2lobes_cingsep_convert.txt (c.f. MRTRIX3:path/to/labelconvertLUTs/f2lobes....txt) # df.drop(columns=[sMRI_nodes_txt],inplace=True) #* Move nodes files #* WMFODs: instead could do on the cluster via ppmi_tckgen_cluster #'msmt_5tt','dhollander' #* Re-run these after fixing the masks #indices_to_rerun: # #* dwi2response # > PPMI has only two b-values, which precludes multi-shell multi-tissue CSD, # > so use 'dwi2response tournier' rather than msmt_5tt or dhollander # > (which will output matrices rather than vectors in the RF txt file) #cmd = ['dwi2response','dhollander','-shells','0,1000',dwi,wm_response_txt,gm_response_txt,csf_response_txt,'-voxels',voxels_mif] #cmd = ['dwi2response','msmt_5tt',dwi,fivett,wm_response_txt,gm_response_txt,csf_response_txt,'-voxels',voxels_mif] #* dwi2fod #cmd = ['dwi2fod','csd',dwi,gm_response_txt,gm_fods_mif,csf_response_txt,csf_fods_mif] #os.system(' '.join(cmd)) #df.at[i,"Path_gm_fod"] = gm_fods_mif #df.at[i,"Path_csf_fod"] = csf_fods_mif #* Check the response functions: # shview wm_response.txt # #* MSMT-only: Check the FODs: convert to RGB image # for i in df.index.tolist(): #     sid = df.at[i,'Subject ID'] #     wm_fods_mif = os.path.join(fods_folder,'%s_wm_fods_%s.mif' % (sid,rf_algo)) #     tissue_rgb_mif = wm_fods_mif.replace('wm_fods','tissueRGB') #     cmd = 'mrconvert %s - -coord 3 0 | mrcat CSF.mif GM.mif - %s -axis 3' % (wm_fods_mif,tissue_rgb_mif) #     if os.path.exists(wm_fods_mif): #         df.at[i,"Path_wm_fod"] = wm_fods_mif # mrconvert WM_FODs.mif - -coord 3 0 | mrcat CSF.mif GM.mif - tissueRGB.mif -axis 3 # mrview tissueRGB.mif -odf.load_sh WM_FODs.mif #*** 7. mtnormalise #indices_to_rerun: #*** 8. Tractography #*** 9. Connectomes #* foreach fail... # batch_cmd = [ #     'foreach', '-%i' % n_cpus, #     'tck_freesurfer/*tcksift4M.tck',':', #'sMRI_freesurfer/*nodes*uint32.mif',':', #     #'echo','sMRI_freesurfer/UNI*uint32.mif' #     'tck2connectome','-assignment_all_voxels','IN','sMRI_freesurfer/UNI_uint32.mif','tck_freesurfer/PRE_connectome.tsv' # ] # if tck2connectome_bool: #     foreach_check = os.system(' '.join(batch_cmd)) #indices_to_rerun: