-
Notifications
You must be signed in to change notification settings - Fork 7
/
generateFileStructure.py
executable file
·237 lines (186 loc) · 9.03 KB
/
generateFileStructure.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
#!/usr/bin/env python
#Code to generate a full shell of diffusion-weighted, eddy distorted images using FSL's possum, along with data that can be used to
#establish a ground truth.
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
# Handle arguments (before slow imports so --help can be fast)
import argparse
def str2bool(v):
#Function allows boolean arguments to take a wider variety of inputs
if v.lower() in ('yes', 'true', 't', 'y', '1'):
return True
elif v.lower() in ('no', 'false', 'f', 'n', '0'):
return False
else:
raise argparse.ArgumentTypeError('Boolean value expected.')
parser = argparse.ArgumentParser(description="Setup all the files required to run the simulations.")
parser.add_argument("possum_dir",help="Path to the possum simulation directory.")
parser.add_argument("output_dir",help="Path to output directory.")
parser.add_argument("bvals",help="Path to bval file.")
parser.add_argument("bvecs",help="Path to bvec file.")
parser.add_argument("--num_images",help='Number of volumes. Defaults to number of entries in bval file.',type=int)
parser.add_argument("--motion_dir",help='Path to directory describing subject motion.')
parser.add_argument("--interleave_factor",help="Interleave factor for slice-acquistion. Default=1. ",type=int,default=1)
parser.add_argument("--brain",help='Path to POSSUM input object.')
parser.add_argument("--brain_diffusion",help='Path to directory containing spherical harmonic coefficients for input object.')
parser.add_argument("--generate_artefact_free",help='Generate datasets without eddy-current and motion artefacts. Default=True.', type=str2bool, nargs='?',const=True,default=True)
parser.add_argument("--generate_distorted",help='Generate datasets with eddy-current and motion artefacts. Default=False', type=str2bool, nargs='?',const=True,default=False)
parser.add_argument("--matlabDir",help='Path to the Matlab run file',type=str,nargs='?',default='/usr/local/bin/matlab')
args=parser.parse_args()
import os
from subprocess import call
from dipy.io import read_bvals_bvecs
import dipy.reconst.shm as shm
import nibabel as nib
from Library import possumLib as pl
import numpy as np
import shutil
#Check arguments and setup paths
simDir = os.path.abspath(args.possum_dir)
outputDir = args.output_dir
bvalDir = args.bvals
bvecDir = args.bvecs
#Check bval entries
bvals, bvecs = read_bvals_bvecs(
args.bvals,
args.bvecs)
#Check num_images
num_bvals = len(bvals)
if args.num_images:
numImages= args.num_images
if numImages > num_bvals:
print("Warning: num_images cannot be greater than the number of entries in the bval file. Setting to this maximum.")
numImages = num_bvals
else:
numImages= num_bvals
if args.motion_dir == None:
motionDir='None'
else:
motionDir = args.motion_dir
if args.brain== None:
segName = 'Files/Segmentations/HCP_seg_clipped.nii.gz'
else:
segName = args.brain
if args.brain_diffusion== None:
sharm_dir = 'Files/SphericalHarmonics'
else:
sharm_dir = args.brain_diffusion
interleaveFactor = args.interleave_factor
normalImages = args.generate_artefact_free
motionAndEddyImages = args.generate_distorted
#Set eddy current distortion parameters
ep =0.006 #default is 0.006
tau = 0.1
basicSettings = [.0]*4
basicSettings[0]=0.001 #Leave a short gap before first gradient pulse /s
basicSettings[1]=0.010 #Pulse width / s
basicSettings[2]=0.025 #Diffusion time/s
basicSettings[3]=0.06 #Diffusion gradient strength T/m
eddyGradients = 'decaying' #Flat or decaying
#Set directories
codeDir = os.path.abspath('.')
#Load in segmentations
print('Loading segmentation...')
segmentedBrain, segmentedBrainData = pl.loadSegData(segName)
print('Finished loading segmentation.')
#Load in spherical harmonic coefficients
order = 8
#Check for problems first
if any(bval > 2500 for bval in bvals):
raise NotImplementedError('bvals > 2000 currently not supported')
if any(bval > 100 and bval < 1500 for bval in bvals):
print('Loading b=1000 spherical harmonics...')
coefficientsNiib1000 = nib.load(os.path.join(sharm_dir,'coefficientsUpsampledb1000n'+str(order)+'.nii.gz'))
coefficientsb1000 = coefficientsNiib1000.get_fdata()
print('Finished loading b=1000 spherical harmonics.')
if any(bval > 1500 and bval < 2500 for bval in bvals):
print('Loading b=2000 spherical harmonics...')
coefficientsNiib2000 = nib.load(os.path.join(sharm_dir,'coefficientsUpsampledb2000n'+str(order)+'.nii.gz'))
coefficientsb2000 = coefficientsNiib2000.get_fdata()
print('Finished loading b=2000 spherical harmonics.')
bvals, bvecs = read_bvals_bvecs(
bvalDir,
bvecDir)
print('Output directory: ' + outputDir)
#Make directory for cluster files
pl.makeFolder(outputDir)
pl.makeFolder(outputDir+"/Results")
pl.makeFolder(outputDir+"/Distortions")
pl.makeFolder(outputDir+"/Distortions/Motion")
pl.makeFolder(outputDir+"/Distortions/Eddy")
if motionAndEddyImages == "on":
shutil.copytree(motionDir, outputDir+"/Distortions/Motion")
shutil.copy(simDir+"/pulse",outputDir+"/Distortions/Eddy")
#Move ref brain for registering
shutil.copy(simDir+"/brainref.nii.gz", outputDir)
for index, bvec in enumerate(bvecs[0:numImages]):
#This workaround lets you carry on if generating is interrupted
if index < 0:
pass
else:
#Make directory for each setting
outputDirDirection = outputDir+"/Direction"+str(index)
pl.makeFolder(outputDirDirection)
#Copy needed files to folder
shutil.copy(simDir+"/MRpar",outputDirDirection)
shutil.copy(simDir+"/motion",outputDirDirection)
shutil.copy(simDir+"/slcprof",outputDirDirection)
shutil.copy(simDir+"/pulse.info",outputDirDirection)
shutil.copy(simDir+"/pulse.readme",outputDirDirection)
shutil.copy(simDir+"/pulse.posx",outputDirDirection)
shutil.copy(simDir+"/pulse.posy",outputDirDirection)
shutil.copy(simDir+"/pulse.posz",outputDirDirection)
shutil.copy(simDir+"/pulse",outputDirDirection)
#Get attenuated segmentations
#First rotate bvec
#bvecRotated = pl.rotateBvecs(bvecs[index], motionParams[index,4:]);
#Workaround: don't rotate bvec
bvecRotated = bvecs[index]
#print bvecs[index]
#print bvecRotated
if bvals[index] < 100:
attenuatedBrainData = segmentedBrainData
else:
#Convert bvecs to angles
x = bvecRotated[0]
y = bvecRotated[1]
z = bvecRotated[2]
r, theta, phi = shm.cart2sphere(x, y, z)
#Make design matrix
B, m, n = shm.real_sym_sh_basis(order, theta, phi)
#Get attenuated data
print('Attenuating volume ' + str(index))
if bvals[index] < 1500:
attenuatedBrainData = pl.attenuateImageSphericalHarmonics (segmentedBrainData, B, coefficientsb1000, bvals[index], 1000)
elif bvals[index] > 1500 and bvals[index] < 2500:
attenuatedBrainData = pl.attenuateImageSphericalHarmonics (segmentedBrainData, B, coefficientsb2000, bvals[index], 2000)
attenuatedBrainNii = nib.Nifti1Image(attenuatedBrainData, segmentedBrain.affine,segmentedBrain.header)
attenuatedBrainNii.to_filename(os.path.join( outputDirDirection,"brain.nii.gz"))
#Register to reference brain to get sizes right
print('Registering volume ' + str(index))
call(["flirt","-in",outputDirDirection+ "/brain.nii.gz","-ref",outputDir+ "/brainref.nii.gz","-applyxfm","-out",outputDirDirection+ "/brain.nii.gz"])
#Apply motion to brain here
if motionAndEddyImages == True:
outputDirDirectionMotionAndEddy = outputDir+"/DirectionMotionAndEddy"+str(index)
shutil.copytree(outputDirDirection,outputDirDirectionMotionAndEddy)
if motionDir is not "None":
shutil.copy(motionDir + "/motion" + str(index) + '.txt', outputDirDirectionMotionAndEddy+ "/motion")
#Read in pulse
if index == 0:
pulse=pl.read_pulse(outputDirDirection+"/pulse")
pulseinfo = np.loadtxt(outputDirDirection+'/pulse.info')
#Make distorted eddy pulse
if eddyGradients=='flat':
pl.generateEddyPulseFromBvecFlat(simDir,codeDir,args.matlabDir,basicSettings,ep, tau,bvals[index], bvec)
else:
new_pulse = pl.addEddyAccordingToBvec(pulse,pulseinfo,basicSettings[0],basicSettings[1],basicSettings[2],basicSettings[3],ep, tau,bvals[index], bvec)
#Interleave
if (interleaveFactor != 1):
new_pulse = pl.interleavePulse(new_pulse,int(pulseinfo[12]),interleaveFactor)
#Save to correct location
pl.write_pulse(outputDirDirectionMotionAndEddy+"/pulse",new_pulse)
#Save copy of distorted pulse
shutil.copy(outputDirDirectionMotionAndEddy+"/pulse", outputDir+"/Distortions/Eddy/pulseEddy"+str(index))
if normalImages == False:
shutil.rmtree(outputDirDirection)