1#!/usr/local/bin/python3.8 2__version__="v2.5 beta6" 3""" 4# Multi-Echo ICA, Version %s 5# See http://dx.doi.org/10.1016/j.neuroimage.2011.12.028 6# 7# Kundu, P., Brenowitz, N.D., Voon, V., Worbe, Y., Vertes, P.E., Inati, S.J., Saad, Z.S., 8# Bandettini, P.A. & Bullmore, E.T. Integrated strategy for improving functional 9# connectivity mapping using multiecho fMRI. PNAS (2013). 10# 11# Kundu, P., Inati, S.J., Evans, J.W., Luh, W.M. & Bandettini, P.A. Differentiating 12# BOLD and non-BOLD signals in fMRI time series using multi-echo EPI. NeuroImage (2011). 13# 14# meica.py version 2.5 (c) 2014 Prantik Kundu 15# PROCEDURE 1c : Preprocess multi-echo datasets and apply multi-echo ICA based on spatial concatenation 16# -Calcluation of functional-anatomical coregistration using EPI gray matter + local Pearson correlation method 17# 18#alignp_mepi_anat.py V.1.0 19#Compute warp parameters for co-registration from anatomical (skullstripped) to S0 and T2* volume (not skullstripped) 20#S0 and T2* volumes should be exactly in register 21# 22""" % (__version__) 23 24import sys 25import commands 26from re import split as resplit 27import re 28from os import system,getcwd,mkdir,chdir,popen 29import os.path 30from string import rstrip,split 31from optparse import OptionParser,OptionGroup 32 33parser=OptionParser() 34parser.add_option('-t',"",dest='t2s',help="T2* volume",default='') 35parser.add_option('-s',"",dest='s0',help="Skull-stripped S0 weighted volume, optional, for masking T2*",default='') 36parser.add_option('-a',"",dest='anat',help="Anatomical volume",default='') 37parser.add_option('-p',"",dest='prefix',help="Alignment matrix prefix" ,default='') 38parser.add_option('',"--cmass",action='store_true',dest='cmass',help="Align cmass before main co-registration" ,default=False) 39parser.add_option('',"--maxrot",dest='maxrot',help="Maximum rotation, default 30" ,default='30') 40parser.add_option('',"--maxshift",dest='maxshf',help="Maximum shift, default 30" ,default='30') 41parser.add_option('',"--maxscl",dest='maxscl',help="Maximum shift, default 30" ,default='1.2') 42(options,args) = parser.parse_args() 43 44""" 45Set up and cd into directory 46""" 47sl = [] 48startdir = os.path.abspath(os.path.curdir) 49walignp_dirname = "alignp.%s" % (options.prefix) 50sl.append("rm -rf %s " % walignp_dirname) 51sl.append("mkdir %s " % walignp_dirname) 52sl.append("cd %s " % walignp_dirname) 53 54def dsprefix(idn): 55 def prefix(datasetname): 56 return split(datasetname,'+')[0] 57 if len(split(idn,'.'))!=0: 58 if split(idn,'.')[-1]=='HEAD' or split(idn,'.')[-1]=='BRIK' or split(idn,'.')[-2:]==['BRIK','gz']: 59 return prefix(idn) 60 elif split(idn,'.')[-1]=='nii' and not split(idn,'.')[-1]=='nii.gz': 61 return '.'.join(split(idn,'.')[:-1]) 62 elif split(idn,'.')[-2:]==['nii','gz']: 63 return '.'.join(split(idn,'.')[:-2]) 64 else: 65 return prefix(idn) 66 else: 67 return prefix(idn) 68 69def dssuffix(idna): 70 suffix = idna.split(dsprefix(idna))[-1] 71 #print suffix 72 spl_suffix=suffix.split('.') 73 #print spl_suffix 74 if len(spl_suffix[0])!=0 and spl_suffix[0][0] == '+': return spl_suffix[0] 75 else: return suffix 76 77def niibrik(nifti): return dsprefix(nifti)+'+orig' 78 79def import_datasets(dsets): 80 """ 81 Base is functional 82 Dset is anatomical 83 """ 84 outnames = [] 85 for dset in dsets: 86 if dset!='' and dset!=None: 87 indir='' 88 if dset[0]!='/': indir = startdir 89 sl.append("3dcopy -overwrite %s/%s %s/%s/%s" % (indir,dset,startdir,walignp_dirname,dsprefix(os.path.basename(dset)))) 90 sl.append("3drefit -view orig `ls %s/%s/%s+*.HEAD`" % (startdir,walignp_dirname,dsprefix(os.path.basename(dset)))) 91 impdset = '%s+orig' % dsprefix(os.path.basename(dset)) 92 outnames.append(niibrik(impdset)) 93 else: 94 outnames.append('') 95 return outnames 96 97def align_centers(base,dset): 98 sl.append("@Align_Centers -base %s+orig -dset %s+orig" % (dsprefix(base),dsprefix(dset))) 99 volname = dsprefix(dset) + '_shft+orig' 100 matname = dsprefix(dset) + '_shft.1D' 101 return volname,matname 102 103def makeflat(dset): 104 flatname = "%s_flat" % (dsprefix(dset)) 105 sl.append("3dcalc -a %s -expr 'step(a)' -prefix %s" % (dset,flatname) ) 106 return niibrik(flatname) 107 108def graywt(t2sname,s0name): 109 basevol='align_base.nii.gz' 110 weightvol='align_weight.nii.gz' 111 if s0name!='': 112 sl.append("3dcalc -overwrite -a %s -b %s -expr 'a*ispositive(a)*step(b)' -prefix t2svm_ss.nii.gz" % (t2sname,s0name) ) 113 t2sname='t2svm_ss.nii.gz' 114 sl.append("3dBrickStat -mask %s -percentile 50 1 50 %s > graywt_thr.1D" % (t2sname,t2sname) ) 115 sl.append("maxthr=`cat graywt_thr.1D`; maxthra=($maxthr); vmax=${maxthra[1]}" ) 116 sl.append("3dcalc -overwrite -a %s -expr \"a*isnegative(a-2*${vmax})\" -prefix t2svm_thr.nii.gz" % (t2sname) ) 117 sl.append("3dUnifize -overwrite -prefix align_base.nii.gz t2svm_thr.nii.gz") 118 sl.append("3dSeg -prefix Segsy.t2svm -anat %s -mask %s" % (basevol,basevol)) 119 sl.append("3dcalc -overwrite -a 'Segsy.t2svm/Posterior+orig[2]' -prefix %s -expr 'a' " % weightvol) 120 return basevol,weightvol 121 122def allineate(sourcevol, weight, targetvol, prefix, maxrot, maxshf,maxscl,do_cmass): 123 """ 124 source should be anatomical 125 target should be T2* 126 """ 127 outvol_prefix = "%s_al" % prefix 128 outmat_prefix = "%s_al_mat" % prefix 129 cmass_opt = '' 130 if do_cmass: cmass_opt = '-cmass' 131 align_opts = "-lpc -weight %s -maxshf %s -maxrot %s -maxscl %s %s" % (weightvol, maxrot, maxshf,maxscl,cmass_opt) 132 sl.append("3dAllineate -overwrite -weight_frac 1.0 -VERB -warp aff -source_automask+2 -master SOURCE -source %s -base %s -prefix ./%s -1Dmatrix_save ./%s %s " \ 133 % (sourcevol,targetvol,outvol_prefix,outmat_prefix,align_opts)) 134 outvol_name = niibrik(outvol_prefix) 135 outmat_name = outmat_prefix + 'blah' 136 return outvol_name, outmat_name 137 138def aligntest(target,source,testname): 139 sl.append("3dresample -overwrite -master %s -inset %s -prefix altestA -rmode NN" % (target,source)) 140 sl.append("echo `3ddot altestA+orig %s` > %s.txt" % (target,testname)) 141 142def runproc(script_prefix, scrlines): 143 script_name = "_walignp_mepi_anat_%s.sh" % script_prefix 144 ofh = open(script_name,'w') 145 ofh.write("\n".join(scrlines)) 146 ofh.close() 147 os.system("bash %s" % script_name) 148 149"""Import datasets""" 150t2s_name,s0_name,anat_name = import_datasets([options.t2s,options.s0,options.anat]) 151 152"""Filter and segment T2* volume to produce T2* base and weight volumes""" 153basevol,weightvol = graywt(t2s_name,s0_name) 154 155"""Run 3dAllineate""" 156allin_volume,allin_matrix = allineate(anat_name,weightvol,basevol,options.prefix,options.maxrot,options.maxshf,options.maxscl,options.cmass) 157 158"""Run procedure""" 159runproc(options.prefix, sl) 160 161 162 163 164 165 166 167 168 169 170