7ce7cd7b
David Mayerich
added structure t...
|
1
2
3
4
5
6
7
|
# -*- coding: utf-8 -*-
"""
Created on Sun Mar 12 21:54:40 2017
@author: david
"""
|
08ada8c2
David Mayerich
updated image to ...
|
8
9
|
import numpy as np
import scipy as sp
|
7ce7cd7b
David Mayerich
added structure t...
|
10
11
12
13
|
import scipy.ndimage
import progressbar
import glob
|
08ada8c2
David Mayerich
updated image to ...
|
14
|
def st2(I, s=1, dtype=np.float32):
|
7ce7cd7b
David Mayerich
added structure t...
|
15
16
17
18
|
#calculates the 2D structure tensor for an image using a gaussian window with standard deviation s
#calculate the gradient
|
08ada8c2
David Mayerich
updated image to ...
|
19
|
dI = np.gradient(I.astype(dtype))
|
a69fee9d
David Mayerich
added python scri...
|
20
|
|
7ce7cd7b
David Mayerich
added structure t...
|
21
22
23
24
|
#calculate the dimensions of the tensor field
field_dims = [dI[0].shape[0], dI[0].shape[1], 3]
#allocate space for the tensor field
|
08ada8c2
David Mayerich
updated image to ...
|
25
|
Tg = np.zeros(field_dims, dtype=dtype)
|
7ce7cd7b
David Mayerich
added structure t...
|
26
27
28
29
30
31
32
33
|
#calculate the gradient components of the tensor
ti = 0
for i in range(2):
for j in range(i + 1):
Tg[:, :, ti] = dI[j] * dI[i]
ti = ti + 1
|
a69fee9d
David Mayerich
added python scri...
|
34
35
36
37
|
#if the user does not want a blurred field
if(s == 0):
return Tg
|
7ce7cd7b
David Mayerich
added structure t...
|
38
|
#blur the tensor field
|
08ada8c2
David Mayerich
updated image to ...
|
39
|
T = np.zeros(field_dims, dtype=dtype)
|
7ce7cd7b
David Mayerich
added structure t...
|
40
41
42
43
44
45
46
|
for i in range(3):
T[:, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, i], [s, s])
return T
|
08ada8c2
David Mayerich
updated image to ...
|
47
48
|
def st3(I, s=1, v=[1, 1, 1], dtype=np.float32):
#calculate the structure tensor field for the 3D input image I given the window size s and voxel size v
|
7ce7cd7b
David Mayerich
added structure t...
|
49
|
#check the format for the window size
|
7ce7cd7b
David Mayerich
added structure t...
|
50
|
|
08ada8c2
David Mayerich
updated image to ...
|
51
|
v = np.array(v)
|
7ce7cd7b
David Mayerich
added structure t...
|
52
|
print("\nCalculating gradient...")
|
08ada8c2
David Mayerich
updated image to ...
|
53
|
dI = np.gradient(I.astype(dtype), v[0], v[1], v[2])
|
7ce7cd7b
David Mayerich
added structure t...
|
54
55
56
57
|
#calculate the dimensions of the tensor field
field_dims = [dI[0].shape[0], dI[0].shape[1], dI[0].shape[2], 6]
#allocate space for the tensor field
|
08ada8c2
David Mayerich
updated image to ...
|
58
|
Tg = np.zeros(field_dims, dtype=np.float32)
|
7ce7cd7b
David Mayerich
added structure t...
|
59
60
61
62
63
64
65
66
67
68
69
70
71
|
#calculate the gradient components of the tensor
ti = 0
print("Calculating tensor components...")
bar = progressbar.ProgressBar()
bar.max_value = 6
for i in range(3):
for j in range(i + 1):
Tg[:, :, :, ti] = dI[j] * dI[i]
ti = ti + 1
bar.update(ti)
#blur the tensor field
|
08ada8c2
David Mayerich
updated image to ...
|
72
|
T = np.zeros(field_dims, dtype=np.float32)
|
7ce7cd7b
David Mayerich
added structure t...
|
73
74
75
76
|
print("\nConvolving tensor field...")
bar = progressbar.ProgressBar()
bar.max_value = 6
|
08ada8c2
David Mayerich
updated image to ...
|
77
78
|
sigma = s / v
print(sigma)
|
7ce7cd7b
David Mayerich
added structure t...
|
79
|
for i in range(6):
|
08ada8c2
David Mayerich
updated image to ...
|
80
|
T[:, :, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, :, i], sigma)
|
7ce7cd7b
David Mayerich
added structure t...
|
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
|
bar.update(i+1)
return T
def st(I, s=1):
if I.ndim == 3:
return st3(I, s)
elif I.ndim == 2:
return st2(I, s)
else:
print("Image must be 2D or 3D")
return
def sym2mat(T):
#Calculate the full symmetric matrix from a list of lower triangular elements.
#The lower triangular components in the input field T are assumed to be the
# final dimension of the input matrix.
# | 1 2 4 7 |
# | 0 3 5 8 |
# | 0 0 6 9 |
# | 0 0 0 10 |
in_s = T.shape
#get the number of tensor elements
n = in_s[T.ndim - 1]
#calculate the dimension of the symmetric matrix
|
08ada8c2
David Mayerich
updated image to ...
|
112
|
d = int(0.5 * (np.sqrt(8. * n + 1.) - 1.))
|
7ce7cd7b
David Mayerich
added structure t...
|
113
114
115
116
117
|
#calculate the dimensions of the output field
out_s = list(in_s)[:-1] + [d] + [d]
#allocate space for the output field
|
08ada8c2
David Mayerich
updated image to ...
|
118
|
R = np.zeros(out_s)
|
7ce7cd7b
David Mayerich
added structure t...
|
119
120
121
122
123
124
125
126
127
128
129
|
ni = 0
for i in range(d):
for j in range(i + 1):
R[..., i, j] = T[..., ni]
if i != j:
R[..., j, i] = T[..., ni]
ni = ni + 1
return R
|
a69fee9d
David Mayerich
added python scri...
|
130
|
def vec(S, vector=0):
|
08ada8c2
David Mayerich
updated image to ...
|
131
|
|
29a37cd9
David Mayerich
fixed matrix_sym ...
|
132
133
134
|
if(S.ndim != 3):
print("ERROR: a 2D slice is expected")
return
|
7ce7cd7b
David Mayerich
added structure t...
|
135
136
137
138
139
140
|
#convert the field to a full rank-2 tensor
T = sym2mat(S);
del(S)
#calculate the eigenvectors and eigenvalues
|
08ada8c2
David Mayerich
updated image to ...
|
141
|
l, v = np.linalg.eig(T)
|
7ce7cd7b
David Mayerich
added structure t...
|
142
143
144
145
146
|
#get the dimension of the tensor field
d = T.shape[2]
#allocate space for the vector field
|
08ada8c2
David Mayerich
updated image to ...
|
147
148
149
|
V = np.zeros([T.shape[0], T.shape[1], 3])
#arrange the indices for each pixel from smallest to largest eigenvalue
|
7ce7cd7b
David Mayerich
added structure t...
|
150
151
152
|
idx = l.argsort()
for di in range(d):
|
08ada8c2
David Mayerich
updated image to ...
|
153
|
b = idx[:, :, -1-vector] == di
|
7ce7cd7b
David Mayerich
added structure t...
|
154
155
156
|
V[b, 0:d] = v[b, :, di]
return V
|
08ada8c2
David Mayerich
updated image to ...
|
157
|
|
7ce7cd7b
David Mayerich
added structure t...
|
158
159
160
161
162
163
164
165
166
167
168
169
170
|
def loadstack(filemask):
#Load an image stack as a 3D grayscale array
#get a list of all files matching the given mask
files = [file for file in glob.glob(filemask)]
#calculate the size of the output stack
I = scipy.misc.imread(files[0])
X = I.shape[0]
Y = I.shape[1]
Z = len(files)
#allocate space for the image stack
|
08ada8c2
David Mayerich
updated image to ...
|
171
|
M = np.zeros([X, Y, Z]).astype('float32')
|
7ce7cd7b
David Mayerich
added structure t...
|
172
173
174
175
176
177
178
179
180
181
182
183
|
#create a progress bar
bar = progressbar.ProgressBar()
bar.max_value = Z
#for each file
for z in range(Z):
#load the file and save it to the image stack
M[:, :, z] = scipy.misc.imread(files[z], flatten="True").astype('float32')
bar.update(z+1)
return M
|
08ada8c2
David Mayerich
updated image to ...
|
184
185
|
#calculate the anisotropy of a structure tensor given the tensor field S
def anisotropy3(S):
|
7ce7cd7b
David Mayerich
added structure t...
|
186
187
188
189
|
Sf = sym2mat(S)
#calculate the eigenvectors and eigenvalues
|
08ada8c2
David Mayerich
updated image to ...
|
190
|
l, v = np.linalg.eig(Sf)
|
7ce7cd7b
David Mayerich
added structure t...
|
191
192
|
#store the sorted eigenvalues
|
08ada8c2
David Mayerich
updated image to ...
|
193
|
ls = np.sort(l)
|
7ce7cd7b
David Mayerich
added structure t...
|
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
l0 = ls[:, :, 0]
l1 = ls[:, :, 1]
l2 = ls[:, :, 2]
#calculate the linear anisotropy
Cl = (l2 - l1)/(l2 + l1 + l0)
#calculate the planar anisotropy
Cp = 2 * (l1 - l0) / (l2 + l1 + l0)
#calculate the spherical anisotropy
Cs = 3 * l0 / (l2 + l1 + l0)
#calculate the fractional anisotropy
l_hat = (l0 + l1 + l2)/3
fa_num = (l2 - l_hat) ** 2 + (l1 - l_hat) ** 2 + (l0 - l_hat) ** 2;
fa_den = l0 ** 2 + l1 ** 2 + l2 ** 2
|
08ada8c2
David Mayerich
updated image to ...
|
211
|
FA = np.sqrt(3./2.) * np.sqrt(fa_num) / np.sqrt(fa_den)
|
7ce7cd7b
David Mayerich
added structure t...
|
212
213
214
|
return FA, Cl, Cp, Cs
|
08ada8c2
David Mayerich
updated image to ...
|
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
|
#calculate the fractional anisotropy
def fa(S):
Sf = sym2mat(S)
#calculate the eigenvectors and eigenvalues
l, v = np.linalg.eig(Sf)
#store the sorted eigenvalues
ls = np.sort(l)
l0 = ls[:, :, 0]
l1 = ls[:, :, 1]
#if this is a 2D tensor, calculate and return the coherence
if(S.shape[2] == 3):
C = ((l0 - l1) / (l0 + l1)) ** 2
return C
#if this is a 3D tensor
elif(S.shape[2] == 6):
l2 = ls[:, :, 2]
#calculate the fractional anisotropy
l_hat = (l0 + l1 + l2)/3
fa_num = (l2 - l_hat) ** 2 + (l1 - l_hat) ** 2 + (l0 - l_hat) ** 2;
fa_den = l0 ** 2 + l1 ** 2 + l2 ** 2
FA = np.sqrt(3./2.) * np.sqrt(fa_num) / np.sqrt(fa_den)
return FA
|
a69fee9d
David Mayerich
added python scri...
|
243
244
245
246
247
248
249
250
251
252
253
254
255
256
|
#calculate the specified eigenvalue for the tensor field
def eigenval(S, ev):
Sf = sym2mat(S)
#calculate the eigenvectors and eigenvalues
l, v = np.linalg.eig(Sf)
#store the sorted eigenvalues
ls = np.sort(l)
evals = ls[:, :, ev]
return evals
def amira(filename, T):
|
7ce7cd7b
David Mayerich
added structure t...
|
257
258
259
260
261
262
263
264
265
266
|
#generates a tensor field that can be imported into Amira
# 0 dx dx ----> 0
# 1 dx dy ----> 1
# 2 dy dy ----> 3
# 3 dx dz ----> 2
# 4 dy dz ----> 4
# 5 dz dz ----> 5
#swap the 2nd and 3rd tensor components
|
08ada8c2
David Mayerich
updated image to ...
|
267
|
A = np.copy(T)
|
95f1e985
David Mayerich
updated structure...
|
268
269
|
A[..., 3] = T[..., 2]
A[..., 2] = T[..., 3]
|
7ce7cd7b
David Mayerich
added structure t...
|
270
271
|
#roll the tensor axis so that it is the leading component
|
14379866
David Mayerich
updated amira out...
|
272
|
#A = numpy.rollaxis(A, A.ndim - 1)
|
7ce7cd7b
David Mayerich
added structure t...
|
273
|
A.tofile(filename)
|
95f1e985
David Mayerich
updated structure...
|
274
275
276
277
278
279
280
281
282
283
284
285
286
287
|
print("\n", A.shape)
def resample3(T, s=2):
#resample a tensor field by an integer factor s
#This function first convolves the field with a box filter and then
# re-samples to create a smaller field
#check the format for the window size
if type(s) is not list:
s = [s] * 3
elif len(s) == 1:
s = s * 3
elif len(s) == 2:
s.insert(1, s[0])
|
08ada8c2
David Mayerich
updated image to ...
|
288
|
s = np.array(s)
|
95f1e985
David Mayerich
updated structure...
|
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
|
bar = progressbar.ProgressBar()
bar.max_value = T.shape[3]
#blur with a uniform box filter of size r
for t in range(T.shape[3]):
T[..., t] = scipy.ndimage.filters.uniform_filter(T[..., t], 2 * s)
bar.update(t+1)
#resample at a rate of r
R = T[::s[0], ::s[1], ::s[2], :]
return R
def color3(prefix, T, vector='largest', aniso=True):
#Saves a stack of color images corresponding to the eigenvector and optionally scaled by anisotropy
bar = progressbar.ProgressBar()
bar.max_value = T.shape[2]
#for each z-axis slice
for z in range(T.shape[2]):
S = T[:, :, z, :] #get the slice
V = st2vec(S, vector='smallest') #calculate the vector
|
08ada8c2
David Mayerich
updated image to ...
|
312
|
C = np.absolute(V) #calculate the absolute value
|
95f1e985
David Mayerich
updated structure...
|
313
314
315
316
317
318
319
320
|
if aniso == True: #if the image is scaled by anisotropy
FA, Cl, Cp, Cs = anisotropy(S) #calculate the anisotropy of the slice
if vector == 'largest':
A = Cl
elif vector == 'smallest':
A = Cp
else: #otherwise just scale by 1
|
08ada8c2
David Mayerich
updated image to ...
|
321
322
|
A = np.ones(T.shape[0], T.shape[1])
image = C * np.expand_dims(A, 3)
|
95f1e985
David Mayerich
updated structure...
|
323
324
325
|
filename = prefix + str(z).zfill(4) + ".bmp"
scipy.misc.imsave(filename, image)
|
08ada8c2
David Mayerich
updated image to ...
|
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
|
bar.update(z + 1)
def st2stack(T, outfile, **kwargs):
eigenvector = False #initialize the colormap flags to false
aniso_color = False
aniso_alpha = False
#image = False
aniso_pwr = 1
cimage_pwr = 1
aimage_pwr = 1
anisostretch = 1 #set the contrast stretch parameter
alpha_channel = False
alpha_image = False
color_image = False
for k,v in kwargs.items(): #for each argument
if(k == "ev"): #if the user wants a colormap based on an eigenvector
eigenvector = True #set the eigenvector flag to true
ev = v #save the desired eigenvector
if(k == "aniso"): #if the user wants to factor in anisotropy
aniso = True #set the anisotropy flag to true
aniso_channel = v #save the anisotropy channel
if(k == "aniso_color"):
aniso_color = v
if(k == "aniso_alpha"):
aniso_alpha = v
if(k == "apwr"): #if the user wants to amplify the anisotropy
aniso_pwr = v #set the anisotropy exponent
if(k == "cipwr"): #if the user specifies the image power
cimage_pwr = v
if(k == "aipwr"):
aimage_pwr = v
if(k == "alphaimage"):
Ia = v
alpha_image = True
if(k == "colorimage"):
Ic = v
color_image = True
if(k == "anisostretch"):
anisostretch = v
if(k == "alpha"):
alpha_channel = v
bar = progressbar.ProgressBar()
bar.max_value = T.shape[2]
for i in range(0, T.shape[2]):
#for i in range(0, 50):
if(alpha_image or alpha_channel):
img = np.ones([T.shape[0], T.shape[1], 4])
else:
img = np.ones([T.shape[0], T.shape[1], 3])
if(eigenvector):
V = st2vec(T[:, :, i], ev) #get the vector field for slice i corresponding to eigenvector ev
img[:, :, 0:3] = V #update the image with the vector field information
if(aniso): #if the user is requesting anisotropy be incorporated into the image
FA, Cl, Cp, Cs = anisotropy(T[:, :, i]) #calculate the anisotropy of the tensors in slice i
if(aniso_channel == "fa"):
A = FA
elif(aniso_channel == "l"):
A = Cl
elif(aniso_channel == "p"):
A = Cp
else:
A = Cs
if(aniso_alpha):
print("rendering anisotropy to the alpha channel")
img[:, :, 3] = A ** aniso_pwr * anisostretch
if(aniso_color):
print("rendering anisotropy to the color channel")
img[:, :, 0:3] = img[:, :, 0:3] * np.expand_dims(A ** aniso_pwr, 3) * anisostretch
if(alpha_image):
img[:, :, 3] = Ia[:, :, i]/255 ** aimage_pwr
if(color_image):
img[:, :, 0:3] = img[:, :, 0:3] * (np.expand_dims(Ic[:, :, i], 3)/255) ** cimage_pwr #multiply the vector field by the image intensity
#outname = outfile + str(i).zfill(3) + ".bmp" #get the file name to be saved
outname = outfile.replace("*", str(i).zfill(3))
sp.misc.imsave(outname, np.ndarray.astype(np.abs(img)*255, "uint8")) #save the output image
bar.update(i+1)
#this function takes a 3D image and turns it into a stack of color images based on the structure tensor
def img2stack(I, outfile, **kwargs):
vs = [1, 1, 1] #set the default voxel size to 1
w = 5
for k,v in kwargs.items(): #for each argument
if(k == "voxelsize"): #if the voxel size is specified
if(len(v) == 1): #if the user just specifies one value
vs = [v] * 3 #assume that the voxels are isotropic, create a list of 3 v's
elif(len(v) == 2): #if the user specifies two values
vs[0] = v[0] #assume that the voxels are isotropic along (x, y) and anisotropic along z
vs[1] = v[0]
vs[2] = v[1]
elif(len(v) == 3):
vs = v
if(k == "window"): #if the user specifies a window size
w = v
T = st3(I, w, vs) #calculate the structure tensor
st2stack(T, outfile, **kwargs)
def stack2stack(infile_mask, outfile, **kwargs):
I = loadstack(infile_mask) #load the file mask
for k,v in kwargs.items(): #for each argument
if(k == "ipwr"):
img = I
img2stack(I, outfile, image=img, **kwargs) #call the function to convert the image to an output ST stack
|