diff --git a/python/structen.py b/python/structen.py new file mode 100644 index 0000000..80728fc --- /dev/null +++ b/python/structen.py @@ -0,0 +1,232 @@ +# -*- coding: utf-8 -*- +""" +Created on Sun Mar 12 21:54:40 2017 + +@author: david +""" + +import numpy +import scipy.ndimage +import progressbar +import glob + +def st2(I, s=1, dtype=numpy.float32): + + + #calculates the 2D structure tensor for an image using a gaussian window with standard deviation s + + #calculate the gradient + dI = numpy.gradient(I) + + #calculate the dimensions of the tensor field + field_dims = [dI[0].shape[0], dI[0].shape[1], 3] + + #allocate space for the tensor field + Tg = numpy.zeros(field_dims, dtype=dtype) + + #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 + T = numpy.zeros(field_dims, dtype=dtype) + + for i in range(3): + T[:, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, i], [s, s]) + + + return T + +def st3(I, s=1): + #calculate the structure tensor field for the 3D input image I given the window size s in 3D + #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]) + + print("\nCalculating gradient...") + dI = numpy.gradient(I) + #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 + Tg = numpy.zeros(field_dims, dtype=numpy.float32) + + #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 + T = numpy.zeros(field_dims, dtype=numpy.float32) + + print("\nConvolving tensor field...") + bar = progressbar.ProgressBar() + bar.max_value = 6 + for i in range(6): + T[:, :, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, :, i], s) + 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 + d = int(0.5 * (numpy.sqrt(8. * n + 1.) - 1.)) + + #calculate the dimensions of the output field + out_s = list(in_s)[:-1] + [d] + [d] + + #allocate space for the output field + R = numpy.zeros(out_s) + + 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 + +def st2vec(S, vector='largest'): + #Create a color image from a 2D or 3D structure tensor slice + + #convert the field to a full rank-2 tensor + T = sym2mat(S); + del(S) + + #calculate the eigenvectors and eigenvalues + l, v = numpy.linalg.eig(T) + + #get the dimension of the tensor field + d = T.shape[2] + + #allocate space for the vector field + V = numpy.zeros([T.shape[0], T.shape[1], 3]) + + idx = l.argsort() + + for di in range(d): + if vector == 'smallest': + b = idx[:, :, 0] == di + elif vector == 'largest': + b = idx[:, :, d-1] == di + else: + b = idx[:, :, 1] == di + V[b, 0:d] = v[b, :, di] + + return V + +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 + M = numpy.zeros([X, Y, Z]).astype('float32') + + #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 + +def anisotropy(S): + + Sf = sym2mat(S) + + #calculate the eigenvectors and eigenvalues + l, v = numpy.linalg.eig(Sf) + + #store the sorted eigenvalues + ls = numpy.sort(l) + 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 + FA = numpy.sqrt(3./2.) * numpy.sqrt(fa_num) / numpy.sqrt(fa_den) + + return FA, Cl, Cp, Cs + +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 + temp = T[..., 3] + T[..., 3] = T[..., 2] + T[..., 2] = temp; + + #roll the tensor axis so that it is the leading component + A = numpy.rollaxis(T, T.ndim - 1) + A.tofile(filename) + print(A.shape) \ No newline at end of file -- libgit2 0.21.4