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