-
Notifications
You must be signed in to change notification settings - Fork 0
/
read_pso_mammoth.py
100 lines (76 loc) · 3.02 KB
/
read_pso_mammoth.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
import nibabel as nib
import numpy as np
from dipy.sims.voxel import MultiTensor
from load_data import get_train_dti
mask = nib.load('wm_mask_hardi_01.nii.gz').get_data()
_, affine, gtab = get_train_dti()
def parse_params(params, NC, iso, fr):
mevals = []
angles = []
fractions = []
for i in range(NC):
mevals.append(np.array([params[4 * i], params[4 * i + 1], params[4 * i + 1]]))
angles.append(np.array([params[4 * i + 2], params[4 * i + 3]]))
if iso:
mevals.append(np.array([params[4 * NC], params[4 * NC], params[4 * NC]]))
angles.append(np.array([0., 0.]))
iso_frac = params[4 * NC + 1]
if ~fr:
if ~iso:
fractions.append(np.ones((1, NC)) * (100. / NC))
else:
fractions.append(np.vstack((np.ones((NC, 1)) * (1. / NC) * (100. - iso_frac), np.array([iso_frac]))).T)
else:
if NC == 1:
frac = 100.
else:
fracc = params[-(NC - 1)]
frac = np.zeros(NC)
frac[0] = fracc[0]
if NC == 2:
frac[1] = 100. - frac[0]
else:
frac[1] = (100. - frac[0]) * fracc[1] / 100.
if NC == 3:
frac[2] = 100. - (frac[0] + frac[1])
if ~iso:
fractions.append(frac)
else:
fractions.append(np.vstack((frac * (100. - iso_frac) / 100., np.array([iso_frac]))).T)
return np.array(mevals), np.array(angles), np.array(fractions).T
f = open('out_pso_dtiPier_hardi_hardiPier', 'r')
# f = open('out__pso_testDataSlice__23', 'r')
# f = open('out__pso_testDataSlice__23', 'r')
line = f.readline()
while(line!=''):
# NC iso fr Npart Niter den
params = np.array(line[:-1].split(' '), dtype=np.int)
sliceXY = np.zeros((50,50,params[0],3))
#first line of data
line = f.readline()
dat = np.array(line[:-1].split(' '),dtype=np.float64)
#number of voxel in slice
slicez = int(dat[2])
nb_to_read = int(mask[:,:,slicez].sum())
#parse first voxel
mevals, angles, fractions = parse_params(dat[3:],params[0], params[1], params[2])
# could keep fractions also
_, sticks = MultiTensor(gtab, mevals, 100, angles, fractions, None)
if params[1]:
sticks = sticks[:-1]
#save sticks
sliceXY[dat[0],dat[1]] = sticks
# read rest of slice
for i in range(nb_to_read-1):
line = f.readline()
dat = np.array(line[:-1].split(' '),dtype=np.float64)
mevals, angles, fractions = parse_params(dat[3:],params[0], params[1], params[2])
_, sticks = MultiTensor(gtab, mevals, 100, angles, fractions, None)
if params[1]:
sticks = sticks[:-1]
sliceXY[dat[0],dat[1]] = sticks
filename = 'sticks_NC={}_iso={}_fr={}_Np={}_Ni={}_snr={}_typ={}_slic={}'.format(params[0],params[1],params[2],params[3],params[4],params[5],params[6],slicez)
nib.save(nib.Nifti1Image(sliceXY, affine), '/media/Data/work/isbi2013/slices_out_pso/' + filename + '.nii.gz')
print('saved! NC = {} iso = {} fr = {} snr = {} typ = {} slic = {}'.format(params[0],params[1],params[2],params[5],params[6],slicez))
line = f.readline()
f.close()