From 3167a07f1af96d11e63336e6294e30c1629ab99e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sun, 12 Aug 2012 19:37:02 +0200
Subject: [PATCH 01/16] add polygon approximation algorithm

---
 CONTRIBUTORS.txt                      |  7 +--
 skimage/measure/__init__.py           |  1 +
 skimage/measure/_polygon.py           | 88 +++++++++++++++++++++++++++
 skimage/measure/tests/test_polygon.py | 27 ++++++++
 4 files changed, 118 insertions(+), 5 deletions(-)
 create mode 100644 skimage/measure/_polygon.py
 create mode 100644 skimage/measure/tests/test_polygon.py

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index 531ab842..4abe9256 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -105,11 +105,8 @@
   Shape views: ``util.shape.view_as_windows`` and ``util.shape.view_as_blocks``
 
 - Johannes Schönberger
-  Polygon, circle and ellipse drawing functions
-  Adaptive thresholding
-  Implementation of Matlab's `regionprops`
-  Estimation of geometric transformation parameters
-  Local binary pattern texture classification
+  Drawing functions, adaptive thresholding, regionprops, geometric
+  transformations, LBPs, polygon approximations, and more.
 
 - Pavel Campr
   Fixes and tests for Histograms of Oriented Gradients.
diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py
index f5109698..d7397b2f 100755
--- a/skimage/measure/__init__.py
+++ b/skimage/measure/__init__.py
@@ -1,3 +1,4 @@
 from .find_contours import find_contours
 from ._regionprops import regionprops, perimeter
 from ._structural_similarity import structural_similarity
+from ._polygon import approximate_polygon
\ No newline at end of file
diff --git a/skimage/measure/_polygon.py b/skimage/measure/_polygon.py
new file mode 100644
index 00000000..0b80599a
--- /dev/null
+++ b/skimage/measure/_polygon.py
@@ -0,0 +1,88 @@
+import numpy as np
+
+
+def approximate_polygon(coords, tolerance):
+    """Approximate a polygonal chain with the specified tolerance.
+
+    It is based on the Douglas-Peucker algorithm.
+
+    Parameters
+    ----------
+    coords : (N, 2) array
+        Coordinate array.
+    tolerance : float
+        Maximum distance from original points of polygon to approximated
+        polygonal chain. If tolerance is 0, the original coordinate array
+        is returned.
+
+    Returns
+    -------
+    coords : (M, 2) array
+        Approximated polygonal chain where M <= N.
+
+    References
+    ----------
+    .. [1] http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
+    """
+    if tolerance == 0:
+        return coords
+
+    chain = np.zeros(coords.shape[0], 'bool')
+    # pre-allocate distance array for all points
+    dists = np.zeros(coords.shape[0])
+    chain[0] = True
+    chain[-1] = True
+    pos_stack = [(0, chain.shape[0] - 1)]
+
+    while 1:
+        start, end = pos_stack.pop()
+        # determine properties of current line segment
+        r0, c0 = coords[start, :]
+        r1, c1 = coords[end, :]
+        dr = r1 - r0
+        dc = c1 - c0
+        segment_angle = - np.arctan2(dr, dc)
+        segment_dist = c0 * np.sin(segment_angle) + r0 * np.cos(segment_angle)
+
+        # select points in-between line segment
+        segment_coords = coords[start + 1:end, :]
+        segment_dists = dists[start + 1:end]
+
+
+        # check whether to take perpendicular or euclidean distance with
+        # inner product of vectors
+
+        # vectors from points -> start and end
+        dr0 = segment_coords[:, 0] - r0
+        dc0 = segment_coords[:, 1] - c0
+        dr1 = segment_coords[:, 0] - r1
+        dc1 = segment_coords[:, 1] - c1
+        # vectors points -> start and end projected on start -> end vector
+        projected_lengths0 = dr0 * dr + dc0 * dc
+        projected_lengths1 = - dr1 * dr - dc1 * dc
+        perp = np.logical_and(projected_lengths0 > 0,
+                              projected_lengths1 > 0)
+        eucl = np.logical_not(perp)
+        segment_dists[perp] = np.abs(
+            segment_coords[perp, 0] * np.cos(segment_angle)
+            + segment_coords[perp, 1] * np.sin(segment_angle)
+            - segment_dist
+        )
+        segment_dists[eucl] = np.minimum(
+            # distance to start point
+            np.sqrt(dc0[eucl] ** 2 + dr0[eucl] ** 2),
+            # distance to end point
+            np.sqrt(dc1[eucl] ** 2 + dr1[eucl] ** 2)
+        )
+
+        if np.any(segment_dists > tolerance):
+            # select point with maximum distance to line
+            new_end = start + np.argmax(segment_dists) + 1
+            pos_stack.append((new_end, end))
+            pos_stack.append((start, new_end))
+            chain[new_end] = True
+
+        if len(pos_stack) == 0:
+            break
+
+    return coords[chain, :]
diff --git a/skimage/measure/tests/test_polygon.py b/skimage/measure/tests/test_polygon.py
new file mode 100644
index 00000000..e539f69c
--- /dev/null
+++ b/skimage/measure/tests/test_polygon.py
@@ -0,0 +1,27 @@
+import numpy as np
+from skimage.measure import approximate_polygon
+
+
+def test_approximate_polygon():
+    square = np.array([
+        [0, 0], [0, 1], [0, 2], [0, 3],
+        [1, 3], [2, 3], [3, 3],
+        [3, 2], [3, 1], [3, 0],
+        [2, 0], [1, 0], [0, 0]
+    ])
+
+    out = approximate_polygon(square, 0.1)
+    np.testing.assert_array_equal(out, square[(0, 3, 6, 9, 12), :])
+
+    out = approximate_polygon(square, 2.2)
+    np.testing.assert_array_equal(out, square[(0, 6, 12), :])
+
+    out = approximate_polygon(square[(0, 1, 3, 4, 5, 6, 7, 9, 11, 12), :], 0.1)
+    np.testing.assert_array_equal(out, square[(0, 3, 6, 9, 12), :])
+
+    out = approximate_polygon(square, -1)
+    np.testing.assert_array_equal(out, square)
+
+
+if __name__ == "__main__":
+    np.testing.run_module_suite()

From ec3ed52da2b0eb6704528d6de79b310370c94b27 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sun, 12 Aug 2012 21:10:57 +0200
Subject: [PATCH 02/16] add example script for polygon approximation

---
 doc/examples/plot_approximate_polygon.py | 35 ++++++++++++++++++++++++
 1 file changed, 35 insertions(+)
 create mode 100644 doc/examples/plot_approximate_polygon.py

diff --git a/doc/examples/plot_approximate_polygon.py b/doc/examples/plot_approximate_polygon.py
new file mode 100644
index 00000000..8c1a126d
--- /dev/null
+++ b/doc/examples/plot_approximate_polygon.py
@@ -0,0 +1,35 @@
+"""
+====================
+Approximate Polygons
+====================
+
+This example shows how to approximate polygonal chains with the Douglas-Peucker
+algorithm.
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from skimage.draw import ellipse
+from skimage.measure import find_contours, approximate_polygon
+from skimage.morphology import erosion
+
+img = np.zeros((800, 800), 'int32')
+
+rr, cc = ellipse(250, 250, 180, 230, img.shape)
+img[rr, cc] = 1
+rr, cc = ellipse(600, 600, 150, 90, img.shape)
+img[rr, cc] = 1
+
+plt.gray()
+plt.imshow(img)
+
+strel = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], 'uint8')
+for contour in find_contours(img, 0):
+    coords = approximate_polygon(contour, 4)
+    plt.plot(coords[:, 1], coords[:, 0], '-r', linewidth=2)
+    coords = approximate_polygon(contour, 40)
+    plt.plot(coords[:, 1], coords[:, 0], '-g', linewidth=2)
+
+plt.axis((0, 800, 0, 800))
+plt.axis('off')
+plt.show()

From 77e4c37e5c318d39feb190dca1a423908c3911a7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sun, 12 Aug 2012 21:12:49 +0200
Subject: [PATCH 03/16] remove unused import and code from polygon
 approximation example

---
 doc/examples/plot_approximate_polygon.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/doc/examples/plot_approximate_polygon.py b/doc/examples/plot_approximate_polygon.py
index 8c1a126d..cec98492 100644
--- a/doc/examples/plot_approximate_polygon.py
+++ b/doc/examples/plot_approximate_polygon.py
@@ -11,7 +11,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 from skimage.draw import ellipse
 from skimage.measure import find_contours, approximate_polygon
-from skimage.morphology import erosion
+
 
 img = np.zeros((800, 800), 'int32')
 
@@ -23,7 +23,6 @@ img[rr, cc] = 1
 plt.gray()
 plt.imshow(img)
 
-strel = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]], 'uint8')
 for contour in find_contours(img, 0):
     coords = approximate_polygon(contour, 4)
     plt.plot(coords[:, 1], coords[:, 0], '-r', linewidth=2)

From 74b358af824232e6419cb40e459aeeaa4805c900 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Sun, 12 Aug 2012 21:14:27 +0200
Subject: [PATCH 04/16] remove colons in comments

---
 doc/examples/plot_shapes.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/doc/examples/plot_shapes.py b/doc/examples/plot_shapes.py
index 8107fc80..9e586209 100644
--- a/doc/examples/plot_shapes.py
+++ b/doc/examples/plot_shapes.py
@@ -19,11 +19,11 @@ import numpy as np
 
 img = np.zeros((500, 500, 3), 'uint8')
 
-#: draw line
+# draw line
 rr, cc = line(120, 123, 20, 400)
 img[rr,cc,0] = 255
 
-#: fill polygon
+# fill polygon
 poly = np.array((
     (300, 300),
     (480, 320),
@@ -34,11 +34,11 @@ poly = np.array((
 rr, cc = polygon(poly[:,0], poly[:,1], img.shape)
 img[rr,cc,1] = 255
 
-#: fill circle
+# fill circle
 rr, cc = circle(200, 200, 100, img.shape)
 img[rr,cc,:] = (255, 255, 0)
 
-#: fill ellipse
+# fill ellipse
 rr, cc = ellipse(300, 300, 100, 200, img.shape)
 img[rr,cc,2] = 255
 

From add94d24cecda037a6b64d3989dde8142731b995 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Mon, 13 Aug 2012 09:30:28 +0200
Subject: [PATCH 05/16] make loop condition more readable

---
 skimage/measure/_polygon.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/skimage/measure/_polygon.py b/skimage/measure/_polygon.py
index 0b80599a..9138faac 100644
--- a/skimage/measure/_polygon.py
+++ b/skimage/measure/_polygon.py
@@ -33,8 +33,9 @@ def approximate_polygon(coords, tolerance):
     chain[0] = True
     chain[-1] = True
     pos_stack = [(0, chain.shape[0] - 1)]
+    end_of_chain = False
 
-    while 1:
+    while not end_of_chain:
         start, end = pos_stack.pop()
         # determine properties of current line segment
         r0, c0 = coords[start, :]
@@ -83,6 +84,6 @@ def approximate_polygon(coords, tolerance):
             chain[new_end] = True
 
         if len(pos_stack) == 0:
-            break
+            end_of_chain = True
 
     return coords[chain, :]

From dcbea756775f9e531861bd0df4a49480012277b3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Mon, 13 Aug 2012 22:21:40 +0200
Subject: [PATCH 06/16] improved polygon approximation example script

made existing example parameters more readable and add another example with
a more complex shape based on a subdivision algorithm.
---
 doc/examples/plot_approximate_polygon.py | 60 +++++++++++++++++++++---
 1 file changed, 53 insertions(+), 7 deletions(-)

diff --git a/doc/examples/plot_approximate_polygon.py b/doc/examples/plot_approximate_polygon.py
index cec98492..ec7e4f4d 100644
--- a/doc/examples/plot_approximate_polygon.py
+++ b/doc/examples/plot_approximate_polygon.py
@@ -9,6 +9,7 @@ algorithm.
 
 import numpy as np
 import matplotlib.pyplot as plt
+from scipy import signal
 from skimage.draw import ellipse
 from skimage.measure import find_contours, approximate_polygon
 
@@ -20,15 +21,60 @@ img[rr, cc] = 1
 rr, cc = ellipse(600, 600, 150, 90, img.shape)
 img[rr, cc] = 1
 
+hand = np.array([[1.64516129, 1.16145833],
+                 [1.64516129, 1.59375   ],
+                 [1.35080645, 1.921875  ],
+                 [1.375     , 2.18229167],
+                 [1.68548387, 1.9375    ],
+                 [1.60887097, 2.55208333],
+                 [1.68548387, 2.69791667],
+                 [1.76209677, 2.56770833],
+                 [1.83064516, 1.97395833],
+                 [1.89516129, 2.75      ],
+                 [1.9516129 , 2.84895833],
+                 [2.01209677, 2.76041667],
+                 [1.99193548, 1.99479167],
+                 [2.11290323, 2.63020833],
+                 [2.2016129 , 2.734375  ],
+                 [2.25403226, 2.60416667],
+                 [2.14919355, 1.953125  ],
+                 [2.30645161, 2.36979167],
+                 [2.39112903, 2.36979167],
+                 [2.41532258, 2.1875    ],
+                 [2.1733871 , 1.703125  ],
+                 [2.07782258, 1.16666667]])
+
+
+def doo_sabin_subdivision(points):
+    new_points = np.zeros((2 * (len(points) - 1), 2))
+    new_points[::2] = 3 / 4. * points[:-1] + 1 / 4. * points[1:]
+    new_points[1::2] = 1 / 4. * points[:-1] + 3 / 4. * points[1:]
+    return new_points
+
+new_hand = hand.copy()
+for _ in range(5):
+    new_hand = doo_sabin_subdivision(new_hand)
+
+appr_hand = approximate_polygon(new_hand, tolerance=0.02)
+
+print "Number of coordinates:", len(hand), len(new_hand), len(appr_hand)
+
+fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(9, 4))
+
+ax1.plot(hand[:, 0], hand[:, 1])
+ax1.plot(new_hand[:, 0], new_hand[:, 1])
+ax1.plot(appr_hand[:, 0], appr_hand[:, 1])
+
 plt.gray()
-plt.imshow(img)
+ax2.imshow(img)
 
 for contour in find_contours(img, 0):
-    coords = approximate_polygon(contour, 4)
-    plt.plot(coords[:, 1], coords[:, 0], '-r', linewidth=2)
-    coords = approximate_polygon(contour, 40)
-    plt.plot(coords[:, 1], coords[:, 0], '-g', linewidth=2)
+    coords = approximate_polygon(contour, tolerance=2.5)
+    ax2.plot(coords[:, 1], coords[:, 0], '-r', linewidth=2)
+    coords2 = approximate_polygon(contour, tolerance=39.5)
+    ax2.plot(coords2[:, 1], coords2[:, 0], '-g', linewidth=2)
+    print "Number of coordinates:", len(contour), len(coords), len(coords2)
+
+ax2.axis((0, 800, 0, 800))
 
-plt.axis((0, 800, 0, 800))
-plt.axis('off')
 plt.show()

From 40c58f33334698d92bcc42e42225ed8c8e6034c5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 15 Aug 2012 08:48:40 +0200
Subject: [PATCH 07/16] add subdivision of polygonal curves using B-Splines

---
 skimage/measure/__init__.py |  2 +-
 skimage/measure/_polygon.py | 64 +++++++++++++++++++++++++++++++++++++
 2 files changed, 65 insertions(+), 1 deletion(-)

diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py
index d7397b2f..1c92d9ee 100755
--- a/skimage/measure/__init__.py
+++ b/skimage/measure/__init__.py
@@ -1,4 +1,4 @@
 from .find_contours import find_contours
 from ._regionprops import regionprops, perimeter
 from ._structural_similarity import structural_similarity
-from ._polygon import approximate_polygon
\ No newline at end of file
+from ._polygon import approximate_polygon, subdivide_polygon
\ No newline at end of file
diff --git a/skimage/measure/_polygon.py b/skimage/measure/_polygon.py
index 9138faac..01bbdd98 100644
--- a/skimage/measure/_polygon.py
+++ b/skimage/measure/_polygon.py
@@ -1,4 +1,5 @@
 import numpy as np
+from scipy import signal
 
 
 def approximate_polygon(coords, tolerance):
@@ -87,3 +88,66 @@ def approximate_polygon(coords, tolerance):
             end_of_chain = True
 
     return coords[chain, :]
+
+
+SUBDIVISION_DEGREES = {
+    # degree: (mask_even, mask_odd)
+    #         extracted from (degree + 2)th row of Pascal's triangle
+    1: ([1, 1], [1, 1]),
+    2: ([3, 1], [1, 3]),
+    3: ([1, 6, 1], [0, 4, 4]),
+    4: ([5, 10, 1], [1, 10, 5]),
+    5: ([1, 15, 15, 1], [0, 6, 20, 6]),
+    6: ([7, 35, 21, 1], [1, 21, 35, 7]),
+    7: ([1, 28, 70, 28, 1], [0, 8, 56, 56, 8]),
+}
+
+
+def subdivide_polygon(coords, degree=2):
+    """Subdivision of polygonal curves using B-Splines.
+
+    Parameters
+    ----------
+    coords : (N, 2) array
+        Coordinate array.
+    method : {1, 2, 3, 4, 5, 6, 7}
+        Degree of B-Spline. Default is 2.
+
+    Returns
+    -------
+    coords : (M, 2) array
+        Subdivided coordinate array.
+
+    References
+    ----------
+    .. [1] http://mrl.nyu.edu/publications/subdiv-course2000/coursenotes00.pdf
+    """
+    if degree not in SUBDIVISION_DEGREES:
+        raise ValueError("Invalid B-Spline degree. Only degree 1 - 7 is "
+                         "supported.")
+
+    circular = np.all(coords[0, :] == coords[-1, :])
+
+    method = 'valid'
+    if circular:
+        coords = coords[:-1, :]
+        method = 'same'
+
+    mask_even, mask_odd = SUBDIVISION_DEGREES[degree]
+    mask_even = np.array(mask_even, 'float') / 2 ** degree
+    mask_odd = np.array(mask_odd, 'float') / 2 ** degree
+
+    even = signal.convolve2d(coords.T, np.atleast_2d(mask_even), mode=method,
+                             boundary='wrap')
+    odd = signal.convolve2d(coords.T, np.atleast_2d(mask_odd), mode=method,
+                            boundary='wrap')
+
+    out = np.zeros((even.shape[1] + odd.shape[1], 2))
+    out[1::2] = even.T
+    out[::2] = odd.T
+
+    if circular:
+        # close polygon
+        out = np.vstack([out, out[0, :]])
+
+    return out

From 339c33eac3e47480d003e7f9b83672fadaeda5a8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 15 Aug 2012 08:50:31 +0200
Subject: [PATCH 08/16] add note about size if resulting polygonal curves

---
 skimage/measure/_polygon.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/skimage/measure/_polygon.py b/skimage/measure/_polygon.py
index 01bbdd98..2d938af7 100644
--- a/skimage/measure/_polygon.py
+++ b/skimage/measure/_polygon.py
@@ -7,6 +7,9 @@ def approximate_polygon(coords, tolerance):
 
     It is based on the Douglas-Peucker algorithm.
 
+    Note that the approximated polygon is always within the convex hull of the
+    original polygon.
+
     Parameters
     ----------
     coords : (N, 2) array
@@ -106,6 +109,9 @@ SUBDIVISION_DEGREES = {
 def subdivide_polygon(coords, degree=2):
     """Subdivision of polygonal curves using B-Splines.
 
+    Note that the resulting curve is always within the convex hull of the
+    original polygon.
+
     Parameters
     ----------
     coords : (N, 2) array

From 76c3574755cf92f59c1cd140caa746d7ec112e0d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 15 Aug 2012 08:55:43 +0200
Subject: [PATCH 09/16] update polygon example with subdivision

---
 ...approximate_polygon.py => plot_polygon.py} | 31 +++++++++----------
 1 file changed, 14 insertions(+), 17 deletions(-)
 rename doc/examples/{plot_approximate_polygon.py => plot_polygon.py} (84%)

diff --git a/doc/examples/plot_approximate_polygon.py b/doc/examples/plot_polygon.py
similarity index 84%
rename from doc/examples/plot_approximate_polygon.py
rename to doc/examples/plot_polygon.py
index ec7e4f4d..1240ff39 100644
--- a/doc/examples/plot_approximate_polygon.py
+++ b/doc/examples/plot_polygon.py
@@ -9,18 +9,11 @@ algorithm.
 
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy import signal
 from skimage.draw import ellipse
-from skimage.measure import find_contours, approximate_polygon
+from skimage.measure import find_contours, approximate_polygon, \
+    subdivide_polygon
 
 
-img = np.zeros((800, 800), 'int32')
-
-rr, cc = ellipse(250, 250, 180, 230, img.shape)
-img[rr, cc] = 1
-rr, cc = ellipse(600, 600, 150, 90, img.shape)
-img[rr, cc] = 1
-
 hand = np.array([[1.64516129, 1.16145833],
                  [1.64516129, 1.59375   ],
                  [1.35080645, 1.921875  ],
@@ -44,17 +37,12 @@ hand = np.array([[1.64516129, 1.16145833],
                  [2.1733871 , 1.703125  ],
                  [2.07782258, 1.16666667]])
 
-
-def doo_sabin_subdivision(points):
-    new_points = np.zeros((2 * (len(points) - 1), 2))
-    new_points[::2] = 3 / 4. * points[:-1] + 1 / 4. * points[1:]
-    new_points[1::2] = 1 / 4. * points[:-1] + 3 / 4. * points[1:]
-    return new_points
-
+# subdivide polygon using 2nd degree B-Splines
 new_hand = hand.copy()
 for _ in range(5):
-    new_hand = doo_sabin_subdivision(new_hand)
+    new_hand = subdivide_polygon(new_hand, degree=2)
 
+# approximate subdivided polygon with Douglas-Peucker algorithm
 appr_hand = approximate_polygon(new_hand, tolerance=0.02)
 
 print "Number of coordinates:", len(hand), len(new_hand), len(appr_hand)
@@ -65,9 +53,18 @@ ax1.plot(hand[:, 0], hand[:, 1])
 ax1.plot(new_hand[:, 0], new_hand[:, 1])
 ax1.plot(appr_hand[:, 0], appr_hand[:, 1])
 
+
+# create two ellipses in image
+img = np.zeros((800, 800), 'int32')
+rr, cc = ellipse(250, 250, 180, 230, img.shape)
+img[rr, cc] = 1
+rr, cc = ellipse(600, 600, 150, 90, img.shape)
+img[rr, cc] = 1
+
 plt.gray()
 ax2.imshow(img)
 
+# approximate / simplify coordinates of the two ellipses
 for contour in find_contours(img, 0):
     coords = approximate_polygon(contour, tolerance=2.5)
     ax2.plot(coords[:, 1], coords[:, 0], '-r', linewidth=2)

From f267f656664ea4040cfe78de7c2da90487ef10c0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 15 Aug 2012 09:04:06 +0200
Subject: [PATCH 10/16] add test case for subdivision of polygons

---
 skimage/measure/tests/test_polygon.py | 24 ++++++++++++++++--------
 1 file changed, 16 insertions(+), 8 deletions(-)

diff --git a/skimage/measure/tests/test_polygon.py b/skimage/measure/tests/test_polygon.py
index e539f69c..4644fa84 100644
--- a/skimage/measure/tests/test_polygon.py
+++ b/skimage/measure/tests/test_polygon.py
@@ -1,15 +1,16 @@
 import numpy as np
-from skimage.measure import approximate_polygon
+from skimage.measure import approximate_polygon, subdivide_polygon
+
+
+square = np.array([
+    [0, 0], [0, 1], [0, 2], [0, 3],
+    [1, 3], [2, 3], [3, 3],
+    [3, 2], [3, 1], [3, 0],
+    [2, 0], [1, 0], [0, 0]
+])
 
 
 def test_approximate_polygon():
-    square = np.array([
-        [0, 0], [0, 1], [0, 2], [0, 3],
-        [1, 3], [2, 3], [3, 3],
-        [3, 2], [3, 1], [3, 0],
-        [2, 0], [1, 0], [0, 0]
-    ])
-
     out = approximate_polygon(square, 0.1)
     np.testing.assert_array_equal(out, square[(0, 3, 6, 9, 12), :])
 
@@ -23,5 +24,12 @@ def test_approximate_polygon():
     np.testing.assert_array_equal(out, square)
 
 
+def test_subdivide_polygon():
+    for degree in range(1, 7):
+        out = subdivide_polygon(square, degree)
+        np.testing.assert_array_equal(out[-1], out[0])
+        np.testing.assert_equal(out.shape[0], 2 * square.shape[0] - 1)
+
+
 if __name__ == "__main__":
     np.testing.run_module_suite()

From 02bf7130172655a4beaf3129fd3c2715eb37d77f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 15 Aug 2012 09:08:20 +0200
Subject: [PATCH 11/16] add some more comments to subdivision of polygons code

---
 skimage/measure/_polygon.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/skimage/measure/_polygon.py b/skimage/measure/_polygon.py
index 2d938af7..c1eaa9d5 100644
--- a/skimage/measure/_polygon.py
+++ b/skimage/measure/_polygon.py
@@ -136,10 +136,13 @@ def subdivide_polygon(coords, degree=2):
 
     method = 'valid'
     if circular:
+        # remove last coordinate because of wrapping
         coords = coords[:-1, :]
+        # circular convolution by wrapping boundaries
         method = 'same'
 
     mask_even, mask_odd = SUBDIVISION_DEGREES[degree]
+    # divide by total weight
     mask_even = np.array(mask_even, 'float') / 2 ** degree
     mask_odd = np.array(mask_odd, 'float') / 2 ** degree
 

From 22a89076c26d4539f7b2947afb0fe23eb08c6bfb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 15 Aug 2012 09:41:19 +0200
Subject: [PATCH 12/16] make subdivision mask code more readable und debuggable

---
 skimage/measure/_polygon.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/skimage/measure/_polygon.py b/skimage/measure/_polygon.py
index c1eaa9d5..06e7a29d 100644
--- a/skimage/measure/_polygon.py
+++ b/skimage/measure/_polygon.py
@@ -93,7 +93,8 @@ def approximate_polygon(coords, tolerance):
     return coords[chain, :]
 
 
-SUBDIVISION_DEGREES = {
+# B-Spline subdivision
+_SUBDIVISION_MASKS = {
     # degree: (mask_even, mask_odd)
     #         extracted from (degree + 2)th row of Pascal's triangle
     1: ([1, 1], [1, 1]),
@@ -128,7 +129,7 @@ def subdivide_polygon(coords, degree=2):
     ----------
     .. [1] http://mrl.nyu.edu/publications/subdiv-course2000/coursenotes00.pdf
     """
-    if degree not in SUBDIVISION_DEGREES:
+    if degree not in _SUBDIVISION_MASKS:
         raise ValueError("Invalid B-Spline degree. Only degree 1 - 7 is "
                          "supported.")
 
@@ -141,10 +142,10 @@ def subdivide_polygon(coords, degree=2):
         # circular convolution by wrapping boundaries
         method = 'same'
 
-    mask_even, mask_odd = SUBDIVISION_DEGREES[degree]
+    mask_even, mask_odd = _SUBDIVISION_MASKS[degree]
     # divide by total weight
-    mask_even = np.array(mask_even, 'float') / 2 ** degree
-    mask_odd = np.array(mask_odd, 'float') / 2 ** degree
+    mask_even = np.array(mask_even, np.float) / (2 ** degree)
+    mask_odd = np.array(mask_odd, np.float) / (2 ** degree)
 
     even = signal.convolve2d(coords.T, np.atleast_2d(mask_even), mode=method,
                              boundary='wrap')

From 50df20f3b71163f62423c812f8415aef49a426cb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 15 Aug 2012 09:41:49 +0200
Subject: [PATCH 13/16] fix wrong parameter name in doc string

---
 skimage/measure/_polygon.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/measure/_polygon.py b/skimage/measure/_polygon.py
index 06e7a29d..07af5f58 100644
--- a/skimage/measure/_polygon.py
+++ b/skimage/measure/_polygon.py
@@ -117,7 +117,7 @@ def subdivide_polygon(coords, degree=2):
     ----------
     coords : (N, 2) array
         Coordinate array.
-    method : {1, 2, 3, 4, 5, 6, 7}
+    degree : {1, 2, 3, 4, 5, 6, 7}
         Degree of B-Spline. Default is 2.
 
     Returns

From 26f6bd4fbafe2a73a5c928a78cea18b25ba6bbed Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Wed, 15 Aug 2012 09:47:34 +0200
Subject: [PATCH 14/16] add test case for non-circular subdivision of polygons

---
 skimage/measure/tests/test_polygon.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/skimage/measure/tests/test_polygon.py b/skimage/measure/tests/test_polygon.py
index 4644fa84..f361bbdc 100644
--- a/skimage/measure/tests/test_polygon.py
+++ b/skimage/measure/tests/test_polygon.py
@@ -1,6 +1,6 @@
 import numpy as np
 from skimage.measure import approximate_polygon, subdivide_polygon
-
+from skimage.measure._polygon import _SUBDIVISION_MASKS
 
 square = np.array([
     [0, 0], [0, 1], [0, 2], [0, 3],
@@ -26,9 +26,14 @@ def test_approximate_polygon():
 
 def test_subdivide_polygon():
     for degree in range(1, 7):
+        # test circular
         out = subdivide_polygon(square, degree)
         np.testing.assert_array_equal(out[-1], out[0])
         np.testing.assert_equal(out.shape[0], 2 * square.shape[0] - 1)
+        # test non-circular
+        out = subdivide_polygon(square[:-1], degree)
+        mask_len = len(_SUBDIVISION_MASKS[degree][0])
+        np.testing.assert_equal(out.shape[0], 2 * (square.shape[0] - mask_len))
 
 
 if __name__ == "__main__":

From 8ff263979be76adc3d5e5288404bc1cb2b107eba Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Thu, 16 Aug 2012 23:03:16 +0200
Subject: [PATCH 15/16] Add possibility to preserve ends and update tests

---
 skimage/measure/_polygon.py           | 12 +++++++---
 skimage/measure/tests/test_polygon.py | 34 ++++++++++++++++++++-------
 2 files changed, 34 insertions(+), 12 deletions(-)

diff --git a/skimage/measure/_polygon.py b/skimage/measure/_polygon.py
index 07af5f58..633a9418 100644
--- a/skimage/measure/_polygon.py
+++ b/skimage/measure/_polygon.py
@@ -107,18 +107,21 @@ _SUBDIVISION_MASKS = {
 }
 
 
-def subdivide_polygon(coords, degree=2):
+def subdivide_polygon(coords, degree=2, preserve_ends=False):
     """Subdivision of polygonal curves using B-Splines.
 
     Note that the resulting curve is always within the convex hull of the
-    original polygon.
+    original polygon. Circular polygons stay closed after subdivision.
 
     Parameters
     ----------
     coords : (N, 2) array
         Coordinate array.
-    degree : {1, 2, 3, 4, 5, 6, 7}
+    degree : {1, 2, 3, 4, 5, 6, 7}, optional
         Degree of B-Spline. Default is 2.
+    preserve_ends : bool, optional
+        Preserve first and last coordinate of non-circular polygon. Default is
+        False.
 
     Returns
     -------
@@ -160,4 +163,7 @@ def subdivide_polygon(coords, degree=2):
         # close polygon
         out = np.vstack([out, out[0, :]])
 
+    if preserve_ends and not circular:
+        out = np.vstack([coords[0, :], out, coords[-1, :]])
+
     return out
diff --git a/skimage/measure/tests/test_polygon.py b/skimage/measure/tests/test_polygon.py
index f361bbdc..f239c735 100644
--- a/skimage/measure/tests/test_polygon.py
+++ b/skimage/measure/tests/test_polygon.py
@@ -25,15 +25,31 @@ def test_approximate_polygon():
 
 
 def test_subdivide_polygon():
-    for degree in range(1, 7):
-        # test circular
-        out = subdivide_polygon(square, degree)
-        np.testing.assert_array_equal(out[-1], out[0])
-        np.testing.assert_equal(out.shape[0], 2 * square.shape[0] - 1)
-        # test non-circular
-        out = subdivide_polygon(square[:-1], degree)
-        mask_len = len(_SUBDIVISION_MASKS[degree][0])
-        np.testing.assert_equal(out.shape[0], 2 * (square.shape[0] - mask_len))
+    new_square1 = square
+    new_square2 = square[:-1]
+    new_square3 = square[:-1]
+    # test iterative subdvision
+    for _ in range(10):
+        square1, square2, square3 = new_square1, new_square2, new_square3
+        # test different B-Spline degrees
+        for degree in range(1, 7):
+            mask_len = len(_SUBDIVISION_MASKS[degree][0])
+            # test circular
+            new_square1 = subdivide_polygon(square1, degree)
+            np.testing.assert_array_equal(new_square1[-1], new_square1[0])
+            np.testing.assert_equal(new_square1.shape[0],
+                                    2 * square1.shape[0] - 1)
+            # test non-circular
+            new_square2 = subdivide_polygon(square2, degree)
+            np.testing.assert_equal(new_square2.shape[0],
+                                    2 * (square2.shape[0] - mask_len + 1))
+            # test non-circular, preserve_ends
+            new_square3 = subdivide_polygon(square3, degree, True)
+            np.testing.assert_equal(new_square3[0], square3[0])
+            np.testing.assert_equal(new_square3[-1], square3[-1])
+            print mask_len
+            np.testing.assert_equal(new_square3.shape[0],
+                                    2 * (square3.shape[0] - mask_len + 2))
 
 
 if __name__ == "__main__":

From f9c81e9ed39cad4616a0c809f1424a7e92e51c31 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <hannesschoenberger@gmail.com>
Date: Thu, 16 Aug 2012 23:06:40 +0200
Subject: [PATCH 16/16] Update subdivision example with preserve_ends parameter

---
 doc/examples/plot_polygon.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/examples/plot_polygon.py b/doc/examples/plot_polygon.py
index 1240ff39..6eaef0bf 100644
--- a/doc/examples/plot_polygon.py
+++ b/doc/examples/plot_polygon.py
@@ -40,7 +40,7 @@ hand = np.array([[1.64516129, 1.16145833],
 # subdivide polygon using 2nd degree B-Splines
 new_hand = hand.copy()
 for _ in range(5):
-    new_hand = subdivide_polygon(new_hand, degree=2)
+    new_hand = subdivide_polygon(new_hand, degree=2, preserve_ends=True)
 
 # approximate subdivided polygon with Douglas-Peucker algorithm
 appr_hand = approximate_polygon(new_hand, tolerance=0.02)