added python script for 3D phantom creation for 2D images

David Mayerich
1 parent 08ada8c2
Showing 2 changed files with 78 additions and 5 deletions Show diff stats
python/phantom.py
python/structen.py
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May 18 15:56:20 2017
+
+@author: david
+"""
+
+import numpy as np
+import scipy as sp
+import scipy.misc
+import scipy.ndimage
+import matplotlib.pyplot as plt
+import os
+
+#calculate a 3D image from a 2D binary mask
+def binary(infile, sigma=2):
+    
+    I = sp.misc.imread(infile).astype(np.bool)
+    
+    #if the image has more than one channel, just keep the first one
+    if(I.ndim == 3):
+        I = I[:, :, 0]
+    
+    L = []
+    while np.count_nonzero(I) != 0:
+        L.append(I)
+        I = sp.ndimage.binary_erosion(I)
+        
+    #create a 3D image representing the new stack    
+    S = np.zeros( (I.shape[0], I.shape[1], len(L) * 2 - 1) )
+    
+    #for each image in the list
+    for i in range(0, len(L)):
+        if(i == 0):
+            S[:, :, len(L) - 1] = L[0]
+        else:
+            S[:, :, len(L) - 1 + i] = L[i]
+            S[:, :, len(L) - 1 - i] = L[i]
+            
+    S = sp.ndimage.filters.gaussian_filter(S, sigma)
+    return S
+
+#generate a 3D image stack from a 2D binary mask
+def binary_stack(infile, outdir, sigma=2):
+    outfile_base = os.path.basename(infile)
+    outfile_prefix, outfile_ext = os.path.splitext(outfile_base)
+    
+    S = binary(infile, sigma)
+
+    zcount = len(str(S.shape[2]))
+    for f in range(0, S.shape[2]):
+        fname = outdir + "/" + outfile_prefix + str(f).zfill(zcount) + outfile_ext
+    
+        sp.misc.imsave(fname, S[:, :, f])
+
+
+
+
@@ -17,9 +17,7 @@ def st2(I, s=1, dtype=np.float32):
  
     #calculate the gradient
     dI = np.gradient(I.astype(dtype))
-    
-    #print(dI[1][20, 164])
-    
+      
     #calculate the dimensions of the tensor field
     field_dims = [dI[0].shape[0], dI[0].shape[1], 3]
  
@@ -33,6 +31,10 @@ def st2(I, s=1, dtype=np.float32):
             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)
  
@@ -125,7 +127,7 @@ def sym2mat(T):
  
     return R   
  
-def st2vec(S, vector=0):
+def vec(S, vector=0):
  
     if(S.ndim != 3):
         print("ERROR: a 2D slice is expected")
@@ -238,7 +240,20 @@ def fa(S):
         FA = np.sqrt(3./2.) * np.sqrt(fa_num) / np.sqrt(fa_den)
         return FA
  
-def st2amira(filename, T):
+#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