From a69fee9dc89a6dab5b4bf4a28b23d945e8182504 Mon Sep 17 00:00:00 2001
From: David Mayerich <mayerich@uh.edu>
Date: Fri, 19 May 2017 13:48:13 -0500
Subject: [PATCH] added python script for 3D phantom creation for 2D images

---
 python/phantom.py  | 58 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 python/structen.py | 25 ++++++++++++++++++++-----
 2 files changed, 78 insertions(+), 5 deletions(-)
 create mode 100644 python/phantom.py

diff --git a/python/phantom.py b/python/phantom.py
new file mode 100644
index 0000000..87e2800
--- /dev/null
+++ b/python/phantom.py
@@ -0,0 +1,58 @@
+# -*- 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])
+
+
+
+
diff --git a/python/structen.py b/python/structen.py
index 672ee5e..cc96fca 100644
--- a/python/structen.py
+++ b/python/structen.py
@@ -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
--
libgit2 0.21.4