python/structen.py

# -*- coding: utf-8 -*-
"""
Created on Sun Mar 12 21:54:40 2017
@author: david
"""
import numpy as np
import scipy as sp
import scipy.ndimage
import progressbar
import glob
def st2(I, s=1, dtype=np.float32):   
    
    #calculates the 2D structure tensor for an image using a gaussian window with standard deviation s
    
    #calculate the gradient
    dI = np.gradient(I.astype(dtype))
      
    #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 = np.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
    
    #if the user does not want a blurred field
    if(s == 0):
        return Tg
        
    #blur the tensor field
    T = np.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, 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
    #check the format for the window size
        
    v = np.array(v)
    print("\nCalculating gradient...")
    dI = np.gradient(I.astype(dtype), v[0], v[1], v[2])
    #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 = np.zeros(field_dims, dtype=np.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 = np.zeros(field_dims, dtype=np.float32)
    
    print("\nConvolving tensor field...")
    bar = progressbar.ProgressBar()
    bar.max_value = 6
    sigma = s / v
    print(sigma)
    for i in range(6):
        T[:, :, :, i] = scipy.ndimage.filters.gaussian_filter(Tg[:, :, :, i], sigma)
        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 * (np.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 = np.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 vec(S, vector=0):
   
    if(S.ndim != 3):
        print("ERROR: a 2D slice is expected")
        return
      
    #convert the field to a full rank-2 tensor
    T = sym2mat(S);
    del(S)
    
    #calculate the eigenvectors and eigenvalues
    l, v = np.linalg.eig(T)
    
    #get the dimension of the tensor field
    d = T.shape[2]
    
    #allocate space for the vector field
    V = np.zeros([T.shape[0], T.shape[1], 3])
    #arrange the indices for each pixel from smallest to largest eigenvalue    
    idx = l.argsort()
    
    for di in range(d):
        b = idx[:, :, -1-vector] == 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 = np.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
#calculate the anisotropy of a structure tensor given the tensor field S
def anisotropy3(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]
    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 = np.sqrt(3./2.) * np.sqrt(fa_num) / np.sqrt(fa_den)
    
    return FA, Cl, Cp, Cs
#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
#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):
    #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
    A = np.copy(T)
    A[..., 3] = T[..., 2]
    A[..., 2] = T[..., 3]
    
    #roll the tensor axis so that it is the leading component
    #A = numpy.rollaxis(A, A.ndim - 1)
    A.tofile(filename)
    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])
    s = np.array(s)
    
    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
        C = np.absolute(V)                       #calculate the absolute value
        
        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
            A = np.ones(T.shape[0], T.shape[1])
        image = C * np.expand_dims(A, 3)
        
        filename = prefix + str(z).zfill(4) + ".bmp"
        scipy.misc.imsave(filename, image)
        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