import Fedor's contribution

CONTRIBUTORS.txt

@@ -154,3 +154,6 @@
 
 - Riaan van den Dool
   skimage.io plugin: GDAL
+
+- Fedor Morozov
+  Drawing: Wu's anti-aliased circle

doc/examples/plot_shapes.py

@@ -48,6 +48,10 @@ img[rr,cc,2] = 255
 rr, cc = circle_perimeter(120, 400, 15)
 img[rr, cc, :] = (255, 0, 0)
 
+# anti-aliased circle
+rr, cc, val = circle_perimeter(120, 400, 70, 'wu')
+img[rr, cc, 1] = val * 255
+
 # ellipses
 rr, cc = ellipse_perimeter(120, 400, 60, 20, orientation=math.pi / 4.)
 img[rr, cc, :] = (255, 0, 255)

skimage/draw/_draw.pyx

@@ -175,9 +175,10 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
         Centre coordinate of circle.
     radius: int
         Radius of circle.
-    method : {'bresenham', 'andres'}, optional
+    method : {'bresenham', 'andres', 'wu'}, optional
         bresenham : Bresenham method (default)
         andres : Andres method
+        wu : Wu's method
 
     Returns
     -------
@@ -192,6 +193,8 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
     circles create a disc whereas Bresenham can make holes. There
     is also less distortions when Andres circles are rotated.
     Bresenham method is also known as midpoint circle algorithm.
+    Wu's method draws anti-aliased circle. This implementation doesn't use
+    lookup table optimization.
 
     References
     ----------
@@ -199,6 +202,7 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
            plotter", 4 (1965) 25-30.
     .. [2] E. Andres, "Discrete circles, rings and spheres",
            18 (1994) 695-706.
+    .. [3] X. Wu, "Fast anti-aliased circle generation", 2 (1995) 446-450.
 
     Examples
     --------
@@ -222,10 +226,15 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
 
     cdef list rr = list()
     cdef list cc = list()
+    cdef list val = list()
 
     cdef Py_ssize_t x = 0
     cdef Py_ssize_t y = radius
     cdef Py_ssize_t d = 0
+
+    cdef double dceil = 0
+    cdef double dceil_prev = 0
+
     cdef char cmethod
     if method == 'bresenham':
         d = 3 - 2 * radius
@@ -233,33 +242,62 @@ def circle_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t radius,
     elif method == 'andres':
         d = radius - 1
         cmethod = 'a'
+    elif method == 'wu':
+        cmethod = 'w'
     else:
         raise ValueError('Wrong method')
 
-    while y >= x:
-        rr.extend([y, -y, y, -y, x, -x, x, -x])
-        cc.extend([x, x, -x, -x, y, y, -y, -y])
+    if cmethod == 'a' or cmethod == 'b':
+        while y >= x:
+            rr.extend([y, -y, y, -y, x, -x, x, -x])
+            cc.extend([x, x, -x, -x, y, y, -y, -y])
 
-        if cmethod == 'b':
-            if d < 0:
-                d += 4 * x + 6
-            else:
-                d += 4 * (x - y) + 10
-                y -= 1
+            if cmethod == 'b':
+                if d < 0:
+                    d += 4 * x + 6
+                else:
+                    d += 4 * (x - y) + 10
+                    y -= 1
+                x += 1
+            elif cmethod == 'a':
+                if d >= 2 * (x - 1):
+                    d = d - 2 * x
+                    x = x + 1
+                elif d <= 2 * (radius - y):
+                    d = d + 2 * y - 1
+                    y = y - 1
+                else:
+                    d = d + 2 * (y - x - 1)
+                    y = y - 1
+                    x = x + 1
+        return (np.array(rr, dtype=np.intp) + cy,
+                np.array(cc, dtype=np.intp) + cx)
+
+    elif cmethod == 'w':
+        dceil_prev = 0
+
+        rr.extend([y, x,  y,  x, -y, -x, -y, -x])
+        cc.extend([x, y, -x, -y,  x,  y, -x, -y])
+        val.extend([1] * 8)
+
+        while y > x + 1:
             x += 1
-        elif cmethod == 'a':
-            if d >= 2 * (x - 1):
-                d = d - 2 * x
-                x = x + 1
-            elif d <= 2 * (radius - y):
-                d = d + 2 * y - 1
-                y = y - 1
-            else:
-                d = d + 2 * (y - x - 1)
-                y = y - 1
-                x = x + 1
+            dceil = math.sqrt(radius**2 - x**2)
+            dceil = math.ceil(dceil) - dceil
+            if dceil < dceil_prev:
+                y -= 1
+            rr.extend([y, y - 1, x, x, y, y - 1, x, x])
+            cc.extend([x, x, y, y - 1, -x, -x, -y, 1 - y])
 
-    return np.array(rr, dtype=np.intp) + cy, np.array(cc, dtype=np.intp) + cx
+            rr.extend([-y, 1 - y, -x, -x, -y, 1 - y, -x, -x])
+            cc.extend([x, x, y, y - 1, -x, -x, -y, 1 - y])
+
+            val.extend([1 - dceil, dceil] * 8)
+            dceil_prev = dceil
+
+        return (np.array(rr, dtype=np.intp) + cy,
+                np.array(cc, dtype=np.intp) + cx,
+                np.array(val, dtype=np.float))
 
 
 def ellipse_perimeter(Py_ssize_t cy, Py_ssize_t cx, Py_ssize_t yradius,