diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 9f3f2aa7..d582ab63 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -129,3 +129,6 @@ - Anders Boesen Lindbo Larsen Dense DAISY feature description, circle perimeter drawing. + +- François Boulogne + Andres Method for circle perimeter, ellipse perimeter drawing. 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 b3b24655..033de8f0 100644 --- a/skimage/draw/_draw.pyx +++ b/skimage/draw/_draw.pyx @@ -126,15 +126,15 @@ def polygon(y, x, shape=None): return np.array(rr), np.array(cc) -def ellipse(double cy, double cx, double b, double a, shape=None): +def ellipse(double cy, double cx, double yradius, double xradius, shape=None): """Generate coordinates of pixels within ellipse. Parameters ---------- cy, cx : double Centre coordinate of ellipse. - b, a: double - Minor and major semi-axes. ``(x/a)**2 + (y/b)**2 = 1``. + yradius, xradius : double + Minor and major semi-axes. ``(x/xradius)**2 + (y/yradius)**2 = 1``. Returns ------- @@ -143,10 +143,10 @@ def ellipse(double cy, double cx, double b, double a, shape=None): May be used to directly index into an array, e.g. ``img[rr, cc] = 1``. """ - cdef int minr = max(0, cy - b) - cdef int maxr = math.ceil(cy + b) - cdef int minc = max(0, cx - a) - cdef int maxc = math.ceil(cx + a) + cdef int minr = max(0, cy - yradius) + cdef int maxr = math.ceil(cy + yradius) + cdef int minc = max(0, cx - xradius) + cdef int maxc = math.ceil(cx + xradius) # make sure output coordinates do not exceed image size if shape is not None: @@ -159,9 +159,9 @@ def ellipse(double cy, double cx, double b, double a, shape=None): cdef list rr = list() cdef list cc = list() - for r in range(minr, maxr + 1): - for c in range(minc, maxc + 1): - if sqrt(((r - cy) / b)**2 + ((c - cx) / a)**2) < 1: + for r in range(minr, maxr+1): + for c in range(minc, maxc+1): + if sqrt(((r - cy) / yradius)**2 + ((c - cx) / xradius)**2) < 1: rr.append(r) cc.append(c) @@ -266,6 +266,92 @@ 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 + + def set_color(img, coords, color): """Set pixel color in the image at the given coordinates. Coordinates that exceed the shape of the image will be ignored. diff --git a/skimage/draw/tests/test_draw.py b/skimage/draw/tests/test_draw.py index b2325015..ef09747c 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(): @@ -155,6 +155,7 @@ def test_circle_perimeter_bresenham(): rr, cc = circle_perimeter(7, 7, 0, method='bresenham') img[rr, cc] = 1 assert(np.sum(img) == 1) + assert(img[7][7] == 1) img = np.zeros((17, 15), 'uint8') rr, cc = circle_perimeter(7, 7, 7, method='bresenham') @@ -185,6 +186,7 @@ def test_circle_perimeter_andres(): rr, cc = circle_perimeter(7, 7, 0, method='andres') img[rr, cc] = 1 assert(np.sum(img) == 1) + assert(img[7][7] == 1) img = np.zeros((17, 15), 'uint8') rr, cc = circle_perimeter(7, 7, 7, method='andres') @@ -210,7 +212,6 @@ def test_circle_perimeter_andres(): ) assert_array_equal(img, img_) - def test_ellipse(): img = np.zeros((15, 15), 'uint8') @@ -237,6 +238,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) + assert(img[15][7] == 1) + + img = np.zeros((30, 15), 'uint8') + rr, cc = ellipse_perimeter(15, 7, 14, 6) + img[rr, cc] = 1 + 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