Merge remote-tracking branch 'JIC-CSB/master'

rosscarter3 · Dec 2, 2015 · bc525d8 · bc525d8
2 parents e1b8ec9 + 729682b
commit bc525d8
Show file tree

Hide file tree

Showing 3 changed files with 207 additions and 83 deletions.
diff --git a/src/gaussproj.py b/src/gaussproj.py
@@ -1,19 +1,17 @@
 #!/usr/bin/env python2.7
 
+"""
+gaussproj - A simple script to project (single channel) three-dimensional
+image stacks onto two dimensional surfaces
+"""
+
 import itertools, sys, os
 import scipy
 import scipy.ndimage as nd
 import numpy as np
 import stackhandle
 import projpp
-
-def numpy_draw_pyplot(r):
-    plt.imshow(r, interpolation='nearest')
-    plt.show()
-
-def numpy_draw_pil(r):
-    i = Image.fromarray(r)
-    i.show()
+import proj
 
 def read_and_conv(filename):
     try:
@@ -32,13 +30,6 @@ def read_and_conv(filename):
 
     return np.amax(na, 2)
 
-def load_png_stack_pattern(impattern, istart, iend):
-    flush_message("Loading images...")
-    ifiles = [impattern % i for i in range(istart, iend)]
-    ma = np.dstack(itertools.imap(read_and_conv, ifiles))
-    print " done, array is", ma.shape
-    return ma
-
 def load_png_stack(ifiles):
     flush_message("Loading images...")
 
@@ -50,55 +41,14 @@ def flush_message(message):
     print message,
     sys.stdout.flush()
 
-def apply_gaussian_filter(ma, sd):
-    flush_message("Applying 3D gaussian filter...")
-    bl = nd.gaussian_filter(ma, sd)
-    print " done"
-    return bl
-
-def get_max(bl, x, y):
-    dp = list(bl[x, y, :])
-    return dp.index(max(dp))
-
-def find_projection_surface(bl):
-    flush_message("Finding projection surface...")
-    xmax, ymax, zmax = bl.shape
-    surface = np.zeros([xmax, ymax], dtype=np.uint8)
-    for x in range(0, xmax):
-        for y in range(0, ymax):
-            z = get_max(bl, x, y)
-            surface[x, y] = z
-    print " done"
-    return surface
-
-def vavg(ma, x, y, z):
-    return ma[x, y, z-3]
-    #return 0.25 * sum(ma[x, y, z-5:z-1])
-
-def projection_from_surface(ma, surface):
-    flush_message("Generating projection from surface...")
-    xmax, ymax, zmax = ma.shape
-    res = np.zeros([xmax, ymax], dtype=np.uint8)
-    for x in range(0, xmax):
-        for y in range(0, ymax):
-            z = surface[x, y]
-            res[x, y] = vavg(ma, x, y, z)
-            #res[x, y] = ma[x, y, z]
-    print " done"
-    return res
-
-def save_numpy_as_png(filename, np):
-    scipy.misc.imsave(filename, np)
 
 def main():
-
     try:
         stackdir = sys.argv[1]
     except IndexError:
         print "Usage: %s stack_dir [output_dir] [sdx] [sdy] [sdz]" % os.path.basename(sys.argv[0])
         sys.exit(1)
 
-    #imgpattern, istart, iend = stackhandle.get_stack_pattern(stackdir)
     ifiles = stackhandle.get_stack_files(stackdir)
 
     try:
@@ -109,53 +59,48 @@ def main():
 
     print "Called with", stackdir, output_dir
 
-    #impattern = 'data/pngstack/ExpID3002_spch4_TL003_plantD_lif_S000_T000_C000_Z0%02d.png'
-    #istart = 0
-    #iend = 92
-
     try:
         sdx, sdy, sdz = int(sys.argv[3]), int(sys.argv[4]), int(sys.argv[5])
     except IndexError:
         print "Using default values for standard deviation"
         sdx, sdy, sdz = 8, 8, 6
 
-    #sds = 10
     sds = 5
 
     print "Using standard deviations: %d, %d, %d" % (sdx, sdy, sdz)
 
-    #ma = load_png_stack(imgpattern, istart, iend)
     ma = load_png_stack(ifiles)
 
+    flush_message("Applying 3D gaussian filter...")
+    bl = nd.gaussian_filter(ma, [sdx, sdy, sdz])
+    print "done"
 
-
-    bl = apply_gaussian_filter(ma, [sdx, sdy, sdz])
-
-    #ps = find_projection_surface(bl)
-    ps = np.argmax(bl, 2)
+    flush_message("Finding projection surface...")
+    ps = proj.max_indices_z(bl)
+    print "done"
 
     sps = nd.gaussian_filter(ps, sds)
     _, _, zmax = ma.shape
     vis_factor = 255 / zmax
     sfilename = os.path.join(output_dir, "surface-g3d-%d-%d-%d-%d.png" % (sdx, sdy, sdz, sds))
-    save_numpy_as_png(sfilename, sps * vis_factor)
+    scipy.misc.imsave(sfilename, sps * vis_factor)
 
-    res = projection_from_surface(ma, sps)
+    flush_message("Generating projection from surface...")
+    res = proj.projection_from_surface_z(ma, sps, dm=3, dp=0) 
+    print "done"
 
     filename = os.path.join(output_dir, "proj-g3d-%d-%d-%d-%d.png" % (sdx, sdy, sdz, sds))
     pmax = np.amax(res)
-
     vis_scale = 255 / pmax
-
     scipy.misc.imsave(filename, res * vis_scale)
 
     flush_message("Post processing...")
     pp = projpp.proj_filter(res * vis_scale, 3, 60, 15)
-    print " done"
+    print "done"
+
     filename = os.path.join(output_dir, 'proj-pp-%d-%d-%d-%d.png' % (sdx, sdy, sdz, sds))
     scipy.misc.imsave(filename, pp)
 
-    #numpy_draw_pil(res)
 
 if __name__ == "__main__":
     main()
diff --git a/src/gaussproj_fromlif_all.py b/src/gaussproj_fromlif_all.py
@@ -14,7 +14,8 @@
 from pylab import imshow, show
 from xml import etree as et
 import matplotlib.pyplot as plt
-from gaussproj import *
+from gaussproj import flush_message
+import proj
 
 
 def parse_xml_metadata(xml_string, array_order='zyx'):
@@ -72,10 +73,10 @@ def main(args):
 		sys.exit(1)
 
 	try:
-        sdx, sdy, sdz = int(sys.argv[2]), int(sys.argv[3]), int(sys.argv[4])
-    except IndexError:
-        print "Using default values for standard deviation"
-        sdx, sdy, sdz = 4, 4, 3
+            sdx, sdy, sdz = int(sys.argv[2]), int(sys.argv[3]), int(sys.argv[4])
+        except IndexError:
+            print "Using default values for standard deviation"
+            sdx, sdy, sdz = 4, 4, 3
 
 	md = bf.get_omexml_metadata(lifdir)
 	mdo = bf.OMEXML(md)	
@@ -119,20 +120,20 @@ def main(args):
 		max_proj = np.amax(ma, axis=2)
 
 		mpfilename = os.path.join(output_dir, "max-proj.png")
-		save_numpy_as_png(mpfilename, max_proj)
+		scipy.misc.imsave(mpfilename, max_proj)
 
-		bl = apply_gaussian_filter(ma, [sdx, sdy, sdz])
+		bl = nd.gaussian_filter(ma, [sdx, sdy, sdz])
 
 		#ps = find_projection_surface(bl)
-		ps = np.argmax(bl, 2)
+		ps = proj.max_indices_z(bl)
 
 		sps = nd.gaussian_filter(ps, sds)
 
 		vis_factor = 255 / z_size
 		sfilename = os.path.join(output_dir, "surface-g3d-%d-%d-%d-%d.png" % (sdx, sdy, sdz, sds))
-		save_numpy_as_png(sfilename, sps * vis_factor)
+		scipy.misc.imsave(sfilename, sps * vis_factor)
 
-		res = projection_from_surface(ma, sps)
+		res = proj.projection_from_surface(ma, sps, dm=3, dp=0)
 
 		filename = os.path.join(output_dir, "proj-g3d-%d-%d-%d-%d.png" % (sdx, sdy, sdz, sds))
 		pmax = np.amax(res)

diff --git a/src/proj.py b/src/proj.py
@@ -0,0 +1,178 @@
+#!/usr/bin/env python2.7
+
+import scipy
+import scipy.ndimage as nd
+import numpy as np
+from math import sqrt
+import numpy.linalg as la
+
+# Methods for simple projection in the z-direction
+
+def max_indices_z(A, clip_bottom=True):
+    """
+    From a 3D matrix representing an image stack, 
+    find the depth of the voxel with maximum intensity in the z-direction.
+
+    Args:
+        A (numpy.ndarray): 3D volumetric image data 
+
+        clip_bottom (Optional[bool]): For columns which have uniform signal,
+                                      whether to return the index of the top or
+                                      bottom (default) of the image stack
+
+    Returns:
+        numpy.ndarray[int]: 2D array containing the depth of the projection surface
+    """
+
+    h = np.argmax(A, 2)
+    if clip_bottom:
+        h[np.logical_and(h==0, A[:,:,0]==A[:,:,-1])] = A.shape[2]-1
+    return h
+
+def threshold_indices_z(A, b, c, clip_bottom=True):
+    """
+    From a 3D matrix representing an image stack, 
+    find the depth of the first voxel which exceeds a threshold intensity
+    (depending on the mean intensity of the pixels in this "pillar")
+
+    Args:
+        A (numpy.ndarray): 3D volumetric image data 
+        b, c (float): Threshold is taken to be b*mean of the pillar + c
+
+        clip_bottom (Optional[bool]): For columns which have uniform signal,
+                                      whether to return the index of the top or
+                                      bottom (default) of the image stack
+
+    Returns:
+        numpy.ndarray[int]: 2D array containing the depth of the projection surface
+    """
+
+    t = b*np.mean(A, axis=2) + c
+    # now wish to find first occurance
+    return max_indices_z(A > t[:,:, np.newaxis], clip_bottom)
+
+
+def sheet_indices_z(A, niter=1000, dt=0.05, a=3.0, b=3.0):
+    """
+    From a 3D matrix representing an image stack, find the projection
+    surface using an active contour. The projection surface starts at the 
+    smallest z-depth (h=0), and moves in the direction of increasing z;
+    this downwards speed is reduced in the presence of signal. A smoothing
+    term is also applied to the surface, which encourages nearby points
+    to be at similar depths.
+
+    Args:
+        A (numpy.ndarray): 3D volumetric image data 
+        niter (int): Number of timesteps for the evolution of the 
+                     projection surface
+        dt (float): Timestep for the surface evolution. 
+                    May need to be limited to maintain stability 
+                    of the numerical method.
+        a (float): Smoothing parameter - larger values penalize differences
+                   in the depth of adjacent points.
+        b (float): Parameter controlling how much the image intensity slows the
+                   downwards movement of the projection surface; larger values
+                   result in the image intensity slowing the projection 
+                   surface more.
+                   
+    Returns:
+        numpy.ndarray[float]: 2D array containing the depth of the projection surface
+    """
+
+    h0 = np.zeros(A.shape[0:2])
+    i, j = np.meshgrid(np.arange(A.shape[0]), np.arange(A.shape[1]), order='ij')
+    d2 = np.zeros(h0.shape)
+    for k in range(niter):
+        d2[1:-1, 1:-1] = h0[:-2, 1:-1]+ h0[2:,1:-1] + h0[1:-1,:-2] + h0[1:-1,2:] - 4*h0[1:-1, 1:-1]
+        A_h = A[i.T, j.T, h0.astype(int)]
+        h0 += dt/(1+b*A_h) + d2*a*dt
+        h0 = np.clip(h0, 0, A.shape[2]-1)
+    return h0    
+
+
+def projection_from_surface_z(ma, surface, dm=0, dp=10, op=np.max):
+    """
+    Project the 3D image stack onto a surface (orthogonally, in the
+    z-direction. This involves applying an operator (maximum, or mean)
+    to the pixels in a band (in the z-direction) about the projection surface.
+
+    Args:
+        ma: 3D image stack.
+        surface: 2D indices of the projection surface.
+        dm: Number of pixels considered "above" (negative z) the surface
+        dp: Number of pixels (-1) considered "below" the surface
+        op: Operator used to calculate the intensity of the projected 
+            image im[i,j] from the band of pixels 
+            ma[i, j, surface-dm:surface+dp]
+
+    Returns:
+        numpy.ndarray: Projected image
+    """
+    xmax, ymax, zmax = ma.shape
+    res = np.zeros([xmax, ymax], dtype=ma.dtype)
+    z0 = np.clip(surface-dm, 0, zmax-1).astype(int)
+    z1 = np.clip(surface+dp, 1, zmax).astype(int)
+    for x in range(0, xmax):
+        for y in range(0, ymax):
+            res[x, y] = op(ma[x, y, z0[x,y]:z1[x,y]])
+    return res
+
+def projection_from_surface_normal_z(ma, surface, dm=0, dp=10, 
+                                        op=np.max, gradient_radius=0.5):
+    """
+    Project the 3D image stack, examining a ray of pixels *normal* to 
+    the projection surface. For each point in 2D, we consider a single ray 
+    passing through this point, normal to the surface. Points are sampled
+    along this ray (from a slightly blurred version of the 3D stack), and
+    an operator (default maximum) applied to these samples to give the 
+    intensity of the 2D image.
+
+    Args:
+        ma: 3D image stack.
+        surface: depth of the projection surface (2D array)
+        dm: Number of pixels considered "above" (negative z) the surface
+        dp: Number of pixels (-1) considered "below" the surface
+        op: Operator used to calculate the intensity of the projected 
+            image im[i,j] from the band of pixels 
+            ma[i, j, surface-dm:surface+dp]
+        gradient_radius: Radius of the gaussian operator used to find the slope
+                         of the projection surface.
+
+    Returns:
+        numpy.ndarray: Projected image
+    """
+
+
+    imax, jmax, kmax = ma.shape
+    res = np.zeros([imax, jmax], dtype=ma.dtype)
+    surface_di = scipy.ndimage.gaussian_filter1d(surface, 0.5, axis=0, order=1)
+    surface_dj = scipy.ndimage.gaussian_filter1d(surface, 0.5, axis=1, order=1)
+
+    def clip(i_,j_,k_):
+        return np.array( (max(0, min(i_, imax-1)), max(0, min(j_, jmax-1)), max(0, min(k_, kmax-1))) )
+
+    for i in range(0, imax):
+        for j in range(0, jmax):
+            si = -surface_di[i,j]
+            sj = -surface_dj[i,j]
+            k = surface[i,j]
+            sk = 1.0
+            h = sqrt(si*si+sj*sj+sk*sk)
+            si = si/h
+            sj = sj/h
+            sk = sk/h
+            p_start = clip(i-dm*si, j-dm*sj, k-dm*sk)
+            p_end = clip(i+dp*si, j+dp*sj, k+dp*sk)
+            ns = int(la.norm(p_end-p_start) +1)
+            tt = np.linspace(0, 1, ns)
+            res[i,j] = op( scipy.ndimage.map_coordinates(ma, 
+                            ([(p_start[0]*(1-t) + p_end[0]*t) for t in tt],
+                             [(p_start[1]*(1-t) + p_end[1]*t) for t in tt],
+                             [(p_start[2]*(1-t) + p_end[2]*t) for t in tt]),
+                                                         order=1 ))
+    return res
+
+
+
+
+