Merge pull request #54 from jni/memreduce

Reduce memory usage and speed up computation

gala/classify.py

@@ -309,7 +309,7 @@ def boundary_overlap_threshold(boundary_idxs, gt, tol_false, tol_true):
 def make_thresholded_boundary_overlap_loss(tol_false, tol_true):
     """Return a merge loss function based on boundary overlaps."""
     def loss(g, n1, n2, gt):
-        boundary_idxs = list(g[n1][n2]['boundary'])
+        boundary_idxs = g[n1][n2]['boundary']
         return \
             boundary_overlap_threshold(boundary_idxs, gt, tol_false, tol_true)
     return loss
@@ -326,7 +326,7 @@ def label_merges(g, merge_history, feature_map_function, gt, loss_function):
         n1, n2 = nodes
         features[i,:] = feature_map_function(g, n1, n2)
         labels[i] = loss_function(g, n1, n2, gt)
-        labeled_image.ravel()[list(g[n1][n2]['boundary'])] = 2+labels[i]
+        labeled_image.ravel()[g[n1][n2]['boundary']] = 2+labels[i]
         g.merge_nodes(n1,n2)
     return features, labels, labeled_image
 

gala/features/contact.pyx

@@ -60,7 +60,7 @@ class Manager(base.Null):
         return feature_vector
 
     def create_edge_cache(self, g, n1, n2):
-        edge_idxs = np.array(list(g[n1][n2]['boundary']))
+        edge_idxs = np.array(g[n1][n2]['boundary'])
         n1_idxs = np.array(list(g.extent(n1)))
         n2_idxs = np.array(list(g.extent(n2)))
         if self.oriented: ar = g.oriented_probabilities_r

gala/features/convex_hull.py

@@ -31,7 +31,7 @@ def write_fm(self, json_fm={}):
     def convex_hull_ind(self, g, n1, n2=None):
         m = np.zeros_like(g.watershed); 
         if n2 is not None:
-            m.ravel()[list(g[n1][n2]['boundary'])]=1
+            m.ravel()[g[n1][n2]['boundary']]=1
         else:
             m.ravel()[list(g.extent(n1))] = 1
         m = m - nd.binary_erosion(m) #Only need border

gala/features/histogram.pyx

@@ -110,7 +110,7 @@ class Manager(base.Null):
         return self.histogram(ar[node_idxs,:])
 
     def create_edge_cache(self, g, n1, n2):
-        edge_idxs = list(g[n1][n2]['boundary'])
+        edge_idxs = g[n1][n2]['boundary']
         if self.oriented:
             ar = g.oriented_probabilities_r
         else:

gala/features/inclusion.py

@@ -32,11 +32,11 @@ def compute_edge_features(self, g, n1, n2, cache=None):
                               for x in g.neighbors(n1)])
         bd_lengths2 = sorted([len(g[n2][x]['boundary'])
                               for x in g.neighbors(n2)])
-        ratios1 = [float(len(g[n1][n2]["boundary"]))/float(sum(bd_lengths1)),
-                   float(len(g[n1][n2]["boundary"]))/float(sum(bd_lengths2))]
+        ratios1 = [float(len(g[n1][n2]['boundary']))/float(sum(bd_lengths1)),
+                   float(len(g[n1][n2]['boundary']))/float(sum(bd_lengths2))]
         ratios1.sort()
-        ratios2 = [float(len(g[n1][n2]["boundary"]))/float(max(bd_lengths1)),
-                   float(len(g[n1][n2]["boundary"]))/float(max(bd_lengths2))]
+        ratios2 = [float(len(g[n1][n2]['boundary']))/float(max(bd_lengths1)),
+                   float(len(g[n1][n2]['boundary']))/float(max(bd_lengths2))]
         ratios2.sort()
         return np.concatenate((ratios1, ratios2))
 

gala/features/moments.pyx

@@ -44,7 +44,7 @@ class Manager(base.Null):
         return self.compute_moment_sums(ar, node_idxs)
 
     def create_edge_cache(self, g, n1, n2):
-        edge_idxs = list(g[n1][n2]['boundary'])
+        edge_idxs = g[n1][n2]['boundary']
         if self.oriented:
             ar = g.oriented_probabilities_r
         else:

gala/features/orientation.py

@@ -36,7 +36,7 @@ def create_node_cache(self, g, n):
     def create_edge_cache(self, g, n1, n2):
         # Get subscripts of extent (morpho.unravel_index was slow)
         M = np.zeros_like(g.watershed); 
-        M.ravel()[list(g[n1][n2]['boundary'])] = 1 
+        M.ravel()[g[n1][n2]['boundary']] = 1 
         ind = np.array(np.nonzero(M)).T
         # Get second moment matrix
         smm = np.cov(ind.T)/float(len(ind))

gala/features/squiggliness.py

@@ -41,5 +41,5 @@ def compute_edge_features(self, g, n1, n2, cache=None):
             cache = g[n1][n2][self.default_cache]
         m, M = cache[:self.ndim], cache[self.ndim:]
         plane_surface = np.sort(M-m)[1:].prod() * (3.0-g.pad_thickness)
-        return np.array([len(g[n1][n2]['boundary']) / plane_surface])
+        return np.array([len(set(g[n1][n2]['boundary'])) / plane_surface])
 

gala/morpho.py

@@ -543,20 +543,80 @@ def build_neighbors_array(ar, connectivity=1):
     idxs = arange(ar.size, dtype=uint32)
     return get_neighbor_idxs(ar, idxs, connectivity)
 
+
+def raveled_steps_to_neighbors(shape, connectivity=1):
+    """Compute the stepsize along all axes for given connectivity and shape.
+
+    Parameters
+    ----------
+    shape : tuple of int
+        The shape of the array along which we are stepping.
+    connectivity : int in {1, 2, ..., ``len(shape)``}
+        The number of orthogonal steps we can take to reach a "neighbor".
+
+    Returns
+    -------
+    steps : array of int64
+        The steps needed to get to neighbors from a particular raveled
+        index.
+
+    Examples
+    --------
+    >>> shape = (5, 4, 9)
+    >>> steps = raveled_steps_to_neighbors(shape)
+    >>> sorted(steps)
+    [-36, -9, -1, 1, 9, 36]
+    >>> steps2 = raveled_steps_to_neighbors(shape, 2)
+    >>> sorted(steps2)
+    [-45, -37, -36, -35, -27, -10, -9, -8, -1, 1, 8, 9, 10, 27, 35, 36, 37, 45]
+    """
+    stepsizes = np.cumprod((1,) + shape[-1:0:-1])[::-1]
+    steps = []
+    steps.extend((stepsizes, -stepsizes))
+    for nhops in range(2, connectivity + 1):
+        prod = np.array(list(it.product(*([[1, -1]] * nhops))))
+        multisteps = np.array(list(it.combinations(stepsizes, nhops))).T
+        steps.append(np.dot(prod, multisteps).ravel())
+    return np.concatenate(steps).astype(np.int64)
+
+
 def get_neighbor_idxs(ar, idxs, connectivity=1):
-    if isscalar(idxs): # in case only a single idx is given
+    """Return indices of neighboring voxels given array, indices, connectivity.
+
+    Parameters
+    ----------
+    ar : ndarray
+        The array in which neighbors are to be found.
+    idxs : int or container of int
+        The indices for which to find neighbors.
+    connectivity : int in {1, 2, ..., ``ar.ndim``}
+        The number of orthogonal steps allowed to be considered a
+        neighbor.
+
+    Returns
+    -------
+    neighbor_idxs : 2D array, shape (nidxs, nneighbors)
+        The neighbor indices for each index passed.
+
+    Examples
+    --------
+    >>> ar = np.arange(16).reshape((4, 4))
+    >>> ar
+    array([[ 0,  1,  2,  3],
+           [ 4,  5,  6,  7],
+           [ 8,  9, 10, 11],
+           [12, 13, 14, 15]])
+    >>> get_neighbor_idxs(ar, [5, 10], connectivity=1)
+    array([[ 9,  6,  1,  4],
+           [14, 11,  6,  9]])
+    >>> get_neighbor_idxs(ar, 9, connectivity=2)
+    array([[13, 10,  5,  8, 14, 12,  6,  4]])
+    """
+    if isscalar(idxs):  # in case only a single idx is given
         idxs = [idxs]
-    idxs = array(idxs) # in case a list or other array-like is given
-    strides = array(ar.strides)/ar.itemsize
-    if connectivity == 1: 
-        steps = (strides, -strides)
-    else:
-        steps = []
-        for i in range(1,connectivity+1):
-            prod = array(list(it.product(*([[1,-1]]*i))))
-            i_strides = array(list(it.combinations(strides,i))).T
-            steps.append(prod.dot(i_strides).ravel())
-    return idxs[:,newaxis] + concatenate(steps).astype(int32)
+    idxs = array(idxs)  # in case a list or other array-like is given
+    steps = raveled_steps_to_neighbors(ar.shape, connectivity)
+    return idxs[:, np.newaxis] + steps
 
 def orphans(a):
     """Find all the segments that do not touch the volume boundary.