diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index aba19f3b..0cda885d 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -78,3 +78,6 @@
 - Christoph Gohlke
   Windows packaging and Python 3 compatibility.
 
+- Neil Yager
+  Skeletonization.
+  
\ No newline at end of file
diff --git a/doc/examples/plot_skeleton.py b/doc/examples/plot_skeleton.py
new file mode 100644
index 00000000..a61bbbdd
--- /dev/null
+++ b/doc/examples/plot_skeleton.py
@@ -0,0 +1,63 @@
+"""
+===========
+Skeletonize
+===========
+
+Skeletonization reduces binary objects to 1 pixel wide representations. This
+can be useful for feature extraction, and/or representing an object's topology.
+
+The algorithm works by making successive passes of the image. On each pass, 
+border pixels are identified and removed on the condition that they do not
+break the connectivity of the corresponding object. 
+
+This module provides an example of calling the routine and displaying the
+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
+import numpy as np
+import matplotlib.pyplot as plt
+
+# an empty image
+image = np.zeros((400, 400))
+
+# foreground object 1
+image[10:-10, 10:100] = 1
+image[-100:-10, 10:-10] = 1
+image[10:-10, -100:-10] = 1
+
+# foreground object 2
+rs, cs = draw.bresenham(250, 150, 10, 280)
+for i in range(10): image[rs+i, cs] = 1 
+rs, cs = draw.bresenham(10, 150, 250, 280)
+for i in range(20): image[rs+i, cs] = 1 
+
+# foreground object 3
+ir, ic = np.indices(image.shape)
+circle1 = (ic - 135)**2 + (ir - 150)**2 < 30**2
+circle2 = (ic - 135)**2 + (ir - 150)**2 < 20**2
+image[circle1] = 1
+image[circle2] = 0
+
+# perform skeletonization
+skeleton = skeletonize(image)
+
+# display results
+plt.figure(figsize=(10,6))
+
+plt.subplot(121)
+plt.imshow(image, cmap=plt.cm.gray)
+plt.axis('off')
+plt.title('original', fontsize=20)
+
+plt.subplot(122)
+plt.imshow(skeleton, cmap=plt.cm.gray)
+plt.axis('off')
+plt.title('skeleton', fontsize=20)
+
+plt.subplots_adjust(wspace=0.02, hspace=0.02, top=0.98, 
+                    bottom=0.02, left=0.02, right=0.98)
+
+plt.show()
diff --git a/scikits/image/data/bw_text.png b/scikits/image/data/bw_text.png
new file mode 100644
index 00000000..fd1a4c7c
Binary files /dev/null and b/scikits/image/data/bw_text.png differ
diff --git a/scikits/image/data/bw_text_skeleton.npy b/scikits/image/data/bw_text_skeleton.npy
new file mode 100644
index 00000000..9492cb64
Binary files /dev/null and b/scikits/image/data/bw_text_skeleton.npy differ
diff --git a/scikits/image/morphology/__init__.py b/scikits/image/morphology/__init__.py
index 03cf49e0..d57cc859 100644
--- a/scikits/image/morphology/__init__.py
+++ b/scikits/image/morphology/__init__.py
@@ -2,3 +2,4 @@ from grey import *
 from selem import *
 from .ccomp import label
 from watershed import watershed, is_local_maximum
+from skeletonize import skeletonize
diff --git a/scikits/image/morphology/skeletonize.py b/scikits/image/morphology/skeletonize.py
new file mode 100644
index 00000000..ea478a9a
--- /dev/null
+++ b/scikits/image/morphology/skeletonize.py
@@ -0,0 +1,141 @@
+"""Use an iterative thinning algorithm to find the skeletons of binary 
+objects in an image.
+
+"""
+
+import numpy as np
+from scipy.ndimage import correlate
+
+def skeletonize(image):
+    """Return the skeleton of a binary image.
+    
+    Thinning is used to reduce each connected component in a binary image
+    to a single-pixel wide skeleton. 
+     
+    Parameters
+    ----------
+    image : numpy.ndarray 
+        A binary image containing the objects to be skeletonized. '1' 
+        represents foreground, and '0' represents background. It 
+        also accepts arrays of boolean values where True is foreground.
+    
+    Returns
+    -------
+    skeleton : ndarray
+        A matrix containing the thinned image.
+    
+    Notes
+    -----
+    The algorithm [1] works by making successive passes of the image, 
+    removing pixels on object borders. This continues until no
+    more pixels can be removed.  The image is correlated with a
+    mask that assigns each pixel a number in the range [0...255]
+    corresponding to each possible pattern of its 8 neighbouring
+    pixels. A look up table is then used to assign the pixels a
+    value of 0, 1, 2 or 3, which are selectively removed during
+    the iterations.
+    
+    Note that this algorithm will give different results than a 
+    medial axis transform, which is also often referred to as
+    "skeletonization".   
+   
+    References
+    ----------
+    .. [1] A fast parallel algorithm for thinning digital patterns, 
+       T. Y. ZHANG and C. Y. SUEN, Communications of the ACM,
+       March 1984, Volume 27, Number 3
+
+    
+    Examples
+    --------
+    >>> X, Y = np.ogrid[0:9, 0:9]
+    >>> ellipse = (1./3 * (X - 4)**2 + (Y - 4)**2 < 3**2).astype(np.uint8)
+    >>> ellipse
+    array([[0, 0, 0, 1, 1, 1, 0, 0, 0],
+           [0, 0, 1, 1, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 1, 0, 0],
+           [0, 0, 1, 1, 1, 1, 1, 0, 0],
+           [0, 0, 0, 1, 1, 1, 0, 0, 0]], dtype=uint8)
+    >>> skel = skeletonize(ellipse)
+    >>> skel
+    array([[0, 0, 0, 0, 0, 0, 0, 0, 0],
+           [0, 0, 0, 0, 0, 0, 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, 1, 0, 0, 0, 0],
+           [0, 0, 0, 0, 1, 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, 0, 0, 0, 0, 0, 0]], dtype=uint8)
+           
+    """
+    # look up table - there is one entry for each of the 2^8=256 possible
+    # combinations of 8 binary neighbours. 1's, 2's and 3's are candidates
+    # for removal at each iteration of the algorithm.
+    lut = [ 0,0,0,1,0,0,1,3,0,0,3,1,1,0,1,3,0,0,0,0,0,0,0,0,2,0,2,0,3,0,3,3,
+            0,0,0,0,0,0,0,0,3,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,2,0,0,0,3,0,2,2,
+            0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+            2,0,0,0,0,0,0,0,2,0,0,0,2,0,0,0,3,0,0,0,0,0,0,0,3,0,0,0,3,0,2,0,
+            0,1,3,1,0,0,1,3,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,
+            3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+            2,3,1,3,0,0,1,3,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+            2,3,0,1,0,0,0,1,0,0,0,0,0,0,0,0,3,3,0,1,0,0,0,0,2,2,0,0,2,0,0,0]
+
+    # convert to unsigned int (this should work for boolean values)
+    skeleton = np.array(image).astype(np.uint8)
+
+    # check some properties of the input image:
+    #  - 2D
+    #  - binary image with only 0's and 1's 
+    if skeleton.ndim != 2:
+        raise ValueError('Skeletonize requires a 2D array')    
+    if not np.all(np.in1d(skeleton.flat, (0, 1))):
+        raise ValueError('Image contains values other than 0 and 1')
+    
+    # create the mask that will assign a unique value based on the
+    #  arrangement of neighbouring pixels
+    mask = np.array([[  1,  2,  4],
+                     [128,  0,  8],
+                     [ 64, 32, 16]], np.uint8)
+
+    pixelRemoved = True
+    while pixelRemoved:
+        pixelRemoved = False;
+
+        # assign each pixel a unique value based on its foreground neighbours
+        neighbours = correlate(skeleton, mask, mode='constant')
+        
+        # ignore background
+        neighbours *= skeleton
+        
+        # use LUT to categorize each foreground pixel as a 0, 1, 2 or 3
+        codes = np.take(lut, neighbours)
+        
+        # pass 1 - remove the 1's and 3's
+        code_mask = (codes == 1)
+        if np.any(code_mask): 
+            pixelRemoved = True
+            skeleton[code_mask] = 0
+        code_mask = (codes == 3)
+        if np.any(code_mask): 
+            pixelRemoved = True
+            skeleton[code_mask] = 0
+
+        # pass 2 - remove the 2's and 3's
+        neighbours = correlate(skeleton, mask, mode='constant')
+        neighbours *= skeleton
+        codes = np.take(lut, neighbours)
+        code_mask = (codes == 2)
+        if np.any(code_mask): 
+            pixelRemoved = True
+            skeleton[code_mask] = 0
+        code_mask = (codes == 3)
+        if np.any(code_mask): 
+            pixelRemoved = True
+            skeleton[code_mask] = 0
+
+    return skeleton
diff --git a/scikits/image/morphology/tests/test_skeletonize.py b/scikits/image/morphology/tests/test_skeletonize.py
new file mode 100644
index 00000000..aada8a46
--- /dev/null
+++ b/scikits/image/morphology/tests/test_skeletonize.py
@@ -0,0 +1,95 @@
+import numpy as np
+from scikits.image.morphology import skeletonize
+import numpy.testing
+from scikits.image.draw import draw
+from scipy.ndimage import correlate
+from scikits.image.io import imread
+from scikits.image import data_dir
+import os.path
+
+class TestSkeletonize():
+    def test_skeletonize_no_foreground(self):
+        im = np.zeros((5,5))
+        result = skeletonize(im)
+        numpy.testing.assert_array_equal(result, np.zeros((5,5)))
+    
+    def test_skeletonize_wrong_dim1(self):
+        im = np.zeros((5))
+        numpy.testing.assert_raises(ValueError, skeletonize, im)    
+
+    def test_skeletonize_wrong_dim2(self):
+        im = np.zeros((5, 5, 5))
+        numpy.testing.assert_raises(ValueError, skeletonize, im)    
+
+    def test_skeletonize_not_binary(self):
+        im = np.zeros((5, 5))
+        im[0, 0] = 1
+        im[0, 1] = 2
+        numpy.testing.assert_raises(ValueError, skeletonize, im)  
+        
+    def test_skeletonize_unexpected_value(self):
+        im = np.zeros((5, 5))
+        im[0, 0] = 2
+        numpy.testing.assert_raises(ValueError, skeletonize, im)          
+        
+    def test_skeletonize_all_foreground(self):
+        im = np.ones((3,4))
+        result = skeletonize(im)
+    
+    def test_skeletonize_single_point(self):
+        im = np.zeros((5, 5), np.uint8)
+        im[3, 3] = 1
+        result = skeletonize(im)
+        numpy.testing.assert_array_equal(result, im)
+    
+    def test_skeletonize_already_thinned(self):
+        im = np.zeros((5, 5), np.uint8)
+        im[3,1:-1] = 1
+        im[2, -1] = 1
+        im[4, 0] = 1
+        result = skeletonize(im)
+        numpy.testing.assert_array_equal(result, im)
+    
+    def test_skeletonize_output(self):
+        im = imread(os.path.join(data_dir, "bw_text.png"), as_grey=True)
+        
+        # make black the foreground
+        im = (im==0)
+        result = skeletonize(im)
+        
+        expected = np.load(os.path.join(data_dir, "bw_text_skeleton.npy"))
+        numpy.testing.assert_array_equal(result, expected)
+        
+    
+    def test_skeletonize_num_neighbours(self):
+        # an empty image
+        image = np.zeros((300, 300))
+        
+        # foreground object 1
+        image[10:-10, 10:100] = 1
+        image[-100:-10, 10:-10] = 1
+        image[10:-10, -100:-10] = 1
+        
+        # foreground object 2
+        rs, cs = draw.bresenham(250, 150, 10, 280)
+        for i in range(10): image[rs+i, cs] = 1 
+        rs, cs = draw.bresenham(10, 150, 250, 280)
+        for i in range(20): image[rs+i, cs] = 1 
+        
+        # foreground object 3
+        ir, ic = np.indices(image.shape)
+        circle1 = (ic - 135)**2 + (ir - 150)**2 < 30**2
+        circle2 = (ic - 135)**2 + (ir - 150)**2 < 20**2
+        image[circle1] = 1
+        image[circle2] = 0
+        result = skeletonize(image)
+        
+        # there should never be a 2x2 block of foreground pixels in a skeleton
+        mask = np.array([[1,  1],
+                         [1,  1]], np.uint8)        
+        blocks = correlate(result, mask, mode='constant')
+        assert not numpy.any(blocks == 4)
+        
+
+if __name__ == '__main__':
+    np.testing.run_module_suite()