From 0f31034a460d73732ec8f28ea200e9c1424aa799 Mon Sep 17 00:00:00 2001
From: emmanuelle <emmanuelle.gouillart@nsup.org>
Date: Sat, 22 Oct 2011 22:38:33 +0200
Subject: [PATCH 1/6] First integration of cell profiler medial axis transform

(skeletonization in cell profiler)
---
 scikits/image/morphology/_cpmorphology2.pyx  | 309 +++++++++++++++++++
 skimage/morphology/__init__.py               |   2 +-
 skimage/morphology/skeletonize.py            | 265 +++++++++++++++-
 skimage/morphology/tests/test_skeletonize.py |  51 ++-
 skimage/transform/setup.py                   |   1 +
 5 files changed, 623 insertions(+), 5 deletions(-)
 create mode 100644 scikits/image/morphology/_cpmorphology2.pyx

diff --git a/scikits/image/morphology/_cpmorphology2.pyx b/scikits/image/morphology/_cpmorphology2.pyx
new file mode 100644
index 00000000..d7b76413
--- /dev/null
+++ b/scikits/image/morphology/_cpmorphology2.pyx
@@ -0,0 +1,309 @@
+'''_cpmorphology2.pyx - support routines for cpmorphology in Cython
+
+Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
+Website: http://www.cellprofiler.org
+
+Copyright (c) 2003-2009 Massachusetts Institute of Technology
+Copyright (c) 2009-2011 Broad Institute
+All rights reserved.
+
+Original author: Lee Kamentsky
+'''
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+cdef extern from "Python.h":
+    ctypedef int Py_intptr_t
+
+cdef extern from "numpy/arrayobject.h":
+    ctypedef class numpy.ndarray [object PyArrayObject]:
+        cdef char *data
+        cdef Py_intptr_t *dimensions
+        cdef Py_intptr_t *strides
+    cdef void import_array()
+    cdef int  PyArray_ITEMSIZE(np.ndarray)
+
+import_array()
+
+@cython.boundscheck(False)
+def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, 
+                                negative_indices=False, mode='c'] result,
+                     np.ndarray[dtype=np.int32_t, ndim=1,
+                                negative_indices=False, mode='c'] i,
+                     np.ndarray[dtype=np.int32_t, ndim=1,
+                                negative_indices=False, mode='c'] j,
+                     np.ndarray[dtype=np.int32_t, ndim=1,
+                                negative_indices=False, mode='c'] order,
+                     np.ndarray[dtype=np.uint8_t, ndim=1,
+                                negative_indices=False, mode='c'] table):
+    '''Inner loop of skeletonize function
+
+    result - on input, the image to be skeletonized, on output the skeletonized
+             image
+    i,j    - the coordinates of each foreground pixel in the image
+    order  - the index of each pixel, in the order of processing
+    table  - the 512-element lookup table of values after transformation
+
+    The loop determines whether each pixel in the image can be removed without
+    changing the Euler number of the image. The pixels are ordered by
+    increasing distance from the background which means a point nearer to
+    the quench-line of the brushfire will be evaluated later than a
+    point closer to the edge.
+    '''
+    cdef:
+        np.int32_t accumulator
+        np.int32_t index,order_index
+        np.int32_t ii,jj
+
+    for index in range(order.shape[0]):
+        accumulator = 16
+        order_index = order[index]
+        ii = i[order_index]
+        jj = j[order_index]
+        if ii > 0:
+            if jj > 0 and result[ii - 1, jj - 1]:
+                accumulator += 1
+            if result[ii - 1, jj]:
+                accumulator += 2
+            if jj < result.shape[1] - 1 and result[ii - 1, jj + 1]:
+                    accumulator += 4
+            if jj > 0 and result[ii, jj - 1]:
+                accumulator += 8
+            if jj < result.shape[1] - 1 and result[ii, jj + 1]:
+                accumulator += 32
+            if ii < result.shape[0]-1:
+                if jj > 0 and result[ii+1,jj-1]:
+                    accumulator += 64
+                if result[ii+1,jj]:
+                    accumulator += 128
+                if jj < result.shape[1]-1 and result[ii+1,jj+1]:
+                    accumulator += 256
+            result[ii,jj] = table[accumulator]
+
+@cython.boundscheck(False)
+def table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2,
+                                  negative_indices=False, mode='c'] image):
+    """
+    Return an index into a table per pixel of a binary image
+
+    Take the sum of true neighborhood pixel values where the neighborhood
+    looks like this:
+     1   2   4
+     8  16  32
+    64 128 256
+
+    This code could be replaced by a convolution with the kernel:
+    256 128 64
+     32  16  8
+      4   2  1
+    but this runs about twice as fast because of inlining and the
+    hardwired kernel.
+    """
+    cdef:
+        np.ndarray[dtype=np.int32_t, ndim=2, 
+                   negative_indices=False, mode='c'] indexer
+        np.int32_t *p_indexer
+        np.uint8_t *p_image
+        np.int32_t i_stride
+        np.int32_t i_shape
+        np.int32_t j_shape
+        np.int32_t i
+        np.int32_t j
+        np.int32_t offset
+
+    i_shape   = image.shape[0]
+    j_shape   = image.shape[1]
+    indexer = np.zeros((i_shape, j_shape), np.int32)
+    p_indexer = <np.int32_t *>indexer.data
+    p_image   = <np.uint8_t *>image.data
+    i_stride  = image.strides[0]
+    assert i_shape >= 3 and j_shape >= 3, \
+        "Please use the slow method for arrays < 3x3"
+
+    for i in range(1, i_shape-1):
+        offset = i_stride* i + 1
+        for j in range(1, j_shape - 1):
+            if p_image[offset]:
+                p_indexer[offset + i_stride + 1] += 1
+                p_indexer[offset + i_stride] += 2
+                p_indexer[offset + i_stride - 1] += 4
+                p_indexer[offset + 1] += 8
+                p_indexer[offset] += 16
+                p_indexer[offset - 1] += 32
+                p_indexer[offset - i_stride + 1] += 64
+                p_indexer[offset - i_stride] += 128
+                p_indexer[offset - i_stride - 1] += 256
+            offset += 1
+    #
+    # Do the corner cases (literally)
+    #
+    if image[0, 0]:
+        indexer[0, 0] += 16
+        indexer[0, 1] += 8
+        indexer[1, 0] += 2
+        indexer[1, 1] += 1
+
+    if image[0, j_shape - 1]:
+        indexer[0, j_shape - 2] += 32
+        indexer[0, j_shape - 1] += 16
+        indexer[1, j_shape - 2] += 4
+        indexer[1, j_shape - 1] += 2
+
+    if image[i_shape - 1, 0]:
+        indexer[i_shape - 2, 0] += 128
+        indexer[i_shape - 2, 1] += 64
+        indexer[i_shape - 1, 0] += 16
+        indexer[i_shape - 1, 1] += 8
+
+    if image[i_shape - 1, j_shape - 1]:
+        indexer[i_shape - 2, j_shape - 2] += 256
+        indexer[i_shape - 2, j_shape - 1] += 128
+        indexer[i_shape - 1, j_shape - 2] += 32
+        indexer[i_shape - 1, j_shape - 1] += 16
+    #
+    # Do the edges
+    #
+    for j in range(1, j_shape - 1):
+        if image[0, j]:
+            indexer[0, j - 1] += 32
+            indexer[0, j]   += 16
+            indexer[0, j + 1] += 8
+            indexer[1, j - 1] += 4
+            indexer[1, j]   += 2
+            indexer[1, j + 1] += 1
+        if image[i_shape - 1, j]:
+            indexer[i_shape - 2, j - 1] += 256
+            indexer[i_shape - 2, j]   += 128
+            indexer[i_shape - 2, j + 1] += 64
+            indexer[i_shape - 1, j - 1] += 32
+            indexer[i_shape - 1, j]   += 16
+            indexer[i_shape - 1, j + 1] += 8
+
+    for i in range(1, i_shape - 1):
+        if image[i, 0]:
+            indexer[i - 1, 0] += 128
+            indexer[i, 0]   += 16
+            indexer[i + 1, 0] += 2
+            indexer[i - 1, 1] += 64
+            indexer[i, 1]   += 8
+            indexer[i + 1, 1] += 1
+        if image[i, j_shape - 1]:
+            indexer[i - 1, j_shape - 2] += 256
+            indexer[i, j_shape - 2]   += 32
+            indexer[i + 1, j_shape - 2] += 4
+            indexer[i - 1, j_shape - 1] += 128
+            indexer[i, j_shape - 1]   += 16
+            indexer[i + 1, j_shape - 1] += 2
+    return indexer
+
+@cython.boundscheck(False)
+def index_lookup(np.ndarray[dtype=np.int32_t, ndim=1, 
+                                negative_indices=False] index_i, 
+                 np.ndarray[dtype=np.int32_t, ndim=1, 
+                                negative_indices=False] index_j, 
+                 np.ndarray[dtype=np.uint32_t, ndim=2, 
+                                negative_indices=False] image,
+                 table_in, 
+                 iterations=None):
+    """
+    Perform a table lookup for only the indexed pixels
+    
+    For morphological operations that only convert 1 to 0, the set of
+    resulting pixels is always a subset of the input set. Therefore, when
+    repeating, it will be faster to operate only on the subsets especially
+    when the results are 1-d or 0-d objects.
+    
+    This function returns a new index_i and index_j array of the pixels
+    that survive the operation. The image is modified in-place to remove
+    the pixels that did not survive.
+    
+    index_i - an array of row indexes into the image.
+    index_j - a similarly-shaped array of column indexes.
+    image - the binary image: *NOTE* add a row and column of border values
+            to the original image to account for pixels on the edge of the
+            image.
+    iterations - # of iterations to do, default is "forever"
+    
+    The idea of index_lookup was taken from
+    http://blogs.mathworks.com/steve/2008/06/13/performance-optimization-for-applylut/
+    which, apparently, is how Matlab achieved its bwmorph speedup.
+    """
+    cdef:
+        np.ndarray[dtype=np.uint8_t, ndim=1, 
+            negative_indices=False] table = table_in.astype(np.uint8)
+        np.uint32_t center, hit_count, idx, indexer
+        np.int32_t idxi, idxj
+
+    if iterations == None:
+        # Worst case - remove one per iteration
+        iterations = len(index_i)
+
+    for i in range(iterations):
+        hit_count = len(index_i)
+        with nogil:
+            #
+            # For integer images (i.e., labels), a neighbor point is
+            # "background" if it doesn't match the central value. This
+            # lets adjacent labeled objects shrink independently of each
+            # other.
+            #
+            for 0 <= idx < hit_count:
+                idxi, idxj = index_i[idx], index_j[idx]
+                center = image[idxi, idxj]
+                indexer =   ((image[idxi - 1, idxj - 1] == center) * 1 +
+                             (image[idxi - 1, idxj]     == center) * 2 +
+                             (image[idxi - 1, idxj + 1] == center) * 4 +
+                             (image[idxi,     idxj - 1] == center) * 8 +
+                                    16 +
+                             (image[idxi,     idxj + 1] == center) * 32 +
+                             (image[idxi + 1, idxj - 1] == center) * 64 +
+                             (image[idxi + 1, idxj]     == center) * 128 +
+                             (image[idxi + 1, idxj + 1] == center) * 256)
+                if table[indexer] == 0:
+                    # mark for deletion
+                    index_i[idx] = -index_i[idx]
+
+            # remove marked pixels
+            for 0 <= idx < hit_count:
+                idxi, idxj = index_i[idx], index_j[idx]
+                if idxi < 0:
+                    image[-idxi, idxj] = 0
+
+        index_j = index_j[index_i >= 0]
+        index_i = index_i[index_i >= 0]
+        if len(index_i) == hit_count:
+            break
+
+    return (index_i, index_j)
+
+def prepare_for_index_lookup(image, border_value):
+    """
+    Return the index arrays of "1" pixels and an image with an added border
+    
+    The routine, index_lookup takes an array of i indexes, an array of
+    j indexes and an image guaranteed to be indexed successfully by 
+    index_<i,j>[:] +/- 1. This routine constructs an image with added border
+    pixels... evilly, the index, 0 - 1, lands on the border because of Python's
+    negative indexing convention.
+    """
+    if np.issubdtype(image.dtype, float):
+        image = image.astype(bool)
+    image_i, image_j = np.argwhere(image.astype(bool)).transpose().\
+                                    astype(np.int32) + 1
+    output_image = (np.ones(np.array(image.shape) + 2, image.dtype) \
+                    if border_value
+                    else np.zeros(np.array(image.shape) + 2, image.dtype))
+    output_image[1:image.shape[0] + 1, 1:image.shape[1] + 1] = image
+    return (image_i, image_j, output_image.astype(np.uint32))
+
+
+def extract_from_image_lookup(orig_image, index_i, index_j):
+    """
+    Extract only one pixel
+    """
+    output = np.zeros(orig_image.shape, orig_image.dtype)
+    output[index_i - 1, index_j - 1] = orig_image[index_i - 1, index_j - 1]
+    return output
+
diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py
index d57cc859..efb92507 100644
--- a/skimage/morphology/__init__.py
+++ b/skimage/morphology/__init__.py
@@ -2,4 +2,4 @@ from grey import *
 from selem import *
 from .ccomp import label
 from watershed import watershed, is_local_maximum
-from skeletonize import skeletonize
+from skeletonize import skeletonize, medial_axis
diff --git a/skimage/morphology/skeletonize.py b/skimage/morphology/skeletonize.py
index ea478a9a..7d76f833 100644
--- a/skimage/morphology/skeletonize.py
+++ b/skimage/morphology/skeletonize.py
@@ -4,7 +4,13 @@ objects in an image.
 """
 
 import numpy as np
-from scipy.ndimage import correlate
+from scipy import ndimage
+
+from _cpmorphology2 import skeletonize_loop, table_lookup_index
+from _cpmorphology2 import extract_from_image_lookup, \
+            prepare_for_index_lookup, index_lookup
+
+# --------- Skeletonization by morphological thinning ---------
 
 def skeletonize(image):
     """Return the skeleton of a binary image.
@@ -24,6 +30,10 @@ def skeletonize(image):
     skeleton : ndarray
         A matrix containing the thinned image.
     
+    See also
+    --------
+    medial_axis
+
     Notes
     -----
     The algorithm [1] works by making successive passes of the image, 
@@ -107,7 +117,7 @@ def skeletonize(image):
         pixelRemoved = False;
 
         # assign each pixel a unique value based on its foreground neighbours
-        neighbours = correlate(skeleton, mask, mode='constant')
+        neighbours = ndimage.correlate(skeleton, mask, mode='constant')
         
         # ignore background
         neighbours *= skeleton
@@ -126,7 +136,7 @@ def skeletonize(image):
             skeleton[code_mask] = 0
 
         # pass 2 - remove the 2's and 3's
-        neighbours = correlate(skeleton, mask, mode='constant')
+        neighbours = ndimage.correlate(skeleton, mask, mode='constant')
         neighbours *= skeleton
         codes = np.take(lut, neighbours)
         code_mask = (codes == 2)
@@ -139,3 +149,252 @@ def skeletonize(image):
             skeleton[code_mask] = 0
 
     return skeleton
+
+# --------- Skeletonization by medial axis transform --------
+
+eight_connect = ndimage.generate_binary_structure(2, 2)
+
+
+def medial_axis(image, mask=None, return_distance=False):
+    """
+    Compute the medial axis transform of a binary image
+
+    Parameters
+    ----------
+
+    image: binary ndarray 
+    
+    mask: binary ndarray, optional
+        If a mask is given, only those elements with a true value in `mask`
+        are used for computing the medial axis.
+
+    return_distance; bool, optional
+        If true, the distance transform is returned as well as the skeleton.
+
+    Returns
+    -------
+
+    out: ndarray of bools
+        Medial axis transform of the image
+
+    dist: ndarray of ints
+        Distance transform of the image (only returned if `return_distance`
+        is True)
+   
+    See also
+    --------
+    skeletonize 
+
+    Notes
+    -----
+    This algorithm computes the medial axis transform of an image
+    as the ridges of its distance transform. First, the distance transform
+    is computed, then the foreground (value of 1) points are ordered by
+    the distance transform. In order to reduce the image to its skeleton, 
+    a point is removed if it has more than one neighbor and if removing it
+    does not change the Euler number (the connectivity).
+
+    Examples
+    --------
+    >>> square = np.zeros((7, 7), dtype=np.uint8)
+    >>> square[1:-1, 2:-2] = 1
+    >>> square
+    array([[0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 0, 0],
+           [0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
+    >>> morphology.medial_axis(square).astype(np.uint8)
+    array([[0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 1, 0, 1, 0, 0],
+           [0, 0, 0, 1, 0, 0, 0],
+           [0, 0, 0, 1, 0, 0, 0],
+           [0, 0, 0, 1, 0, 0, 0],
+           [0, 0, 1, 0, 1, 0, 0],
+           [0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
+
+    """
+    global eight_connect
+    if mask is None:
+        masked_image = image.astype(np.bool)
+    else:
+        masked_image = image.astype(bool).copy()
+        masked_image[~mask] = False
+    #
+    # Lookup table - start with only positive pixels.
+    # Keep if # pixels in neighborhood is 2 or less
+    # Keep if removing the pixel results in a different connectivity
+    # table is independent of image
+    table = (_make_table(True, 
+                        np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], bool),
+                        np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], bool)) &
+            (np.array([ndimage.label(_pattern_of(index), eight_connect)[1] !=
+                        ndimage.label(_pattern_of(index & ~ 2**4),
+                                    eight_connect)[1]
+                        for index in range(512)]) |
+        np.array([np.sum(_pattern_of(index)) < 3 for index in range(512)])))
+    distance = ndimage.distance_transform_edt(masked_image)
+    if return_distance:
+        store_distance = distance.copy()
+    #
+        # The processing order along the edge is critical to the shape of the
+    # resulting skeleton: if you process a corner first, that corner will
+    # be eroded and the skeleton will miss the arm from that corner. Pixels
+    # with fewer neighbors are more "cornery" and should be processed last.
+    #
+    cornerness_table = np.array([9 - np.sum(_pattern_of(index))
+                                 for index in range(512)])
+    corner_score = _table_lookup(masked_image, cornerness_table, False, 1)
+    i, j = np.mgrid[0:image.shape[0], 0:image.shape[1]]
+    result = masked_image.copy()
+    distance = distance[result]
+    i = np.ascontiguousarray(i[result], np.int32)
+    j = np.ascontiguousarray(j[result], np.int32)
+    result = np.ascontiguousarray(result, np.uint8)
+    #
+    # We use a random # for tiebreaking. Assign each pixel in the image a
+    # predictable, random # so that masking doesn't affect arbitrary choices
+    # of skeletons
+    #
+    # Why fix the seed? Should we pass a random number generator instead?
+    np.random.seed(0)
+    tiebreaker = np.random.permutation(np.arange(masked_image.sum()))
+    order = np.lexsort((tiebreaker,
+                        corner_score[masked_image],
+                        distance))
+    order = np.ascontiguousarray(order, np.int32)
+    table = np.ascontiguousarray(table, np.uint8)
+    # Remove pixels not belonging to the medial axis
+    skeletonize_loop(result, i, j, order, table)
+
+    result = result.astype(bool)
+    if not mask is None:
+        result[~mask] = image[~mask]
+    if return_distance:
+        return result, store_distance
+    else:
+        return result
+
+def _pattern_of(index):
+    """
+    Return the pattern represented by an index value
+    Byte decomposition of index
+    """
+    return np.array([[index & 2**0,index & 2**1,index & 2**2],
+                     [index & 2**3,index & 2**4,index & 2**5],
+                     [index & 2**6,index & 2**7,index & 2**8]], bool)
+
+
+def _table_lookup(image, table, border_value, iterations = None):
+    """
+    Perform a morphological transform on an image, directed by its 
+    neighbors
+    
+    Parameters
+    ----------
+    image - a binary image
+    table - a 512-element table giving the transform of each pixel given
+            the values of that pixel and its 8-connected neighbors.
+    border_value - the value of pixels beyond the border of the image.
+                   This should test as True or False.
+    
+    Returns
+    -------
+    result: ndarray of same shape as `image`
+        Transformed image
+
+    Notes
+    -----
+    The pixels are numbered like this:
+    
+    0 1 2
+    3 4 5
+    6 7 8
+    The index at a pixel is the sum of 2**<pixel-number> for pixels
+    that evaluate to true. 
+    """
+    #
+    # Test for a table that never transforms a zero into a one:
+    #
+    center_is_zero = np.array([(x & 2**4) == 0 for x in range(2**9)])
+    use_index_trick = False
+    if (not np.any(table[center_is_zero]) and
+        (np.issubdtype(image.dtype, bool) or np.issubdtype(image.dtype, int))):
+        # Use the index trick
+        use_index_trick = True
+        invert = False
+    elif (np.all(table[~center_is_zero]) and np.issubdtype(image.dtype, bool)):
+        # All ones stay ones, invert the table and the image and do the trick
+        use_index_trick = True
+        invert = True
+        image = ~ image
+        # table index 0 -> 511 and the output is reversed
+        table = ~ table[511-np.arange(512)]
+        border_value = not border_value
+    if use_index_trick:
+        orig_image = image
+        index_i, index_j, image = prepare_for_index_lookup(image, border_value)
+        index_i, index_j = index_lookup(index_i, index_j,
+                                        image, table, iterations)
+        image = extract_from_image_lookup(orig_image, index_i, index_j)
+        if invert:
+            image = ~ image
+        return image
+    print(use_index_trick)
+    counter = 0
+    while counter != iterations:
+        counter += 1
+        #
+        # We accumulate into the indexer to get the index into the table
+        # at each point in the image
+        #
+        if image.shape[0] < 3 or image.shape[1] < 3:
+            image = image.astype(bool)
+            indexer = np.zeros(image.shape,int)
+            indexer[1:, 1:]   += image[:-1, :-1] * 2**0
+            indexer[1:, :]    += image[:-1, :] * 2**1
+            indexer[1:, :-1]  += image[:-1, 1:] * 2**2
+            
+            indexer[:, 1:]    += image[:, :-1] * 2**3
+            indexer[:, :]     += image[:, :] * 2**4
+            indexer[:, :-1]   += image[:, 1:] * 2**5
+        
+            indexer[:-1, 1:]  += image[1:, :-1] * 2**6
+            indexer[:-1, :]   += image[1:, :] * 2**7
+            indexer[:-1, :-1] += image[1:, 1:] * 2**8
+        else:
+            indexer = table_lookup_index(np.ascontiguousarray(image, np.uint8))
+        if border_value:
+            indexer[0,:]   |= 2**0 + 2**1 + 2**2
+            indexer[-1,:]  |= 2**6 + 2**7 + 2**8
+            indexer[:,0]   |= 2**0 + 2**3 + 2**6
+            indexer[:,-1]  |= 2**2 + 2**5 + 2**8
+        new_image = table[indexer]
+        if np.all(new_image == image):
+            break
+        image = new_image
+    return image
+
+def _make_table(value, pattern, care=np.ones((3,3),bool)):
+    '''Return a table suitable for table_lookup
+    
+    value - set all table entries matching "pattern" to "value", all others
+            to not "value"
+    pattern - a 3x3 boolean array with the pattern to match
+    care    - a 3x3 boolean array where each value is true if the pattern
+              must match at that position and false if we don't care if
+              the pattern matches at that position.
+    '''
+    def fn(index, p, i, j):
+        '''Return true if bit position "p" in index matches pattern'''
+        return ((((index & 2**p) > 0) == pattern[i, j]) or not care[i, j])
+    return np.array([value     
+                     if (fn(i, 0, 0, 0) and fn(i, 1, 0, 1) and fn(i, 2, 0, 2) 
+                     and fn(i, 3, 1, 0) and fn(i, 4, 1, 1) and fn(i, 5, 1, 2) 
+                     and fn(i, 6, 2, 0) and fn(i, 7, 2, 1) and fn(i, 8, 2, 2))
+                     else not value
+                     for i in range(512)], bool)
+
+
diff --git a/skimage/morphology/tests/test_skeletonize.py b/skimage/morphology/tests/test_skeletonize.py
index 8d900cd7..c1b69119 100644
--- a/skimage/morphology/tests/test_skeletonize.py
+++ b/skimage/morphology/tests/test_skeletonize.py
@@ -1,5 +1,5 @@
 import numpy as np
-from skimage.morphology import skeletonize
+from skimage.morphology import skeletonize, medial_axis
 import numpy.testing
 from skimage.draw import draw
 from scipy.ndimage import correlate
@@ -90,6 +90,55 @@ class TestSkeletonize():
         blocks = correlate(result, mask, mode='constant')
         assert not numpy.any(blocks == 4)
         
+class TestMedialAxis():
+    def test_00_00_zeros(self):
+        '''Test skeletonize on an array of all zeros'''
+        result = medial_axis(np.zeros((10, 10), bool))
+        assert np.all(result == False)
+
+    def test_00_01_zeros_masked(self):
+        '''Test skeletonize on an array that is completely masked'''
+        result = medial_axis(np.zeros((10, 10), bool),
+                                   np.zeros((10, 10), bool))
+        assert np.all(result == False)
+    
+    def test_01_01_rectangle(self):
+        '''Test skeletonize on a rectangle'''
+        image = np.zeros((9, 15), bool)
+        image[1:-1, 1:-1] = True
+        #
+        # The result should be four diagonals from the
+        # corners, meeting in a horizontal line
+        #
+        expected = np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
+                             [0,1,0,0,0,0,0,0,0,0,0,0,0,1,0],
+                             [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
+                             [0,0,0,1,0,0,0,0,0,0,0,1,0,0,0],
+                             [0,0,0,0,1,1,1,1,1,1,1,0,0,0,0],
+                             [0,0,0,1,0,0,0,0,0,0,0,1,0,0,0],
+                             [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
+                             [0,1,0,0,0,0,0,0,0,0,0,0,0,1,0],
+                             [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]], bool)
+        result = medial_axis(image)
+        assert np.all(result == expected)
+
+    def test_01_02_hole(self):
+        '''Test skeletonize on a rectangle with a hole in the middle'''
+        image = np.zeros((9, 15), bool)
+        image[1:-1, 1:-1] = True
+        image[4, 4:-4] = False
+        expected = np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
+                             [0,1,0,0,0,0,0,0,0,0,0,0,0,1,0],
+                             [0,0,1,1,1,1,1,1,1,1,1,1,1,0,0],
+                             [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
+                             [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
+                             [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0],
+                             [0,0,1,1,1,1,1,1,1,1,1,1,1,0,0],
+                             [0,1,0,0,0,0,0,0,0,0,0,0,0,1,0],
+                             [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]],bool)
+        result = medial_axis(image)
+        assert np.all(result == expected)
+
 
 if __name__ == '__main__':
     np.testing.run_module_suite()
diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py
index 58cfaf32..34b76832 100644
--- a/skimage/transform/setup.py
+++ b/skimage/transform/setup.py
@@ -23,6 +23,7 @@ def configuration(parent_package='', top_path=None):
     config.add_extension('_project', sources=['_project.c'],
                          include_dirs=[get_numpy_include_dirs()])
 
+
     return config
 
 if __name__ == '__main__':

From 736802bad4400eec4e9a1f5c22725f39bde01091 Mon Sep 17 00:00:00 2001
From: emmanuelle <emmanuelle.gouillart@nsup.org>
Date: Sun, 23 Oct 2011 11:06:56 +0200
Subject: [PATCH 2/6] Medial axis transform (skeletonization) from cell
 profiler

 * medial_axis function in morphology.skeletonize
 * plot_medial_axis.py example for the gallery

Compared to the skeletonize algorithm in morphology.skeletonize,
medial_axis is faster because it processes all pixel in one pass.
However, the resulting skeleton has more branches, that may be unwanted
depending on the application.
---
 doc/examples/plot_medial_transform.py        |  70 ++++++++
 scikits/image/morphology/_cpmorphology2.pyx  | 153 ++++------------
 skimage/morphology/skeletonize.py            | 176 ++++++++-----------
 skimage/morphology/tests/test_skeletonize.py |   9 +
 4 files changed, 183 insertions(+), 225 deletions(-)
 create mode 100644 doc/examples/plot_medial_transform.py

diff --git a/doc/examples/plot_medial_transform.py b/doc/examples/plot_medial_transform.py
new file mode 100644
index 00000000..47fcd5bd
--- /dev/null
+++ b/doc/examples/plot_medial_transform.py
@@ -0,0 +1,70 @@
+"""
+===========================
+Medial axis skeletonization
+===========================
+
+The medial axis of an object is the set of all points having more than one 
+closest point on the object's boundary. It is often called the **topological
+skeleton**, because it is a 1-pixel wide skeleton of the object, with the same
+connectivity as the original object.
+
+Here, we use the medial axis transform to compute the width of the foreground
+objects. As the function ``medial_axis`` (``skimage.morphology.medial_axis``)
+returns the distance transform in addition to the medial axis (with the keyword
+argument ``return_distance=True``), it is possible to compute the distance to
+the background for all points of the medial axis with this function. This gives
+an estimate of the local width of the objects.
+
+For a skeleton with less branches, there exists another skeletonization 
+algorithm in ``skimage``: ``skimage.morphology.skeletonize``, that computes
+a skeleton by iterative morphological thinnings.
+"""
+
+import numpy as np
+from scipy import ndimage
+from scikits.image.morphology import medial_axis
+import matplotlib.pyplot as plt
+
+
+def microstructure(l=256):
+    """
+    Synthetic binary data: binary microstructure with blobs.
+
+    Parameters
+    ----------
+
+    l: int, optional 
+        linear size of the returned image
+
+    """
+    n = 5
+    x, y = np.ogrid[0:l, 0:l]
+    mask_outer = (x - l/2)**2 + (y - l/2)**2 < (l/2)**2
+    mask = np.zeros((l, l))
+    generator = np.random.RandomState(1)
+    points = l * generator.rand(2, n**2)
+    mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1
+    mask = ndimage.gaussian_filter(mask, sigma=l/(4.*n))
+    return mask > mask.mean()
+
+data = microstructure(l=64)
+
+# Compute the medial axis (skeleton) and the distance transform
+skel, distance = medial_axis(data, return_distance=True)
+
+# Distance to the background for pixels of the skeleton
+dist_on_skel = distance * skel
+
+plt.figure(figsize=(8, 4))
+plt.subplot(121)
+plt.imshow(data, cmap=plt.cm.gray, interpolation='nearest')
+plt.axis('off')
+plt.subplot(122)
+plt.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest')
+plt.contour(data, [0.5], colors='w')
+plt.axis('off')
+
+plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0,
+                    right=1)
+plt.show()
+
diff --git a/scikits/image/morphology/_cpmorphology2.pyx b/scikits/image/morphology/_cpmorphology2.pyx
index d7b76413..23f9736f 100644
--- a/scikits/image/morphology/_cpmorphology2.pyx
+++ b/scikits/image/morphology/_cpmorphology2.pyx
@@ -38,30 +38,49 @@ def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
                                 negative_indices=False, mode='c'] order,
                      np.ndarray[dtype=np.uint8_t, ndim=1,
                                 negative_indices=False, mode='c'] table):
-    '''Inner loop of skeletonize function
+    """
+    Inner loop of skeletonize function
 
-    result - on input, the image to be skeletonized, on output the skeletonized
-             image
-    i,j    - the coordinates of each foreground pixel in the image
-    order  - the index of each pixel, in the order of processing
-    table  - the 512-element lookup table of values after transformation
+    Parameters
+    ----------
+
+    result: ndarray of uint8
+        On input, the image to be skeletonized, on output the skeletonized
+        image.
+    i, j: ndarrays
+        The coordinates of each foreground pixel in the image
+    order: ndarray
+        The index of each pixel, in the order of processing (order[0] is
+        the first pixel to process, etc.)
+    table: ndarray
+        The 512-element lookup table of values after transformation 
+        (whether to keep or not each configuration in a binary 3x3 array)
+
+    Notes
+    -----
 
     The loop determines whether each pixel in the image can be removed without
     changing the Euler number of the image. The pixels are ordered by
     increasing distance from the background which means a point nearer to
     the quench-line of the brushfire will be evaluated later than a
     point closer to the edge.
-    '''
+
+    Note that the neighbourhood of a pixel may evolve before the loop 
+    arrives at this pixel. This is why it is possible to compute the 
+    skeleton in only one pass, thanks to an adapted ordering of the
+    pixels.
+    """
     cdef:
         np.int32_t accumulator
-        np.int32_t index,order_index
-        np.int32_t ii,jj
+        np.int32_t index, order_index
+        np.int32_t ii, jj
 
     for index in range(order.shape[0]):
         accumulator = 16
         order_index = order[index]
         ii = i[order_index]
         jj = j[order_index]
+        # Compute the configuration around the pixel
         if ii > 0:
             if jj > 0 and result[ii - 1, jj - 1]:
                 accumulator += 1
@@ -74,13 +93,14 @@ def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
             if jj < result.shape[1] - 1 and result[ii, jj + 1]:
                 accumulator += 32
             if ii < result.shape[0]-1:
-                if jj > 0 and result[ii+1,jj-1]:
+                if jj > 0 and result[ii + 1, jj - 1]:
                     accumulator += 64
-                if result[ii+1,jj]:
+                if result[ii + 1, jj]:
                     accumulator += 128
-                if jj < result.shape[1]-1 and result[ii+1,jj+1]:
+                if jj < result.shape[1] - 1 and result[ii + 1, jj + 1]:
                     accumulator += 256
-            result[ii,jj] = table[accumulator]
+            # Assign the value of table corresponding to the configuration
+            result[ii, jj] = table[accumulator]
 
 @cython.boundscheck(False)
 def table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2,
@@ -198,112 +218,7 @@ def table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2,
             indexer[i + 1, j_shape - 1] += 2
     return indexer
 
-@cython.boundscheck(False)
-def index_lookup(np.ndarray[dtype=np.int32_t, ndim=1, 
-                                negative_indices=False] index_i, 
-                 np.ndarray[dtype=np.int32_t, ndim=1, 
-                                negative_indices=False] index_j, 
-                 np.ndarray[dtype=np.uint32_t, ndim=2, 
-                                negative_indices=False] image,
-                 table_in, 
-                 iterations=None):
-    """
-    Perform a table lookup for only the indexed pixels
-    
-    For morphological operations that only convert 1 to 0, the set of
-    resulting pixels is always a subset of the input set. Therefore, when
-    repeating, it will be faster to operate only on the subsets especially
-    when the results are 1-d or 0-d objects.
-    
-    This function returns a new index_i and index_j array of the pixels
-    that survive the operation. The image is modified in-place to remove
-    the pixels that did not survive.
-    
-    index_i - an array of row indexes into the image.
-    index_j - a similarly-shaped array of column indexes.
-    image - the binary image: *NOTE* add a row and column of border values
-            to the original image to account for pixels on the edge of the
-            image.
-    iterations - # of iterations to do, default is "forever"
-    
-    The idea of index_lookup was taken from
-    http://blogs.mathworks.com/steve/2008/06/13/performance-optimization-for-applylut/
-    which, apparently, is how Matlab achieved its bwmorph speedup.
-    """
-    cdef:
-        np.ndarray[dtype=np.uint8_t, ndim=1, 
-            negative_indices=False] table = table_in.astype(np.uint8)
-        np.uint32_t center, hit_count, idx, indexer
-        np.int32_t idxi, idxj
-
-    if iterations == None:
-        # Worst case - remove one per iteration
-        iterations = len(index_i)
-
-    for i in range(iterations):
-        hit_count = len(index_i)
-        with nogil:
-            #
-            # For integer images (i.e., labels), a neighbor point is
-            # "background" if it doesn't match the central value. This
-            # lets adjacent labeled objects shrink independently of each
-            # other.
-            #
-            for 0 <= idx < hit_count:
-                idxi, idxj = index_i[idx], index_j[idx]
-                center = image[idxi, idxj]
-                indexer =   ((image[idxi - 1, idxj - 1] == center) * 1 +
-                             (image[idxi - 1, idxj]     == center) * 2 +
-                             (image[idxi - 1, idxj + 1] == center) * 4 +
-                             (image[idxi,     idxj - 1] == center) * 8 +
-                                    16 +
-                             (image[idxi,     idxj + 1] == center) * 32 +
-                             (image[idxi + 1, idxj - 1] == center) * 64 +
-                             (image[idxi + 1, idxj]     == center) * 128 +
-                             (image[idxi + 1, idxj + 1] == center) * 256)
-                if table[indexer] == 0:
-                    # mark for deletion
-                    index_i[idx] = -index_i[idx]
-
-            # remove marked pixels
-            for 0 <= idx < hit_count:
-                idxi, idxj = index_i[idx], index_j[idx]
-                if idxi < 0:
-                    image[-idxi, idxj] = 0
-
-        index_j = index_j[index_i >= 0]
-        index_i = index_i[index_i >= 0]
-        if len(index_i) == hit_count:
-            break
-
-    return (index_i, index_j)
-
-def prepare_for_index_lookup(image, border_value):
-    """
-    Return the index arrays of "1" pixels and an image with an added border
-    
-    The routine, index_lookup takes an array of i indexes, an array of
-    j indexes and an image guaranteed to be indexed successfully by 
-    index_<i,j>[:] +/- 1. This routine constructs an image with added border
-    pixels... evilly, the index, 0 - 1, lands on the border because of Python's
-    negative indexing convention.
-    """
-    if np.issubdtype(image.dtype, float):
-        image = image.astype(bool)
-    image_i, image_j = np.argwhere(image.astype(bool)).transpose().\
-                                    astype(np.int32) + 1
-    output_image = (np.ones(np.array(image.shape) + 2, image.dtype) \
-                    if border_value
-                    else np.zeros(np.array(image.shape) + 2, image.dtype))
-    output_image[1:image.shape[0] + 1, 1:image.shape[1] + 1] = image
-    return (image_i, image_j, output_image.astype(np.uint32))
 
 
-def extract_from_image_lookup(orig_image, index_i, index_j):
-    """
-    Extract only one pixel
-    """
-    output = np.zeros(orig_image.shape, orig_image.dtype)
-    output[index_i - 1, index_j - 1] = orig_image[index_i - 1, index_j - 1]
-    return output
+
 
diff --git a/skimage/morphology/skeletonize.py b/skimage/morphology/skeletonize.py
index 7d76f833..748d1b80 100644
--- a/skimage/morphology/skeletonize.py
+++ b/skimage/morphology/skeletonize.py
@@ -1,14 +1,11 @@
-"""Use an iterative thinning algorithm to find the skeletons of binary 
-objects in an image.
-
+"""
+Algorithms for computing the skeleton of a binary image
 """
 
 import numpy as np
 from scipy import ndimage
 
 from _cpmorphology2 import skeletonize_loop, table_lookup_index
-from _cpmorphology2 import extract_from_image_lookup, \
-            prepare_for_index_lookup, index_lookup
 
 # --------- Skeletonization by morphological thinning ---------
 
@@ -168,7 +165,7 @@ def medial_axis(image, mask=None, return_distance=False):
         If a mask is given, only those elements with a true value in `mask`
         are used for computing the medial axis.
 
-    return_distance; bool, optional
+    return_distance: bool, optional
         If true, the distance transform is returned as well as the skeleton.
 
     Returns
@@ -188,11 +185,24 @@ def medial_axis(image, mask=None, return_distance=False):
     Notes
     -----
     This algorithm computes the medial axis transform of an image
-    as the ridges of its distance transform. First, the distance transform
-    is computed, then the foreground (value of 1) points are ordered by
-    the distance transform. In order to reduce the image to its skeleton, 
-    a point is removed if it has more than one neighbor and if removing it
-    does not change the Euler number (the connectivity).
+    as the ridges of its distance transform. 
+    
+    The different steps of the algorithm are as follows
+     * A lookup table is used, that assigns 0 or 1 to each configuration of 
+       the 3x3 binary square, whether the central pixel should be removed
+       or kept. We want a point to be removed if it has more than one neighbor
+       and if removing it does not change the number of connected components.
+
+     * The distance transform to the background is computed, as well as
+       the cornerness of the pixel.
+
+     * The foreground (value of 1) points are ordered by
+       the distance transform, then the cornerness. 
+       
+     * A cython function is called to reduce the image to its skeleton. It 
+       processes pixels in the order determined at the previous step, and 
+       removes or not a pixel according to the lookup table. Because of the
+       ordering, it is possible to process all pixels in only one pass.
 
     Examples
     --------
@@ -223,48 +233,60 @@ def medial_axis(image, mask=None, return_distance=False):
         masked_image = image.astype(bool).copy()
         masked_image[~mask] = False
     #
-    # Lookup table - start with only positive pixels.
-    # Keep if # pixels in neighborhood is 2 or less
-    # Keep if removing the pixel results in a different connectivity
-    # table is independent of image
-    table = (_make_table(True, 
-                        np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], bool),
-                        np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], bool)) &
+    # Build lookup table - three conditions 
+    # 1. Keep only positive pixels (center_is_foreground array).
+    # AND
+    # 2. Keep if removing the pixel results in a different connectivity
+    # (if the number of connected components is different with and
+    # without the central pixel)
+    # OR
+    # 3. Keep if # pixels in neighbourhood is 2 or less
+    # Note that table is independent of image
+    center_is_foreground = (np.arange(512) & 2**4).astype(bool)
+    table = (center_is_foreground & 
             (np.array([ndimage.label(_pattern_of(index), eight_connect)[1] !=
                         ndimage.label(_pattern_of(index & ~ 2**4),
                                     eight_connect)[1]
                         for index in range(512)]) |
         np.array([np.sum(_pattern_of(index)) < 3 for index in range(512)])))
+    
+    # Build distance transform
     distance = ndimage.distance_transform_edt(masked_image)
     if return_distance:
         store_distance = distance.copy()
-    #
-        # The processing order along the edge is critical to the shape of the
+    
+    # Corners
+    # The processing order along the edge is critical to the shape of the
     # resulting skeleton: if you process a corner first, that corner will
     # be eroded and the skeleton will miss the arm from that corner. Pixels
     # with fewer neighbors are more "cornery" and should be processed last.
-    #
+    # We use a cornerness_table lookup table where the score of a 
+    # configuration is the number of background (0-value) pixels in the
+    # 3x3 neighbourhood
     cornerness_table = np.array([9 - np.sum(_pattern_of(index))
                                  for index in range(512)])
-    corner_score = _table_lookup(masked_image, cornerness_table, False, 1)
+    corner_score = _table_lookup(masked_image, cornerness_table)
+    
+    # Define arrays for inner loop
     i, j = np.mgrid[0:image.shape[0], 0:image.shape[1]]
     result = masked_image.copy()
     distance = distance[result]
     i = np.ascontiguousarray(i[result], np.int32)
     j = np.ascontiguousarray(j[result], np.int32)
     result = np.ascontiguousarray(result, np.uint8)
-    #
+    
+    # Determine the order in which pixels are processed.
     # We use a random # for tiebreaking. Assign each pixel in the image a
     # predictable, random # so that masking doesn't affect arbitrary choices
     # of skeletons
     #
-    # Why fix the seed? Should we pass a random number generator instead?
-    np.random.seed(0)
-    tiebreaker = np.random.permutation(np.arange(masked_image.sum()))
+    generator = np.random.RandomState(0)
+    tiebreaker = generator.permutation(np.arange(masked_image.sum()))
     order = np.lexsort((tiebreaker,
                         corner_score[masked_image],
                         distance))
     order = np.ascontiguousarray(order, np.int32)
+
     table = np.ascontiguousarray(table, np.uint8)
     # Remove pixels not belonging to the medial axis
     skeletonize_loop(result, i, j, order, table)
@@ -282,12 +304,12 @@ def _pattern_of(index):
     Return the pattern represented by an index value
     Byte decomposition of index
     """
-    return np.array([[index & 2**0,index & 2**1,index & 2**2],
-                     [index & 2**3,index & 2**4,index & 2**5],
-                     [index & 2**6,index & 2**7,index & 2**8]], bool)
+    return np.array([[index & 2**0, index & 2**1, index & 2**2],
+                     [index & 2**3, index & 2**4, index & 2**5],
+                     [index & 2**6, index & 2**7, index & 2**8]], bool)
 
 
-def _table_lookup(image, table, border_value, iterations = None):
+def _table_lookup(image, table):
     """
     Perform a morphological transform on an image, directed by its 
     neighbors
@@ -316,85 +338,27 @@ def _table_lookup(image, table, border_value, iterations = None):
     that evaluate to true. 
     """
     #
-    # Test for a table that never transforms a zero into a one:
+    # We accumulate into the indexer to get the index into the table
+    # at each point in the image
     #
-    center_is_zero = np.array([(x & 2**4) == 0 for x in range(2**9)])
-    use_index_trick = False
-    if (not np.any(table[center_is_zero]) and
-        (np.issubdtype(image.dtype, bool) or np.issubdtype(image.dtype, int))):
-        # Use the index trick
-        use_index_trick = True
-        invert = False
-    elif (np.all(table[~center_is_zero]) and np.issubdtype(image.dtype, bool)):
-        # All ones stay ones, invert the table and the image and do the trick
-        use_index_trick = True
-        invert = True
-        image = ~ image
-        # table index 0 -> 511 and the output is reversed
-        table = ~ table[511-np.arange(512)]
-        border_value = not border_value
-    if use_index_trick:
-        orig_image = image
-        index_i, index_j, image = prepare_for_index_lookup(image, border_value)
-        index_i, index_j = index_lookup(index_i, index_j,
-                                        image, table, iterations)
-        image = extract_from_image_lookup(orig_image, index_i, index_j)
-        if invert:
-            image = ~ image
-        return image
-    print(use_index_trick)
-    counter = 0
-    while counter != iterations:
-        counter += 1
-        #
-        # We accumulate into the indexer to get the index into the table
-        # at each point in the image
-        #
-        if image.shape[0] < 3 or image.shape[1] < 3:
-            image = image.astype(bool)
-            indexer = np.zeros(image.shape,int)
-            indexer[1:, 1:]   += image[:-1, :-1] * 2**0
-            indexer[1:, :]    += image[:-1, :] * 2**1
-            indexer[1:, :-1]  += image[:-1, 1:] * 2**2
-            
-            indexer[:, 1:]    += image[:, :-1] * 2**3
-            indexer[:, :]     += image[:, :] * 2**4
-            indexer[:, :-1]   += image[:, 1:] * 2**5
+    if image.shape[0] < 3 or image.shape[1] < 3:
+        image = image.astype(bool)
+        indexer = np.zeros(image.shape,int)
+        indexer[1:, 1:]   += image[:-1, :-1] * 2**0
+        indexer[1:, :]    += image[:-1, :] * 2**1
+        indexer[1:, :-1]  += image[:-1, 1:] * 2**2
         
-            indexer[:-1, 1:]  += image[1:, :-1] * 2**6
-            indexer[:-1, :]   += image[1:, :] * 2**7
-            indexer[:-1, :-1] += image[1:, 1:] * 2**8
-        else:
-            indexer = table_lookup_index(np.ascontiguousarray(image, np.uint8))
-        if border_value:
-            indexer[0,:]   |= 2**0 + 2**1 + 2**2
-            indexer[-1,:]  |= 2**6 + 2**7 + 2**8
-            indexer[:,0]   |= 2**0 + 2**3 + 2**6
-            indexer[:,-1]  |= 2**2 + 2**5 + 2**8
-        new_image = table[indexer]
-        if np.all(new_image == image):
-            break
-        image = new_image
+        indexer[:, 1:]    += image[:, :-1] * 2**3
+        indexer[:, :]     += image[:, :] * 2**4
+        indexer[:, :-1]   += image[:, 1:] * 2**5
+    
+        indexer[:-1, 1:]  += image[1:, :-1] * 2**6
+        indexer[:-1, :]   += image[1:, :] * 2**7
+        indexer[:-1, :-1] += image[1:, 1:] * 2**8
+    else:
+        indexer = table_lookup_index(np.ascontiguousarray(image, np.uint8))
+    image = table[indexer]
     return image
 
-def _make_table(value, pattern, care=np.ones((3,3),bool)):
-    '''Return a table suitable for table_lookup
-    
-    value - set all table entries matching "pattern" to "value", all others
-            to not "value"
-    pattern - a 3x3 boolean array with the pattern to match
-    care    - a 3x3 boolean array where each value is true if the pattern
-              must match at that position and false if we don't care if
-              the pattern matches at that position.
-    '''
-    def fn(index, p, i, j):
-        '''Return true if bit position "p" in index matches pattern'''
-        return ((((index & 2**p) > 0) == pattern[i, j]) or not care[i, j])
-    return np.array([value     
-                     if (fn(i, 0, 0, 0) and fn(i, 1, 0, 1) and fn(i, 2, 0, 2) 
-                     and fn(i, 3, 1, 0) and fn(i, 4, 1, 1) and fn(i, 5, 1, 2) 
-                     and fn(i, 6, 2, 0) and fn(i, 7, 2, 1) and fn(i, 8, 2, 2))
-                     else not value
-                     for i in range(512)], bool)
 
 
diff --git a/skimage/morphology/tests/test_skeletonize.py b/skimage/morphology/tests/test_skeletonize.py
index c1b69119..539e0ad8 100644
--- a/skimage/morphology/tests/test_skeletonize.py
+++ b/skimage/morphology/tests/test_skeletonize.py
@@ -121,6 +121,8 @@ class TestMedialAxis():
                              [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]], bool)
         result = medial_axis(image)
         assert np.all(result == expected)
+        result, distance = medial_axis(image, return_distance=True)
+        assert distance.max() == 4
 
     def test_01_02_hole(self):
         '''Test skeletonize on a rectangle with a hole in the middle'''
@@ -139,6 +141,13 @@ class TestMedialAxis():
         result = medial_axis(image)
         assert np.all(result == expected)
 
+    def test_narrow_image(self):
+        """Test skeletonize on a 1-pixel thin strip"""
+        image = np.zeros((1, 5), bool)
+        image[:, 1:-1] = True
+        result = medial_axis(image)
+        assert np.all(result == image)
+
 
 if __name__ == '__main__':
     np.testing.run_module_suite()

From ed1f05427e741175c17467fabbc20f0def6b18aa Mon Sep 17 00:00:00 2001
From: emmanuelle <emmanuelle.gouillart@nsup.org>
Date: Mon, 24 Oct 2011 21:31:24 +0200
Subject: [PATCH 3/6] ENH address comments on PR #65 made by Stefan and Neil

---
 doc/examples/plot_medial_transform.py         |  2 +-
 .../{_cpmorphology2.pyx => _skeletonize.pyx}  | 34 +++++------
 skimage/morphology/setup.py                   |  4 ++
 skimage/morphology/skeletonize.py             | 58 +++++++++++--------
 4 files changed, 52 insertions(+), 46 deletions(-)
 rename scikits/image/morphology/{_cpmorphology2.pyx => _skeletonize.pyx} (92%)

diff --git a/doc/examples/plot_medial_transform.py b/doc/examples/plot_medial_transform.py
index 47fcd5bd..af0b4c38 100644
--- a/doc/examples/plot_medial_transform.py
+++ b/doc/examples/plot_medial_transform.py
@@ -15,7 +15,7 @@ argument ``return_distance=True``), it is possible to compute the distance to
 the background for all points of the medial axis with this function. This gives
 an estimate of the local width of the objects.
 
-For a skeleton with less branches, there exists another skeletonization 
+For a skeleton with fewer branches, there exists another skeletonization 
 algorithm in ``skimage``: ``skimage.morphology.skeletonize``, that computes
 a skeleton by iterative morphological thinnings.
 """
diff --git a/scikits/image/morphology/_cpmorphology2.pyx b/scikits/image/morphology/_skeletonize.pyx
similarity index 92%
rename from scikits/image/morphology/_cpmorphology2.pyx
rename to scikits/image/morphology/_skeletonize.pyx
index 23f9736f..0b2f8a8c 100644
--- a/scikits/image/morphology/_cpmorphology2.pyx
+++ b/scikits/image/morphology/_skeletonize.pyx
@@ -14,21 +14,9 @@ import numpy as np
 cimport numpy as np
 cimport cython
 
-cdef extern from "Python.h":
-    ctypedef int Py_intptr_t
-
-cdef extern from "numpy/arrayobject.h":
-    ctypedef class numpy.ndarray [object PyArrayObject]:
-        cdef char *data
-        cdef Py_intptr_t *dimensions
-        cdef Py_intptr_t *strides
-    cdef void import_array()
-    cdef int  PyArray_ITEMSIZE(np.ndarray)
-
-import_array()
 
 @cython.boundscheck(False)
-def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, 
+def _skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, 
                                 negative_indices=False, mode='c'] result,
                      np.ndarray[dtype=np.int32_t, ndim=1,
                                 negative_indices=False, mode='c'] i,
@@ -44,15 +32,18 @@ def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
     Parameters
     ----------
 
-    result: ndarray of uint8
+    result : ndarray of uint8
         On input, the image to be skeletonized, on output the skeletonized
         image.
-    i, j: ndarrays
+
+    i, j : ndarrays
         The coordinates of each foreground pixel in the image
-    order: ndarray
+   
+   order : ndarray
         The index of each pixel, in the order of processing (order[0] is
         the first pixel to process, etc.)
-    table: ndarray
+   
+   table : ndarray
         The 512-element lookup table of values after transformation 
         (whether to keep or not each configuration in a binary 3x3 array)
 
@@ -103,21 +94,24 @@ def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
             result[ii, jj] = table[accumulator]
 
 @cython.boundscheck(False)
-def table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2,
+def _table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2,
                                   negative_indices=False, mode='c'] image):
     """
     Return an index into a table per pixel of a binary image
 
     Take the sum of true neighborhood pixel values where the neighborhood
-    looks like this:
+    looks like this::
+
      1   2   4
      8  16  32
     64 128 256
 
-    This code could be replaced by a convolution with the kernel:
+    This code could be replaced by a convolution with the kernel::
+
     256 128 64
      32  16  8
       4   2  1
+    
     but this runs about twice as fast because of inlining and the
     hardwired kernel.
     """
diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py
index 9d8b5527..46010989 100644
--- a/skimage/morphology/setup.py
+++ b/skimage/morphology/setup.py
@@ -14,6 +14,7 @@ def configuration(parent_package='', top_path=None):
     cython(['ccomp.pyx'], working_path=base_path)
     cython(['cmorph.pyx'], working_path=base_path)
     cython(['_watershed.pyx'], working_path=base_path)
+    cython(['_skeletonize.pyx'], working_path=base_path)
 
     config.add_extension('ccomp', sources=['ccomp.c'],
                          include_dirs=[get_numpy_include_dirs()])
@@ -21,6 +22,9 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_watershed', sources=['_watershed.c'],
                          include_dirs=[get_numpy_include_dirs()])
+    config.add_extension('_skeletonize', sources=['_skeletonize.c'],
+                         include_dirs=[get_numpy_include_dirs()])
+
 
     return config
 
diff --git a/skimage/morphology/skeletonize.py b/skimage/morphology/skeletonize.py
index 748d1b80..b7c54c9e 100644
--- a/skimage/morphology/skeletonize.py
+++ b/skimage/morphology/skeletonize.py
@@ -5,7 +5,7 @@ Algorithms for computing the skeleton of a binary image
 import numpy as np
 from scipy import ndimage
 
-from _cpmorphology2 import skeletonize_loop, table_lookup_index
+from _skeletonize import _skeletonize_loop, _table_lookup_index
 
 # --------- Skeletonization by morphological thinning ---------
 
@@ -149,7 +149,7 @@ def skeletonize(image):
 
 # --------- Skeletonization by medial axis transform --------
 
-eight_connect = ndimage.generate_binary_structure(2, 2)
+_eight_connect = ndimage.generate_binary_structure(2, 2)
 
 
 def medial_axis(image, mask=None, return_distance=False):
@@ -159,22 +159,22 @@ def medial_axis(image, mask=None, return_distance=False):
     Parameters
     ----------
 
-    image: binary ndarray 
+    image : binary ndarray 
     
-    mask: binary ndarray, optional
+    mask : binary ndarray, optional
         If a mask is given, only those elements with a true value in `mask`
         are used for computing the medial axis.
 
-    return_distance: bool, optional
+    return_distance : bool, optional
         If true, the distance transform is returned as well as the skeleton.
 
     Returns
     -------
 
-    out: ndarray of bools
+    out : ndarray of bools
         Medial axis transform of the image
 
-    dist: ndarray of ints
+    dist : ndarray of ints
         Distance transform of the image (only returned if `return_distance`
         is True)
    
@@ -201,8 +201,9 @@ def medial_axis(image, mask=None, return_distance=False):
        
      * A cython function is called to reduce the image to its skeleton. It 
        processes pixels in the order determined at the previous step, and 
-       removes or not a pixel according to the lookup table. Because of the
-       ordering, it is possible to process all pixels in only one pass.
+       removes or maintains a pixel according to the lookup table. Because 
+       of the ordering, it is possible to process all pixels in only one 
+       pass.
 
     Examples
     --------
@@ -226,7 +227,7 @@ def medial_axis(image, mask=None, return_distance=False):
            [0, 0, 0, 0, 0, 0, 0]], dtype=uint8)
 
     """
-    global eight_connect
+    global _eight_connect
     if mask is None:
         masked_image = image.astype(np.bool)
     else:
@@ -243,12 +244,17 @@ def medial_axis(image, mask=None, return_distance=False):
     # 3. Keep if # pixels in neighbourhood is 2 or less
     # Note that table is independent of image
     center_is_foreground = (np.arange(512) & 2**4).astype(bool)
-    table = (center_is_foreground & 
-            (np.array([ndimage.label(_pattern_of(index), eight_connect)[1] !=
+    table = (center_is_foreground  # condition 1.
+                &
+            (np.array([ndimage.label(_pattern_of(index), _eight_connect)[1] !=
                         ndimage.label(_pattern_of(index & ~ 2**4),
-                                    eight_connect)[1]
-                        for index in range(512)]) |
-        np.array([np.sum(_pattern_of(index)) < 3 for index in range(512)])))
+                                    _eight_connect)[1]
+                        for index in range(512)]) # condition 2 
+                | 
+        np.array([np.sum(_pattern_of(index)) < 3 for index in range(512)]))
+        # condition 3
+            )
+        
     
     # Build distance transform
     distance = ndimage.distance_transform_edt(masked_image)
@@ -289,7 +295,7 @@ def medial_axis(image, mask=None, return_distance=False):
 
     table = np.ascontiguousarray(table, np.uint8)
     # Remove pixels not belonging to the medial axis
-    skeletonize_loop(result, i, j, order, table)
+    _skeletonize_loop(result, i, j, order, table)
 
     result = result.astype(bool)
     if not mask is None:
@@ -316,20 +322,22 @@ def _table_lookup(image, table):
     
     Parameters
     ----------
-    image - a binary image
-    table - a 512-element table giving the transform of each pixel given
-            the values of that pixel and its 8-connected neighbors.
-    border_value - the value of pixels beyond the border of the image.
-                   This should test as True or False.
+    image : ndarray
+        A binary image
+    table : ndarray
+        A 512-element table giving the transform of each pixel given
+        the values of that pixel and its 8-connected neighbors.
+    border_value : bool
+        The value of pixels beyond the border of the image.
     
     Returns
     -------
-    result: ndarray of same shape as `image`
+    result : ndarray of same shape as `image`
         Transformed image
 
     Notes
     -----
-    The pixels are numbered like this:
+    The pixels are numbered like this::
     
     0 1 2
     3 4 5
@@ -343,7 +351,7 @@ def _table_lookup(image, table):
     #
     if image.shape[0] < 3 or image.shape[1] < 3:
         image = image.astype(bool)
-        indexer = np.zeros(image.shape,int)
+        indexer = np.zeros(image.shape, int)
         indexer[1:, 1:]   += image[:-1, :-1] * 2**0
         indexer[1:, :]    += image[:-1, :] * 2**1
         indexer[1:, :-1]  += image[:-1, 1:] * 2**2
@@ -356,7 +364,7 @@ def _table_lookup(image, table):
         indexer[:-1, :]   += image[1:, :] * 2**7
         indexer[:-1, :-1] += image[1:, 1:] * 2**8
     else:
-        indexer = table_lookup_index(np.ascontiguousarray(image, np.uint8))
+        indexer = _table_lookup_index(np.ascontiguousarray(image, np.uint8))
     image = table[indexer]
     return image
 

From 25e986d545b712b482a7b67e370aebd18e306a70 Mon Sep 17 00:00:00 2001
From: emmanuelle <emmanuelle.gouillart@nsup.org>
Date: Mon, 24 Oct 2011 23:27:29 +0200
Subject: [PATCH 4/6] Moved cython file that had not been automaticall moved by
 rename scikits/image -> skimage

---
 {scikits/image => skimage}/morphology/_skeletonize.pyx | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename {scikits/image => skimage}/morphology/_skeletonize.pyx (100%)

diff --git a/scikits/image/morphology/_skeletonize.pyx b/skimage/morphology/_skeletonize.pyx
similarity index 100%
rename from scikits/image/morphology/_skeletonize.pyx
rename to skimage/morphology/_skeletonize.pyx

From 39a0e3577c7206e0ebb3aee44dfa9a06a57828fa Mon Sep 17 00:00:00 2001
From: emmanuelle <emmanuelle.gouillart@nsup.org>
Date: Mon, 24 Oct 2011 23:28:21 +0200
Subject: [PATCH 5/6] scikits/image -> skimage in examples

---
 doc/examples/plot_medial_transform.py | 2 +-
 doc/examples/plot_skeleton.py         | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/examples/plot_medial_transform.py b/doc/examples/plot_medial_transform.py
index af0b4c38..1421b4ce 100644
--- a/doc/examples/plot_medial_transform.py
+++ b/doc/examples/plot_medial_transform.py
@@ -22,7 +22,7 @@ a skeleton by iterative morphological thinnings.
 
 import numpy as np
 from scipy import ndimage
-from scikits.image.morphology import medial_axis
+from skimage.morphology import medial_axis
 import matplotlib.pyplot as plt
 
 
diff --git a/doc/examples/plot_skeleton.py b/doc/examples/plot_skeleton.py
index a61bbbdd..7cb33d95 100644
--- a/doc/examples/plot_skeleton.py
+++ b/doc/examples/plot_skeleton.py
@@ -15,8 +15,8 @@ results. The input is a 2D ndarray, with either boolean or integer elements.
 In the case of boolean, 'True' indicates foreground, and for integer arrays,
 the foreground is 1's.  
 """
-from scikits.image.morphology import skeletonize
-from scikits.image.draw import draw
+from skimage.morphology import skeletonize
+from skimage.draw import draw
 import numpy as np
 import matplotlib.pyplot as plt
 

From f8a51a4e3e3284922db011d19ada0768a5ab5110 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Mon, 24 Oct 2011 17:23:39 -0700
Subject: [PATCH 6/6] DOC: Minor spacing fixes in docstrings.

---
 skimage/morphology/_skeletonize.pyx | 40 ++++++++++++-----------------
 skimage/morphology/skeletonize.py   |  9 +++----
 2 files changed, 21 insertions(+), 28 deletions(-)

diff --git a/skimage/morphology/_skeletonize.pyx b/skimage/morphology/_skeletonize.pyx
index 0b2f8a8c..ff5fcdf2 100644
--- a/skimage/morphology/_skeletonize.pyx
+++ b/skimage/morphology/_skeletonize.pyx
@@ -1,5 +1,4 @@
-'''_cpmorphology2.pyx - support routines for cpmorphology in Cython
-
+'''
 Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
 Website: http://www.cellprofiler.org
 
@@ -39,11 +38,11 @@ def _skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
     i, j : ndarrays
         The coordinates of each foreground pixel in the image
    
-   order : ndarray
+    order : ndarray
         The index of each pixel, in the order of processing (order[0] is
         the first pixel to process, etc.)
    
-   table : ndarray
+    table : ndarray
         The 512-element lookup table of values after transformation 
         (whether to keep or not each configuration in a binary 3x3 array)
 
@@ -102,15 +101,15 @@ def _table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2,
     Take the sum of true neighborhood pixel values where the neighborhood
     looks like this::
 
-     1   2   4
-     8  16  32
-    64 128 256
+       1   2   4
+       8  16  32
+      64 128 256
 
     This code could be replaced by a convolution with the kernel::
 
-    256 128 64
-     32  16  8
-      4   2  1
+      256 128 64
+       32  16  8
+        4   2  1
     
     but this runs about twice as fast because of inlining and the
     hardwired kernel.
@@ -182,37 +181,32 @@ def _table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2,
     for j in range(1, j_shape - 1):
         if image[0, j]:
             indexer[0, j - 1] += 32
-            indexer[0, j]   += 16
+            indexer[0, j] += 16
             indexer[0, j + 1] += 8
             indexer[1, j - 1] += 4
-            indexer[1, j]   += 2
+            indexer[1, j] += 2
             indexer[1, j + 1] += 1
         if image[i_shape - 1, j]:
             indexer[i_shape - 2, j - 1] += 256
-            indexer[i_shape - 2, j]   += 128
+            indexer[i_shape - 2, j] += 128
             indexer[i_shape - 2, j + 1] += 64
             indexer[i_shape - 1, j - 1] += 32
-            indexer[i_shape - 1, j]   += 16
+            indexer[i_shape - 1, j] += 16
             indexer[i_shape - 1, j + 1] += 8
 
     for i in range(1, i_shape - 1):
         if image[i, 0]:
             indexer[i - 1, 0] += 128
-            indexer[i, 0]   += 16
+            indexer[i, 0] += 16
             indexer[i + 1, 0] += 2
             indexer[i - 1, 1] += 64
-            indexer[i, 1]   += 8
+            indexer[i, 1] += 8
             indexer[i + 1, 1] += 1
         if image[i, j_shape - 1]:
             indexer[i - 1, j_shape - 2] += 256
-            indexer[i, j_shape - 2]   += 32
+            indexer[i, j_shape - 2] += 32
             indexer[i + 1, j_shape - 2] += 4
             indexer[i - 1, j_shape - 1] += 128
-            indexer[i, j_shape - 1]   += 16
+            indexer[i, j_shape - 1] += 16
             indexer[i + 1, j_shape - 1] += 2
     return indexer
-
-
-
-
-
diff --git a/skimage/morphology/skeletonize.py b/skimage/morphology/skeletonize.py
index b7c54c9e..5ddf2c4c 100644
--- a/skimage/morphology/skeletonize.py
+++ b/skimage/morphology/skeletonize.py
@@ -339,9 +339,10 @@ def _table_lookup(image, table):
     -----
     The pixels are numbered like this::
     
-    0 1 2
-    3 4 5
-    6 7 8
+      0 1 2
+      3 4 5
+      6 7 8
+
     The index at a pixel is the sum of 2**<pixel-number> for pixels
     that evaluate to true. 
     """
@@ -368,5 +369,3 @@ def _table_lookup(image, table):
     image = table[indexer]
     return image
 
-
-