Merge branch 'master' of github.com:jni/skan

.bumpversion.cfg

@@ -0,0 +1,20 @@
+[bumpversion]
+commit = True
+tag = True
+current_version = 0.2.0-dev
+parse = (?P<major>\d+)\.(?P<minor>\d+)(\.(?P<patch>\d+))?(\-(?P<release>\w+))?
+serialize = 
+	{major}.{minor}.{patch}-{release}
+	{major}.{minor}.{patch}
+	{major}.{minor}
+
+[bumpversion:file:setup.py]
+
+[bumpversion:file:skan/__init__.py]
+
+[bumpversion:part:release]
+optional_value = final
+values = 
+	dev
+	final
+

.travis.yml

@@ -32,7 +32,7 @@ script:
   - py.test
 
 after_success:
-  - if [[ $TRAVIS_PYTHON_VERSION == 3.5 ]]; then make numba-clean; py.test --cov-report term-missing --cov .; coveralls; fi
+  - if [[ $TRAVIS_PYTHON_VERSION == 3.5 ]]; then make numba-clean; export NUMBA_DISABLE_JIT=1; py.test --cov-report term-missing --cov .; coveralls; fi
 
 cache:
     directories:

MANIFEST.in

@@ -0,0 +1,2 @@
+include LICENSE
+recursive-include skan/test *.py

benchmarks/bench_skan.py

@@ -27,6 +27,10 @@ def bench_suite():
         g, indices, degrees = csr.skeleton_to_csgraph(skeleton,
                                                       spacing=2.24826)
     times['build graph'] = t_build_graph[0]
+    with timer() as t_build_graph2:
+        g, indices, degrees = csr.skeleton_to_csgraph(skeleton,
+                                                      spacing=2.24826)
+    times['build graph again'] = t_build_graph2[0]
     with timer() as t_stats:
         stats = csr.branch_statistics(g, indices, degrees)
     times['compute statistics'] = t_stats[0]

requirements.txt

@@ -1,5 +1,5 @@
 numpy>=1.11
-scipy>=1.18
+scipy>=0.18
 networkx>=1.11
 numba>=0.29
 pandas>=0.18

setup.py

@@ -11,10 +11,10 @@
 LONG_DESCRIPTION    = descr
 MAINTAINER          = 'Juan Nunez-Iglesias'
 MAINTAINER_EMAIL    = 'juan.n@unimelb.edu.au'
-URL                 = 'https://github.com/jni/skeleton-analysis'
+URL                 = 'https://github.com/jni/skan'
 LICENSE             = 'BSD 3-clause'
-DOWNLOAD_URL        = 'https://github.com/jni/skeleton-analysis'
-VERSION             = '0.1-dev'
+DOWNLOAD_URL        = 'https://github.com/jni/skan'
+VERSION             = '0.2.0-dev'
 PYTHON_VERSION      = (3, 5)
 INST_DEPENDENCIES   = {} 
 

skan/__init__.py

@@ -1,5 +1,7 @@
 from .csr import skeleton_to_csgraph, branch_statistics, summarise
 
+__version__ = '0.2.0-dev'
+
 __all__ = ['skeleton_to_csgraph',
            'branch_statistics',
-           'summarise']
+           'summarise']

skan/_testdata.py

@@ -45,3 +45,5 @@
                         [0, 1, 0, 0, 1],
                         [1, 0, 1, 0, 1],
                         [0, 0, 0, 0, 1]]], dtype=bool)
+
+topograph1d = np.array([3., 2., 3.])

skan/csr.py

@@ -36,6 +36,19 @@ def neighbors(self, row):
         return self.indices[loc:stop]
 
 
+def _pixel_graph(image, steps, distances, num_edges, height=None):
+    row = np.empty(num_edges, dtype=int)
+    col = np.empty(num_edges, dtype=int)
+    data = np.empty(num_edges, dtype=float)
+    if height is None:
+        _write_pixel_graph(image, steps, distances, row, col, data)
+    else:
+        _write_pixel_graph_height(image, height, steps, distances,
+                                  row, col, data)
+    graph = sparse.coo_matrix((data, (row, col))).tocsr()
+    return graph
+
+
 @numba.jit(nopython=True, cache=True, nogil=True)
 def _write_pixel_graph(image, steps, distances, row, col, data):
     """Step over `image` to build a graph of nonzero pixel neighbors.
@@ -81,6 +94,60 @@ def _write_pixel_graph(image, steps, distances, row, col, data):
                     data[k] = distances[j]
                     k += 1
 
+@numba.jit(nopython=True, cache=True, nogil=True)
+def _write_pixel_graph_height(image, height, steps, distances, row, col, data):
+    """Step over `image` to build a graph of nonzero pixel neighbors.
+
+    Parameters
+    ----------
+    image : int array
+        The input image.
+    height : float array, same shape as `image`
+        This is taken to be a height map along an additional
+        dimension (in addition to the image dimensions), so the distance
+        between two neighbors `i` and `n` separated by `j` is given by:
+
+             `np.sqrt(distances[j]**2 + (height[i] - height[n])**2)`
+
+    steps : int array, shape (N,)
+        The raveled index steps to find a pixel's neighbors in `image`.
+    distances : float array, shape (N,)
+        The euclidean distance from a pixel to its corresponding
+        neighbor in `steps`.
+    row : int array
+        Output array to be filled with the "center" pixel IDs.
+    col : int array
+        Output array to be filled with the "neighbor" pixel IDs.
+    data : float array
+        Output array to be filled with the distances from center to
+        neighbor pixels.
+
+    Notes
+    -----
+    No size or bounds checking is performed. Users should ensure that
+    - No index in `indices` falls on any edge of `image` (or the
+      neighbor computation will fail or segfault).
+    - The `steps` and `distances` arrays have the same shape.
+    - The `row`, `col`, `data` are long enough to hold all of the
+      edges.
+    """
+    image = image.ravel()
+    height = height.ravel()
+    n_neighbors = steps.size
+    start_idx = np.max(steps)
+    end_idx = image.size + np.min(steps)
+    k = 0
+    for i in range(start_idx, end_idx + 1):
+        if image[i] != 0:
+            for j in range(n_neighbors):
+                n = steps[j] + i
+                if image[n] != 0:
+                    row[k] = image[i]
+                    col[k] = image[n]
+                    data[k] = np.sqrt(distances[j] ** 2 +
+                                      (height[i] - height[n]) ** 2)
+                    k += 1
+
 
 def skeleton_to_csgraph(skel, *, spacing=1):
     """Convert a skeleton image of thin lines to a graph of neighbor pixels.
@@ -112,6 +179,10 @@ def skeleton_to_csgraph(skel, *, spacing=1):
         An image where each pixel value contains the degree of its
         corresponding node in `graph`. This is useful to classify nodes.
     """
+    if np.issubdtype(skel.dtype, float):  # interpret float skels as height
+        height = pad(skel, 0.)
+    else:
+        height = None
     skel = skel.astype(bool)  # ensure we have a bool image
                               # since we later use it for bool indexing
     spacing = np.ones(skel.ndim, dtype=float) * spacing
@@ -127,12 +198,9 @@ def skeleton_to_csgraph(skel, *, spacing=1):
     degree_image = ndi.convolve(skel.astype(int), degree_kernel,
                                 mode='constant') * skel
     num_edges = np.sum(degree_image)  # *2, which is how many we need to store
-    row, col = np.zeros(num_edges, dtype=int), np.zeros(num_edges, dtype=int)
-    data = np.zeros(num_edges, dtype=float)
     steps, distances = raveled_steps_to_neighbors(skelint.shape, ndim,
                                                   spacing=spacing)
-    _write_pixel_graph(skelint, steps, distances, row, col, data)
-    graph = sparse.coo_matrix((data, (row, col))).tocsr()
+    graph = _pixel_graph(skelint, steps, distances, num_edges, height)
     return graph, pixel_indices, degree_image
 
 
@@ -283,9 +351,9 @@ def submatrix(M, idxs):
     --------
     >>> Md = np.arange(16).reshape((4, 4))
     >>> M = sparse.csr_matrix(Md)
-    >>> submatrix(M, [0, 2]).toarray()
-    array([[ 0,  2],
-           [ 8, 10]], dtype=int64)
+    >>> print(submatrix(M, [0, 2]).toarray())
+    [[ 0  2]
+     [ 8 10]]
     """
     Msub = M[idxs, :][:, idxs]
     return Msub

skan/test/test_csr.py

@@ -7,7 +7,7 @@
 sys.path.append(rundir)
 
 from skan._testdata import (tinycycle, tinyline, skeleton1, skeleton2,
-                            skeleton3d)
+                            skeleton3d, topograph1d)
 
 
 def test_tiny_cycle():
@@ -80,3 +80,10 @@ def test_3d_spacing():
     assert_equal(stats.shape, (5, 4))
     assert_almost_equal(stats[0], [1, 11, 2 * np.sqrt(27), 1])
     assert_equal(np.unique(stats[:, 3].astype(int)), [1, 2, 3])
+
+
+def test_topograph():
+    g, idxs, degimg = csr.skeleton_to_csgraph(topograph1d)
+    stats = csr.branch_statistics(g, idxs, degimg)
+    assert stats.shape == (1, 4)
+    assert_almost_equal(stats[0], [1, 3, 2 * np.sqrt(2), 0])

test-environment.yml

@@ -5,8 +5,8 @@ dependencies:
 - numba>=0.29
 - numpy>=1.11*
 - numpydoc>=0.5*
-- scipy>=0.17*
-- networkx>=1.10*
+- scipy>=0.18*
+- networkx>=1.11*
 - pytest>=3*
 - coverage>=4.0
 - pytest-cov>=2.2