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__':