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
20
21
|
dI = np.gradient(I.astype(dtype))
#print(dI[1][20, 164])
|
7ce7cd7b
David Mayerich
added structure t...
|
22
23
24
25
26
|
#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 ...
|
27
|
Tg = np.zeros(field_dims, dtype=dtype)
|
7ce7cd7b
David Mayerich
added structure t...
|
28
29
30
31
32
33
34
35
36
|
#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
#blur the tensor field
|
08ada8c2
David Mayerich
updated image to ...
|
37
|
T = np.zeros(field_dims, dtype=dtype)
|
7ce7cd7b
David Mayerich
added structure t...
|
38
39
40
41
42
43
44
|
for i in range(3):
T[:, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, i], [s, s])
return T
|
08ada8c2
David Mayerich
updated image to ...
|
45
46
|
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...
|
47
|
#check the format for the window size
|
7ce7cd7b
David Mayerich
added structure t...
|
48
|
|
08ada8c2
David Mayerich
updated image to ...
|
49
|
v = np.array(v)
|
7ce7cd7b
David Mayerich
added structure t...
|
50
|
print("\nCalculating gradient...")
|
08ada8c2
David Mayerich
updated image to ...
|
51
|
dI = np.gradient(I.astype(dtype), v[0], v[1], v[2])
|
7ce7cd7b
David Mayerich
added structure t...
|
52
53
54
55
|
#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 ...
|
56
|
Tg = np.zeros(field_dims, dtype=np.float32)
|
7ce7cd7b
David Mayerich
added structure t...
|
57
58
59
60
61
62
63
64
65
66
67
68
69
|
#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 ...
|
70
|
T = np.zeros(field_dims, dtype=np.float32)
|
7ce7cd7b
David Mayerich
added structure t...
|
71
72
73
74
|
print("\nConvolving tensor field...")
bar = progressbar.ProgressBar()
bar.max_value = 6
|
08ada8c2
David Mayerich
updated image to ...
|
75
76
|
sigma = s / v
print(sigma)
|
7ce7cd7b
David Mayerich
added structure t...
|
77
|
for i in range(6):
|
08ada8c2
David Mayerich
updated image to ...
|
78
|
T[:, :, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, :, i], sigma)
|
7ce7cd7b
David Mayerich
added structure t...
|
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
|
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 ...
|
110
|
d = int(0.5 * (np.sqrt(8. * n + 1.) - 1.))
|
7ce7cd7b
David Mayerich
added structure t...
|
111
112
113
114
115
|
#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 ...
|
116
|
R = np.zeros(out_s)
|
7ce7cd7b
David Mayerich
added structure t...
|
117
118
119
120
121
122
123
124
125
126
127
|
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
|
08ada8c2
David Mayerich
updated image to ...
|
128
129
|
def st2vec(S, vector=0):
|
29a37cd9
David Mayerich
fixed matrix_sym ...
|
130
131
132
|
if(S.ndim != 3):
print("ERROR: a 2D slice is expected")
return
|
7ce7cd7b
David Mayerich
added structure t...
|
133
134
135
136
137
138
|
#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 ...
|
139
|
l, v = np.linalg.eig(T)
|
7ce7cd7b
David Mayerich
added structure t...
|
140
141
142
143
144
|
#get the dimension of the tensor field
d = T.shape[2]
#allocate space for the vector field
|
08ada8c2
David Mayerich
updated image to ...
|
145
146
147
|
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...
|
148
149
150
|
idx = l.argsort()
for di in range(d):
|
08ada8c2
David Mayerich
updated image to ...
|
151
|
b = idx[:, :, -1-vector] == di
|
7ce7cd7b
David Mayerich
added structure t...
|
152
153
154
|
V[b, 0:d] = v[b, :, di]
return V
|
08ada8c2
David Mayerich
updated image to ...
|
155
|
|
7ce7cd7b
David Mayerich
added structure t...
|
156
157
158
159
160
161
162
163
164
165
166
167
168
|
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 ...
|
169
|
M = np.zeros([X, Y, Z]).astype('float32')
|
7ce7cd7b
David Mayerich
added structure t...
|
170
171
172
173
174
175
176
177
178
179
180
181
|
#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 ...
|
182
183
|
#calculate the anisotropy of a structure tensor given the tensor field S
def anisotropy3(S):
|
7ce7cd7b
David Mayerich
added structure t...
|
184
185
186
187
|
Sf = sym2mat(S)
#calculate the eigenvectors and eigenvalues
|
08ada8c2
David Mayerich
updated image to ...
|
188
|
l, v = np.linalg.eig(Sf)
|
7ce7cd7b
David Mayerich
added structure t...
|
189
190
|
#store the sorted eigenvalues
|
08ada8c2
David Mayerich
updated image to ...
|
191
|
ls = np.sort(l)
|
7ce7cd7b
David Mayerich
added structure t...
|
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
|
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 ...
|
209
|
FA = np.sqrt(3./2.) * np.sqrt(fa_num) / np.sqrt(fa_den)
|
7ce7cd7b
David Mayerich
added structure t...
|
210
211
212
|
return FA, Cl, Cp, Cs
|
08ada8c2
David Mayerich
updated image to ...
|
213
214
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
|
#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
|
7ce7cd7b
David Mayerich
added structure t...
|
241
242
243
244
245
246
247
248
249
250
251
|
def st2amira(filename, T):
#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 ...
|
252
|
A = np.copy(T)
|
95f1e985
David Mayerich
updated structure...
|
253
254
|
A[..., 3] = T[..., 2]
A[..., 2] = T[..., 3]
|
7ce7cd7b
David Mayerich
added structure t...
|
255
256
|
#roll the tensor axis so that it is the leading component
|
14379866
David Mayerich
updated amira out...
|
257
|
#A = numpy.rollaxis(A, A.ndim - 1)
|
7ce7cd7b
David Mayerich
added structure t...
|
258
|
A.tofile(filename)
|
95f1e985
David Mayerich
updated structure...
|
259
260
261
262
263
264
265
266
267
268
269
270
271
272
|
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 ...
|
273
|
s = np.array(s)
|
95f1e985
David Mayerich
updated structure...
|
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
|
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 ...
|
297
|
C = np.absolute(V) #calculate the absolute value
|
95f1e985
David Mayerich
updated structure...
|
298
299
300
301
302
303
304
305
|
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 ...
|
306
307
|
A = np.ones(T.shape[0], T.shape[1])
image = C * np.expand_dims(A, 3)
|
95f1e985
David Mayerich
updated structure...
|
308
309
310
|
filename = prefix + str(z).zfill(4) + ".bmp"
scipy.misc.imsave(filename, image)
|
08ada8c2
David Mayerich
updated image to ...
|
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
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
|
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
|