updated structure tensor functions

David Mayerich
1 parent 7ce7cd7b
Showing 2 changed files with 60 additions and 7 deletions Show diff stats
python/structen.py
stim/parser/filename.h
@@ -10,8 +10,7 @@ import scipy.ndimage
 import progressbar
 import glob
  
-def st2(I, s=1, dtype=numpy.float32):
-    
+def st2(I, s=1, dtype=numpy.float32):   
  
     #calculates the 2D structure tensor for an image using a gaussian window with standard deviation s
  
@@ -222,11 +221,63 @@ def st2amira(filename, T):
     #   5    dz dz   ----> 5
  
     #swap the 2nd and 3rd tensor components
-    temp = T[..., 3]
-    T[..., 3] = T[..., 2]
-    T[..., 2] = temp;
+    A = numpy.copy(T)
+    A[..., 3] = T[..., 2]
+    A[..., 2] = T[..., 3]
  
     #roll the tensor axis so that it is the leading component
-    A = numpy.rollaxis(T, T.ndim - 1)
+    A = numpy.rollaxis(A, A.ndim - 1)
     A.tofile(filename)
-    print(A.shape)
 \ No newline at end of file
+    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 = numpy.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 = numpy.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 = numpy.ones(T.shape[0], T.shape[1])
+        image = C * numpy.expand_dims(A, 3)
+        
+        filename = prefix + str(z).zfill(4) + ".bmp"
+        scipy.misc.imsave(filename, image)
+        bar.update(z + 1)
 \ No newline at end of file
@@ -89,6 +89,8 @@ protected:
 				absolute.push_back(relative[i]);
 			}
 		}
+		else
+			absolute = relative;
 	}
  
 	/// Parses a directory string into a drive (NULL if not Windows) and list of hierarchical directories