diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index aba19f3b..dbcb17d3 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -78,3 +78,5 @@
 - 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
index 68958c6d..a61bbbdd 100644
--- a/doc/examples/plot_skeleton.py
+++ b/doc/examples/plot_skeleton.py
@@ -3,7 +3,17 @@
 Skeletonize
 ===========
 
-An example of thinning a binary image using 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
@@ -32,7 +42,7 @@ image[circle1] = 1
 image[circle2] = 0
 
 # perform skeletonization
-skeleton = skeletonize.skeletonize(image)
+skeleton = skeletonize(image)
 
 # display results
 plt.figure(figsize=(10,6))
@@ -50,4 +60,4 @@ 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()
\ No newline at end of file
+plt.show()
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
index 853aba3e..dc04e7fa 100644
--- a/scikits/image/morphology/skeletonize.py
+++ b/scikits/image/morphology/skeletonize.py
@@ -1,7 +1,5 @@
 """skeletonize.py - Use an iterative thinning algorithm to find the
                     skeletons of binary objects in an image.
-
-Original author: Neil Yager
 """
 
 import numpy as np
@@ -11,16 +9,16 @@ from .. import util
 def skeletonize(image):
     """
     Return a single pixel wide skeleton of all connected components 
-     in a binary image. 
+    in a binary image. 
      
     The algorithm 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. 
+    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. 
 
     Parameters
     ----------
@@ -34,9 +32,9 @@ def skeletonize(image):
     -----
     
     This implementation gives different results than a medial 
-     axis transforrmation, which can be can be implemented using  
-     morphological operations. This implementation is generally much
-     faster.
+    axis transformation, which can be can be implemented using  
+    morphological operations. This implementation is generally much
+    faster.
     
     Returns 
     -------
@@ -55,7 +53,9 @@ def skeletonize(image):
     --------
     """
     
-    # look up table 
+    # 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,
@@ -73,9 +73,8 @@ def skeletonize(image):
     #  - binary image with only 0's and 1's 
     if skeleton.ndim != 2:
         raise ValueError('Skeletonize requires a 2D array')    
-    for val in np.unique(skeleton):
-        if val not in [0, 1]:
-            raise ValueError('Invalid value in the image: %d'%(val))
+    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
@@ -87,26 +86,36 @@ def skeletonize(image):
     while pixelRemoved:
         pixelRemoved = False;
 
-        # pass 1 - remove the 1's and 3's
+        # assign each pixel a unique value based on its foreground neighbours
         neighbours = correlate(skeleton, mask, mode='constant')
+        
+        # ignore background
         neighbours[skeleton == 0] = 0
+        
+        # use LUT to categorize each foreground pixel as a 0, 1, 2 or 3
         codes = np.take(lut, neighbours)
-        if np.any(codes == 1): 
+        
+        # pass 1 - remove the 1's and 3's
+        code_mask = (codes == 1)
+        if np.any(code_mask): 
             pixelRemoved = True
-            skeleton[codes == 1] = 0
-        if np.any(codes == 3): 
+            skeleton[code_mask] = 0
+        code_mask = (codes == 3)
+        if np.any(code_mask): 
             pixelRemoved = True
-            skeleton[codes == 3] = 0
+            skeleton[code_mask] = 0
 
         # pass 2 - remove the 2's and 3's
         neighbours = correlate(skeleton, mask, mode='constant')
         neighbours[skeleton == 0] = 0
         codes = np.take(lut, neighbours)
-        if np.any(codes == 2): 
+        code_mask = (codes == 2)
+        if np.any(code_mask): 
             pixelRemoved = True
-            skeleton[codes == 2] = 0
-        if np.any(codes == 3): 
+            skeleton[code_mask] = 0
+        code_mask = (codes == 3)
+        if np.any(code_mask): 
             pixelRemoved = True
-            skeleton[codes == 3] = 0
+            skeleton[code_mask] = 0
 
     return skeleton
diff --git a/scikits/image/morphology/tests/test_skeletonize.py b/scikits/image/morphology/tests/test_skeletonize.py
index 841fcfbe..aada8a46 100644
--- a/scikits/image/morphology/tests/test_skeletonize.py
+++ b/scikits/image/morphology/tests/test_skeletonize.py
@@ -1,4 +1,3 @@
-import unittest
 import numpy as np
 from scikits.image.morphology import skeletonize
 import numpy.testing
@@ -8,39 +7,39 @@ from scikits.image.io import imread
 from scikits.image import data_dir
 import os.path
 
-class TestSkeletonize(unittest.TestCase):
+class TestSkeletonize():
     def test_skeletonize_no_foreground(self):
         im = np.zeros((5,5))
-        result = skeletonize.skeletonize(im)
+        result = skeletonize(im)
         numpy.testing.assert_array_equal(result, np.zeros((5,5)))
     
     def test_skeletonize_wrong_dim1(self):
         im = np.zeros((5))
-        self.assertRaises(ValueError, skeletonize.skeletonize, im)    
+        numpy.testing.assert_raises(ValueError, skeletonize, im)    
 
     def test_skeletonize_wrong_dim2(self):
         im = np.zeros((5, 5, 5))
-        self.assertRaises(ValueError, skeletonize.skeletonize, im)    
+        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
-        self.assertRaises(ValueError, skeletonize.skeletonize, im)  
+        numpy.testing.assert_raises(ValueError, skeletonize, im)  
         
     def test_skeletonize_unexpected_value(self):
         im = np.zeros((5, 5))
         im[0, 0] = 2
-        self.assertRaises(ValueError, skeletonize.skeletonize, im)          
+        numpy.testing.assert_raises(ValueError, skeletonize, im)          
         
     def test_skeletonize_all_foreground(self):
         im = np.ones((3,4))
-        result = skeletonize.skeletonize(im)
+        result = skeletonize(im)
     
     def test_skeletonize_single_point(self):
         im = np.zeros((5, 5), np.uint8)
         im[3, 3] = 1
-        result = skeletonize.skeletonize(im)
+        result = skeletonize(im)
         numpy.testing.assert_array_equal(result, im)
     
     def test_skeletonize_already_thinned(self):
@@ -48,7 +47,7 @@ class TestSkeletonize(unittest.TestCase):
         im[3,1:-1] = 1
         im[2, -1] = 1
         im[4, 0] = 1
-        result = skeletonize.skeletonize(im)
+        result = skeletonize(im)
         numpy.testing.assert_array_equal(result, im)
     
     def test_skeletonize_output(self):
@@ -56,7 +55,7 @@ class TestSkeletonize(unittest.TestCase):
         
         # make black the foreground
         im = (im==0)
-        result = skeletonize.skeletonize(im)
+        result = skeletonize(im)
         
         expected = np.load(os.path.join(data_dir, "bw_text_skeleton.npy"))
         numpy.testing.assert_array_equal(result, expected)
@@ -83,15 +82,14 @@ class TestSkeletonize(unittest.TestCase):
         circle2 = (ic - 135)**2 + (ir - 150)**2 < 20**2
         image[circle1] = 1
         image[circle2] = 0
-        result = skeletonize.skeletonize(image)
+        result = skeletonize(image)
         
-        # there should never be a 2x2 block of foreground pixels
-        #  in a skeleton
+        # 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')
-        self.assertFalse(numpy.any(blocks == 4))
+        assert not numpy.any(blocks == 4)
         
 
 if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+    np.testing.run_module_suite()