From 075c43109abeec1dfd35e1a7af7365143a665649 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Fran=C3=A7ois=20Boulogne?= <fboulogne sciunto org>
Date: Sun, 13 Jan 2013 16:39:02 +0100
Subject: [PATCH] new method: ellipse_perimeter

---
 skimage/draw/__init__.py        |  2 +-
 skimage/draw/_draw.pyx          | 84 +++++++++++++++++++++++++++++++++
 skimage/draw/tests/test_draw.py | 47 +++++++++++++++++-
 3 files changed, 130 insertions(+), 3 deletions(-)

diff --git a/skimage/draw/__init__.py b/skimage/draw/__init__.py
index 4b30ea18..67692f1b 100644
--- a/skimage/draw/__init__.py
+++ b/skimage/draw/__init__.py
@@ -1,2 +1,2 @@
-from ._draw import line, polygon, ellipse, circle, circle_perimeter, set_color
+from ._draw import line, polygon, ellipse, circle, circle_perimeter, ellipse_perimeter, set_color
 bresenham = line
diff --git a/skimage/draw/_draw.pyx b/skimage/draw/_draw.pyx
index 66b41d2c..73c9a341 100644
--- a/skimage/draw/_draw.pyx
+++ b/skimage/draw/_draw.pyx
@@ -263,6 +263,90 @@ def circle_perimeter(int cy, int cx, int radius, method='bresenham'):
 
     return np.array(rr) + cy, np.array(cc) + cx
 
+def ellipse_perimeter(int cy, int cx, int yradius, int xradius):
+    """Generate ellipse perimeter coordinates.
+
+    Parameters
+    ----------
+    cy, cx : int
+        Centre coordinate of ellipse.
+    yradius, xradius: int
+        Main radial values.
+
+    Returns
+    -------
+    rr, cc : (N,) ndarray of int
+        Indices of pixels that belong to the circle perimeter.
+        May be used to directly index into an array, e.g.
+        ``img[rr, cc] = 1``.
+
+    References
+    ----------
+    .. [1] J. Kennedy "A fast Bresenham type algorithm for
+        drawing ellipses".
+    """
+    # If both radii == 0, return the center
+    # to avoid infinite loop in 2nd set
+    if (xradius == 0 and yradius == 0):
+        return np.array(cy), np.array(cx)
+
+    # a and b are xradius an yradius
+    # compute 2a^2 and 2b^2
+    cdef int twoasquared = 2 * xradius * xradius
+    cdef int twobsquared = 2 * yradius * yradius
+
+    # Pixels
+    cdef list px = list()
+    cdef list py = list()
+
+    # First set of points:
+    # start at the top
+    cdef int x = xradius
+    cdef int y = 0
+
+    cdef int err = 0
+    cdef int xstop = twobsquared * xradius
+    cdef int ystop = 0
+    cdef int xchange = yradius * yradius * (1 - 2 * xradius)
+    cdef int ychange = xradius * xradius
+
+    while(xstop > ystop):
+        px.extend([x, -x, -x, x])
+        py.extend([y, y, -y, -y])
+        y += 1
+        ystop += twoasquared
+        err += ychange
+        ychange += twoasquared
+        if ((2 * err + xchange) > 0):
+            x -= 1
+            xstop -= twobsquared
+            err += xchange
+            xchange += twobsquared
+
+    # Second set of points:
+    x = 0
+    y = yradius
+
+    err = 0
+    xstop = 0
+    ystop = twoasquared * yradius
+    xchange = yradius * yradius
+    ychange = xradius * xradius * (1 - 2 * yradius)
+
+    while(xstop <= ystop):
+        px.extend([x, -x, -x, x])
+        py.extend([y, y, -y, -y])
+        x += 1
+        xstop += twobsquared
+        err += xchange
+        xchange += twobsquared
+        if ((2 * err + ychange) > 0):
+            y -= 1
+            ystop -= twoasquared
+            err += ychange
+            ychange += twobsquared
+
+    return np.array(py) + cy, np.array(px) + cx
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
diff --git a/skimage/draw/tests/test_draw.py b/skimage/draw/tests/test_draw.py
index b2325015..7b4edbe4 100644
--- a/skimage/draw/tests/test_draw.py
+++ b/skimage/draw/tests/test_draw.py
@@ -1,7 +1,7 @@
 from numpy.testing import assert_array_equal
 import numpy as np
 
-from skimage.draw import line, polygon, circle, circle_perimeter, ellipse
+from skimage.draw import line, polygon, circle, circle_perimeter, ellipse, ellipse_perimeter
 
 
 def test_line_horizontal():
@@ -210,7 +210,6 @@ def test_circle_perimeter_andres():
     )
     assert_array_equal(img, img_)
 
-
 def test_ellipse():
     img = np.zeros((15, 15), 'uint8')
 
@@ -237,6 +236,50 @@ def test_ellipse():
 
     assert_array_equal(img, img_)
 
+def test_ellipse_perimeter():
+    img = np.zeros((30, 15), 'uint8')
+    rr, cc = ellipse_perimeter(15, 7, 0, 0)
+    img[rr, cc] = 1
+    assert(np.sum(img) == 1)
+
+    img = np.zeros((30, 15), 'uint8')
+    rr, cc = ellipse_perimeter(15, 7, 14, 6)
+    img[rr, cc] = 1
+    print(img)
+    img_ = np.array(
+      [[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+       [0, 0, 0, 0, 0, 0, 1, 1, 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, 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, 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, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
+       [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 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, 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, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+       [0, 0, 0, 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, 1, 1, 0, 0, 0, 0, 0, 0]]
+    )
+
+    assert_array_equal(img, img_)
 
 if __name__ == "__main__":
     from numpy.testing import run_module_suite